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 Exact spectral flattening of a translation-invariant AO torus in primitive Bloch blocks.
10 : !> \author Thomas D. Kuehne, tkuehne@cp2k.org
11 : ! **************************************************************************************************
12 : MODULE spectral_localizer_bloch
13 : USE cp_dbcsr_api, ONLY: dbcsr_add, dbcsr_create, dbcsr_distribution_get, &
14 : dbcsr_distribution_type, dbcsr_finalize, dbcsr_get_info, &
15 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
16 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_put_block, dbcsr_type
17 : USE cp_dbcsr_contrib, ONLY: dbcsr_frobenius_norm
18 : USE ieee_arithmetic, ONLY: ieee_is_finite
19 : #if defined(__HAS_IEEE_EXCEPTIONS)
20 : USE ieee_exceptions, ONLY: IEEE_ALL, ieee_get_halting_mode, ieee_set_halting_mode
21 : #endif
22 : USE kinds, ONLY: dp
23 : USE kpoint_lattice_fft, ONLY: k_grid_to_cell_fft
24 : USE message_passing, ONLY: mp_para_env_type
25 : USE spectral_localizer_dbcsr, ONLY: localizer_dbcsr_check, localizer_dbcsr_release, localizer_matrix_type
26 : #include "./base/base_uses.f90"
27 :
28 : IMPLICIT NONE
29 : PRIVATE
30 : PUBLIC :: localizer_flatten_bloch
31 :
32 : CONTAINS
33 :
34 : ! **************************************************************************************************
35 : !> \brief Flatten all Bloch bands, then return to the original distributed real-cell/spin AO basis.
36 : !> \param h Original complex Hamiltonian in spin/cell/atom order, unchanged
37 : !> \param s Original complex metric in the same order, unchanged
38 : !> \param grid Full Gamma-centered torus dimensions
39 : !> \param nspin Number of physical spin blocks (one or two)
40 : !> \param energy Absolute query energy
41 : !> \param scale Positive flattened energy
42 : !> \param tolerance Structure, covariance and eigensystem tolerance
43 : !> \param metric_tolerance Minimum global overlap eigenvalue ratio
44 : !> \param gap_tolerance Minimum original electronic distance from the query
45 : !> \param para_env Communicator; independent eigensystems are distributed across its ranks
46 : !> \param flattened Newly allocated covariant shifted Hamiltonian only on success
47 : !> \param gap Original electronic gap, not a flattened localizer gap
48 : !> \param metric_ratio Global minimum/maximum overlap eigenvalue ratio across all k-points
49 : !> \param residuals Covariance, metric orthonormality, eigenproblem and final Hermiticity residuals
50 : !> \param status 0 success, 1 unresolved gap, -1 invalid layout/data, -2 metric,
51 : !> -4 eigensolver, -6 excessive residual, -7 broken translation covariance
52 : !> \note No occupied-state projection, symmetry reduction or modification of coordinate operators.
53 : ! **************************************************************************************************
54 62 : SUBROUTINE localizer_flatten_bloch(h, s, grid, nspin, energy, scale, tolerance, metric_tolerance, &
55 : gap_tolerance, para_env, flattened, gap, metric_ratio, residuals, status)
56 : TYPE(localizer_matrix_type), INTENT(INOUT) :: h, s
57 : INTEGER, INTENT(IN) :: grid(3), nspin
58 : REAL(dp), INTENT(IN) :: energy, scale, tolerance, &
59 : metric_tolerance, gap_tolerance
60 : TYPE(mp_para_env_type), POINTER :: para_env
61 : TYPE(localizer_matrix_type), INTENT(OUT) :: flattened
62 : REAL(dp), INTENT(OUT) :: gap, metric_ratio, residuals(4)
63 : INTEGER, INTENT(OUT) :: status
64 :
65 62 : COMPLEX(dp), ALLOCATABLE :: fk(:, :, :), hc(:, :, :), hk(:, :, :), &
66 62 : sc(:, :, :), sk(:, :, :)
67 : INTEGER :: ik, nb, ncell, ns, nscalar, nsmall, &
68 : substatus
69 62 : INTEGER, ALLOCATABLE :: cell_index(:), cells(:, :), starts(:)
70 62 : INTEGER, POINTER :: sizes(:), ssize(:)
71 : LOGICAL :: valid
72 : REAL(dp) :: all_bounds(2), bounds(2), difference, &
73 : errors(2), gap_threshold, local_gap, &
74 : threshold
75 62 : REAL(dp), ALLOCATABLE :: weights(:), xkp(:, :)
76 : TYPE(localizer_matrix_type) :: reconstructed
77 :
78 62 : status = -1
79 62 : gap = 0.0_dp
80 62 : metric_ratio = 0.0_dp
81 310 : residuals = HUGE(1.0_dp)
82 242 : IF (ANY(grid < 1) .OR. (nspin /= 1 .AND. nspin /= 2)) RETURN
83 360 : IF (.NOT. ALL(ieee_is_finite([energy, scale, tolerance, metric_tolerance, gap_tolerance]))) RETURN
84 60 : IF (scale <= 0.0_dp .OR. tolerance <= 0.0_dp .OR. tolerance >= 1.0_dp) RETURN
85 60 : IF (metric_tolerance <= 0.0_dp .OR. metric_tolerance >= 1.0_dp .OR. gap_tolerance < 0.0_dp) RETURN
86 60 : CALL dbcsr_get_info(h%real_part, nblkrows_total=nb, row_blk_size=sizes)
87 60 : CALL dbcsr_get_info(s%real_part, nblkrows_total=ns, row_blk_size=ssize)
88 60 : IF (nb < 1 .OR. nb /= ns) RETURN
89 460 : IF (ANY(sizes /= ssize)) RETURN
90 460 : IF (ANY(sizes <= 0)) RETURN
91 : ! Bound the products before interpreting the repeated atom-block layout.
92 : ncell = 1
93 240 : DO ik = 1, 3
94 180 : IF (grid(ik) > nb/ncell) RETURN
95 240 : ncell = ncell*grid(ik)
96 : END DO
97 60 : IF (nspin > nb/ncell) RETURN
98 60 : IF (MOD(nb, nspin*ncell) /= 0) RETURN
99 60 : ns = nb/(nspin*ncell)
100 460 : DO ik = 1, nb
101 460 : IF (sizes(ik) /= sizes(1 + MODULO(ik - 1, ns))) RETURN
102 : END DO
103 60 : CALL localizer_dbcsr_check(h, para_env, tolerance, valid)
104 60 : IF (.NOT. valid) RETURN
105 60 : CALL localizer_dbcsr_check(s, para_env, tolerance, valid)
106 60 : IF (.NOT. valid) RETURN
107 146 : nscalar = SUM(sizes(:ns))
108 60 : nsmall = nscalar*nspin
109 600 : ALLOCATE (starts(nb), cell_index(nb), cells(3, ncell), xkp(3, ncell), weights(ncell))
110 460 : DO ik = 1, nb
111 488 : starts(ik) = 1 + SUM(sizes(:MODULO(ik - 1, ns))) + ((ik - 1)/(ns*ncell))*nscalar
112 460 : cell_index(ik) = 1 + MODULO((ik - 1)/ns, ncell)
113 : END DO
114 262 : DO ik = 1, ncell
115 : cells(:, ik) = [MODULO(ik - 1, grid(1)), MODULO((ik - 1)/grid(1), grid(2)), &
116 808 : (ik - 1)/(grid(1)*grid(2))]
117 868 : xkp(:, ik) = REAL(cells(:, ik), dp)/REAL(grid, dp)
118 : END DO
119 : ALLOCATE (hc(nsmall, nsmall, ncell), sc(nsmall, nsmall, ncell), &
120 1020 : hk(nsmall, nsmall, ncell), sk(nsmall, nsmall, ncell), fk(nsmall, nsmall, ncell))
121 60 : CALL torus_blocks(h, starts, cell_index, para_env, hc)
122 60 : CALL torus_blocks(s, starts, cell_index, para_env, sc)
123 60 : residuals = 0.0_dp
124 : ! Check the complete matrix, including missing blocks, before using one cell as representative.
125 60 : CALL expand_blocks(hc, h, grid, starts, cell_index, para_env, reconstructed)
126 60 : CALL matrix_difference(reconstructed, h, difference)
127 60 : CALL localizer_dbcsr_release(reconstructed)
128 60 : residuals(1) = difference
129 60 : CALL expand_blocks(sc, s, grid, starts, cell_index, para_env, reconstructed)
130 60 : CALL matrix_difference(reconstructed, s, difference)
131 60 : CALL localizer_dbcsr_release(reconstructed)
132 60 : residuals(1) = MAX(residuals(1), difference)
133 60 : status = -7
134 60 : IF (.NOT. ieee_is_finite(residuals(1)) .OR. residuals(1) > tolerance) RETURN
135 : ! B(R) -> B(k) has positive phase and no normalization; the inverse has negative phase / Nk.
136 208 : weights = 1.0_dp
137 688 : CALL k_grid_to_cell_fft(hc, -xkp, grid, cells, hk, weights=weights)
138 688 : CALL k_grid_to_cell_fft(sc, -xkp, grid, cells, sk, weights=weights)
139 48 : fk = 0.0_dp
140 48 : all_bounds = [HUGE(1.0_dp), 0.0_dp]
141 48 : gap = HUGE(1.0_dp)
142 48 : gap_threshold = gap_tolerance
143 48 : status = 0
144 208 : DO ik = 1, ncell
145 160 : IF (MODULO(ik - 1, para_env%num_pe) /= para_env%mepos) CYCLE
146 : CALL flatten_block(hk(:, :, ik), sk(:, :, ik), energy, scale, tolerance, metric_tolerance, fk(:, :, ik), &
147 80 : local_gap, bounds, threshold, errors, substatus)
148 80 : IF (substatus /= 0) THEN
149 22 : status = MIN(status, substatus)
150 22 : CYCLE
151 : END IF
152 174 : all_bounds = [MIN(all_bounds(1), bounds(1)), MAX(all_bounds(2), bounds(2))]
153 58 : gap = MIN(gap, local_gap)
154 58 : gap_threshold = MAX(gap_threshold, threshold)
155 302 : residuals(2:3) = MAX(residuals(2:3), errors)
156 : END DO
157 48 : CALL para_env%min(status)
158 48 : IF (status /= 0) RETURN
159 34 : CALL para_env%min(all_bounds(1))
160 34 : CALL para_env%max(all_bounds(2))
161 34 : CALL para_env%min(gap)
162 34 : CALL para_env%max(gap_threshold)
163 34 : CALL para_env%max(residuals)
164 34 : metric_ratio = all_bounds(1)/all_bounds(2)
165 34 : status = -2
166 34 : IF (metric_ratio <= metric_tolerance) RETURN
167 32 : status = 1
168 32 : IF (gap <= gap_threshold) RETURN
169 18 : status = -6
170 162 : IF (MAXVAL(residuals) > tolerance .OR. .NOT. ALL(ieee_is_finite(residuals))) RETURN
171 18 : CALL para_env%sum(fk)
172 18 : CALL k_grid_to_cell_fft(fk, xkp, grid, cells, hc)
173 18 : CALL expand_blocks(hc, h, grid, starts, cell_index, para_env, flattened)
174 18 : CALL localizer_dbcsr_check(flattened, para_env, tolerance, valid, residuals(4))
175 18 : IF (.NOT. valid) THEN
176 0 : CALL localizer_dbcsr_release(flattened)
177 0 : RETURN
178 : END IF
179 18 : status = 0
180 442 : END SUBROUTINE localizer_flatten_bloch
181 :
182 : ! **************************************************************************************************
183 : !> \brief Extract the distributed first bra cell for every physical spin; no full AO gathering.
184 : !> \param matrix Original matrix
185 : !> \param starts Primitive spin/AO offsets for each atomic block
186 : !> \param cell_index Cell of each block
187 : !> \param para_env Communicator
188 : !> \param blocks One primitive spin matrix for every cell displacement
189 : ! **************************************************************************************************
190 120 : SUBROUTINE torus_blocks(matrix, starts, cell_index, para_env, blocks)
191 : TYPE(localizer_matrix_type), INTENT(INOUT) :: matrix
192 : INTEGER, INTENT(IN) :: starts(:), cell_index(:)
193 : TYPE(mp_para_env_type), POINTER :: para_env
194 : COMPLEX(dp), INTENT(OUT) :: blocks(:, :, :)
195 :
196 120 : CALL extract(matrix%real_part, CMPLX(1.0_dp, 0.0_dp, dp), .TRUE.)
197 120 : CALL extract(matrix%imaginary_part, CMPLX(0.0_dp, 1.0_dp, dp), .FALSE.)
198 149216 : CALL para_env%sum(blocks)
199 :
200 : CONTAINS
201 :
202 : ! **************************************************************************************************
203 : !> \brief Accumulate one real component into the complex block row.
204 : !> \param component Real DBCSR component
205 : !> \param factor Real or imaginary unit
206 : !> \param initialize Reset output before accumulation
207 : ! **************************************************************************************************
208 480 : SUBROUTINE extract(component, factor, initialize)
209 : TYPE(dbcsr_type), INTENT(INOUT) :: component
210 : COMPLEX(dp), INTENT(IN) :: factor
211 : LOGICAL, INTENT(IN) :: initialize
212 :
213 : INTEGER :: c, col, i, j, row
214 : LOGICAL :: transposed
215 240 : REAL(dp), POINTER :: block(:, :)
216 : TYPE(dbcsr_iterator_type) :: iterator
217 :
218 74668 : IF (initialize) blocks = 0.0_dp
219 240 : CALL dbcsr_iterator_start(iterator, component)
220 5552 : DO WHILE (dbcsr_iterator_blocks_left(iterator))
221 5312 : CALL dbcsr_iterator_next_block(iterator, row, col, block, transposed=transposed)
222 5312 : CPASSERT(.NOT. transposed)
223 5312 : IF (cell_index(row) /= 1) CYCLE
224 1800 : i = starts(row)
225 1800 : j = starts(col)
226 1800 : c = cell_index(col)
227 : blocks(i:i + SIZE(block, 1) - 1, j:j + SIZE(block, 2) - 1, c) = &
228 58592 : blocks(i:i + SIZE(block, 1) - 1, j:j + SIZE(block, 2) - 1, c) + factor*block
229 : END DO
230 240 : CALL dbcsr_iterator_stop(iterator)
231 240 : END SUBROUTINE extract
232 : END SUBROUTINE torus_blocks
233 :
234 : ! **************************************************************************************************
235 : !> \brief Expand one block row into the original DBCSR distribution without changing spin ordering.
236 : !> \param blocks Primitive spin matrices at each displacement
237 : !> \param template Original full-storage AO layout
238 : !> \param grid Torus dimensions
239 : !> \param starts Primitive AO starts
240 : !> \param cell_index Cell of each atomic block
241 : !> \param para_env Communicator
242 : !> \param matrix Newly allocated full torus matrix
243 : ! **************************************************************************************************
244 138 : SUBROUTINE expand_blocks(blocks, template, grid, starts, cell_index, para_env, matrix)
245 : COMPLEX(dp), INTENT(IN) :: blocks(:, :, :)
246 : TYPE(localizer_matrix_type), INTENT(IN) :: template
247 : INTEGER, INTENT(IN) :: grid(3), starts(:), cell_index(:)
248 : TYPE(mp_para_env_type), POINTER :: para_env
249 : TYPE(localizer_matrix_type), INTENT(OUT) :: matrix
250 :
251 : INTEGER :: row, col, i, j, c, d(3), t, u
252 138 : INTEGER, POINTER :: sizes(:), rows(:), cols(:), pgrid(:, :)
253 : TYPE(dbcsr_distribution_type) :: distribution
254 :
255 138 : CALL dbcsr_get_info(template%real_part, distribution=distribution, row_blk_size=sizes)
256 138 : CALL dbcsr_distribution_get(distribution, row_dist=rows, col_dist=cols, pgrid=pgrid)
257 138 : CALL dbcsr_create(matrix%real_part, template=template%real_part)
258 138 : CALL dbcsr_create(matrix%imaginary_part, template=template%imaginary_part)
259 1086 : DO row = 1, SIZE(starts)
260 948 : t = cell_index(row) - 1
261 948 : i = starts(row)
262 8694 : DO col = 1, SIZE(starts)
263 7608 : IF (pgrid(rows(row), cols(col)) /= para_env%mepos) CYCLE
264 3804 : u = cell_index(col) - 1
265 3804 : j = starts(col)
266 : d = MODULO([MODULO(u, grid(1)) - MODULO(t, grid(1)), &
267 : MODULO(u/grid(1), grid(2)) - MODULO(t/grid(1), grid(2)), &
268 15216 : u/(grid(1)*grid(2)) - t/(grid(1)*grid(2))], grid)
269 3804 : c = 1 + d(1) + grid(1)*(d(2) + grid(2)*d(3))
270 : CALL dbcsr_put_block(matrix%real_part, row, col, &
271 296964 : REAL(blocks(i:i + sizes(row) - 1, j:j + sizes(col) - 1, c), dp))
272 : CALL dbcsr_put_block(matrix%imaginary_part, row, col, &
273 301716 : AIMAG(blocks(i:i + sizes(row) - 1, j:j + sizes(col) - 1, c)))
274 : END DO
275 : END DO
276 138 : CALL dbcsr_finalize(matrix%real_part)
277 138 : CALL dbcsr_finalize(matrix%imaginary_part)
278 138 : END SUBROUTINE expand_blocks
279 :
280 : ! **************************************************************************************************
281 : !> \brief Relative complex Frobenius difference; leaves the trial as the difference.
282 : !> \param trial Reconstructed matrix
283 : !> \param reference Original matrix
284 : !> \param error Relative difference
285 : ! **************************************************************************************************
286 120 : SUBROUTINE matrix_difference(trial, reference, error)
287 : TYPE(localizer_matrix_type), INTENT(INOUT) :: trial, reference
288 : REAL(dp), INTENT(OUT) :: error
289 :
290 : REAL(dp) :: norm
291 :
292 120 : norm = HYPOT(dbcsr_frobenius_norm(reference%real_part), dbcsr_frobenius_norm(reference%imaginary_part))
293 120 : CALL dbcsr_add(trial%real_part, reference%real_part, 1.0_dp, -1.0_dp)
294 120 : CALL dbcsr_add(trial%imaginary_part, reference%imaginary_part, 1.0_dp, -1.0_dp)
295 : error = HYPOT(dbcsr_frobenius_norm(trial%real_part), dbcsr_frobenius_norm(trial%imaginary_part))/ &
296 120 : MAX(norm, TINY(1.0_dp))
297 120 : END SUBROUTINE matrix_difference
298 :
299 : ! **************************************************************************************************
300 : !> \brief Full generalized Bloch eigensystem, with metric/gap and numerical residual diagnostics.
301 : !> \param h Hermitian primitive Bloch Hamiltonian
302 : !> \param s Hermitian positive metric
303 : !> \param energy Absolute query energy
304 : !> \param scale Flattened band energy
305 : !> \param tolerance Structure tolerance
306 : !> \param metric_tolerance Minimum local overlap eigenvalue ratio
307 : !> \param flat Covariant shifted sign operator
308 : !> \param gap Distance of the original spectrum from energy
309 : !> \param bounds Smallest/largest original metric eigenvalues
310 : !> \param threshold Roundoff gap threshold
311 : !> \param errors Metric orthonormality and relative eigenproblem residual
312 : !> \param status 0 success, -1 invalid data, -2 metric, -4 eigensolver
313 : ! **************************************************************************************************
314 80 : SUBROUTINE flatten_block(h, s, energy, scale, tolerance, metric_tolerance, flat, gap, bounds, threshold, errors, status)
315 : COMPLEX(dp), INTENT(IN) :: h(:, :), s(:, :)
316 : REAL(dp), INTENT(IN) :: energy, scale, tolerance, metric_tolerance
317 : COMPLEX(dp), INTENT(OUT) :: flat(:, :)
318 : REAL(dp), INTENT(OUT) :: gap, bounds(2), threshold, errors(2)
319 : INTEGER, INTENT(OUT) :: status
320 :
321 80 : COMPLEX(dp), ALLOCATABLE :: a(:, :), b(:, :), work(:), sc(:, :), residual(:, :)
322 80 : REAL(dp), ALLOCATABLE :: eigenvalues(:), rwork(:)
323 : REAL(dp) :: hnorm, snorm, herror, serror
324 : INTEGER :: n, info, i
325 : #if defined(__HAS_IEEE_EXCEPTIONS)
326 : LOGICAL :: halt(SIZE(IEEE_ALL))
327 : #endif
328 :
329 80 : status = -1
330 80 : n = SIZE(h, 1)
331 36582 : IF (.NOT. ALL(ieee_is_finite(REAL(h, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(h)))) RETURN
332 36560 : IF (.NOT. ALL(ieee_is_finite(REAL(s, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(s)))) RETURN
333 18280 : hnorm = MAX(MAXVAL(ABS(h)), TINY(1.0_dp))
334 18280 : snorm = MAX(MAXVAL(ABS(s)), TINY(1.0_dp))
335 18280 : herror = MAXVAL(ABS(h - CONJG(TRANSPOSE(h))))
336 18280 : serror = MAXVAL(ABS(s - CONJG(TRANSPOSE(s))))
337 80 : IF (herror > tolerance*hnorm .OR. serror > tolerance*snorm) RETURN
338 960 : ALLOCATE (a(n, n), b(n, n), work(MAX(1, 2*n)), eigenvalues(n), rwork(MAX(1, 3*n)))
339 18280 : a(:, :) = s
340 : #if defined(__HAS_IEEE_EXCEPTIONS)
341 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
342 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
343 : #endif
344 80 : CALL zheev('N', 'U', n, a, n, eigenvalues, work, SIZE(work), rwork, info)
345 : #if defined(__HAS_IEEE_EXCEPTIONS)
346 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
347 : #endif
348 80 : status = -4
349 760 : IF (info /= 0 .OR. .NOT. ALL(ieee_is_finite(eigenvalues))) RETURN
350 1600 : bounds = [MINVAL(eigenvalues), MAXVAL(eigenvalues)]
351 80 : status = -2
352 80 : IF (bounds(2) <= 0.0_dp .OR. bounds(1) <= metric_tolerance*bounds(2)) RETURN
353 17766 : a(:, :) = h - energy*s
354 17766 : b(:, :) = s
355 : #if defined(__HAS_IEEE_EXCEPTIONS)
356 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
357 : #endif
358 58 : CALL zhegv(1, 'V', 'U', n, a, n, b, n, eigenvalues, work, SIZE(work), rwork, info)
359 : #if defined(__HAS_IEEE_EXCEPTIONS)
360 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
361 : #endif
362 58 : status = -4
363 654 : IF (info /= 0 .OR. .NOT. ALL(ieee_is_finite(eigenvalues))) RETURN
364 654 : gap = MINVAL(ABS(eigenvalues))
365 : threshold = REAL(n, dp)*(herror + (ABS(energy) + MAXVAL(ABS(eigenvalues)))*serror)/bounds(1) + &
366 : 100.0_dp*EPSILON(1.0_dp)*REAL(n, dp)* &
367 18362 : MAX(MAXVAL(ABS(eigenvalues)), MAXVAL(ABS(h - energy*s))/bounds(1))
368 348 : ALLOCATE (sc(n, n), residual(n, n))
369 153254 : sc(:, :) = MATMUL(s, a)
370 164278 : residual(:, :) = MATMUL(CONJG(TRANSPOSE(a)), sc)
371 654 : DO i = 1, n
372 654 : residual(i, i) = residual(i, i) - 1.0_dp
373 : END DO
374 17766 : errors(1) = SQRT(SUM(ABS(residual)**2))
375 188674 : residual(:, :) = MATMUL(h - energy*s, a) - sc*SPREAD(eigenvalues, 1, n)
376 : errors(2) = SQRT(SUM(ABS(residual)**2))/MAX(TINY(1.0_dp), &
377 53182 : SQRT(SUM(ABS(h - energy*s)**2))*SQRT(SUM(ABS(a)**2)))
378 189266 : flat = MATMUL(sc*SPREAD(SIGN(scale, eigenvalues), 1, n), CONJG(TRANSPOSE(sc)))
379 58 : status = 0
380 80 : END SUBROUTINE flatten_block
381 : END MODULE spectral_localizer_bloch
|
