# Solve a very large linear system (question about a library linear algebra to do this)

I need to solve a very large linear system (coming from finite element method). I'm currently using the Intel MKL library, but the system has been delayed more than 20 hours.

The matrix of the system is sparse and it is written in Morse (CRS)

I'm programming in Fortran 90, but I guess the language of the library is not important.

• Large of the matrix: 5636087x5636087 or bigger.
• Type: sparse (CRS), with real numbers, non symmetric.
• Number of nonzero elements of the matrix: 202516857
• The matrix is non-symmetric, and I suppose that is well conditioned (coming from finite element method).
• Can you give an order of magnitude for how large the matrices are? Also is there any structure to the matrices? Are the symmetric positive definite, are the matrices well conditioned, are they real or complex etc? Also when you say invert the matrix do you mean you need the actually inverse or just need to solve a linear system. The more info you can give the better. Commented Jan 21, 2016 at 5:43
• Which routine from the MKL are you using? I don't know what kind of machine you have (especially how much memory), but I wouldn't even try to touch such a matrix with a direct solver... Commented Jan 21, 2016 at 14:12
• I'm using pardiso on MKL.My machine has 256GB of RAM with 96 processors Intel(R) Xeon(R) CPU E7-4860 v2 @ 2.60GHz Commented Jan 21, 2016 at 14:27
• Are you using the out of core option in MKL Pardiso? Commented Jan 21, 2016 at 14:28

Since the matrix is considerably sparse and well-conditioned (if it is true), I suggest you can try to use Krylov subspace method, which can only use the information of matrix-vector product; such as GMRES(m), CGS, BiCGSTAB, and IDR(s). With aspect of memory and well-conditioned matrix, maybe BiCGSTAB, CGS and IDR(s) are preferable. They are many FORTRAN codes of these above four solvers online. Or even, you can try to use the AGMG method, which have a FORTRAN code in website http://homepages.ulb.ac.be/~ynotay/AGMG/ and other method is http://www.hsl.rl.ac.uk/catalogue/hsl_mi20.html

LSQR uses only 4 N-vectors of MEMORY. It might be your choice if storage must be minimized. It is fairly good otherwise, too, but probably not the best one (unless your matrix is non-square; then the ones below do not work bot LSQR does). LSQR is CG for $$A^*A$$ but numerically more robust etc. 1DOTs/MatVec, 3AXPYs/MatVec (2 MatVec). https://web.stanford.edu/group/SOL/software/lsqr/lsqr-toms82a.pdf

Slide 3 of [1] suggests GMRES only if you have a good preconditioner (and lots of memory!). Then GMRES usually converges quickly in terms of MatVec (Ax matrix-vector products) (not in terms of DOTs and AXPYs) but uses m+3 MEMORY (N-vectors) if you restart at m.

If the matrix is NOT strongly indefinite (re $$x^*Ax$$ is usually roughly $$\ge0$$ or usually roughly $$\le0$$, where "usually" means "for most random x"), use BiCGSTAB (or BiCGSTAB(l) if its eigenvalues have large imaginary parts).

Else use IDR (or IDRstab if your matrix is real and its eigenvalues have large imaginary parts - I guess that there is no similar quick test to this?).

Of course, a good preconditioner is recommendable in any case.

IDR(s): the best code is IDR(s)BiOrtho, i.e., "Algorithm 913". Try s=4. Increase s, if it does not converge quickly enough. Decrease s if it does. Or start with s=1 if you lack memory.

IDR(s)BiO = $$s+\frac2{s+1}$$ DOTs/MatVec, $$2s+\frac2{s+1}$$ AXPYs/MatVec, $$3s+4$$ MEMORY (N-vectors), by page 5:11 of [https://dl.acm.org/doi/10.1145/2049662.2049667] That also contains the Matlab code (and the algorithm). More can be found here: [https://web.archive.org/web/20190522043534/http://ta.twi.tudelft.nl/nw/users/gijzen/IDR.html]

IDRstab (IDR(s)stab(L)): The best current Matlab code (and algorithm) for is probably that by Martin Neuenhofen (according to van Gijzen it is better than his and the best that he knows). In the experiments by Neuenhofen, $$L\le2$$ is always enough.

IDR(s)stab(2) = $$s-1+\frac9{2(s+1)}$$ DOTs/MatVec, $$\frac52s+2$$ AXPYs/MatVec, $$5(s+1)-1$$ MEMORY (N-vectors), by pages 35-37: https://arxiv.org/pdf/1708.01926.pdf

Matlab code: http://www.martinneuenhofen.de/GMstab/GMstab.html It also contains Mstab which is IDRstab for the case where you solve for multiple $$b$$ with the same $$A$$.

All these codes work for general complex A and b. The problems with IDR(s) (vs. IDRstab) are probably not likely unless A is real but its eigenvalues are far from real.

For random matrices, GMRES tends to converge linearly (i.e., all too slow) and the others even more slowly. Fortunately, you rarely meet random matrices, even if not preconditioned.