LU_decomposition Subroutine

public subroutine LU_decomposition(A, L, U)

LU decomposition of a matrix A This subroutine performs LU decomposition of a given matrix A, where L is a lower triangular 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)) :: U

Source Code

    subroutine LU_decomposition(A, L, U)

        real(dp), dimension(:, :), intent(in) :: A
        real(dp), dimension(size(A, 1), size(A, 1)), intent(out) :: L, U
        integer :: i, j, N

        N = size(A, 1)

        L = 0.d0
        U = 0.d0

        do j = 1, N
            L(j, j) = 1.d0

            do i = 1, j
                U(i, j) = A(i, j) - dot_product(L(i, 1:i - 1), U(1:i - 1, j))
            end do

            do i = j + 1, N
                L(i, j) = (A(i, j) - dot_product(L(i, 1:j - 1), U(1:j - 1, j))) / U(j, j)
            end do
        end do

    end subroutine LU_decomposition