linear algebra - Calculation mistakes of SparseLU in Eigen (C++)

Question

Recently I found that in some certain circumstances, the SparseLU solver will get the wrong answer.

I give an example in the attachment, where "A.csv" is the coefficient matrix of the equation Ax=b, and "b.csv" is the right part of the equation.

The codes are like this:

// solve the linear system: Ax = b
Eigen::VectorXd sparselu_x, pardisolu_x;
SparseLU<spmat> solver_splu;
PardisoLU<spmat> solver_pdslu;

// SparseLU results
solver_splu.compute(A);
sparselu_x = solver_splu.solve(b);

// PardisoLU results
solver_pdslu.compute(A);
pardisolu_x = solver_pdslu.solve(b);


// print some of the difference
cout << "sparselu results - pardisolu results (last 5 items):

" << (sparselu_x - pardisolu_x).tail(5) << endl;

cout << "



A*x - b (sparselu results, last 5 items):

" << (A * sparselu_x - b).tail(5) << endl;
cout << "



A*x - b (pardisolu results, last 5 items):

" << (A * pardisolu_x - b).tail(5) << endl;

And the output:

sparselu results - pardisolu results (last 5 items):

 6.25516e-10
 6.36814e-08
-1.45986e-05
  -0.0170633
   -0.437249




A*x - b (sparselu results, last 5 items):

           0
 1.81899e-12
-1.81899e-12
     18.0459
     621.588




A*x - b (pardisolu results, last 5 items):

1.81899e-12
          0
          0
          0
          0

Full data and codes about reading the csv file and solving the system are uploaded at:

https://github.com/cpc-tyym/SparseLU-eg

Because of the efficency, I really prefer the sparselu solver, but this mistake makes me choose the pardiso solver, which is 3 times slower than the sparselu solver. So could anyone find out whether the sparselu solver is wrong in somewhere, or did I make any mistakes when solving the equation?

question from:https://stackoverflow.com/questions/65912200/calculation-mistakes-of-sparselu-in-eigen-c

Answer 1

Waitting for answers