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 computes preconditioners, and implements methods to apply them
10 : !> currently used in qs_ot
11 : !> \par History
12 : !> - [UB] 2009-05-13 Adding stable approximate inverse (full and sparse)
13 : !> \author Joost VandeVondele (09.2002)
14 : ! **************************************************************************************************
15 : MODULE preconditioner_apply
16 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm
17 : USE cp_cfm_types, ONLY: cp_cfm_create,&
18 : cp_cfm_get_info,&
19 : cp_cfm_release,&
20 : cp_cfm_to_fm,&
21 : cp_cfm_type,&
22 : cp_fm_to_cfm
23 : USE cp_dbcsr_api, ONLY: &
24 : dbcsr_add, dbcsr_copy, dbcsr_get_info, dbcsr_init_p, dbcsr_iterator_blocks_left, &
25 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
26 : dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
27 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 : copy_dbcsr_to_local,&
29 : copy_fm_to_dbcsr,&
30 : copy_local_to_dbcsr,&
31 : cp_dbcsr_m_by_n_from_template
32 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
33 : cp_fm_scale_and_add
34 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_restore
35 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
36 : cp_fm_struct_release,&
37 : cp_fm_struct_type
38 : USE cp_fm_types, ONLY: cp_fm_create,&
39 : cp_fm_get_info,&
40 : cp_fm_release,&
41 : cp_fm_to_fm,&
42 : cp_fm_type
43 : USE input_constants, ONLY: &
44 : ot_low_rank_base_lattice_fft, ot_low_rank_base_overlap, ot_precond_fermi_low_rank, &
45 : ot_precond_full_all, ot_precond_full_all_covariant, ot_precond_full_kinetic, &
46 : ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_s_inverse, &
47 : ot_precond_solver_chebyshev, ot_precond_solver_direct, ot_precond_solver_inv_chol, &
48 : ot_precond_solver_update
49 : USE kinds, ONLY: dp
50 : USE lattice_preconditioner_operator, ONLY: apply_lattice_inverse_dense,&
51 : apply_lattice_state_inverse_dense
52 : USE local_gemm_api, ONLY: LOCAL_GEMM_PU_GPU,&
53 : local_gemm_ctxt_type
54 : USE low_rank_preconditioner_model, ONLY: low_rank_inverse_weight
55 : USE mathconstants, ONLY: z_one,&
56 : z_zero
57 : USE parallel_gemm_api, ONLY: parallel_gemm
58 : USE preconditioner_types, ONLY: lattice_rhs_distribution_type,&
59 : preconditioner_p_type,&
60 : preconditioner_type,&
61 : release_lattice_rhs_distribution
62 :
63 : !$ USE OMP_LIB, ONLY: omp_get_max_threads
64 : #include "./base/base_uses.f90"
65 :
66 : IMPLICIT NONE
67 : PRIVATE
68 :
69 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner_apply'
70 :
71 : PUBLIC :: apply_preconditioner_cfm_complex, apply_preconditioner_dbcsr_complex, &
72 : apply_preconditioner_fm, apply_preconditioner_dbcsr, apply_preconditioner_complex_batch
73 :
74 : ! Reused across channels in each worker's batch; no FM/CFM lifecycle in the local path.
75 : TYPE complex_preconditioner_work_type
76 : COMPLEX(dp), ALLOCATABLE :: buffer(:, :, :)
77 : REAL(dp), ALLOCATABLE :: packed(:, :, :)
78 : TYPE(local_gemm_ctxt_type) :: gemm_ctx
79 : END TYPE complex_preconditioner_work_type
80 :
81 : CONTAINS
82 :
83 : ! **************************************************************************************************
84 : !> \brief Apply single-rank complex preconditioners in bounded OpenMP batches.
85 : !> Packing and DBCSR writeback stay outside OpenMP; each worker owns its dense scratch.
86 : !> \param preconditioners channel preconditioners
87 : !> \param matrix_in real/imaginary input pointers, indexed by component and channel
88 : !> \param matrix_out real/imaginary output pointers, indexed by component and channel
89 : !> \param handled false when a channel requires the distributed implementation
90 : !> \param output_scale optional real scale per channel, applied before writeback
91 : ! **************************************************************************************************
92 6363 : SUBROUTINE apply_preconditioner_complex_batch(preconditioners, matrix_in, matrix_out, handled, output_scale)
93 : TYPE(preconditioner_p_type), INTENT(IN) :: preconditioners(:)
94 : TYPE(dbcsr_p_type), INTENT(IN) :: matrix_in(:, :), matrix_out(:, :)
95 : LOGICAL, INTENT(OUT) :: handled
96 : REAL(dp), INTENT(IN), OPTIONAL :: output_scale(:)
97 :
98 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_preconditioner_complex_batch'
99 :
100 : INTEGER :: first, handle, i, j, k, n, nbatch, &
101 : nchannels, nworkers
102 : TYPE(complex_preconditioner_work_type), &
103 6363 : ALLOCATABLE :: work(:)
104 : TYPE(preconditioner_type), POINTER :: p
105 :
106 6363 : handled = .FALSE.
107 6363 : nchannels = SIZE(preconditioners)
108 6363 : IF (nchannels == 0) RETURN
109 9943 : DO i = 1, nchannels
110 7718 : p => preconditioners(i)%preconditioner
111 7718 : IF (.NOT. ASSOCIATED(p)) RETURN
112 7690 : IF (.NOT. ASSOCIATED(p%complex_fm)) RETURN
113 7690 : IF (p%complex_fm%matrix_struct%para_env%num_pe /= 1) RETURN
114 2225 : SELECT CASE (p%in_use)
115 : CASE (ot_precond_full_all, ot_precond_full_all_covariant, ot_precond_fermi_low_rank, &
116 : ot_precond_full_single, ot_precond_full_single_inverse, &
117 : ot_precond_full_kinetic, ot_precond_s_inverse)
118 : CASE DEFAULT
119 3580 : RETURN
120 : END SELECT
121 : END DO
122 :
123 2225 : CALL timeset(routineN, handle)
124 2225 : nworkers = 1
125 2225 : !$ nworkers = MIN(nchannels, omp_get_max_threads())
126 8900 : ALLOCATE (work(nworkers))
127 4450 : DO j = 1, nworkers
128 4450 : CALL work(j)%gemm_ctx%create(LOCAL_GEMM_PU_GPU, timing=.FALSE.)
129 : END DO
130 5805 : DO first = 1, nchannels, nworkers
131 3580 : nbatch = MIN(nworkers, nchannels - first + 1)
132 7160 : DO j = 1, nbatch
133 3580 : i = first + j - 1
134 3580 : p => preconditioners(i)%preconditioner
135 3580 : CALL dbcsr_get_info(matrix_in(1, i)%matrix, nfullrows_total=n, nfullcols_total=k)
136 3580 : CPASSERT(n == p%complex_fm%matrix_struct%nrow_global)
137 3580 : IF (ALLOCATED(work(j)%buffer)) THEN
138 1355 : IF (SIZE(work(j)%buffer, 1) /= n .OR. SIZE(work(j)%buffer, 2) /= k) THEN
139 54 : DEALLOCATE (work(j)%buffer, work(j)%packed)
140 : END IF
141 : END IF
142 3580 : IF (.NOT. ALLOCATED(work(j)%buffer)) THEN
143 18232 : ALLOCATE (work(j)%buffer(n, k, 2), work(j)%packed(n, k, 2))
144 : END IF
145 3580 : CALL copy_dbcsr_to_local(matrix_in(1, i)%matrix, work(j)%packed(:, :, 1))
146 10740 : CALL copy_dbcsr_to_local(matrix_in(2, i)%matrix, work(j)%packed(:, :, 2))
147 : END DO
148 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) NUM_THREADS(nbatch) IF(nbatch > 1) &
149 3580 : !$OMP SHARED(work,preconditioners,first,nbatch,output_scale) PRIVATE(i)
150 : DO j = 1, nbatch
151 : i = first + j - 1
152 : CALL apply_preconditioner_complex_local(preconditioners(i)%preconditioner, work(j))
153 : IF (PRESENT(output_scale)) THEN
154 : IF (output_scale(i) /= 1.0_dp) work(j)%packed(:, :, :) = output_scale(i)*work(j)%packed(:, :, :)
155 : END IF
156 : END DO
157 : !$OMP END PARALLEL DO
158 9385 : DO j = 1, nbatch
159 3580 : i = first + j - 1
160 3580 : CALL copy_local_to_dbcsr(matrix_out(1, i)%matrix, work(j)%packed(:, :, 1))
161 7160 : CALL copy_local_to_dbcsr(matrix_out(2, i)%matrix, work(j)%packed(:, :, 2))
162 : END DO
163 : END DO
164 4450 : DO j = 1, nworkers
165 4450 : CALL work(j)%gemm_ctx%destroy()
166 : END DO
167 2225 : handled = .TRUE.
168 2225 : CALL timestop(handle)
169 10813 : END SUBROUTINE apply_preconditioner_complex_batch
170 :
171 : ! **************************************************************************************************
172 : !> \brief Apply the existing complex preconditioner with private local GEMM storage.
173 : !> \param p channel preconditioner; its matrices are read only
174 : !> \param w worker storage
175 : ! **************************************************************************************************
176 3580 : SUBROUTINE apply_preconditioner_complex_local(p, w)
177 : TYPE(preconditioner_type), INTENT(IN) :: p
178 : TYPE(complex_preconditioner_work_type), &
179 : INTENT(INOUT), TARGET :: w
180 :
181 3580 : COMPLEX(dp), CONTIGUOUS, POINTER :: basis(:, :), RESULT(:, :), rhs(:, :), &
182 3580 : rotation(:, :)
183 : INTEGER :: k, n
184 :
185 3580 : NULLIFY (rotation)
186 3580 : n = SIZE(w%buffer, 1)
187 3580 : k = SIZE(w%buffer, 2)
188 0 : IF (n == 0 .OR. k == 0) RETURN
189 3580 : basis => p%complex_fm%local_data
190 3580 : rhs => w%buffer(:, :, 1)
191 3580 : RESULT => w%buffer(:, :, 2)
192 1087988 : rhs(:, :) = CMPLX(w%packed(:, :, 1), w%packed(:, :, 2), KIND=dp)
193 4876 : SELECT CASE (p%in_use)
194 : CASE (ot_precond_full_all, ot_precond_full_all_covariant, ot_precond_fermi_low_rank)
195 1296 : CPASSERT(ASSOCIATED(p%full_evals))
196 1296 : CPASSERT(ASSOCIATED(p%occ_evals))
197 1296 : CPASSERT(n == SIZE(p%full_evals))
198 1296 : IF (p%in_use /= ot_precond_fermi_low_rank) THEN
199 1296 : CPASSERT(k == SIZE(p%occ_evals))
200 : END IF
201 1296 : IF (p%in_use == ot_precond_full_all_covariant) THEN
202 108 : CPASSERT(ASSOCIATED(p%occ_rotation_complex))
203 108 : rotation => p%occ_rotation_complex%local_data
204 : CALL w%gemm_ctx%gemm('N', 'N', n, k, k, z_one, rhs, n, &
205 108 : rotation, SIZE(rotation, 1), z_zero, RESULT, n)
206 108 : rhs => w%buffer(:, :, 2)
207 108 : RESULT => w%buffer(:, :, 1)
208 : END IF
209 : CALL w%gemm_ctx%gemm('C', 'N', n, k, n, z_one, basis, SIZE(basis, 1), &
210 1296 : rhs, n, z_zero, RESULT, n)
211 : CALL scale_complex_spectral(p, RESULT, p%complex_fm%matrix_struct%row_indices, &
212 1296 : p%complex_fm%matrix_struct%col_indices)
213 : CALL w%gemm_ctx%gemm('N', 'N', n, k, n, z_one, basis, SIZE(basis, 1), &
214 1296 : RESULT, n, z_zero, rhs, n)
215 1296 : IF (p%in_use == ot_precond_full_all_covariant) THEN
216 : CALL w%gemm_ctx%gemm('N', 'C', n, k, k, z_one, rhs, n, &
217 108 : rotation, SIZE(rotation, 1), z_zero, RESULT, n)
218 : ELSE
219 1188 : RESULT => rhs
220 : END IF
221 : CASE DEFAULT
222 : CALL w%gemm_ctx%gemm('N', 'N', n, k, n, z_one, basis, SIZE(basis, 1), &
223 3580 : rhs, n, z_zero, RESULT, n)
224 : END SELECT
225 1087988 : w%packed(:, :, 1) = REAL(RESULT, KIND=dp)
226 1087988 : w%packed(:, :, 2) = AIMAG(RESULT)
227 3580 : END SUBROUTINE apply_preconditioner_complex_local
228 :
229 : ! **************************************************************************************************
230 : !> \brief Scale spectral coefficients identically for local and distributed applications.
231 : !> \param p spectral preconditioner
232 : !> \param DATA local spectral coefficients
233 : !> \param rows global row indices
234 : !> \param cols global column indices
235 : ! **************************************************************************************************
236 12471 : SUBROUTINE scale_complex_spectral(p, DATA, rows, cols)
237 : TYPE(preconditioner_type), INTENT(IN) :: p
238 : COMPLEX(dp), INTENT(INOUT) :: data(:, :)
239 : INTEGER, INTENT(IN) :: rows(:), cols(:)
240 :
241 : INTEGER :: i, j, nocc, rank_used
242 : REAL(dp) :: scale
243 :
244 12471 : IF (p%in_use == ot_precond_fermi_low_rank) THEN
245 252 : nocc = SIZE(p%occ_evals)
246 252 : rank_used = MIN(p%spectral_rank, SIZE(p%full_evals) - nocc)
247 1498 : DO i = 1, SIZE(DATA, 1)
248 1246 : scale = p%spectral_base_scale
249 1246 : IF (rows(i) > nocc .AND. rows(i) <= nocc + rank_used) THEN
250 196 : scale = low_rank_inverse_weight(p%full_evals(rows(i)), p%spectral_reference, p%energy_gap)
251 : END IF
252 11620 : DATA(i, :) = scale*DATA(i, :)
253 : END DO
254 : ELSE
255 214958 : DO j = 1, SIZE(DATA, 2)
256 6938897 : DO i = 1, SIZE(DATA, 1)
257 6723939 : scale = 1.0_dp/MAX(p%energy_gap, p%full_evals(rows(i)) - p%occ_evals(cols(j)))
258 6926678 : DATA(i, j) = scale*DATA(i, j)
259 : END DO
260 : END DO
261 : END IF
262 12471 : END SUBROUTINE scale_complex_spectral
263 :
264 : ! **************************************************************************************************
265 : !> \brief applies a previously created preconditioner to a full matrix
266 : !> \param preconditioner_env ...
267 : !> \param matrix_in ...
268 : !> \param matrix_out ...
269 : ! **************************************************************************************************
270 48956 : SUBROUTINE apply_preconditioner_fm(preconditioner_env, matrix_in, matrix_out)
271 :
272 : TYPE(preconditioner_type) :: preconditioner_env
273 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
274 :
275 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_fm'
276 :
277 : INTEGER :: handle
278 :
279 48956 : CALL timeset(routineN, handle)
280 :
281 48956 : SELECT CASE (preconditioner_env%in_use)
282 : CASE (0)
283 0 : CPABORT("No preconditioner in use")
284 : CASE (ot_precond_full_single)
285 1352 : CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
286 : CASE (ot_precond_full_all)
287 27458 : CALL apply_full_all(preconditioner_env, matrix_in, matrix_out)
288 : CASE (ot_precond_fermi_low_rank)
289 6 : CALL apply_low_rank_fm(preconditioner_env, matrix_in, matrix_out)
290 : CASE (ot_precond_full_all_covariant)
291 8 : CALL apply_full_all_covariant(preconditioner_env, matrix_in, matrix_out)
292 : CASE (ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
293 20132 : IF (preconditioner_env%lattice_fft_active) THEN
294 0 : CALL apply_lattice_operator_fm(preconditioner_env, matrix_in, matrix_out)
295 : ELSE
296 40264 : SELECT CASE (preconditioner_env%solver)
297 : CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
298 20132 : CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
299 : CASE (ot_precond_solver_direct)
300 0 : CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
301 : CASE (ot_precond_solver_chebyshev)
302 0 : CALL apply_chebyshev_fm(preconditioner_env, matrix_in, matrix_out)
303 : CASE DEFAULT
304 20132 : CPABORT("Solver not implemented")
305 : END SELECT
306 : END IF
307 : CASE DEFAULT
308 48956 : CPABORT("Unknown preconditioner")
309 : END SELECT
310 :
311 48956 : CALL timestop(handle)
312 :
313 48956 : END SUBROUTINE apply_preconditioner_fm
314 :
315 : ! **************************************************************************************************
316 : !> \brief ...
317 : !> \param preconditioner_env ...
318 : !> \param matrix_in ...
319 : !> \param matrix_out ...
320 : ! **************************************************************************************************
321 80475 : SUBROUTINE apply_preconditioner_dbcsr(preconditioner_env, matrix_in, matrix_out)
322 :
323 : TYPE(preconditioner_type) :: preconditioner_env
324 : TYPE(dbcsr_type) :: matrix_in, matrix_out
325 :
326 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_dbcsr'
327 :
328 : INTEGER :: handle
329 :
330 80475 : CALL timeset(routineN, handle)
331 :
332 80475 : SELECT CASE (preconditioner_env%in_use)
333 : CASE (0)
334 0 : CPABORT("No preconditioner in use")
335 : CASE (ot_precond_full_single)
336 202 : CALL apply_single(preconditioner_env, matrix_in, matrix_out)
337 : CASE (ot_precond_full_all)
338 25376 : CALL apply_all(preconditioner_env, matrix_in, matrix_out)
339 : CASE (ot_precond_fermi_low_rank)
340 0 : CALL apply_low_rank_dbcsr(preconditioner_env, matrix_in, matrix_out)
341 : CASE (ot_precond_full_all_covariant)
342 0 : CALL apply_all_covariant(preconditioner_env, matrix_in, matrix_out)
343 : CASE (ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
344 54897 : IF (preconditioner_env%lattice_fft_active) THEN
345 156 : CALL apply_lattice_operator_dbcsr(preconditioner_env, matrix_in, matrix_out)
346 : ELSE
347 109482 : SELECT CASE (preconditioner_env%solver)
348 : CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
349 54741 : CALL apply_single(preconditioner_env, matrix_in, matrix_out)
350 : CASE (ot_precond_solver_direct)
351 0 : CPABORT("Apply_full_direct not supported with ot")
352 : !CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
353 : CASE (ot_precond_solver_chebyshev)
354 0 : CALL apply_chebyshev_dbcsr(preconditioner_env, matrix_in, matrix_out)
355 : CASE DEFAULT
356 54741 : CPABORT("Wrong solver")
357 : END SELECT
358 : END IF
359 : CASE DEFAULT
360 80475 : CPABORT("Wrong preconditioner")
361 : END SELECT
362 :
363 80475 : CALL timestop(handle)
364 :
365 80475 : END SUBROUTINE apply_preconditioner_dbcsr
366 :
367 : ! **************************************************************************************************
368 : !> \brief Apply a complex k-point orbital preconditioner.
369 : !> \param preconditioner_env complex preconditioner storage
370 : !> \param matrix_in complex input channel
371 : !> \param matrix_out complex output channel
372 : ! **************************************************************************************************
373 45111 : SUBROUTINE apply_preconditioner_cfm_complex(preconditioner_env, matrix_in, matrix_out)
374 :
375 : TYPE(preconditioner_type) :: preconditioner_env
376 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_in
377 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_out
378 :
379 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_cfm_complex'
380 :
381 : INTEGER :: handle, k, n, npre
382 : TYPE(cp_cfm_type) :: canonical_in, canonical_out
383 :
384 15037 : CALL timeset(routineN, handle)
385 :
386 15037 : SELECT CASE (preconditioner_env%in_use)
387 : CASE (ot_precond_fermi_low_rank, ot_precond_full_all, ot_precond_full_all_covariant, &
388 : ot_precond_full_single, ot_precond_full_single_inverse, &
389 : ot_precond_full_kinetic, ot_precond_s_inverse)
390 : CASE DEFAULT
391 15037 : CPABORT("Unsupported complex K-point OT preconditioner")
392 : END SELECT
393 15037 : CPASSERT(ASSOCIATED(preconditioner_env%complex_fm))
394 15037 : CALL cp_cfm_get_info(matrix_in, nrow_global=n, ncol_global=k)
395 15037 : CALL cp_cfm_get_info(preconditioner_env%complex_fm, nrow_global=npre)
396 15037 : CPASSERT(n == npre)
397 :
398 26128 : SELECT CASE (preconditioner_env%in_use)
399 : CASE (ot_precond_full_all, ot_precond_fermi_low_rank)
400 11091 : CALL apply_complex_spectral(preconditioner_env, matrix_in, matrix_out)
401 : CASE (ot_precond_full_all_covariant)
402 84 : CPASSERT(ASSOCIATED(preconditioner_env%occ_rotation_complex))
403 : CALL cp_cfm_create(canonical_in, matrix_in%matrix_struct, &
404 84 : name='complex covariant FULL_ALL input')
405 : CALL cp_cfm_create(canonical_out, matrix_out%matrix_struct, &
406 84 : name='complex covariant FULL_ALL output')
407 : CALL cp_cfm_gemm('N', 'N', n, k, k, z_one, matrix_in, &
408 84 : preconditioner_env%occ_rotation_complex, z_zero, canonical_in)
409 84 : CALL apply_complex_spectral(preconditioner_env, canonical_in, canonical_out)
410 : CALL cp_cfm_gemm('N', 'C', n, k, k, z_one, canonical_out, &
411 84 : preconditioner_env%occ_rotation_complex, z_zero, matrix_out)
412 84 : CALL cp_cfm_release(canonical_out)
413 84 : CALL cp_cfm_release(canonical_in)
414 : CASE DEFAULT
415 : CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, preconditioner_env%complex_fm, &
416 15121 : matrix_in, z_zero, matrix_out)
417 : END SELECT
418 :
419 15037 : CALL timestop(handle)
420 :
421 15037 : END SUBROUTINE apply_preconditioner_cfm_complex
422 :
423 : ! **************************************************************************************************
424 : !> \brief Apply the spectral inverse in its canonical complex basis.
425 : !> \param preconditioner_env complex spectral preconditioner
426 : !> \param matrix_in input orbital-gradient matrix
427 : !> \param matrix_out preconditioned orbital-gradient matrix
428 : ! **************************************************************************************************
429 44700 : SUBROUTINE apply_complex_spectral(preconditioner_env, matrix_in, matrix_out)
430 :
431 : TYPE(preconditioner_type) :: preconditioner_env
432 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_in
433 : TYPE(cp_cfm_type), INTENT(INOUT) :: matrix_out
434 :
435 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
436 11175 : POINTER :: local_data
437 : INTEGER :: k, n, ncol_local, nrow_local
438 11175 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
439 : TYPE(cp_cfm_type) :: matrix_spectral
440 :
441 0 : CPASSERT(ASSOCIATED(preconditioner_env%full_evals))
442 11175 : CPASSERT(ASSOCIATED(preconditioner_env%occ_evals))
443 11175 : CALL cp_cfm_get_info(matrix_in, nrow_global=n, ncol_global=k)
444 11175 : CPASSERT(n == SIZE(preconditioner_env%full_evals))
445 11175 : IF (preconditioner_env%in_use /= ot_precond_fermi_low_rank) THEN
446 10923 : CPASSERT(k == SIZE(preconditioner_env%occ_evals))
447 : END IF
448 : CALL cp_cfm_create(matrix_spectral, matrix_in%matrix_struct, &
449 11175 : name='complex spectral input')
450 : CALL cp_cfm_gemm('C', 'N', n, k, n, z_one, preconditioner_env%complex_fm, &
451 11175 : matrix_in, z_zero, matrix_spectral)
452 : CALL cp_cfm_get_info(matrix_spectral, nrow_local=nrow_local, ncol_local=ncol_local, &
453 11175 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
454 : CALL scale_complex_spectral(preconditioner_env, local_data(1:nrow_local, 1:ncol_local), &
455 11175 : row_indices, col_indices)
456 : CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, preconditioner_env%complex_fm, &
457 11175 : matrix_spectral, z_zero, matrix_out)
458 11175 : CALL cp_cfm_release(matrix_spectral)
459 :
460 11175 : END SUBROUTINE apply_complex_spectral
461 :
462 : ! **************************************************************************************************
463 : !> \brief Apply a complex orbital preconditioner to paired real/imaginary DBCSR matrices.
464 : !> \param preconditioner_env complex k-point preconditioner
465 : !> \param matrix_in_re real input channel
466 : !> \param matrix_in_im imaginary input channel
467 : !> \param matrix_out_re real output channel
468 : !> \param matrix_out_im imaginary output channel
469 : ! **************************************************************************************************
470 4075 : SUBROUTINE apply_preconditioner_dbcsr_complex(preconditioner_env, matrix_in_re, matrix_in_im, &
471 : matrix_out_re, matrix_out_im)
472 :
473 : TYPE(preconditioner_type), TARGET :: preconditioner_env
474 : TYPE(dbcsr_type), TARGET :: matrix_in_re, matrix_in_im, &
475 : matrix_out_re, matrix_out_im
476 :
477 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_dbcsr_complex'
478 :
479 : INTEGER :: handle, k, n
480 : LOGICAL :: handled
481 : TYPE(cp_cfm_type) :: matrix_in, matrix_out
482 : TYPE(cp_fm_type) :: matrix_in_im_fm, matrix_in_re_fm, &
483 : matrix_out_im_fm, matrix_out_re_fm
484 32600 : TYPE(dbcsr_p_type) :: input(2, 1), output(2, 1)
485 8150 : TYPE(preconditioner_p_type) :: preconditioners(1)
486 :
487 4075 : preconditioners(1)%preconditioner => preconditioner_env
488 4075 : input(1, 1)%matrix => matrix_in_re
489 4075 : input(2, 1)%matrix => matrix_in_im
490 4075 : output(1, 1)%matrix => matrix_out_re
491 4075 : output(2, 1)%matrix => matrix_out_im
492 4075 : CALL apply_preconditioner_complex_batch(preconditioners, input, output, handled)
493 4075 : IF (handled) RETURN
494 :
495 3290 : CALL timeset(routineN, handle)
496 :
497 3290 : CALL dbcsr_get_info(matrix_in_re, nfullrows_total=n, nfullcols_total=k)
498 3290 : CPASSERT(ASSOCIATED(preconditioner_env%complex_fm))
499 :
500 : CALL cp_fm_create(matrix_in_re_fm, preconditioner_env%complex_fm%matrix_struct, &
501 3290 : nrow=n, ncol=k, name='complex preconditioner input real')
502 : CALL cp_fm_create(matrix_in_im_fm, preconditioner_env%complex_fm%matrix_struct, &
503 3290 : nrow=n, ncol=k, name='complex preconditioner input imaginary')
504 : CALL cp_fm_create(matrix_out_re_fm, preconditioner_env%complex_fm%matrix_struct, &
505 3290 : nrow=n, ncol=k, name='complex preconditioner output real')
506 : CALL cp_fm_create(matrix_out_im_fm, preconditioner_env%complex_fm%matrix_struct, &
507 3290 : nrow=n, ncol=k, name='complex preconditioner output imaginary')
508 3290 : CALL copy_dbcsr_to_fm(matrix_in_re, matrix_in_re_fm)
509 3290 : CALL copy_dbcsr_to_fm(matrix_in_im, matrix_in_im_fm)
510 :
511 : CALL cp_cfm_create(matrix_in, matrix_in_re_fm%matrix_struct, &
512 3290 : name='complex preconditioner input')
513 : CALL cp_cfm_create(matrix_out, matrix_in_re_fm%matrix_struct, &
514 3290 : name='complex preconditioner output')
515 3290 : CALL cp_fm_to_cfm(matrix_in_re_fm, matrix_in_im_fm, matrix_in)
516 3290 : CALL apply_preconditioner_cfm_complex(preconditioner_env, matrix_in, matrix_out)
517 3290 : CALL cp_cfm_to_fm(matrix_out, matrix_out_re_fm, matrix_out_im_fm)
518 3290 : CALL dbcsr_set(matrix_out_re, 0.0_dp)
519 3290 : CALL dbcsr_set(matrix_out_im, 0.0_dp)
520 3290 : CALL copy_fm_to_dbcsr(matrix_out_re_fm, matrix_out_re)
521 3290 : CALL copy_fm_to_dbcsr(matrix_out_im_fm, matrix_out_im)
522 :
523 3290 : CALL cp_cfm_release(matrix_out)
524 3290 : CALL cp_cfm_release(matrix_in)
525 3290 : CALL cp_fm_release(matrix_out_im_fm)
526 3290 : CALL cp_fm_release(matrix_out_re_fm)
527 3290 : CALL cp_fm_release(matrix_in_im_fm)
528 3290 : CALL cp_fm_release(matrix_in_re_fm)
529 :
530 3290 : CALL timestop(handle)
531 :
532 4075 : END SUBROUTINE apply_preconditioner_dbcsr_complex
533 :
534 : ! **************************************************************************************************
535 : !> \brief apply to full matrix, complete inversion has already been done
536 : !> \param preconditioner_env ...
537 : !> \param matrix_in ...
538 : !> \param matrix_out ...
539 : ! **************************************************************************************************
540 42968 : SUBROUTINE apply_full_single(preconditioner_env, matrix_in, matrix_out)
541 :
542 : TYPE(preconditioner_type) :: preconditioner_env
543 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
544 :
545 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_single'
546 :
547 : INTEGER :: handle, k, n
548 :
549 21484 : CALL timeset(routineN, handle)
550 :
551 21484 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
552 : CALL parallel_gemm('N', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
553 21484 : matrix_in, 0.0_dp, matrix_out)
554 21484 : CALL timestop(handle)
555 :
556 21484 : END SUBROUTINE apply_full_single
557 :
558 : ! **************************************************************************************************
559 : !> \brief apply to dbcsr matrix, complete inversion has already been done
560 : !> \param preconditioner_env ...
561 : !> \param matrix_in ...
562 : !> \param matrix_out ...
563 : ! **************************************************************************************************
564 54943 : SUBROUTINE apply_single(preconditioner_env, matrix_in, matrix_out)
565 :
566 : TYPE(preconditioner_type) :: preconditioner_env
567 : TYPE(dbcsr_type) :: matrix_in, matrix_out
568 :
569 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_single'
570 :
571 : INTEGER :: handle
572 :
573 54943 : CALL timeset(routineN, handle)
574 :
575 54943 : IF (.NOT. ASSOCIATED(preconditioner_env%dbcsr_matrix)) THEN
576 0 : CPABORT("NOT ASSOCIATED preconditioner_env%dbcsr_matrix")
577 : END IF
578 : CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, matrix_in, &
579 54943 : 0.0_dp, matrix_out)
580 :
581 54943 : CALL timestop(handle)
582 :
583 54943 : END SUBROUTINE apply_single
584 :
585 : ! **************************************************************************************************
586 : !> \brief preconditioner contains the factorization, application done by
587 : !> solving the linear system
588 : !> \param preconditioner_env ...
589 : !> \param matrix_in ...
590 : !> \param matrix_out ...
591 : ! **************************************************************************************************
592 0 : SUBROUTINE apply_full_direct(preconditioner_env, matrix_in, matrix_out)
593 :
594 : TYPE(preconditioner_type) :: preconditioner_env
595 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
596 :
597 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_direct'
598 :
599 : INTEGER :: handle, k, n
600 : TYPE(cp_fm_type) :: work
601 :
602 0 : CALL timeset(routineN, handle)
603 :
604 0 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
605 0 : CALL cp_fm_create(work, matrix_in%matrix_struct, name="apply_full_single")
606 : CALL cp_fm_cholesky_restore(matrix_in, k, preconditioner_env%fm, work,&
607 0 : & "SOLVE", transa="T")
608 : CALL cp_fm_cholesky_restore(work, k, preconditioner_env%fm, matrix_out,&
609 0 : & "SOLVE", transa="N")
610 0 : CALL cp_fm_release(work)
611 :
612 0 : CALL timestop(handle)
613 :
614 0 : END SUBROUTINE apply_full_direct
615 :
616 : ! **************************************************************************************************
617 : !> \brief Apply a Chebyshev approximation to the inverse of the stored dense SPD operator.
618 : !> \param preconditioner_env ...
619 : !> \param matrix_in ...
620 : !> \param matrix_out ...
621 : ! **************************************************************************************************
622 0 : SUBROUTINE apply_chebyshev_fm(preconditioner_env, matrix_in, matrix_out)
623 :
624 : TYPE(preconditioner_type) :: preconditioner_env
625 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
626 :
627 : INTEGER :: degree, iteration, k, n
628 : REAL(KIND=dp) :: delta, rho, rho_previous, sigma, theta
629 : TYPE(cp_fm_type) :: direction, residual
630 :
631 0 : CPASSERT(ASSOCIATED(preconditioner_env%fm))
632 0 : degree = preconditioner_env%polynomial_degree
633 0 : CPASSERT(degree >= 1)
634 0 : theta = 0.5_dp*(preconditioner_env%polynomial_max + preconditioner_env%polynomial_min)
635 0 : delta = 0.5_dp*(preconditioner_env%polynomial_max - preconditioner_env%polynomial_min)
636 0 : sigma = theta/delta
637 0 : rho_previous = 1.0_dp/sigma
638 0 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
639 0 : CALL cp_fm_create(direction, matrix_in%matrix_struct, name='Chebyshev direction')
640 0 : CALL cp_fm_create(residual, matrix_in%matrix_struct, name='Chebyshev residual')
641 0 : CALL cp_fm_to_fm(matrix_in, direction)
642 0 : CALL cp_fm_scale(1.0_dp/theta, direction)
643 0 : CALL cp_fm_to_fm(direction, matrix_out)
644 0 : DO iteration = 2, degree
645 0 : CALL cp_fm_to_fm(matrix_in, residual)
646 : CALL parallel_gemm('N', 'N', n, k, n, -1.0_dp, preconditioner_env%fm, &
647 0 : matrix_out, 1.0_dp, residual)
648 0 : rho = 1.0_dp/(2.0_dp*sigma - rho_previous)
649 0 : CALL cp_fm_scale(rho*rho_previous, direction)
650 0 : CALL cp_fm_scale_and_add(1.0_dp, direction, 2.0_dp*rho/delta, residual)
651 0 : CALL cp_fm_scale_and_add(1.0_dp, matrix_out, 1.0_dp, direction)
652 0 : rho_previous = rho
653 : END DO
654 0 : CALL cp_fm_release(residual)
655 0 : CALL cp_fm_release(direction)
656 :
657 0 : END SUBROUTINE apply_chebyshev_fm
658 :
659 : ! **************************************************************************************************
660 : !> \brief Apply a Chebyshev approximation to the inverse of the stored sparse SPD operator.
661 : !> \param preconditioner_env ...
662 : !> \param matrix_in ...
663 : !> \param matrix_out ...
664 : ! **************************************************************************************************
665 0 : SUBROUTINE apply_chebyshev_dbcsr(preconditioner_env, matrix_in, matrix_out)
666 :
667 : TYPE(preconditioner_type) :: preconditioner_env
668 : TYPE(dbcsr_type) :: matrix_in, matrix_out
669 :
670 : INTEGER :: degree, iteration
671 : REAL(KIND=dp) :: delta, rho, rho_previous, sigma, theta
672 : TYPE(dbcsr_type) :: direction, residual
673 :
674 0 : CPASSERT(ASSOCIATED(preconditioner_env%dbcsr_matrix))
675 0 : degree = preconditioner_env%polynomial_degree
676 0 : CPASSERT(degree >= 1)
677 0 : theta = 0.5_dp*(preconditioner_env%polynomial_max + preconditioner_env%polynomial_min)
678 0 : delta = 0.5_dp*(preconditioner_env%polynomial_max - preconditioner_env%polynomial_min)
679 0 : sigma = theta/delta
680 0 : rho_previous = 1.0_dp/sigma
681 0 : CALL dbcsr_copy(direction, matrix_in, name='Chebyshev direction')
682 0 : CALL dbcsr_copy(residual, matrix_in, name='Chebyshev residual')
683 0 : CALL dbcsr_set(matrix_out, 0.0_dp)
684 0 : CALL dbcsr_add(matrix_out, direction, 1.0_dp, 1.0_dp/theta)
685 0 : CALL dbcsr_set(direction, 0.0_dp)
686 0 : CALL dbcsr_add(direction, matrix_in, 1.0_dp, 1.0_dp/theta)
687 0 : DO iteration = 2, degree
688 0 : CALL dbcsr_set(residual, 0.0_dp)
689 0 : CALL dbcsr_add(residual, matrix_in, 1.0_dp, 1.0_dp)
690 : CALL dbcsr_multiply('N', 'N', -1.0_dp, preconditioner_env%dbcsr_matrix, &
691 0 : matrix_out, 1.0_dp, residual)
692 0 : rho = 1.0_dp/(2.0_dp*sigma - rho_previous)
693 0 : CALL dbcsr_add(direction, residual, rho*rho_previous, 2.0_dp*rho/delta)
694 0 : CALL dbcsr_add(matrix_out, direction, 1.0_dp, 1.0_dp)
695 0 : rho_previous = rho
696 : END DO
697 0 : CALL dbcsr_release(residual)
698 0 : CALL dbcsr_release(direction)
699 :
700 0 : END SUBROUTINE apply_chebyshev_dbcsr
701 :
702 : ! **************************************************************************************************
703 : !> \brief Applies the overlap-inverse base and a bounded low-rank spectral correction.
704 : !> \param preconditioner_env ...
705 : !> \param matrix_in ...
706 : !> \param matrix_out ...
707 : ! **************************************************************************************************
708 12 : SUBROUTINE apply_low_rank_fm(preconditioner_env, matrix_in, matrix_out)
709 :
710 : TYPE(preconditioner_type) :: preconditioner_env
711 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
712 :
713 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_low_rank_fm'
714 :
715 : INTEGER :: handle, i, j, k, n, ncol_local, nocc, &
716 : nrow_local, rank_used
717 6 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
718 : REAL(KIND=dp) :: correction
719 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
720 6 : POINTER :: local_data
721 : TYPE(cp_fm_type) :: matrix_tmp
722 :
723 6 : CALL timeset(routineN, handle)
724 :
725 6 : CPASSERT(ASSOCIATED(preconditioner_env%fm))
726 6 : CPASSERT(ASSOCIATED(preconditioner_env%full_evals))
727 6 : CPASSERT(ASSOCIATED(preconditioner_env%occ_evals))
728 :
729 6 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
730 6 : nocc = SIZE(preconditioner_env%occ_evals)
731 6 : CPASSERT(n == SIZE(preconditioner_env%full_evals))
732 6 : rank_used = MIN(preconditioner_env%spectral_rank, n - nocc)
733 :
734 : ! With the exact base V V^T = S^-1. The optional lattice base replaces
735 : ! S^-1 by the inverse of its block-circulant translational projection:
736 : ! P = B^-1/window + V_r diag(1/max(gap, epsilon_a-mu) - 1/window) V_r^T.
737 : ! The common reference mu makes the operator right-covariant for any number
738 : ! of orbital columns. The lattice path retains only non-negative corrections.
739 12 : SELECT CASE (preconditioner_env%low_rank_base)
740 : CASE (ot_low_rank_base_overlap)
741 6 : CPASSERT(ASSOCIATED(preconditioner_env%base_fm))
742 : CALL parallel_gemm('N', 'N', n, k, n, preconditioner_env%spectral_base_scale, &
743 6 : preconditioner_env%base_fm, matrix_in, 0.0_dp, matrix_out)
744 : CASE (ot_low_rank_base_lattice_fft)
745 : CALL apply_lattice_operator_fm(preconditioner_env, matrix_in, matrix_out, &
746 0 : preconditioner_env%spectral_base_scale)
747 : CASE DEFAULT
748 6 : CPABORT("Unknown FERMI_LOW_RANK base operator during application")
749 : END SELECT
750 :
751 6 : IF (rank_used > 0) THEN
752 : CALL cp_fm_create(matrix_tmp, matrix_in%matrix_struct, name=routineN, &
753 4 : nrow=rank_used, ncol=k)
754 : CALL cp_fm_get_info(matrix_tmp, nrow_local=nrow_local, ncol_local=ncol_local, &
755 4 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
756 :
757 : CALL parallel_gemm('T', 'N', rank_used, k, n, 1.0_dp, preconditioner_env%fm, &
758 4 : matrix_in, 0.0_dp, matrix_tmp, a_first_col=nocc + 1)
759 :
760 16 : DO j = 1, ncol_local
761 25 : DO i = 1, nrow_local
762 : correction = low_rank_inverse_weight( &
763 : preconditioner_env%full_evals(nocc + row_indices(i)), &
764 : preconditioner_env%spectral_reference, &
765 : preconditioner_env%energy_gap) &
766 9 : - preconditioner_env%spectral_base_scale
767 21 : local_data(i, j) = local_data(i, j)*correction
768 : END DO
769 : END DO
770 :
771 : CALL parallel_gemm('N', 'N', n, k, rank_used, 1.0_dp, preconditioner_env%fm, &
772 4 : matrix_tmp, 1.0_dp, matrix_out, a_first_col=nocc + 1)
773 4 : CALL cp_fm_release(matrix_tmp)
774 : END IF
775 :
776 6 : CALL timestop(handle)
777 :
778 6 : END SUBROUTINE apply_low_rank_fm
779 :
780 : ! **************************************************************************************************
781 : !> \brief Apply a replicated-cell inverse operator to a distributed full matrix.
782 : !> \param preconditioner_env ...
783 : !> \param matrix_in ...
784 : !> \param matrix_out ...
785 : !> \param scale_factor ...
786 : ! **************************************************************************************************
787 156 : SUBROUTINE apply_lattice_operator_fm(preconditioner_env, matrix_in, matrix_out, scale_factor)
788 :
789 : TYPE(preconditioner_type) :: preconditioner_env
790 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
791 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: scale_factor
792 :
793 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_lattice_operator_fm'
794 :
795 : INTEGER :: fft_ok, handle, i
796 : LOGICAL :: used_fft
797 : REAL(KIND=dp) :: imaginary_residual, output_max, scale
798 156 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: base_input, coarse_projected, &
799 156 : coarse_residual, coarse_rhs, &
800 156 : local_input, local_output
801 :
802 156 : CALL timeset(routineN, handle)
803 156 : CPASSERT(ALLOCATED(preconditioner_env%lattice_inverse_k))
804 156 : scale = 1.0_dp
805 156 : IF (PRESENT(scale_factor)) scale = scale_factor
806 : CALL distribute_lattice_rhs(preconditioner_env, matrix_in, &
807 156 : preconditioner_env%lattice_rhs_distribution, local_input)
808 624 : ALLOCATE (local_output(SIZE(local_input, 1), SIZE(local_input, 2)))
809 156 : IF (ALLOCATED(preconditioner_env%lattice_local_indices)) THEN
810 0 : CPASSERT(ALLOCATED(preconditioner_env%lattice_local_operator))
811 0 : CPASSERT(ALLOCATED(preconditioner_env%lattice_local_coarse_inverse))
812 : ALLOCATE (base_input(SIZE(local_input, 1), SIZE(local_input, 2)), &
813 : coarse_projected(SIZE(preconditioner_env%lattice_local_indices), SIZE(local_input, 2)), &
814 0 : coarse_rhs(SIZE(preconditioner_env%lattice_local_indices), SIZE(local_input, 2)))
815 0 : DO i = 1, SIZE(preconditioner_env%lattice_local_indices)
816 0 : coarse_projected(i, :) = local_input(preconditioner_env%lattice_local_indices(i), :)
817 : END DO
818 0 : coarse_rhs(:, :) = MATMUL(preconditioner_env%lattice_local_coarse_inverse, &
819 0 : coarse_projected)
820 : base_input(:, :) = local_input - &
821 0 : MATMUL(preconditioner_env%lattice_local_operator, coarse_rhs)
822 : ELSE
823 13728 : ALLOCATE (base_input(SIZE(local_input, 1), SIZE(local_input, 2)), source=local_input)
824 : END IF
825 156 : used_fft = .TRUE.
826 156 : imaginary_residual = 0.0_dp
827 156 : IF (SIZE(local_input, 2) > 0) THEN
828 : CALL apply_lattice_inverse_dense(preconditioner_env%lattice_inverse_k, &
829 : preconditioner_env%lattice_dims, base_input, local_output, &
830 156 : used_fft, imaginary_residual)
831 : END IF
832 156 : IF (ALLOCATED(preconditioner_env%lattice_local_indices)) THEN
833 0 : ALLOCATE (coarse_residual(SIZE(coarse_rhs, 1), SIZE(coarse_rhs, 2)))
834 0 : coarse_projected(:, :) = MATMUL(TRANSPOSE(preconditioner_env%lattice_local_operator), &
835 0 : local_output)
836 0 : coarse_residual(:, :) = MATMUL(preconditioner_env%lattice_local_coarse_inverse, &
837 0 : coarse_projected)
838 : local_output(preconditioner_env%lattice_local_indices, :) = &
839 0 : local_output(preconditioner_env%lattice_local_indices, :) + coarse_rhs - coarse_residual
840 0 : DEALLOCATE (coarse_projected, coarse_residual, coarse_rhs)
841 : END IF
842 156 : fft_ok = MERGE(1, 0, used_fft)
843 156 : output_max = 0.0_dp
844 13728 : IF (SIZE(local_output) > 0) output_max = MAXVAL(ABS(local_output))
845 156 : CALL preconditioner_env%para_env%min(fft_ok)
846 156 : CALL preconditioner_env%para_env%max(imaginary_residual)
847 156 : CALL preconditioner_env%para_env%max(output_max)
848 156 : IF (fft_ok == 0) CPABORT("Lattice preconditioner inverse failed to use FFTs")
849 156 : IF (imaginary_residual > 1.0E-10_dp*MAX(1.0_dp, output_max)) THEN
850 0 : CPABORT("Lattice preconditioner inverse produced a non-real Gamma-point result")
851 : END IF
852 : CALL collect_lattice_rhs(preconditioner_env, preconditioner_env%lattice_rhs_distribution, &
853 156 : local_output, matrix_out, scale)
854 156 : DEALLOCATE (base_input, local_input, local_output)
855 156 : CALL timestop(handle)
856 :
857 312 : END SUBROUTINE apply_lattice_operator_fm
858 :
859 : ! **************************************************************************************************
860 : !> \brief Apply state-dependent replicated-cell inverse operators to a distributed full matrix.
861 : !> \param preconditioner_env ...
862 : !> \param matrix_in ...
863 : !> \param matrix_out ...
864 : ! **************************************************************************************************
865 0 : SUBROUTINE apply_lattice_state_operator_fm(preconditioner_env, matrix_in, matrix_out)
866 :
867 : TYPE(preconditioner_type) :: preconditioner_env
868 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
869 :
870 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_lattice_state_operator_fm'
871 :
872 : INTEGER :: fft_ok, global_column, handle, i, j, k, &
873 : n, state
874 : LOGICAL :: used_fft
875 : REAL(KIND=dp) :: imaginary_residual, output_max
876 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: base_input, coarse_projected, &
877 0 : coarse_residual, coarse_rhs, complement_input, local_input, local_output, &
878 0 : occupied_component, occupied_update
879 :
880 0 : CALL timeset(routineN, handle)
881 0 : CPASSERT(ALLOCATED(preconditioner_env%lattice_state_inverse_k))
882 0 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
883 0 : CPASSERT(SIZE(preconditioner_env%lattice_state_inverse_k, 4) == k)
884 : CALL distribute_lattice_rhs(preconditioner_env, matrix_in, &
885 0 : preconditioner_env%lattice_rhs_distribution, local_input)
886 0 : CPASSERT(ALLOCATED(preconditioner_env%lattice_occ_vectors))
887 0 : CPASSERT(ALLOCATED(preconditioner_env%lattice_occ_dual))
888 : ALLOCATE (complement_input(n, SIZE(local_input, 2)), &
889 : local_output(n, SIZE(local_input, 2)), &
890 : occupied_component(k, SIZE(local_input, 2)), &
891 0 : occupied_update(n, SIZE(local_input, 2)))
892 0 : occupied_component(:, :) = MATMUL(TRANSPOSE(preconditioner_env%lattice_occ_vectors), &
893 0 : local_input)
894 : complement_input(:, :) = local_input - &
895 0 : MATMUL(preconditioner_env%lattice_occ_dual, occupied_component)
896 0 : ALLOCATE (base_input(n, SIZE(local_input, 2)), source=complement_input)
897 0 : IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
898 0 : CPASSERT(ALLOCATED(preconditioner_env%lattice_state_local_operator))
899 0 : CPASSERT(ALLOCATED(preconditioner_env%lattice_state_local_coarse_inverse))
900 : ALLOCATE (coarse_rhs(SIZE(preconditioner_env%lattice_state_local_indices, 1), &
901 : SIZE(local_input, 2)), &
902 : coarse_projected(SIZE(preconditioner_env%lattice_state_local_indices, 1), &
903 0 : SIZE(local_input, 2)))
904 0 : DO j = 1, SIZE(local_input, 2)
905 0 : state = preconditioner_env%lattice_rhs_distribution%global_columns(j)
906 0 : DO i = 1, SIZE(preconditioner_env%lattice_state_local_indices, 1)
907 : coarse_projected(i, j) = &
908 0 : complement_input(preconditioner_env%lattice_state_local_indices(i, state), j)
909 : END DO
910 : coarse_rhs(:, j) = MATMUL( &
911 : preconditioner_env%lattice_state_local_coarse_inverse(:, :, state), &
912 0 : coarse_projected(:, j))
913 : base_input(:, j) = complement_input(:, j) - &
914 : MATMUL(preconditioner_env%lattice_state_local_operator(:, :, state), &
915 0 : coarse_rhs(:, j))
916 : END DO
917 : END IF
918 0 : used_fft = .TRUE.
919 0 : imaginary_residual = 0.0_dp
920 0 : IF (SIZE(local_input, 2) > 0) THEN
921 : CALL apply_lattice_state_inverse_dense(preconditioner_env%lattice_state_inverse_k, &
922 : preconditioner_env%lattice_dims, &
923 : base_input, local_output, used_fft, imaginary_residual, &
924 0 : preconditioner_env%lattice_rhs_distribution%global_columns)
925 : END IF
926 0 : IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
927 0 : ALLOCATE (coarse_residual(SIZE(coarse_rhs, 1), SIZE(coarse_rhs, 2)))
928 0 : DO j = 1, SIZE(local_input, 2)
929 0 : state = preconditioner_env%lattice_rhs_distribution%global_columns(j)
930 : coarse_projected(:, j) = &
931 : MATMUL(TRANSPOSE(preconditioner_env%lattice_state_local_operator(:, :, state)), &
932 0 : local_output(:, j))
933 : coarse_residual(:, j) = MATMUL( &
934 : preconditioner_env%lattice_state_local_coarse_inverse(:, :, state), &
935 0 : coarse_projected(:, j))
936 : local_output(preconditioner_env%lattice_state_local_indices(:, state), j) = &
937 : local_output(preconditioner_env%lattice_state_local_indices(:, state), j) + &
938 0 : coarse_rhs(:, j) - coarse_residual(:, j)
939 : END DO
940 0 : DEALLOCATE (coarse_projected, coarse_residual, coarse_rhs)
941 : END IF
942 0 : fft_ok = MERGE(1, 0, used_fft)
943 0 : output_max = 0.0_dp
944 0 : IF (SIZE(local_output) > 0) output_max = MAXVAL(ABS(local_output))
945 0 : CALL preconditioner_env%para_env%min(fft_ok)
946 0 : CALL preconditioner_env%para_env%max(imaginary_residual)
947 0 : CALL preconditioner_env%para_env%max(output_max)
948 0 : IF (fft_ok == 0) THEN
949 0 : CPABORT("State-dependent lattice preconditioner inverse failed to use FFTs")
950 : END IF
951 0 : IF (imaginary_residual > 1.0E-10_dp*MAX(1.0_dp, output_max)) THEN
952 0 : CPABORT("State-dependent lattice preconditioner produced a non-real Gamma-point result")
953 : END IF
954 0 : occupied_component(:, :) = MATMUL(TRANSPOSE(preconditioner_env%lattice_occ_dual), &
955 0 : local_output)
956 0 : occupied_update(:, :) = MATMUL(preconditioner_env%lattice_occ_vectors, occupied_component)
957 0 : local_output(:, :) = local_output - occupied_update
958 0 : occupied_component(:, :) = MATMUL(TRANSPOSE(preconditioner_env%lattice_occ_vectors), &
959 0 : local_input)
960 0 : DO j = 1, SIZE(local_input, 2)
961 0 : global_column = preconditioner_env%lattice_rhs_distribution%global_columns(j)
962 0 : DO i = 1, k
963 : occupied_component(i, j) = occupied_component(i, j)/MAX( &
964 : preconditioner_env%energy_gap, &
965 : preconditioner_env%full_evals(i) - &
966 0 : preconditioner_env%occ_evals(global_column))
967 : END DO
968 : END DO
969 0 : occupied_update(:, :) = MATMUL(preconditioner_env%lattice_occ_vectors, occupied_component)
970 0 : local_output(:, :) = local_output + occupied_update
971 : CALL collect_lattice_rhs(preconditioner_env, preconditioner_env%lattice_rhs_distribution, &
972 0 : local_output, matrix_out)
973 0 : DEALLOCATE (base_input, complement_input, local_input, local_output, occupied_component, occupied_update)
974 0 : CALL timestop(handle)
975 :
976 0 : END SUBROUTINE apply_lattice_state_operator_fm
977 :
978 : ! **************************************************************************************************
979 : !> \brief Route complete orbital columns to distinct ranks for lattice-FFT application.
980 : !> \param preconditioner_env ...
981 : !> \param matrix_in distributed input matrix
982 : !> \param distribution communication map retained for the reverse route
983 : !> \param local_input complete columns owned by this rank
984 : ! **************************************************************************************************
985 156 : SUBROUTINE distribute_lattice_rhs(preconditioner_env, matrix_in, distribution, local_input)
986 :
987 : TYPE(preconditioner_type) :: preconditioner_env
988 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in
989 : TYPE(lattice_rhs_distribution_type), INTENT(INOUT) :: distribution
990 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
991 : INTENT(OUT) :: local_input
992 :
993 : INTEGER :: destination, i, ientry, j, k, n, &
994 : ncol_local, nowned, nrow_local, &
995 : num_pe, slot
996 156 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cursor, send_global_columns, &
997 156 : send_global_rows
998 156 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
999 : LOGICAL :: rebuild
1000 156 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: recv_values, send_values
1001 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1002 156 : POINTER :: input_data
1003 :
1004 156 : num_pe = preconditioner_env%para_env%num_pe
1005 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k, &
1006 : nrow_local=nrow_local, ncol_local=ncol_local, &
1007 156 : row_indices=row_indices, col_indices=col_indices, local_data=input_data)
1008 : rebuild = distribution%nrow_global /= n .OR. distribution%ncol_global /= k .OR. &
1009 : distribution%nrow_local /= nrow_local .OR. distribution%ncol_local /= ncol_local .OR. &
1010 156 : distribution%num_pe /= num_pe
1011 : IF (.NOT. rebuild) THEN
1012 : rebuild = .NOT. ALLOCATED(distribution%local_rows) .OR. &
1013 154 : .NOT. ALLOCATED(distribution%local_columns)
1014 : END IF
1015 : IF (.NOT. rebuild) THEN
1016 : rebuild = ANY(distribution%local_rows /= row_indices) .OR. &
1017 2926 : ANY(distribution%local_columns /= col_indices)
1018 : END IF
1019 154 : IF (rebuild) THEN
1020 2 : CALL release_lattice_rhs_distribution(distribution)
1021 2 : distribution%nrow_global = n
1022 2 : distribution%ncol_global = k
1023 2 : distribution%nrow_local = nrow_local
1024 2 : distribution%ncol_local = ncol_local
1025 2 : distribution%num_pe = num_pe
1026 10 : ALLOCATE (distribution%local_rows(nrow_local), distribution%local_columns(ncol_local))
1027 18 : distribution%local_rows(:) = row_indices
1028 22 : distribution%local_columns(:) = col_indices
1029 0 : ALLOCATE (distribution%send_counts(num_pe), distribution%send_displacements(num_pe), &
1030 0 : distribution%recv_counts(num_pe), distribution%recv_displacements(num_pe), &
1031 14 : cursor(num_pe))
1032 6 : distribution%send_counts = 0
1033 22 : DO j = 1, ncol_local
1034 20 : destination = MODULO(col_indices(j) - 1, num_pe) + 1
1035 : distribution%send_counts(destination) = &
1036 22 : distribution%send_counts(destination) + nrow_local
1037 : END DO
1038 : CALL preconditioner_env%para_env%alltoall(distribution%send_counts, &
1039 2 : distribution%recv_counts, 1)
1040 2 : distribution%send_displacements(1) = 0
1041 2 : distribution%recv_displacements(1) = 0
1042 4 : DO i = 2, num_pe
1043 : distribution%send_displacements(i) = distribution%send_displacements(i - 1) + &
1044 2 : distribution%send_counts(i - 1)
1045 : distribution%recv_displacements(i) = distribution%recv_displacements(i - 1) + &
1046 4 : distribution%recv_counts(i - 1)
1047 : END DO
1048 :
1049 : ALLOCATE (send_global_rows(SUM(distribution%send_counts)), &
1050 : send_global_columns(SUM(distribution%send_counts)), &
1051 0 : distribution%send_local_rows(SUM(distribution%send_counts)), &
1052 0 : distribution%send_local_columns(SUM(distribution%send_counts)), &
1053 0 : distribution%recv_global_rows(SUM(distribution%recv_counts)), &
1054 50 : distribution%recv_global_columns(SUM(distribution%recv_counts)))
1055 6 : cursor(:) = distribution%send_displacements
1056 22 : DO j = 1, ncol_local
1057 20 : destination = MODULO(col_indices(j) - 1, num_pe) + 1
1058 182 : DO i = 1, nrow_local
1059 160 : cursor(destination) = cursor(destination) + 1
1060 160 : ientry = cursor(destination)
1061 160 : send_global_rows(ientry) = row_indices(i)
1062 160 : send_global_columns(ientry) = col_indices(j)
1063 160 : distribution%send_local_rows(ientry) = i
1064 180 : distribution%send_local_columns(ientry) = j
1065 : END DO
1066 : END DO
1067 : CALL preconditioner_env%para_env%alltoall( &
1068 : send_global_rows, distribution%send_counts, distribution%send_displacements, &
1069 2 : distribution%recv_global_rows, distribution%recv_counts, distribution%recv_displacements)
1070 : CALL preconditioner_env%para_env%alltoall( &
1071 : send_global_columns, distribution%send_counts, distribution%send_displacements, &
1072 2 : distribution%recv_global_columns, distribution%recv_counts, distribution%recv_displacements)
1073 :
1074 2 : nowned = 0
1075 2 : DO j = preconditioner_env%para_env%mepos + 1, k, num_pe
1076 10 : nowned = nowned + 1
1077 : END DO
1078 6 : ALLOCATE (distribution%global_columns(nowned))
1079 2 : nowned = 0
1080 2 : DO j = preconditioner_env%para_env%mepos + 1, k, num_pe
1081 10 : nowned = nowned + 1
1082 10 : distribution%global_columns(nowned) = j
1083 : END DO
1084 2 : DEALLOCATE (cursor, send_global_columns, send_global_rows)
1085 : END IF
1086 :
1087 156 : nowned = SIZE(distribution%global_columns)
1088 624 : ALLOCATE (local_input(n, nowned), source=0.0_dp)
1089 : ALLOCATE (send_values(SUM(distribution%send_counts)), &
1090 1404 : recv_values(SUM(distribution%recv_counts)))
1091 12636 : DO ientry = 1, SIZE(send_values)
1092 : send_values(ientry) = input_data(distribution%send_local_rows(ientry), &
1093 12636 : distribution%send_local_columns(ientry))
1094 : END DO
1095 : CALL preconditioner_env%para_env%alltoall( &
1096 : send_values, distribution%send_counts, distribution%send_displacements, &
1097 156 : recv_values, distribution%recv_counts, distribution%recv_displacements)
1098 12636 : DO ientry = 1, SIZE(recv_values)
1099 12480 : slot = (distribution%recv_global_columns(ientry) - 1)/num_pe + 1
1100 12636 : local_input(distribution%recv_global_rows(ientry), slot) = recv_values(ientry)
1101 : END DO
1102 :
1103 156 : DEALLOCATE (recv_values, send_values)
1104 :
1105 312 : END SUBROUTINE distribute_lattice_rhs
1106 :
1107 : ! **************************************************************************************************
1108 : !> \brief Return rank-owned lattice-FFT columns to their original full-matrix layout.
1109 : !> \param preconditioner_env ...
1110 : !> \param distribution communication map from distribute_lattice_rhs
1111 : !> \param local_output complete preconditioned columns owned by this rank
1112 : !> \param matrix_out distributed output matrix
1113 : !> \param scale_factor optional final scaling
1114 : ! **************************************************************************************************
1115 156 : SUBROUTINE collect_lattice_rhs(preconditioner_env, distribution, local_output, matrix_out, scale_factor)
1116 :
1117 : TYPE(preconditioner_type) :: preconditioner_env
1118 : TYPE(lattice_rhs_distribution_type), INTENT(IN) :: distribution
1119 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: local_output
1120 : TYPE(cp_fm_type), INTENT(IN) :: matrix_out
1121 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: scale_factor
1122 :
1123 : INTEGER :: ientry, num_pe, slot
1124 : REAL(KIND=dp) :: scale
1125 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: recv_values, send_values
1126 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1127 156 : POINTER :: output_data
1128 :
1129 156 : num_pe = preconditioner_env%para_env%num_pe
1130 0 : CPASSERT(SIZE(local_output, 1) == distribution%nrow_global)
1131 156 : CPASSERT(SIZE(local_output, 2) == SIZE(distribution%global_columns))
1132 156 : scale = 1.0_dp
1133 156 : IF (PRESENT(scale_factor)) scale = scale_factor
1134 : ALLOCATE (send_values(SUM(distribution%recv_counts)), &
1135 1404 : recv_values(SUM(distribution%send_counts)))
1136 12636 : DO ientry = 1, SIZE(send_values)
1137 12480 : slot = (distribution%recv_global_columns(ientry) - 1)/num_pe + 1
1138 12636 : send_values(ientry) = local_output(distribution%recv_global_rows(ientry), slot)
1139 : END DO
1140 : CALL preconditioner_env%para_env%alltoall( &
1141 : send_values, distribution%recv_counts, distribution%recv_displacements, &
1142 156 : recv_values, distribution%send_counts, distribution%send_displacements)
1143 :
1144 156 : CALL cp_fm_get_info(matrix_out, local_data=output_data)
1145 156 : output_data = 0.0_dp
1146 12636 : DO ientry = 1, SIZE(recv_values)
1147 : output_data(distribution%send_local_rows(ientry), &
1148 12636 : distribution%send_local_columns(ientry)) = scale*recv_values(ientry)
1149 : END DO
1150 156 : DEALLOCATE (recv_values, send_values)
1151 :
1152 156 : END SUBROUTINE collect_lattice_rhs
1153 :
1154 : ! **************************************************************************************************
1155 : !> \brief Bridge a DBCSR OT gradient to the full-matrix lattice-FFT operator.
1156 : !> \param preconditioner_env ...
1157 : !> \param matrix_in ...
1158 : !> \param matrix_out ...
1159 : ! **************************************************************************************************
1160 468 : SUBROUTINE apply_lattice_operator_dbcsr(preconditioner_env, matrix_in, matrix_out)
1161 :
1162 : TYPE(preconditioner_type) :: preconditioner_env
1163 : TYPE(dbcsr_type) :: matrix_in, matrix_out
1164 :
1165 : INTEGER :: k, n
1166 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1167 : TYPE(cp_fm_type) :: input_fm, output_fm
1168 :
1169 156 : NULLIFY (fm_struct)
1170 156 : CALL dbcsr_get_info(matrix_in, nfullrows_total=n, nfullcols_total=k)
1171 : CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=k, &
1172 : context=preconditioner_env%ctxt, &
1173 156 : para_env=preconditioner_env%para_env)
1174 156 : CALL cp_fm_create(input_fm, fm_struct, name="lattice preconditioner input")
1175 156 : CALL cp_fm_create(output_fm, fm_struct, name="lattice preconditioner output")
1176 156 : CALL cp_fm_struct_release(fm_struct)
1177 156 : CALL copy_dbcsr_to_fm(matrix_in, input_fm)
1178 156 : CALL apply_lattice_operator_fm(preconditioner_env, input_fm, output_fm)
1179 156 : CALL copy_fm_to_dbcsr(output_fm, matrix_out)
1180 156 : CALL cp_fm_release(input_fm)
1181 156 : CALL cp_fm_release(output_fm)
1182 :
1183 156 : END SUBROUTINE apply_lattice_operator_dbcsr
1184 :
1185 : ! **************************************************************************************************
1186 : !> \brief DBCSR variant of the bounded low-rank spectral correction.
1187 : !> \param preconditioner_env ...
1188 : !> \param matrix_in ...
1189 : !> \param matrix_out ...
1190 : ! **************************************************************************************************
1191 0 : SUBROUTINE apply_low_rank_dbcsr(preconditioner_env, matrix_in, matrix_out)
1192 :
1193 : TYPE(preconditioner_type) :: preconditioner_env
1194 : TYPE(dbcsr_type) :: matrix_in, matrix_out
1195 :
1196 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_low_rank_dbcsr'
1197 :
1198 : INTEGER :: col, col_offset, col_size, handle, i, j, &
1199 : k, n, nocc, rank_used, row, &
1200 : row_offset, row_size
1201 : REAL(KIND=dp) :: correction
1202 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: DATA
1203 : TYPE(dbcsr_iterator_type) :: iter
1204 : TYPE(dbcsr_type) :: matrix_tmp
1205 :
1206 0 : CALL timeset(routineN, handle)
1207 :
1208 0 : IF (preconditioner_env%low_rank_base == ot_low_rank_base_lattice_fft) THEN
1209 0 : CALL apply_low_rank_dbcsr_via_fm(preconditioner_env, matrix_in, matrix_out)
1210 0 : CALL timestop(handle)
1211 0 : RETURN
1212 : END IF
1213 :
1214 0 : CPASSERT(ASSOCIATED(preconditioner_env%base_dbcsr_matrix))
1215 0 : CPASSERT(ASSOCIATED(preconditioner_env%dbcsr_matrix))
1216 0 : CPASSERT(ASSOCIATED(preconditioner_env%full_evals))
1217 0 : CPASSERT(ASSOCIATED(preconditioner_env%occ_evals))
1218 :
1219 0 : CALL dbcsr_get_info(matrix_in, nfullrows_total=n, nfullcols_total=k)
1220 0 : nocc = SIZE(preconditioner_env%occ_evals)
1221 0 : CPASSERT(n == SIZE(preconditioner_env%full_evals))
1222 0 : rank_used = MIN(preconditioner_env%spectral_rank, n - nocc)
1223 :
1224 : ! DBCSR form of the same overlap-inverse base plus positive low-rank update.
1225 : CALL dbcsr_multiply('N', 'N', preconditioner_env%spectral_base_scale, &
1226 0 : preconditioner_env%base_dbcsr_matrix, matrix_in, 0.0_dp, matrix_out)
1227 :
1228 0 : IF (rank_used > 0) THEN
1229 0 : CALL dbcsr_copy(matrix_tmp, matrix_in, name=routineN)
1230 0 : CALL dbcsr_set(matrix_tmp, 0.0_dp)
1231 : CALL dbcsr_multiply('T', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
1232 : matrix_in, 0.0_dp, matrix_tmp, &
1233 0 : first_row=nocc + 1, last_row=nocc + rank_used)
1234 :
1235 0 : CALL dbcsr_iterator_start(iter, matrix_tmp)
1236 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1237 : CALL dbcsr_iterator_next_block(iter, row, col, DATA, &
1238 : row_size=row_size, col_size=col_size, &
1239 0 : row_offset=row_offset, col_offset=col_offset)
1240 0 : DO j = 1, col_size
1241 0 : DO i = 1, row_size
1242 0 : IF (row_offset + i - 1 >= nocc + 1 .AND. &
1243 0 : row_offset + i - 1 <= nocc + rank_used) THEN
1244 : correction = low_rank_inverse_weight( &
1245 : preconditioner_env%full_evals(row_offset + i - 1), &
1246 : preconditioner_env%spectral_reference, &
1247 : preconditioner_env%energy_gap) &
1248 0 : - preconditioner_env%spectral_base_scale
1249 0 : DATA(i, j) = DATA(i, j)*correction
1250 : END IF
1251 : END DO
1252 : END DO
1253 : END DO
1254 0 : CALL dbcsr_iterator_stop(iter)
1255 :
1256 : CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
1257 : matrix_tmp, 1.0_dp, matrix_out, &
1258 0 : first_k=nocc + 1, last_k=nocc + rank_used)
1259 0 : CALL dbcsr_release(matrix_tmp)
1260 : END IF
1261 :
1262 0 : CALL timestop(handle)
1263 :
1264 0 : END SUBROUTINE apply_low_rank_dbcsr
1265 :
1266 : ! **************************************************************************************************
1267 : !> \brief Bridge DBCSR OT gradients to the full-matrix lattice-FFT implementation.
1268 : !> \param preconditioner_env ...
1269 : !> \param matrix_in ...
1270 : !> \param matrix_out ...
1271 : ! **************************************************************************************************
1272 0 : SUBROUTINE apply_low_rank_dbcsr_via_fm(preconditioner_env, matrix_in, matrix_out)
1273 :
1274 : TYPE(preconditioner_type) :: preconditioner_env
1275 : TYPE(dbcsr_type) :: matrix_in, matrix_out
1276 :
1277 : INTEGER :: k, n
1278 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1279 : TYPE(cp_fm_type) :: input_fm, output_fm
1280 :
1281 0 : NULLIFY (fm_struct)
1282 0 : CALL dbcsr_get_info(matrix_in, nfullrows_total=n, nfullcols_total=k)
1283 : CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=k, &
1284 : context=preconditioner_env%ctxt, &
1285 0 : para_env=preconditioner_env%para_env)
1286 0 : CALL cp_fm_create(input_fm, fm_struct, name="lattice low-rank input")
1287 0 : CALL cp_fm_create(output_fm, fm_struct, name="lattice low-rank output")
1288 0 : CALL cp_fm_struct_release(fm_struct)
1289 0 : CALL copy_dbcsr_to_fm(matrix_in, input_fm)
1290 0 : CALL apply_low_rank_fm(preconditioner_env, input_fm, output_fm)
1291 0 : CALL copy_fm_to_dbcsr(output_fm, matrix_out)
1292 0 : CALL cp_fm_release(input_fm)
1293 0 : CALL cp_fm_release(output_fm)
1294 :
1295 0 : END SUBROUTINE apply_low_rank_dbcsr_via_fm
1296 :
1297 : ! **************************************************************************************************
1298 : !> \brief full all to a full matrix
1299 : !> \param preconditioner_env ...
1300 : !> \param matrix_in ...
1301 : !> \param matrix_out ...
1302 : ! **************************************************************************************************
1303 109864 : SUBROUTINE apply_full_all(preconditioner_env, matrix_in, matrix_out)
1304 :
1305 : TYPE(preconditioner_type) :: preconditioner_env
1306 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
1307 :
1308 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_all'
1309 :
1310 : INTEGER :: handle, i, j, k, n, ncol_local, &
1311 : nrow_local
1312 27466 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1313 : REAL(KIND=dp) :: dum
1314 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1315 27466 : POINTER :: local_data
1316 : TYPE(cp_fm_type) :: matrix_tmp
1317 :
1318 27466 : CALL timeset(routineN, handle)
1319 :
1320 27466 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
1321 :
1322 27466 : CALL cp_fm_create(matrix_tmp, matrix_in%matrix_struct, name="apply_full_all")
1323 : CALL cp_fm_get_info(matrix_tmp, nrow_local=nrow_local, ncol_local=ncol_local, &
1324 27466 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
1325 :
1326 : !
1327 : CALL parallel_gemm('T', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
1328 27466 : matrix_in, 0.0_dp, matrix_tmp)
1329 :
1330 : ! do the right scaling
1331 218394 : DO j = 1, ncol_local
1332 2750195 : DO i = 1, nrow_local
1333 : dum = 1.0_dp/MAX(preconditioner_env%energy_gap, &
1334 2531801 : preconditioner_env%full_evals(row_indices(i)) - preconditioner_env%occ_evals(col_indices(j)))
1335 2722729 : local_data(i, j) = local_data(i, j)*dum
1336 : END DO
1337 : END DO
1338 :
1339 : ! mult back
1340 : CALL parallel_gemm('N', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
1341 27466 : matrix_tmp, 0.0_dp, matrix_out)
1342 :
1343 27466 : CALL cp_fm_release(matrix_tmp)
1344 :
1345 27466 : CALL timestop(handle)
1346 :
1347 27466 : END SUBROUTINE apply_full_all
1348 :
1349 : ! **************************************************************************************************
1350 : !> \brief full all to a dbcsr matrix
1351 : !> \param preconditioner_env ...
1352 : !> \param matrix_in ...
1353 : !> \param matrix_out ...
1354 : ! **************************************************************************************************
1355 50752 : SUBROUTINE apply_all(preconditioner_env, matrix_in, matrix_out)
1356 :
1357 : TYPE(preconditioner_type) :: preconditioner_env
1358 : TYPE(dbcsr_type) :: matrix_in, matrix_out
1359 :
1360 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_all'
1361 :
1362 : INTEGER :: col, col_offset, col_size, handle, i, j, &
1363 : row, row_offset, row_size
1364 : REAL(KIND=dp) :: dum
1365 25376 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: DATA
1366 : TYPE(dbcsr_iterator_type) :: iter
1367 : TYPE(dbcsr_type) :: matrix_tmp
1368 :
1369 25376 : CALL timeset(routineN, handle)
1370 :
1371 25376 : CALL dbcsr_copy(matrix_tmp, matrix_in, name="apply_full_all")
1372 : CALL dbcsr_multiply('T', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
1373 25376 : matrix_in, 0.0_dp, matrix_tmp)
1374 : ! do the right scaling
1375 25376 : CALL dbcsr_iterator_start(iter, matrix_tmp)
1376 72572 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1377 : CALL dbcsr_iterator_next_block(iter, row, col, DATA, &
1378 : row_size=row_size, col_size=col_size, &
1379 47196 : row_offset=row_offset, col_offset=col_offset)
1380 487270 : DO j = 1, col_size
1381 3712637 : DO i = 1, row_size
1382 : dum = 1.0_dp/MAX(preconditioner_env%energy_gap, &
1383 : preconditioner_env%full_evals(row_offset + i - 1) &
1384 3250743 : - preconditioner_env%occ_evals(col_offset + j - 1))
1385 3665441 : DATA(i, j) = DATA(i, j)*dum
1386 : END DO
1387 : END DO
1388 : END DO
1389 25376 : CALL dbcsr_iterator_stop(iter)
1390 :
1391 : ! mult back
1392 : CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
1393 25376 : matrix_tmp, 0.0_dp, matrix_out)
1394 25376 : CALL dbcsr_release(matrix_tmp)
1395 25376 : CALL timestop(handle)
1396 :
1397 25376 : END SUBROUTINE apply_all
1398 :
1399 : ! **************************************************************************************************
1400 : !> \brief Apply the state-selective inverse in a private canonical occupied gauge.
1401 : !> \param preconditioner_env ...
1402 : !> \param matrix_in ...
1403 : !> \param matrix_out ...
1404 : ! **************************************************************************************************
1405 32 : SUBROUTINE apply_full_all_covariant(preconditioner_env, matrix_in, matrix_out)
1406 :
1407 : TYPE(preconditioner_type) :: preconditioner_env
1408 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
1409 :
1410 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_all_covariant'
1411 :
1412 : INTEGER :: handle, k, n
1413 : TYPE(cp_fm_type) :: canonical_in, canonical_out
1414 :
1415 8 : CALL timeset(routineN, handle)
1416 8 : CPASSERT(ASSOCIATED(preconditioner_env%occ_rotation))
1417 8 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
1418 :
1419 8 : CALL cp_fm_create(canonical_in, matrix_in%matrix_struct, name=routineN)
1420 8 : CALL cp_fm_create(canonical_out, matrix_out%matrix_struct, name=routineN)
1421 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, matrix_in, &
1422 8 : preconditioner_env%occ_rotation, 0.0_dp, canonical_in)
1423 8 : IF (preconditioner_env%lattice_fft_active) THEN
1424 0 : CALL apply_lattice_state_operator_fm(preconditioner_env, canonical_in, canonical_out)
1425 : ELSE
1426 8 : CALL apply_full_all(preconditioner_env, canonical_in, canonical_out)
1427 : END IF
1428 : CALL parallel_gemm('N', 'T', n, k, k, 1.0_dp, canonical_out, &
1429 8 : preconditioner_env%occ_rotation, 0.0_dp, matrix_out)
1430 8 : CALL cp_fm_release(canonical_in)
1431 8 : CALL cp_fm_release(canonical_out)
1432 8 : CALL timestop(handle)
1433 :
1434 8 : END SUBROUTINE apply_full_all_covariant
1435 :
1436 : ! **************************************************************************************************
1437 : !> \brief DBCSR form of the rotation-covariant state-selective inverse.
1438 : !> \param preconditioner_env ...
1439 : !> \param matrix_in ...
1440 : !> \param matrix_out ...
1441 : ! **************************************************************************************************
1442 0 : SUBROUTINE apply_all_covariant(preconditioner_env, matrix_in, matrix_out)
1443 :
1444 : TYPE(preconditioner_type) :: preconditioner_env
1445 : TYPE(dbcsr_type) :: matrix_in, matrix_out
1446 :
1447 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_all_covariant'
1448 :
1449 : INTEGER :: handle, k, n
1450 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1451 : TYPE(cp_fm_type) :: input_fm, output_fm
1452 : TYPE(dbcsr_type) :: canonical_in, canonical_out
1453 :
1454 0 : CALL timeset(routineN, handle)
1455 0 : CPASSERT(ASSOCIATED(preconditioner_env%occ_rotation))
1456 0 : CALL dbcsr_get_info(matrix_in, nfullrows_total=n, nfullcols_total=k)
1457 :
1458 0 : IF (preconditioner_env%lattice_fft_active) THEN
1459 0 : NULLIFY (fm_struct)
1460 : CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=k, &
1461 : context=preconditioner_env%ctxt, &
1462 0 : para_env=preconditioner_env%para_env)
1463 0 : CALL cp_fm_create(input_fm, fm_struct, name=routineN)
1464 0 : CALL cp_fm_create(output_fm, fm_struct, name=routineN)
1465 0 : CALL cp_fm_struct_release(fm_struct)
1466 0 : CALL copy_dbcsr_to_fm(matrix_in, input_fm)
1467 0 : CALL apply_full_all_covariant(preconditioner_env, input_fm, output_fm)
1468 0 : CALL copy_fm_to_dbcsr(output_fm, matrix_out)
1469 0 : CALL cp_fm_release(input_fm)
1470 0 : CALL cp_fm_release(output_fm)
1471 0 : CALL timestop(handle)
1472 0 : RETURN
1473 : END IF
1474 :
1475 0 : IF (.NOT. ASSOCIATED(preconditioner_env%occ_rotation_dbcsr)) THEN
1476 0 : CALL dbcsr_init_p(preconditioner_env%occ_rotation_dbcsr)
1477 : CALL cp_dbcsr_m_by_n_from_template(preconditioner_env%occ_rotation_dbcsr, &
1478 : template=matrix_in, m=k, n=k, &
1479 0 : sym=dbcsr_type_no_symmetry)
1480 : CALL copy_fm_to_dbcsr(preconditioner_env%occ_rotation, &
1481 0 : preconditioner_env%occ_rotation_dbcsr)
1482 : END IF
1483 :
1484 0 : CALL dbcsr_copy(canonical_in, matrix_in, name=routineN)
1485 0 : CALL dbcsr_set(canonical_in, 0.0_dp)
1486 0 : CALL dbcsr_copy(canonical_out, matrix_out, name=routineN)
1487 0 : CALL dbcsr_set(canonical_out, 0.0_dp)
1488 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_in, &
1489 0 : preconditioner_env%occ_rotation_dbcsr, 0.0_dp, canonical_in)
1490 0 : CALL apply_all(preconditioner_env, canonical_in, canonical_out)
1491 : CALL dbcsr_multiply('N', 'T', 1.0_dp, canonical_out, &
1492 0 : preconditioner_env%occ_rotation_dbcsr, 0.0_dp, matrix_out)
1493 0 : CALL dbcsr_release(canonical_in)
1494 0 : CALL dbcsr_release(canonical_out)
1495 0 : CALL timestop(handle)
1496 :
1497 0 : END SUBROUTINE apply_all_covariant
1498 :
1499 0 : END MODULE preconditioner_apply
|
