Hello,
I'm testing some Finite Volume stuff with Eigen and I am having some weird problems with the sparse implementation of BiCGSTAB: Sometimes when I try to solve a linear system with a diagonal Matrix (I can't tell whether the matrix will be diagonal or even symmetric, therefore I don't use any more simple solver) it just gives me nans.

I reproduced the behavior in the following simple code. Using Eigen 3.1 (same with 3.1.1) for me this yields nans as result.
Furthermore strange things happen if I change some values:
-Changing b(0) to any other value yields the correct value
-Changing all diagonal elements to 1.7 will yield nans independent of b(0)

Code: Select all

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <Eigen/IterativeLinearSolvers>

using namespace Eigen;
using namespace std;

typedef SparseMatrix<double, RowMajor> SpMat;

int main(int argc, char *argv[])
{   
    BiCGSTAB <SpMat>  solver;
    VectorXd b=VectorXd::Constant(4,0.6);
    SpMat A(4,4);
    for(int i=0;i<4;i++)
        A.insert(i,i)=1.6;
    A.makeCompressed();
    solver.compute(A);
    cout << A << endl;
    cout << b << endl;
    cout << solver.solve(b)<<endl;
    return 0;
}

Hello,
BiCGSTAB uses the diagonal preconditioner by default.
Since you have a diagonal matrix, the method should converge in one iteration.
We didn't consider that case where the preconditioner is exact.
It is now modified in the devel branch.