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 Block-circulant inverse operators for replicated-cell OT preconditioners.
10 : ! **************************************************************************************************
11 : MODULE lattice_preconditioner_operator
12 : USE fft_tools, ONLY: BWFFT,&
13 : FFT_RADIX_NEXT,&
14 : FWFFT,&
15 : fft_1d_many,&
16 : fft_alloc,&
17 : fft_dealloc,&
18 : fft_radix_operations
19 : USE ieee_arithmetic, ONLY: ieee_is_finite
20 : USE kinds, ONLY: dp
21 : USE kpoint_lattice_fft, ONLY: cell_to_k_grid_fft
22 : #include "./base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 : PRIVATE
26 :
27 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lattice_preconditioner_operator'
28 :
29 : INTEGER, PARAMETER, PUBLIC :: lattice_fft_auto_min_cells = 8
30 : REAL(KIND=dp), PARAMETER, PUBLIC :: lattice_fft_auto_max_error = 1.0E-3_dp
31 : REAL(KIND=dp), PARAMETER, PUBLIC :: lattice_fft_auto_max_cost_ratio = 0.8_dp
32 : REAL(KIND=dp), PARAMETER, PUBLIC :: lattice_fft_auto_max_storage_ratio = 4.0_dp
33 :
34 : PUBLIC :: apply_lattice_inverse_dense, &
35 : apply_lattice_state_inverse_dense, &
36 : build_lattice_inverse, &
37 : build_lattice_inverse_from_blocks, &
38 : build_lattice_local_correction, &
39 : lattice_fft_cost_ratio, &
40 : lattice_fft_auto_select, &
41 : lattice_fft_storage_ratio, &
42 : lattice_local_cost_ratio, &
43 : lattice_local_storage_ratio, &
44 : lattice_grid
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief Enumerate a regular direct/reciprocal lattice with the first direction varying fastest.
50 : !> \param lattice_dims number of cells in every lattice direction
51 : !> \param index_to_cell integer direct-lattice coordinates
52 : !> \param xkp fractional reciprocal-lattice coordinates
53 : ! **************************************************************************************************
54 54 : SUBROUTINE lattice_grid(lattice_dims, index_to_cell, xkp)
55 :
56 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
57 : INTEGER, DIMENSION(:, :), INTENT(OUT) :: index_to_cell
58 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: xkp
59 :
60 : INTEGER :: icell, ix, iy, iz
61 :
62 216 : IF (ANY(lattice_dims <= 0)) CPABORT("Invalid lattice dimensions")
63 : IF (SIZE(index_to_cell, 1) < 3 .OR. SIZE(xkp, 1) < 3 .OR. &
64 432 : SIZE(index_to_cell, 2) /= PRODUCT(lattice_dims) .OR. &
65 : SIZE(xkp, 2) /= PRODUCT(lattice_dims)) THEN
66 0 : CPABORT("Invalid lattice-grid output shape")
67 : END IF
68 :
69 54 : icell = 0
70 116 : DO iz = 0, lattice_dims(3) - 1
71 206 : DO iy = 0, lattice_dims(2) - 1
72 424 : DO ix = 0, lattice_dims(1) - 1
73 272 : icell = icell + 1
74 1088 : index_to_cell(:, icell) = [ix, iy, iz]
75 : xkp(:, icell) = [REAL(ix, KIND=dp)/REAL(lattice_dims(1), KIND=dp), &
76 : REAL(iy, KIND=dp)/REAL(lattice_dims(2), KIND=dp), &
77 1178 : REAL(iz, KIND=dp)/REAL(lattice_dims(3), KIND=dp)]
78 : END DO
79 : END DO
80 : END DO
81 :
82 54 : END SUBROUTINE lattice_grid
83 :
84 : ! **************************************************************************************************
85 : !> \brief Project a cell-major SPD operator onto lattice translations and invert its Fourier blocks.
86 : !> \param operator_matrix dense SPD matrix with identical cell blocks stored consecutively
87 : !> \param lattice_dims number of replicated cells in every lattice direction
88 : !> \param inverse_k inverse Hermitian operator block at every reciprocal-grid point
89 : !> \param projection_error relative Frobenius error of the block-circulant projection
90 : !> \param used_fft reports whether construction used the FFT path
91 : !> \param info zero on success; otherwise the failing Cholesky factorization index
92 : ! **************************************************************************************************
93 8 : SUBROUTINE build_lattice_inverse(operator_matrix, lattice_dims, inverse_k, projection_error, &
94 : used_fft, info)
95 :
96 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: operator_matrix
97 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
98 : COMPLEX(KIND=dp), ALLOCATABLE, &
99 : DIMENSION(:, :, :), INTENT(OUT) :: inverse_k
100 : REAL(KIND=dp), INTENT(OUT) :: projection_error
101 : LOGICAL, INTENT(OUT) :: used_fft
102 : INTEGER, INTENT(OUT) :: info
103 :
104 : INTEGER :: n, ncell
105 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: blocks
106 :
107 8 : IF (SIZE(operator_matrix, 1) /= SIZE(operator_matrix, 2)) THEN
108 0 : CPABORT("Lattice preconditioner operator must be square")
109 : END IF
110 32 : IF (ANY(lattice_dims <= 0)) CPABORT("Lattice dimensions must be positive")
111 8 : n = SIZE(operator_matrix, 1)
112 32 : ncell = PRODUCT(lattice_dims)
113 8 : IF (MODULO(n, ncell) /= 0) THEN
114 0 : CPABORT("AO dimension does not match replicated cells")
115 : END IF
116 8 : CALL project_lattice_operator(operator_matrix, lattice_dims, blocks, projection_error)
117 8 : CALL build_lattice_inverse_from_blocks(blocks, lattice_dims, inverse_k, used_fft, info)
118 8 : DEALLOCATE (blocks)
119 :
120 8 : END SUBROUTINE build_lattice_inverse
121 :
122 : ! **************************************************************************************************
123 : !> \brief Transform and invert already projected real-space lattice blocks.
124 : !> \param blocks translation-averaged real-space blocks
125 : !> \param lattice_dims number of replicated cells in every lattice direction
126 : !> \param inverse_k inverse Hermitian operator block at every reciprocal-grid point
127 : !> \param used_fft reports whether construction used the FFT path
128 : !> \param info zero on success; otherwise the failing Cholesky factorization index
129 : ! **************************************************************************************************
130 22 : SUBROUTINE build_lattice_inverse_from_blocks(blocks, lattice_dims, inverse_k, used_fft, info)
131 :
132 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: blocks
133 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
134 : COMPLEX(KIND=dp), ALLOCATABLE, &
135 : DIMENSION(:, :, :), INTENT(OUT) :: inverse_k
136 : LOGICAL, INTENT(OUT) :: used_fft
137 : INTEGER, INTENT(OUT) :: info
138 :
139 : INTEGER :: block_size, ncell
140 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell
141 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp
142 :
143 88 : IF (ANY(lattice_dims <= 0)) CPABORT("Lattice dimensions must be positive")
144 88 : ncell = PRODUCT(lattice_dims)
145 22 : block_size = SIZE(blocks, 1)
146 22 : IF (block_size <= 0 .OR. SIZE(blocks, 2) /= block_size .OR. SIZE(blocks, 3) /= ncell) THEN
147 0 : CPABORT("Projected lattice blocks have inconsistent dimensions")
148 : END IF
149 :
150 110 : ALLOCATE (index_to_cell(3, ncell), xkp(3, ncell))
151 22 : CALL lattice_grid(lattice_dims, index_to_cell, xkp)
152 110 : ALLOCATE (inverse_k(block_size, block_size, ncell))
153 22 : CALL cell_to_k_grid_fft(blocks, index_to_cell, xkp, lattice_dims, inverse_k, used_fft)
154 22 : CALL invert_hermitian_blocks(inverse_k, info)
155 22 : DEALLOCATE (index_to_cell, xkp)
156 :
157 22 : END SUBROUTINE build_lattice_inverse_from_blocks
158 :
159 : ! **************************************************************************************************
160 : !> \brief Build a balanced local correction for selected non-circulant cell blocks.
161 : !> \param operator_matrix dense SPD operator
162 : !> \param lattice_dims number of replicated cells in every lattice direction
163 : !> \param local_cell_count number of cell blocks retained in the correction
164 : !> \param correction_indices selected global AO indices
165 : !> \param operator_columns columns A U for the selected coordinate basis U
166 : !> \param coarse_inverse inverse of U^T A U
167 : !> \param corrected_projection_error residual projection error after removing selected rows/columns
168 : !> \param info zero on success; otherwise the failing Cholesky factorization index
169 : ! **************************************************************************************************
170 8 : SUBROUTINE build_lattice_local_correction(operator_matrix, lattice_dims, local_cell_count, &
171 : correction_indices, operator_columns, coarse_inverse, &
172 : corrected_projection_error, info)
173 :
174 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: operator_matrix
175 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
176 : INTEGER, INTENT(IN) :: local_cell_count
177 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: correction_indices
178 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
179 : INTENT(OUT) :: operator_columns, coarse_inverse
180 : REAL(KIND=dp), INTENT(OUT) :: corrected_projection_error
181 : INTEGER, INTENT(OUT) :: info
182 :
183 : INTEGER :: block_size, cell, chosen, i, index, &
184 : local_rank, n, ncell, row_first, &
185 : row_last
186 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: selected_cells
187 8 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: selected
188 : REAL(KIND=dp) :: denominator
189 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: scores
190 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: projected_operator, residual
191 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: blocks
192 :
193 8 : IF (SIZE(operator_matrix, 1) /= SIZE(operator_matrix, 2)) THEN
194 0 : CPABORT("Lattice local-correction operator must be square")
195 : END IF
196 32 : IF (ANY(lattice_dims <= 0)) CPABORT("Lattice dimensions must be positive")
197 8 : n = SIZE(operator_matrix, 1)
198 32 : ncell = PRODUCT(lattice_dims)
199 8 : IF (MODULO(n, ncell) /= 0) CPABORT("AO dimension does not match replicated cells")
200 8 : IF (local_cell_count <= 0 .OR. local_cell_count >= ncell) THEN
201 0 : CPABORT("Lattice local correction requires between one and ncell-1 cells")
202 : END IF
203 8 : block_size = n/ncell
204 8 : local_rank = local_cell_count*block_size
205 :
206 : CALL project_lattice_operator(operator_matrix, lattice_dims, blocks, &
207 8 : corrected_projection_error, projected_operator)
208 8 : DEALLOCATE (blocks)
209 32 : ALLOCATE (residual(n, n))
210 776 : residual(:, :) = operator_matrix - projected_operator
211 8 : DEALLOCATE (projected_operator)
212 56 : ALLOCATE (scores(ncell), selected(ncell), selected_cells(local_cell_count))
213 8 : selected = .FALSE.
214 16 : DO chosen = 1, local_cell_count
215 40 : scores = -1.0_dp
216 40 : DO cell = 1, ncell
217 32 : IF (selected(cell)) CYCLE
218 32 : row_first = (cell - 1)*block_size + 1
219 32 : row_last = cell*block_size
220 : scores(cell) = SUM(residual(row_first:row_last, :)**2) + &
221 : SUM(residual(:, row_first:row_last)**2) - &
222 2056 : SUM(residual(row_first:row_last, row_first:row_last)**2)
223 : END DO
224 48 : index = MAXLOC(scores, DIM=1, MASK=.NOT. selected)
225 8 : selected(index) = .TRUE.
226 8 : selected_cells(chosen) = index
227 8 : row_first = (index - 1)*block_size + 1
228 8 : row_last = index*block_size
229 200 : residual(row_first:row_last, :) = 0.0_dp
230 160 : residual(:, row_first:row_last) = 0.0_dp
231 : END DO
232 :
233 24 : ALLOCATE (correction_indices(local_rank))
234 8 : index = 0
235 16 : DO chosen = 1, local_cell_count
236 8 : cell = selected_cells(chosen)
237 32 : DO i = 1, block_size
238 16 : index = index + 1
239 24 : correction_indices(index) = (cell - 1)*block_size + i
240 : END DO
241 : END DO
242 :
243 776 : denominator = SUM(operator_matrix**2)
244 776 : corrected_projection_error = SQRT(SUM(residual**2)/MAX(denominator, TINY(denominator)))
245 :
246 56 : ALLOCATE (operator_columns(n, local_rank), coarse_inverse(local_rank, local_rank))
247 152 : operator_columns(:, :) = operator_matrix(:, correction_indices)
248 56 : coarse_inverse(:, :) = operator_matrix(correction_indices, correction_indices)
249 8 : CALL dpotrf('U', local_rank, coarse_inverse, local_rank, info)
250 8 : IF (info == 0) CALL dpotri('U', local_rank, coarse_inverse, local_rank, info)
251 8 : IF (info == 0) THEN
252 16 : DO i = 1, local_rank - 1
253 24 : coarse_inverse(i + 1:local_rank, i) = coarse_inverse(i, i + 1:local_rank)
254 : END DO
255 : END IF
256 :
257 8 : DEALLOCATE (residual, scores, selected, selected_cells)
258 :
259 8 : END SUBROUTINE build_lattice_local_correction
260 :
261 : ! **************************************************************************************************
262 : !> \brief Project a dense cell-major operator onto its block-circulant translation average.
263 : !> \param operator_matrix dense square operator
264 : !> \param lattice_dims number of replicated cells in every lattice direction
265 : !> \param blocks translation-averaged real-space blocks
266 : !> \param projection_error relative Frobenius projection error
267 : !> \param projected_operator optional dense block-circulant projection
268 : ! **************************************************************************************************
269 16 : SUBROUTINE project_lattice_operator(operator_matrix, lattice_dims, blocks, projection_error, &
270 : projected_operator)
271 :
272 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: operator_matrix
273 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
274 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
275 : INTENT(OUT) :: blocks
276 : REAL(KIND=dp), INTENT(OUT) :: projection_error
277 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
278 : INTENT(OUT), OPTIONAL :: projected_operator
279 :
280 : INTEGER :: alpha, beta, block_size, delta_index, &
281 : icell, jcell, n, ncell
282 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell
283 : INTEGER, DIMENSION(3) :: delta
284 : REAL(KIND=dp) :: denominator, difference, numerator
285 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp
286 :
287 16 : IF (SIZE(operator_matrix, 1) /= SIZE(operator_matrix, 2)) THEN
288 0 : CPABORT("Lattice preconditioner operator must be square")
289 : END IF
290 64 : IF (ANY(lattice_dims <= 0)) CPABORT("Lattice dimensions must be positive")
291 16 : n = SIZE(operator_matrix, 1)
292 64 : ncell = PRODUCT(lattice_dims)
293 16 : IF (MODULO(n, ncell) /= 0) CPABORT("AO dimension does not match replicated cells")
294 16 : block_size = n/ncell
295 :
296 80 : ALLOCATE (index_to_cell(3, ncell), xkp(3, ncell))
297 16 : CALL lattice_grid(lattice_dims, index_to_cell, xkp)
298 80 : ALLOCATE (blocks(block_size, block_size, ncell), source=0.0_dp)
299 88 : DO jcell = 1, ncell
300 536 : DO icell = 1, ncell
301 : ! With X(k)=sum_R exp(+ikR) X(R), the convolution block is indexed by R-T.
302 1792 : delta = MODULO(index_to_cell(:, icell) - index_to_cell(:, jcell), lattice_dims)
303 448 : delta_index = cell_linear_index(delta, lattice_dims)
304 1416 : DO beta = 1, block_size
305 3136 : DO alpha = 1, block_size
306 : blocks(alpha, beta, delta_index) = blocks(alpha, beta, delta_index) + &
307 : operator_matrix((icell - 1)*block_size + alpha, &
308 : (jcell - 1)*block_size + beta)/ &
309 2688 : REAL(ncell, KIND=dp)
310 : END DO
311 : END DO
312 : END DO
313 : END DO
314 :
315 40 : IF (PRESENT(projected_operator)) ALLOCATE (projected_operator(n, n))
316 : numerator = 0.0_dp
317 : denominator = 0.0_dp
318 88 : DO jcell = 1, ncell
319 536 : DO icell = 1, ncell
320 1792 : delta = MODULO(index_to_cell(:, icell) - index_to_cell(:, jcell), lattice_dims)
321 448 : delta_index = cell_linear_index(delta, lattice_dims)
322 1416 : DO beta = 1, block_size
323 3136 : DO alpha = 1, block_size
324 1792 : IF (PRESENT(projected_operator)) THEN
325 : projected_operator((icell - 1)*block_size + alpha, &
326 : (jcell - 1)*block_size + beta) = &
327 704 : blocks(alpha, beta, delta_index)
328 : END IF
329 : difference = operator_matrix((icell - 1)*block_size + alpha, &
330 : (jcell - 1)*block_size + beta) - &
331 1792 : blocks(alpha, beta, delta_index)
332 1792 : numerator = numerator + difference*difference
333 : denominator = denominator + &
334 : operator_matrix((icell - 1)*block_size + alpha, &
335 2688 : (jcell - 1)*block_size + beta)**2
336 : END DO
337 : END DO
338 : END DO
339 : END DO
340 16 : projection_error = SQRT(numerator/MAX(denominator, TINY(denominator)))
341 :
342 16 : DEALLOCATE (index_to_cell, xkp)
343 :
344 16 : END SUBROUTINE project_lattice_operator
345 :
346 : ! **************************************************************************************************
347 : !> \brief Estimate lattice-FFT application work relative to a dense inverse application.
348 : !> \param lattice_dims number of explicitly replicated cells
349 : !> \param block_size number of orbitals in the reference cell
350 : !> \param rhs_count number of orbital columns transformed together
351 : !> \param num_pe number of message-passing ranks sharing the operator
352 : !> \return estimated lattice-to-dense work ratio
353 : ! **************************************************************************************************
354 30 : PURE REAL(KIND=dp) FUNCTION lattice_fft_cost_ratio(lattice_dims, block_size, rhs_count, num_pe) &
355 : RESULT(ratio)
356 :
357 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
358 : INTEGER, INTENT(IN) :: block_size, rhs_count, num_pe
359 :
360 : INTEGER :: ncell, nrhs
361 : REAL(KIND=dp) :: block_size_real, communication_work, &
362 : dense_work, fft_work, matrix_work, &
363 : ncell_real, nrhs_real, setup_work
364 :
365 120 : ncell = MAX(1, PRODUCT(lattice_dims))
366 30 : nrhs = MAX(1, rhs_count)
367 30 : IF (block_size <= 0) THEN
368 30 : ratio = HUGE(ratio)
369 : RETURN
370 : END IF
371 30 : ncell_real = REAL(ncell, KIND=dp)
372 30 : nrhs_real = REAL(nrhs, KIND=dp)
373 30 : block_size_real = REAL(block_size, KIND=dp)
374 :
375 : ! This deliberately simple model counts real dense multiply-adds as the reference.
376 : ! Four scalar operations per FFT butterfly account for the two complex transforms;
377 : ! the fixed packing term prevents AUTO from selecting tiny batches. The logarithmic
378 : ! communication term represents the input/output all-to-all routing.
379 30 : dense_work = (ncell_real*block_size_real)**2*nrhs_real
380 30 : matrix_work = ncell_real*block_size_real**2*nrhs_real
381 : fft_work = 4.0_dp*ncell_real*block_size_real*nrhs_real* &
382 30 : LOG(REAL(MAX(2, ncell), KIND=dp))/LOG(2.0_dp)
383 30 : setup_work = 8.0_dp*ncell_real*block_size_real
384 : communication_work = 2.0_dp*ncell_real*block_size_real*nrhs_real* &
385 30 : LOG(REAL(MAX(2, num_pe + 1), KIND=dp))/LOG(2.0_dp)
386 30 : ratio = (matrix_work + fft_work + setup_work + communication_work)/MAX(dense_work, TINY(dense_work))
387 :
388 30 : END FUNCTION lattice_fft_cost_ratio
389 :
390 : ! **************************************************************************************************
391 : !> \brief Estimate reciprocal-block storage relative to one dense real operator.
392 : !> \param lattice_dims number of explicitly replicated cells
393 : !> \param operator_count number of distinct state-dependent reciprocal operators
394 : !> \return estimated lattice-to-dense storage ratio
395 : ! **************************************************************************************************
396 26 : PURE REAL(KIND=dp) FUNCTION lattice_fft_storage_ratio(lattice_dims, operator_count) RESULT(ratio)
397 :
398 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
399 : INTEGER, INTENT(IN) :: operator_count
400 :
401 : INTEGER :: ncell
402 :
403 104 : ncell = MAX(1, PRODUCT(lattice_dims))
404 : ! Reciprocal blocks are complex; the block-size dependence cancels against
405 : ! the reference dense matrix. A state-independent operator has count one.
406 26 : ratio = 2.0_dp*REAL(MAX(1, operator_count), KIND=dp)/REAL(ncell, KIND=dp)
407 :
408 26 : END FUNCTION lattice_fft_storage_ratio
409 :
410 : ! **************************************************************************************************
411 : !> \brief Estimate balanced-local-correction work relative to a dense inverse application.
412 : !> \param lattice_dims number of replicated cells
413 : !> \param local_cell_count number of cell blocks retained in the correction
414 : !> \return estimated correction-to-dense work ratio
415 : ! **************************************************************************************************
416 20 : PURE REAL(KIND=dp) FUNCTION lattice_local_cost_ratio(lattice_dims, local_cell_count) &
417 : RESULT(ratio)
418 :
419 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
420 : INTEGER, INTENT(IN) :: local_cell_count
421 :
422 : INTEGER :: ncell
423 : REAL(KIND=dp) :: local_fraction
424 :
425 80 : ncell = MAX(1, PRODUCT(lattice_dims))
426 20 : IF (local_cell_count <= 0) THEN
427 20 : ratio = 0.0_dp
428 : RETURN
429 : END IF
430 6 : local_fraction = REAL(local_cell_count, KIND=dp)/REAL(ncell, KIND=dp)
431 : ! Both corrected and dense applications scale linearly with the right-hand-side count.
432 6 : ratio = 2.0_dp*(local_fraction + local_fraction**2)
433 :
434 6 : END FUNCTION lattice_local_cost_ratio
435 :
436 : ! **************************************************************************************************
437 : !> \brief Estimate balanced-local-correction storage relative to one dense real operator.
438 : !> \param lattice_dims number of replicated cells
439 : !> \param operator_count number of distinct state-dependent local corrections
440 : !> \param local_cell_count number of cell blocks retained in each correction
441 : !> \return estimated correction-to-dense storage ratio
442 : ! **************************************************************************************************
443 20 : PURE REAL(KIND=dp) FUNCTION lattice_local_storage_ratio(lattice_dims, operator_count, &
444 : local_cell_count) RESULT(ratio)
445 :
446 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
447 : INTEGER, INTENT(IN) :: operator_count, local_cell_count
448 :
449 : INTEGER :: ncell
450 : REAL(KIND=dp) :: local_fraction
451 :
452 80 : ncell = MAX(1, PRODUCT(lattice_dims))
453 20 : IF (local_cell_count <= 0) THEN
454 20 : ratio = 0.0_dp
455 : RETURN
456 : END IF
457 6 : local_fraction = REAL(local_cell_count, KIND=dp)/REAL(ncell, KIND=dp)
458 6 : ratio = REAL(MAX(1, operator_count), KIND=dp)*(local_fraction + local_fraction**2)
459 :
460 6 : END FUNCTION lattice_local_storage_ratio
461 :
462 : ! **************************************************************************************************
463 : !> \brief Conservative error-, work-, and storage-aware AUTO selector.
464 : !> \param lattice_dims number of explicitly replicated cells
465 : !> \param block_size number of orbitals in the reference cell
466 : !> \param rhs_count number of orbital columns transformed together
467 : !> \param operator_count number of distinct reciprocal operators stored
468 : !> \param num_pe number of message-passing ranks sharing the operator
469 : !> \param projection_error relative Frobenius projection error
470 : !> \param used_fft whether construction used the FFT path
471 : !> \param info zero when all projected Fourier blocks were positive definite
472 : !> \return true when AUTO may use the lattice inverse
473 : ! **************************************************************************************************
474 16 : PURE LOGICAL FUNCTION lattice_fft_auto_select(lattice_dims, block_size, rhs_count, operator_count, &
475 : num_pe, projection_error, used_fft, info) &
476 : RESULT(selected)
477 :
478 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
479 : INTEGER, INTENT(IN) :: block_size, rhs_count, operator_count, &
480 : num_pe
481 : REAL(KIND=dp), INTENT(IN) :: projection_error
482 : LOGICAL, INTENT(IN) :: used_fft
483 : INTEGER, INTENT(IN) :: info
484 :
485 : selected = PRODUCT(lattice_dims) >= lattice_fft_auto_min_cells .AND. &
486 68 : projection_error <= lattice_fft_auto_max_error .AND. used_fft .AND. info == 0
487 : selected = selected .AND. &
488 : lattice_fft_cost_ratio(lattice_dims, block_size, rhs_count, num_pe) <= &
489 : lattice_fft_auto_max_cost_ratio .AND. &
490 : lattice_fft_storage_ratio(lattice_dims, operator_count) <= &
491 8 : lattice_fft_auto_max_storage_ratio
492 :
493 16 : END FUNCTION lattice_fft_auto_select
494 :
495 : ! **************************************************************************************************
496 : !> \brief Apply inverse Fourier blocks to a dense cell-major real matrix.
497 : !> \param inverse_k inverse overlap blocks on a complete reciprocal grid
498 : !> \param lattice_dims number of replicated cells in every lattice direction
499 : !> \param matrix_in matrix to precondition
500 : !> \param matrix_out preconditioned matrix
501 : !> \param used_fft reports whether both transforms used the FFT path
502 : !> \param imaginary_residual largest discarded imaginary component
503 : ! **************************************************************************************************
504 298 : SUBROUTINE apply_lattice_inverse_dense(inverse_k, lattice_dims, matrix_in, matrix_out, &
505 : used_fft, imaginary_residual)
506 :
507 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: inverse_k
508 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
509 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix_in
510 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: matrix_out
511 : LOGICAL, INTENT(OUT) :: used_fft
512 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: imaginary_residual
513 :
514 298 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: input_k, input_rs, output_k, output_rs
515 : INTEGER :: block_size, icell, ncell
516 : LOGICAL :: forward_fft, inverse_fft
517 :
518 1192 : IF (ANY(lattice_dims <= 0)) CPABORT("Lattice dimensions must be positive")
519 1192 : ncell = PRODUCT(lattice_dims)
520 298 : block_size = SIZE(inverse_k, 1)
521 298 : IF (SIZE(inverse_k, 2) /= block_size .OR. SIZE(inverse_k, 3) /= ncell) THEN
522 0 : CPABORT("Inconsistent inverse blocks in lattice preconditioner")
523 : END IF
524 894 : IF (SIZE(matrix_in, 1) /= block_size*ncell .OR. &
525 : ANY(SHAPE(matrix_out) /= SHAPE(matrix_in))) THEN
526 0 : CPABORT("Inconsistent dense matrix shape in lattice preconditioner")
527 : END IF
528 :
529 1490 : ALLOCATE (input_rs(block_size, SIZE(matrix_in, 2), ncell))
530 1506 : DO icell = 1, ncell
531 : input_rs(:, :, icell) = CMPLX( &
532 43818 : matrix_in((icell - 1)*block_size + 1:icell*block_size, :), 0.0_dp, KIND=dp)
533 : END DO
534 :
535 1192 : ALLOCATE (input_k(block_size, SIZE(matrix_in, 2), ncell))
536 1192 : ALLOCATE (output_k(block_size, SIZE(matrix_in, 2), ncell))
537 1192 : ALLOCATE (output_rs(block_size, SIZE(matrix_in, 2), ncell))
538 298 : CALL batched_lattice_fft(input_rs, lattice_dims, BWFFT, input_k, forward_fft)
539 1506 : DO icell = 1, ncell
540 625642 : output_k(:, :, icell) = MATMUL(inverse_k(:, :, icell), input_k(:, :, icell))
541 : END DO
542 298 : CALL batched_lattice_fft(output_k, lattice_dims, FWFFT, output_rs, inverse_fft)
543 298 : used_fft = forward_fft .AND. inverse_fft
544 43516 : IF (PRESENT(imaginary_residual)) imaginary_residual = MAXVAL(ABS(AIMAG(output_rs)))
545 :
546 1506 : DO icell = 1, ncell
547 43818 : matrix_out((icell - 1)*block_size + 1:icell*block_size, :) = REAL(output_rs(:, :, icell), KIND=dp)
548 : END DO
549 :
550 298 : DEALLOCATE (input_k, input_rs, output_k, output_rs)
551 :
552 298 : END SUBROUTINE apply_lattice_inverse_dense
553 :
554 : ! **************************************************************************************************
555 : !> \brief Apply one inverse Fourier operator to every column of a dense real matrix.
556 : !> \param inverse_k state-dependent inverse blocks on a complete reciprocal grid
557 : !> \param lattice_dims number of replicated cells in every lattice direction
558 : !> \param matrix_in matrix to precondition, with one state per column
559 : !> \param matrix_out preconditioned matrix
560 : !> \param used_fft reports whether both transforms used the FFT path
561 : !> \param imaginary_residual largest discarded imaginary component
562 : !> \param state_indices optional global operator index for every input column
563 : ! **************************************************************************************************
564 28 : SUBROUTINE apply_lattice_state_inverse_dense(inverse_k, lattice_dims, matrix_in, matrix_out, &
565 14 : used_fft, imaginary_residual, state_indices)
566 :
567 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :), &
568 : INTENT(IN) :: inverse_k
569 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
570 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix_in
571 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: matrix_out
572 : LOGICAL, INTENT(OUT) :: used_fft
573 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: imaginary_residual
574 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: state_indices
575 :
576 14 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: input_k, input_rs, output_k, output_rs
577 : INTEGER :: block_size, icell, j, ncell
578 : LOGICAL :: forward_fft, inverse_fft
579 :
580 56 : IF (ANY(lattice_dims <= 0)) CPABORT("Lattice dimensions must be positive")
581 56 : ncell = PRODUCT(lattice_dims)
582 14 : block_size = SIZE(inverse_k, 1)
583 14 : IF (SIZE(inverse_k, 2) /= block_size .OR. SIZE(inverse_k, 3) /= ncell) THEN
584 0 : CPABORT("Inconsistent state-dependent inverse blocks in lattice preconditioner")
585 : END IF
586 14 : IF (PRESENT(state_indices)) THEN
587 : IF (SIZE(state_indices) /= SIZE(matrix_in, 2) .OR. &
588 36 : ANY(state_indices < 1) .OR. ANY(state_indices > SIZE(inverse_k, 4))) THEN
589 0 : CPABORT("Invalid state indices in lattice preconditioner")
590 : END IF
591 2 : ELSE IF (SIZE(inverse_k, 4) /= SIZE(matrix_in, 2)) THEN
592 0 : CPABORT("Inconsistent state-dependent inverse blocks in lattice preconditioner")
593 : END IF
594 42 : IF (SIZE(matrix_in, 1) /= block_size*ncell .OR. &
595 : ANY(SHAPE(matrix_out) /= SHAPE(matrix_in))) THEN
596 0 : CPABORT("Inconsistent dense state matrix shape in lattice preconditioner")
597 : END IF
598 :
599 70 : ALLOCATE (input_rs(block_size, SIZE(matrix_in, 2), ncell))
600 54 : DO icell = 1, ncell
601 : input_rs(:, :, icell) = CMPLX( &
602 270 : matrix_in((icell - 1)*block_size + 1:icell*block_size, :), 0.0_dp, KIND=dp)
603 : END DO
604 :
605 56 : ALLOCATE (input_k(block_size, SIZE(matrix_in, 2), ncell))
606 56 : ALLOCATE (output_k(block_size, SIZE(matrix_in, 2), ncell))
607 56 : ALLOCATE (output_rs(block_size, SIZE(matrix_in, 2), ncell))
608 14 : CALL batched_lattice_fft(input_rs, lattice_dims, BWFFT, input_k, forward_fft)
609 54 : DO icell = 1, ncell
610 126 : DO j = 1, SIZE(matrix_in, 2)
611 112 : IF (PRESENT(state_indices)) THEN
612 : output_k(:, j, icell) = &
613 216 : MATMUL(inverse_k(:, :, icell, state_indices(j)), input_k(:, j, icell))
614 : ELSE
615 432 : output_k(:, j, icell) = MATMUL(inverse_k(:, :, icell, j), input_k(:, j, icell))
616 : END IF
617 : END DO
618 : END DO
619 14 : CALL batched_lattice_fft(output_k, lattice_dims, FWFFT, output_rs, inverse_fft)
620 14 : used_fft = forward_fft .AND. inverse_fft
621 284 : IF (PRESENT(imaginary_residual)) imaginary_residual = MAXVAL(ABS(AIMAG(output_rs)))
622 :
623 54 : DO icell = 1, ncell
624 270 : matrix_out((icell - 1)*block_size + 1:icell*block_size, :) = REAL(output_rs(:, :, icell), KIND=dp)
625 : END DO
626 :
627 14 : DEALLOCATE (input_k, input_rs, output_k, output_rs)
628 :
629 14 : END SUBROUTINE apply_lattice_state_inverse_dense
630 :
631 : ! **************************************************************************************************
632 : !> \brief Transform a complete replicated-cell batch with three batched one-dimensional FFTs.
633 : !> \param values_in matrices on the complete direct or reciprocal lattice
634 : !> \param lattice_dims logical lattice dimensions
635 : !> \param fsign transform direction, BWFFT for cell-to-k and FWFFT for k-to-cell
636 : !> \param values_out transformed matrices on the complete lattice
637 : !> \param used_fft reports whether the batched FFT path succeeded
638 : ! **************************************************************************************************
639 624 : SUBROUTINE batched_lattice_fft(values_in, lattice_dims, fsign, values_out, used_fft)
640 :
641 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: values_in
642 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
643 : INTEGER, INTENT(IN) :: fsign
644 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: values_out
645 : LOGICAL, INTENT(OUT) :: used_fft
646 :
647 : COMPLEX(KIND=dp), ALLOCATABLE, &
648 624 : DIMENSION(:, :, :, :) :: grid
649 : INTEGER :: attempt, batch_index, d, i, icell, ix, &
650 : iy, iz, j, radix_length, stat
651 : INTEGER, DIMENSION(3) :: input_stride, nfft, output_stride, stride
652 : LOGICAL :: compatible
653 : REAL(KIND=dp) :: normalization
654 :
655 6240 : IF (ANY(lattice_dims <= 0) .OR. SIZE(values_in, 3) /= PRODUCT(lattice_dims) .OR. &
656 : ANY(SHAPE(values_out) /= SHAPE(values_in))) THEN
657 0 : CPABORT("Invalid shape in batched lattice FFT")
658 : END IF
659 624 : IF (fsign /= FWFFT .AND. fsign /= BWFFT) THEN
660 0 : CPABORT("Invalid direction in batched lattice FFT")
661 : END IF
662 624 : nfft = lattice_dims
663 624 : compatible = .TRUE.
664 2496 : DO d = 1, 3
665 5760 : compatible = .FALSE.
666 5760 : DO attempt = 0, 15
667 5760 : IF (nfft(d) >= 3) THEN
668 2952 : CALL fft_radix_operations(nfft(d), radix_length, FFT_RADIX_NEXT)
669 2952 : IF (radix_length == nfft(d)) THEN
670 : compatible = .TRUE.
671 : EXIT
672 : END IF
673 : END IF
674 3888 : nfft(d) = nfft(d) + lattice_dims(d)
675 : END DO
676 2496 : IF (.NOT. compatible) EXIT
677 : END DO
678 624 : IF (.NOT. compatible) THEN
679 0 : values_out = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
680 0 : used_fft = .FALSE.
681 0 : RETURN
682 : END IF
683 :
684 2496 : stride = nfft/lattice_dims
685 2496 : input_stride = 1
686 2496 : output_stride = 1
687 624 : IF (fsign == BWFFT) input_stride = stride
688 624 : IF (fsign == FWFFT) output_stride = stride
689 3744 : ALLOCATE (grid(SIZE(values_in, 1)*SIZE(values_in, 2), nfft(1), nfft(2), nfft(3)))
690 624 : grid = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
691 624 : icell = 0
692 1320 : DO iz = 0, lattice_dims(3) - 1
693 2184 : DO iy = 0, lattice_dims(2) - 1
694 4056 : DO ix = 0, lattice_dims(1) - 1
695 2496 : icell = icell + 1
696 17536 : DO j = 1, SIZE(values_in, 2)
697 87552 : DO i = 1, SIZE(values_in, 1)
698 70880 : batch_index = i + SIZE(values_in, 1)*(j - 1)
699 : grid(batch_index, ix*input_stride(1) + 1, iy*input_stride(2) + 1, &
700 85056 : iz*input_stride(3) + 1) = values_in(i, j, icell)
701 : END DO
702 : END DO
703 : END DO
704 : END DO
705 : END DO
706 :
707 624 : stat = 0
708 2496 : DO d = 1, 3
709 1872 : CALL transform_grid_dimension(grid, d, fsign, stat)
710 2496 : IF (stat /= 0) EXIT
711 : END DO
712 624 : IF (stat /= 0) THEN
713 0 : values_out = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
714 0 : used_fft = .FALSE.
715 0 : DEALLOCATE (grid)
716 0 : RETURN
717 : END IF
718 :
719 624 : normalization = 1.0_dp
720 624 : IF (fsign == FWFFT) THEN
721 2184 : normalization = REAL(PRODUCT(nfft), KIND=dp)/REAL(PRODUCT(lattice_dims), KIND=dp)
722 : END IF
723 624 : icell = 0
724 1320 : DO iz = 0, lattice_dims(3) - 1
725 2184 : DO iy = 0, lattice_dims(2) - 1
726 4056 : DO ix = 0, lattice_dims(1) - 1
727 2496 : icell = icell + 1
728 17536 : DO j = 1, SIZE(values_out, 2)
729 87552 : DO i = 1, SIZE(values_out, 1)
730 70880 : batch_index = i + SIZE(values_out, 1)*(j - 1)
731 : values_out(i, j, icell) = normalization* &
732 : grid(batch_index, ix*output_stride(1) + 1, iy*output_stride(2) + 1, &
733 85056 : iz*output_stride(3) + 1)
734 : END DO
735 : END DO
736 : END DO
737 : END DO
738 : END DO
739 624 : used_fft = .TRUE.
740 624 : DEALLOCATE (grid)
741 :
742 0 : END SUBROUTINE batched_lattice_fft
743 :
744 : ! **************************************************************************************************
745 : !> \brief Transform every line of one dimension of a four-dimensional FFT batch.
746 : !> \param grid batched three-dimensional grid
747 : !> \param direction grid dimension to transform
748 : !> \param fsign transform direction
749 : !> \param stat zero on success
750 : ! **************************************************************************************************
751 1872 : SUBROUTINE transform_grid_dimension(grid, direction, fsign, stat)
752 :
753 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :), &
754 : INTENT(INOUT) :: grid
755 : INTEGER, INTENT(IN) :: direction, fsign
756 : INTEGER, INTENT(OUT) :: stat
757 :
758 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
759 1872 : POINTER :: line_in, line_out
760 : INTEGER :: ibatch, iline, ix, iy, iz, m, n
761 : REAL(KIND=dp) :: scale
762 :
763 1872 : n = SIZE(grid, direction + 1)
764 9360 : m = SIZE(grid)/n
765 5616 : CALL fft_alloc(line_in, [n, m])
766 5616 : CALL fft_alloc(line_out, [n, m])
767 1872 : iline = 0
768 624 : SELECT CASE (direction)
769 : CASE (1)
770 3120 : DO iz = 1, SIZE(grid, 4)
771 13104 : DO iy = 1, SIZE(grid, 3)
772 440896 : DO ibatch = 1, SIZE(grid, 1)
773 428416 : iline = iline + 1
774 2152064 : line_in(:, iline) = grid(ibatch, :, iy, iz)
775 : END DO
776 : END DO
777 : END DO
778 : CASE (2)
779 3120 : DO iz = 1, SIZE(grid, 4)
780 13104 : DO ix = 1, SIZE(grid, 2)
781 440896 : DO ibatch = 1, SIZE(grid, 1)
782 428416 : iline = iline + 1
783 2152064 : line_in(:, iline) = grid(ibatch, ix, :, iz)
784 : END DO
785 : END DO
786 : END DO
787 : CASE (3)
788 3120 : DO iy = 1, SIZE(grid, 3)
789 13104 : DO ix = 1, SIZE(grid, 2)
790 440896 : DO ibatch = 1, SIZE(grid, 1)
791 428416 : iline = iline + 1
792 2152064 : line_in(:, iline) = grid(ibatch, ix, iy, :)
793 : END DO
794 : END DO
795 : END DO
796 : CASE DEFAULT
797 1872 : CPABORT("Invalid batched FFT grid direction")
798 : END SELECT
799 :
800 1872 : scale = 1.0_dp
801 1872 : IF (fsign == FWFFT) scale = 1.0_dp/REAL(n, KIND=dp)
802 : CALL fft_1d_many(fsign, n, m, .FALSE., .FALSE., n, n, &
803 1872 : line_in, line_out, scale, stat)
804 1872 : IF (stat == 0) THEN
805 1872 : iline = 0
806 624 : SELECT CASE (direction)
807 : CASE (1)
808 3120 : DO iz = 1, SIZE(grid, 4)
809 13104 : DO iy = 1, SIZE(grid, 3)
810 440896 : DO ibatch = 1, SIZE(grid, 1)
811 428416 : iline = iline + 1
812 2152064 : grid(ibatch, :, iy, iz) = line_out(:, iline)
813 : END DO
814 : END DO
815 : END DO
816 : CASE (2)
817 3120 : DO iz = 1, SIZE(grid, 4)
818 13104 : DO ix = 1, SIZE(grid, 2)
819 440896 : DO ibatch = 1, SIZE(grid, 1)
820 428416 : iline = iline + 1
821 2152064 : grid(ibatch, ix, :, iz) = line_out(:, iline)
822 : END DO
823 : END DO
824 : END DO
825 : CASE (3)
826 4368 : DO iy = 1, SIZE(grid, 3)
827 13104 : DO ix = 1, SIZE(grid, 2)
828 440896 : DO ibatch = 1, SIZE(grid, 1)
829 428416 : iline = iline + 1
830 2152064 : grid(ibatch, ix, iy, :) = line_out(:, iline)
831 : END DO
832 : END DO
833 : END DO
834 : END SELECT
835 : END IF
836 1872 : CALL fft_dealloc(line_in)
837 1872 : CALL fft_dealloc(line_out)
838 :
839 1872 : END SUBROUTINE transform_grid_dimension
840 :
841 : ! **************************************************************************************************
842 : !> \brief Convert zero-based lattice coordinates to the cell-major linear index.
843 : !> \param coordinate zero-based lattice coordinate
844 : !> \param lattice_dims number of cells in every direction
845 : !> \return one-based linear index
846 : ! **************************************************************************************************
847 896 : PURE INTEGER FUNCTION cell_linear_index(coordinate, lattice_dims) RESULT(index)
848 :
849 : INTEGER, DIMENSION(3), INTENT(IN) :: coordinate, lattice_dims
850 :
851 : index = coordinate(1) + lattice_dims(1)*(coordinate(2) + &
852 896 : lattice_dims(2)*coordinate(3)) + 1
853 :
854 896 : END FUNCTION cell_linear_index
855 :
856 : ! **************************************************************************************************
857 : !> \brief Hermitize, Cholesky factorize, and invert each small complex block.
858 : !> \param blocks matrix blocks, replaced by their inverse
859 : !> \param info zero on success; otherwise the one-based failing block index
860 : ! **************************************************************************************************
861 22 : SUBROUTINE invert_hermitian_blocks(blocks, info)
862 :
863 : COMPLEX(KIND=dp), DIMENSION(:, :, :), &
864 : INTENT(INOUT) :: blocks
865 : INTEGER, INTENT(OUT) :: info
866 :
867 : COMPLEX(KIND=dp) :: value
868 22 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: solution, work
869 22 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: factor, inverse
870 : INTEGER :: i, iblock, j, k, m, n
871 : REAL(KIND=dp) :: diagonal
872 :
873 22 : n = SIZE(blocks, 1)
874 22 : CPASSERT(SIZE(blocks, 2) == n)
875 198 : ALLOCATE (factor(n, n), inverse(n, n), solution(n), work(n))
876 22 : info = 0
877 134 : DO iblock = 1, SIZE(blocks, 3)
878 112 : factor = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
879 : inverse(:, :) = 0.5_dp*(blocks(:, :, iblock) + &
880 5736 : CONJG(TRANSPOSE(blocks(:, :, iblock))))
881 :
882 : ! Small-block upper Cholesky factorization, inverse=U^H U.
883 516 : DO j = 1, n
884 404 : diagonal = REAL(inverse(j, j), KIND=dp)
885 2812 : DO m = 1, j - 1
886 2812 : diagonal = diagonal - ABS(factor(m, j))**2
887 : END DO
888 404 : IF (.NOT. ieee_is_finite(diagonal) .OR. diagonal <= 0.0_dp) THEN
889 0 : info = iblock
890 0 : DEALLOCATE (factor, inverse, solution, work)
891 : RETURN
892 : END IF
893 404 : factor(j, j) = CMPLX(SQRT(diagonal), 0.0_dp, KIND=dp)
894 2924 : DO k = j + 1, n
895 2408 : value = inverse(j, k)
896 24520 : DO m = 1, j - 1
897 24520 : value = value - CONJG(factor(m, j))*factor(m, k)
898 : END DO
899 2812 : factor(j, k) = value/factor(j, j)
900 : END DO
901 : END DO
902 :
903 : ! Solve U^H U x_j=e_j for every inverse column. A small explicit
904 : ! triangular solve keeps the complete small-block kernel self-contained.
905 516 : DO j = 1, n
906 404 : work = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
907 404 : solution = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
908 5624 : DO i = 1, n
909 5220 : value = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
910 5220 : IF (i == j) value = CMPLX(1.0_dp, 0.0_dp, KIND=dp)
911 76372 : DO m = 1, i - 1
912 76372 : value = value - CONJG(factor(m, i))*work(m)
913 : END DO
914 5624 : work(i) = value/CONJG(factor(i, i))
915 : END DO
916 5624 : DO i = n, 1, -1
917 5220 : value = work(i)
918 76372 : DO m = i + 1, n
919 76372 : value = value - factor(i, m)*solution(m)
920 : END DO
921 5624 : solution(i) = value/factor(i, i)
922 : END DO
923 5736 : inverse(:, j) = solution
924 : END DO
925 5758 : blocks(:, :, iblock) = 0.5_dp*(inverse + CONJG(TRANSPOSE(inverse)))
926 : END DO
927 22 : DEALLOCATE (factor, inverse, solution, work)
928 :
929 : END SUBROUTINE invert_hermitian_blocks
930 :
931 : END MODULE lattice_preconditioner_operator
|
