Skip to content

Commit ed0a4e0

Browse files
loiseaujcloiseaujc
committed
Jose's fix for the expm tests.
1 parent 5eb6b47 commit ed0a4e0

File tree

1 file changed

+91
-68
lines changed

1 file changed

+91
-68
lines changed

‎test/linalg/test_linalg_expm.fypp‎

Lines changed: 91 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,10 @@
1+
12
#:include "common.fypp"
23
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
3-
! Test Schur decomposition
4+
! Test the matrix exponential.
45
module test_linalg_expm
56
use testdrive, only: error_type, check, new_unittest, unittest_type
7+
use stdlib_constants
68
use stdlib_linalg_constants
79
use stdlib_linalg, only: expm, eye, norm, matrix_exp
810
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
@@ -14,99 +16,120 @@ module test_linalg_expm
1416

1517
contains
1618

17-
!> schur decomposition tests
19+
! gcc-15 bugfix utility
20+
subroutine add_test(tests,new_test)
21+
type(unittest_type), allocatable, intent(inout) :: tests(:)
22+
type(unittest_type), intent(in) :: new_test
23+
24+
integer :: n
25+
type(unittest_type), allocatable :: new_tests(:)
26+
27+
if (allocated(tests)) then
28+
n = size(tests)
29+
else
30+
n = 0
31+
end if
32+
33+
allocate(new_tests(n+1))
34+
if (n>0) new_tests(1:n) = tests(1:n)
35+
new_tests(1+n) = new_test
36+
call move_alloc(from=new_tests,to=tests)
37+
38+
end subroutine add_test
39+
40+
!> Exponent of matrix tests
1841
subroutine test_expm_computation(tests)
1942
!> Collection of tests
2043
type(unittest_type), allocatable, intent(out) :: tests(:)
2144

22-
allocate(tests(0))
23-
24-
#:for rk,rt,ri in RC_KINDS_TYPES
25-
tests = [tests, new_unittest("expm_${ri}$",test_expm_${ri}$)]
26-
tests = [tests, new_unittest("Error-handling expm_${ri}$",test_error_handling_expm_${ri}$)]
27-
#:endfor
45+
call add_test(tests,new_unittest("expm",test_expm))
46+
call add_test(tests,new_unittest("error_handling_expm",test_error_handling_expm))
2847

2948
end subroutine test_expm_computation
3049

3150
!> Matrix exponential with analytic expression.
32-
#:for rk,rt,ri in RC_KINDS_TYPES
33-
subroutine test_expm_${ri}$(error)
51+
subroutine test_expm(error)
3452
type(error_type), allocatable, intent(out) :: error
3553
! Problem dimension.
3654
integer(ilp), parameter :: n = 5, m = 6
3755
! Test matrix.
38-
${rt}$ :: A(n, n), E(n, n), Eref(n, n)
39-
real(${rk}$) :: err
4056
integer(ilp) :: i, j
4157

42-
! Initialize matrix.
43-
A = 0.0_${rk}$
44-
do i = 1, n-1
45-
A(i, i+1) = m*1.0_${rk}$
46-
enddo
47-
48-
! Reference with analytical exponential.
49-
Eref = eye(n, mold=1.0_${rk}$)
50-
do i = 1, n-1
51-
do j = 1, n-i
52-
Eref(i, i+j) = Eref(i, i+j-1)*m/j
58+
#:for rk,rt,ri in RC_KINDS_TYPES
59+
block
60+
${rt}$ :: A(n, n), E(n, n), Eref(n, n)
61+
real(${rk}$) :: err
62+
63+
! Initialize matrix.
64+
A = zero_${rk}$
65+
do i = 1, n-1
66+
A(i, i+1) = m*one_${rk}$
67+
enddo
68+
69+
! Reference with analytical exponential
70+
Eref = eye(n, mold=one_${rk}$)
71+
do i = 1, n-1
72+
do j = 1, n-i
73+
Eref(i, i+j) = Eref(i, i+j-1)*m/j
74+
enddo
5375
enddo
54-
enddo
5576

