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 : INTEGER :: info, n, negative, pf, pf_reference, &
431 : positive, ldl_positive, ldl_negative
432 : #if defined(__HAS_IEEE_EXCEPTIONS)
433 : LOGICAL :: halt(SIZE(IEEE_ALL))
434 : #endif
435 : REAL(KIND=dp) :: metric_min, residual, threshold
436 251 : REAL(KIND=dp), ALLOCATABLE :: eigenvalues(:), rwork(:), skew(:, :)
437 :
438 251 : status = -1
439 251 : invariant = 0
440 251 : gap = 0.0_dp
441 251 : metric_ratio = 0.0_dp
442 251 : n = SIZE(s, 1)
443 251 : IF (n < 1 .OR. SIZE(s, 2) /= n) RETURN
444 753 : IF (ANY(SHAPE(localizer) /= [2*n, 2*n])) RETURN
445 1004 : IF (.NOT. ALL(ieee_is_finite([matrix_tol, metric_tol, gap_tol]))) RETURN
446 251 : IF (matrix_tol <= 0.0_dp .OR. matrix_tol >= 1.0_dp) RETURN
447 251 : IF (metric_tol <= 0.0_dp .OR. metric_tol >= 1.0_dp .OR. gap_tol < 0.0_dp) RETURN
448 251 : IF (.NOT. hermitian(s, matrix_tol) .OR. .NOT. hermitian(localizer, matrix_tol)) RETURN
449 251 : IF (aii) THEN
450 30 : status = -3
451 30 : IF (MOD(n, 2) /= 0) RETURN
452 30 : IF (.NOT. self_dual(s, matrix_tol)) RETURN
453 : END IF
454 :
455 2510 : ALLOCATE (a(n, n), eigenvalues(n), work(MAX(1, 2*n)), rwork(MAX(1, 3*n)))
456 6267843 : a(:, :) = s
457 : ! Follow CP2K's eigensolvers: LAPACK may use exceptional IEEE intermediates internally.
458 : #if defined(__HAS_IEEE_EXCEPTIONS)
459 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
460 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
461 : #endif
462 251 : CALL zheev('N', 'U', n, a, n, eigenvalues, work, SIZE(work), rwork, info)
463 : #if defined(__HAS_IEEE_EXCEPTIONS)
464 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
465 : #endif
466 251 : status = ABS(info) + 2
467 32031 : IF (info /= 0 .OR. .NOT. ALL(ieee_is_finite(eigenvalues))) RETURN
468 251 : status = -2
469 32031 : IF (MAXVAL(eigenvalues) <= 0.0_dp) RETURN
470 63811 : metric_ratio = MINVAL(eigenvalues)/MAXVAL(eigenvalues)
471 251 : IF (metric_ratio <= metric_tol) RETURN
472 31175 : metric_min = MINVAL(eigenvalues)
473 243 : DEALLOCATE (a, eigenvalues, work, rwork)
474 :
475 1458 : ALLOCATE (a(2*n, 2*n), b(2*n, 2*n))
476 1701 : ALLOCATE (eigenvalues(2*n), work(4*n), rwork(6*n))
477 243 : status = -1
478 48954582 : IF (.NOT. ALL(ieee_is_finite(REAL(localizer, dp))) .OR. &
479 : .NOT. ALL(ieee_is_finite(AIMAG(localizer)))) RETURN
480 24477291 : a(:, :) = localizer
481 243 : b(:, :) = CMPLX(0.0_dp, 0.0_dp, dp)
482 6134971 : b(:n, :n) = s
483 6134971 : b(n + 1:, n + 1:) = s
484 : ! Generalized eigenvalues, unlike eigenvalues of raw AO matrices, are basis invariant.
485 : #if defined(__HAS_IEEE_EXCEPTIONS)
486 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
487 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
488 : #endif
489 243 : CALL zhegv(1, 'N', 'U', 2*n, a, 2*n, b, 2*n, eigenvalues, work, SIZE(work), rwork, info)
490 : #if defined(__HAS_IEEE_EXCEPTIONS)
491 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
492 : #endif
493 243 : status = ABS(info) + 2
494 62107 : IF (info /= 0 .OR. .NOT. ALL(ieee_is_finite(eigenvalues))) RETURN
495 62107 : gap = MINVAL(ABS(eigenvalues))
496 : ! Bound roundoff and triangle inconsistencies after conversion to the orthonormal metric.
497 : residual = REAL(2*n, dp)*MAXVAL(ABS(localizer - CONJG(TRANSPOSE(localizer))))/metric_min + &
498 30674126 : REAL(n, dp)*MAXVAL(ABS(eigenvalues))*MAXVAL(ABS(s - CONJG(TRANSPOSE(s))))/metric_min
499 : threshold = MAX(gap_tol, residual + 100.0_dp*EPSILON(1.0_dp)*REAL(2*n, dp)* &
500 24539155 : MAX(MAXVAL(ABS(eigenvalues)), MAXVAL(ABS(localizer))/metric_min))
501 243 : status = 1
502 243 : IF (gap <= threshold) RETURN
503 57685 : positive = COUNT(eigenvalues > threshold)
504 57685 : negative = COUNT(eigenvalues < -threshold)
505 209 : CALL localizer_inertia(localizer, ldl_positive, ldl_negative, info)
506 209 : status = -4
507 209 : IF (info /= 0) RETURN
508 209 : IF (positive /= ldl_positive .OR. negative /= ldl_negative) RETURN
509 209 : IF (.NOT. aii) THEN
510 181 : invariant = (ldl_positive - ldl_negative)/2
511 : ELSE
512 28 : status = -3
513 28 : IF (positive /= negative) RETURN
514 112 : ALLOCATE (skew(2*n, 2*n))
515 28 : CALL localizer_aii_skew(localizer, skew, matrix_tol, residual, info)
516 28 : IF (info /= 0) RETURN
517 28 : status = 1
518 28 : IF (gap <= threshold + REAL(2*n, dp)*residual/metric_min) RETURN
519 28 : CALL localizer_pfaffian_sign(skew, pf, info)
520 28 : status = -4
521 28 : IF (info /= 0) RETURN
522 : ! Fix the orientation against the positive atomic Hamiltonian in the SAME AO metric.
523 28 : a(:, :) = CMPLX(0.0_dp, 0.0_dp, dp)
524 1237854 : a(:n, :n) = s
525 1237854 : a(n + 1:, n + 1:) = -s
526 28 : CALL localizer_aii_skew(a, skew, matrix_tol, residual, info)
527 28 : IF (info /= 0) RETURN
528 28 : IF (REAL(2*n, dp)*residual >= metric_min) RETURN
529 28 : CALL localizer_pfaffian_sign(skew, pf_reference, info)
530 28 : IF (info /= 0) RETURN
531 112 : invariant = (1 - pf*pf_reference)/2
532 : END IF
533 209 : status = 0
534 251 : END SUBROUTINE localizer_matrix_evaluate
535 :
536 : ! **************************************************************************************************
537 : !> \brief Count inertia by pivoted Hermitian LDL^dagger, without orthogonalizing the AO basis.
538 : !> \param matrix Hermitian matrix; input is not modified
539 : !> \param positive Number of positive eigenvalues
540 : !> \param negative Number of negative eigenvalues
541 : !> \param status 0 resolved, 1 singular, -1 invalid data, -4 factorization failure
542 : !> \note Pivot magnitudes are not eigenvalues or a metric-aware protection gap.
543 : ! **************************************************************************************************
544 1303 : SUBROUTINE localizer_inertia_complex(matrix, positive, negative, status)
545 : COMPLEX(KIND=dp), INTENT(IN) :: matrix(:, :)
546 : INTEGER, INTENT(OUT) :: positive, negative, status
547 :
548 1303 : COMPLEX(KIND=dp), ALLOCATABLE :: a(:, :), work(:)
549 : INTEGER :: i, info, n
550 1303 : INTEGER, ALLOCATABLE :: pivots(:)
551 : REAL(KIND=dp) :: scale
552 : #if defined(__HAS_IEEE_EXCEPTIONS)
553 : LOGICAL :: halt(SIZE(IEEE_ALL))
554 : #endif
555 :
556 1303 : positive = 0
557 1303 : negative = 0
558 1303 : status = -1
559 1303 : n = SIZE(matrix, 1)
560 1303 : IF (n < 1 .OR. SIZE(matrix, 2) /= n) RETURN
561 1303 : IF (.NOT. hermitian(matrix, 1.e-10_dp)) RETURN
562 : ! Never discard a small imaginary term: only exactly real matrices use LDL^T.
563 133123 : IF (ALL(AIMAG(matrix) == 0.0_dp)) THEN
564 112556 : CALL localizer_inertia_real(REAL(matrix, dp), positive, negative, status)
565 568 : RETURN
566 : END IF
567 22933613 : scale = MAXVAL(ABS(matrix))
568 0 : status = 1
569 733 : IF (scale <= TINY(1.0_dp)) RETURN
570 5864 : ALLOCATE (a(n, n), pivots(n), work(MAX(1, 64*n)))
571 22933613 : a(:, :) = matrix/scale
572 : #if defined(__HAS_IEEE_EXCEPTIONS)
573 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
574 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
575 : #endif
576 733 : CALL zhetrf('L', n, a, n, pivots, work, SIZE(work), info)
577 : #if defined(__HAS_IEEE_EXCEPTIONS)
578 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
579 : #endif
580 733 : IF (info > 0) RETURN
581 731 : status = -4
582 731 : IF (info < 0) RETURN
583 45867142 : IF (.NOT. ALL(ieee_is_finite(REAL(a, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(a)))) RETURN
584 : CALL block_inertia([(REAL(a(i, i), dp), i=1, n)], [(ABS(a(i + 1, i)), i=1, n - 1)], &
585 244357 : pivots, positive, negative, status)
586 1303 : END SUBROUTINE localizer_inertia_complex
587 :
588 : ! **************************************************************************************************
589 : !> \brief Native real symmetric LDL^T inertia, without converting to complex storage.
590 : !> \param matrix Real symmetric matrix; input is not modified
591 : !> \param positive Number of positive eigenvalues
592 : !> \param negative Number of negative eigenvalues
593 : !> \param status 0 resolved, 1 singular, -1 invalid data, -4 factorization failure
594 : ! **************************************************************************************************
595 1114 : SUBROUTINE localizer_inertia_real(matrix, positive, negative, status)
596 : REAL(dp), INTENT(IN) :: matrix(:, :)
597 : INTEGER, INTENT(OUT) :: positive, negative, status
598 1114 : REAL(dp), ALLOCATABLE :: a(:, :), work(:)
599 1114 : INTEGER, ALLOCATABLE :: pivots(:)
600 : INTEGER :: i, info, n
601 : REAL(dp) :: scale
602 : #if defined(__HAS_IEEE_EXCEPTIONS)
603 : LOGICAL :: halt(SIZE(IEEE_ALL))
604 : #endif
605 :
606 1114 : positive = 0
607 1114 : negative = 0
608 1114 : status = -1
609 1114 : n = SIZE(matrix, 1)
610 1114 : IF (n < 1 .OR. SIZE(matrix, 2) /= n) RETURN
611 157910 : IF (.NOT. ALL(ieee_is_finite(matrix))) RETURN
612 157910 : scale = MAXVAL(ABS(matrix))
613 157910 : IF (MAXVAL(ABS(matrix - TRANSPOSE(matrix))) > 1.e-10_dp*MAX(TINY(1.0_dp), scale)) RETURN
614 1112 : status = 1
615 1112 : IF (scale <= TINY(1.0_dp)) RETURN
616 8880 : ALLOCATE (a(n, n), pivots(n), work(MAX(1, 64*n)))
617 157882 : a(:, :) = matrix/scale
618 : #if defined(__HAS_IEEE_EXCEPTIONS)
619 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
620 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
621 : #endif
622 1110 : CALL dsytrf('L', n, a, n, pivots, work, SIZE(work), info)
623 : #if defined(__HAS_IEEE_EXCEPTIONS)
624 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
625 : #endif
626 1110 : IF (info > 0) RETURN
627 1110 : status = -4
628 157882 : IF (info < 0 .OR. .NOT. ALL(ieee_is_finite(a))) RETURN
629 : CALL block_inertia([(a(i, i), i=1, n)], [(ABS(a(i + 1, i)), i=1, n - 1)], &
630 36298 : pivots, positive, negative, status)
631 1114 : END SUBROUTINE localizer_inertia_real
632 :
633 : ! **************************************************************************************************
634 : !> \brief Count signs of lower Bunch-Kaufman D blocks, shared by real and complex factors.
635 : !> \param diagonal Real diagonal of the factored matrix
636 : !> \param coupling Magnitudes of first subdiagonal entries
637 : !> \param pivots LAPACK lower-triangle pivot encoding
638 : !> \param positive Positive inertia count
639 : !> \param negative Negative inertia count
640 : !> \param status 0 resolved, 1 singular, -4 invalid block structure
641 : !> \note A Hermitian 2x2 block has the same eigenvalues when its offdiagonal is replaced by its modulus.
642 : ! **************************************************************************************************
643 1841 : SUBROUTINE block_inertia(diagonal, coupling, pivots, positive, negative, status)
644 : REAL(dp), INTENT(IN) :: diagonal(:), coupling(:)
645 : INTEGER, INTENT(IN) :: pivots(:)
646 : INTEGER, INTENT(OUT) :: positive, negative, status
647 :
648 : INTEGER :: i, n
649 : REAL(dp) :: values(2)
650 : INTERFACE
651 : SUBROUTINE dlae2(a, b, c, rt1, rt2)
652 : IMPORT :: dp
653 : REAL(dp), INTENT(IN) :: a, b, c
654 : REAL(dp), INTENT(OUT) :: rt1, rt2
655 : END SUBROUTINE dlae2
656 : END INTERFACE
657 :
658 1841 : positive = 0
659 1841 : negative = 0
660 1841 : status = -4
661 1841 : n = SIZE(pivots)
662 1841 : i = 1
663 52590 : DO WHILE (i <= n)
664 50749 : IF (pivots(i) == 0) RETURN
665 50749 : IF (pivots(i) > 0) THEN
666 92622 : values = [diagonal(i), 0.0_dp]
667 30874 : IF (values(1) == 0.0_dp) THEN
668 0 : status = 1
669 0 : RETURN
670 : END IF
671 30874 : i = i + 1
672 : ELSE
673 19875 : IF (i == n) RETURN
674 19875 : IF (pivots(i + 1) /= pivots(i)) RETURN
675 : ! A 2x2 pivot can have zero diagonal and still contribute one sign of each kind.
676 19875 : CALL dlae2(diagonal(i), coupling(i), diagonal(i + 1), values(1), values(2))
677 59625 : IF (ANY(values == 0.0_dp)) THEN
678 0 : status = 1
679 0 : RETURN
680 : END IF
681 19875 : i = i + 2
682 : END IF
683 152247 : positive = positive + COUNT(values > 0.0_dp)
684 154088 : negative = negative + COUNT(values < 0.0_dp)
685 : END DO
686 1841 : status = 0
687 : END SUBROUTINE block_inertia
688 :
689 : ! **************************************************************************************************
690 : !> \brief Test finiteness and Hermiticity before passing only one triangle to LAPACK.
691 : !> \param a Square matrix
692 : !> \param tolerance Relative tolerance
693 : !> \return Whether the matrix is Hermitian and finite
694 : ! **************************************************************************************************
695 2846 : LOGICAL FUNCTION hermitian(a, tolerance)
696 : COMPLEX(KIND=dp), INTENT(IN) :: a(:, :)
697 : REAL(KIND=dp), INTENT(IN) :: tolerance
698 :
699 2846 : hermitian = .FALSE.
700 142566540 : IF (.NOT. ALL(ieee_is_finite(REAL(a, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(a)))) RETURN
701 142542992 : hermitian = MAXVAL(ABS(a - CONJG(TRANSPOSE(a)))) <= tolerance*MAX(TINY(1.0_dp), MAXVAL(ABS(a)))
702 2844 : END FUNCTION hermitian
703 :
704 : ! **************************************************************************************************
705 : !> \brief Check J conj(A) J^dagger=A for J=(0,I;-I,0), without assuming spin conservation.
706 : !> \param a Even-dimensional matrix in the fixed up/down ordering
707 : !> \param tolerance Relative tolerance
708 : !> \return Whether physical fermionic time reversal is respected
709 : ! **************************************************************************************************
710 657 : LOGICAL FUNCTION self_dual(a, tolerance)
711 : COMPLEX(KIND=dp), INTENT(IN) :: a(:, :)
712 : REAL(KIND=dp), INTENT(IN) :: tolerance
713 :
714 : INTEGER :: n
715 : REAL(KIND=dp) :: scale
716 :
717 657 : n = SIZE(a, 1)/2
718 13202191 : scale = MAX(TINY(1.0_dp), MAXVAL(ABS(a)))
719 : self_dual = MAXVAL(ABS(a(:n, :n) - CONJG(a(n + 1:, n + 1:)))) <= tolerance*scale .AND. &
720 6644430 : MAXVAL(ABS(a(:n, n + 1:) + CONJG(a(n + 1:, :n)))) <= tolerance*scale
721 657 : END FUNCTION self_dual
722 :
723 : ! **************************************************************************************************
724 : !> \brief Convert the class AII localizer to a real skew matrix using Q=(I-i C)/sqrt(2).
725 : !> \param localizer Hermitian localizer; auxiliary index outside physical up/down index
726 : !> \param skew Real part of i Q^dagger L Q
727 : !> \param tolerance Relative imaginary/skew residual tolerance
728 : !> \param residual Entrywise bound for the discarded real and symmetric components
729 : !> \param status Zero on success; negative if the AII structure is violated
730 : ! **************************************************************************************************
731 160 : SUBROUTINE localizer_aii_skew(localizer, skew, tolerance, residual, status)
732 : COMPLEX(KIND=dp), INTENT(IN) :: localizer(:, :)
733 : REAL(KIND=dp), INTENT(OUT) :: skew(:, :)
734 : REAL(KIND=dp), INTENT(IN) :: tolerance
735 : REAL(KIND=dp), INTENT(OUT) :: residual
736 : INTEGER, INTENT(OUT) :: status
737 :
738 : COMPLEX(KIND=dp) :: value
739 : INTEGER :: i, j, m, n, pi, pj, si, sj
740 : REAL(KIND=dp) :: antisymmetric, scale, skew_residual
741 :
742 160 : m = SIZE(localizer, 1)
743 160 : status = -1
744 160 : residual = HUGE(1.0_dp)
745 160 : IF (m < 4 .OR. MOD(m, 4) /= 0) RETURN
746 402 : IF (SIZE(localizer, 2) /= m .OR. ANY(SHAPE(skew) /= [m, m])) RETURN
747 134 : IF (.NOT. ieee_is_finite(tolerance)) RETURN
748 134 : IF (tolerance <= 0.0_dp .OR. tolerance >= 1.0_dp) RETURN
749 19923468 : IF (.NOT. ALL(ieee_is_finite(REAL(localizer, dp))) .OR. &
750 : .NOT. ALL(ieee_is_finite(AIMAG(localizer)))) RETURN
751 134 : n = m/4
752 134 : residual = 0.0_dp
753 : ! C is a signed permutation: Q^dagger L Q = (L + CLC + i(CL-LC))/2.
754 24134 : DO j = 1, m
755 24000 : pj = (3 - (j - 1)/n)*n + MODULO(j - 1, n) + 1
756 24000 : sj = MERGE(1, -1, j <= n .OR. j > 3*n)
757 9961734 : DO i = 1, m
758 9937600 : pi = (3 - (i - 1)/n)*n + MODULO(i - 1, n) + 1
759 9937600 : si = MERGE(1, -1, i <= n .OR. i > 3*n)
760 : value = 0.5_dp*(localizer(i, j) + si*sj*localizer(pi, pj) + &
761 9937600 : CMPLX(0.0_dp, 1.0_dp, dp)*(si*localizer(pi, j) - sj*localizer(i, pj)))
762 9937600 : IF (.NOT. ieee_is_finite(REAL(value, dp)) .OR. .NOT. ieee_is_finite(AIMAG(value))) THEN
763 0 : residual = HUGE(1.0_dp)
764 0 : RETURN
765 : END IF
766 9937600 : residual = MAX(residual, ABS(REAL(value, dp)))
767 9961600 : skew(i, j) = -AIMAG(value)
768 : END DO
769 : END DO
770 9961734 : scale = MAX(TINY(1.0_dp), MAXVAL(ABS(localizer)))
771 134 : IF (residual > tolerance*scale) RETURN
772 : skew_residual = 0.0_dp
773 23380 : DO j = 1, m
774 4990712 : DO i = 1, j
775 4990604 : skew_residual = MAX(skew_residual, ABS(0.5_dp*skew(i, j) + 0.5_dp*skew(j, i)))
776 : END DO
777 : END DO
778 108 : IF (skew_residual > 0.5_dp*tolerance*scale) RETURN
779 : ! Project only validated roundoff; callers include this bound in the gap test.
780 108 : residual = residual + skew_residual
781 23380 : DO j = 1, m
782 4990712 : DO i = 1, j
783 4967332 : antisymmetric = 0.5_dp*skew(i, j) - 0.5_dp*skew(j, i)
784 4967332 : skew(i, j) = antisymmetric
785 4990604 : skew(j, i) = -antisymmetric
786 : END DO
787 : END DO
788 108 : status = 0
789 : END SUBROUTINE localizer_aii_skew
790 :
791 : ! **************************************************************************************************
792 : !> \brief Dense pivoted skew elimination, accumulating only the sign of the Pfaffian.
793 : !> \param matrix Real skew-symmetric matrix; input is not modified
794 : !> \param pf_sign Sign, or zero if singular/unresolved
795 : !> \param status Zero on success, -1 invalid data, -2 numerically singular factorization
796 : !> \note Each simultaneous row/column interchange changes the Pfaffian sign, not its square.
797 : ! **************************************************************************************************
798 704 : SUBROUTINE localizer_pfaffian_sign(matrix, pf_sign, status)
799 : REAL(KIND=dp), INTENT(IN) :: matrix(:, :)
800 : INTEGER, INTENT(OUT) :: pf_sign, status
801 :
802 : INTEGER :: i, j, k, n, p, sign_accum
803 : REAL(KIND=dp) :: pivot, scale
804 704 : REAL(KIND=dp), ALLOCATABLE :: a(:, :), row(:)
805 :
806 704 : n = SIZE(matrix, 1)
807 704 : pf_sign = 0
808 704 : status = -1
809 704 : IF (n < 2 .OR. MOD(n, 2) /= 0 .OR. SIZE(matrix, 2) /= n) RETURN
810 9904056 : IF (.NOT. ALL(ieee_is_finite(matrix))) RETURN
811 9904056 : scale = MAXVAL(ABS(matrix))
812 9904056 : IF (MAXVAL(ABS(matrix + TRANSPOSE(matrix))) > 100.0_dp*EPSILON(1.0_dp)*scale) RETURN
813 704 : status = -2
814 704 : IF (scale <= TINY(1.0_dp)) RETURN
815 4212 : ALLOCATE (a(n, n), row(n))
816 9904014 : a(:, :) = matrix/scale
817 702 : sign_accum = 1
818 13220 : DO k = 1, n - 1, 2
819 2494608 : p = k + MAXLOC(ABS(a(k, k + 1:)), DIM=1)
820 12520 : IF (p /= k + 1) THEN
821 2807856 : row(:) = a(k + 1, :)
822 2807856 : a(k + 1, :) = a(p, :)
823 2807856 : a(p, :) = row
824 2807856 : row(:) = a(:, k + 1)
825 2807856 : a(:, k + 1) = a(:, p)
826 2807856 : a(:, p) = row
827 6256 : sign_accum = -sign_accum
828 : END IF
829 12520 : pivot = a(k, k + 1)
830 773628416 : IF (ABS(pivot) <= 100.0_dp*EPSILON(1.0_dp)*MAXVAL(ABS(a(k:, k:)))) RETURN
831 12518 : IF (pivot < 0.0_dp) sign_accum = -sign_accum
832 2469566 : DO j = k + 2, n
833 381868810 : DO i = k + 2, j - 1
834 379399244 : a(i, j) = a(i, j) + (a(k + 1, i)*a(k, j) - a(k, i)*a(k + 1, j))/pivot
835 381856292 : a(j, i) = -a(i, j)
836 : END DO
837 : END DO
838 773629102 : IF (.NOT. ALL(ieee_is_finite(a(k:, k:)))) RETURN
839 : END DO
840 700 : pf_sign = sign_accum
841 700 : status = 0
842 704 : END SUBROUTINE localizer_pfaffian_sign
843 :
844 : END MODULE spectral_localizer
|
