An incomplete Cholesky (IC) factorization can be used as a preconditioning matrix for the conjugate gradient algorithm, and is often the best choice if
the matrix a is large, sparse (many 0 entries) and symmetric-positive definite (A=A^T, x^TAx > 0). The IC factorization essentially just ignores the 0 entries of the matrix
A and thus the preconditioner matrix M is able to maintain A's sparse structure. One can use the form M=LL^T where the lower triangular matrix L stores the
square root of the diagonal of A, or one can use the root free form M=LDL^T where D stores the diagonal of A and L has 1 diagonal. Below we give a simple Fortran subroutine for computing the M=LL^T form, for more information on this version there is
an excellent book which provides this algorithm "Matrix Computations" by Golub and van Loan. For a parallel version, one can use the M=LDL^T version with a
red black ordering and the algorithm provided in Ortega's book "Introduction to Parallel and Vector Solution of Linear Systems".
SUBROUTINE IC(A,N) ! Subroutine overwrites the sparse SPD matrix A ! with L from its IC factorization LL^T Implicit None Integer(4) :: i,j,k,N Real(8), Dimension(N,N) :: A Do k=1,N A(k,k)=sqrt(A(k,k)) Do i=k+1,N If (A(i,k).ne.0) Then A(i,k)=A(i,k)/A(k,k) End If End Do Do j=k+1,N Do i=j,N If (A(i,j).ne.0) Then A(i,j)=A(i,j)-A(i,k)*A(j,k) End If End Do End Do End Do END SUBROUTINE