Hello,
I am trying to use the BiCGSTAB solver to solve a linear problem but it is not working.
The help on your website says:
VectorXd x(n), b(n);
SparseMatrix<double> A(n,n);
BiCGSTAB<SparseMatrix<double> > solver;
solver(A);
x = solver.solve(b);
The line solver(A) gives me an error when trying to launch the program. I gave a look at other solvers and noticed that the line is usually solver.compute(A). I tried but then the result is wrong.
I also would like to know if there is a way to check if a Sparse matrix is positive definite without converting it to a dense matrix.
Thank you for the help

Make sure you are using latest Eigen version: https://bitbucket.org/eigen/eigen/get/3.2.tar.gz
Do you have a self-contained example?
How did you check the result?
What about solver.info(), solver.iterations(), solver.error()?

I solved my system of equations using the SparseLU solver and then the BiCGSTAB solver.
My system is a FEM problem. When I use the SparseLU solver, the result seems correct because it does verify my boundary conditions and the values are physically possible. The BiCGSTAB solver gives totally wrong results not even good on imposed boundary conditions. I tried on a very simple system of equations and the BiCGSTAB solver is working but the correct way is solver.compute(A). If I use solver(A), it cannot work because solver is not a function nor a vector.
I checked solver.info() and solver.iterations(). It gives me 2 and something around 1500.
Because of imposed boundary conditions, some entries of the matrix are 1e100 which is far bigger than other normal entries. Could this be the problem ?

Cannot you impose your boundary condition by moving the fixed variable to the right hand side? Moreover, by default BiCGSTAB use a Jacobi preconditioner which makes sense only if the matrix is diagonally dominant. Depending on how you built your matrix, the huge coefficients might not be on the diagonal which is a big problem for the Jacobi preconditioner. Finally, when using BiCGSTAB, it is very often required to use a more advanced preconditioner like ILUT:

BiCGSTAB<SparseMatrix<double>, IncompleteLUT<double> >

My big values are only on the diagonal and nowhere else.
I did more tests and I have to admit I am a bit lost. I have 498 nodes, so it makes 498*3 = 1494 degrees of freedom. That will be the size of my stiffness matrix.
I tried to build my stiffness matrix without applying any boundary. This way I just have a sparse matrix with values of the same range, no more big values.
Using this matrix, I tried to solve a system KU = F, K is the stiffness, U the unknown and F the force vector. To do my test, I took every component of F to 100. Since my matrix is not too big, I can convert it to a dense matrix and make sure the determinant is not 0.
To find U, I used 3 methods:
- Sparse LU solver
- BiCGSTAB solver
- convert K to a dense matrix and use K.inverse()

For those 3 methods, I computed the norm of the residual, KU - F. Here is the part I do not understand: no matter what method I used, the norm is still very big.
For sparse LU the norm is 1e5, for BICGSTAB 1e111 (completly wrong it seems) and 2e6 for the last method.
Even if Sparse LU is the best, it is still far from being good. I also tried BiCGSTAB<SparseMatrix<double>, IncompleteLUT<double> > and the result gives me #IND, so it does not work.
Any idea where I might have done something wrong ?

"K.inverse()" ouch! never ever inverse a matrix unless you know it's closed to unitary. Also, don't forget to normalize the error measure by the norm of f. If you want a reference solution do:

Code: Select all

MatrixXd A = ...;
x = A.colPivHouseholderQr().solve(f);
cout << (A*x-f).norm()/b.norm() << endl;


Then, it might also be interesting to look at the condition number of A:

Code: Select all

JacobiSVD<MatrixXd> svd(A); // will take times
cout << svd.singularValues()(A.cols()-1)/svd.singularValues()(0) << endl;


Finally, note that the BiCGSTAB algorithm is not guaranteed. You can also export your matrix so that we can have a look to see if there is anything wrong on our side. See http://eigen.tuxfamily.org/dox-devel/gr ... tml#title3 for a code snippet.

I tried to use the preconditioner IncompleteLUT and my system of 1494 unknowns can be solved using BiCGSTAB.
I also imposed my boundary conditions by moving them to the right hand side, therefore my matrix had no more big numbers. But the normal BiCGSTAB was still not working.
Finally, I tried the preconditioner with bigger problems, and it stops instantly, probably because the sparse matrix is too big. It was around 100 000 unknowns. For such problem, it is not possible to use the SparseLU solver, it takes too much time.
Is there anyway to use the preconditioner with very big ( and very sparse) problem ?

What do you mean by stops instantly? memory overflow? If so, you can play with solver.setFillfactor(N) before calling solver.compute(). E.g., try with small N (e.g., in the [2:10] range).

I tried to use setFillfactor, but it does not exist for any type of solver I tried. Did I forget to add something in my code to make it available ?

sorry, that's a method of the preconditioner:

solver.preconditioner().setFillfactor(N);