Incomplete Cholesky decomposition of a matrix A This subroutine performs incomplete Cholesky decomposition of a given matrix A, where L is a lower triangular matrix and U is an upper triangular matrix.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=dp), | intent(in), | DIMENSION(:, :) | :: | A | ||
real(kind=dp), | intent(out), | DIMENSION(SIZE(A, 1), SIZE(A, 1)) | :: | L |
SUBROUTINE Incomplete_Cholesky_decomposition(A, L) REAL(dp), DIMENSION(:, :), INTENT(IN) :: A REAL(dp), DIMENSION(SIZE(A, 1), SIZE(A, 1)), INTENT(OUT) :: L LOGICAL, DIMENSION(SIZE(A, 1), SIZE(A, 1)) :: S INTEGER :: i, j, N N = SIZE(A, 1) L = Identity_n(N) S = A /= 0 DO i = 1, N DO j = 1, i-1 IF (S(i,j)) L(i,j) = (A(i,j) - DOT_PRODUCT(L(i, 1:j-1), L(j, 1:j-1))) / L(j,j) END DO IF (S(i,i)) L(i,i) = SQRT(A(i,i) - DOT_PRODUCT(L(i, 1:i-1), L(i, 1:i-1))) END DO END SUBROUTINE Incomplete_Cholesky_decomposition