Résolution d'un système clairsemé: Eigen vs MATLAB

Question

J'ai un système linéaire clairsemé Ax = b. Dans mon application, Aest une matrice clairsemée symétrique avec une taille typique d'environ 2,500,000 x 2,500,000, avec des non-zéros sur la diagonale principale et sur une autre diagonale (plus le symétrique de celui-ci). Cela en fait 2-3 non-zéros par ligne / colonne. 

Pour tester mon code, je compare MATLAB et Eigen. J'ai créé une matrice clairsemée de 1 000 000 x 1 000 000 A. Dans MATLAB, je l'utilise simplement x = A\bet cela prend environ 8 secondes. Dans Eigen, j'ai essayé plusieurs solveurs. SuperLU prend environ 150 s. SimplicialCholesky prend environ 300 secondes. UmfPackLU prend environ 490 s. Ces temps sont trop longs pour moi; sur des données réelles, cela prend trop de temps pour être utile. D'autres solveurs donnent des résultats complètement différents par rapport à MATLAB, les solveurs itératifs ont pris trop de temps. SimplicialCholesky, SuperLU et UmfPackLU donnent des résultats similaires (ils diffèrent aux décimales), donc j'espère que cela compte comme la même chose. Code propre: 

// 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";

Peut-être devrais-je mentionner que le calcul fonctionne sur les 8 cœurs, donc quand je regarde l'heure, j'obtiens 8 fois, ce que je résume. De plus, le calcul est (jusqu'à présent) encapsulé dans la bibliothèque .dll .cu, il sera parallélisé via CUDA à l'étape suivante. J'ai mesuré les temps pour toutes les méthodes séparément afin d'éviter certains chevauchements de comptage.

J'ai trouvé les solutions possibles suivantes pour accélérer le calcul:

• Utilisez lu normal , ne fonctionne pas pour un système clairsemé;  
• En lien avec la bibliothèque BLAS / LAPACK , je pense avoir fait cela.
• essayez différents solveurs ,  ou wrappers , d'autres solveurs n'ont pas donné les mêmes résultats que MATLAB; les réponses ici étaient trop spécifiques à chaque cas;
• multithreading, utilisez le compilateur avec les optimisations activées effectuées (compilateur - optimisations maximales, vitesse de faveur), toujours très lent;
• utilisez UmfPack, comme le fait MATLAB, pour obtenir des performances similaires - c'est encore plus lent que SimlicialCholesky
• liste d'autres bibliothèques possibles travaillant avec des matrices , mais je ne sais pas comment elles traiteraient mon cas

Est-ce que je peux faire quelque chose pour accélérer les calculs en utilisant Eigen, donc cela prend le même temps que MATLAB? Est-ce que j'utilise le bon solveur, en ce qui concerne la taille et la rareté de la matrice? Est-ce que j'utilise correctement les solveurs actuels? Dois-je effectuer une configuration supplémentaire, inclure d'autres bibliothèques? Si ce n'est pas possible, y a-t-il d'autres bibliothèques que je pourrais utiliser? 

Je travaille sur Windows 10, machine 64 bits. J'ai Visual Studio 2019. 

Answer 1

J'ai récemment essayé de nombreux solveurs linéaires pour mon solveur de collocation spectrale, et j'ai trouvé que le "tatou" est le plus rapide qui résout Ax = b dense basé sur la bibliothèque openblas. Eigen3.3 est très lent même avec "setNumbthreads", je ne trouve toujours pas la raison. Si vous voulez le résoudre avec Cuda ou OpenMP. Je vous suggère fortement d'utiliser la bibliothèque de paralution. cela fonctionne bien pour mon problème. Salutations

http://www.paralution.com/