boost::uBlas: Incomplete LU für dünn besetzte Matrix

Danke für den Link, das schaue ich mir mal an. Meine Systemmatrix ist leider nicht richtig sparse, nur fast. viele Einträge sind sehr klein, vernachlässigen will ich sie aber nicht. Mit Matlab habe ich mal ausprobiert, alle Einträge mit Betrag kleiner Epsilon zu Null zu setzen, dann eine ILU zu machen und dann GMRES auf die volle Matrix mit dem spärlichen Präkonditionierer loszulassen. Das hat sehr fein funktioniert (sogar überraschend gut: Meine Systematrix ist ansonsten total garstig; iterative Löser tun alles, nur nicht in endlicher Zeit ohne Präkonditionierer zu konvergieren)

Über weitere Tipps und Hinweise freue ich mich natürlich

darkfate

Ist sie wenigstens symmetrisch? Sonst macht ja ein Sparse-Algorithmus keinen Sinn.
Je nachdem wie stark du filterst, bekommst du ein rauschen in deine Lösung. Unter bestimmten Umständen könnte deine Determinante sogar Null sein und mit der Vernachlässigung ist sie dann doch wieder ungleich Null und sehr klein. Damit produzierst du mMn falsche Ergebnisse.

Hat deine Matrix ein bestimmtes Muster? Ich kann mir schlecht vorstellen was in deinem Solver vorgeht

Die Matrix ist reell und nicht symmetrisch. Sie entsteht in einem BEM-Code. Grob kann man sie sich als Interaktion vorstellen von Oberflächenpaneelen und Quellpunkten. Sind Paneel und Quellpunkt weit entfernt, wird der entsprechende Matrixeintrag klein. Vernachlässigen tue ich Elemente nur bei der Berechnung des Präkonditionierers:

Algorithmus (in etwa)

Erstelle vollbesetzte Matrix A und rechte Seite b
Erstelle dünn besetzte Matrix B durch Entfernen kleiner Einträge aus A
[L,U] = ilu(B)
x = gmres(A,b, ..., L,U)

Im Moment geht es mir erstmal unr darum, zu testen, wie gut mein BEM-Code mit iterativen Lösern umgehen kann.

Warum muss die Matrix symmetrisch sein? Wie gesagt, ich habe das ganz ja schon unter Matlab ausprobiert und es funktioniert anständig. Ich möchte es bloß in meinem eigentlichen Code einbauen, und dafür brauche ich eine anständige C++ ILU implementierung, die (am einfachsten für mich) mit boost::ublas zusammenarbeiten kann.

darkfate

Mups schrieb:

Warum muss die Matrix symmetrisch sein? Wie gesagt, ich habe das ganz ja schon unter Matlab ausprobiert und es funktioniert anständig. Ich möchte es bloß in meinem eigentlichen Code einbauen, und dafür brauche ich eine anständige C++ ILU implementierung, die (am einfachsten für mich) mit boost::ublas zusammenarbeiten kann.

Nein muss sie natürlich nicht, du hast nur im Topic von dünn besetzen Matrizen gesprochen. Ich weiß überhaupt nicht ob es in deinem Fall einen Leistungsschub bringt. Wie hoch ist die Prozentzahll der Nulleinträge im Schnitt?

Kennst du LAPACK? Ist zwar Fortran... soll aber auch C Ports geben.
Ich benutze daraus manche Funktionen übersetze sie selber in C obwohl ich nicht begabt bin klappt es ganz gut. Wenn das nicht klappt, dann JScience in Java... da kann man eventuell auch ein wenig Code übernehmen.

LU ist auch nicht wirklich ein schwerer Algorithmuss.. ich schätze mal er wird nicht länger als 50-100 Zeilen lang sein.

Dann gibt es noch "Numerical Recipies in C" leider nur als Buch mit CD-ROM

http://www.amazon.com/s/ref=nb_ss?url=search-alias%3Dstripbooks&field-keywords=numerical+recipes&x=0&y=0

Es gibt zwar viele offene Implementierungen. Diesen sollte man aber nicht blind vertrauen weil sie eben nicht sicher bei Optimierungen laufen. Zumindest ist der Großteil der Mathematik-Programmierer der Meinung.

Ich benutze zur Zeit Lapack, um das volle System zu lösen und möchte nun auf iterative Löser umsteigen (mit denen kann ich hoffentlich irgendwann die Struktur der vollen Matrix besser ausnutzen - ganz anderes Thema). Mir geht es nicht um eine LU (das kann Lapack), sondern um eine ILU: *incomplete* LU (so heisst auch das Topic).
Sowas selber implementieren ist doof. Feritge, optimierte Löser sind in der Regel um Größenordnungen schneller als eigener Fricklkram (habe ich sehr deutlich beim testen von ATLAS oder der MKL gemerkt).
Und was den Leistungsschub angeht: Das habe ich ja schon mit Matlab schon getestet. Das geht es anständig, und vor allem geht es mir ja gerade darum, vom direkten auf einen iterativen Löser umzusteigen. Dafür der Präkonditionierer. Dafür die incomplete LU.
Das man *irgendwelchen* Implementierungen nicht vertrauen soll, ist mir klar. Deswegen frage ich ja hier

Also zurück zum Topic: Kann jemand von euch eine ILU implementierung (die am besten mit ublas zuammenarbeitet) empfehlen oder kennt zumindest eine?

Wir die ILU auch in "Numerical Recipies in C" beschrieben?

darkfate

Ich bin mir fast sicher das du bei PARDISO zwischen iterativ und direkt umschalten kannst. MKL Pardiso hat es in jedem Fall.

Walli

Mups schrieb:

iterative Löser tun alles, nur nicht in endlicher Zeit ohne Präkonditionierer zu konvergieren

GMRES konvergiert schon wenn die Matrix vollen Rang hat, denn GMRES ist letztlich ein direkter Löser.

boost::ublas ist meiner Erfahrung nach ziemlich lahm. Ich benutze PETSc für solche Sachen. Man hat viele Möglichkeiten und man kann den relativ objektorientierten C-Kram auch schön hinter einem C++ Wrapper verstecken, wenn man ihn nicht direkt benutzen mag.

Walli schrieb:

Mups schrieb:

iterative Löser tun alles, nur nicht in endlicher Zeit ohne Präkonditionierer zu konvergieren

GMRES konvergiert schon wenn die Matrix vollen Rang hat, denn GMRES ist letztlich ein direkter Löser.

Das habe ich mit "in endlicher Zeit" andeuten wollen. Mit jedem Iterationsschritt wird GMRES langsamer und verbraucht mehr Speicher, was in der Praxis die Anwendbarkeit als direkten Löser minimiert.

Walli

Naja, ein ordentlicher Vorkonditionierer und ein restarted GMRES sind schon in der Regel deutlich performanter .