Sparse-System lösen: Eigen vs. MATLAB

Question

Ich habe ein spärliches lineares System Ax = b. In meiner Anwendung Ahandelt es sich um eine symmetrische Matrix mit geringer Dichte und einer typischen Größe von etwa 2.500.000 x 2.500.000, mit Nicht-Nullen auf der Hauptdiagonale und auf einer anderen Diagonale (plus der Symmetrie zu dieser). Dies macht es 2-3 Nicht-Nullen pro Zeile / Spalte. 

Um meinen Code zu testen, vergleiche ich MATLAB und Eigen. Ich habe eine 1.000.000 x 1.000.000 spärliche Matrix erstellt A. In MATLAB benutze ich einfach x = A\bund es dauert ungefähr 8 Sekunden. In Eigen habe ich mehrere Löser ausprobiert. SuperLU dauert ca. 150 s. SimplicialCholesky dauert ungefähr 300 Sekunden. UmfPackLU dauert ca. 490 s. Diese Zeiten sind zu lang für mich; Bei realen Daten dauert es einfach zu lange, um nützlich zu sein. Andere Löser liefern völlig andere Ergebnisse als MATLAB, iterative Löser haben zu lange gedauert. SimplicialCholesky, SuperLU und UmfPackLU geben ähnliche Werte an (sie unterscheiden sich an den Dezimalstellen), daher hoffe ich, dass dies gleich zählt. Eigencode: 

// prepare sparse matrix A
    std::vector<T> tripletList; // I am leaving filling the triplet list out
    Eigen::SparseMatrix<float> A(k, k); // k is usually around 2500000, in the test case I described here it is 1000000
    A.setFromTriplets(tripletList.begin(), tripletList.end());
    A.makeCompressed();

// prepare vector b
    Eigen::Map<Eigen::VectorXf> b; // vector b is filled with values

// calculate A x = b and measure time - for SimplicialCholesky
    t1 = std::chrono::steady_clock::now();
    Eigen::SimplicialCholesky<Eigen::SparseMatrix<float>> solver_chol(A);
    x = solver_chol.solve(b);
    t2 = std::chrono::steady_clock::now();
    log_file << "SimlicialCholeskytime: t2 - t1 = " << std::chrono::duration_cast<std::chrono::seconds>(t2 - t1).count() << " s \n";

// calculate A x = b and measure time - for SparseLU
    t1 = std::chrono::steady_clock::now();
    Eigen::SparseLU<Eigen::SparseMatrix<float>> solver_slu(A);
    x = solver_slu.solve(b);
    t2 = std::chrono::steady_clock::now();
    log_file << "SparseLU time: t2 - t1 = " << std::chrono::duration_cast<std::chrono::seconds>(t2 - t1).count() << " s \n";

// calculate A x = b and measure time - for UmfPackLU - here I had to convert to double.
    Eigen::SparseMatrix<double> Ad = A.cast <double>();
    Ad.makeCompressed();
    Eigen::VectorXd bd = b.cast <double>();
    t1 = std::chrono::steady_clock::now();
    Eigen::UmfPackLU<Eigen::SparseMatrix<double>> solver(Ad);
    Eigen::VectorXd xd = solver.solve(bd);
    t2 = std::chrono::steady_clock::now();
    log_file << "UmfPackLU time: t2 - t1 = " << std::chrono::duration_cast<std::chrono::seconds>(t2 - t1).count() << " s \n";

Vielleicht sollte ich erwähnen, dass die Berechnung auf allen 8 Kernen ausgeführt wird. Wenn ich also die Zeit anschaue, erhalte ich 8 Mal, was ich zusammenfasse. Außerdem wird die Berechnung (bisher) in die DLL-Bibliothek .cu eingeschlossen und im nächsten Schritt über CUDA parallelisiert. Ich habe die Zeiten für alle Methoden separat gemessen, um Überschneidungen zu vermeiden.

Ich habe die folgenden möglichen Lösungen gefunden, um die Berechnung zu beschleunigen:

• Verwenden Sie normales lu , funktioniert nicht für spärliches System;  
• Ich glaube, ich habe dies getan, indem ich auf die BLAS / LAPACK-Bibliothek verlinkt habe.
• versuchen , verschiedene Solver ,  oder Wrapper , andere Löser nicht gleiche Ergebnisse wie MATLAB gegeben haben; Die Antworten hier waren zu fallspezifisch.
• Multithreading, Compiler mit aktivierten Optimierungen verwenden (Compiler - maximale Optimierungen, Geschwindigkeit bevorzugen), immer noch sehr langsam;
• Verwenden Sie UmfPack wie MATLAB, um eine ähnliche Leistung zu erzielen - es ist sogar langsamer als SimlicialCholesky
• Liste anderer möglicher Bibliotheken, die mit Matrizen arbeiten , aber ich weiß nicht, wie sie mit meinem Fall umgehen würden

Kann ich irgendetwas tun, um die Berechnungen mit Eigen zu beschleunigen, sodass es eine ähnliche Zeit wie bei MATLAB dauert? Benutze ich den richtigen Löser in Bezug auf die Größe und Sparsamkeit der Matrix? Benutze ich aktuelle Löser richtig? Muss ich zusätzliche Einstellungen vornehmen und einige andere Bibliotheken hinzufügen? Wenn es nicht möglich ist, gibt es einige andere Bibliotheken, die ich verwenden könnte? 

Ich arbeite an Windows 10, 64-Bit-Maschine. Ich habe Visual Studio 2019. 

Answer 1

Ich habe kürzlich viele lineare Löser für meinen spektralen Kollokationslöser ausprobiert und festgestellt, dass das "Gürteltier" das schnelle ist, das dichtes Ax = b basierend auf der Openblas-Bibliothek löst. Eigen3.3 ist selbst mit "setNumbthreads" sehr langsam, ich kann den Grund immer noch nicht finden. Wenn Sie es mit Cuda oder OpenMP lösen möchten. Ich empfehle Ihnen dringend, die Paralutionsbibliothek zu verwenden. es funktioniert gut für mein Problem. Grüße

http://www.paralution.com/