# Trying to understand splitting-based iterative method for 2D Laplace problem

I am trying to understand the theory behind a splitting based iterative method which uses the incomplete Cholesky factorization. Before giving the specific details, let me first give the problem background:

I have discretized a linear system $$Au = f$$ for the following PDE problem:

$$-\frac{\partial u}{\partial x^{2}} -\frac{\partial u}{\partial y^{2}} = f$$

If for convenience we eliminate the boundary conditions and have the usual lexicographic ordering while using finite differences, then $$A$$ will have the same shape as described on page 4 here. Note that it is symmetric and positive definite. We will assume that it is of size $$n \times n$$.

Since it is a large sparse matrix, direct solvers are not a good method. So, I want to employ splitting based iterative solvers, i.e. we let $$A = M - N$$ which defines the following recursion:

$$u_{k+1} = M^{-1}Nu_{k} + M^{-1}f$$ where the starting vector is the zero vector.

My first question: If we let $$L$$ denote the incomplete Cholesky factor with zero fill-in, then since $$L \cdot L^{T}$$ is an approximation of A, can we take $$M = L \cdot L^{T}$$ in the above method?

Second question: If the previous answer is yes, how can I theoretically prove that the above method converges? I know that I need to prove that the spectral radius(i.e. maximum absolute eigenvalue) $$\rho (B_{M}) < 1$$, where $$B_{M} = I - M^{-1}A$$ is the iteration matrix. Especially since it is not clear what form $$M^{-1}$$ takes here, I am unsure of how to proceed.

The third question is closely related to the previous one. Given that the stopping criterion is $$\frac{\left\lVert f - Au_{k}\right\rVert _{2}}{\left\lVert f\right\rVert _{2}} \leq 10^{-10}$$, where $$f = -12x^{2}y^{5} - 20x^{4}y^{3}$$, I am also interested in quantifying the required amount of iterations for convergence. All I know is that the smaller $$\rho (B_{M})$$ is the faster it converges, but how to translate this into a minimum number of iterations?

The final question has to do with the amount of computational work per iteration, if we look at the above recursion then it is clear that at the very least a linear system solve with $$M$$ is required, but it is unclear to me what amount of work this requires.

Edit: I have indeed looked into the alternative direct methods suggested in the comments, but for me it's important to gain a good grasp for relatively simpler methods like the one described here. Especially for comparative purposes I'd like to be able to answer the above questions, so that I can better consider its weaknesses and see what kind of alternative methods work best. I still do appreciate the suggestions though!

New contributor
Drake Mannis is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
• Just a note (nothing to do with an answer to your question): sparse matrix direct solvers are now very efficient and widely used; often, outperforming iterative solvers (certainly, case-dependent). – Anton Menshov Jan 13 at 3:23
• In order to add something to @AntonMenshov comment, I suggest you look at MUMPS (MUltifrontal Massively Parallel sparse direct Solver), which I would really like to see if any iterative solver could beat it, cause I think it is really unbeatable! – Alone Programmer Jan 13 at 3:38
• I was under the impression that mumps performance suffers for above 1000 cores or so. Which is part of why iterative solvers are still used for HPC applications – EMP Jan 13 at 5:49
• Let's not turn it into a debate about direct vs iterative :) – Anton Menshov Jan 13 at 6:05
• I second @AntonMenshov's suggestion and just stick with the question at hand :-) For reference, one a single core, direct solvers are generally the faster choice for 2d problems until around 100,000 unknowns (plus minus a factor of 2 or 3 maybe). – Wolfgang Bangerth Jan 13 at 17:31