ILU_decomposition Subroutine

public subroutine ILU_decomposition(A, L, U, level)

Incomplete LU decomposition of a matrix A This subroutine performs incomplete 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
integer, intent(in), optional :: level

Called by

proc~~ilu_decomposition~~CalledByGraph proc~ilu_decomposition ILU_decomposition proc~calculate_ilu_preconditioner Calculate_ILU_preconditioner proc~calculate_ilu_preconditioner->proc~ilu_decomposition

Source Code

    subroutine ILU_decomposition(A, L, U, level)

        real(dp), dimension(:, :), intent(in) :: A
        real(dp), dimension(size(A, 1), size(A, 1)), intent(out) :: L, U
        integer, optional, intent(in) :: level
        integer :: N, i, j
        integer, dimension(size(A, 1), size(A, 1)) :: fill_level
        logical, dimension(size(A, 1), size(A, 1)) :: S

        N = size(A, 1)

        L = 0.d0
        U = 0.d0

        if (present(level)) then
            call compute_fill_pattern_ILU(A, fill_level, level, N)
            S = (fill_level <= level)
        else
            S = A /= 0
        end if

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

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

            if (abs(U(j, j)) < 1.0e-12_dp) then
                print*,"Warning: Near-zero pivot at row ", j, ", value =", U(j, j)
                print*,"Replacing with small value, value =", sign(1.0e-12_dp, U(j, j))
                U(j, j) = sign(1.0e-12_dp, U(j, j))
            end if

            do i = j + 1, N
                if (S(i, j)) 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 ILU_decomposition