|
1 | 1 | #:include "common.fypp"
|
2 | 2 | #:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
|
3 |
| -#:set KINDS_TYPES = R_KINDS_TYPES |
| 3 | +#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) |
| 4 | +#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES |
4 | 5 | module test_specialmatrices
|
5 | 6 | use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test
|
6 | 7 | use stdlib_kinds
|
| 8 | + use stdlib_constants |
7 | 9 | use stdlib_linalg, only: hermitian
|
8 | 10 | use stdlib_linalg_state, only: linalg_state_type
|
9 | 11 | use stdlib_math, only: all_close
|
@@ -41,17 +43,33 @@ contains
|
41 | 43 | integer, parameter :: n = 5
|
42 | 44 | type(tridiagonal_${s1}$_type) :: A
|
43 | 45 | ${t1}$, allocatable :: Amat(:,:), dl(:), dv(:), du(:)
|
| 46 | + #: if t1.startswith('complex') |
| 47 | + real(wp), allocatable :: data(:, :) |
| 48 | + #:endif |
44 | 49 | ${t1}$, allocatable :: x(:)
|
45 | 50 | ${t1}$, allocatable :: y1(:), y2(:)
|
46 | 51 |
|
47 | 52 | ! Initialize matrix.
|
48 | 53 | allocate(dl(n-1), dv(n), du(n-1))
|
| 54 | + #:if t1.startswith('real') |
49 | 55 | call random_number(dl) ; call random_number(dv) ; call random_number(du)
|
| 56 | + #:else |
| 57 | + allocate(data(n, 2)) |
| 58 | + call random_number(data) ; dl%re = data(:n-1, 1) ; dl%im = data(:n-1, 2) |
| 59 | + call random_number(data) ; dv%re = data(:n, 1) ; dv%im = data(:n, 2) |
| 60 | + call random_number(data) ; du%re = data(:n-1, 1) ; du%im = data(:n-1, 2) |
| 61 | + #:endif |
50 | 62 | A = tridiagonal(dl, dv, du) ; Amat = dense(A)
|
51 | 63 |
|
52 | 64 | ! Random vectors.
|
| 65 | + #:if t1.startswith('real') |
53 | 66 | allocate(x(n), source = 0.0_wp) ; call random_number(x)
|
54 | 67 | allocate(y1(n), source = 0.0_wp) ; allocate(y2(n), source=0.0_wp)
|
| 68 | + #:else |
| 69 | + allocate(x(n), source=zero_c${k1}$) |
| 70 | + call random_number(data) ; x%re = data(:, 1) ; x%im = data(:, 2) |
| 71 | + allocate(y1(n), source = zero_c${k1}$) ; allocate(y2(n), source=zero_c${k1}$) |
| 72 | + #:endif |
55 | 73 |
|
56 | 74 | ! Test y = A @ x
|
57 | 75 | y1 = matmul(Amat, x) ; call spmv(A, x, y2)
|
@@ -117,33 +135,51 @@ contains
|
117 | 135 | integer, parameter :: n = 5
|
118 | 136 | type(symtridiagonal_${s1}$_type) :: A
|
119 | 137 | ${t1}$, allocatable :: Amat(:,:), dv(:), ev(:)
|
| 138 | + #:if t1.startswith('complex') |
| 139 | + real(wp), allocatable :: data(:, :) |
| 140 | + #:endif |
120 | 141 | ${t1}$, allocatable :: x(:)
|
121 | 142 | ${t1}$, allocatable :: y1(:), y2(:)
|
122 | 143 |
|
123 | 144 | ! Initialize matrix.
|
124 | 145 | allocate(ev(n-1), dv(n))
|
| 146 | + #:if t1.startswith('real') |
125 | 147 | call random_number(dv) ; call random_number(ev)
|
| 148 | + #:else |
| 149 | + allocate(data(n, 2), source=0.0_wp) |
| 150 | + call random_number(data) ; dv%re = data(:, 1) ; dv%im = data(:, 2) |
| 151 | + call random_number(data) ; ev%re = data(:n-1, 1) ; ev%im = data(:n-1, 2) |
| 152 | + #:endif |
126 | 153 | A = symtridiagonal(dv, ev) ; Amat = dense(A)
|
127 | 154 |
|
128 | 155 | ! Random vectors.
|
| 156 | + #:if t1.startswith('real') |
129 | 157 | allocate(x(n), source = 0.0_wp) ; call random_number(x)
|
130 | 158 | allocate(y1(n), source = 0.0_wp) ; allocate(y2(n), source=0.0_wp)
|
| 159 | + #:else |
| 160 | + allocate(x(n), source=zero_c${k1}$) |
| 161 | + call random_number(data) ; x%re = data(:, 1) ; x%im = data(:, 2) |
| 162 | + allocate(y1(n), source = zero_c${k1}$) ; allocate(y2(n), source=zero_c${k1}$) |
| 163 | + #:endif |
131 | 164 |
|
132 | 165 | ! Test y = A @ x
|
133 | 166 | y1 = matmul(Amat, x) ; call spmv(A, x, y2)
|
| 167 | + print *, "matvec :", all_close(y1, y2) |
134 | 168 | call check(error, all_close(y1, y2), .true.)
|
135 | 169 | if (allocated(error)) return
|
136 | 170 |
|
137 | 171 | ! Test y = A.T @ x
|
138 | 172 | y1 = 0.0_wp ; y2 = 0.0_wp
|
139 | 173 | y1 = matmul(transpose(Amat), x) ; call spmv(A, x, y2, op="T")
|
| 174 | + print *, "tranposed matvec :", all_close(y1, y2) |
140 | 175 | call check(error, all_close(y1, y2), .true.)
|
141 | 176 | if (allocated(error)) return
|
142 | 177 |
|
143 | 178 | #:if t1.startswith('complex')
|
144 | 179 | ! Test y = A.H @ x
|
145 | 180 | y1 = 0.0_wp ; y2 = 0.0_wp
|
146 | 181 | y1 = matmul(hermitian(Amat), x) ; call spmv(A, x, y2, op="H")
|
| 182 | + print *, "hermitian matvec :", all_close(y1, y2) |
147 | 183 | call check(error, all_close(y1, y2), .true.)
|
148 | 184 | if (allocated(error)) return
|
149 | 185 | #:endif
|
|
0 commit comments