test_la_determinant.fypp

#:include "common.fypp"
! Test matrix determinant
module test_linalg_determinant
    use linear_algebra

    implicit none (type,external)

    contains

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

        real :: t0,t1

        call cpu_time(t0)

        #:for rk,rt,ri in ALL_KINDS_TYPES
        call test_${ri}$_eye_determinant(error)
        if (error) return

        call test_${ri}$_eye_multiple(error)
        if (error) return
        #: endfor

        #:for ck,ct,ci in CMPL_KINDS_TYPES
        call test_${ci}$_complex_determinant(error)
        if (error) return
        #: endfor

        call cpu_time(t1)

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

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

    end subroutine test_matrix_determinant

    !> Determinant of identity matrix
    #:for rk,rt,ri in ALL_KINDS_TYPES
    subroutine test_${ri}$_eye_determinant(error)
        logical, intent(out) :: error

        type(la_state) :: state

        integer(ilp) :: i
        integer(ilp), parameter :: n = 128_ilp

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

        a = eye(n)

        !> Determinant function
        deta = det(a,err=state)
        error = state%error() .or. .not.abs(deta-1.0_${rk}$)<tiny(0.0_${rk}$)
        if (error) return

    end subroutine test_${ri}$_eye_determinant

    !> Determinant of identity matrix multiplier
    subroutine test_${ri}$_eye_multiple(error)
        logical, intent(out) :: error

        type(la_state) :: state

        integer(ilp), parameter :: n = 10_ilp
        real(${rk}$), parameter :: coef = 0.0001_${rk}$
        integer(ilp) :: i,j
        ${rt}$ :: a(n,n),deta

        !> Multiply eye by a very small number
        a = eye(n)
        do concurrent (i=1:n)
          a(i,i) = coef
        end do

        !> Determinant: small, but a is not singular, because it is a multiple of the identity.
        deta = det(a,err=state)
        error = state%error() .or. .not.abs(deta-coef**n)<max(tiny(0.0_${rk}$),epsilon(0.0_${rk}$)*coef**n)
        if (error) return

    end subroutine test_${ri}$_eye_multiple

    #:endfor

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

        type(la_state) :: state

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

        ${ct}$, parameter :: res(nmax) = [${ct}$::(1,1),(0,2),(-2,2),(-4,0),(-4,-4), &
                                                  (0,-8),(8,-8),(16,0),(16,16),(0,32)]

        ${ct}$, allocatable :: a(:,:)
        ${ct}$ :: deta(nmax)

        !> Test determinant for all sizes, 1:nmax
        matrix_size: do n=1,nmax

           ! Put 1+i on each diagonal element
           a = eye(n)
           do concurrent (i=1:n)
             a(i,i) = (1.0_${ck}$,1.0_${ck}$)
           end do

           ! Expected result
           deta(n) = det(a,err=state)

           deallocate(a)
           if (state%error()) exit matrix_size

        end do matrix_size

        error = state%error() .or. any(.not.abs(res-deta)<=tiny(0.0_${ck}$))

    end subroutine test_${ci}$_complex_determinant

    #:endfor

end module test_linalg_determinant