Solving a large (n>1000000) linear system of equations

Question

I have reached my limit on the following problem:

As part of my FEA code (in MATLAB) I need to find x, x=A\b. Both A and b are sparse, complex, double precision matrix and vector respectively. The size of A is (n,n) and b is (n,1) where n is 850000 and can increase to up 2000000. In addition, A is symmetric and mostly diagonal.

I have two HPC clusters to my disposal, one with 5.7TB of RAM and the other with 1.5TB (but faster cores). Technically speaking I can solve the system as is, and just wait for approximately 15days. However, I need to perform solve the system of equations up to 10x per simulations. Therefore, 15days is not an acceptable amount of time.

I tried iterative methods however these did not yield the same results as the backslash operator. Also in some cases convergence was not obtained.

I have converted the x=A\b part into a mexa64-format to potentially reduce the time. However, I fear it will still take days.

Any suggestions on how to tackle this problem? Are there any commercial codes that can do this faster/more efficient? And how commercial FEA packages solve this problem when a model has over 1m nodes?

Any help in much appreciated! Thanks in advance.

Answer 1

If your matrix is Hermitian (Aij = Aji*) and positive definite, then you can't go wrong with a conjugate gradient approach. That is almost always the fastest method if you have a single right hand side and is guaranteed to converge. Additionally, this allows you to use a guess solution (for instance, from a previous time step) which can be beneficial. While the solution is different from a direct decomposition (how could it not be with a gazillion coefficients), its probably not that different.

If it's literally a symmetric complex matrix (Aij = Aji), then that is very weird but maybe still better than a general complex matrix. In either case, try using a routine that only requires you to pass half of the matrix (I don't think matlab supports this functionality). That will cut your memory usage in half and might make things a bit faster (moving 16 GB worth of numbers is slower than 8GB).

Finally, I can't find an answer on whether Matlab supports the Intel MKL Inspector-executor sparse BLAS routines. In my experience, using MKL in my bottlenecks always increases performance. It's closed source and you have no idea whats happening in there -- probably magic.

Answer 2

The sparse matrix computation is a big area itself. Only thing you can do is to get more specific in it and use a more customizable solution, as Matlab do not allow you to do so and see it's source code and see if you can do any improvement. As i know mainly there are two methods for solving a sparse linear system, Direct and Iterative. For each type many algorithms are there. For example Cholesky Decomposition have a very good performance and low memory usage, but only works on SPD (symmetric Positive Definite) matrices. other methods are LU, LDL, QR etc. Each one meet specific matrix types to solve.

Matlab also uses one of these methods behind the scene (not sure but I've read somewhere that it have used SuiteSparse code underneath the A/b operator). Based on type of your matrix (whether it is Symmetric Positive Definite or Non-symetric etc.) you need to choose one schema (i mean either Iterative or Direct). I think Direct methods suits better for such large system (because of rounding errors). Also there are GPU based solvers available for both direct and iterative methods and highly parallelized, but maybe not suitable for your big problem where matrix being factorized is more than 8GB itself.

Steps you should do:

• Identify the type of your matrix (SPD, Symmetric or general).
• Choose best Direct method for you matrix.
• Find a free or commercial package that can handle such factorization (examples are Eigen, Intel MKL, SuitSpare) and have optimized code also support parallel computation.
• Make your hand a little dirty, write a very small program that takes the matrix A and vector b from a file and find x using one of these packages with for example C++ code (they have ready to use example codes).

Also pls keep in mind that each method and package can change the ordering of matrix before starting computation and choosing a good ordering can have a big effect on speed or even accuracy.

This link should be also interesting:

Do not invert that matrix

Also two books on sparse matrices:

• Sparse Matrix Technology by Sergio Pissanetzky
• Direct Methods for Sparse Linear Systems by Timothy A. Davis