Matrix exponentiell einer realen asymmetrischen Matrix mit Fortran 95 und LAPACK

10

Ich habe kürzlich eine ähnliche Frage für schief-hermitische Matrizen gestellt. Inspiriert vom Erfolg dieser Frage und nachdem ich ein paar Stunden lang meinen Kopf gegen eine Wand geschlagen habe, betrachte ich die Exponentialmatrix realer asymmetrischer Matrizen. Der Weg zum Finden der Eigenwerte und Eigenvektoren scheint ziemlich kompliziert zu sein, und ich fürchte, ich habe mich verlaufen.

Hintergrund: Vor einiger Zeit habe ich diese Frage zur theoretischen Physik SE gestellt. Das Ergebnis ermöglicht es mir, Master-Gleichungen als echte asymmetrische Matrizen zu formulieren. Im zeitunabhängigen Fall wird die Hauptgleichung durch Exponentiierung dieser Matrix gelöst. Im zeitabhängigen Fall ist eine Integration erforderlich. Im Moment geht es mir nur um Zeitunabhängigkeit.

Wenn ich mir die verschiedenen Subroutinen ansehe , die ich aufrufen sollte ( ? Gehrd , ? Orghr , ? Hseqr ...), ist es unklar, ob es einfacher wäre, die Matrix von real*8nach zu gießen complex*16und mit den komplexen Doppelversionen dieser Routinen fortzufahren. oder bleiben Sie dabei real*8und versuchen Sie, die Anzahl meiner Arrays zu verdoppeln und später eine komplexe Matrix daraus zu erstellen.

Also, welche Routinen sollte ich aufrufen (und in welcher Reihenfolge) und sollte ich die echten Doppelversionen oder die komplexen Doppelversionen verwenden? Im Folgenden finden Sie einen Versuch, dies mit echten Doppelversionen zu tun. Ich bin festgefahren, die Eigenwerte und Eigenvektoren von zu finden 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
qubyte
quelle

Antworten:

8

Ich würde zuerst wirklich gründlich darüber nachdenken, ob die Matrix wirklich völlig willkürlich ist oder nicht: Gibt es eine Transformation, die sie hermitisch machen würde? Garantiert die Physik, dass die Matrix diagonalisierbar sein sollte (mit einer vernünftig konditionierten Eigenvektormatrix)?

Wenn sich herausstellt, dass es wirklich keine Symmetrie gibt, die ausgenutzt werden kann, sollten Sie zunächst Neunzehn zweifelhafte Methoden zur Berechnung des Matrixexponentials lesen , das die Standardreferenz ist (und vom Autor von MATLAB und dem Mitautor von G & vL geschrieben wurde). .

Jack Poulson
quelle
1
Das Beste, was ich tun kann, ist, daraus eine Blockdiagonalmatrix mit asymmetrischen Blöcken zu machen. Das an sich ist allerdings sehr interessant. Die meisten dieser Blöcke sind , und ich kann sie nur analytisch lösen. Es gibt jedoch einen verbleibenden Block ohne auszunutzende Symmetrien. 2×24×4
Qubyte
1
Ich mag diese Antwort; Der unsymmetrische Fall weist genügend Fallstricke auf, die es wert sind, in Betracht gezogen zu werden, wenn es eine Formulierung Ihres Problems gibt, die zu symmetrischen Matrizen anstelle von unsymmetrischen führt.
JM
@ MarkS.Everitt: Du scheinst fast da zu sein ... wie groß sind die Matrizen? Wieder ~ 36 x 36?
Jack Poulson
In diesem Fall , aber es besteht die Möglichkeit, dass es bis zu . 36 × 3616×1636×36
Qubyte
2
@ MarkS.Everitt: Ihr Problem ist also jetzt effektiv, wie man 4x4-Matrizen potenziert. Dies ist ausreichend klein, damit eine asymptotische Analyse irrelevant ist, sodass die Antwort vollständig von den Werten abhängt. Ich kann es nicht mehr wirklich sagen, es sei denn, Sie übersetzen Ihren verknüpften Physikbeitrag in lineare Algebra (was ist ein Superoperator?!?).
Jack Poulson
7

Um auf dem aufzubauen, was Jack gesagt hat, ist der Standardansatz, der in Software verwendet zu werden scheint (wie EXPOKIT, der in Ihrer früheren Frage erwähnt wurde), Skalieren und Quadrieren, gefolgt von Padé-Approximation (Methoden 2 und 3) oder Krylov-Subraummethoden (Methode) 20). Insbesondere wenn Sie Exponentialintegratoren betrachten, sollten Sie die Krylov-Subraummethoden betrachten und Artikel über Exponentialintegratoren betrachten (einige Referenzen werden zusammen mit Methode 20 im Moler & van Loan-Artikel erwähnt).

