Exponentialmatrix einer Skew-Hermitian-Matrix mit Fortran 95 und LAPACK

11

Ich stecke mich gerade für einige quantenmechanische Simulationen in fortran 95. Ehrlich gesagt wurde ich von Octave verwöhnt, daher habe ich die Matrixexponentiation als selbstverständlich angesehen. Was ist bei einer (kleinen, ) Schräg- Hermitianischen Matrix der Größe der effizienteste Weg, um dieses Problem mit LAPACK zu lösen? Ich verwende nicht den LAPACK95-Wrapper, sondern nur direkte Anrufe an LAPACK. $n\leq 36$ $n\times n$

fortran lapack qubyte
quelle

2

Benötigen Sie die Exponentialmatrix für sich oder benötigen Sie die Exponentialmatrix multipliziert mit einem Vektor?

Paul

@ Paul: Sorry, habe das vorher nicht gesehen. Nein, ich brauche die gesamte Matrix.

Qubyte

Warum sollte jemand diese Frage ablehnen? Wenn Sie abstimmen, hinterlassen Sie bitte einen Grund in den Kommentaren! Vielleicht kann die Frage auf diese Weise verbessert werden.

Qubyte

Wir verlassen uns auf DGPADM , aber bei Jack Poulson gibt es einen besseren Weg.

Mike Dunlavey

16

Matrixexponentiale von schief-hermitischen Matrizen sind billig zu berechnen:

Angenommen, ist Ihre schief- Matrix, dann ist , und über zheevd und Freunde können Sie die Zerlegung erhalten $A$ $iA$

i A = U Λ U^{H},

$iA = U \Lambda U^H,$

Dabei ist die einheitliche Eigenvektormatrix und ist real und diagonal. Dann, trivial, $U$ $\Lambda$

A = U (- i Λ) U^{H} .

$A = U (-i \Lambda) U^H.$

Sobald Sie und , ist es einfach zu berechnen $U$ $\Lambda$

\exp (A) = \exp (U (- i Λ) U^{H}) = U \exp (- i Λ) U^{H}

$\exp(A)=\exp(U (-i\Lambda) U^H)=U \exp(-i\Lambda) U^H$

indem Sie zuerst die Eigenwerte potenzieren , über zcopy setzen , ausführen, indem Sie zscal für jede Spalte mit einem potenzierten Eigenwert ausführen und schließlich Ihr Ergebnis auf setzen $B := U$ $B := B \exp(-i \Lambda)$

\exp (A) := B U^{H}

$\exp(A) := B U^H$

über zgemm .

Jack Poulson
quelle

Vielen Dank! Ich habe dort einen offensichtlichen Trick mit dem verpasst . Sie haben mich auf die spezifischen LAPACK-Unterprogramme gesetzt, die ich benötige. Vielen Dank dafür. Ich werde dies noch nicht als korrekt markieren (möchte es zuerst testen).

i

$i$

Qubyte

1

Keine Eile. Ich habe es tatsächlich schon einmal implementiert, also bin ich ziemlich zuversichtlich :-)

Jack Poulson

Dies wird einer dieser magischen Code-Teile sein, die ich überall verwende. Für das, was es wert ist, werde ich mich auch in einer Kommentarzeile bedanken, die wahrscheinlich niemand sonst jemals sehen wird.

Qubyte

2

@ JackPoulson: Gut gespielt, Sir. Das bekomme ich, wenn ich einen Major auswähle, der nicht an imaginäre Zahlen glaubt (außer an Eigenwerte).

Geoff Oxberry

1

@ JackPoulson: Es funktioniert wunderbar. Nochmals vielen Dank dafür. Besonders das zscal-Bit. Ich hatte den größten Teil des restlichen Codes in einer anderen Unterroutine, aber das war etwas, das ich übersehen hatte.

Qubyte

5

Da ich auf meinem Handy bin, kann ich die Dinge nicht einfach verknüpfen und werde später Links hinzufügen. Sie sollten sich wahrscheinlich das Papier "19 zweifelhafte Methoden zur Berechnung des Matrixexponentials", die Fortran-Bibliothek EXPOKIT, Jitse Niesens Artikel über Krylov-Methoden zur Berechnung des Matrixexponentials und einige der jüngsten Artikel von Nick Higham über Matrixexponentiale ansehen. Es ist üblicher, das Produkt einer Exponentialmatrix und eines Vektors als das Exponential der Matrix allein zu benötigen, und hier können Krylov-Methoden sehr hilfreich sein. Für kleinere, dichte Matrizen wie die von Ihnen beschriebenen sind Padé-Methoden vielleicht besser, aber ich habe mit Krylov-Methoden viel Erfolg gehabt, wenn sie in exponentiellen Methoden zur numerischen Integration von ODEs verwendet wurden.

