Skip to content

Commit fb52445

Browse files
loiseaujcloiseaujc
committed
Add tests for expm.
1 parent ead54b1 commit fb52445

File tree

2 files changed

+161
-0
lines changed

2 files changed

+161
-0
lines changed

‎test/linalg/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@ set(
1717
"test_linalg_sparse.fypp"
1818
"test_linalg_specialmatrices.fypp"
1919
"test_linalg_cholesky.fypp"
20+
"test_linalg_expm.fypp"
2021
)
2122

2223
# Preprocessed files to contain preprocessor directives -> .F90
@@ -32,6 +33,7 @@ ADDTEST(linalg)
3233
ADDTEST(linalg_cholesky)
3334
ADDTEST(linalg_determinant)
3435
ADDTEST(linalg_eigenvalues)
36+
ADDTEST(linalg_expm)
3537
ADDTEST(linalg_matrix_property_checks)
3638
ADDTEST(linalg_inverse)
3739
ADDTEST(linalg_pseudoinverse)

‎test/linalg/test_linalg_expm.fypp‎

Lines changed: 159 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,159 @@
1+
#:include "common.fypp"
2+
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
3+
! Test the matrix exponential.
4+
module test_linalg_expm
5+
use testdrive, only: error_type, check, new_unittest, unittest_type
6+
use stdlib_constants
7+
use stdlib_linalg_constants
8+
use stdlib_linalg, only: expm, eye, norm, matrix_exp
9+
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
10+
LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
11+
12+
implicit none (type,external)
13+
14+
public :: test_expm_computation
15+
16+
contains
17+
18+
! gcc-15 bugfix utility
19+
subroutine add_test(tests,new_test)
20+
type(unittest_type), allocatable, intent(inout) :: tests(:)
21+
type(unittest_type), intent(in) :: new_test
22+
23+
integer :: n
24+
type(unittest_type), allocatable :: new_tests(:)
25+
26+
if (allocated(tests)) then
27+
n = size(tests)
28+
else
29+
n = 0
30+
end if
31+
32+
allocate(new_tests(n+1))
33+
if (n>0) new_tests(1:n) = tests(1:n)
34+
new_tests(1+n) = new_test
35+
call move_alloc(from=new_tests,to=tests)
36+
37+
end subroutine add_test
38+
39+
!> Exponent of matrix tests
40+
subroutine test_expm_computation(tests)
41+
!> Collection of tests
42+
type(unittest_type), allocatable, intent(out) :: tests(:)
43+
44+
call add_test(tests,new_unittest("expm",test_expm))
45+
call add_test(tests,new_unittest("error_handling_expm",test_error_handling_expm))
46+
47+
end subroutine test_expm_computation
48+
49+
!> Matrix exponential with analytic expression.
50+
subroutine test_expm(error)
51+
type(error_type), allocatable, intent(out) :: error
52+
! Problem dimension.
53+
integer(ilp), parameter :: n = 5, m = 6
54+
! Test matrix.
55+
integer(ilp) :: i, j
56+
57+
#:for rk,rt,ri in RC_KINDS_TYPES
58+
block
59+
${rt}$ :: A(n, n), E(n, n), Eref(n, n)
60+
real(${rk}$) :: err
61+
62+
! Initialize matrix.
63+
A = zero_${rk}$
64+
do i = 1, n-1
65+
A(i, i+1) = m*one_${rk}$
66+
enddo
67+
68+
! Reference with analytical exponential
69+
Eref = eye(n, mold=one_${rk}$)
70+
do i = 1, n-1
71+
do j = 1, n-i
72+
Eref(i, i+j) = Eref(i, i+j-1)*m/j
73+
enddo
74+
enddo
75+
76+
! Compute matrix exponential.
77+
E = expm(A)
78+
79+
! Check result.
80+
err = norm(Eref - E, "inf")
81+
print *, err , (n**2)*epsilon(1.0_${rk}$)
82+
call check(error, err < (n**2)*epsilon(1.0_${rk}$), "Analytical matrix exponential.")
83+
if (allocated(error)) return
84+
end block
85+
#:endfor
86+
end subroutine test_expm
87+
88+
!> Test error handler.
89+
subroutine test_error_handling_expm(error)
90+
type(error_type), allocatable, intent(out) :: error
91+
! Problem dimension.
92+
integer(ilp), parameter :: n = 5, m = 6
93+
! Test matrix.
94+
95+
type(linalg_state_type) :: err
96+
integer(ilp) :: i
97+
98+
#:for rk,rt,ri in RC_KINDS_TYPES
99+
block
100+
${rt}$ :: A(n, n), E(n, n)
101+
! Initialize matrix.
102+
A = zero_${rk}$
103+
do i = 1, n-1
104+
A(i, i+1) = m*one_${rk}$
105+
enddo
106+
107+
! Compute matrix exponential.
108+
call matrix_exp(A, E, order=-1, err=err)
109+
! Check result.
110+
call check(error, err%error(), "Negative Pade order")
111+
if (allocated(error)) return
112+
113+
call matrix_exp(A, order=-1, err=err)
114+
! Check result.
115+
call check(error, err%error(), "Negative Pade order")
116+
if (allocated(error)) return
117+
118+
! Compute matrix exponential.
119+
call matrix_exp(A, E(:n, :n-1), err=err)
120+
! Check result.
121+
call check(error, err%error(), "Invalid matrix size")
122+
if (allocated(error)) return
123+
124+
! Compute matrix exponential.
125+
call matrix_exp(A(:n, :n-1), err=err)
126+
! Check result.
127+
call check(error, err%error(), "Invalid matrix size")
128+
if (allocated(error)) return
129+
end block
130+
#:endfor
131+
end subroutine test_error_handling_expm
132+
133+
end module test_linalg_expm
134+
135+
program test_expm
136+
use, intrinsic :: iso_fortran_env, only : error_unit
137+
use testdrive, only : run_testsuite, new_testsuite, testsuite_type
138+
use test_linalg_expm, only : test_expm_computation
139+
implicit none
140+
integer :: stat, is
141+
type(testsuite_type), allocatable :: testsuites(:)
142+
character(len=*), parameter :: fmt = '("#", *(1x, a))'
143+
144+
stat = 0
145+
146+
testsuites = [ &
147+
new_testsuite("linalg_expm", test_expm_computation) &
148+
]
149+
150+
do is = 1, size(testsuites)
151+
write(error_unit, fmt) "Testing:", testsuites(is)%name
152+
call run_testsuite(testsuites(is)%collect, error_unit, stat)
153+
end do
154+
155+
if (stat > 0) then
156+
write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
157+
error stop
158+
end if
159+
end program test_expm

0 commit comments

Comments
 (0)