LDU_decomposition Subroutine

public subroutine LDU_decomposition(A, L, D, U)

LDU decomposition of a matrix A This subroutine performs LDU decomposition of a given matrix A, where L is a lower triangular matrix, D is a diagonal matrix, and U is an upper triangular matrix.

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), DIMENSION(:, :) :: A
real(kind=dp), intent(out), DIMENSION(SIZE(A, 1),SIZE(A, 1)) :: L
real(kind=dp), intent(out), DIMENSION(SIZE(A, 1),SIZE(A, 1)) :: D
real(kind=dp), intent(out), DIMENSION(SIZE(A, 1),SIZE(A, 1)) :: U

Source Code

    SUBROUTINE LDU_decomposition(A, L, D, U)

        REAL(dp),DIMENSION(:, :), INTENT(IN) :: A
        REAL(dp),DIMENSION(SIZE(A, 1),SIZE(A, 1)), INTENT(OUT) :: L, U, D
        INTEGER :: i, j, k, N

        N = SIZE(A, 1)

        L = 0.d0
        D = 0.d0
        U = 0.d0

        DO j = 1, N
            L(j, j) = 1.d0
            U(j, j) = 1.d0
            
            DO i = 1, j-1
                U(i, j) = (A(i, j) - DOT_PRODUCT(L(i, 1:i-1), U(1:i-1, j) * [ (D(k,k), k = 1, i-1) ])) / D(i, i)
            END DO

            i = j
            D(j, j) = A(j, j) - DOT_PRODUCT(L(j, 1:j-1), U(1:j-1, j) * [ (D(k,k), k = 1, j-1) ])

            DO i = j+1, N
                L(i, j) = (A(i, j) - DOT_PRODUCT(L(i, 1:j-1), U(1:j-1, j) * [ (D(k,k), k = 1, j-1) ])) / D(j, j)
            END DO
        END DO

    END SUBROUTINE LDU_decomposition