## Incomplete Cholesky Factorization

Author: Austen Duffy, Florida State University

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

``````