The solver

For large sparse matrices, direct methods that build a full decomposition of the matrix are usually unprofitable: during the decomposition, new nonzero elements appear in the positions of the original zeros — the so-called fill-in. As a result, the main advantage of a sparse matrix is lost — a small memory footprint and a small number of arithmetic operations.

We will solve system (3.1) by the iterative preconditioned conjugate gradient method (PCG). The method applies to symmetric positive definite matrices — exactly those that arise in FEM problems with elliptic operators and in implicit heat conduction schemes. Let us write out the algorithm step by step.

Let x0x_0 =0= 0 be the initial guess. The initial residual takes the form

r0\displaystyle r_0 =b\displaystyle = b −A⋅x0.\displaystyle - A \cdot x_0.
(3.2)

The initial search direction and the preconditioning vector have the form

P⋅z0\displaystyle P \cdot z_0 =r0,p0\displaystyle = r_0, \qquad p_0 =z0,\displaystyle = z_0,
(3.3)

where PP is the preconditioner matrix: it approximates AA, but systems with it are much cheaper to solve. Then at every iteration kk =0,1,2,…= 0, 1, 2, \ldots we compute

αk\displaystyle \alpha_k =rkT⋅zkpkT⋅A⋅pk,\displaystyle = \frac{r_k^T \cdot z_k}{p_k^T \cdot A \cdot p_k},
(3.4)
xk+1\displaystyle x_{k+1} =xk\displaystyle = x_k +αk⋅pk,\displaystyle + \alpha_k \cdot p_k,
(3.5)
rk+1\displaystyle r_{k+1} =rk\displaystyle = r_k −αk⋅A⋅pk,\displaystyle - \alpha_k \cdot A \cdot p_k,
(3.6)
P⋅zk+1\displaystyle P \cdot z_{k+1} =rk+1,\displaystyle = r_{k+1},
(3.7)
βk\displaystyle \beta_k =rk+1T⋅zk+1rkT⋅zk,\displaystyle = \frac{r_{k+1}^T \cdot z_{k+1}}{r_k^T \cdot z_k},
(3.8)
pk+1\displaystyle p_{k+1} =zk+1\displaystyle = z_{k+1} +βk⋅pk.\displaystyle + \beta_k \cdot p_k.
(3.9)

The iterations continue until the squared norm of the residual drops below a given threshold

rk+1T⋅rk+1\displaystyle r_{k+1}^T \cdot r_{k+1} <ε.\displaystyle < \varepsilon.
(3.10)

The main computational cost of each iteration is one multiplication of the sparse matrix by a vector, A⋅pkA \cdot p_k. This is precisely why the matrix is converted from COO to CSR before the iterations start: in this format the multiplication of a row by the vector runs only over the nonzero elements.

Flowchart of one iteration of the preconditioned conjugate gradient method
Fig. 3.2. One PCG iteration: the step A⋅pkA \cdot p_k is the most expensive one.

Next, let us talk about preconditioners.