lizhe
Registered Member

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 builtin 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 
ggael
Moderator

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.

nicolo
Registered Member

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?
Thanks, Nicolo 
ggael
Moderator

That's easy:
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/doxdevel/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! 
nicolo
Registered Member

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. 
nicolo
Registered Member

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:
The solution to this sparse system is x=(1,2,3,4). Only case b) will solve this.
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:
then both cases a) and b) will make it. It seems fair to have a properly conditioned equation system as a requirement.) 
ggael
Moderator

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., 1e4). But if the default precondioner works you, then go for it!! Also, since you are running a iterative nonlinear solver, I suggest you the solveWithGuess function: x = bicg.solveWithGuess(b, x0); 
nicolo
Registered Member

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 offtridiagonal 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! 
Registered users: Baidu [Spider], Bing [Bot], Exabot [Bot], Google [Bot], Sogou [Bot], vlas, Yahoo [Bot]