rozwiąż

12

Przenoszę istniejący kod z MATLAB do C ++ i mam liniowy system do rozwiązania (zamiast bardziej typowej formy ) $xA=b$ $Ax=b$

Matryca jest gęsta i ma ogólną formę, ale nie jest większa niż 1000 x 1000. Tak więc w MATLAB rozwiązaniem jest funkcja lub notacja ukośnika do przodu $A$ mrdivide(b,A)x = b/A;

Jak mam rozwiązać ten problem w kodzie C ++ przy użyciu procedur BLAS i LAPACK?

Znam procedurę LAPACK, DGESVktóra rozwiązuje dla . $Ax=b$ $x$

Jedną z moich myśli było wykonanie pewnych manipulacji przy użyciu macierzy transponujących tożsamość:

$(xA)^T=b^T$

$A^T x^T = b^T$

$x^T = (A^T)^{-1} b^T$

Następnie rozwiązać ostateczny kształt za pomocą DGESVdziałających na transponowana . (więc koszt transpozycji i koszt rozwiązania systemu) $A^T$ $A$

Czy istnieje podejście bardziej wydajne lub w inny sposób lepsze ?

Pracuję z klasami macierzy i wektorów, a także implementacją BLAS z biblioteki BOOST uBLAS, a także powiązaniami z procedurami biblioteki LAPACK. Korzystałem z tej konfiguracji z powodzeniem do innych operacji i mam nadzieję znaleźć rozwiązanie ograniczone do tych bibliotek.

Powinienem również zauważyć, że wykonuję ten typ operacji tylko kilka razy podczas konfiguracji kodu, więc wydajność nie jest krytyczna.

Może to MATLAB dokumentacja na mrdivideto pomocne dla innych.

linear-solver c++ lapack

— NoahR
źródło

10

Trywialna odpowiedź dla kwadratu : zastosowanie, które rozwiązuje również dla kiedy . $A$ dgesvx $A^T x = b$ TRANS = 'T'

Należy pamiętać, że w przypadku BLAS lub LAPACK prawie nie trzeba transponować (zamieniać elementów w pamięci) macierzy: większość podprogramów ma TRANSargument umożliwiający działanie na macierzy transpozycji lub na macierzy przechowywanej z innym układem pamięci. (Transpozycja jest równoważna zmianie układu pamięci ciągłej Fortran na C-contiguos i viceversa.)

— Stefano M.
źródło

Dzięki za odpowiedź i wyjaśnienie! Z LAPACK wykonałem bardzo mało pracy i teraz wiem, jak szukać opcji TRANS. Mam problem z przepracowaniem argumentu TRANS boost::numeric::bindings::lapack::gesvx(), ale to nie jest część mojego pytania tutaj. Jeśli mi się uda, wrócę z notatką, jak to zrobić.

— NoahR

Mam działające rozwiązanie gesvx(), ale nie bez przeszkód po drodze. Gdy argumentem TRANS jest „T”, dokumentacja LAPACK mówi, że gesvxrozwiązuje

, ale tak naprawdę rozwiązuje

ponieważ oczekuje się, że postać argumentów wejściowych

i

nie będzie transponowana Formularz. Tak więc argument

jest transponowany, podczas gdy

i

A^{T} X = B

$A^T X = B$

A^{T} X^{T} = B^{T}

$A^T X^T = B^T$

X

$X$

B

$B$

A

$A$

X

$X$

B

$B$ nie są. Świetnie, to wygodniejsze. Jeśli ktoś inny natknie się na to, próbując użyć powiązań numerycznych boost, powiem, że nie byłem w stanie uzyskać interfejsu transpozycji używanego w tym rozwiązaniu. pracować przez powiązania.

— NoahR

Ok, znalazłem sposób na skorzystanie z gesvxformularza transpozycji boost::numeric::bindings. Zawiń transponowaną macierz

w funkcji. To identyfikuje parametr jako typ transpozycji :

A^{T}

$A^T$ trans()

boost::numeric::bindings::lapack::gesvx( FACT, boost::numeric::bindings::trans(Atransposed), af, ipiv, equed, r, c, b, x, rcond, ferr, berr );

— NoahR

0

Możesz obliczyć pseudo odwrotność , używając powiedzmy dekompozycji SVD, które są zawarte w LAPACK. $A$

x ZA = b x Q R = b x = b R^{- 1} Q^{T.}

$\begin{equation} xA=b\\ xQR=b\\ x=bR^{-1}Q^T \end{equation}$

$A$

— Gil
źródło

3

A

$A$

R

$R$

R^{- 1}

$R^{-1}$