56-
! Compute matrix exponential.
57-
E = expm(A)
77+
! Compute matrix exponential.
78+
E = expm(A)
5879

59-
! Check result.
60-
err = norm(Eref - E, "inf")
61-
call check(error, err < (n**2)*epsilon(1.0_${rk}$), "Analytical matrix exponential.")
62-
if (allocated(error)) return
63-
return
64-
end subroutine test_expm_${ri}$
65-
#:endfor
80+
! Check result.
81+
err = norm(Eref - E, "inf")
82+
print *, err , (n**2)*epsilon(1.0_${rk}$)
83+
call check(error, err < (n**2)*epsilon(1.0_${rk}$), "Analytical matrix exponential.")
84+
if (allocated(error)) return
85+
end block
86+
#:endfor
87+
end subroutine test_expm
6688

6789
!> Test error handler.
68-
#:for rk,rt,ri in RC_KINDS_TYPES
69-
subroutine test_error_handling_expm_${ri}$(error)
90+
subroutine test_error_handling_expm(error)
7091
type(error_type), allocatable, intent(out) :: error
7192
! Problem dimension.
7293
integer(ilp), parameter :: n = 5, m = 6
7394
! Test matrix.
74-
${rt}$ :: A(n, n), E(n, n)
95+
7596
type(linalg_state_type) :: err
7697
integer(ilp) :: i
7798

78-
! Initialize matrix.
79-
A = 0.0_${rk}$
80-
do i = 1, n-1
81-
A(i, i+1) = m*1.0_${rk}$
82-
enddo
83-
84-
! Compute matrix exponential.
85-
call matrix_exp(A, E, order=-1, err=err)
86-
! Check result.
87-
call check(error, err%error(), "Negative Pade order")
88-
if (allocated(error)) return
89-
90-
call matrix_exp(A, order=-1, err=err)
91-
! Check result.
92-
call check(error, err%error(), "Negative Pade order")
93-
if (allocated(error)) return
94-
95-
! Compute matrix exponential.
96-
call matrix_exp(A, E(:n, :n-1), err=err)
97-
! Check result.
98-
call check(error, err%error(), "Invalid matrix size")
99-
if (allocated(error)) return
100-
101-
! Compute matrix exponential.
102-
call matrix_exp(A(:n, :n-1), err=err)
103-
! Check result.
104-
call check(error, err%error(), "Invalid matrix size")
105-
if (allocated(error)) return
106-
107-
return
108-
end subroutine test_error_handling_expm_${ri}$
109-
#:endfor
99+
#:for rk,rt,ri in RC_KINDS_TYPES
100+
block
101+
${rt}$ :: A(n, n), E(n, n)
102+
! Initialize matrix.
103+
A = zero_${rk}$
104+
do i = 1, n-1
105+
A(i, i+1) = m*one_${rk}$
106+
enddo
107+
108+
! Compute matrix exponential.
109+
call matrix_exp(A, E, order=-1, err=err)
110+
! Check result.
111+
call check(error, err%error(), "Negative Pade order")
112+
if (allocated(error)) return
113+
114+
call matrix_exp(A, order=-1, err=err)
115+
! Check result.
116+
call check(error, err%error(), "Negative Pade order")
117+
if (allocated(error)) return
118+
119+
! Compute matrix exponential.
120+
call matrix_exp(A, E(:n, :n-1), err=err)
121+
! Check result.
122+
call check(error, err%error(), "Invalid matrix size")
123+
if (allocated(error)) return
124+
125+
! Compute matrix exponential.
126+
call matrix_exp(A(:n, :n-1), err=err)
127+
! Check result.
128+
call check(error, err%error(), "Invalid matrix size")
129+
if (allocated(error)) return
130+
end block
131+
#:endfor
132+
end subroutine test_error_handling_expm
110133

111134
end module test_linalg_expm
112135

0 commit comments

Comments
 (0)