Incomplete Cholesky Factorization

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 
         


 


 

Contact Info

E.mail: aduffy@math.fsu.edu