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