Ich verwende Numerische Bibliotheksbindungen für Boost UBlas, um ein einfaches lineares System zu lösen. Das Folgende funktioniert gut, außer dass es auf die Handhabung von Matrizen A (m x m) für relativ beschränkt ist kleines 'm'.
In der Praxis habe ich eine viel größere Matrix mit der Dimension m = 10 ^ 6 (bis zu 10 ^ 7).
Gibt es einen C ++ - Ansatz zum Lösen von Ax = b, der Speicher effizient verwendet.
Kurze Antwort: Verwenden Sie nicht Boosts LAPACK -Bindungen, diese wurden für dichte Matrizen entworfen,
Nicht spärliche Matrizen, verwenden Sie stattdessen UMFPACK .
Lange Antwort: UMFPACK ist eine der besten Bibliotheken zum Lösen von Ax = b, wenn A groß und dünn ist.
Unten ist ein Beispielcode (basierend auf umfpack_simple.c ), der eine einfache A und b generiert.
und löst Ax = b .
Die Funktion generate_sparse_matrix_problem erstellt die Matrix A und die
rechte Seite b . Die Matrix wird zuerst in Triplettform konstruiert. Das
Vektoren Ti, Tj und Tx beschreiben vollständig A. Die Tripletform ist einfach zu erstellen
Effiziente Sparse-Matrix-Methoden erfordern das Komprimierte Sparse-Spaltenformat. Umwandlung
wird mit umfpack_di_triplet_to_col ausgeführt.
Eine symbolische Faktorisierung wird mit umfpack_di_symbolic durchgeführt. Ein spärlicher
Die LU-Zerlegung von A erfolgt mit umfpack_di_numeric .
Die unteren und oberen Dreieckslösungen werden mit umfpack_di_solve durchgeführt.
Mit n als 500.000 auf meinem Computer dauert das gesamte Programm etwa eine Sekunde.
Valgrind berichtet, dass 369.239.649 Bytes (knapp über 352 MB) zugewiesen wurden.
Beachten Sie, dass Seite die Unterstützung von Boost für dünn besetzte Matrizen diskutiert in Triplett (Koordinate)
und komprimiertes Format. Wenn Sie möchten, können Sie Routinen schreiben, um diese Boost-Objekte zu konvertieren
zu den einfachen Arrays UMFPACK benötigt als Eingabe.
Angenommen, Ihre riesigen Matrizen sind spärlich, und ich hoffe, sie sind in dieser Größe, sehen Sie sich das PARDISO -Projekt an ein spärlicher linearer Löser, das ist, was Sie brauchen, wenn Sie Matrizen so groß wie Sie gesagt haben wollen. Ermöglicht eine effiziente Speicherung von nur Werten ungleich Null und ist viel schneller als das Lösen des gleichen Systems von dichten Matrizen.
Ich nehme an, dass deine Matrix dicht ist. Wenn es spärlich ist, können Sie zahlreiche spezialisierte Algorithmen finden, wie bereits von DeusAduro und Duffymo .
Wenn Ihnen kein (ausreichend großes) Cluster zur Verfügung steht, sollten Sie sich die Algorithmen ansehen, die außerhalb des Kerns liegen. ScaLAPACK hat einige Out-of-Core-Solver als Teil seiner Prototyp-Paket , siehe die Dokumentation hier und Google für weitere Informationen. Wenn Sie im Internet nach "Out-of-Core LU / (Matrix) Solvern / Paketen" suchen, erhalten Sie Links zu einer Fülle weiterer Algorithmen und Tools. Ich bin kein Experte für diese.
Für dieses Problem würden die meisten Leute jedoch einen Cluster verwenden. Das Paket, das Sie auf fast jedem Cluster finden, ist wieder ScaLAPACK. Darüber hinaus gibt es in der Regel zahlreiche andere Pakete auf dem typischen Cluster, so dass Sie auswählen können, was zu Ihrem Problem passt (Beispiele hier) und hier ).
Bevor Sie mit der Programmierung beginnen, möchten Sie wahrscheinlich schnell überprüfen, wie lange es dauern wird, um Ihr Problem zu lösen. Ein typischer Löser benötigt etwa 0 (3 * N ^ 3) Flops (N ist die Dimension der Matrix). Wenn N = 100000 ist, betrachten Sie daher 3000000 Gflops. Unter der Annahme, dass Ihr In-Memory-Solver 10 Gflops / s pro Kern leistet, betrachten Sie 3 1/2 Tage auf einem einzelnen Kern. Wenn die Algorithmen gut skalieren, sollte die Erhöhung der Anzahl der Kerne die Zeit in der Nähe von linear reduzieren. Hinzu kommt die I / O.
Sie sind sich nicht sicher über C ++ - Implementierungen, aber es gibt mehrere Dinge, die Sie tun können, wenn der Speicher abhängig vom Typ der Matrix, mit der Sie es zu tun haben, ein Problem ist:
Werfen Sie einen Blick auf die Liste frei verfügbarer Software zur Lösung von Problemen der linearen Algebra , zusammengestellt von Jack Dongarra und Hatem Ltaief.
Ich denke, dass Sie für die Problemgröße, die Sie betrachten, wahrscheinlich einen iterativen Algorithmus benötigen. Wenn Sie die Matrix A nicht in einem Sparse-Format speichern möchten, können Sie eine matrixfreie Implementierung verwenden. Iterative Algorithmen müssen typischerweise nicht auf einzelne Einträge der Matrix A zugreifen, sie müssen nur Matrix-Vektor-Produkte Av (und manchmal A ^ T v, das Produkt der transponierten Matrix mit dem Vektor) berechnen. Wenn die Bibliothek also gut gestaltet ist, sollte es ausreichen, wenn Sie eine Klasse übergeben, die Matrix-Vektor-Produkte zu beherrschen weiß.
Wie die angenommene Antwort nahelegt, gibt es UMFPACK. Aber wenn Sie BOOST verwenden, können Sie immer noch die kompakten Matrizen in BOOST verwenden und UMFPACK verwenden, um das System zu lösen. Es gibt eine Bindung, die es einfach macht:
Es ist ungefähr zwei Jahre datiert, aber es ist nur eine Bindung (zusammen mit ein paar anderen).
Siehe verwandte Frage: UMFPACK und BOOSTs uBLAS Sparse Matrix
Tags und Links c++ boost linear-algebra lapack umfpack