FORTRAN90: Matrizen-System loesen (LAPACK)

[zahor]

Hi DPler,

tut mir Leid wegem der Textmenge, aber mein Problem ist folgendes:
Und zwar uebersetzte ich gerade ein Script von Matlab nach Fortran.
Von mir aus koennte es ja Matlab bleiben, aber ich programmiere das ja nicht fuer mich,
sondern fuer ein physikalisches Projekt an nem Forschungsinstitut.
Die Matlab-Stelle sieht so aus:

Code:

A=B\M

Dabei sind sowohl A, B als auch Matrizen vom Typ Double Complex (Komplexe Zahlen mit erhoehter Genauigkeit).
Ihre Groesse sei n x n.
Ich hab mir dann gedacht - klar, ist ja logisch, a=b\m heisst soviel wie a=INVERSE(b)*m - so ist der Matlab-Operator \ definiert. also invertiere ich zuerst Matrix B und multipliziere sie dann mit M. Fuer die Multiplikation hat ifort (der verwendete Fortran-compiler) eine eingebaute Funktion. Fuer die Invertierung habe ich folgende LAPACK-Routinen gefunden, und auch versucht, sie zu verwenden.

Code:

! LU-Factorisation of M
call zgetrf(n, n, b, n, ipiv, info)
PRINT *, "Successfully LU-Factorised Matrix 2", info

! Invert it
call zgetri(n, b, n, ipiv, work, lwork, info)
PRINT *, "Successfully inverted Matrix 2!", info

A=MATMUL(b,m)
PRINT *,"Successfully multiplied matrices"

dabei hat IPIV - ein Integer-Array - die Dimension n, und lwork (typ integer) den Wert 1594, was gleichzeitig auch die Dimension von Work ist - das ist uebrigens 2n. info ist ein Ausgabewert (integer), sozusagen der Fehlercode.
nun gibt info aber bei den beiden oberen Werten 28 zurueck.
Die Dokumentation zu den obigen befehlen findet ihr hier fuer zgetrf und hier fuer zgetri. Die Links sind zwar fuer die Variante mit Single Precision, das macht bei den Parametern aber keinen Unterschied.
Fehlercode 28 sagt mir jetzt also, dass meine Matrix ungeeignet ist. Die Frage ist nun - wie aendere ich die Matrix so, dass die Funktion sie annimmt?
Der Matlab-Dokumentation (klick) ist zu entnehmen, dass es dazu die Routinen ZLANGE, ZGESV, ZGECON verwendet. ZGESV mag meine Matrix aber auch nicht und gibt auch 28 als Fehlercode zurueck. mein Problem mit ZLANGE und ZGECON ist allerdings, dass ich nicht weiss, welche der verschiedenen Anwendungsmoeglichkeiten die richtige ist! Oder sind das gar nicht die "richtigen" Funktionen? sollte ich vllt. was ganz anderes probieren? jede hilfe ist willkommen!

[calculon]

Bei Matrizen die nahe singulär sind nimmt man zum Lösen von linearen Gleichungssystemen auch gerne die Singulärwertzerlegung. Vielleicht ist die in deinem Pack ja auch enthalten.

Btw.: Als ich 16 war, hab' ich glaub' ich gerade gelernt wie man Polynomdivision macht; wenn überhaupt ;-)

Gruß
--

[zahor]

joa, ueber SVD koennte man es auch berechnen. Danke fuer die Idee! Ich hab mal gegoogelt, wie ich das mit LAPACK anstellen muesste, und bin auch fuendig geworden. Das wuerde dann aber ziemlich umstaendlich:
1. ZGEBRD aufrufen, um die Matrix in bidiagonale form zu bringen
2. ZUNGBR aufrufen (2 mal), um Q sowie P^H (P**H) zu kriegen
3. Dann das ganze durch ZBDSQR jagen um die Singulaerwerte zu erhalten
4. Letzendlich das Pseudinverse mithilfe von ZGELSS und den Singulaerwerten erhalten
5. Dieses Ergebnis mithilfe von MATMUL mit der anderen Matrix multiplizieren.
Da wird man doch doof davon, vor allem weil das ja ungemein viel Zeit beansprucht! Da wird ja alles moegliche berechnet, was ich weiter dann gar nicht brauche. Da ist ZGESV wohl doch schneller. Der Fakt, das MATLAB es fuer seine Berechnungen benutzt, will ja doch was heissen. Die werden ja alle Varianten durchdacht haben.
Ich wuerde schon lieber bei ZGESV oder notfalls ZGERTF / ZGETRI / ZGETRS bleiben... Danke trotzdem
Meine Frage war eher so gedacht: Wie mache ich ZGESV meine Matrizen schmackhaft?

[calculon]

Zitat:

Zitat von MatLab

The specific algorithm used for solving the simultaneous linear equations denoted by X = A\B and X = B/A depends upon the structure of the coefficient matrix A.

Welchen Algorithmus MatLab nimmt, weißt du eigentlich nicht genau, da MatLab das erst mal durch Analyse der Matrix selbst entscheiden muss. Selbst bei einer singulären Matrix wird noch versucht das Gleichungssystem zu lösen.
Daher solltest du nicht zu sehr kucken wie MatLab das eigentlich macht, sondern selber kucken wie du das Lösen kannst (ist die Matrix symmetrisch?, singulär?, sparse?, etc.) und meine Erfahrung ist, dass man mit Hilfe der Singulärwertzerlegung einen sehr universellen Löser hat. Aber du hast Recht, es ist sehr lahm.

Gruß
--

[zahor]

also die matrix ist weder symmetrisch noch sparse. aber faktor 28 ist laut lapack singulaer.
deshalb geht ja auch die "normale" routine nicht -.- und ich will versuchen, sie zu desingularisieren (heisst das so?^^)
oder eben das ganze ueber das pseudoinverse loesen... oder per cgesv, was aber anscheinend wieder eine nicht-singulaere matrix will...
und ich bin den "Algorithm"-Teil der matlab-doku-seite einfach von oben nach unten durchgegangen und hab die bedingungen geprueft...
im anhang mal ein Screenshot von

Code:

spy(M)

und

Code:

spy(B)

[calculon]

Also wenn die Matrizen nicht sparse sind, dann weiß' ich auch nicht ;-)

Gruß
--

[zahor]

die ist nicht als sparse erstellt, sondern als volle matrix. von den werten her ist sie natuerlich sparse.
das problem ist uebrigens geloest. die matrizen waren zu gross, irgendwo im code wurde von der variable - nennen wirsie mal n -, die vorher fuer die dimension verwendet wird (vergleichbaer mit setlength(matrix,n) - nur halt als matrix und nich als vektor) um 2 verringert, aber die groesse der matrix bleibt gleich. dann wird n als parameter, der die groesse beschreibt, verwendet, was natuerlich nicht stinmmt. jetzt geht alles.
dennoch danke fuer euere hilfe :-)

[calculon]

Zitat:

Zitat von zahor

[..] dennoch danke fuer euere hilfe :-)

Bitte gern geschehen, aber du musst mich deshalb nicht in der dritten Person anreden :P

Gruß
--

[zahor]

ja stimmt eigentlich ;-)