test_la_inverse.fypp

#:include "common.fypp"
! Test inverse matrix
module test_linalg_inverse
    use linear_algebra

    implicit none (type,external)

    contains

    !> Matrix inversion tests
    subroutine test_inverse_matrix(error)
        logical, intent(out) :: error

        real :: t0,t1

        call cpu_time(t0)

        #:for rk,rt,ri in REAL_KINDS_TYPES
        call test_${ri}$_eye_inverse(error)
        if (error) return
        #: endfor
        #:for ck,ct,ci in CMPL_KINDS_TYPES
        call test_${ci}$_eye_inverse(error)
        if (error) return
        #: endfor

        call cpu_time(t1)

        print 1, 1000*(t1-t0), merge('SUCCESS','ERROR  ',.not.error)

        1 format('Inverse matrix tests completed in ',f9.4,' milliseconds, result=',a)

    end subroutine test_inverse_matrix

    !> Invert identity matrix
    #:for rk,rt,ri in REAL_KINDS_TYPES
    subroutine test_${ri}$_eye_inverse(error)
        logical, intent(out) :: error

        type(la_state) :: state

        integer(ilp) :: i,j
        integer(ilp), parameter :: n = 250_ilp

        ${rt}$ :: a(n,n),inva(n,n)

        do concurrent (i=1:n,j=1:n)
          a(i,j) = merge(1.0_${rk}$,0.0_${rk}$,i==j)
        end do

        !> Invert function
        inva = inv(a,err=state)
        error = state%error() .or. .not.all(abs(a-inva)<tiny(0.0_${rk}$))
        if (error) return

        !> Inverse subroutine
        call invert(a,err=state)
        error = state%error() .or. .not.all(abs(a-inva)<tiny(0.0_${rk}$))

    end subroutine test_${ri}$_eye_inverse

    #:endfor

    !> Invert identity matrix
    #:for ck,ct,ci in CMPL_KINDS_TYPES
    subroutine test_${ci}$_eye_inverse(error)
        logical, intent(out) :: error

        type(la_state) :: state

        integer(ilp) :: i,j,failed
        integer(ilp), parameter :: n = 250_ilp

        ${ct}$ :: a(n,n),copya(n,n),inva(n,n)

        do concurrent (i=1:n,j=1:n)
          a(i,j) = merge((1.0_${ck}$,1.0_${ck}$),(0.0_${ck}$,0.0_${ck}$),i==j)
        end do
        copya = a

        !> The inverse of a complex diagonal matrix has conjg(z_ii)/abs(z_ii)^2 on the diagonal
        inva = inv(a,err=state)

        failed = 0
        do i=1,n
            do j=1,n
                if (.not.is_diagonal_inverse(a(i,j),inva(i,j),i,j)) failed = failed+1
            end do
        end do

        error = state%error() .or. .not.failed==0
        if (error) return

        !> Inverse subroutine
        call invert(copya,err=state)

        failed = 0
        do i=1,n
            do j=1,n
                if (.not.is_diagonal_inverse(a(i,j),copya(i,j),i,j)) failed = failed+1
            end do
        end do

        error = state%error() .or. .not.failed==0

        contains

           elemental logical function is_diagonal_inverse(aij,invaij,i,j)
               ${ct}$, intent(in) :: aij,invaij
               integer(ilp), intent(in) :: i,j
               if (i/=j) then
                  is_diagonal_inverse = max(abs(aij),abs(invaij))<tiny(0.0_${ck}$)
               else
                  ! Product should return the real identity
                  is_diagonal_inverse = abs(aij*invaij - (1.0_${ck}$,0.0_${ck}$))<tiny(0.0_${ck}$)
               endif
           end function is_diagonal_inverse

    end subroutine test_${ci}$_eye_inverse

    #:endfor

end module test_linalg_inverse