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 :
8 : ! **************************************************************************************************
9 : !> \brief Metric-aware quadratic pseudospectra for finite AO operators.
10 : !> \author Thomas D. Kuehne, tkuehne@cp2k.org
11 : ! **************************************************************************************************
12 : MODULE quadratic_pseudospectrum
13 : USE ieee_arithmetic, ONLY: ieee_is_finite
14 : #if defined(__HAS_IEEE_EXCEPTIONS)
15 : USE ieee_exceptions, ONLY: IEEE_ALL, ieee_get_halting_mode, ieee_set_halting_mode
16 : #endif
17 : USE kinds, ONLY: dp
18 : IMPLICIT NONE
19 : PRIVATE
20 : PUBLIC :: quadratic_dense, quadratic_solve, quadratic_result_type
21 :
22 : TYPE quadratic_result_type
23 : REAL(dp), ALLOCATABLE :: eigenvalues(:), residuals(:), errors(:, :), means(:, :)
24 : REAL(dp), ALLOCATABLE :: projected_errors(:, :), leakage(:, :)
25 : COMPLEX(dp), ALLOCATABLE :: states(:, :), complex_means(:, :)
26 : INTEGER :: iterations = 0
27 : END TYPE quadratic_result_type
28 :
29 : INTERFACE quadratic_dense
30 : MODULE PROCEDURE quadratic_dense_real, quadratic_dense_complex
31 : END INTERFACE
32 : INTERFACE quadratic_solve
33 : MODULE PROCEDURE quadratic_solve_real, quadratic_solve_complex
34 : END INTERFACE
35 :
36 : ABSTRACT INTERFACE
37 : ! **************************************************************************************************
38 : !> \brief Apply a covariant AO observable or the overlap to a block of coefficients.
39 : !> \param OPERATOR Zero for overlap, otherwise the requested observable index
40 : !> \param v Input coefficients
41 : !> \param y Covariant output vectors
42 : !> \param status Zero on success
43 : ! **************************************************************************************************
44 : SUBROUTINE operator_apply(OPERATOR, v, y, status)
45 : IMPORT dp
46 : INTEGER, INTENT(IN) :: operator
47 : COMPLEX(dp), INTENT(IN) :: v(:, :)
48 : COMPLEX(dp), INTENT(OUT) :: y(:, :)
49 : INTEGER, INTENT(OUT) :: status
50 : END SUBROUTINE operator_apply
51 : ! **************************************************************************************************
52 : !> \brief Solve the positive AO metric for a block of covariant right hand sides.
53 : !> \param v Right hand sides
54 : !> \param y Contravariant solutions
55 : !> \param status Zero on success
56 : ! **************************************************************************************************
57 : SUBROUTINE metric_solve(v, y, status)
58 : IMPORT dp
59 : COMPLEX(dp), INTENT(IN) :: v(:, :)
60 : COMPLEX(dp), INTENT(OUT) :: y(:, :)
61 : INTEGER, INTENT(OUT) :: status
62 : END SUBROUTINE metric_solve
63 : END INTERFACE
64 :
65 : TYPE, ABSTRACT :: quadratic_operator_type
66 : CONTAINS
67 : PROCEDURE(operator_method), DEFERRED :: apply
68 : PROCEDURE(operator_method), DEFERRED :: adjoint
69 : PROCEDURE(metric_method), DEFERRED :: solve
70 : END TYPE quadratic_operator_type
71 :
72 : ABSTRACT INTERFACE
73 : ! **************************************************************************************************
74 : !> \brief ...
75 : !> \param this ...
76 : !> \param OPERATOR ...
77 : !> \param v ...
78 : !> \param y ...
79 : !> \param status ...
80 : ! **************************************************************************************************
81 : SUBROUTINE operator_method(this, OPERATOR, v, y, status)
82 : IMPORT dp, quadratic_operator_type
83 : CLASS(quadratic_operator_type), INTENT(INOUT) :: this
84 : INTEGER, INTENT(IN) :: operator
85 : COMPLEX(dp), INTENT(IN) :: v(:, :)
86 : COMPLEX(dp), INTENT(OUT) :: y(:, :)
87 : INTEGER, INTENT(OUT) :: status
88 : END SUBROUTINE operator_method
89 :
90 : ! **************************************************************************************************
91 : !> \brief ...
92 : !> \param this ...
93 : !> \param v ...
94 : !> \param y ...
95 : !> \param status ...
96 : ! **************************************************************************************************
97 : SUBROUTINE metric_method(this, v, y, status)
98 : IMPORT dp, quadratic_operator_type
99 : CLASS(quadratic_operator_type), INTENT(INOUT) :: this
100 : COMPLEX(dp), INTENT(IN) :: v(:, :)
101 : COMPLEX(dp), INTENT(OUT) :: y(:, :)
102 : INTEGER, INTENT(OUT) :: status
103 : END SUBROUTINE metric_method
104 : END INTERFACE
105 :
106 : TYPE, EXTENDS(quadratic_operator_type) :: dense_operator_type
107 : COMPLEX(dp), POINTER :: operators(:, :, :) => NULL(), metric(:, :) => NULL()
108 : COMPLEX(dp), ALLOCATABLE :: factor(:, :)
109 : CONTAINS
110 : PROCEDURE :: apply => dense_apply
111 : PROCEDURE :: adjoint => dense_adjoint
112 : PROCEDURE :: solve => dense_solve
113 : END TYPE dense_operator_type
114 :
115 : TYPE, EXTENDS(quadratic_operator_type) :: callback_operator_type
116 : PROCEDURE(operator_apply), POINTER, NOPASS :: apply_callback => NULL()
117 : PROCEDURE(operator_apply), POINTER, NOPASS :: adjoint_callback => NULL()
118 : PROCEDURE(metric_solve), POINTER, NOPASS :: solve_callback => NULL()
119 : CONTAINS
120 : PROCEDURE :: apply => callback_apply
121 : PROCEDURE :: adjoint => callback_adjoint
122 : PROCEDURE :: solve => callback_solve
123 : END TYPE callback_operator_type
124 :
125 : CONTAINS
126 :
127 : ! **************************************************************************************************
128 : !> \brief Dense AO oracle, also providing a dense-metric reference for the iterative solver.
129 : !> \param operators Hermitian observables in a covariant AO basis
130 : !> \param metric Positive AO overlap
131 : !> \param shift Query value for each observable
132 : !> \param weights Nonnegative observable scales
133 : !> \param nstates Number of lowest eigenpairs
134 : !> \param iterative Use the matrix-free Ritz iteration rather than full diagonalization
135 : !> \param max_subspace Maximum trial-space dimension
136 : !> \param max_iter Maximum Ritz iterations
137 : !> \param tolerance Relative eigenpair residual threshold
138 : !> \param metric_tolerance Minimum overlap eigenvalue ratio
139 : !> \param RESULT Eigenpairs, separate observable residuals and expectations
140 : !> \param status 0 success, 1 unconverged, -1 input, -2 metric, -4 numerical failure
141 : ! **************************************************************************************************
142 53 : SUBROUTINE quadratic_dense_real(operators, metric, shift, weights, nstates, iterative, max_subspace, &
143 : max_iter, tolerance, metric_tolerance, RESULT, status)
144 : COMPLEX(dp), INTENT(IN) :: operators(:, :, :), metric(:, :)
145 : REAL(dp), INTENT(IN) :: shift(:), weights(:)
146 : INTEGER, INTENT(IN) :: nstates
147 : LOGICAL, INTENT(IN) :: iterative
148 : INTEGER, INTENT(IN) :: max_subspace, max_iter
149 : REAL(dp), INTENT(IN) :: tolerance, metric_tolerance
150 : TYPE(quadratic_result_type), INTENT(OUT) :: result
151 : INTEGER, INTENT(OUT) :: status
152 :
153 : INTEGER :: j
154 :
155 53 : status = -1
156 53 : IF (SIZE(operators, 1) < 1 .OR. SIZE(operators, 1) /= SIZE(operators, 2)) RETURN
157 313 : DO j = 1, SIZE(operators, 3)
158 313 : IF (.NOT. hermitian(operators(:, :, j))) RETURN
159 : END DO
160 : CALL quadratic_dense_complex(operators, metric, CMPLX(shift, 0.0_dp, dp), weights, nstates, iterative, &
161 309 : max_subspace, max_iter, tolerance, metric_tolerance, RESULT, status)
162 : END SUBROUTINE quadratic_dense_real
163 :
164 : ! **************************************************************************************************
165 : !> \brief Dense oracle for general complex queries and possibly non-Hermitian AO operators.
166 : !> \param operators Covariant AO operators, not necessarily normal or Hermitian
167 : !> \param metric Positive AO overlap
168 : !> \param shift Complex query value for each operator
169 : !> \param weights Nonnegative residual scales
170 : !> \param nstates Number of lowest eigenpairs
171 : !> \param iterative Whether to use the bounded Ritz iteration
172 : !> \param max_subspace Maximum trial-space dimension
173 : !> \param max_iter Maximum Ritz iterations
174 : !> \param tolerance Relative eigenpair residual threshold
175 : !> \param metric_tolerance Minimum overlap eigenvalue ratio
176 : !> \param RESULT Eigenpairs, residuals and complex expectations; means holds their real parts
177 : !> \param status 0 success, 1 unconverged, -1 input, -2 metric, -4 numerical failure
178 : !> \param unitary Optional mask for compressions of continuum unitary operators, whose full normal Gram is S
179 : ! **************************************************************************************************
180 1972 : SUBROUTINE quadratic_dense_complex(operators, metric, shift, weights, nstates, iterative, max_subspace, &
181 986 : max_iter, tolerance, metric_tolerance, RESULT, status, unitary)
182 : COMPLEX(dp), INTENT(IN), TARGET :: operators(:, :, :), metric(:, :)
183 : COMPLEX(dp), INTENT(IN) :: shift(:)
184 : REAL(dp), INTENT(IN) :: weights(:), tolerance, metric_tolerance
185 : INTEGER, INTENT(IN) :: nstates, max_subspace, max_iter
186 : LOGICAL, INTENT(IN) :: iterative
187 : TYPE(quadratic_result_type), INTENT(OUT) :: result
188 : INTEGER, INTENT(OUT) :: status
189 : LOGICAL, INTENT(IN), OPTIONAL :: unitary(:)
190 986 : TYPE(dense_operator_type) :: operator_backend
191 986 : REAL(dp), ALLOCATABLE :: eigenvalues(:)
192 : INTEGER :: n
193 : #if defined(__HAS_IEEE_EXCEPTIONS)
194 : LOGICAL :: halt(SIZE(IEEE_ALL))
195 : #endif
196 :
197 986 : status = -1
198 986 : n = SIZE(metric, 1)
199 986 : IF (n < 1 .OR. SIZE(metric, 2) /= n) RETURN
200 986 : IF (SIZE(shift) < 1 .OR. SIZE(weights) /= SIZE(shift)) RETURN
201 3928 : IF (ANY(SHAPE(operators) /= [n, n, SIZE(shift)])) RETURN
202 982 : IF (.NOT. ieee_is_finite(metric_tolerance)) RETURN
203 982 : IF (metric_tolerance <= 0.0_dp .OR. metric_tolerance >= 1.0_dp) RETURN
204 982 : IF (.NOT. hermitian(metric)) RETURN
205 465812 : IF (.NOT. ALL(ieee_is_finite(REAL(operators, dp))) .OR. &
206 : .NOT. ALL(ieee_is_finite(AIMAG(operators)))) RETURN
207 982 : operator_backend%operators => operators
208 982 : operator_backend%metric => metric
209 5892 : ALLOCATE (operator_backend%factor(n, n), eigenvalues(n))
210 86254 : operator_backend%factor(:, :) = metric
211 982 : CALL diagonalize(operator_backend%factor, eigenvalues, status)
212 982 : IF (status /= 0) RETURN
213 982 : status = -2
214 982 : IF (eigenvalues(n) <= 0.0_dp) RETURN
215 982 : IF (eigenvalues(1) <= metric_tolerance*eigenvalues(n)) RETURN
216 79690 : operator_backend%factor(:, :) = metric
217 : #if defined(__HAS_IEEE_EXCEPTIONS)
218 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
219 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
220 : #endif
221 978 : CALL zpotrf('L', n, operator_backend%factor, n, status)
222 : #if defined(__HAS_IEEE_EXCEPTIONS)
223 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
224 : #endif
225 978 : IF (status /= 0) THEN
226 0 : status = -2
227 0 : RETURN
228 : END IF
229 : CALL quadratic_solve_core(n, shift, weights, nstates, iterative, max_subspace, max_iter, tolerance, &
230 1943 : operator_backend, RESULT, status, unitary)
231 1968 : END SUBROUTINE quadratic_dense_complex
232 :
233 : ! **************************************************************************************************
234 : !> \brief Apply a dense covariant operator or the metric.
235 : !> \param this Dense operator backend
236 : !> \param OPERATOR Zero for metric, otherwise an observable index
237 : !> \param v Input block
238 : !> \param y Output block
239 : !> \param status Return status
240 : ! **************************************************************************************************
241 21590 : SUBROUTINE dense_apply(this, OPERATOR, v, y, status)
242 : CLASS(dense_operator_type), INTENT(INOUT) :: this
243 : INTEGER, INTENT(IN) :: operator
244 : COMPLEX(dp), INTENT(IN) :: v(:, :)
245 : COMPLEX(dp), INTENT(OUT) :: y(:, :)
246 : INTEGER, INTENT(OUT) :: status
247 :
248 21590 : status = 0
249 21590 : IF (OPERATOR == 0) THEN
250 23045016 : y = MATMUL(this%metric, v)
251 : ELSE
252 15360908 : y = MATMUL(this%operators(:, :, OPERATOR), v)
253 : END IF
254 21590 : END SUBROUTINE dense_apply
255 :
256 : ! **************************************************************************************************
257 : !> \brief Apply the actual dense covariant adjoint, without assuming a normal operator.
258 : !> \param this Dense operator backend
259 : !> \param OPERATOR Observable index
260 : !> \param v Input block
261 : !> \param y Adjoint output block
262 : !> \param status Return status
263 : ! **************************************************************************************************
264 4834 : SUBROUTINE dense_adjoint(this, OPERATOR, v, y, status)
265 : CLASS(dense_operator_type), INTENT(INOUT) :: this
266 : INTEGER, INTENT(IN) :: operator
267 : COMPLEX(dp), INTENT(IN) :: v(:, :)
268 : COMPLEX(dp), INTENT(OUT) :: y(:, :)
269 : INTEGER, INTENT(OUT) :: status
270 :
271 11442079 : y = MATMUL(CONJG(TRANSPOSE(this%operators(:, :, OPERATOR))), v)
272 4834 : status = 0
273 4834 : END SUBROUTINE dense_adjoint
274 :
275 : ! **************************************************************************************************
276 : !> \brief Reuse one dense Cholesky factorization for AO-metric solutions.
277 : !> \param this Dense operator backend
278 : !> \param v Right hand sides
279 : !> \param y Solutions
280 : !> \param status Return status
281 : ! **************************************************************************************************
282 8978 : SUBROUTINE dense_solve(this, v, y, status)
283 : CLASS(dense_operator_type), INTENT(INOUT) :: this
284 : COMPLEX(dp), INTENT(IN) :: v(:, :)
285 : COMPLEX(dp), INTENT(OUT) :: y(:, :)
286 : INTEGER, INTENT(OUT) :: status
287 :
288 : INTEGER :: n
289 : #if defined(__HAS_IEEE_EXCEPTIONS)
290 : LOGICAL :: halt(SIZE(IEEE_ALL))
291 : #endif
292 :
293 8978 : n = SIZE(this%factor, 1)
294 590367 : y = v
295 : #if defined(__HAS_IEEE_EXCEPTIONS)
296 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
297 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
298 : #endif
299 8978 : CALL zpotrs('L', n, SIZE(v, 2), this%factor, n, y, n, status)
300 : #if defined(__HAS_IEEE_EXCEPTIONS)
301 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
302 : #endif
303 8978 : IF (status /= 0) status = -4
304 8978 : END SUBROUTINE dense_solve
305 :
306 : ! **************************************************************************************************
307 : !> \brief Solve sum_j w_j^2 A_j S^-1 A_j c = lambda S c without squaring sparse operators.
308 : !> \param n AO/spinor dimension
309 : !> \param shift Query values for the observables
310 : !> \param weights Nonnegative observable scales
311 : !> \param nstates Number of lowest eigenpairs
312 : !> \param iterative Whether to use a bounded restarted block Ritz space
313 : !> \param max_subspace Maximum trial-space size
314 : !> \param max_iter Maximum Ritz iterations
315 : !> \param tolerance Relative S^-1 eigenpair residual threshold
316 : !> \param apply Apply covariant observables or metric (operator zero)
317 : !> \param solve Apply inverse metric using already validated factors
318 : !> \param RESULT Normalized eigenstates and separated physical residuals
319 : !> \param status 0 success, 1 unconverged, -1 invalid input, -4 numerical failure
320 : ! **************************************************************************************************
321 0 : SUBROUTINE quadratic_solve_real(n, shift, weights, nstates, iterative, max_subspace, max_iter, tolerance, &
322 : apply, solve, RESULT, status)
323 : INTEGER, INTENT(IN) :: n
324 : REAL(dp), INTENT(IN) :: shift(:), weights(:)
325 : INTEGER, INTENT(IN) :: nstates
326 : LOGICAL, INTENT(IN) :: iterative
327 : INTEGER, INTENT(IN) :: max_subspace, max_iter
328 : REAL(dp), INTENT(IN) :: tolerance
329 :
330 : PROCEDURE(operator_apply) :: apply
331 : PROCEDURE(metric_solve) :: solve
332 : TYPE(quadratic_result_type), INTENT(OUT) :: result
333 : INTEGER, INTENT(OUT) :: status
334 :
335 : CALL quadratic_solve_complex(n, CMPLX(shift, 0.0_dp, dp), weights, nstates, iterative, max_subspace, &
336 0 : max_iter, tolerance, apply, solve, RESULT, status, apply)
337 0 : END SUBROUTINE quadratic_solve_real
338 :
339 : ! **************************************************************************************************
340 : !> \brief Solve sum_j w_j^2 (A_j-z_j S)^dagger S^-1 (A_j-z_j S) c = lambda S c.
341 : !> \param n AO/spinor dimension
342 : !> \param shift Complex queries for arbitrary, possibly nonnormal operators
343 : !> \param weights Nonnegative residual scales
344 : !> \param nstates Number of lowest eigenpairs
345 : !> \param iterative Whether to use a bounded restarted block Ritz space
346 : !> \param max_subspace Maximum trial-space size
347 : !> \param max_iter Maximum Ritz iterations
348 : !> \param tolerance Relative S^-1 eigenpair residual threshold
349 : !> \param apply Apply covariant operators or the positive metric (operator zero)
350 : !> \param solve Solve with the already validated positive metric
351 : !> \param RESULT Eigenpairs, residuals and complex expectations; means holds their real parts
352 : !> \param status 0 success, 1 unconverged, -1 input, -4 numerical failure
353 : !> \param adjoint Apply A_j^dagger, not the coefficient-space S-adjoint; positive operator indices
354 : !> \param unitary True for compressions of continuum unitary operators; retain full norms and projection leakage
355 : ! **************************************************************************************************
356 0 : SUBROUTINE quadratic_solve_complex(n, shift, weights, nstates, iterative, max_subspace, max_iter, tolerance, &
357 0 : apply, solve, RESULT, status, adjoint, unitary)
358 : INTEGER, INTENT(IN) :: n
359 : COMPLEX(dp), INTENT(IN) :: shift(:)
360 : REAL(dp), INTENT(IN) :: weights(:)
361 : INTEGER, INTENT(IN) :: nstates
362 : LOGICAL, INTENT(IN) :: iterative
363 : INTEGER, INTENT(IN) :: max_subspace, max_iter
364 : REAL(dp), INTENT(IN) :: tolerance
365 :
366 : PROCEDURE(operator_apply) :: apply, adjoint
367 : PROCEDURE(metric_solve) :: solve
368 : TYPE(quadratic_result_type), INTENT(OUT) :: result
369 : INTEGER, INTENT(OUT) :: status
370 : LOGICAL, INTENT(IN), OPTIONAL :: unitary(:)
371 :
372 : TYPE(callback_operator_type) :: operator_backend
373 :
374 0 : operator_backend%apply_callback => apply
375 0 : operator_backend%adjoint_callback => adjoint
376 0 : operator_backend%solve_callback => solve
377 : CALL quadratic_solve_core(n, shift, weights, nstates, iterative, max_subspace, max_iter, tolerance, &
378 0 : operator_backend, RESULT, status, unitary)
379 0 : END SUBROUTINE quadratic_solve_complex
380 :
381 : ! **************************************************************************************************
382 : !> \brief Shared quadratic solver using an explicit operator backend.
383 : !> \param n AO/spinor dimension
384 : !> \param shift Complex query values
385 : !> \param weights Nonnegative observable scales
386 : !> \param nstates Number of lowest eigenpairs
387 : !> \param iterative Whether to use a bounded restarted block Ritz space
388 : !> \param max_subspace Maximum trial-space size
389 : !> \param max_iter Maximum Ritz iterations
390 : !> \param tolerance Relative S^-1 eigenpair residual threshold
391 : !> \param operator_backend Operator application and metric-solve backend
392 : !> \param RESULT Eigenpairs, residuals and expectations
393 : !> \param status 0 success, 1 unconverged, -1 input, -4 numerical failure
394 : !> \param unitary Optional continuum-unitary mask
395 : ! **************************************************************************************************
396 1956 : SUBROUTINE quadratic_solve_core(n, shift, weights, nstates, iterative, max_subspace, max_iter, tolerance, &
397 978 : operator_backend, RESULT, status, unitary)
398 : INTEGER, INTENT(IN) :: n, nstates, max_subspace, max_iter
399 : COMPLEX(dp), INTENT(IN) :: shift(:)
400 : REAL(dp), INTENT(IN) :: weights(:), tolerance
401 : LOGICAL, INTENT(IN) :: iterative
402 : CLASS(quadratic_operator_type), INTENT(INOUT) :: operator_backend
403 : TYPE(quadratic_result_type), INTENT(OUT) :: result
404 : INTEGER, INTENT(OUT) :: status
405 : LOGICAL, INTENT(IN), OPTIONAL :: unitary(:)
406 978 : COMPLEX(dp), ALLOCATABLE :: v(:, :), sv(:, :), qv(:, :), projected(:, :), sc(:, :), qc(:, :), &
407 978 : r(:, :), z(:, :), trial(:, :), strial(:, :), work(:, :), invwork(:, :), &
408 978 : restarted(:, :)
409 978 : REAL(dp), ALLOCATABLE :: eigenvalues(:)
410 : REAL(dp) :: square, scale, norm_before, norm_after, full_square, lost_square
411 : COMPLEX(dp) :: coefficient
412 : INTEGER :: capacity, m, old, i, j, p, iteration, seed, initial, pass, keep
413 : LOGICAL :: accepted
414 978 : LOGICAL, ALLOCATABLE :: continuum(:)
415 :
416 978 : status = -1
417 978 : IF (n < 1 .OR. nstates < 1 .OR. nstates > n .OR. max_iter < 1) RETURN
418 978 : IF (SIZE(shift) < 1 .OR. SIZE(weights) /= SIZE(shift)) RETURN
419 2934 : ALLOCATE (continuum(SIZE(shift)))
420 978 : continuum = .FALSE.
421 978 : IF (PRESENT(unitary)) THEN
422 13 : IF (SIZE(unitary) /= SIZE(shift)) RETURN
423 41 : continuum(:) = unitary
424 : END IF
425 14556 : IF (.NOT. ALL(ieee_is_finite([REAL(shift, dp), AIMAG(shift), weights, tolerance]))) RETURN
426 6156 : IF (ANY(weights < 0.0_dp) .OR. MAXVAL(weights) <= 0.0_dp .OR. tolerance <= 0.0_dp) RETURN
427 978 : capacity = n
428 978 : IF (iterative) capacity = MIN(n, max_subspace)
429 978 : IF (capacity < nstates .OR. (capacity < n .AND. capacity < 2*nstates)) RETURN
430 : ALLOCATE (v(n, capacity), sv(n, capacity), qv(n, capacity), sc(n, nstates), qc(n, nstates), &
431 23472 : r(n, nstates), z(n, nstates), trial(n, 1), strial(n, 1), work(n, nstates), invwork(n, nstates))
432 : ALLOCATE (result%states(n, nstates), result%eigenvalues(nstates), result%residuals(nstates), &
433 : result%errors(SIZE(shift), nstates), result%means(SIZE(shift), nstates), &
434 : result%complex_means(SIZE(shift), nstates), &
435 17604 : result%projected_errors(SIZE(shift), nstates), result%leakage(SIZE(shift), nstates))
436 16357 : result%leakage = 0.0_dp
437 978 : m = 0
438 978 : seed = 0
439 978 : initial = n
440 978 : IF (iterative) initial = MIN(capacity, MAX(nstates + 2, 2*nstates))
441 6344 : DO j = 1, initial
442 5366 : trial = CMPLX(0.0_dp, 0.0_dp, dp)
443 5366 : IF (.NOT. iterative) THEN
444 5274 : trial(j, 1) = 1.0_dp
445 : ELSE
446 92 : seed = seed + 1
447 92 : CALL seed_vector()
448 : END IF
449 5366 : CALL append_vector(accepted)
450 5366 : IF (status /= 0) RETURN
451 11710 : IF (.NOT. accepted) THEN
452 0 : status = -4
453 0 : RETURN
454 : END IF
455 : END DO
456 978 : CALL apply_q(v(:, :m), qv(:, :m), status)
457 978 : IF (status /= 0) RETURN
458 1210 : DO iteration = 1, max_iter
459 1210 : result%iterations = iteration
460 7260 : ALLOCATE (projected(m, m), eigenvalues(m))
461 3976490 : projected(:, :) = MATMUL(CONJG(TRANSPOSE(v(:, :m))), qv(:, :m))
462 103802 : scale = MAX(1.0_dp, MAXVAL(ABS(projected)))
463 103802 : IF (MAXVAL(ABS(projected - CONJG(TRANSPOSE(projected)))) > 1.e-8_dp*scale) THEN
464 0 : status = -4
465 0 : RETURN
466 : END IF
467 206394 : projected(:, :) = 0.5_dp*(projected + CONJG(TRANSPOSE(projected)))
468 1210 : CALL diagonalize(projected, eigenvalues, status)
469 1210 : IF (status /= 0) RETURN
470 6725 : result%eigenvalues(:) = eigenvalues(:nstates)
471 4338822 : result%states(:, :) = MATMUL(v(:, :m), projected(:, :nstates))
472 4265019 : sc(:, :) = MATMUL(sv(:, :m), projected(:, :nstates))
473 2165643 : qc(:, :) = MATMUL(qv(:, :m), projected(:, :nstates))
474 3162 : DO pass = 1, 2
475 2186 : IF (pass == 2) THEN
476 : ! Cached metric/operator images can drift through repeated Ritz restarts.
477 976 : CALL operator_backend%apply(0, result%states, sc, status)
478 976 : IF (status /= 0) RETURN
479 5943 : DO j = 1, nstates
480 58605 : square = REAL(DOT_PRODUCT(result%states(:, j), sc(:, j)), dp)
481 4967 : IF (.NOT. ieee_is_finite(square) .OR. square <= 0.0_dp) THEN
482 0 : status = -4
483 0 : RETURN
484 : END IF
485 59581 : result%states(:, j) = result%states(:, j)/SQRT(square)
486 : END DO
487 976 : CALL operator_backend%apply(0, result%states, sc, status)
488 976 : IF (status /= 0) RETURN
489 976 : CALL apply_q(result%states, qc, status)
490 976 : IF (status /= 0) RETURN
491 : END IF
492 12668 : DO j = 1, nstates
493 134594 : r(:, j) = qc(:, j) - result%eigenvalues(j)*sc(:, j)
494 : END DO
495 2186 : CALL operator_backend%solve(r, z, status)
496 2186 : IF (status /= 0) RETURN
497 12668 : DO j = 1, nstates
498 132408 : square = REAL(DOT_PRODUCT(r(:, j), z(:, j)), dp)
499 10482 : IF (.NOT. ieee_is_finite(square) .OR. square < -EPSILON(square)*scale) THEN
500 0 : status = -4
501 0 : RETURN
502 : END IF
503 12668 : result%residuals(j) = SQRT(MAX(0.0_dp, square))
504 : END DO
505 15338 : IF (ANY(result%residuals > tolerance*MAX(1.0_dp, ABS(result%eigenvalues)))) EXIT
506 : END DO
507 6233 : IF (ALL(result%residuals <= tolerance*MAX(1.0_dp, ABS(result%eigenvalues)))) THEN
508 976 : DEALLOCATE (projected, eigenvalues)
509 976 : EXIT
510 : END IF
511 234 : IF (.NOT. iterative .OR. m == n .OR. iteration == max_iter) THEN
512 2 : status = 1
513 2 : RETURN
514 : END IF
515 : ! Fill remaining slots with a partial residual block before restarting.
516 232 : IF (m == capacity) THEN
517 : ! Retain nearby Ritz vectors across restarts, important for clustered spectra.
518 64 : keep = MIN(m, MAX(nstates, capacity/2))
519 256 : ALLOCATE (restarted(n, keep))
520 348464 : restarted(:, :) = MATMUL(v(:, :m), projected(:, :keep))
521 11184 : v(:, :keep) = restarted
522 348464 : restarted(:, :) = MATMUL(sv(:, :m), projected(:, :keep))
523 11184 : sv(:, :keep) = restarted
524 177264 : restarted(:, :) = MATMUL(qv(:, :m), projected(:, :keep))
525 11184 : qv(:, :keep) = restarted
526 64 : DEALLOCATE (restarted)
527 64 : m = keep
528 : END IF
529 232 : DEALLOCATE (projected, eigenvalues)
530 232 : old = m
531 692 : DO j = 1, nstates
532 526 : IF (m == capacity) EXIT
533 12818 : trial(:, 1) = z(:, j)
534 460 : CALL append_vector(accepted)
535 692 : IF (status /= 0) RETURN
536 : END DO
537 232 : IF (m == old) THEN
538 0 : seed = seed + 1
539 0 : CALL seed_vector()
540 0 : CALL append_vector(accepted)
541 0 : IF (status /= 0) RETURN
542 0 : IF (.NOT. accepted) THEN
543 0 : status = 1
544 0 : RETURN
545 : END IF
546 : END IF
547 232 : CALL apply_q(v(:, old + 1:m), qv(:, old + 1:m), status)
548 1442 : IF (status /= 0) RETURN
549 : END DO
550 10906 : IF (ANY(result%eigenvalues < -tolerance*MAX(1.0_dp, MAXVAL(ABS(result%eigenvalues))))) THEN
551 2 : status = -4
552 2 : RETURN
553 : END IF
554 5937 : result%eigenvalues = MAX(0.0_dp, result%eigenvalues)
555 3062 : DO p = 1, SIZE(shift)
556 2088 : CALL operator_backend%apply(p, result%states, work, status)
557 2088 : IF (status /= 0) RETURN
558 12462 : DO j = 1, nstates
559 123022 : result%complex_means(p, j) = DOT_PRODUCT(result%states(:, j), work(:, j))
560 12462 : result%means(p, j) = REAL(result%complex_means(p, j), dp)
561 : END DO
562 125110 : work(:, :) = work - shift(p)*sc
563 2088 : CALL operator_backend%solve(work, invwork, status)
564 2088 : IF (status /= 0) RETURN
565 17612 : DO j = 1, nstates
566 123022 : square = REAL(DOT_PRODUCT(work(:, j), invwork(:, j)), dp)
567 10374 : IF (.NOT. ieee_is_finite(square) .OR. square < -tolerance) THEN
568 0 : status = -4
569 0 : RETURN
570 : END IF
571 10374 : result%projected_errors(p, j) = SQRT(MAX(0.0_dp, square))
572 10374 : IF (continuum(p)) THEN
573 : full_square = (1.0_dp + ABS(shift(p))**2)*REAL(DOT_PRODUCT(result%states(:, j), sc(:, j)), dp) &
574 375 : - 2.0_dp*REAL(CONJG(shift(p))*result%complex_means(p, j), dp)
575 29 : lost_square = full_square - square
576 29 : scale = MAX(1.0_dp, ABS(full_square), ABS(square))
577 29 : IF (.NOT. ieee_is_finite(full_square) .OR. full_square < -tolerance*scale .OR. &
578 : lost_square < -tolerance*scale) THEN
579 0 : status = -4
580 0 : RETURN
581 : END IF
582 29 : result%leakage(p, j) = SQRT(MAX(0.0_dp, lost_square))
583 29 : square = full_square
584 : END IF
585 12462 : result%errors(p, j) = SQRT(MAX(0.0_dp, square))
586 : END DO
587 : END DO
588 5937 : DO j = 1, nstates
589 15337 : square = SUM((weights*result%errors(:, j))**2)
590 5937 : IF (ABS(square - result%eigenvalues(j)) > 10.0_dp*tolerance*MAX(1.0_dp, square)) THEN
591 0 : status = -4
592 0 : RETURN
593 : END IF
594 : END DO
595 978 : status = 0
596 : CONTAINS
597 : ! **************************************************************************************************
598 : !> \brief Deterministic complex trial block, independent of compiler random-number generators.
599 : ! **************************************************************************************************
600 92 : SUBROUTINE seed_vector()
601 : INTEGER :: k
602 :
603 1952 : DO k = 1, n
604 : trial(k, 1) = CMPLX(SIN(REAL(k, dp)*REAL(seed + 1, dp)*1.4142135623730951_dp), &
605 1952 : COS(REAL(k, dp)*REAL(seed + 2, dp)*1.7320508075688772_dp), dp)
606 : END DO
607 92 : END SUBROUTINE seed_vector
608 : ! **************************************************************************************************
609 : !> \brief Twice-reorthogonalized S-metric expansion, rejecting linearly dependent vectors.
610 : !> \param accepted Whether a new basis vector was appended
611 : ! **************************************************************************************************
612 5826 : SUBROUTINE append_vector(accepted)
613 : LOGICAL, INTENT(OUT) :: accepted
614 :
615 5826 : accepted = .FALSE.
616 5826 : CALL operator_backend%apply(0, trial, strial, status)
617 5826 : IF (status /= 0) RETURN
618 81850 : norm_before = REAL(DOT_PRODUCT(trial(:, 1), strial(:, 1)), dp)
619 5826 : IF (.NOT. ieee_is_finite(norm_before) .OR. norm_before < 0.0_dp) THEN
620 0 : status = -4
621 0 : RETURN
622 : END IF
623 17478 : DO pass = 1, 2
624 83418 : DO i = 1, m
625 1452532 : coefficient = DOT_PRODUCT(v(:, i), strial(:, 1))
626 1452532 : trial(:, 1) = trial(:, 1) - coefficient*v(:, i)
627 1464184 : strial(:, 1) = strial(:, 1) - coefficient*sv(:, i)
628 : END DO
629 : END DO
630 81850 : norm_after = REAL(DOT_PRODUCT(trial(:, 1), strial(:, 1)), dp)
631 5826 : IF (.NOT. ieee_is_finite(norm_after)) THEN
632 0 : status = -4
633 0 : RETURN
634 : END IF
635 5826 : IF (norm_after <= 1.e-24_dp*MAX(TINY(1.0_dp), norm_before)) RETURN
636 5826 : m = m + 1
637 81850 : v(:, m) = trial(:, 1)/SQRT(norm_after)
638 81850 : sv(:, m) = strial(:, 1)/SQRT(norm_after)
639 5826 : accepted = .TRUE.
640 : END SUBROUTINE append_vector
641 : ! **************************************************************************************************
642 : !> \brief Matrix-free covariant normal products, reusing the positive metric factorization.
643 : !> \param block Input AO coefficients
644 : !> \param image Covariant Q times the input
645 : !> \param status Return status
646 : ! **************************************************************************************************
647 2186 : SUBROUTINE apply_q(block, image, status)
648 : COMPLEX(dp), INTENT(IN) :: block(:, :)
649 : COMPLEX(dp), INTENT(OUT) :: image(:, :)
650 : INTEGER, INTENT(OUT) :: status
651 :
652 : COMPLEX(dp), ALLOCATABLE :: ab(:, :), sb(:, :), second(:, :), &
653 : solved(:, :), ss(:, :)
654 : INTEGER :: nrhs, op
655 :
656 2186 : nrhs = SIZE(block, 2)
657 26232 : ALLOCATE (sb(n, nrhs), ab(n, nrhs), solved(n, nrhs), second(n, nrhs), ss(n, nrhs))
658 2186 : CALL operator_backend%apply(0, block, sb, status)
659 2186 : IF (status /= 0) RETURN
660 142641 : image = CMPLX(0.0_dp, 0.0_dp, dp)
661 7020 : DO op = 1, SIZE(shift)
662 4834 : IF (weights(op) == 0.0_dp) CYCLE
663 4834 : CALL operator_backend%apply(op, block, ab, status)
664 4834 : IF (status /= 0) RETURN
665 4834 : IF (continuum(op)) THEN
666 130 : CALL operator_backend%adjoint(op, block, second, status)
667 130 : IF (status /= 0) RETURN
668 : image = image + weights(op)**2*((1.0_dp + ABS(shift(op))**2)*sb &
669 5585 : - CONJG(shift(op))*ab - shift(op)*second)
670 : CYCLE
671 : END IF
672 330663 : ab(:, :) = ab - shift(op)*sb
673 4704 : CALL operator_backend%solve(ab, solved, status)
674 4704 : IF (status /= 0) RETURN
675 4704 : CALL operator_backend%adjoint(op, solved, second, status)
676 4704 : IF (status /= 0) RETURN
677 4704 : CALL operator_backend%apply(0, solved, ss, status)
678 4704 : IF (status /= 0) RETURN
679 342257 : image = image + weights(op)**2*(second - CONJG(shift(op))*ss)
680 : END DO
681 2186 : END SUBROUTINE apply_q
682 : END SUBROUTINE quadratic_solve_core
683 :
684 : ! **************************************************************************************************
685 : !> \brief Forward an operator application to the public callback interface.
686 : !> \param this ...
687 : !> \param OPERATOR ...
688 : !> \param v ...
689 : !> \param y ...
690 : !> \param status ...
691 : ! **************************************************************************************************
692 0 : SUBROUTINE callback_apply(this, OPERATOR, v, y, status)
693 : CLASS(callback_operator_type), INTENT(INOUT) :: this
694 : INTEGER, INTENT(IN) :: operator
695 : COMPLEX(dp), INTENT(IN) :: v(:, :)
696 : COMPLEX(dp), INTENT(OUT) :: y(:, :)
697 : INTEGER, INTENT(OUT) :: status
698 :
699 0 : CALL this%apply_callback(OPERATOR, v, y, status)
700 0 : END SUBROUTINE callback_apply
701 :
702 : ! **************************************************************************************************
703 : !> \brief Forward an adjoint application to the public callback interface.
704 : !> \param this ...
705 : !> \param OPERATOR ...
706 : !> \param v ...
707 : !> \param y ...
708 : !> \param status ...
709 : ! **************************************************************************************************
710 0 : SUBROUTINE callback_adjoint(this, OPERATOR, v, y, status)
711 : CLASS(callback_operator_type), INTENT(INOUT) :: this
712 : INTEGER, INTENT(IN) :: operator
713 : COMPLEX(dp), INTENT(IN) :: v(:, :)
714 : COMPLEX(dp), INTENT(OUT) :: y(:, :)
715 : INTEGER, INTENT(OUT) :: status
716 :
717 0 : CALL this%adjoint_callback(OPERATOR, v, y, status)
718 0 : END SUBROUTINE callback_adjoint
719 :
720 : ! **************************************************************************************************
721 : !> \brief Forward a metric solve to the public callback interface.
722 : !> \param this ...
723 : !> \param v ...
724 : !> \param y ...
725 : !> \param status ...
726 : ! **************************************************************************************************
727 0 : SUBROUTINE callback_solve(this, v, y, status)
728 : CLASS(callback_operator_type), INTENT(INOUT) :: this
729 : COMPLEX(dp), INTENT(IN) :: v(:, :)
730 : COMPLEX(dp), INTENT(OUT) :: y(:, :)
731 : INTEGER, INTENT(OUT) :: status
732 :
733 0 : CALL this%solve_callback(v, y, status)
734 0 : END SUBROUTINE callback_solve
735 :
736 : ! **************************************************************************************************
737 : !> \brief Hermiticity and finiteness check for the dense reference.
738 : !> \param matrix Input square matrix
739 : !> \return Whether the matrix is finite Hermitian
740 : ! **************************************************************************************************
741 1244 : FUNCTION hermitian(matrix) RESULT(valid)
742 : COMPLEX(dp), INTENT(IN) :: matrix(:, :)
743 : LOGICAL :: valid
744 :
745 1244 : valid = .FALSE.
746 399536 : IF (.NOT. ALL(ieee_is_finite(REAL(matrix, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(matrix)))) RETURN
747 399536 : valid = MAXVAL(ABS(matrix - CONJG(TRANSPOSE(matrix)))) <= 1.e-10_dp*MAX(TINY(1.0_dp), MAXVAL(ABS(matrix)))
748 1244 : END FUNCTION hermitian
749 :
750 : ! **************************************************************************************************
751 : !> \brief Small dense Hermitian eigensolve shared by the oracle and Ritz projection.
752 : !> \param matrix Overwritten with eigenvectors
753 : !> \param values Ascending eigenvalues
754 : !> \param status Return status
755 : ! **************************************************************************************************
756 2192 : SUBROUTINE diagonalize(matrix, values, status)
757 : COMPLEX(dp), INTENT(INOUT) :: matrix(:, :)
758 : REAL(dp), INTENT(OUT) :: values(:)
759 : INTEGER, INTENT(OUT) :: status
760 2192 : COMPLEX(dp), ALLOCATABLE :: workspace(:)
761 2192 : REAL(dp), ALLOCATABLE :: rwork(:)
762 : INTEGER :: n
763 : #if defined(__HAS_IEEE_EXCEPTIONS)
764 : LOGICAL :: halt(SIZE(IEEE_ALL))
765 : #endif
766 2192 : n = SIZE(values)
767 10960 : ALLOCATE (workspace(MAX(1, 2*n)), rwork(MAX(1, 3*n - 2)))
768 : #if defined(__HAS_IEEE_EXCEPTIONS)
769 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
770 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
771 : #endif
772 2192 : CALL zheev('V', 'L', n, matrix, n, values, workspace, SIZE(workspace), rwork, status)
773 : #if defined(__HAS_IEEE_EXCEPTIONS)
774 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
775 : #endif
776 15916 : IF (status /= 0 .OR. .NOT. ALL(ieee_is_finite(values))) status = -4
777 2192 : END SUBROUTINE diagonalize
778 0 : END MODULE quadratic_pseudospectrum
|
