Skip to content

Commit e9a3f5b

Browse files
perazzperazz
committed
add tests
1 parent 92b1ac5 commit e9a3f5b

File tree

2 files changed

+166
-0
lines changed

2 files changed

+166
-0
lines changed

‎test/linalg/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ set(
33
"test_linalg.fypp"
44
"test_blas_lapack.fypp"
55
"test_linalg_solve.fypp"
6+
"test_linalg_inverse.fypp"
67
"test_linalg_lstsq.fypp"
78
"test_linalg_determinant.fypp"
89
"test_linalg_matrix_property_checks.fypp"
@@ -12,6 +13,7 @@ fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
1213
ADDTEST(linalg)
1314
ADDTEST(linalg_determinant)
1415
ADDTEST(linalg_matrix_property_checks)
16+
ADDTEST(linalg_inverse)
1517
ADDTEST(linalg_solve)
1618
ADDTEST(linalg_lstsq)
1719
ADDTEST(blas_lapack)

‎test/linalg/test_linalg_inverse.fypp

Lines changed: 164 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,164 @@
1+
#:include "common.fypp"
2+
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
3+
! Test inverse matrix operator
4+
module test_linalg_inverse
5+
use testdrive, only: error_type, check, new_unittest, unittest_type
6+
use stdlib_linalg_constants
7+
use stdlib_linalg, only: inv,invert,operator(.inv.)
8+
use stdlib_linalg_state, only: linalg_state_type
9+
10+
implicit none (type,external)
11+
private
12+
13+
public :: test_inverse_matrix
14+
15+
contains
16+
17+
!> Matrix inversion tests
18+
subroutine test_inverse_matrix(tests)
19+
!> Collection of tests
20+
type(unittest_type), allocatable, intent(out) :: tests(:)
21+
22+
allocate(tests(0))
23+
24+
#:for rk,rt,ri in REAL_KINDS_TYPES
25+
#:if rk!="xdp"
26+
tests = [tests,new_unittest("inverse_${ri}$_eye_inverse",test_${ri}$_eye_inverse)]
27+
#:endif
28+
#:endfor
29+
30+
end subroutine test_inverse_matrix
31+
32+
!> Invert real identity matrix
33+
#:for rk,rt,ri in REAL_KINDS_TYPES
34+
#:if rk!="xdp"
35+
subroutine test_${ri}$_eye_inverse(error)
36+
type(error_type), allocatable, intent(out) :: error
37+
38+
type(linalg_state_type) :: state
39+
40+
integer(ilp), parameter :: n = 25_ilp
41+
integer(ilp) :: i,j
42+
${rt}$ :: a(n,n),inva(n,n)
43+
44+
do concurrent (i=1:n,j=1:n)
45+
a(i,j) = merge(1.0_${rk}$,0.0_${rk}$,i==j)
46+
end do
47+
48+
!> Inverse function
49+
inva = inv(a,err=state)
50+
call check(error,state%ok(),'inverse_${ri}$_eye (function): '//state%print())
51+
if (allocated(error)) return
52+
53+
call check(error,all(abs(a-inva)<epsilon(0.0_${rk}$)),'inverse_${ri}$_eye (function): data converged')
54+
if (allocated(error)) return
55+
56+
!> Inverse subroutine
57+
call invert(a,err=state)
58+
59+
call check(error,state%ok(),'inverse_${ri}$_eye (subroutine): '//state%print())
60+
if (allocated(error)) return
61+
62+
call check(error,all(abs(a-inva)<epsilon(0.0_${rk}$)),'inverse_${ri}$_eye (subroutine): data converged')
63+
if (allocated(error)) return
64+
65+
end subroutine test_${ri}$_eye_inverse
66+
67+
#:endif
68+
#:endfor
69+
70+
!> Invert complex identity matrix
71+
#:for ck,ct,ci in CMPLX_KINDS_TYPES
72+
#:if ck!="xdp"
73+
subroutine test_${ci}$_eye_inverse(error)
74+
type(error_type), allocatable, intent(out) :: error
75+
76+
type(linalg_state_type) :: state
77+
78+
integer(ilp) :: i,j,failed
79+
integer(ilp), parameter :: n = 25_ilp
80+
81+
${ct}$ :: a(n,n),copya(n,n),inva(n,n)
82+
83+
do concurrent (i=1:n,j=1:n)
84+
a(i,j) = merge((1.0_${ck}$,1.0_${ck}$),(0.0_${ck}$,0.0_${ck}$),i==j)
85+
end do
86+
copya = a
87+
88+
!> The inverse of a complex diagonal matrix has conjg(z_ii)/abs(z_ii)^2 on the diagonal
89+
inva = inv(a,err=state)
90+
91+
call check(error,state%ok(),'inverse_${ci}$_eye (function): '//state%print())
92+
if (allocated(error)) return
93+
94+
failed = 0
95+
do i=1,n
96+
do j=1,n
97+
if (.not.is_diagonal_inverse(a(i,j),inva(i,j),i,j)) failed = failed+1
98+
end do
99+
end do
100+
101+
call check(error,failed==0,'inverse_${ci}$_eye (function): data converged')
102+
if (allocated(error)) return
103+
104+
!> Inverse subroutine
105+
call invert(copya,err=state)
106+
107+
call check(error,state%ok(),'inverse_${ci}$_eye (subroutine): '//state%print())
108+
if (allocated(error)) return
109+
110+
failed = 0
111+
do i=1,n
112+
do j=1,n
113+
if (.not.is_diagonal_inverse(a(i,j),copya(i,j),i,j)) failed = failed+1
114+
end do
115+
end do
116+
117+
call check(error,failed==0,'inverse_${ci}$_eye (subroutine): data converged')
118+
if (allocated(error)) return
119+
120+
contains
121+
122+
elemental logical function is_diagonal_inverse(aij,invaij,i,j)
123+
${ct}$, intent(in) :: aij,invaij
124+
integer(ilp), intent(in) :: i,j
125+
if (i/=j) then
126+
is_diagonal_inverse = max(abs(aij),abs(invaij))<epsilon(0.0_${ck}$)
127+
else
128+
! Product should return the real identity
129+
is_diagonal_inverse = abs(aij*invaij - (1.0_${ck}$,0.0_${ck}$))<epsilon(0.0_${ck}$)
130+
endif
131+
end function is_diagonal_inverse
132+
133+
end subroutine test_${ci}$_eye_inverse
134+
135+
#:endif
136+
#:endfor
137+
138+
end module test_linalg_inverse
139+
140+
program test_inv
141+
use, intrinsic :: iso_fortran_env, only : error_unit
142+
use testdrive, only : run_testsuite, new_testsuite, testsuite_type
143+
use test_linalg_inverse, only : test_inverse_matrix
144+
implicit none
145+
integer :: stat, is
146+
type(testsuite_type), allocatable :: testsuites(:)
147+
character(len=*), parameter :: fmt = '("#", *(1x, a))'
148+
149+
stat = 0
150+
151+
testsuites = [ &
152+
new_testsuite("linalg_inverse", test_inverse_matrix) &
153+
]
154+
155+
do is = 1, size(testsuites)
156+
write(error_unit, fmt) "Testing:", testsuites(is)%name
157+
call run_testsuite(testsuites(is)%collect, error_unit, stat)
158+
end do
159+
160+
if (stat > 0) then
161+
write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
162+
error stop
163+
end if
164+
end program test_inv

0 commit comments

Comments
 (0)