Skip to content

Commit 2d58c68

Browse files
perazzperazz
committed
add tests
1 parent 6ce5172 commit 2d58c68

File tree

3 files changed

+214
-0
lines changed

3 files changed

+214
-0
lines changed

‎src/stdlib_linalg.fypp

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,10 @@ module stdlib_linalg
1919
public :: det
2020
public :: operator(.det.)
2121
public :: diag
22+
public :: eig
23+
public :: eigh
24+
public :: eigvals
25+
public :: eigvalsh
2226
public :: eye
2327
public :: lstsq
2428
public :: lstsq_space

‎test/linalg/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@ set(
22
fppFiles
33
"test_linalg.fypp"
44
"test_blas_lapack.fypp"
5+
"test_linalg_eigenvalues.fypp"
56
"test_linalg_solve.fypp"
67
"test_linalg_lstsq.fypp"
78
"test_linalg_determinant.fypp"
@@ -11,6 +12,7 @@ fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
1112

1213
ADDTEST(linalg)
1314
ADDTEST(linalg_determinant)
15+
ADDTEST(linalg_eigenvalues)
1416
ADDTEST(linalg_matrix_property_checks)
1517
ADDTEST(linalg_solve)
1618
ADDTEST(linalg_lstsq)

‎test/linalg/test_linalg_eigenvalues.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 eigenvalues and eigendecompositions
3+
module test_linalg_eigenvalues
4+
use stdlib_linalg_constants
5+
use stdlib_linalg_state
6+
use stdlib_linalg, only: eig, eigh, eigvals, eigvalsh, diag
7+
use testdrive, only: error_type, check, new_unittest, unittest_type
8+
9+
implicit none (type,external)
10+
private
11+
12+
public :: test_eig_eigh
13+
14+
contains
15+
16+
!> SVD tests
17+
subroutine test_eig_eigh(tests)
18+
!> Collection of tests
19+
type(unittest_type), allocatable, intent(out) :: tests(:)
20+
21+
allocate(tests(0))
22+
23+
#:for rk,rt,ri in REAL_KINDS_TYPES
24+
#:if rk!="xdp"
25+
tests = [tests,new_unittest("test_eig_real_${ri}$",test_eig_real_${ri}$), &
26+
new_unittest("test_eigh_real_${ri}$",test_eigh_real_${ri}$)]
27+
#:endif
28+
#: endfor
29+
30+
#:for ck,ct,ci in CMPLX_KINDS_TYPES
31+
tests = [tests,new_unittest("test_eig_complex_${ci}$",test_eig_complex_${ci}$)]
32+
#: endfor
33+
34+
end subroutine test_eig_eigh
35+
36+
!> Simple real matrix eigenvalues
37+
#:for rk,rt,ri in REAL_KINDS_TYPES
38+
#:if rk!="xdp"
39+
subroutine test_eig_real_${ri}$(error)
40+
type(error_type), allocatable, intent(out) :: error
41+
42+
!> Reference solution
43+
real(${rk}$), parameter :: zero = 0.0_${rk}$
44+
real(${rk}$), parameter :: two = 2.0_${rk}$
45+
real(${rk}$), parameter :: sqrt2o2 = sqrt(two)*0.5_${rk}$
46+
real(${rk}$), parameter :: tol = sqrt(epsilon(zero))
47+
48+
!> Local variables
49+
type(linalg_state_type) :: state
50+
${rt}$ :: A(3,3),B(2,2)
51+
complex(${rk}$) :: lambda(3),Bvec(2,2),Bres(2,2)
52+
53+
!> Matrix with real eigenvalues
54+
A = reshape([1,0,0, &
55+
0,2,0, &
56+
0,0,3],[3,3])
57+
58+
call eig(A,lambda,err=state)
59+
60+
call check(error,state%ok(),state%print())
61+
if (allocated(error)) return
62+
63+
call check(error, all(aimag(lambda)==zero.and.real(lambda,kind=${rk}$)==[1,2,3]),'expected results')
64+
if (allocated(error)) return
65+
66+
!> Matrix with complex eigenvalues
67+
B = transpose(reshape([1, -1, &
68+
1, 1],[2,2]))
69+
70+
!> Expected right eigenvectors
71+
Bres(1,1:2) = sqrt2o2
72+
Bres(2,1) = cmplx(zero,-sqrt2o2,kind=${rk}$)
73+
Bres(2,2) = cmplx(zero,+sqrt2o2,kind=${rk}$)
74+
75+
call eig(B,lambda,right=Bvec,err=state)
76+
77+
call check(error,state%ok(),state%print())
78+
if (allocated(error)) return
79+
80+
call check(error, all(abs(Bres-Bvec)<=tol),'expected results')
81+
if (allocated(error)) return
82+
83+
end subroutine test_eig_real_${ri}$
84+
85+
! Symmetric matrix eigenvalues
86+
subroutine test_eigh_real_${ri}$(error)
87+
type(error_type), allocatable, intent(out) :: error
88+
89+
!> Reference solution
90+
real(${rk}$), parameter :: zero = 0.0_${rk}$
91+
real(${rk}$), parameter :: tol = sqrt(epsilon(zero))
92+
real(${rk}$), parameter :: A(4,4) = reshape([6,3,1,5, &
93+
3,0,5,1, &
94+
1,5,6,2, &
95+
5,1,2,2],[4,4])
96+
97+
!> Local variables
98+
real(${rk}$) :: Amat(4,4),lambda(4),vect(4,4),Av(4,4),lv(4,4)
99+
type(linalg_state_type) :: state
100+
101+
Amat = A
102+
103+
call eigh(Amat,lambda,vect,err=state)
104+
105+
Av = matmul(A,vect)
106+
lv = matmul(vect,diag(lambda))
107+
108+
call check(error,state%ok(),state%print())
109+
if (allocated(error)) return
110+
111+
call check(error, all(abs(Av-lv)<=tol*abs(Av)),'expected results')
112+
if (allocated(error)) return
113+
114+
!> Test functional versions: no state interface
115+
lambda = eigvalsh(Amat)
116+
117+
!> State interface
118+
lambda = eigvalsh(Amat,err=state)
119+
120+
call check(error,state%ok(),state%print())
121+
if (allocated(error)) return
122+
123+
!> Functional version, lower A
124+
Amat = A
125+
lambda = eigvalsh(Amat,upper_a=.false.,err=state)
126+
127+
call check(error,state%ok(),state%print())
128+
if (allocated(error)) return
129+
130+
end subroutine test_eigh_real_${ri}$
131+
132+
#:endif
133+
#:endfor
134+
135+
!> Simple complex matrix eigenvalues
136+
#:for ck,ct,ci in CMPLX_KINDS_TYPES
137+
#:if ck!="xdp"
138+
subroutine test_eig_complex_${ci}$(error)
139+
type(error_type), allocatable, intent(out) :: error
140+
141+
!> Reference solution
142+
real(${ck}$), parameter :: zero = 0.0_${ck}$
143+
real(${ck}$), parameter :: two = 2.0_${ck}$
144+
real(${ck}$), parameter :: sqrt2o2 = sqrt(two)*0.5_${ck}$
145+
real(${ck}$), parameter :: tol = sqrt(epsilon(zero))
146+
${ct}$, parameter :: cone = (1.0_${ck}$,0.0_${ck}$)
147+
${ct}$, parameter :: cimg = (0.0_${ck}$,1.0_${ck}$)
148+
${ct}$, parameter :: czero = (0.0_${ck}$,0.0_${ck}$)
149+
150+
!> Local vaciables
151+
type(linalg_state_type) :: state
152+
${ct}$ :: A(2,2),lambda(2),Avec(2,2),Ares(2,2),lres(2)
153+
154+
!> Matcix with real eigenvalues
155+
A = transpose(reshape([ cone, cimg, &
156+
-cimg, cone], [2,2]))
157+
158+
call eig(A,lambda,right=Avec,err=state)
159+
160+
!> Expected eigenvalues and eigenvectors
161+
lres(1) = two
162+
lres(2) = zero
163+
164+
!> Eigenvectors may vary: do not use for error
165+
Ares(1,1) = cmplx(zero,sqrt2o2,kind=${ck}$)
166+
Ares(1,2) = cmplx(sqrt2o2,zero,kind=${ck}$)
167+
Ares(2,1) = cmplx(sqrt2o2,zero,kind=${ck}$)
168+
Ares(2,2) = cmplx(zero,sqrt2o2,kind=${ck}$)
169+
170+
call check(error,state%ok(),state%print())
171+
if (allocated(error)) return
172+
173+
call check(error, all(abs(lambda-lres)<=tol), 'results match expected')
174+
if (allocated(error)) return
175+
176+
end subroutine test_eig_complex_${ci}$
177+
178+
#:endif
179+
#:endfor
180+
181+
182+
end module test_linalg_eigenvalues
183+
184+
program test_eigenvalues
185+
use, intrinsic :: iso_fortran_env, only : error_unit
186+
use testdrive, only : run_testsuite, new_testsuite, testsuite_type
187+
use test_linalg_eigenvalues, only : test_eig_eigh
188+
implicit none
189+
integer :: stat, is
190+
type(testsuite_type), allocatable :: testsuites(:)
191+
character(len=*), parameter :: fmt = '("#", *(1x, a))'
192+
193+
stat = 0
194+
195+
testsuites = [ &
196+
new_testsuite("linalg_eigenvalues", test_eig_eigh) &
197+
]
198+
199+
do is = 1, size(testsuites)
200+
write(error_unit, fmt) "Testing:", testsuites(is)%name
201+
call run_testsuite(testsuites(is)%collect, error_unit, stat)
202+
end do
203+
204+
if (stat > 0) then
205+
write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
206+
error stop
207+
end if
208+
end program test_eigenvalues

0 commit comments

Comments
 (0)