下三角行列なのかジャッジ

行列 A を引数にして call judge_lower_triangle(A) で下三角行列が判定するサブルーチン。

下三角成分に 0 があると Ax=b を解くときにランク落ちになってしまうので、下三角成分がノンゼロという厳しめの制約を課した。

subroutine judge_lower_triangle( A )
implicit none

    real(8), intent(in) :: A(:,:)
    integer :: i, j, flag=0

    if( size(A,1) /= size(A,2) ) stop "input matrix is not a square"

    do i = 1, size(A,1)
        do j = 1, size(A,2)
            if( j <= i ) then
                if( A(i,j)==0.d0 ) then
                    print*, "BAD component was found : ", i, j, A(i,j)
                    flag=1
                end if
            else
                if( A(i,j)/=0.d0 ) then
                    print*, "BAD component was found : ", i, j, A(i,j)
                    flag=1
                end if
            end if
        end do
    end do

    if( flag==1 ) stop "input matrix is NOT lower triangle"

end subroutine judge_lower_triangle