Hello everyone,

I'am using Eigen as the library for some calculations with large and sparse matrix and vectors, where the core job is to solve the linear system of equations:

A * x = b (1)

with A a square, sparse, unsymmetric matrix.

I've succeeded in using Eigen to solve a linear system of equations with a symmetric and sparse A. So, here, thanks all the authors of Eigen.

But, it seems that the solvers built-in all need a sysmmetric matrix A and that perhaps with the external library SuperLU one could solve such an equation like (1). If so, can anyone offer one or some examples for solving (1)? And is there any guide that presents how to install and use SuperLU with Eigen?


Thanks a lot

Zhe LI

Hi, you can either use the BICGSTAB<> iterative solver with the default Jacobi or, probably better, the IncompleteLUT<> preconditioner. The devel branch includes a SparseLU solver inspired from SuperLU but but zero dependency and even faster performance.

ggael wrote:Hi, you can either use the BICGSTAB<> iterative solver with the default Jacobi or, probably better, the IncompleteLUT<> preconditioner. The devel branch includes a SparseLU solver inspired from SuperLU but but zero dependency and even faster performance.

Hello team,
may I add a follow up to this question and the comment on preconditioning for BICGSTAB? How, exactly, would one precondition the sparse matrix?
Is it fair to ask if someone could illustrate this for the code pasted below? Here, the sparse Matrix SparseA, intentionally kept simple and small, is 4x4 dimensioned. The solution vector b is given such that a solution x for A*x=b will be x=(1,2,3,4). BICGSTAB will easily solve this.
However, would it be possible to implement preconditioning into this simple code, just to illustrate the required steps and syntax?

Code: Select all


#include <D:\VC2005\Eigen\LU>// provides LU decomposition
#include <D:\VC2005\Eigen\Array>// provides random matrices

#include "stdafx.h"
#include <iostream>
#include <vector>
#include <D:\VC2005\Eigen\Dense>
#include <D:\VC2005\Eigen\Sparse>
#include <ctime>
using namespace std;
using namespace Eigen;
using Eigen::MatrixXd;
typedef Eigen::SparseMatrix<double> SpMat; // declares a column-major sparse matrix type of double
typedef Eigen::Triplet<double> T;

// * * * * * * * * * * * * * * * * * * * * * * * * * * *
int main()
// * * * * * * * * * * * * * * * * * * * * * * * * * * *
{
//-----------------------------------------------------------
// We will considr a linear system with 4 variables and 4 equations
//-----------------------------------------------------------
   int n = 4; // size of the image
   std::vector<T> coefficients; // list of non-zeros coefficients
   Eigen::VectorXd b(n); // the rbight hand side-vector
//-----------------------------------------------------------
// Create dense Matrix A
//-----------------------------------------------------------
   MatrixXd MtxA = MatrixXd::Random(4, 4);
   MtxA << 5,2,4,0, 1,5,0,4, 3,0,5,2, 0,3,1,5;

//------------------------------------------
// Now create/populate an according Sparse Matrix
//------------------------------------------
   for(int j=0; j<n; ++j)
   {
      for(int i=0; i<n; ++i)
      {
           if (MtxA(j,i)*MtxA(j,i)>1e-5) {
         coefficients.push_back(T(j,i,MtxA(j,i)));
         } //if..
      } // ..for i
   }
   SpMat SparseA(n,n);
   SparseA.setFromTriplets(coefficients.begin(), coefficients.end());
//------------------------------------------
   cout << "Here is the Sparse matrix SparseA:\n"<< SparseA << endl;
   cin.get();   
   
//-----------------------------------------------------------
//   A right hand side vector (RHS Vector)is declared. I chose these values such that
//   the solution will be x1=1 x2=2 x3=3 x4=4
//-----------------------------------------------------------
   b << 21,27,26,29; //this is RHS for solution vector x=1,2,3,4
   //b << 26,28,29,29; //this is RHS for solution vector x=2,2,3,4
   cout << "Here is the Vector b:\n"<< b << endl;
   
//*************************************************************
//   SOLVE: Now we solve the NLS for RHS b using BICGSTAB  method.
//*************************************************************
   Eigen::BiCGSTAB<SparseMatrix<double> >  BCGST;
   BCGST.compute(SparseA);
   Eigen::VectorXd X;
   X = BCGST.solve(b);
   cout << "Here is the Vector x using BICGSTAB :\n"<< X << endl;
   cin.get();
}


Thanks,
Nicolo

That's easy:

Code: Select all

Eigen::BiCGSTAB<SparseMatrix<double>, Eigen::IncompleteLUT<double> >  BCGST;


that's all! Then you can access and tweak the ILUT before calling .compute(), e.g.:

BCGST.preconditioner().setDroptol(0.001);
BCGST.compute(SparseA);

See also: http://eigen.tuxfamily.org/dox-devel/cl ... teLUT.html

Finally, let me add that in the devel branch (soon 3.2) there is a SparseLU direct solver. Might be worth giving it a try!

