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 Dense reference spectral localizers in a nonorthogonal finite basis.
10 : !> \author Thomas D. Kuehne, tkuehne@cp2k.org
11 : ! **************************************************************************************************
12 : MODULE spectral_localizer
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 :
19 : IMPLICIT NONE
20 : PRIVATE
21 : PUBLIC :: localizer_evaluate, localizer_matrix_evaluate, localizer_inertia, localizer_pfaffian_sign
22 : PUBLIC :: localizer_periodic_evaluate, localizer_aii_skew
23 : PUBLIC :: localizer_3d_evaluate, localizer_3d_matrix_evaluate
24 :
25 : INTERFACE localizer_inertia
26 : MODULE PROCEDURE localizer_inertia_complex, localizer_inertia_real
27 : END INTERFACE localizer_inertia
28 :
29 : CONTAINS
30 :
31 : ! **************************************************************************************************
32 : !> \brief Finite three-dimensional class AII localizer, using all three Cartesian moments.
33 : !> \param h Covariant physical-spin Hamiltonian, all up AOs followed by all down AOs
34 : !> \param s Positive physical-spin AO metric
35 : !> \param r Three covariant Hermitian position operators
36 : !> \param point Cartesian query coordinates
37 : !> \param energy Absolute energy query
38 : !> \param kappa Positive energy/length scale
39 : !> \param matrix_tol Relative matrix-structure tolerance
40 : !> \param metric_tol Minimum overlap eigenvalue ratio
41 : !> \param gap_tol Absolute unresolved-gap threshold
42 : !> \param invariant Strong local Z2 index, valid only for status zero
43 : !> \param gap Smallest metric-normalized singular value of the chiral block
44 : !> \param metric_ratio Smallest/largest physical metric eigenvalue
45 : !> \param status Same convention as localizer_evaluate
46 : ! **************************************************************************************************
47 90 : SUBROUTINE localizer_3d_evaluate(h, s, r, point, energy, kappa, matrix_tol, metric_tol, gap_tol, &
48 : invariant, gap, metric_ratio, status)
49 : COMPLEX(dp), INTENT(IN) :: h(:, :), s(:, :), r(:, :, :)
50 : REAL(dp), INTENT(IN) :: point(3), energy, kappa, matrix_tol, &
51 : metric_tol, gap_tol
52 : INTEGER, INTENT(OUT) :: invariant
53 : REAL(dp), INTENT(OUT) :: gap, metric_ratio
54 : INTEGER, INTENT(OUT) :: status
55 :
56 : COMPLEX(dp), PARAMETER :: ii = CMPLX(0.0_dp, 1.0_dp, dp)
57 :
58 90 : COMPLEX(dp), ALLOCATABLE :: block(:, :), shifted(:, :, :)
59 : INTEGER :: j, n
60 :
61 90 : invariant = 0
62 90 : gap = 0.0_dp
63 90 : metric_ratio = 0.0_dp
64 90 : status = -1
65 90 : n = SIZE(h, 1)
66 90 : IF (n < 2 .OR. MOD(n, 2) /= 0 .OR. SIZE(h, 2) /= n) RETURN
67 540 : IF (ANY(SHAPE(s) /= [n, n]) .OR. ANY(SHAPE(r) /= [n, n, 3])) RETURN
68 1078 : IF (.NOT. ALL(ieee_is_finite([point, energy, kappa, matrix_tol, metric_tol, gap_tol]))) RETURN
69 88 : IF (kappa <= 0.0_dp .OR. matrix_tol <= 0.0_dp .OR. matrix_tol >= 1.0_dp) RETURN
70 88 : IF (metric_tol <= 0.0_dp .OR. metric_tol >= 1.0_dp .OR. gap_tol < 0.0_dp) RETURN
71 86 : IF (.NOT. hermitian(h, matrix_tol) .OR. .NOT. hermitian(s, matrix_tol)) RETURN
72 336 : DO j = 1, 3
73 336 : IF (.NOT. hermitian(r(:, :, j), matrix_tol)) RETURN
74 : END DO
75 84 : status = -3
76 84 : IF (.NOT. self_dual(h, matrix_tol) .OR. .NOT. self_dual(s, matrix_tol)) RETURN
77 328 : DO j = 1, 3
78 328 : IF (.NOT. self_dual(r(:, :, j), matrix_tol)) RETURN
79 : END DO
80 656 : ALLOCATE (block(2*n, 2*n), shifted(n, n, 3))
81 328 : DO j = 1, 3
82 2896240 : shifted(:, :, j) = kappa*(r(:, :, j) - point(j)*s)
83 : END DO
84 : ! A = I_aux (H-E S) - i sum_j sigma_j (kappa (X_j-x_j S)).
85 965386 : block(:n, :n) = h - energy*s - ii*shifted(:, :, 3)
86 965386 : block(n + 1:, n + 1:) = h - energy*s + ii*shifted(:, :, 3)
87 965386 : block(:n, n + 1:) = -ii*shifted(:, :, 1) - shifted(:, :, 2)
88 965386 : block(n + 1:, :n) = -ii*shifted(:, :, 1) + shifted(:, :, 2)
89 82 : CALL localizer_3d_matrix_evaluate(block, s, matrix_tol, metric_tol, gap_tol, invariant, gap, metric_ratio, status)
90 82 : END SUBROUTINE localizer_3d_evaluate
91 :
92 : ! **************************************************************************************************
93 : !> \brief Metric gap and real determinant sign of a three-dimensional AII chiral block.
94 : !> \param block Covariant A=I_aux(H-E S)-i kappa sigma.(X-x S); auxiliary index outermost
95 : !> \param s Positive physical-spin AO metric
96 : !> \param matrix_tol Relative structure tolerance
97 : !> \param metric_tol Minimum overlap eigenvalue ratio
98 : !> \param gap_tol Absolute unresolved-gap threshold
99 : !> \param invariant Zero for positive determinant, one for negative; valid only on success
100 : !> \param gap Minimum absolute generalized eigenvalue of [0,A;A^dagger,0]
101 : !> \param metric_ratio Smallest/largest physical metric eigenvalue
102 : !> \param status Same convention as localizer_evaluate
103 : !> \note A positive atomic Hamiltonian gives positive determinant in every nonsingular AO basis.
104 : !> Unlike the 2D Pfaffian, no basis-orientation correction is required for this determinant.
105 : ! **************************************************************************************************
106 93 : SUBROUTINE localizer_3d_matrix_evaluate(block, s, matrix_tol, metric_tol, gap_tol, &
107 : invariant, gap, metric_ratio, status)
108 : COMPLEX(dp), INTENT(IN) :: block(:, :), s(:, :)
109 : REAL(dp), INTENT(IN) :: matrix_tol, metric_tol, gap_tol
110 : INTEGER, INTENT(OUT) :: invariant, status
111 : REAL(dp), INTENT(OUT) :: gap, metric_ratio
112 93 : COMPLEX(dp), ALLOCATABLE :: pencil(:, :), metric(:, :)
113 : COMPLEX(dp) :: value
114 93 : REAL(dp), ALLOCATABLE :: real_block(:, :)
115 : REAL(dp) :: residual, scale, metric_scale, error_budget
116 : INTEGER :: n, m, q, i, j, pi, pj, si, sj, signature, determinant
117 : #if defined(__HAS_IEEE_EXCEPTIONS)
118 : LOGICAL :: halt(SIZE(IEEE_ALL))
119 : #endif
120 :
121 93 : status = -1
122 93 : invariant = 0
123 93 : gap = 0.0_dp
124 93 : metric_ratio = 0.0_dp
125 93 : n = SIZE(s, 1)
126 93 : IF (n < 2 .OR. MOD(n, 2) /= 0 .OR. SIZE(s, 2) /= n) RETURN
127 93 : m = 2*n
128 279 : IF (ANY(SHAPE(block) /= [m, m])) RETURN
129 372 : IF (.NOT. ALL(ieee_is_finite([matrix_tol, metric_tol, gap_tol]))) RETURN
130 93 : IF (matrix_tol <= 0.0_dp .OR. matrix_tol >= 1.0_dp) RETURN
131 93 : IF (metric_tol <= 0.0_dp .OR. metric_tol >= 1.0_dp .OR. gap_tol < 0.0_dp) RETURN
132 93 : IF (.NOT. hermitian(s, matrix_tol)) RETURN
133 8469482 : IF (.NOT. ALL(ieee_is_finite(REAL(block, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(block)))) RETURN
134 93 : status = -3
135 93 : IF (.NOT. self_dual(s, matrix_tol)) RETURN
136 372 : ALLOCATE (real_block(m, m))
137 93 : q = n/2
138 93 : residual = 0.0_dp
139 : ! Q=(I-i C)/sqrt(2), C=J_aux tensor J_spin: Q^dagger A Q is REAL, not skew.
140 19253 : DO j = 1, m
141 19160 : pj = (3 - (j - 1)/q)*q + MODULO(j - 1, q) + 1
142 19160 : sj = MERGE(1, -1, j <= q .OR. j > 3*q)
143 4234741 : DO i = 1, m
144 4215488 : pi = (3 - (i - 1)/q)*q + MODULO(i - 1, q) + 1
145 4215488 : si = MERGE(1, -1, i <= q .OR. i > 3*q)
146 : value = 0.5_dp*(block(i, j) + si*sj*block(pi, pj) + &
147 4215488 : CMPLX(0.0_dp, 1.0_dp, dp)*(si*block(pi, j) - sj*block(i, pj)))
148 4215488 : residual = MAX(residual, ABS(AIMAG(value)))
149 4234648 : real_block(i, j) = REAL(value, dp)
150 : END DO
151 : END DO
152 4234741 : scale = MAX(TINY(1.0_dp), MAXVAL(ABS(block)))
153 4234741 : IF (.NOT. ALL(ieee_is_finite(real_block)) .OR. .NOT. ieee_is_finite(residual)) RETURN
154 93 : IF (residual > matrix_tol*scale) RETURN
155 :
156 651 : ALLOCATE (pencil(2*m, 2*m), metric(m, m))
157 93 : pencil = 0.0_dp
158 4234741 : pencil(:m, m + 1:) = block
159 4234741 : pencil(m + 1:, :m) = CONJG(TRANSPOSE(block))
160 93 : metric = 0.0_dp
161 1063545 : metric(:n, :n) = s
162 1063545 : metric(n + 1:, n + 1:) = s
163 : CALL localizer_matrix_evaluate(pencil, metric, .FALSE., matrix_tol, metric_tol, gap_tol, &
164 93 : signature, gap, metric_ratio, status)
165 93 : IF (status /= 0) RETURN
166 85 : status = -4
167 85 : IF (signature /= 0) RETURN
168 : ! ||S||_max <= lambda_max(S), hence this lower bound is conservative.
169 969361 : metric_scale = metric_ratio*MAXVAL(ABS(s))
170 85 : status = 1
171 85 : IF (gap <= REAL(m, dp)*(residual + 100.0_dp*EPSILON(1.0_dp)*scale)/metric_scale) RETURN
172 85 : DEALLOCATE (pencil, metric)
173 3859757 : real_block = real_block/scale
174 : ! Frobenius residuals bound the spectral perturbation after metric normalization.
175 85 : error_budget = gap*metric_scale/scale - REAL(m, dp)*(residual/scale + 100.0_dp*EPSILON(1.0_dp))
176 : #if defined(__HAS_IEEE_EXCEPTIONS)
177 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
178 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
179 : #endif
180 85 : CALL real_determinant_sign(real_block, error_budget, determinant, status)
181 : #if defined(__HAS_IEEE_EXCEPTIONS)
182 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
183 : #endif
184 85 : IF (status /= 0) RETURN
185 85 : invariant = (1 - determinant)/2
186 85 : status = 0
187 178 : END SUBROUTINE localizer_3d_matrix_evaluate
188 :
189 : ! **************************************************************************************************
190 : !> \brief LU determinant sign with a QR fallback when pivot growth exhausts the error budget.
191 : !> \param matrix Scaled real nonsingular chiral block
192 : !> \param error_budget Maximum admissible Frobenius backward error from the independent metric gap
193 : !> \param determinant Sign of the determinant, only valid on success
194 : !> \param status Zero on success, one for unresolved numerical error, minus four for LAPACK failure
195 : !> \note The caller disables IEEE halting while LAPACK is active and restores it afterwards.
196 : ! **************************************************************************************************
197 85 : SUBROUTINE real_determinant_sign(matrix, error_budget, determinant, status)
198 : REAL(dp), INTENT(IN) :: matrix(:, :), error_budget
199 : INTEGER, INTENT(OUT) :: determinant, status
200 :
201 : INTEGER :: i, info, method, n
202 85 : INTEGER, ALLOCATABLE :: pivots(:)
203 : REAL(dp) :: error, roundoff
204 85 : REAL(dp), ALLOCATABLE :: factors(:, :), permuted(:, :), &
205 85 : PRODUCT(:, :), row(:), tau(:), work(:)
206 :
207 85 : n = SIZE(matrix, 1)
208 1275 : ALLOCATE (factors(n, n), PRODUCT(n, n), permuted(n, n), row(n), pivots(n), tau(n), work(64*n))
209 85 : determinant = 0
210 85 : status = -4
211 93 : DO method = 1, 2
212 4063893 : factors(:, :) = matrix
213 4063893 : permuted(:, :) = matrix
214 93 : IF (method == 1) THEN
215 85 : CALL dgetrf(n, n, factors, n, pivots, info)
216 : ELSE
217 8 : CALL dgeqrf(n, n, factors, n, tau, work, SIZE(work), info)
218 : END IF
219 4063893 : IF (info /= 0 .OR. .NOT. ALL(ieee_is_finite(factors))) CYCLE
220 4063893 : PRODUCT(:, :) = factors
221 18645 : DO i = 1, n
222 2031993 : PRODUCT(i + 1:, i) = 0.0_dp
223 : END DO
224 : ! Include rounding of both the triangular product and the residual evaluation.
225 4063893 : roundoff = 8.0_dp*EPSILON(1.0_dp)*REAL(n, dp)*SQRT(SUM(product**2))
226 93 : IF (method == 1) THEN
227 85 : CALL dtrmm('L', 'L', 'N', 'U', n, n, 1.0_dp, factors, n, product, n)
228 17517 : DO i = 1, n
229 3859672 : row(:) = permuted(i, :)
230 3859672 : permuted(i, :) = permuted(pivots(i), :)
231 3859757 : permuted(pivots(i), :) = row
232 : END DO
233 : ELSE
234 8 : CALL dormqr('L', 'N', n, n, n, factors, n, tau, product, n, work, SIZE(work), info)
235 8 : IF (info /= 0) CYCLE
236 : END IF
237 4063893 : IF (.NOT. ALL(ieee_is_finite(product))) CYCLE
238 4063893 : error = SQRT(SUM((product - permuted)**2))
239 93 : IF (method == 1) THEN
240 : ! Bound cancellation in the reconstructed product by |L| |U|, not by |L U|.
241 3859757 : permuted(:, :) = ABS(factors)
242 3859757 : PRODUCT(:, :) = permuted
243 17517 : DO i = 1, n
244 1929921 : PRODUCT(i + 1:, i) = 0.0_dp
245 : END DO
246 85 : CALL dtrmm('L', 'L', 'N', 'U', n, n, 1.0_dp, permuted, n, product, n)
247 3859757 : roundoff = 8.0_dp*EPSILON(1.0_dp)*REAL(n, dp)*SQRT(SUM(product**2))
248 : END IF
249 93 : error = error + roundoff
250 93 : status = 1
251 93 : IF (.NOT. ieee_is_finite(error) .OR. error >= error_budget) CYCLE
252 85 : determinant = 1
253 17517 : DO i = 1, n
254 17432 : IF (factors(i, i) == 0.0_dp) EXIT
255 17432 : IF (factors(i, i) < 0.0_dp) determinant = -determinant
256 17517 : IF (method == 1) THEN
257 16312 : IF (pivots(i) /= i) determinant = -determinant
258 : ELSE
259 1120 : IF (tau(i) /= 0.0_dp) determinant = -determinant
260 : END IF
261 : END DO
262 85 : IF (i <= n) CYCLE
263 85 : status = 0
264 93 : RETURN
265 : END DO
266 0 : determinant = 0
267 85 : END SUBROUTINE real_determinant_sign
268 :
269 : ! **************************************************************************************************
270 : !> \brief Two-dimensional periodic localizer from covariant trigonometric position operators.
271 : !> \param h Periodic Hamiltonian, including all boundary-crossing matrix elements
272 : !> \param s Positive definite overlap
273 : !> \param cosine Covariant cos(G_j.r) matrices, not cos of projected AO positions
274 : !> \param sine Covariant sin(G_j.r) matrices
275 : !> \param origin Two dimensionless phase offsets G_j.r0
276 : !> \param energy Absolute reference energy
277 : !> \param eta Positive energy scale multiplying the torus coordinates
278 : !> \param aii Use physical spinful time reversal
279 : !> \param matrix_tol Relative structure tolerance
280 : !> \param metric_tol Minimum overlap eigenvalue ratio
281 : !> \param gap_tol Absolute unresolved-gap threshold in energy units
282 : !> \param invariant Index in the same orientation convention as localizer_evaluate
283 : !> \param gap Minimum absolute generalized eigenvalue, in energy units
284 : !> \param metric_ratio Smallest/largest overlap eigenvalue ratio
285 : !> \param status Same status convention as localizer_evaluate
286 : !> \note Adapted from Doll, Loring, Schulz-Baldes (2025), Eq. (1): minus eta times
287 : !> their localizer, conjugated by diag(I,-I), to retain our (H,X,Y) orientation.
288 : !> A resolved finite-matrix signature is not a certificate of bulk convergence.
289 : ! **************************************************************************************************
290 60 : SUBROUTINE localizer_periodic_evaluate(h, s, cosine, sine, origin, energy, eta, aii, &
291 : matrix_tol, metric_tol, gap_tol, invariant, gap, metric_ratio, status)
292 : COMPLEX(dp), INTENT(IN) :: h(:, :), s(:, :), cosine(:, :, :), &
293 : sine(:, :, :)
294 : REAL(dp), INTENT(IN) :: origin(2), energy, eta
295 : LOGICAL, INTENT(IN) :: aii
296 : REAL(dp), INTENT(IN) :: matrix_tol, metric_tol, gap_tol
297 : INTEGER, INTENT(OUT) :: invariant
298 : REAL(dp), INTENT(OUT) :: gap, metric_ratio
299 : INTEGER, INTENT(OUT) :: status
300 :
301 60 : COMPLEX(dp), ALLOCATABLE :: l(:, :), mass(:, :), &
302 60 : shifted_sine(:, :, :)
303 : INTEGER :: j, n
304 :
305 60 : status = -1
306 60 : invariant = 0
307 60 : gap = 0.0_dp
308 60 : metric_ratio = 0.0_dp
309 60 : n = SIZE(h, 1)
310 60 : IF (n < 1 .OR. SIZE(h, 2) /= n) RETURN
311 180 : IF (ANY(SHAPE(s) /= [n, n])) RETURN
312 420 : IF (ANY(SHAPE(cosine) /= [n, n, 2]) .OR. ANY(SHAPE(sine) /= [n, n, 2])) RETURN
313 480 : IF (.NOT. ALL(ieee_is_finite([origin, energy, eta, matrix_tol]))) RETURN
314 60 : IF (eta <= 0.0_dp .OR. matrix_tol <= 0.0_dp .OR. matrix_tol >= 1.0_dp) RETURN
315 54 : IF (.NOT. hermitian(h, matrix_tol) .OR. .NOT. hermitian(s, matrix_tol)) RETURN
316 150 : DO j = 1, 2
317 102 : IF (.NOT. hermitian(cosine(:, :, j), matrix_tol)) RETURN
318 144 : IF (.NOT. hermitian(sine(:, :, j), matrix_tol)) RETURN
319 : END DO
320 48 : IF (aii) THEN
321 12 : status = -3
322 12 : IF (MOD(n, 2) /= 0) RETURN
323 12 : IF (.NOT. self_dual(h, matrix_tol) .OR. .NOT. self_dual(s, matrix_tol)) RETURN
324 36 : DO j = 1, 2
325 24 : IF (.NOT. self_dual(cosine(:, :, j), matrix_tol)) RETURN
326 36 : IF (.NOT. self_dual(sine(:, :, j), matrix_tol)) RETURN
327 : END DO
328 : END IF
329 528 : ALLOCATE (mass(n, n), shifted_sine(n, n, 2), l(2*n, 2*n))
330 1383984 : mass(:, :) = h - energy*s
331 144 : DO j = 1, 2
332 2767968 : mass(:, :) = mass - eta*(s - COS(origin(j))*cosine(:, :, j) - SIN(origin(j))*sine(:, :, j))
333 2768016 : shifted_sine(:, :, j) = eta*(COS(origin(j))*sine(:, :, j) - SIN(origin(j))*cosine(:, :, j))
334 : END DO
335 1383984 : l(:n, :n) = mass
336 1383984 : l(n + 1:, n + 1:) = -mass
337 1383984 : l(:n, n + 1:) = shifted_sine(:, :, 1) - CMPLX(0.0_dp, 1.0_dp, dp)*shifted_sine(:, :, 2)
338 2767920 : l(n + 1:, :n) = CONJG(TRANSPOSE(l(:n, n + 1:)))
339 : CALL localizer_matrix_evaluate(l, s, aii, matrix_tol, metric_tol, gap_tol, &
340 48 : invariant, gap, metric_ratio, status)
341 60 : END SUBROUTINE localizer_periodic_evaluate
342 :
343 : ! **************************************************************************************************
344 : !> \brief Evaluate the 2D class A or AII localizer without projecting onto occupied states.
345 : !> \param h Covariant Hamiltonian matrix
346 : !> \param s Positive definite AO overlap
347 : !> \param r Covariant Cartesian position matrices in the two oriented directions
348 : !> \param point Position in the same coordinates and length units as r
349 : !> \param energy Absolute reference energy in the units of h
350 : !> \param kappa Positive energy/length scale
351 : !> \param aii Use spinful time reversal, with all up AOs followed by all down AOs
352 : !> \param matrix_tol Relative Hermiticity and time-reversal tolerance
353 : !> \param metric_tol Minimum allowed smallest/largest overlap eigenvalue ratio
354 : !> \param gap_tol Absolute unresolved-gap threshold in energy units
355 : !> \param invariant Half-signature (A), or 0 trivial / 1 nontrivial (AII); valid only if status=0
356 : !> \param gap Minimum absolute generalized localizer eigenvalue
357 : !> \param metric_ratio Smallest/largest overlap eigenvalue ratio
358 : !> \param status 0 resolved; 1 unresolved gap; -1 invalid data; -2 bad metric; -3 broken TR;
359 : !> -4 factorization failure; positive values above 1 indicate LAPACK failure
360 : !> \note This is a finite-basis diagnostic, not a certificate of basis or volume convergence.
361 : ! **************************************************************************************************
362 58 : SUBROUTINE localizer_evaluate(h, s, r, point, energy, kappa, aii, matrix_tol, metric_tol, &
363 : gap_tol, invariant, gap, metric_ratio, status)
364 : COMPLEX(KIND=dp), INTENT(IN) :: h(:, :), s(:, :), r(:, :, :)
365 : REAL(KIND=dp), INTENT(IN) :: point(2), energy, kappa
366 : LOGICAL, INTENT(IN) :: aii
367 : REAL(KIND=dp), INTENT(IN) :: matrix_tol, metric_tol, gap_tol
368 : INTEGER, INTENT(OUT) :: invariant
369 : REAL(KIND=dp), INTENT(OUT) :: gap, metric_ratio
370 : INTEGER, INTENT(OUT) :: status
371 :
372 : COMPLEX(KIND=dp), ALLOCATABLE :: localizer(:, :)
373 : INTEGER :: i, n
374 :
375 58 : status = -1
376 58 : invariant = 0
377 58 : gap = 0.0_dp
378 58 : metric_ratio = 0.0_dp
379 58 : n = SIZE(h, 1)
380 58 : IF (n < 1 .OR. SIZE(h, 2) /= n) RETURN
381 348 : IF (ANY(SHAPE(s) /= [n, n]) .OR. ANY(SHAPE(r) /= [n, n, 2])) RETURN
382 580 : IF (.NOT. ALL(ieee_is_finite([point, energy, kappa, matrix_tol, metric_tol, gap_tol]))) RETURN
383 58 : IF (kappa <= 0.0_dp .OR. matrix_tol <= 0.0_dp .OR. matrix_tol >= 1.0_dp) RETURN
384 58 : IF (metric_tol <= 0.0_dp .OR. metric_tol >= 1.0_dp .OR. gap_tol < 0.0_dp) RETURN
385 58 : IF (.NOT. hermitian(h, matrix_tol) .OR. .NOT. hermitian(s, matrix_tol)) RETURN
386 162 : DO i = 1, 2
387 162 : IF (.NOT. hermitian(r(:, :, i), matrix_tol)) RETURN
388 : END DO
389 54 : IF (aii) THEN
390 14 : status = -3
391 14 : IF (MOD(n, 2) /= 0) RETURN
392 14 : IF (.NOT. self_dual(h, matrix_tol) .OR. .NOT. self_dual(s, matrix_tol)) RETURN
393 36 : DO i = 1, 2
394 36 : IF (.NOT. self_dual(r(:, :, i), matrix_tol)) RETURN
395 : END DO
396 : END IF
397 :
398 208 : ALLOCATE (localizer(2*n, 2*n))
399 503172 : localizer(:n, :n) = h - energy*s
400 503172 : localizer(n + 1:, n + 1:) = -localizer(:n, :n)
401 : localizer(:n, n + 1:) = kappa*((r(:, :, 1) - point(1)*s) - &
402 503172 : CMPLX(0.0_dp, 1.0_dp, dp)*(r(:, :, 2) - point(2)*s))
403 1006292 : localizer(n + 1:, :n) = CONJG(TRANSPOSE(localizer(:n, n + 1:)))
404 : CALL localizer_matrix_evaluate(localizer, s, aii, matrix_tol, metric_tol, gap_tol, &
405 52 : invariant, gap, metric_ratio, status)
406 52 : END SUBROUTINE localizer_evaluate
407 :
408 : ! **************************************************************************************************
409 : !> \brief Dense reference evaluation of an assembled localizer with metric diag(S,S).
410 : !> \param localizer Covariant Hermitian localizer
411 : !> \param s Positive definite physical AO/spinor overlap
412 : !> \param aii Whether to use the physical spinful time-reversal structure
413 : !> \param matrix_tol Relative structure tolerance
414 : !> \param metric_tol Minimum overlap eigenvalue ratio
415 : !> \param gap_tol Absolute unresolved-gap tolerance
416 : !> \param invariant Resolved class A or AII index
417 : !> \param gap Minimum absolute generalized eigenvalue
418 : !> \param metric_ratio Smallest/largest overlap eigenvalue ratio
419 : !> \param status Same status convention as localizer_evaluate
420 : ! **************************************************************************************************
421 251 : SUBROUTINE localizer_matrix_evaluate(localizer, s, aii, matrix_tol, metric_tol, gap_tol, &
422 : invariant, gap, metric_ratio, status)
423 : COMPLEX(KIND=dp), INTENT(IN) :: localizer(:, :), s(:, :)
424 : LOGICAL, INTENT(IN) :: aii
425 : REAL(KIND=dp), INTENT(IN) :: matrix_tol, metric_tol, gap_tol
426 : INTEGER, INTENT(OUT) :: invariant, status
427 : REAL(KIND=dp), INTENT(OUT) :: gap, metric_ratio
428 :
429 251 : COMPLEX(KIND=dp), ALLOCATABLE :: a(:, :), b(:, :), work(:)
430 251 : INTEGER, ALLOCATABLE :: iwork(:)
431 : INTEGER :: info, n, negative, pf, pf_reference, &
432 : positive, ldl_positive, ldl_negative
433 : #if defined(__HAS_IEEE_EXCEPTIONS)
434 : LOGICAL :: halt(SIZE(IEEE_ALL))
435 : #endif
436 : REAL(KIND=dp) :: metric_min, residual, threshold
437 251 : REAL(KIND=dp), ALLOCATABLE :: eigenvalues(:), rwork(:), skew(:, :)
438 :
439 251 : status = -1
440 251 : invariant = 0
441 251 : gap = 0.0_dp
442 251 : metric_ratio = 0.0_dp
443 251 : n = SIZE(s, 1)
444 251 : IF (n < 1 .OR. SIZE(s, 2) /= n) RETURN
445 753 : IF (ANY(SHAPE(localizer) /= [2*n, 2*n])) RETURN
446 1004 : IF (.NOT. ALL(ieee_is_finite([matrix_tol, metric_tol, gap_tol]))) RETURN
447 251 : IF (matrix_tol <= 0.0_dp .OR. matrix_tol >= 1.0_dp) RETURN
448 251 : IF (metric_tol <= 0.0_dp .OR. metric_tol >= 1.0_dp .OR. gap_tol < 0.0_dp) RETURN
449 251 : IF (.NOT. hermitian(s, matrix_tol) .OR. .NOT. hermitian(localizer, matrix_tol)) RETURN
450 251 : IF (aii) THEN
451 30 : status = -3
452 30 : IF (MOD(n, 2) /= 0) RETURN
453 30 : IF (.NOT. self_dual(s, matrix_tol)) RETURN
454 : END IF
455 :
456 2510 : ALLOCATE (a(n, n), eigenvalues(n), work(MAX(1, 2*n)), rwork(MAX(1, 3*n)), iwork(1))
457 6267843 : a(:, :) = s
458 : ! Follow CP2K's eigensolvers: LAPACK may use exceptional IEEE intermediates internally.
459 : #if defined(__HAS_IEEE_EXCEPTIONS)
460 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
461 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
462 : #endif
463 251 : CALL zheevd('N', 'U', n, a, n, eigenvalues, work, SIZE(work), rwork, SIZE(rwork), iwork, SIZE(iwork), info)
464 : #if defined(__HAS_IEEE_EXCEPTIONS)
465 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
466 : #endif
467 251 : status = ABS(info) + 2
468 32031 : IF (info /= 0 .OR. .NOT. ALL(ieee_is_finite(eigenvalues))) RETURN
469 251 : status = -2
470 32031 : IF (MAXVAL(eigenvalues) <= 0.0_dp) RETURN
471 63811 : metric_ratio = MINVAL(eigenvalues)/MAXVAL(eigenvalues)
472 251 : IF (metric_ratio <= metric_tol) RETURN
473 31175 : metric_min = MINVAL(eigenvalues)
474 243 : DEALLOCATE (a, eigenvalues, work, rwork, iwork)
475 :
476 1458 : ALLOCATE (a(2*n, 2*n), b(2*n, 2*n))
477 1701 : ALLOCATE (eigenvalues(2*n), work(4*n), rwork(6*n), iwork(1))
478 243 : status = -1
479 48954582 : IF (.NOT. ALL(ieee_is_finite(REAL(localizer, dp))) .OR. &
480 : .NOT. ALL(ieee_is_finite(AIMAG(localizer)))) RETURN
481 24477291 : a(:, :) = localizer
482 243 : b(:, :) = CMPLX(0.0_dp, 0.0_dp, dp)
483 6134971 : b(:n, :n) = s
484 6134971 : b(n + 1:, n + 1:) = s
485 : ! Generalized eigenvalues, unlike eigenvalues of raw AO matrices, are basis invariant.
486 : #if defined(__HAS_IEEE_EXCEPTIONS)
487 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
488 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
489 : #endif
490 : CALL zhegvd(1, 'N', 'U', 2*n, a, 2*n, b, 2*n, eigenvalues, work, SIZE(work), rwork, SIZE(rwork), &
491 243 : iwork, SIZE(iwork), info)
492 : #if defined(__HAS_IEEE_EXCEPTIONS)
493 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
494 : #endif
495 243 : status = ABS(info) + 2
496 62107 : IF (info /= 0 .OR. .NOT. ALL(ieee_is_finite(eigenvalues))) RETURN
497 62107 : gap = MINVAL(ABS(eigenvalues))
498 : ! Bound roundoff and triangle inconsistencies after conversion to the orthonormal metric.
499 : residual = REAL(2*n, dp)*MAXVAL(ABS(localizer - CONJG(TRANSPOSE(localizer))))/metric_min + &
500 30674126 : REAL(n, dp)*MAXVAL(ABS(eigenvalues))*MAXVAL(ABS(s - CONJG(TRANSPOSE(s))))/metric_min
501 : threshold = MAX(gap_tol, residual + 100.0_dp*EPSILON(1.0_dp)*REAL(2*n, dp)* &
502 24539155 : MAX(MAXVAL(ABS(eigenvalues)), MAXVAL(ABS(localizer))/metric_min))
503 243 : status = 1
504 243 : IF (gap <= threshold) RETURN
505 57685 : positive = COUNT(eigenvalues > threshold)
506 57685 : negative = COUNT(eigenvalues < -threshold)
507 209 : CALL localizer_inertia(localizer, ldl_positive, ldl_negative, info)
508 209 : status = -4
509 209 : IF (info /= 0) RETURN
510 209 : IF (positive /= ldl_positive .OR. negative /= ldl_negative) RETURN
511 209 : IF (.NOT. aii) THEN
512 181 : invariant = (ldl_positive - ldl_negative)/2
513 : ELSE
514 28 : status = -3
515 28 : IF (positive /= negative) RETURN
516 112 : ALLOCATE (skew(2*n, 2*n))
517 28 : CALL localizer_aii_skew(localizer, skew, matrix_tol, residual, info)
518 28 : IF (info /= 0) RETURN
519 28 : status = 1
520 28 : IF (gap <= threshold + REAL(2*n, dp)*residual/metric_min) RETURN
521 28 : CALL localizer_pfaffian_sign(skew, pf, info)
522 28 : status = -4
523 28 : IF (info /= 0) RETURN
524 : ! Fix the orientation against the positive atomic Hamiltonian in the SAME AO metric.
525 28 : a(:, :) = CMPLX(0.0_dp, 0.0_dp, dp)
526 1237854 : a(:n, :n) = s
527 1237854 : a(n + 1:, n + 1:) = -s
528 28 : CALL localizer_aii_skew(a, skew, matrix_tol, residual, info)
529 28 : IF (info /= 0) RETURN
530 28 : IF (REAL(2*n, dp)*residual >= metric_min) RETURN
531 28 : CALL localizer_pfaffian_sign(skew, pf_reference, info)
532 28 : IF (info /= 0) RETURN
533 112 : invariant = (1 - pf*pf_reference)/2
534 : END IF
535 209 : status = 0
536 251 : END SUBROUTINE localizer_matrix_evaluate
537 :
538 : ! **************************************************************************************************
539 : !> \brief Count inertia by pivoted Hermitian LDL^dagger, without orthogonalizing the AO basis.
540 : !> \param matrix Hermitian matrix; input is not modified
541 : !> \param positive Number of positive eigenvalues
542 : !> \param negative Number of negative eigenvalues
543 : !> \param status 0 resolved, 1 singular, -1 invalid data, -4 factorization failure
544 : !> \note Pivot magnitudes are not eigenvalues or a metric-aware protection gap.
545 : ! **************************************************************************************************
546 1303 : SUBROUTINE localizer_inertia_complex(matrix, positive, negative, status)
547 : COMPLEX(KIND=dp), INTENT(IN) :: matrix(:, :)
548 : INTEGER, INTENT(OUT) :: positive, negative, status
549 :
550 1303 : COMPLEX(KIND=dp), ALLOCATABLE :: a(:, :), work(:)
551 : INTEGER :: i, info, n
552 1303 : INTEGER, ALLOCATABLE :: pivots(:)
553 : REAL(KIND=dp) :: scale
554 : #if defined(__HAS_IEEE_EXCEPTIONS)
555 : LOGICAL :: halt(SIZE(IEEE_ALL))
556 : #endif
557 :
558 1303 : positive = 0
559 1303 : negative = 0
560 1303 : status = -1
561 1303 : n = SIZE(matrix, 1)
562 1303 : IF (n < 1 .OR. SIZE(matrix, 2) /= n) RETURN
563 1303 : IF (.NOT. hermitian(matrix, 1.e-10_dp)) RETURN
564 : ! Never discard a small imaginary term: only exactly real matrices use LDL^T.
565 133123 : IF (ALL(AIMAG(matrix) == 0.0_dp)) THEN
566 112556 : CALL localizer_inertia_real(REAL(matrix, dp), positive, negative, status)
567 568 : RETURN
568 : END IF
569 22933613 : scale = MAXVAL(ABS(matrix))
570 0 : status = 1
571 733 : IF (scale <= TINY(1.0_dp)) RETURN
572 5864 : ALLOCATE (a(n, n), pivots(n), work(MAX(1, 64*n)))
573 22933613 : a(:, :) = matrix/scale
574 : #if defined(__HAS_IEEE_EXCEPTIONS)
575 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
576 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
577 : #endif
578 733 : CALL zhetrf('L', n, a, n, pivots, work, SIZE(work), info)
579 : #if defined(__HAS_IEEE_EXCEPTIONS)
580 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
581 : #endif
582 733 : IF (info > 0) RETURN
583 731 : status = -4
584 731 : IF (info < 0) RETURN
585 45867142 : IF (.NOT. ALL(ieee_is_finite(REAL(a, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(a)))) RETURN
586 : CALL block_inertia([(REAL(a(i, i), dp), i=1, n)], [(ABS(a(i + 1, i)), i=1, n - 1)], &
587 244357 : pivots, positive, negative, status)
588 1303 : END SUBROUTINE localizer_inertia_complex
589 :
590 : ! **************************************************************************************************
591 : !> \brief Native real symmetric LDL^T inertia, without converting to complex storage.
592 : !> \param matrix Real symmetric matrix; input is not modified
593 : !> \param positive Number of positive eigenvalues
594 : !> \param negative Number of negative eigenvalues
595 : !> \param status 0 resolved, 1 singular, -1 invalid data, -4 factorization failure
596 : ! **************************************************************************************************
597 1114 : SUBROUTINE localizer_inertia_real(matrix, positive, negative, status)
598 : REAL(dp), INTENT(IN) :: matrix(:, :)
599 : INTEGER, INTENT(OUT) :: positive, negative, status
600 1114 : REAL(dp), ALLOCATABLE :: a(:, :), work(:)
601 1114 : INTEGER, ALLOCATABLE :: pivots(:)
602 : INTEGER :: i, info, n
603 : REAL(dp) :: scale
604 : #if defined(__HAS_IEEE_EXCEPTIONS)
605 : LOGICAL :: halt(SIZE(IEEE_ALL))
606 : #endif
607 :
608 1114 : positive = 0
609 1114 : negative = 0
610 1114 : status = -1
611 1114 : n = SIZE(matrix, 1)
612 1114 : IF (n < 1 .OR. SIZE(matrix, 2) /= n) RETURN
613 157910 : IF (.NOT. ALL(ieee_is_finite(matrix))) RETURN
614 157910 : scale = MAXVAL(ABS(matrix))
615 157910 : IF (MAXVAL(ABS(matrix - TRANSPOSE(matrix))) > 1.e-10_dp*MAX(TINY(1.0_dp), scale)) RETURN
616 1112 : status = 1
617 1112 : IF (scale <= TINY(1.0_dp)) RETURN
618 8880 : ALLOCATE (a(n, n), pivots(n), work(MAX(1, 64*n)))
619 157882 : a(:, :) = matrix/scale
620 : #if defined(__HAS_IEEE_EXCEPTIONS)
621 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
622 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
623 : #endif
624 1110 : CALL dsytrf('L', n, a, n, pivots, work, SIZE(work), info)
625 : #if defined(__HAS_IEEE_EXCEPTIONS)
626 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
627 : #endif
628 1110 : IF (info > 0) RETURN
629 1110 : status = -4
630 157882 : IF (info < 0 .OR. .NOT. ALL(ieee_is_finite(a))) RETURN
631 : CALL block_inertia([(a(i, i), i=1, n)], [(ABS(a(i + 1, i)), i=1, n - 1)], &
632 36298 : pivots, positive, negative, status)
633 1114 : END SUBROUTINE localizer_inertia_real
634 :
635 : ! **************************************************************************************************
636 : !> \brief Count signs of lower Bunch-Kaufman D blocks, shared by real and complex factors.
637 : !> \param diagonal Real diagonal of the factored matrix
638 : !> \param coupling Magnitudes of first subdiagonal entries
639 : !> \param pivots LAPACK lower-triangle pivot encoding
640 : !> \param positive Positive inertia count
641 : !> \param negative Negative inertia count
642 : !> \param status 0 resolved, 1 singular, -4 invalid block structure
643 : !> \note A Hermitian 2x2 block has the same eigenvalues when its offdiagonal is replaced by its modulus.
644 : ! **************************************************************************************************
645 1841 : SUBROUTINE block_inertia(diagonal, coupling, pivots, positive, negative, status)
646 : REAL(dp), INTENT(IN) :: diagonal(:), coupling(:)
647 : INTEGER, INTENT(IN) :: pivots(:)
648 : INTEGER, INTENT(OUT) :: positive, negative, status
649 :
650 : INTEGER :: i, n
651 : REAL(dp) :: values(2)
652 : INTERFACE
653 : SUBROUTINE dlae2(a, b, c, rt1, rt2)
654 : IMPORT :: dp
655 : REAL(dp), INTENT(IN) :: a, b, c
656 : REAL(dp), INTENT(OUT) :: rt1, rt2
657 : END SUBROUTINE dlae2
658 : END INTERFACE
659 :
660 1841 : positive = 0
661 1841 : negative = 0
662 1841 : status = -4
663 1841 : n = SIZE(pivots)
664 1841 : i = 1
665 52590 : DO WHILE (i <= n)
666 50749 : IF (pivots(i) == 0) RETURN
667 50749 : IF (pivots(i) > 0) THEN
668 92622 : values = [diagonal(i), 0.0_dp]
669 30874 : IF (values(1) == 0.0_dp) THEN
670 0 : status = 1
671 0 : RETURN
672 : END IF
673 30874 : i = i + 1
674 : ELSE
675 19875 : IF (i == n) RETURN
676 19875 : IF (pivots(i + 1) /= pivots(i)) RETURN
677 : ! A 2x2 pivot can have zero diagonal and still contribute one sign of each kind.
678 19875 : CALL dlae2(diagonal(i), coupling(i), diagonal(i + 1), values(1), values(2))
679 59625 : IF (ANY(values == 0.0_dp)) THEN
680 0 : status = 1
681 0 : RETURN
682 : END IF
683 19875 : i = i + 2
684 : END IF
685 152247 : positive = positive + COUNT(values > 0.0_dp)
686 154088 : negative = negative + COUNT(values < 0.0_dp)
687 : END DO
688 1841 : status = 0
689 : END SUBROUTINE block_inertia
690 :
691 : ! **************************************************************************************************
692 : !> \brief Test finiteness and Hermiticity before passing only one triangle to LAPACK.
693 : !> \param a Square matrix
694 : !> \param tolerance Relative tolerance
695 : !> \return Whether the matrix is Hermitian and finite
696 : ! **************************************************************************************************
697 2846 : LOGICAL FUNCTION hermitian(a, tolerance)
698 : COMPLEX(KIND=dp), INTENT(IN) :: a(:, :)
699 : REAL(KIND=dp), INTENT(IN) :: tolerance
700 :
701 2846 : hermitian = .FALSE.
702 142566540 : IF (.NOT. ALL(ieee_is_finite(REAL(a, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(a)))) RETURN
703 142542992 : hermitian = MAXVAL(ABS(a - CONJG(TRANSPOSE(a)))) <= tolerance*MAX(TINY(1.0_dp), MAXVAL(ABS(a)))
704 2844 : END FUNCTION hermitian
705 :
706 : ! **************************************************************************************************
707 : !> \brief Check J conj(A) J^dagger=A for J=(0,I;-I,0), without assuming spin conservation.
708 : !> \param a Even-dimensional matrix in the fixed up/down ordering
709 : !> \param tolerance Relative tolerance
710 : !> \return Whether physical fermionic time reversal is respected
711 : ! **************************************************************************************************
712 657 : LOGICAL FUNCTION self_dual(a, tolerance)
713 : COMPLEX(KIND=dp), INTENT(IN) :: a(:, :)
714 : REAL(KIND=dp), INTENT(IN) :: tolerance
715 :
716 : INTEGER :: n
717 : REAL(KIND=dp) :: scale
718 :
719 657 : n = SIZE(a, 1)/2
720 13202191 : scale = MAX(TINY(1.0_dp), MAXVAL(ABS(a)))
721 : self_dual = MAXVAL(ABS(a(:n, :n) - CONJG(a(n + 1:, n + 1:)))) <= tolerance*scale .AND. &
722 6644430 : MAXVAL(ABS(a(:n, n + 1:) + CONJG(a(n + 1:, :n)))) <= tolerance*scale
723 657 : END FUNCTION self_dual
724 :
725 : ! **************************************************************************************************
726 : !> \brief Convert the class AII localizer to a real skew matrix using Q=(I-i C)/sqrt(2).
727 : !> \param localizer Hermitian localizer; auxiliary index outside physical up/down index
728 : !> \param skew Real part of i Q^dagger L Q
729 : !> \param tolerance Relative imaginary/skew residual tolerance
730 : !> \param residual Entrywise bound for the discarded real and symmetric components
731 : !> \param status Zero on success; negative if the AII structure is violated
732 : ! **************************************************************************************************
733 160 : SUBROUTINE localizer_aii_skew(localizer, skew, tolerance, residual, status)
734 : COMPLEX(KIND=dp), INTENT(IN) :: localizer(:, :)
735 : REAL(KIND=dp), INTENT(OUT) :: skew(:, :)
736 : REAL(KIND=dp), INTENT(IN) :: tolerance
737 : REAL(KIND=dp), INTENT(OUT) :: residual
738 : INTEGER, INTENT(OUT) :: status
739 :
740 : COMPLEX(KIND=dp) :: value
741 : INTEGER :: i, j, m, n, pi, pj, si, sj
742 : REAL(KIND=dp) :: antisymmetric, scale, skew_residual
743 :
744 160 : m = SIZE(localizer, 1)
745 160 : status = -1
746 160 : residual = HUGE(1.0_dp)
747 160 : IF (m < 4 .OR. MOD(m, 4) /= 0) RETURN
748 402 : IF (SIZE(localizer, 2) /= m .OR. ANY(SHAPE(skew) /= [m, m])) RETURN
749 134 : IF (.NOT. ieee_is_finite(tolerance)) RETURN
750 134 : IF (tolerance <= 0.0_dp .OR. tolerance >= 1.0_dp) RETURN
751 19923468 : IF (.NOT. ALL(ieee_is_finite(REAL(localizer, dp))) .OR. &
752 : .NOT. ALL(ieee_is_finite(AIMAG(localizer)))) RETURN
753 134 : n = m/4
754 134 : residual = 0.0_dp
755 : ! C is a signed permutation: Q^dagger L Q = (L + CLC + i(CL-LC))/2.
756 24134 : DO j = 1, m
757 24000 : pj = (3 - (j - 1)/n)*n + MODULO(j - 1, n) + 1
758 24000 : sj = MERGE(1, -1, j <= n .OR. j > 3*n)
759 9961734 : DO i = 1, m
760 9937600 : pi = (3 - (i - 1)/n)*n + MODULO(i - 1, n) + 1
761 9937600 : si = MERGE(1, -1, i <= n .OR. i > 3*n)
762 : value = 0.5_dp*(localizer(i, j) + si*sj*localizer(pi, pj) + &
763 9937600 : CMPLX(0.0_dp, 1.0_dp, dp)*(si*localizer(pi, j) - sj*localizer(i, pj)))
764 9937600 : IF (.NOT. ieee_is_finite(REAL(value, dp)) .OR. .NOT. ieee_is_finite(AIMAG(value))) THEN
765 0 : residual = HUGE(1.0_dp)
766 0 : RETURN
767 : END IF
768 9937600 : residual = MAX(residual, ABS(REAL(value, dp)))
769 9961600 : skew(i, j) = -AIMAG(value)
770 : END DO
771 : END DO
772 9961734 : scale = MAX(TINY(1.0_dp), MAXVAL(ABS(localizer)))
773 134 : IF (residual > tolerance*scale) RETURN
774 : skew_residual = 0.0_dp
775 23380 : DO j = 1, m
776 4990712 : DO i = 1, j
777 4990604 : skew_residual = MAX(skew_residual, ABS(0.5_dp*skew(i, j) + 0.5_dp*skew(j, i)))
778 : END DO
779 : END DO
780 108 : IF (skew_residual > 0.5_dp*tolerance*scale) RETURN
781 : ! Project only validated roundoff; callers include this bound in the gap test.
782 108 : residual = residual + skew_residual
783 23380 : DO j = 1, m
784 4990712 : DO i = 1, j
785 4967332 : antisymmetric = 0.5_dp*skew(i, j) - 0.5_dp*skew(j, i)
786 4967332 : skew(i, j) = antisymmetric
787 4990604 : skew(j, i) = -antisymmetric
788 : END DO
789 : END DO
790 108 : status = 0
791 : END SUBROUTINE localizer_aii_skew
792 :
793 : ! **************************************************************************************************
794 : !> \brief Dense pivoted skew elimination, accumulating only the sign of the Pfaffian.
795 : !> \param matrix Real skew-symmetric matrix; input is not modified
796 : !> \param pf_sign Sign, or zero if singular/unresolved
797 : !> \param status Zero on success, -1 invalid data, -2 numerically singular factorization
798 : !> \note Each simultaneous row/column interchange changes the Pfaffian sign, not its square.
799 : ! **************************************************************************************************
800 704 : SUBROUTINE localizer_pfaffian_sign(matrix, pf_sign, status)
801 : REAL(KIND=dp), INTENT(IN) :: matrix(:, :)
802 : INTEGER, INTENT(OUT) :: pf_sign, status
803 :
804 : INTEGER :: i, j, k, n, p, sign_accum
805 : REAL(KIND=dp) :: pivot, scale
806 704 : REAL(KIND=dp), ALLOCATABLE :: a(:, :), row(:)
807 :
808 704 : n = SIZE(matrix, 1)
809 704 : pf_sign = 0
810 704 : status = -1
811 704 : IF (n < 2 .OR. MOD(n, 2) /= 0 .OR. SIZE(matrix, 2) /= n) RETURN
812 9904056 : IF (.NOT. ALL(ieee_is_finite(matrix))) RETURN
813 9904056 : scale = MAXVAL(ABS(matrix))
814 9904056 : IF (MAXVAL(ABS(matrix + TRANSPOSE(matrix))) > 100.0_dp*EPSILON(1.0_dp)*scale) RETURN
815 704 : status = -2
816 704 : IF (scale <= TINY(1.0_dp)) RETURN
817 4212 : ALLOCATE (a(n, n), row(n))
818 9904014 : a(:, :) = matrix/scale
819 702 : sign_accum = 1
820 13220 : DO k = 1, n - 1, 2
821 2494608 : p = k + MAXLOC(ABS(a(k, k + 1:)), DIM=1)
822 12520 : IF (p /= k + 1) THEN
823 2807856 : row(:) = a(k + 1, :)
824 2807856 : a(k + 1, :) = a(p, :)
825 2807856 : a(p, :) = row
826 2807856 : row(:) = a(:, k + 1)
827 2807856 : a(:, k + 1) = a(:, p)
828 2807856 : a(:, p) = row
829 6256 : sign_accum = -sign_accum
830 : END IF
831 12520 : pivot = a(k, k + 1)
832 773628416 : IF (ABS(pivot) <= 100.0_dp*EPSILON(1.0_dp)*MAXVAL(ABS(a(k:, k:)))) RETURN
833 12518 : IF (pivot < 0.0_dp) sign_accum = -sign_accum
834 2469566 : DO j = k + 2, n
835 381868810 : DO i = k + 2, j - 1
836 379399244 : a(i, j) = a(i, j) + (a(k + 1, i)*a(k, j) - a(k, i)*a(k + 1, j))/pivot
837 381856292 : a(j, i) = -a(i, j)
838 : END DO
839 : END DO
840 773629102 : IF (.NOT. ALL(ieee_is_finite(a(k:, k:)))) RETURN
841 : END DO
842 700 : pf_sign = sign_accum
843 700 : status = 0
844 704 : END SUBROUTINE localizer_pfaffian_sign
845 :
846 : END MODULE spectral_localizer
|
