test_linalg_expm.fypp

#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
! Test Schur decomposition
module test_linalg_expm
    use testdrive, only: error_type, check, new_unittest, unittest_type
    use stdlib_linalg_constants
    use stdlib_linalg, only: expm, eye, mnorm

    implicit none (type,external)
    
    public :: test_expm_computation

    contains

    !> schur decomposition tests
    subroutine test_expm_computation(tests)
        !> Collection of tests
        type(unittest_type), allocatable, intent(out) :: tests(:)
        
        allocate(tests(0))
        
        #:for rk,rt,ri in RC_KINDS_TYPES
        tests = [tests, new_unittest("expm_${ri}$",test_expm_${ri}$)]
        #:endfor

    end subroutine test_expm_computation

    !> Matrix exponential with analytic expression.
    #:for rk,rt,ri in RC_KINDS_TYPES
    subroutine test_expm_${ri}$(error)
        type(error_type), allocatable, intent(out) :: error
        ! Problem dimension.
        integer(ilp), parameter :: n = 5, m = 6
        ! Test matrix.
        ${rt}$ :: A(n, n), E(n, n), Eref(n, n)
        real(${rk}$) :: err
        integer(ilp) :: i, j

        ! Initialize matrix.
        A = 0.0_${rk}$
        do i = 1, n-1
            A(i, i+1) = m*1.0_${rk}$
        enddo

        ! Reference with analytical exponential.
        Eref = eye(n, mold=1.0_${rk}$)
        do i = 1, n-1
            do j = 1, n-i
                Eref(i, i+j) = Eref(i, i+j-1)*m/j
            enddo
        enddo

        ! Compute matrix exponential.
        E = expm(A)

        ! Check result.
        err = mnorm(Eref - E, "inf")
        call check(error, err < (n**2)*epsilon(1.0_${rk}$), "Analytical matrix exponential.")
        if (allocated(error)) return
        return        
    end subroutine test_expm_${ri}$    
    #:endfor    

end module test_linalg_expm

program test_expm
    use, intrinsic :: iso_fortran_env, only : error_unit
    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
    use test_linalg_expm, only : test_expm_computation
    implicit none
    integer :: stat, is
    type(testsuite_type), allocatable :: testsuites(:)
    character(len=*), parameter :: fmt = '("#", *(1x, a))'

    stat = 0

    testsuites = [ &
        new_testsuite("linalg_expm", test_expm_computation) &
        ]

    do is = 1, size(testsuites)
        write(error_unit, fmt) "Testing:", testsuites(is)%name
        call run_testsuite(testsuites(is)%collect, error_unit, stat)
    end do

    if (stat > 0) then
        write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
        error stop
    end if
end program test_expm