Skip to content

Commit 9d35084

Browse files
perazzperazz
committed
copy in tests
1 parent f5de00f commit 9d35084

File tree

2 files changed

+211
-1
lines changed

2 files changed

+211
-1
lines changed

‎test/linalg/CMakeLists.txt

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ set(
66
"test_linalg_eigenvalues.fypp"
77
"test_linalg_solve.fypp"
88
"test_linalg_inverse.fypp"
9+
"test_linalg_pseudoinverse.fypp"
910
"test_linalg_lstsq.fypp"
1011
"test_linalg_norm.fypp"
1112
"test_linalg_determinant.fypp"
@@ -22,10 +23,11 @@ ADDTEST(linalg_determinant)
2223
ADDTEST(linalg_eigenvalues)
2324
ADDTEST(linalg_matrix_property_checks)
2425
ADDTEST(linalg_inverse)
26+
ADDTEST(linalg_pseudoinverse)
2527
ADDTEST(linalg_norm)
2628
ADDTEST(linalg_solve)
2729
ADDTEST(linalg_lstsq)
2830
ADDTEST(linalg_qr)
2931
ADDTEST(linalg_svd)
3032
ADDTEST(blas_lapack)
31-
ADDTEST(linalg_sparse)
33+
ADDTEST(linalg_sparse)

‎test/linalg/test_linalg_pseudoinverse.fypp

Lines changed: 208 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,208 @@
1+
#:include "common.fypp"
2+
! Test inverse matrix
3+
module test_linalg_pseudoinverse
4+
use stdlib_linalg_interface
5+
6+
implicit none (type,external)
7+
8+
contains
9+
10+
!> Matrix inversion tests
11+
subroutine test_pseudoinverse_matrix(error)
12+
logical, intent(out) :: error
13+
14+
real :: t0,t1
15+
16+
call cpu_time(t0)
17+
18+
#:for rk,rt,ri in REAL_KINDS_TYPES
19+
call test_${ri}$_eye_pseudoinverse(error)
20+
if (error) return
21+
#:endfor
22+
#:for ck,ct,ci in ALL_KINDS_TYPES
23+
call test_${ci}$_square_pseudoinverse(error)
24+
if (error) return
25+
call test_${ci}$_tall_pseudoinverse(error)
26+
if (error) return
27+
call test_${ci}$_wide_pseudoinverse(error)
28+
if (error) return
29+
call test_${ci}$_singular_pseudoinverse(error)
30+
if (error) return
31+
32+
#:endfor
33+
34+
call cpu_time(t1)
35+
36+
print 1, 1000*(t1-t0), merge('SUCCESS','ERROR ',.not.error)
37+
38+
1 format('Pseudo-Inverse matrix tests completed in ',f9.4,' milliseconds, result=',a)
39+
40+
end subroutine test_pseudoinverse_matrix
41+
42+
!> Invert identity matrix
43+
#:for rk,rt,ri in REAL_KINDS_TYPES
44+
subroutine test_${ri}$_eye_pseudoinverse(error)
45+
logical, intent(out) :: error
46+
47+
type(linalg_state) :: state
48+
49+
integer(ilp) :: i,j
50+
integer(ilp), parameter :: n = 15_ilp
51+
real(${rk}$), parameter :: tol = sqrt(epsilon(0.0_${rk}$))
52+
53+
${rt}$ :: a(n,n),inva(n,n)
54+
55+
do concurrent (i=1:n,j=1:n)
56+
a(i,j) = merge(1.0_${rk}$,0.0_${rk}$,i==j)
57+
end do
58+
59+
!> Invert function
60+
inva = pinv(a,err=state)
61+
error = state%error() .or. .not.all(abs(a-inva)<tol)
62+
if (error) return
63+
64+
!> Inverse subroutine
65+
call pseudoinvert(a,inva,err=state)
66+
error = state%error() .or. .not.all(abs(a-inva)<tol)
67+
68+
!> Operator
69+
inva = .pinv.a
70+
error = .not.all(abs(a-inva)<tol)
71+
72+
end subroutine test_${ri}$_eye_pseudoinverse
73+
74+
#:endfor
75+
76+
#:for ck,ct,ci in ALL_KINDS_TYPES
77+
78+
!> Test edge case: square matrix
79+
subroutine test_${ci}$_square_pseudoinverse(error)
80+
logical, intent(out) :: error
81+
82+
type(linalg_state) :: state
83+
84+
integer(ilp) :: failed
85+
integer(ilp), parameter :: n = 10
86+
real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$))
87+
${ct}$ :: a(n, n), inva(n, n)
88+
#:if ct.startswith('complex')
89+
real(${ck}$) :: rea(n, n, 2)
90+
91+
call random_number(rea)
92+
a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${ck}$)
93+
#:else
94+
95+
call random_number(a)
96+
#:endif
97+
98+
inva = pinv(a, err=state)
99+
error = state%error(); if (error) return
100+
101+
failed = count(abs(a - matmul(a, matmul(inva, a))) > tol)
102+
error = failed > 0; if (error) return
103+
104+
failed = count(abs(inva - matmul(inva, matmul(a, inva))) > tol)
105+
error = failed > 0; if (error) return
106+
107+
end subroutine test_${ci}$_square_pseudoinverse
108+
109+
!> Test edge case: tall matrix
110+
subroutine test_${ci}$_tall_pseudoinverse(error)
111+
logical, intent(out) :: error
112+
113+
type(linalg_state) :: state
114+
115+
integer(ilp) :: failed
116+
integer(ilp), parameter :: m = 20, n = 10
117+
real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$))
118+
${ct}$ :: a(m, n), inva(n, m)
119+
#:if ct.startswith('complex')
120+
real(${ck}$) :: rea(m, n, 2)
121+
122+
call random_number(rea)
123+
a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${ck}$)
124+
#:else
125+
126+
call random_number(a)
127+
#:endif
128+
129+
inva = pinv(a, err=state)
130+
error = state%error(); if (error) return
131+
132+
failed = count(abs(a - matmul(a, matmul(inva, a))) > tol)
133+
error = failed > 0; if (error) return
134+
135+
failed = count(abs(inva - matmul(inva, matmul(a, inva))) > tol)
136+
error = failed > 0; if (error) return
137+
138+
end subroutine test_${ci}$_tall_pseudoinverse
139+
140+
!> Test edge case: wide matrix
141+
subroutine test_${ci}$_wide_pseudoinverse(error)
142+
logical, intent(out) :: error
143+
144+
type(linalg_state) :: state
145+
146+
integer(ilp) :: failed
147+
integer(ilp), parameter :: m = 10, n = 20
148+
real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$))
149+
${ct}$ :: a(m, n), inva(n, m)
150+
#:if ct.startswith('complex')
151+
real(${ck}$) :: rea(m, n, 2)
152+
153+
call random_number(rea)
154+
a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${ck}$)
155+
#:else
156+
157+
call random_number(a)
158+
#:endif
159+
160+
inva = pinv(a, err=state)
161+
error = state%error(); if (error) return
162+
163+
failed = count(abs(a - matmul(a, matmul(inva, a))) > tol)
164+
error = failed > 0; if (error) return
165+
166+
failed = count(abs(inva - matmul(inva, matmul(a, inva))) > tol)
167+
error = failed > 0; if (error) return
168+
169+
end subroutine test_${ci}$_wide_pseudoinverse
170+
171+
!> Test edge case: singular matrix
172+
subroutine test_${ci}$_singular_pseudoinverse(error)
173+
logical, intent(out) :: error
174+
175+
type(linalg_state) :: state
176+
177+
integer(ilp) :: failed
178+
integer(ilp), parameter :: n = 10
179+
real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$))
180+
${ct}$ :: a(n, n), inva(n, n)
181+
#:if ct.startswith('complex')
182+
real(${ck}$) :: rea(n, n, 2)
183+
184+
call random_number(rea)
185+
a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${ck}$)
186+
#:else
187+
188+
call random_number(a)
189+
#:endif
190+
191+
! Make the matrix singular
192+
a(:, 1) = a(:, 2)
193+
194+
inva = pinv(a, err=state)
195+
196+
failed = count(abs(a - matmul(a, matmul(inva, a))) > tol)
197+
error = failed > 0; if (error) return
198+
199+
failed = count(abs(inva - matmul(inva, matmul(a, inva))) > tol)
200+
error = failed > 0; if (error) return
201+
202+
end subroutine test_${ci}$_singular_pseudoinverse
203+
204+
#:endfor
205+
206+
end module test_linalg_pseudoinverse
207+
208+

0 commit comments

Comments
 (0)