Geoff Oxberry
quelle

Vielen Dank. Ich kenne 19 zweifelhafte Wege und auch Expokit, aber einige der Leute, mit denen ich zusammenarbeite, sind in der Industrie, deshalb möchte ich dies aus urheberrechtlichen Gründen vermeiden. Ich bin sehr daran interessiert, es mit LAPACK / BLAS zu implementieren, da ich bereits auf diese Bibliotheken verlinke. Eines jedoch; Ich brauche die Matrix exponentiell selbst. Ich arbeite an einer Variante der Quantenprozesstomographie, und der betreffende Prozess wird durch die Matrix verkörpert. Später werde ich mich mit einem Integrator in Kombination mit dieser Exponentialmatrix befassen, wenn es wirklich interessant wird!

Qubyte

1

Der komplexe Eigenlösungsansatz ist mathematisch korrekt, macht aber mehr Arbeit als nötig. Leider kann der verbesserte Ansatz, den ich beschreiben werde, nicht mit LAPACK-Aufrufen implementiert werden.

Schauen Sie sich RC Ward und LJ Gray, ACM Trans an. Mathematik. Sanft. 4, 278 (1978). Dies beschreibt die Software, die im TOMS-Algorithmus 530 verfügbar ist und die Sie von netlib herunterladen können. Dies beschreibt, wie die symmetrische Versatzmatrix als faktorisiert wird $X$

X = U D U^{T}

$X = U D U^T$

wobei real orthogonal ist und real schiefsymmetrisch und blockdiagonal ist. Die diagonalen Unterblöcke sind entweder oder . Da es sich um eine Blockdiagonale handelt, können Sie jeden Unterblock separat potenzieren. Die Blöcke sind Null und , also sind diese trivial. Die Unterblöcke sind fertig $U$ $D$ $2\times 2$ $1\times 1$ $1\times 1$ $\exp(0)=1$ $2\times 2$

\exp (\begin{matrix} 0 & - t \\ t & 0 \end{matrix}) = (\begin{matrix} \cos t & - \sin t \\ \sin t & \cos t \end{matrix})

$\exp \begin{pmatrix} 0 & -t \\ t & 0 \end{pmatrix} = \begin{pmatrix} \cos t & -\sin t \\ \sin t & \cos t \end{pmatrix}$

Die gewünschte Exponentialmatrix ist dann gegeben durch

\exp (X) = U \exp (D) U^{T}

$\exp(X) = U \exp(D) U^T$

Ich habe diesen Ansatz in meinen quantenchemischen Codes seit mehreren Jahrzehnten verwendet und hatte nie Probleme mit der beteiligten Software.

Ron Shepard
quelle

Hallo @ Ron Shepard, und willkommen bei Computational Exchange SE. Können Sie Ihre zweite und dritte Gleichung bearbeiten? Sie sind etwas schwer zu verstehen.

Nicoguaro

0

Wenn Sie nur die Exponentialmatrix multipliziert mit einem Vektor benötigen, kann diese fortran-Subroutine für Sie von Nutzen sein. Es berechnet:

$(e^A)v$

Dabei ist v ein Vektor und A eine reguläre Einsiedlermatrix. Es ist eine Unterroutine aus der EXPOKIT- Bibliothek

Andernfalls möchten Sie möglicherweise diese Unterroutine berücksichtigen , die für jede allgemeine komplexe Matrix A funktioniert.

Paul
quelle

Das sieht nicht nach einem Verweis auf Fortran-Bibliotheken aus.

Geoff Oxberry

@ GeoffOxberry: Ich habe es umgeschrieben, um die fortran-Subroutinen einzuschließen

Paul

@ Paul: Nein gut, fürchte ich. Was ich mache, ist eine All-Matrix-Variante der Prozesstomographie. Außerdem Skew- Hermitian!

Qubyte

Ich weiß es zu schätzen, dass Sie Ihre Antwort umgeschrieben haben, aber basierend auf dem Bearbeitungspfad haben Sie anscheinend Ihre Antwort vollständig geändert, Elemente meiner chronologisch früheren Antwort übernommen und Links hinzugefügt.

Geoff Oxberry

@GeoffOxberry: Im Gegenteil ... Meine Ergebnisse kamen unabhängig von Ihren, aber Sie haben gepostet, bevor ich die Möglichkeit hatte, meine Antwort neu zu schreiben :)

Paul

Exponentialmatrix einer Skew-Hermitian-Matrix mit Fortran 95 und LAPACK

Antworten: