Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 12 : PROGRAM quadratic_pseudospectrum_unittest
8 2 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
9 : #if defined(__HAS_IEEE_EXCEPTIONS)
10 : USE ieee_exceptions, ONLY: IEEE_ALL, ieee_get_halting_mode, ieee_set_halting_mode
11 : #endif
12 : USE kinds, ONLY: dp
13 : USE libcp2k, ONLY: cp2k_init, cp2k_finalize
14 : USE quadratic_pseudospectrum, ONLY: quadratic_dense, &
15 : quadratic_result_type
16 :
17 : IMPLICIT NONE
18 : COMPLEX(dp) :: op(2, 2, 4), s(2, 2), change(2, 2), rotated(2, 2, 4), metric(2, 2)
19 2 : COMPLEX(dp), ALLOCATABLE :: defect(:, :, :), overlap(:, :), projected(:, :)
20 : REAL(dp) :: shift(4), weights(4), expected, center
21 2 : TYPE(quadratic_result_type) :: dense, iterative, transformed
22 : INTEGER :: status, i, n, io_unit
23 :
24 2 : CALL cp2k_init()
25 2 : io_unit = cp_logger_get_default_io_unit()
26 2 : s = CMPLX(0.0_dp, 0.0_dp, dp)
27 2 : s(1, 1) = 1.0_dp
28 2 : s(2, 2) = 1.0_dp
29 2 : op = CMPLX(0.0_dp, 0.0_dp, dp)
30 2 : op(1, 1, 1) = 1.0_dp
31 2 : op(2, 2, 1) = -1.0_dp
32 2 : op(1, 2, 2) = CMPLX(0.0_dp, 1.0_dp, dp)
33 2 : op(2, 1, 2) = CMPLX(0.0_dp, -1.0_dp, dp)
34 2 : shift = [0.3_dp, 0.2_dp, 0.0_dp, 0.0_dp]
35 2 : weights = [1.0_dp, 0.7_dp, 0.7_dp, 0.7_dp]
36 : expected = 1.0_dp + shift(1)**2 + weights(2)**2*(1.0_dp + shift(2)**2) &
37 2 : - 2.0_dp*SQRT(shift(1)**2 + weights(2)**4*shift(2)**2)
38 2 : CALL quadratic_dense(op, s, shift, weights, 1, .FALSE., 2, 200, 1.e-10_dp, 1.e-10_dp, dense, status)
39 2 : IF (status /= 0) ERROR STOP 'Quadratic Pauli oracle failed'
40 2 : IF (ABS(dense%eigenvalues(1) - expected) > 1.e-12_dp) ERROR STOP 'Analytic quadratic eigenvalue mismatch'
41 2 : CALL quadratic_dense(op, s, shift, weights, 1, .TRUE., 2, 200, 1.e-10_dp, 1.e-10_dp, iterative, status)
42 2 : IF (status /= 0) ERROR STOP 'Quadratic Pauli iteration failed'
43 4 : IF (MAXVAL(ABS(dense%eigenvalues - iterative%eigenvalues)) > 1.e-12_dp) ERROR STOP 'Pauli Ritz mismatch'
44 : change = RESHAPE([CMPLX(1.2_dp, 0.1_dp, dp), CMPLX(0.15_dp, -0.1_dp, dp), &
45 2 : CMPLX(0.1_dp, 0.2_dp, dp), CMPLX(0.8_dp, 0.0_dp, dp)], [2, 2])
46 30 : metric = MATMUL(CONJG(TRANSPOSE(change)), change)
47 10 : DO i = 1, 4
48 282 : rotated(:, :, i) = MATMUL(CONJG(TRANSPOSE(change)), MATMUL(op(:, :, i), change))
49 : END DO
50 : CALL quadratic_dense(rotated, metric, shift, weights, 1, .FALSE., 2, 200, 1.e-10_dp, 1.e-10_dp, &
51 2 : transformed, status)
52 2 : IF (status /= 0) ERROR STOP 'Nonorthogonal complex basis failed'
53 4 : IF (MAXVAL(ABS(dense%eigenvalues - transformed%eigenvalues)) > 1.e-12_dp) ERROR STOP 'Basis covariance failed'
54 12 : IF (MAXVAL(ABS(dense%errors - transformed%errors)) > 1.e-12_dp) ERROR STOP 'Observable residual covariance failed'
55 12 : IF (MAXVAL(ABS(dense%means - transformed%means)) > 1.e-12_dp) ERROR STOP 'Observable mean covariance failed'
56 10 : DO i = 1, 4
57 56 : op(:, :, i) = op(:, :, i) + REAL(i, dp)*s
58 10 : shift(i) = shift(i) + REAL(i, dp)
59 : END DO
60 2 : CALL quadratic_dense(op, s, shift, weights, 1, .FALSE., 2, 200, 1.e-10_dp, 1.e-10_dp, transformed, status)
61 2 : IF (status /= 0) ERROR STOP 'Shifted quadratic oracle failed'
62 4 : IF (MAXVAL(ABS(dense%eigenvalues - transformed%eigenvalues)) > 1.e-12_dp) ERROR STOP 'Origin invariance failed'
63 12 : IF (MAXVAL(ABS(dense%errors - transformed%errors)) > 1.e-12_dp) ERROR STOP 'Residual origin invariance failed'
64 :
65 2 : n = 40
66 2 : center = 20.0_dp
67 2 : ALLOCATE (defect(n, n, 4), overlap(n, n), projected(3, 3))
68 2 : defect = CMPLX(0.0_dp, 0.0_dp, dp)
69 2 : overlap = CMPLX(0.0_dp, 0.0_dp, dp)
70 82 : DO i = 1, n
71 80 : overlap(i, i) = 1.0_dp
72 80 : defect(i, i, 1) = 5.0_dp
73 80 : defect(i, i, 2) = REAL(i, dp) - center
74 82 : IF (i < n) THEN
75 78 : defect(i, i + 1, 1) = CMPLX(-0.1_dp, 0.02_dp, dp)
76 78 : defect(i + 1, i, 1) = CONJG(defect(i, i + 1, 1))
77 : END IF
78 : END DO
79 2 : defect(20, 20, 1) = 0.0_dp
80 2 : shift = 0.0_dp
81 2 : weights = [1.0_dp, 0.2_dp, 0.2_dp, 0.2_dp]
82 2 : CALL quadratic_dense(defect, overlap, shift, weights, 3, .FALSE., n, 500, 1.e-9_dp, 1.e-10_dp, dense, status)
83 2 : IF (status /= 0) ERROR STOP 'Defect oracle failed'
84 2 : IF (ABS(dense%states(20, 1))**2 < 0.95_dp) ERROR STOP 'Defect state not localized'
85 2 : CALL quadratic_dense(defect, overlap, shift, weights, 3, .TRUE., 20, 500, 1.e-9_dp, 1.e-10_dp, iterative, status)
86 2 : IF (status /= 0) ERROR STOP 'Defect iteration failed'
87 8 : IF (MAXVAL(ABS(dense%eigenvalues - iterative%eigenvalues)) > 1.e-8_dp) ERROR STOP 'Defect Ritz eigenvalue mismatch'
88 2906 : projected(:, :) = MATMUL(CONJG(TRANSPOSE(dense%states)), iterative%states)
89 26 : IF (SUM(ABS(projected)**2) < 3.0_dp - 1.e-7_dp) ERROR STOP 'Defect invariant subspace mismatch'
90 2 : CALL quadratic_dense(defect, overlap, shift, weights, 3, .TRUE., 12, 1, 1.e-12_dp, 1.e-10_dp, iterative, status)
91 2 : IF (status /= 1) ERROR STOP 'Unconverged quadratic solve accepted'
92 : ! Exactly degenerate joint states: compare the projector, not individual eigenvectors.
93 2 : defect = CMPLX(0.0_dp, 0.0_dp, dp)
94 82 : DO i = 1, n
95 82 : defect(i, i, 1) = 2.0_dp
96 : END DO
97 2 : defect(19, 19, 1) = 0.0_dp
98 2 : defect(20, 20, 1) = 0.0_dp
99 2 : CALL quadratic_dense(defect, overlap, shift, weights, 2, .TRUE., 8, 200, 1.e-10_dp, 1.e-10_dp, iterative, status)
100 2 : IF (status /= 0) ERROR STOP 'Degenerate quadratic iteration failed'
101 6 : IF (MAXVAL(ABS(iterative%eigenvalues)) > 1.e-12_dp) ERROR STOP 'Degenerate eigenvalues failed'
102 14 : IF (ABS(SUM(ABS(iterative%states(19:20, :))**2) - 2.0_dp) > 1.e-12_dp) ERROR STOP 'Degenerate projector failed'
103 2 : overlap(1, 1) = -1.0_dp
104 2 : CALL quadratic_dense(defect, overlap, shift, weights, 3, .FALSE., n, 200, 1.e-9_dp, 1.e-10_dp, dense, status)
105 2 : IF (status /= -2) ERROR STOP 'Indefinite metric accepted'
106 2 : overlap(1, 1) = 1.e-14_dp
107 2 : CALL quadratic_dense(defect, overlap, shift, weights, 3, .FALSE., n, 200, 1.e-9_dp, 1.e-10_dp, dense, status)
108 2 : IF (status /= -2) ERROR STOP 'Ill-conditioned metric accepted'
109 2 : DEALLOCATE (defect, overlap, projected)
110 2 : CALL periodic_tests()
111 2 : CALL translation_tests()
112 2 : CALL folded_translation_tests()
113 2 : CALL compressed_translation_tests()
114 2 : CALL cp2k_finalize()
115 : CONTAINS
116 : ! **************************************************************************************************
117 : !> \brief Compression of a three-state unitary to two AOs retains the missing normal Gram exactly.
118 : ! **************************************************************************************************
119 2 : SUBROUTINE compressed_translation_tests()
120 : COMPLEX(dp) :: change(2, 2), changed(2, 2, 2), &
121 : metric(2, 2), op(2, 2, 2), query(2), &
122 : transformed(2, 2)
123 : INTEGER :: i, status
124 : REAL(dp) :: center, diagonal(2), err, expected(2), &
125 : split, weights(2)
126 2 : TYPE(quadratic_result_type) :: full, projected, rotated
127 :
128 2 : metric = 0
129 2 : metric(1, 1) = 1
130 2 : metric(2, 2) = 1
131 2 : op = 0
132 2 : op(1, 1, 1) = 0.2_dp
133 2 : op(2, 2, 1) = -0.4_dp
134 2 : op(1, 2, 2) = 1
135 2 : query = [CMPLX(0.1_dp, 0.0_dp, dp), EXP(CMPLX(0.0_dp, 0.73_dp, dp))]
136 2 : weights = [1.0_dp, 0.7_dp]
137 6 : diagonal = [0.1_dp**2, 0.5_dp**2] + 2*weights(2)**2
138 6 : center = SUM(diagonal)/2
139 2 : split = SQRT(((diagonal(1) - diagonal(2))/2)**2 + weights(2)**4)
140 6 : expected = [center - split, center + split]
141 : CALL quadratic_dense(op, metric, query, weights, 2, .FALSE., 2, 100, 1.e-11_dp, 1.e-10_dp, &
142 2 : full, status, unitary=[.FALSE., .TRUE.])
143 2 : IF (status /= 0) ERROR STOP 'Compressed unitary solve failed'
144 6 : err = MAXVAL(ABS(full%eigenvalues - expected))
145 2 : IF (err > 1.e-12_dp) ERROR STOP 'Full translation normal Gram differs from analytic result'
146 14 : IF (MAXVAL(ABS(full%errors**2 - full%projected_errors**2 - full%leakage**2)) > 1.e-12_dp) THEN
147 0 : ERROR STOP 'Continuum/projected residual decomposition failed'
148 : END IF
149 6 : IF (MAXVAL(ABS(full%leakage(2, :)**2 - ABS(full%states(1, :))**2)) > 1.e-12_dp) THEN
150 0 : ERROR STOP 'Translation leakage differs from discarded ambient component'
151 : END IF
152 2 : CALL quadratic_dense(op, metric, query, weights, 2, .FALSE., 2, 100, 1.e-11_dp, 1.e-10_dp, projected, status)
153 2 : IF (status /= 0) ERROR STOP 'Projected unitary comparison failed'
154 6 : IF (MAXVAL(ABS(projected%eigenvalues - full%eigenvalues)) < 0.1_dp) ERROR STOP 'Projection negative control failed'
155 : change = RESHAPE([CMPLX(1.2_dp, 0.1_dp, dp), CMPLX(0.15_dp, -0.1_dp, dp), &
156 2 : CMPLX(0.1_dp, 0.2_dp, dp), CMPLX(0.8_dp, 0.0_dp, dp)], [2, 2])
157 30 : transformed = MATMUL(CONJG(TRANSPOSE(change)), change)
158 6 : DO i = 1, 2
159 142 : changed(:, :, i) = MATMUL(CONJG(TRANSPOSE(change)), MATMUL(op(:, :, i), change))
160 : END DO
161 : CALL quadratic_dense(changed, transformed, query, weights, 2, .TRUE., 2, 100, 1.e-11_dp, 1.e-10_dp, &
162 2 : rotated, status, unitary=[.FALSE., .TRUE.])
163 2 : IF (status /= 0) ERROR STOP 'Compressed unitary metric iteration failed'
164 6 : IF (MAXVAL(ABS(rotated%eigenvalues - full%eigenvalues)) > 1.e-12_dp) ERROR STOP 'Unitary metric covariance failed'
165 14 : IF (MAXVAL(ABS(rotated%leakage - full%leakage)) > 1.e-12_dp) ERROR STOP 'Leakage metric covariance failed'
166 14 : op(:, :, 2) = 2*metric
167 : CALL quadratic_dense(op, metric, query, weights, 2, .FALSE., 2, 100, 1.e-11_dp, 1.e-10_dp, &
168 2 : rotated, status, unitary=[.FALSE., .TRUE.])
169 2 : IF (status /= -4) ERROR STOP 'Invalid expanding unitary compression accepted'
170 2 : IF (io_unit > 0) WRITE (io_unit, '(A,ES15.6)') 'Compressed unitary full-normal analytic error: ', err
171 4 : END SUBROUTINE compressed_translation_tests
172 :
173 : ! **************************************************************************************************
174 : !> \brief A trimerized periodic chain must equal its independently assembled three-band Bloch blocks.
175 : ! **************************************************************************************************
176 2 : SUBROUTINE folded_translation_tests()
177 : INTEGER, PARAMETER :: cells = 6, n = 3*cells
178 :
179 : COMPLEX(dp) :: bloch(3, 3, 2), identity(3, 3), &
180 : metric(n, n), operators(n, n, 2), &
181 : phase, query(2)
182 : INTEGER :: first, i, ie, ik, iq, j, loc, status
183 : REAL(dp) :: energy(5), expected(n), hopping, &
184 : joined(n), k, maximum_error, q, &
185 : weights(2)
186 2 : TYPE(quadratic_result_type) :: full, sector
187 :
188 2 : operators = 0
189 2 : metric = 0
190 2 : identity = 0
191 8 : DO i = 1, 3
192 8 : identity(i, i) = 1
193 : END DO
194 38 : DO i = 1, n
195 36 : metric(i, i) = 1
196 36 : j = MODULO(i, n) + 1
197 36 : operators(i, j, 2) = 1
198 36 : hopping = 0.8_dp
199 36 : IF (MODULO(i, 3) == 0) hopping = 1.2_dp
200 36 : operators(i, j, 1) = -hopping
201 38 : operators(j, i, 1) = -hopping
202 : END DO
203 2 : energy = [-2.2_dp, -1.0_dp, 0.0_dp, 0.7_dp, 2.0_dp]
204 2 : weights = [1.0_dp, 0.7_dp]
205 2 : maximum_error = 0
206 28 : DO ik = 0, 12
207 26 : k = -ACOS(-1.0_dp) + 2*ACOS(-1.0_dp)*REAL(ik, dp)/13
208 158 : DO ie = 1, SIZE(energy)
209 390 : query = [CMPLX(energy(ie), 0.0_dp, dp), EXP(CMPLX(0.0_dp, k, dp))]
210 : CALL quadratic_dense(operators, metric, query, weights, n, .FALSE., n, 100, &
211 130 : 1.e-10_dp, 1.e-10_dp, full, status)
212 130 : IF (status /= 0) ERROR STOP 'Full trimerized quadratic scan failed'
213 910 : DO iq = 0, cells - 1
214 780 : q = 2*ACOS(-1.0_dp)*REAL(iq, dp)/cells
215 780 : phase = EXP(CMPLX(0.0_dp, q, dp))
216 780 : bloch = 0
217 780 : bloch(1, 2, 2) = 1
218 780 : bloch(2, 3, 2) = 1
219 780 : bloch(3, 1, 2) = phase
220 780 : bloch(1, 2, 1) = -0.8_dp
221 780 : bloch(2, 3, 1) = -0.8_dp
222 780 : bloch(3, 1, 1) = -1.2_dp*phase
223 19500 : bloch(:, :, 1) = bloch(:, :, 1) + CONJG(TRANSPOSE(bloch(:, :, 1)))
224 : CALL quadratic_dense(bloch, identity, query, weights, 3, .FALSE., 3, 100, &
225 780 : 1.e-10_dp, 1.e-10_dp, sector, status)
226 780 : IF (status /= 0) ERROR STOP 'Bloch-block quadratic scan failed'
227 780 : first = 3*iq + 1
228 3250 : joined(first:first + 2) = sector%eigenvalues
229 : END DO
230 2470 : DO i = 1, n
231 44460 : loc = MINLOC(joined, DIM=1)
232 2340 : expected(i) = joined(loc)
233 2470 : joined(loc) = HUGE(1.0_dp)
234 : END DO
235 2496 : maximum_error = MAX(maximum_error, MAXVAL(ABS(full%eigenvalues - expected)))
236 : END DO
237 : END DO
238 2 : IF (maximum_error > 1.e-11_dp) ERROR STOP 'Full/Bloch trimerized quadratic spectra differ'
239 3 : IF (io_unit > 0) WRITE (io_unit, '(A,I0,A,ES15.6)') 'Trimerized full/Bloch spectra: ', 13*SIZE(energy)*n, &
240 2 : ' eigenvalues; maximum error ', maximum_error
241 2 : END SUBROUTINE folded_translation_tests
242 :
243 : ! **************************************************************************************************
244 : !> \brief Exact Bloch phases, nonnormal open translations and an independent generalized Gram oracle.
245 : ! **************************************************************************************************
246 2 : SUBROUTINE translation_tests()
247 : INTEGER, PARAMETER :: n = 18
248 :
249 : COMPLEX(dp) :: changed(n, n, 2), identity(n, n), metric(n, n), offset(2), &
250 : operators(n, n, 2), projector(n, n), query(2), split(n, n, 3), transform(n, n)
251 : INTEGER :: i, j, status
252 : REAL(dp) :: covariance_error, expected(n), &
253 : gap_error, hopping, mean_error, &
254 : nonnormal_error, phase, weights(2)
255 2 : TYPE(quadratic_result_type) :: altered, incorrect, iterative, original
256 :
257 2 : identity = 0
258 2 : operators = 0
259 38 : DO i = 1, n
260 36 : identity(i, i) = 1
261 38 : operators(i, MODULO(i, n) + 1, 2) = 1
262 : END DO
263 1370 : operators(:, :, 1) = -operators(:, :, 2) - CONJG(TRANSPOSE(operators(:, :, 2)))
264 2 : phase = 2*ACOS(-1.0_dp)*5/REAL(n, dp)
265 6 : query = [CMPLX(-2*COS(phase) + 0.031_dp, 0.0_dp, dp), EXP(CMPLX(0.0_dp, phase, dp))]
266 2 : weights = [1.0_dp, 1.3_dp]
267 : CALL quadratic_dense(operators, identity, query, weights, 1, .FALSE., n, 200, &
268 2 : 1.e-10_dp, 1.e-10_dp, original, status)
269 2 : IF (status /= 0) ERROR STOP 'Unitary translation quadratic solve failed'
270 2 : IF (ABS(original%eigenvalues(1) - 0.031_dp**2) > 1.e-12_dp) ERROR STOP 'Exact Bloch energy/phase failed'
271 2 : IF (ABS(original%complex_means(2, 1) - query(2)) > 1.e-12_dp) ERROR STOP 'Bloch translation mean lost phase'
272 2 : IF (original%errors(2, 1) > 1.e-12_dp) ERROR STOP 'Bloch translation residual failed'
273 : CALL quadratic_dense(operators, identity, REAL(query, dp), weights, 1, .FALSE., n, 200, &
274 6 : 1.e-10_dp, 1.e-10_dp, altered, status)
275 2 : IF (status /= -1) ERROR STOP 'Hermitian interface accepted non-Hermitian translation'
276 :
277 : ! The truncated translation is nonnormal. Its adjoint cannot be replaced by itself
278 : ! or by an independent sum of squares of its Hermitian real/imaginary parts.
279 2 : operators(n, 1, 2) = 0
280 686 : operators(:, :, 1) = 0
281 36 : DO i = 1, n - 1
282 34 : hopping = 0.8_dp
283 34 : IF (MODULO(i, 3) == 0) hopping = 1.2_dp
284 34 : operators(i, i + 1, 1) = -hopping
285 36 : operators(i + 1, i, 1) = -hopping
286 : END DO
287 2 : operators(4, 4, 1) = 0.17_dp
288 2 : query = [CMPLX(0.33_dp, 0.0_dp, dp), EXP(CMPLX(0.0_dp, 0.73_dp, dp))]
289 : CALL quadratic_dense(operators, identity, query, weights, n, .FALSE., n, 200, &
290 2 : 1.e-10_dp, 1.e-10_dp, original, status)
291 2 : IF (status /= 0) ERROR STOP 'Open trimerized translation solve failed'
292 2 : CALL gram_reference(operators, identity, query, weights, expected)
293 38 : gap_error = MAXVAL(ABS(original%eigenvalues - expected))
294 2 : IF (gap_error > 1.e-11_dp) ERROR STOP 'General quadratic Gram oracle mismatch'
295 : CALL quadratic_dense(operators, identity, query, weights, 2, .TRUE., 10, 2000, &
296 2 : 1.e-10_dp, 1.e-10_dp, iterative, status)
297 2 : IF (status /= 0) ERROR STOP 'Nonnormal restarted iteration failed'
298 6 : IF (MAXVAL(ABS(iterative%eigenvalues - expected(:2))) > 1.e-10_dp) ERROR STOP 'Nonnormal Ritz spectrum mismatch'
299 2 : projector = MATMUL(iterative%states, CONJG(TRANSPOSE(iterative%states))) - &
300 6094 : MATMUL(original%states(:, :2), CONJG(TRANSPOSE(original%states(:, :2))))
301 686 : IF (MAXVAL(ABS(projector)) > 1.e-8_dp) ERROR STOP 'Nonnormal invariant-subspace mismatch'
302 686 : split(:, :, 1) = operators(:, :, 1)
303 686 : split(:, :, 2) = 0.5_dp*(operators(:, :, 2) + CONJG(TRANSPOSE(operators(:, :, 2))))
304 686 : split(:, :, 3) = CMPLX(0.0_dp, -0.5_dp, dp)*(operators(:, :, 2) - CONJG(TRANSPOSE(operators(:, :, 2))))
305 : CALL quadratic_dense(split, identity, [REAL(query(1), dp), REAL(query(2), dp), AIMAG(query(2))], &
306 : [weights(1), weights(2), weights(2)], n, .FALSE., n, 200, &
307 14 : 1.e-10_dp, 1.e-10_dp, incorrect, status)
308 2 : IF (status /= 0) ERROR STOP 'Hermitian-split comparison solve failed'
309 38 : nonnormal_error = MAXVAL(ABS(original%eigenvalues - incorrect%eigenvalues))
310 2 : IF (nonnormal_error < 0.01_dp) ERROR STOP 'Open-translation counterexample lost nonnormality'
311 :
312 2 : transform = identity
313 36 : DO i = 1, n - 1
314 34 : transform(i, i) = 1 + 0.02_dp*i
315 36 : transform(i, i + 1) = CMPLX(0.21_dp, -0.12_dp, dp)
316 : END DO
317 12350 : metric = MATMUL(CONJG(TRANSPOSE(transform)), transform)
318 6 : DO j = 1, 2
319 50766 : changed(:, :, j) = MATMUL(CONJG(TRANSPOSE(transform)), MATMUL(operators(:, :, j), transform))
320 : END DO
321 : CALL quadratic_dense(changed, metric, query, weights, n, .FALSE., n, 200, &
322 2 : 1.e-10_dp, 1.e-10_dp, altered, status)
323 2 : IF (status /= 0) ERROR STOP 'Nonnormal complex AO basis failed'
324 2 : CALL gram_reference(changed, metric, query, weights, expected)
325 38 : covariance_error = MAXVAL(ABS(original%eigenvalues - altered%eigenvalues))
326 2 : IF (covariance_error > 1.e-11_dp) ERROR STOP 'Nonnormal metric covariance failed'
327 38 : IF (MAXVAL(ABS(altered%eigenvalues - expected)) > 1.e-11_dp) ERROR STOP 'Nonorthogonal Gram oracle mismatch'
328 110 : mean_error = MAXVAL(ABS(original%complex_means - altered%complex_means))
329 2 : IF (mean_error > 1.e-10_dp) ERROR STOP 'Complex expectation covariance failed'
330 110 : IF (MAXVAL(ABS(original%errors - altered%errors)) > 1.e-10_dp) ERROR STOP 'Complex residual covariance failed'
331 2 : offset = [CMPLX(0.21_dp, -0.12_dp, dp), CMPLX(-0.43_dp, 0.31_dp, dp)]
332 6 : DO j = 1, 2
333 1374 : changed(:, :, j) = changed(:, :, j) + offset(j)*metric
334 : END DO
335 : CALL quadratic_dense(changed, metric, query + offset, weights, n, .FALSE., n, 200, &
336 6 : 1.e-10_dp, 1.e-10_dp, altered, status)
337 2 : IF (status /= 0) ERROR STOP 'Complex operator/query origin solve failed'
338 38 : IF (MAXVAL(ABS(original%eigenvalues - altered%eigenvalues)) > 1.e-11_dp) ERROR STOP 'Complex origin covariance failed'
339 6 : DO j = 1, 2
340 78 : IF (MAXVAL(ABS(altered%complex_means(j, :) - original%complex_means(j, :) - offset(j))) > 1.e-10_dp) THEN
341 0 : ERROR STOP 'Complex expectation origin covariance failed'
342 : END IF
343 : END DO
344 3 : IF (io_unit > 0) WRITE (io_unit, '(A,4ES15.6)') 'Translation Gram/metric/complex-mean/split-error: ', &
345 2 : gap_error, covariance_error, mean_error, nonnormal_error
346 25382 : changed(:, :, 1) = MATMUL(CONJG(TRANSPOSE(transform)), MATMUL(operators(:, :, 1), transform))
347 686 : changed(:, :, 2) = metric
348 2 : query(2) = 1
349 : CALL quadratic_dense(changed, metric, query, weights, 2, .TRUE., 10, 2000, &
350 2 : 1.e-10_dp, 1.e-10_dp, iterative, status, unitary=[.FALSE., .TRUE.])
351 2 : IF (status /= 0) ERROR STOP 'Restarted identity translation failed'
352 2 : IF (iterative%iterations < 2) ERROR STOP 'Identity translation did not exercise restarts'
353 6 : IF (MAXVAL(iterative%errors(2, :)) > 5.e-8_dp) ERROR STOP 'Identity full residual exceeds roundoff'
354 6 : IF (MAXVAL(iterative%projected_errors(2, :)) > 1.e-12_dp) ERROR STOP 'Identity projected residual drift'
355 6 : IF (MAXVAL(iterative%leakage(2, :)) > 5.e-8_dp) ERROR STOP 'Identity projection leakage exceeds roundoff'
356 2 : END SUBROUTINE translation_tests
357 :
358 : ! **************************************************************************************************
359 : !> \brief Independent explicit normal-product reference using LU solves and a generalized eigensolve.
360 : !> \param operators Covariant operators
361 : !> \param metric Positive AO overlap
362 : !> \param query Complex operator queries
363 : !> \param weights Nonnegative residual scales
364 : !> \param values Complete ascending generalized spectrum
365 : ! **************************************************************************************************
366 4 : SUBROUTINE gram_reference(operators, metric, query, weights, values)
367 : COMPLEX(dp), INTENT(IN) :: operators(:, :, :), metric(:, :), query(:)
368 : REAL(dp), INTENT(IN) :: weights(:)
369 : REAL(dp), INTENT(OUT) :: values(:)
370 4 : COMPLEX(dp), ALLOCATABLE :: gram(:, :), shifted(:, :), solved(:, :), lu(:, :), work(:)
371 4 : REAL(dp), ALLOCATABLE :: rwork(:)
372 4 : INTEGER, ALLOCATABLE :: pivots(:)
373 : INTEGER :: n, i, status
374 : #if defined(__HAS_IEEE_EXCEPTIONS)
375 : LOGICAL :: halt(SIZE(IEEE_ALL))
376 : #endif
377 :
378 4 : n = SIZE(metric, 1)
379 64 : ALLOCATE (gram(n, n), shifted(n, n), solved(n, n), lu(n, n), pivots(n), work(4*n), rwork(4*n))
380 4 : gram = 0
381 12 : DO i = 1, SIZE(query)
382 2744 : lu(:, :) = metric
383 2744 : shifted(:, :) = operators(:, :, i) - query(i)*metric
384 2744 : solved(:, :) = shifted
385 : #if defined(__HAS_IEEE_EXCEPTIONS)
386 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
387 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
388 : #endif
389 8 : CALL zgesv(n, n, lu, n, pivots, solved, n, status)
390 : #if defined(__HAS_IEEE_EXCEPTIONS)
391 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
392 : #endif
393 8 : IF (status /= 0) ERROR STOP 'Explicit Gram metric solve failed'
394 52164 : gram(:, :) = gram + weights(i)**2*MATMUL(CONJG(TRANSPOSE(shifted)), solved)
395 : END DO
396 1372 : IF (MAXVAL(ABS(gram - CONJG(TRANSPOSE(gram)))) > 1.e-11_dp) ERROR STOP 'Explicit Gram is not Hermitian'
397 1372 : lu(:, :) = metric
398 : #if defined(__HAS_IEEE_EXCEPTIONS)
399 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
400 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
401 : #endif
402 4 : CALL zhegv(1, 'N', 'L', n, gram, n, lu, n, values, work, SIZE(work), rwork, status)
403 : #if defined(__HAS_IEEE_EXCEPTIONS)
404 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
405 : #endif
406 4 : IF (status /= 0) ERROR STOP 'Explicit generalized Gram spectrum failed'
407 4 : END SUBROUTINE gram_reference
408 :
409 : ! **************************************************************************************************
410 : !> \brief Three-dimensional trigonometric coordinates, basis covariance and torus-query invariance.
411 : ! **************************************************************************************************
412 2 : SUBROUTINE periodic_tests()
413 : INTEGER, PARAMETER :: n = 8
414 :
415 : COMPLEX(dp) :: changed(n, n, 7), identity(n, n), &
416 : metric(n, n), operators(n, n, 7), &
417 : transform(n, n)
418 : INTEGER :: i, j, k, status
419 : REAL(dp) :: diagonal(n), expected(n), move(3), &
420 : radius(3), shift(7), theta(3), value, &
421 : weights(7)
422 2 : TYPE(quadratic_result_type) :: altered, iterative, original
423 :
424 2 : identity = 0
425 2 : operators = 0
426 2 : theta = [0.21_dp, -0.37_dp, 1.43_dp]
427 2 : radius = [0.7_dp, 1.2_dp, 0.9_dp]
428 2 : shift(1) = 0.73_dp
429 2 : weights(1) = 1
430 8 : DO j = 1, 3
431 18 : shift(2*j:2*j + 1) = [COS(theta(j)), SIN(theta(j))]
432 20 : weights(2*j:2*j + 1) = radius(j)
433 : END DO
434 18 : DO i = 1, n
435 16 : identity(i, i) = 1
436 16 : operators(i, i, 1) = 0.2_dp*i
437 64 : DO j = 1, 3
438 48 : value = 2*ACOS(-1.0_dp)*MODULO(REAL(i*(2*j - 1), dp)/REAL(n + j, dp), 1.0_dp)
439 48 : operators(i, i, 2*j) = COS(value)
440 64 : operators(i, i, 2*j + 1) = SIN(value)
441 : END DO
442 130 : diagonal(i) = SUM((weights*ABS(operators(i, i, :) - shift))**2)
443 : END DO
444 18 : DO i = 1, n
445 144 : k = MINLOC(diagonal, DIM=1)
446 16 : expected(i) = diagonal(k)
447 18 : diagonal(k) = HUGE(1.0_dp)
448 : END DO
449 2 : CALL quadratic_dense(operators, identity, shift, weights, n, .FALSE., n, 100, 1.e-10_dp, 1.e-10_dp, original, status)
450 2 : IF (status /= 0) ERROR STOP 'Periodic quadratic diagonal oracle failed'
451 18 : IF (MAXVAL(ABS(original%eigenvalues - expected)) > 1.e-12_dp) ERROR STOP 'Seven-observable analytic spectrum failed'
452 16 : DO i = 1, n - 1
453 14 : operators(i, i + 1, 1) = CMPLX(-0.12_dp, 0.04_dp, dp)
454 16 : operators(i + 1, i, 1) = CONJG(operators(i, i + 1, 1))
455 : END DO
456 : ! The six initial vectors leave only one or two slots for the three-vector residual block.
457 6 : DO k = n - 1, n
458 : CALL quadratic_dense(operators(:k, :k, :), identity(:k, :k), shift, weights, 3, .FALSE., k, 100, &
459 4 : 1.e-10_dp, 1.e-10_dp, original, status)
460 4 : IF (status /= 0) ERROR STOP 'Partial residual block reference failed'
461 : CALL quadratic_dense(operators(:k, :k, :), identity(:k, :k), shift, weights, 3, .TRUE., k, 2, &
462 4 : 1.e-10_dp, 1.e-10_dp, iterative, status)
463 4 : IF (status /= 0) ERROR STOP 'Partial residual block did not fill the available space'
464 4 : IF (iterative%iterations /= 2) ERROR STOP 'Partial residual block iteration count mismatch'
465 18 : IF (MAXVAL(ABS(original%eigenvalues - iterative%eigenvalues)) > 1.e-12_dp) THEN
466 0 : ERROR STOP 'Partial residual block eigenvalue mismatch'
467 : END IF
468 : END DO
469 : ! Reuse the full-size reference from the final partial-residual-block check.
470 2 : transform = identity
471 16 : DO i = 1, n - 1
472 16 : transform(i, i + 1) = CMPLX(0.13_dp, -0.08_dp, dp)
473 : END DO
474 1170 : metric = MATMUL(CONJG(TRANSPOSE(transform)), transform)
475 16 : DO j = 1, 7
476 17376 : changed(:, :, j) = MATMUL(CONJG(TRANSPOSE(transform)), MATMUL(operators(:, :, j), transform))
477 : END DO
478 2 : CALL quadratic_dense(changed, metric, shift, weights, 3, .FALSE., n, 100, 1.e-10_dp, 1.e-10_dp, altered, status)
479 2 : IF (status /= 0) ERROR STOP 'Periodic complex metric failed'
480 8 : IF (MAXVAL(ABS(original%eigenvalues - altered%eigenvalues)) > 1.e-12_dp) ERROR STOP 'Periodic metric covariance failed'
481 2 : CALL quadratic_dense(changed, metric, shift, weights, 3, .TRUE., n, 2, 1.e-10_dp, 1.e-10_dp, iterative, status)
482 2 : IF (status /= 0) ERROR STOP 'Partial residual block with complex metric failed'
483 8 : IF (MAXVAL(ABS(altered%eigenvalues - iterative%eigenvalues)) > 1.e-12_dp) THEN
484 0 : ERROR STOP 'Partial residual block metric covariance failed'
485 : END IF
486 2 : move = [0.43_dp, -0.67_dp, 0.31_dp]
487 2 : changed = operators
488 8 : DO j = 1, 3
489 438 : changed(:, :, 2*j) = COS(move(j))*operators(:, :, 2*j) - SIN(move(j))*operators(:, :, 2*j + 1)
490 438 : changed(:, :, 2*j + 1) = SIN(move(j))*operators(:, :, 2*j) + COS(move(j))*operators(:, :, 2*j + 1)
491 20 : shift(2*j:2*j + 1) = [COS(theta(j) + move(j)), SIN(theta(j) + move(j))]
492 : END DO
493 2 : CALL quadratic_dense(changed, identity, shift, weights, 3, .FALSE., n, 100, 1.e-10_dp, 1.e-10_dp, altered, status)
494 2 : IF (status /= 0) ERROR STOP 'Shifted periodic coordinates failed'
495 8 : IF (MAXVAL(ABS(original%eigenvalues - altered%eigenvalues)) > 1.e-12_dp) ERROR STOP 'Periodic coordinate-origin covariance failed'
496 8 : DO j = 1, 3
497 6 : value = theta(j) + 2*ACOS(-1.0_dp)*(j - 2)
498 20 : shift(2*j:2*j + 1) = [COS(value), SIN(value)]
499 : END DO
500 2 : CALL quadratic_dense(operators, identity, shift, weights, 3, .FALSE., n, 100, 1.e-10_dp, 1.e-10_dp, altered, status)
501 2 : IF (status /= 0) ERROR STOP 'Periodic query wrapping failed'
502 8 : IF (MAXVAL(ABS(original%eigenvalues - altered%eigenvalues)) > 1.e-12_dp) ERROR STOP 'Periodic query invariance failed'
503 2 : CALL quadratic_dense(operators, identity, shift, weights(:6), 1, .FALSE., n, 100, 1.e-10_dp, 1.e-10_dp, altered, status)
504 2 : IF (status /= -1) ERROR STOP 'Mismatched observable weights accepted'
505 : CALL quadratic_dense(operators(:, :, :0), identity, shift(:0), weights(:0), 1, .FALSE., n, 100, &
506 2 : 1.e-10_dp, 1.e-10_dp, altered, status)
507 2 : IF (status /= -1) ERROR STOP 'Empty observable tuple accepted'
508 2 : IF (io_unit > 0) WRITE (io_unit, *) 'Periodic quadratic tuple, metric, iteration and origin checks passed.'
509 2 : END SUBROUTINE periodic_tests
510 : END PROGRAM quadratic_pseudospectrum_unittest
|
