Macierz wykładniczy prawdziwej macierzy asymetrycznej z Fortran 95 i LAPACK

Niedawno zadałem pytanie w tym samym stylu dla matryc skośno-hermitowskich. Zainspirowany sukcesem tego pytania i po kilku godzinach uderzenia głową o ścianę patrzę na wykładniczą macierz prawdziwych matryc asymetrycznych. Droga do znalezienia wartości własnych i wektorów własnych wydaje się dość skomplikowana i obawiam się, że się zgubiłem.

Tło: Jakiś czas temu zadałem to pytanie na temat fizyki teoretycznej SE. Wynik pozwala mi sformułować równania wzorcowe jako rzeczywiste macierze asymetryczne. W przypadku niezależnym od czasu równanie wzorcowe rozwiązuje się przez wykładnik tej macierzy. W przypadku zależnym od czasu będzie wymagało integracji. W tej chwili interesuje mnie tylko niezależność czasowa.

Po patrząc na różnych podprogramów myślę, że powinno być wywołanie ( ? Gehrd , ? Orghr , ? Hseqr ...) nie jest jasne, czy byłoby prostsze do oddania matrycę od real*8celu complex*16i postępować ze skomplikowanymi podwójnymi wersje tych funkcji, lub trzymaj się real*8i podbij trafienie podwojenia liczby moich tablic i zrobienia ich złożonej macierzy później.

Więc do jakich procedur mam dzwonić (i w jakiej kolejności) i czy powinienem używać prawdziwych podwójnych wersji czy złożonych podwójnych wersji? Poniżej znajduje się próba zrobienia tego z prawdziwymi podwójnymi wersjami. Utknąłem, szukając wartości własnych i wektorów własnych L*t.

function time_indep_master(s,L,t)
  ! s is the length of a side of L, which is square.
  ! L is a real*8, asymmetric square matrix.
  ! t is a real*8 value corresponding to time.
  ! This function (will) compute expm(L*t).

  integer, intent(in)    :: s
  real*8,  intent(in)    :: L(s,s), t
  real*8                 :: tau(s-1), work(s), wr(s), wi(s), vl
  real*8, dimension(s,s) :: time_indep_master, A, H, vr
  integer                :: info, m, ifaill(2*s), ifailr(2*s)
  logical                :: sel(s)

  A = L*t
  sel = .true.

  call dgehrd(s,1,s,A,s,tau,work,s,info)
  H = A
  call dorghr(s,1,s,A,s,tau,work,s,info)
  call dhseqr('e','v',s,1,s,H,s,wr,wi,A,s,work,s,info)
  call dhsein('r','q','n',sel,H,s,wr,wi,vl,1,vr,s,2*s,m,work,ifaill,ifailr,info)

  ! Confused now...

end function

Najpierw zastanowiłbym się, czy matryca jest naprawdę całkowicie dowolna: czy jest jakaś transformacja, która uczyniłaby ją pustelnikiem? Czy fizyka gwarantuje, że matryca powinna być diagonalizowana (z odpowiednio uwarunkowaną matrycą wektorów własnych)?

Jeśli okaże się, że tak naprawdę nie ma żadnej symetrii do wykorzystania, powinieneś zacząć od przeczytania Dziewiętnastu podejrzanych sposobów obliczania wykładniczej macierzy , która jest standardowym odniesieniem (i została napisana przez autora MATLAB i współautora G & vL) .

Opierając się na tym, co powiedział Jack, standardowym podejściem, które wydaje się być stosowane w oprogramowaniu (np. EXPOKIT, wspomnianym we wcześniejszym pytaniu) jest skalowanie i kwadratowanie, a następnie aproksymacja Padé (metody 2 i 3) lub metody podprzestrzeni Kryłowa (metoda 20). W szczególności, jeśli patrzysz na integratory wykładnicze, powinieneś rozważyć metody podprzestrzeni Kryłowa i spojrzeć na artykuły na temat integratorów wykładniczych (niektóre odniesienia są wymienione wraz z Metodą 20 w pracy Moler & Van Loan).

Jeśli zależy ci na użyciu wektorów własnych, rozważ użycie trójkątnych układów wektorów własnych (Metoda 15); ponieważ macierz może być niediagonalizowalna, to podejście może nie być najlepsze, ale jest lepsze niż próba bezpośredniego obliczenia wektorów własnych i wartości własnych (tj. Metoda 14).

Redukcja do formy Hessenberga nie jest dobrym pomysłem (Metoda 13).

Nie jest dla mnie oczywiste, czy lepiej byłoby byś posługiwał się rzeczywistą czy złożoną arytmetyką, ponieważ złożona arytmetyka Fortrana jest szybka, ale może przepełniać / niedopełniać (patrz „O ile lepsze są naprawdę kompilatory Fortran?” ).

Możesz bezpiecznie zignorować Metody 5-7 (metody oparte na rozwiązaniu ODE są nieefektywne), Metody 8-13 (drogie), Metoda 14 (obliczanie wektorów własnych dużych macierzy jest trudne bez specjalnej struktury i podatne na błędy numeryczne w przypadkach źle uwarunkowanych) oraz Metoda 16 (obliczanie rozkładu Jordana macierzy jest niestabilne numerycznie). Metody 17–19 są trudniejsze do wdrożenia; w szczególności metody 17 i 18 wymagałyby więcej czytania. Metoda 1 jest rezerwową opcją skalowania i kwadratowania, jeśli przybliżenia Padé nie działają dobrze.

$B_{j}$

$\gamma_{j}$$j$$E_{j}$

Mam prosty podprogram Fortran, który oblicza wykładnik dowolnej macierzy. Sprawdziłem to pod kątem polecenia Matlaba i jest w porządku. Opiera się na skalowaniu i kwadratowaniu. Napisałem to kilka lat temu.

Chciałbym sfinalizować inną procedurę, taką jak te, które pobieram z gams.nist.gov. Ale jeszcze nie ma szczęścia.