# Matlab preconditioned conjugate gradient on big matrix

I have a sparse $5\,656\,236 * 5\,656\,236$ matrix $A$ with $166\,526\,888$ non-zero elements. The matrix comes from using the finite element method on a linear elasticity problem and is positive semi-definite. I'm trying to use the preconditioned conjugate gradient method to solve it, particularly, the pcg() function in MATLAB.

Can I expect this to converge at all? I tried

L = ichol(A);
[u,flag,res,iter,resvec] = pcg(A,F,1e-6,max_ter,L,L');


But it does not seem to converge at all. The residual starts by increasing a lot and after 2000 iterations (which takes some hours) the relative residual is up at $10^4$.

So I guess my question is, how would you attack a system of this size? I'm sure there are a lot of ways to lure the method into converging faster. Any hints and tips are welcome!

• What is $A$? Where does it come from? Do you have to use MATLAB? Do you have to use CG? – Geoff Oxberry Nov 18 '13 at 6:25
• Sorry, I should have specified that the system comes from using linear finite elments on a linear elasticity problem. I do not have to use MATLAB or CG, no. – burk Nov 18 '13 at 8:13
• What do you get when you simply do A\F? – Rody Oldenhuis Nov 19 '13 at 10:22
• When doing A\F, MATLABs memory consumption passes 60GB, and it did not converge overnight. From what I heard, this is not the way to do it with matrices of this size, what are your experiences? – burk Nov 19 '13 at 10:31

The MATLAB function ichol computes by default the zero fill-in variant of the incomplete Cholesky factorisation. You could try to allow more fill-in (and thus try to improve the preconditioner) by using a drop tolerance, e.g.:

L = ichol(A,struct('type','ict','droptol',1e-03,'michol','on'));


Also note, that single level preconditioners such as ILU (or ICHOL), sparse inverses, simple splittings etc. do not scale well for large problem sizes. So if the above won't work, you could try to use algebraic multigrid (AMG). There are AMG packages which are well-suited for solving linear elasticity, e.g., BoomerAMG from Hypre or ML in Trilinos (as far as I know, ML has also a Matlab interface). If you don't like C/C++, you might consider trying PyAMG which is in Python (build upon NumPy and SciPy).

• PETSc has a MATLAB interface for serial problems, and interfaces with hypre (and thus, BoomerAMG) and ML. I'm not sure if you can call a hypre solver or ML from MATLAB, but it's worth looking into. – Geoff Oxberry Nov 18 '13 at 19:20
• @Geoff Oxberry: Thanks for the reminder, I've totally forgot about PetSc :) I don't think Hypre has a MATLAB interface but I think ML has something like that (at least it seems from the ML user's guide but I have no experience with that). Also, ML has at least a Python interface via PyTrilinos. – Algebraic Pavel Nov 18 '13 at 22:31
• I mean, you may or may not be able to call hypre and ML via PETSc's MATLAB interface. In any case, PETSc also has a Python interface, which can certainly be used to call hypre or ML. – Geoff Oxberry Nov 18 '13 at 22:58
• OK guys, thanks for the help. I did not manage to solve the original system before the project was due, but I will accept this answer. MATLAB can only get you so far, but I did not have time to try any of the other solvers you suggested. – burk Nov 27 '13 at 13:10

If your residual increases, then you probably have a matrix that is either not symmetric or at least not positive definite. In those cases, CG will not work, with or without preconditioner.

There are two ways to find out:

• Try a smaller system that results from the same program but using a much coarser mesh and see if CG converges for that.
• Inspect your code to see if it really does what it's supposed to do when assembling the matrix.
• I actually though about this. I have checked that the matrix is both symmetric and positive definite for smaller problems (using chol). Weird thing is that I just solved a 3500718 x 3500718 system with 111551450 non-zero elements, and it took about 5 hours. The original system has been running for over a day. I always seem to get a bit of a jump in the residual at the beginning of the method. – burk Nov 19 '13 at 7:36
• @burk That's possible; CG minimises the $A^{-1}$-norm of the residual instead of the Euclidean one. – Algebraic Pavel Nov 20 '13 at 16:25
• Ahh, of course, I knew that.. Thanks for pointing it out :-) – burk Nov 20 '13 at 20:37