Wenn Sie unbedingt Eigenvektoren verwenden möchten, sollten Sie dreieckige Eigenvektorsysteme verwenden (Methode 15). Da Ihre Matrix möglicherweise nicht diagonalisierbar ist, ist dieser Ansatz möglicherweise nicht der beste, aber es ist besser, als zu versuchen, die Eigenvektoren und Eigenwerte direkt zu berechnen (dh Methode 14).

Eine Reduktion auf die Hessenberg-Form ist keine gute Idee (Methode 13).

Mir ist nicht klar, ob Sie mit realer oder komplexer Arithmetik besser bedient werden können, da die komplexe Fortran-Arithmetik schnell ist, aber möglicherweise über- oder unterläuft (siehe "Wie viel besser sind Fortran-Compiler wirklich?" ).

Sie können die Methoden 5-7 (ODE-Solver-basierte Methoden sind ineffizient), die Methoden 8-13 (teuer) und die Methode 14 ignorieren (die Berechnung von Eigenvektoren großer Matrizen ist ohne spezielle Struktur schwierig und in schlecht konditionierten Fällen anfällig für numerische Fehler). und Methode 16 (die Berechnung der Jordan-Zerlegung einer Matrix ist numerisch instabil). Die Methoden 17-19 sind schwieriger zu implementieren. Insbesondere würden die Methoden 17 und 18 mehr Lesen erfordern. Methode 1 ist eine Ersatzoption zum Skalieren und Quadrieren, wenn Padé-Näherungen nicht gut funktionieren.

Edit # 1 : Basierend auf den Kommentaren als Antwort auf Jacks Antwort scheint die Blockdiagonalisierung eine Option zu sein. In diesem Fall ist etwas wie Methode 18 (Blockdreieck-Diagonalisierung) eine sehr gute Methode. Ich habe zu Beginn gezögert, es zu empfehlen, da Ihre Frage diese Struktur nicht spezifiziert hat. Wenn Sie jedoch eine Transformation haben, die Ihre Matrix diagonalisiert, wird der Ansatz größtenteils komplexer. Sie sollten nur sicherstellen, dass Sie GW Stewarts Trick verwenden, jede Blockdiagonalmatrix in zu zerlegenBj

Bj=γjI+Ej,

Dabei ist der Durchschnitt der Eigenwerte der ten Blockdiagonalmatrix. Diese Zerlegung macht nahezu nullpotent, was die Genauigkeit Ihrer Exponentialberechnungen für die Matrix verbessert. Dieser Trick wird auf Seite 26 der Version des von Jack verlinkten Moler & van Loan-Papiers erläutert. j E jγjjEj

Geoff Oxberry
quelle
1
Upvoted, aber die Leute, die LAPACK implementiert haben, sind nicht naiv in Bezug auf komplexe Arithmetik, insbesondere angesichts der Zeit, die Kahan in die Analyse investiert hat. Er sollte jedoch immer noch in reeller Arithmetik rechnen, wenn dies Flops spart. Die Kosten für die spätere Umwandlung betragen nur gegenüber der gesamten -dichten linearen Algebra. O ( n 3 )O(n2)O(n3)
Jack Poulson
Zweifellos wissen sie, was sie tun; Ich mache mir keine Sorgen um die Implementierung von LAPACK. Ich bin mehr überrascht über das Verhalten des Fortran-Compilers.
Geoff Oxberry
2
Ja, der Compiler ist möglicherweise eher ein Problem als gut geschriebenes LAPACK. Es kann beunruhigend sein, festzustellen, dass Ihr Programm nur fehlgeschlagen ist, weil die vom Compiler verwendeten Implementierungen für Absolutwert und Division verpfuscht wurden ...
JM
-1

Ich habe eine einfache Fortran-Subroutine, die den Exponenten einer beliebigen Matrix berechnet. Ich habe es mit dem Matlab-Befehl verglichen und es ist in Ordnung. Es basiert auf Skalierung und Quadrierung. Ich habe es vor ein paar Jahren geschrieben.

Ich hätte gerne eine andere Unterroutine, wie die, die ich von gams.nist.gov heruntergeladen habe. Aber noch kein Glück.

Freejack
quelle