2

I want to solve the Ax=b equation where A is a sparse matrix (1,964,568 x 1,964,568 nnz=75256446) and b is also sparse (1,964,568 x 1,964,568 nnz= 25354926) by using Eigen library in C++.

At first I was trying to use the Eigen Sparse LU to solve my problem, after a few hours I ran out of memory (I have 128GB RAM). After this I included the INTEL MKL library with Pardiso solver. Even with this I cant solve my problem. Maybe anyone has some tips to solve my problem?

#define EIGEN_USE_MKL_ALL

#include <Eigen/Sparse>
#include <unsupported/Eigen/SparseExtra>
#include <iostream>
#include <Eigen/OrderingMethods>
#include <Eigen/PardisoSupport>


typedef Eigen::SparseMatrix<double>SpMat; 
typedef Eigen::COLAMDOrdering<int>Order;

int main()
{

    SpMat A;
    SpMat B;
    SpMat X;

    Eigen::loadMarket(A, "MatK.mtx");
    Eigen::loadMarket(B, "MatM.txt");

    A.makeCompressed();
    B.makeCompressed();

    Eigen::PardisoLU<SpMat>solver;

    solver.analyzePattern(A);

    solver.factorize(A);

    X = solver.solve(B);


}

I can compile my code and run it. I just need better performance and less memory.

Vex994
  • 49
  • 1
  • 6
  • 2
    PARDISO is a _direct_ solver and likely requires a large amout of additional memory. Try to use some _iterative_ solver. These might also require some memory, e.g., for preconditioning and storing Krylov subspace vectors, however typically not as much as direct solvers. Both Eigen and Intel MKL include some iterative solvers. BTW, do you use 64-bit binary and libraries? – Daniel Langr Aug 29 '19 at 07:59
  • So you are trying to solve `A * X = B` for `X` where both `A` and `B` are huge, sparse matrices? And you assume that the solution `X` is also a sparse matrix? The usual situation is that `x` and `b` are vectors. Not sure if it is possible to obtain a sparse solution matrix `X` with usual solvers. I would probably rather solve for each column `b_i` of `B` and assemble the results `x_i` to form a matrix `X`. – RHertel Aug 29 '19 at 12:52
  • Yes, I noticed that I have to slice my B matrix. The Problem I want so solve is: X=inv(A)*B. Since I dont want to calculate the inverse I wrote it as AX=B and solve this. – Vex994 Aug 30 '19 at 10:31
  • The amount of fill-in (and therefore the amount of memory required) strongly depends on the structure of your matrices (and the ordering method, which it seems can't be configured for `Eigen::ParadisoLU`). Does any of your matrices have some special structure which could be exploited? For a definitive answer you should not only provide the source, but all data required to reproduce the problem (see [mre]). – chtz Aug 30 '19 at 14:14
  • I don't see any reason to assume that `X` is sparse, and even if it is, I doubt that a solver would return a `SparseMatrix`. – RHertel Aug 30 '19 at 17:32
  • I tried the BiCGSTAB now. I am seperating the matrix B and start solving. Afterwards I assemble the result matrix X. It is working like this but each Iteration (seperating, solving, assembly) takes 2-3s and this about 2 million times. Any suggestion to speed it up? – Vex994 Sep 04 '19 at 08:34
  • @Vex Yes. You can speed it up by using a SparseMatrix with RowMajor ordering (instead of the default ColMajor). Eigen has a multi-threaded version of BiCGSTAB, but it requires a SparseMatrix with RowMajor ordering. Don't forget to compile with the required OpenMP flag (depends on your compiler). – RHertel Sep 04 '19 at 10:11
  • @RHertel: Thank you for your advice. It is running now at 93% CPU and each iteration takes around 0.25s. In my case that's still a long time to wait :) – Vex994 Sep 04 '19 at 12:25

1 Answers1

2

Since I haven't found any method to solve this. I tried to change my RHS matrix, because the huge amount of columns was a big problem. I saw in my algorithm that my RHS gets multiplied by another matrix, which reduces the columns to 10. With this it was possible to solve with Intel Pardiso LDLT solver instead of LU. The iterative solver took too many iterations to converge.

Vex994
  • 49
  • 1
  • 6