0

最近我发现在某些特定情况下，SparseLU 求解器会得到错误的答案。

我在附件中举了一个例子，其中“A.csv”是方程Ax=b的系数矩阵，“b.csv”是方程的右边部分。

代码是这样的：

// 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):\n\n" << (sparselu_x - pardisolu_x).tail(5) << endl;

cout << "\n\n\n\nA*x - b (sparselu results, last 5 items):\n\n" << (A * sparselu_x - b).tail(5) << endl;
cout << "\n\n\n\nA*x - b (pardisolu results, last 5 items):\n\n" << (A * pardisolu_x - b).tail(5) << endl;

和输出：

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

关于读取csv文件和解决系统的完整数据和代码上传到：

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

因为效率，我真的更喜欢sparselu求解器，但是这个错误让我选择了pardiso求解器，它比sparselu求解器慢3倍。那么任何人都可以找出 sparselu 求解器是否在某个地方出错，或者我在求解方程时是否犯了任何错误？

4

0 回答 0