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 2 : PROGRAM kpoint_lattice_fft_unittest
8 2 : USE fft_tools, ONLY: finalize_fft,&
9 : init_fft
10 : USE kinds, ONLY: dp
11 : USE kpoint_lattice_fft, ONLY: cell_to_k_grid_fft,&
12 : k_grid_to_cell_execute,&
13 : k_grid_to_cell_fft,&
14 : k_grid_to_cell_prepare,&
15 : k_grid_to_cell_release,&
16 : k_grid_to_cell_work_type,&
17 : regular_kpoint_grid
18 : USE lattice_preconditioner_operator, ONLY: &
19 : apply_lattice_inverse_dense, apply_lattice_state_inverse_dense, build_lattice_inverse, &
20 : build_lattice_local_correction, lattice_fft_auto_select, lattice_fft_cost_ratio, &
21 : lattice_fft_storage_ratio, lattice_grid, lattice_local_cost_ratio, &
22 : lattice_local_storage_ratio
23 : USE mathconstants, ONLY: gaussi,&
24 : twopi,&
25 : z_zero
26 : USE message_passing, ONLY: mp_comm_type,&
27 : mp_world_finalize,&
28 : mp_world_init
29 :
30 : IMPLICIT NONE
31 :
32 : INTEGER, PARAMETER :: base_block = 2, base_cell = 8, base_n = 16, &
33 : base_ncol = 3, ncell = 31, nkp = 64, &
34 : nround = 16, nselected = 3
35 : INTEGER, PARAMETER :: large_block = 1000, large_ncol = 1000, &
36 : large_num_pe = 16
37 2 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: inverse_k
38 2 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: state_inverse_k
39 : COMPLEX(KIND=dp), DIMENSION(2, 3, nkp) :: direct, transformed
40 : COMPLEX(KIND=dp), DIMENSION(2, 3, nround) :: reconstructed, round_k
41 : COMPLEX(KIND=dp), DIMENSION(2, 3, nselected) :: selected_transformed
42 : INTEGER :: d, i, icell, info, j, jcell, ik, ix, iy, iz
43 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: correction_indices
44 : INTEGER, DIMENSION(3) :: base_dims, delta, large_dims
45 : INTEGER, DIMENSION(3, base_cell) :: base_index
46 : INTEGER, DIMENSION(3) :: nkp_grid
47 : INTEGER, DIMENSION(3, ncell) :: index_to_cell
48 : INTEGER, DIMENSION(3, nround) :: round_index
49 : INTEGER, DIMENSION(nselected) :: selected_kpoints
50 : LOGICAL :: used_fft
51 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
52 : REAL(KIND=dp) :: base_residual_norm, corrected_projection_error, &
53 : corrected_residual_norm, cost_ratio, &
54 : expected_cost_ratio, imaginary_residual, &
55 : large_block_real, large_ncell_real, &
56 : large_ncol_real, projection_error
57 : REAL(KIND=dp), DIMENSION(3, base_cell) :: base_xkp
58 : REAL(KIND=dp), DIMENSION(base_n, base_n) :: base_overlap
59 : REAL(KIND=dp), DIMENSION(base_n, base_n) :: defect_operator
60 : REAL(KIND=dp), DIMENSION(base_n, base_ncol) :: base_input, base_input_rotated, &
61 : base_output, base_output_rotated, &
62 : base_state_output, defect_residual
63 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coarse_inverse, defect_base_output, &
64 2 : defect_corrected_output, local_operator, &
65 2 : selected_output, selected_rhs
66 : REAL(KIND=dp), DIMENSION(base_ncol, base_ncol) :: column_rotation
67 : REAL(KIND=dp), DIMENSION(3, nkp) :: xkp, xkp_bad
68 : REAL(KIND=dp), DIMENSION(2, 3, ncell) :: values_rs
69 : REAL(KIND=dp), DIMENSION(2, 3, nround) :: round_rs
70 : TYPE(mp_comm_type) :: mp_comm
71 8 : TYPE(k_grid_to_cell_work_type) :: inverse_work
72 : REAL(KIND=dp), DIMENSION(9) :: subset_weights
73 : COMPLEX(KIND=dp), DIMENSION(2, 3, ncell) :: weighted_reference, weighted_result
74 :
75 2 : CALL mp_world_init(mp_comm)
76 : CALL init_fft("FFTSG", alltoall=.FALSE., pool_limit=4, &
77 2 : wisdom_file="", plan_style=1)
78 :
79 2 : nkp_grid = [4, 4, 4]
80 : hmat = RESHAPE([4.0_dp, 0.1_dp, 0.2_dp, &
81 : 0.0_dp, 5.0_dp, 0.3_dp, &
82 2 : 0.0_dp, 0.0_dp, 6.0_dp], [3, 3])
83 :
84 2 : ik = 0
85 10 : DO iz = 0, nkp_grid(3) - 1
86 42 : DO iy = 0, nkp_grid(2) - 1
87 168 : DO ix = 0, nkp_grid(1) - 1
88 128 : ik = ik + 1
89 : xkp(:, ik) = [(REAL(ix, dp) + 0.5_dp)/REAL(nkp_grid(1), dp), &
90 : REAL(iy, dp)/REAL(nkp_grid(2), dp), &
91 544 : (REAL(iz, dp) + 0.25_dp)/REAL(nkp_grid(3), dp)]
92 : END DO
93 : END DO
94 : END DO
95 2 : IF (.NOT. regular_kpoint_grid(xkp, nkp_grid)) ERROR STOP "Regular k-point grid not recognized"
96 :
97 64 : DO icell = 1, ncell
98 : index_to_cell(:, icell) = [MODULO(3*icell + 1, 11) - 5, &
99 : MODULO(5*icell + 2, 9) - 4, &
100 248 : MODULO(7*icell + 3, 7) - 3]
101 250 : DO d = 1, 3
102 186 : values_rs(1, d, icell) = SIN(0.13_dp*REAL(icell + 2*d, dp))
103 248 : values_rs(2, d, icell) = COS(0.17_dp*REAL(2*icell - d, dp))
104 : END DO
105 : END DO
106 :
107 2 : CALL cell_to_k_grid_fft(values_rs, index_to_cell, xkp, nkp_grid, transformed, used_fft)
108 2 : IF (.NOT. used_fft) ERROR STOP "Complete regular grid did not use lattice FFT"
109 2 : CALL direct_reference(values_rs, index_to_cell, xkp, direct)
110 1282 : IF (MAXVAL(ABS(transformed - direct)) > 2.0E-12_dp) THEN
111 0 : ERROR STOP "Shifted lattice FFT differs from direct transform"
112 : END IF
113 :
114 2 : selected_kpoints = [2, 31, 64]
115 : CALL cell_to_k_grid_fft(values_rs, index_to_cell, xkp, nkp_grid, selected_transformed, &
116 2 : used_fft, selected_kpoints=selected_kpoints)
117 2 : IF (.NOT. used_fft) ERROR STOP "Selected distributed output did not use lattice FFT"
118 62 : IF (MAXVAL(ABS(selected_transformed - direct(:, :, selected_kpoints))) > 2.0E-12_dp) THEN
119 0 : ERROR STOP "Selected lattice-FFT output differs from the full transform"
120 : END IF
121 :
122 2 : CALL cell_to_k_grid_fft(values_rs, index_to_cell, xkp, nkp_grid, transformed, used_fft, 2, hmat)
123 2 : IF (.NOT. used_fft) ERROR STOP "Derivative did not use lattice FFT"
124 2 : CALL direct_reference(values_rs, index_to_cell, xkp, direct, 2, hmat)
125 1282 : IF (MAXVAL(ABS(transformed - direct)) > 5.0E-12_dp) THEN
126 0 : ERROR STOP "Lattice-FFT derivative differs from direct transform"
127 : END IF
128 :
129 2 : nkp_grid = [4, 4, 1]
130 2 : ik = 0
131 10 : DO iy = 0, nkp_grid(2) - 1
132 42 : DO ix = 0, nkp_grid(1) - 1
133 32 : ik = ik + 1
134 : xkp(:, ik) = [(REAL(ix, dp) + 0.5_dp)/REAL(nkp_grid(1), dp), &
135 136 : REAL(iy, dp)/REAL(nkp_grid(2), dp), 0.0_dp]
136 : END DO
137 : END DO
138 : CALL cell_to_k_grid_fft(values_rs, index_to_cell, xkp(:, 1:ik), nkp_grid, &
139 2 : transformed(:, :, 1:ik), used_fft)
140 2 : IF (.NOT. used_fft) ERROR STOP "Singleton grid direction did not use padded lattice FFT"
141 2 : CALL direct_reference(values_rs, index_to_cell, xkp(:, 1:ik), direct(:, :, 1:ik))
142 322 : IF (MAXVAL(ABS(transformed(:, :, 1:ik) - direct(:, :, 1:ik))) > 2.0E-12_dp) THEN
143 0 : ERROR STOP "Padded lattice FFT differs from direct transform"
144 : END IF
145 :
146 : icell = 0
147 10 : DO iy = 0, nkp_grid(2) - 1
148 42 : DO ix = 0, nkp_grid(1) - 1
149 32 : icell = icell + 1
150 128 : round_index(:, icell) = [ix, iy, 0]
151 328 : round_rs(:, :, icell) = values_rs(:, :, icell)
152 : END DO
153 : END DO
154 2 : CALL cell_to_k_grid_fft(round_rs, round_index, xkp(:, 1:ik), nkp_grid, round_k, used_fft)
155 2 : IF (.NOT. used_fft) ERROR STOP "Round-trip forward transform did not use lattice FFT"
156 2 : CALL k_grid_to_cell_fft(round_k, xkp(:, 1:ik), nkp_grid, round_index, reconstructed, used_fft)
157 2 : IF (.NOT. used_fft) ERROR STOP "Round-trip inverse transform did not use lattice FFT"
158 322 : IF (MAXVAL(ABS(reconstructed - CMPLX(round_rs, KIND=dp))) > 3.0E-12_dp) THEN
159 0 : ERROR STOP "Forward/inverse lattice-FFT round trip is inconsistent"
160 : END IF
161 :
162 : ! An incomplete shifted grid uses the supplied weights exactly, including
163 : ! negative translations. Reuse the preparation for separate matrix batches.
164 2 : subset_weights = [(REAL(ik, dp)/45.0_dp, ik=1, 9)]
165 2 : weighted_reference = z_zero
166 64 : DO icell = 1, ncell
167 622 : DO ik = 1, 9
168 : weighted_reference(:, :, icell) = weighted_reference(:, :, icell) + subset_weights(ik)* &
169 7316 : EXP(-gaussi*twopi*SUM(xkp(:, ik)*REAL(index_to_cell(:, icell), dp)))*round_k(:, :, ik)
170 : END DO
171 : END DO
172 : CALL k_grid_to_cell_prepare(inverse_work, xkp(:, 1:9), nkp_grid, index_to_cell, &
173 2 : weights=subset_weights, allow_incomplete=.TRUE.)
174 6 : DO i = 1, 2
175 4 : CALL k_grid_to_cell_execute(inverse_work, round_k(i:i, :, 1:9), weighted_result(i:i, :, :), used_fft)
176 6 : IF (.NOT. used_fft) ERROR STOP "Prepared incomplete grid did not use FFT"
177 : END DO
178 622 : IF (MAXVAL(ABS(weighted_result - weighted_reference)) > 3.0E-12_dp) THEN
179 0 : ERROR STOP "Prepared weighted inverse differs from explicit phase sum"
180 : END IF
181 : CALL k_grid_to_cell_fft(round_k(:, :, 1:9), xkp(:, 1:9), nkp_grid, index_to_cell, weighted_result, &
182 2 : used_fft, weights=subset_weights, allow_incomplete=.TRUE.)
183 622 : IF (MAXVAL(ABS(weighted_result - weighted_reference)) > 3.0E-12_dp) THEN
184 0 : ERROR STOP "One-shot and prepared inverse transforms disagree"
185 : END IF
186 : CALL k_grid_to_cell_prepare(inverse_work, xkp(:, 1:9), nkp_grid, index_to_cell, &
187 20 : weights=2.0_dp*subset_weights, allow_incomplete=.TRUE.)
188 2 : CALL k_grid_to_cell_execute(inverse_work, round_k(:, :, 1:9), weighted_result, used_fft)
189 622 : IF (MAXVAL(ABS(weighted_result - 2.0_dp*weighted_reference)) > 6.0E-12_dp) THEN
190 0 : ERROR STOP "Repreparing inverse transform retained stale weights"
191 : END IF
192 2 : CALL k_grid_to_cell_release(inverse_work)
193 2 : CALL k_grid_to_cell_release(inverse_work)
194 :
195 2 : nkp_grid = [4, 4, 4]
196 :
197 2 : xkp_bad = xkp
198 2 : xkp_bad(1, nkp) = xkp_bad(1, nkp) + 0.03125_dp
199 2 : CALL cell_to_k_grid_fft(values_rs, index_to_cell, xkp_bad, nkp_grid, transformed, used_fft)
200 2 : IF (used_fft) ERROR STOP "Irregular k-point list incorrectly used lattice FFT"
201 2 : CALL direct_reference(values_rs, index_to_cell, xkp_bad, direct)
202 1282 : IF (MAXVAL(ABS(transformed - direct)) > 1.0E-14_dp) THEN
203 0 : ERROR STOP "Direct lattice-transform fallback is inconsistent"
204 : END IF
205 :
206 : ! Verify the block-circulant inverse against direct matrix application and right covariance.
207 2 : base_dims = [4, 2, 1]
208 2 : CALL lattice_grid(base_dims, base_index, base_xkp)
209 2 : base_overlap = 0.0_dp
210 18 : DO jcell = 1, base_cell
211 146 : DO icell = 1, base_cell
212 512 : delta = MODULO(base_index(:, icell) - base_index(:, jcell), base_dims)
213 384 : DO i = 1, base_block
214 512 : IF (ALL(delta == 0)) THEN
215 : base_overlap((icell - 1)*base_block + i, (jcell - 1)*base_block + i) = &
216 32 : 2.8_dp + 0.2_dp*REAL(i, KIND=dp)
217 224 : ELSE IF ((delta(1) == 1 .OR. delta(1) == base_dims(1) - 1) .AND. &
218 : delta(2) == 0) THEN
219 64 : base_overlap((icell - 1)*base_block + i, (jcell - 1)*base_block + i) = -0.12_dp
220 160 : ELSE IF (delta(1) == 0 .AND. delta(2) == 1) THEN
221 32 : base_overlap((icell - 1)*base_block + i, (jcell - 1)*base_block + i) = -0.08_dp
222 : END IF
223 : END DO
224 208 : IF (ALL(delta == 0)) THEN
225 16 : base_overlap((icell - 1)*base_block + 1, (jcell - 1)*base_block + 2) = 0.15_dp
226 16 : base_overlap((icell - 1)*base_block + 2, (jcell - 1)*base_block + 1) = 0.15_dp
227 : END IF
228 : END DO
229 : END DO
230 8 : DO j = 1, base_ncol
231 104 : DO i = 1, base_n
232 102 : base_input(i, j) = SIN(0.11_dp*REAL(2*i + 3*j, KIND=dp))
233 : END DO
234 : END DO
235 2 : CALL build_lattice_inverse(base_overlap, base_dims, inverse_k, projection_error, used_fft, info)
236 2 : IF (info /= 0) ERROR STOP "Positive lattice overlap failed block Cholesky inversion"
237 2 : IF (.NOT. used_fft) ERROR STOP "Block-circulant overlap construction did not use lattice FFT"
238 2 : IF (projection_error > 2.0E-15_dp) ERROR STOP "Exact block-circulant overlap was changed by projection"
239 2 : IF (lattice_fft_auto_select(base_dims, base_block, base_ncol, 1, 1, &
240 : projection_error, used_fft, info)) THEN
241 0 : ERROR STOP "Cost-aware lattice-FFT AUTO accepted a tiny eight-cell operator"
242 : END IF
243 2 : IF (.NOT. lattice_fft_auto_select(base_dims, 6, 8, 1, 2, &
244 : projection_error, used_fft, info)) THEN
245 0 : ERROR STOP "Cost-aware lattice-FFT AUTO rejected an economical eight-cell operator"
246 : END IF
247 2 : IF (lattice_fft_auto_select([4, 1, 1], 6, 8, 1, 2, &
248 : projection_error, used_fft, info)) THEN
249 0 : ERROR STOP "Conservative lattice-FFT AUTO accepted fewer than eight cells"
250 : END IF
251 2 : IF (lattice_fft_auto_select(base_dims, 6, 8, 1, 2, 2.0E-3_dp, used_fft, info)) THEN
252 0 : ERROR STOP "Conservative lattice-FFT AUTO accepted a large projection error"
253 : END IF
254 2 : IF (lattice_fft_auto_select(base_dims, 6, 8, 17, 2, &
255 : projection_error, used_fft, info)) THEN
256 0 : ERROR STOP "Storage-aware lattice-FFT AUTO accepted too many state operators"
257 : END IF
258 2 : IF (lattice_fft_cost_ratio(base_dims, base_block, base_ncol, 1) <= 0.8_dp) THEN
259 0 : ERROR STOP "Lattice-FFT cost estimate missed the tiny-batch overhead"
260 : END IF
261 2 : large_dims = [100, 100, 1]
262 8 : large_ncell_real = REAL(PRODUCT(large_dims), KIND=dp)
263 2 : large_block_real = REAL(large_block, KIND=dp)
264 2 : large_ncol_real = REAL(large_ncol, KIND=dp)
265 : expected_cost_ratio = ( &
266 : large_ncell_real*large_block_real**2*large_ncol_real + &
267 : 4.0_dp*large_ncell_real*large_block_real*large_ncol_real* &
268 : LOG(large_ncell_real)/LOG(2.0_dp) + &
269 : 8.0_dp*large_ncell_real*large_block_real + &
270 : 2.0_dp*large_ncell_real*large_block_real*large_ncol_real* &
271 : LOG(REAL(large_num_pe + 1, KIND=dp))/LOG(2.0_dp))/ &
272 2 : ((large_ncell_real*large_block_real)**2*large_ncol_real)
273 2 : cost_ratio = lattice_fft_cost_ratio(large_dims, large_block, large_ncol, large_num_pe)
274 2 : IF (ABS(cost_ratio - expected_cost_ratio) > 1.0E-14_dp*expected_cost_ratio) THEN
275 0 : ERROR STOP "Lattice-FFT cost estimate overflowed for a large problem"
276 : END IF
277 2 : IF (ABS(lattice_fft_storage_ratio(base_dims, 1) - 0.25_dp) > 1.0E-14_dp) THEN
278 0 : ERROR STOP "Lattice-FFT storage estimate is inconsistent"
279 : END IF
280 2 : IF (ABS(lattice_local_cost_ratio(base_dims, 1) - 0.28125_dp) > 1.0E-14_dp) THEN
281 0 : ERROR STOP "Lattice local-correction cost estimate is inconsistent"
282 : END IF
283 2 : IF (ABS(lattice_local_storage_ratio(base_dims, 3, 1) - 0.421875_dp) > 1.0E-14_dp) THEN
284 0 : ERROR STOP "Lattice local-correction storage estimate is inconsistent"
285 : END IF
286 : CALL apply_lattice_inverse_dense(inverse_k, base_dims, base_input, base_output, &
287 2 : used_fft, imaginary_residual)
288 2 : IF (.NOT. used_fft) ERROR STOP "Block-circulant inverse application did not use lattice FFT"
289 2 : IF (imaginary_residual > 3.0E-13_dp) ERROR STOP "Lattice inverse did not reconstruct a real matrix"
290 1742 : IF (MAXVAL(ABS(MATMUL(base_overlap, base_output) - base_input)) > 2.0E-11_dp) THEN
291 0 : ERROR STOP "FFT block-circulant inverse differs from the dense operator inverse"
292 : END IF
293 :
294 2 : ALLOCATE (state_inverse_k(base_block, base_block, base_cell, base_ncol))
295 8 : DO j = 1, base_ncol
296 344 : state_inverse_k(:, :, :, j) = (0.5_dp + 0.25_dp*REAL(j, KIND=dp))*inverse_k
297 : END DO
298 : CALL apply_lattice_state_inverse_dense(state_inverse_k, base_dims, base_input, &
299 2 : base_state_output, used_fft, imaginary_residual)
300 2 : IF (.NOT. used_fft) ERROR STOP "State-dependent inverse application did not use lattice FFT"
301 8 : DO j = 1, base_ncol
302 102 : IF (MAXVAL(ABS(base_state_output(:, j) - &
303 2 : (0.5_dp + 0.25_dp*REAL(j, KIND=dp))*base_output(:, j))) > 3.0E-12_dp) THEN
304 0 : ERROR STOP "State-dependent lattice inverse used the wrong column operator"
305 : END IF
306 : END DO
307 2 : DEALLOCATE (state_inverse_k)
308 :
309 2 : column_rotation = 0.0_dp
310 2 : column_rotation(1, 1) = COS(0.31_dp)
311 2 : column_rotation(1, 2) = -SIN(0.31_dp)
312 2 : column_rotation(2, 1) = SIN(0.31_dp)
313 2 : column_rotation(2, 2) = COS(0.31_dp)
314 2 : column_rotation(3, 3) = 1.0_dp
315 314 : base_input_rotated = MATMUL(base_input, column_rotation)
316 : CALL apply_lattice_inverse_dense(inverse_k, base_dims, base_input_rotated, &
317 2 : base_output_rotated, used_fft)
318 416 : IF (MAXVAL(ABS(base_output_rotated - MATMUL(base_output, column_rotation))) > 3.0E-12_dp) THEN
319 0 : ERROR STOP "FFT lattice base is not right-covariant under orbital rotations"
320 : END IF
321 :
322 : ! Add one localized positive defect and verify the balanced cell correction.
323 2 : defect_operator = base_overlap
324 2 : defect_operator(5, 5) = defect_operator(5, 5) + 0.90_dp
325 2 : defect_operator(6, 6) = defect_operator(6, 6) + 0.60_dp
326 2 : CALL build_lattice_inverse(defect_operator, base_dims, inverse_k, projection_error, used_fft, info)
327 2 : IF (info /= 0 .OR. .NOT. used_fft) ERROR STOP "Defective lattice operator failed construction"
328 2 : IF (projection_error <= 1.0E-3_dp) ERROR STOP "Localized defect did not break lattice translations"
329 2 : IF (lattice_fft_auto_select(base_dims, 6, 8, 1, 2, projection_error, used_fft, info)) THEN
330 0 : ERROR STOP "Lattice-FFT AUTO accepted a localized defect"
331 : END IF
332 : CALL build_lattice_local_correction( &
333 : defect_operator, base_dims, 1, correction_indices, local_operator, coarse_inverse, &
334 2 : corrected_projection_error, info)
335 2 : IF (info /= 0) ERROR STOP "Local correction coarse inverse failed"
336 6 : IF (ANY(correction_indices /= [5, 6])) ERROR STOP "Local correction selected the wrong defect cell"
337 2 : IF (corrected_projection_error >= projection_error) THEN
338 0 : ERROR STOP "Local correction did not reduce the non-circulant residual"
339 : END IF
340 :
341 2 : ALLOCATE (defect_base_output(base_n, base_ncol), defect_corrected_output(base_n, base_ncol))
342 2 : CALL apply_lattice_inverse_dense(inverse_k, base_dims, base_input, defect_base_output, used_fft)
343 : CALL apply_balanced_local_inverse(inverse_k, base_dims, correction_indices, local_operator, &
344 2 : coarse_inverse, base_input, defect_corrected_output)
345 3280 : defect_residual = MATMUL(defect_operator, defect_base_output) - base_input
346 104 : base_residual_norm = NORM2(defect_residual)
347 1742 : defect_residual = MATMUL(defect_operator, defect_corrected_output) - base_input
348 104 : corrected_residual_norm = NORM2(defect_residual)
349 2 : IF (corrected_residual_norm >= base_residual_norm) THEN
350 0 : ERROR STOP "Balanced local correction did not improve inverse application"
351 : END IF
352 104 : IF (SUM(base_input*defect_corrected_output) <= 0.0_dp) THEN
353 0 : ERROR STOP "Balanced local correction lost positive definiteness"
354 : END IF
355 :
356 0 : ALLOCATE (selected_rhs(base_n, SIZE(correction_indices)), &
357 8 : selected_output(base_n, SIZE(correction_indices)))
358 70 : selected_rhs(:, :) = local_operator
359 : CALL apply_balanced_local_inverse(inverse_k, base_dims, correction_indices, local_operator, &
360 2 : coarse_inverse, selected_rhs, selected_output)
361 2 : selected_rhs = 0.0_dp
362 6 : DO i = 1, SIZE(correction_indices)
363 6 : selected_rhs(correction_indices(i), i) = 1.0_dp
364 : END DO
365 70 : IF (MAXVAL(ABS(selected_output - selected_rhs)) > 3.0E-11_dp) THEN
366 0 : ERROR STOP "Balanced local correction is not exact in its selected subspace"
367 : END IF
368 0 : DEALLOCATE (coarse_inverse, correction_indices, defect_base_output, defect_corrected_output, &
369 2 : local_operator, selected_output, selected_rhs)
370 2 : DEALLOCATE (inverse_k)
371 :
372 2 : CALL finalize_fft(mp_comm, "")
373 2 : CALL mp_world_finalize()
374 :
375 : CONTAINS
376 :
377 : ! **************************************************************************************************
378 : !> \brief Apply the balanced two-level inverse used by the production local correction.
379 : !> \param inverse_k block-circulant base inverse
380 : !> \param lattice_dims replicated-cell dimensions
381 : !> \param correction_indices selected coordinate basis
382 : !> \param operator_columns columns A U
383 : !> \param coarse_inverse inverse of U^T A U
384 : !> \param input right-hand sides
385 : !> \param output corrected inverse action
386 : ! **************************************************************************************************
387 20 : SUBROUTINE apply_balanced_local_inverse(inverse_k, lattice_dims, correction_indices, &
388 4 : operator_columns, coarse_inverse, input, output)
389 :
390 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: inverse_k
391 : INTEGER, DIMENSION(3), INTENT(IN) :: lattice_dims
392 : INTEGER, DIMENSION(:), INTENT(IN) :: correction_indices
393 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: operator_columns, coarse_inverse, input
394 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: output
395 :
396 : INTEGER :: i
397 : LOGICAL :: used_fft
398 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: base_input, coarse_input, &
399 4 : coarse_projected, coarse_residual, &
400 : coarse_rhs
401 :
402 : ALLOCATE (base_input(SIZE(input, 1), SIZE(input, 2)), &
403 : coarse_input(SIZE(correction_indices), SIZE(input, 2)), &
404 : coarse_projected(SIZE(correction_indices), SIZE(input, 2)), &
405 : coarse_residual(SIZE(correction_indices), SIZE(input, 2)), &
406 52 : coarse_rhs(SIZE(correction_indices), SIZE(input, 2)))
407 12 : DO i = 1, SIZE(correction_indices)
408 32 : coarse_input(i, :) = input(correction_indices(i), :)
409 : END DO
410 74 : coarse_rhs(:, :) = MATMUL(coarse_inverse, coarse_input)
411 528 : base_input(:, :) = input - MATMUL(operator_columns, coarse_rhs)
412 4 : CALL apply_lattice_inverse_dense(inverse_k, lattice_dims, base_input, output, used_fft)
413 4 : IF (.NOT. used_fft) ERROR STOP "Balanced correction base did not use lattice FFT"
414 352 : coarse_projected(:, :) = MATMUL(TRANSPOSE(operator_columns), output)
415 74 : coarse_residual(:, :) = MATMUL(coarse_inverse, coarse_projected)
416 34 : output(correction_indices, :) = output(correction_indices, :) + coarse_rhs - coarse_residual
417 4 : DEALLOCATE (base_input, coarse_input, coarse_projected, coarse_residual, coarse_rhs)
418 :
419 2 : END SUBROUTINE apply_balanced_local_inverse
420 :
421 : ! **************************************************************************************************
422 : !> \brief ...
423 : !> \param values_rs ...
424 : !> \param index_to_cell ...
425 : !> \param xkp ...
426 : !> \param values_k ...
427 : !> \param deriv_direction ...
428 : !> \param hmat ...
429 : ! **************************************************************************************************
430 8 : SUBROUTINE direct_reference(values_rs, index_to_cell, xkp, values_k, deriv_direction, hmat)
431 :
432 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: values_rs
433 : INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
434 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp
435 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: values_k
436 : INTEGER, INTENT(IN), OPTIONAL :: deriv_direction
437 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
438 : OPTIONAL :: hmat
439 :
440 : COMPLEX(KIND=dp) :: factor
441 : INTEGER :: icell, ik
442 : REAL(KIND=dp), DIMENSION(3) :: cell_index_real, cell_vector
443 :
444 4168 : values_k = z_zero
445 424 : DO ik = 1, SIZE(xkp, 2)
446 13320 : DO icell = 1, SIZE(values_rs, 3)
447 51584 : factor = EXP(gaussi*twopi*SUM(xkp(:, ik)*REAL(index_to_cell(:, icell), KIND=dp)))
448 12896 : IF (PRESENT(deriv_direction)) THEN
449 15872 : cell_index_real = REAL(index_to_cell(:, icell), KIND=dp)
450 51584 : cell_vector = MATMUL(hmat, cell_index_real)
451 3968 : factor = factor*gaussi*cell_vector(deriv_direction)
452 : END IF
453 129376 : values_k(:, :, ik) = values_k(:, :, ik) + factor*values_rs(:, :, icell)
454 : END DO
455 : END DO
456 :
457 8 : END SUBROUTINE direct_reference
458 :
459 : END PROGRAM kpoint_lattice_fft_unittest
|