Thanks a lot. I implemented your code snippet and it works. This will also help me understand more about the inner mechanics of the Eigen library.
I will need to get a feeling for the DropTolerance though. My first impression is that the conditioning works well for relatively small (hardly sparse at all) test systems where I have up to 10 variables to solve. Here, of course, BICGSTAB is an overkill and standard Gauss approach would solve that easily, too.
However, as I move to my real work the results get bad. Depending on segmentation (my work relates to finite element problems) I need to solve for 100 , 1000 and even 10000 variables and here is where results get bad. Let me mention that my matrix A is a Jacobi matrix and my vector x is a correction vector dx. I am solving the dx vector in the context of a Newton Raphson iteration scheme. Using either self written Gauss algorithm or Eigen's LUFactorization, though slow, the Newton Raphson iterations converge. Using Eigen's BICGSTAB with conditioning they only converge for extremely small problem statements (low segmentation, low number of cells, low resolution). Since my numerical skills hardly go beyond the Gauss/LUFactorization I am somewhat helpless on how to make best use of BICGSTAB and its options (such as conditioning). What other settings may I need to learn/understand to properly use that feature?
Another approach may be to identify submatrices and solve those sequentially to derive the entire matrix solution, but before going that way I like to find out about the capabilities of sparse solvers.
Just in case you may have additional guidance, be it links or existing samples, that would be appreciated. I may be asking for too much and I already appreciate your earlier response.
The direct SparseLU solver sounds promising and I will certainly have a take on it once v3.2 gets released.

In addition to my post above let me also provide an example that somehow plays into this. The following equation system Ax=b will (a) fail to solve using preconditioning and it will (b) solve using Bicgstab without preconditioning. Declare A and b as follows:

Code: Select all

MtxA << 1,2,0,1,  3,4,0,0,  0,0,2,3,  1,0,4,5;
 b << 9,11,18,33 ; //this is RHS for solution vector x=1,2,3,4
 

The solution to this sparse system is x=(1,2,3,4). Only case b) will solve this.

Code: Select all

//-------------------------------------------------------
//   a) solve using preconditioned BICGSTAB method
//-------------------------------------------------------
   Eigen::BiCGSTAB<SparseMatrix<double> , Eigen::IncompleteLUT<double>>  BCGSTCOND;
   BCGSTCOND.preconditioner().setDroptol(.001);
   BCGSTCOND.compute(SparseA);
   Eigen::VectorXd X;
   X = BCGSTCOND.solve(b);
   cout << "Here is the Vector x using BICGSTAB Preconditioned:\n"<< X << endl;
   cin.get();

//-------------------------------------------------------
//   b) solve using non-preconditioned BICGSTAB method
//-------------------------------------------------------
   Eigen::BiCGSTAB<SparseMatrix<double> >  BCGST;
   BCGST.compute(SparseA);
   Eigen::VectorXd Xx;
   Xx = BCGST.solve(b);   
   cout << "Here is the Vector x using BICGSTAB :\n"<< Xx << endl;
   cin.get();


Understanding why case a) fails may be helpful in dealing with Bicgstab. Is this system "playing against the Bigcstab rules"?

(As an afterthought, if we introduce a dependency into 3rd row for variable 3 and then adjust the solution vector b accordingly, that is:

Code: Select all

MtxA << 1,2,0,1,  3,4,1,0,  0,0,2,3,  1,0,4,5;
 b << 9,14,18,33 ;
 

then both cases a) and b) will make it. It seems fair to have a properly conditioned equation system as a requirement.)

Please, try the devel branch: http://bitbucket.org/eigen/eigen/get/default.tar.bz2. It should become the 3.2 stable release next monday and it contains some numerical fixes in the BiCGSTAB.

The default preconditioner of "BiCGSTAB<SparseMatrix<double> >" is a Jacobi preconditioner (inverse of the diagonal). Usually, ILUT helps the convergence of BICGSTAB, but bad settings might also kill the convergence. If I remember correctly, it is usually better to set a large Fillfactor (100) but a quite large DropTol (e.g., 1e-4). But if the default precondioner works you, then go for it!!

Also, since you are running a iterative non-linear solver, I suggest you the solveWithGuess function:

x = bicg.solveWithGuess(b, x0);

I must admit I did not really progress much with the SparseLU solver but I may simply not be using it correctly.

On another note I took my own initiative and considering that the equation systems that I need to solve refers to a finite element problem. The resulting equation system for such kind of problems can be formulated as a Tridiagonal Block Matrix, that is: the square coefficient matrix can be described as a tridiagonal matrix where each matrix element itself represents a submatrix. Of course, this is a special kind of a sparse matrix.

Solving standard tridiagonal matrices can be achieved easily (and without matrix operations) by means of the Thomas Algorithm. Solving Tridiagonal Block Matrices can be achieved by means of a modified Thomas algorithm (now with matrix operations) where divisions are carried out through matrix inverse operations. The algorithm is extremely fast compared to LU factorization since it only processes the tridiagonal elements and completely ignores all off-tridiagonal elements. Unless I oversaw something here, my understanding is that Eigen currently does not have the Block Thomas Algorithm. It may be a nice addition to the existing library. If you are interested I can forward my code to the Eigen development team. Just let me know and we can take this offline.

Regards!