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 48970 : 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 48970 : CALL timeset(routineN, handle)
85 :
86 48970 : SELECT CASE (preconditioner_env%in_use)
87 : CASE (0)
88 0 : CPABORT("No preconditioner in use")
89 : CASE (ot_precond_full_single)
90 1352 : CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
91 : CASE (ot_precond_full_all)
92 27458 : 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 20148 : IF (preconditioner_env%lattice_fft_active) THEN
99 0 : CALL apply_lattice_operator_fm(preconditioner_env, matrix_in, matrix_out)
100 : ELSE
101 40296 : SELECT CASE (preconditioner_env%solver)
102 : CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
103 20148 : 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 20148 : CPABORT("Solver not implemented")
110 : END SELECT
111 : END IF
112 : CASE DEFAULT
113 48970 : CPABORT("Unknown preconditioner")
114 : END SELECT
115 :
116 48970 : CALL timestop(handle)
117 :
118 48970 : END SUBROUTINE apply_preconditioner_fm
119 :
120 : ! **************************************************************************************************
121 : !> \brief ...
122 : !> \param preconditioner_env ...
123 : !> \param matrix_in ...
124 : !> \param matrix_out ...
125 : ! **************************************************************************************************
126 80061 : 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 80061 : CALL timeset(routineN, handle)
136 :
137 80061 : 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 24958 : CALL apply_all(preconditioner_env, matrix_in, matrix_out)
144 : CASE (ot_precond_fermi_low_rank)
145 0 : CALL apply_low_rank_dbcsr(preconditioner_env, matrix_in, matrix_out)
146 : CASE (ot_precond_full_all_covariant)
147 0 : 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 54901 : IF (preconditioner_env%lattice_fft_active) THEN
150 156 : CALL apply_lattice_operator_dbcsr(preconditioner_env, matrix_in, matrix_out)
151 : ELSE
152 109490 : SELECT CASE (preconditioner_env%solver)
153 : CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
154 54745 : 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 0 : CALL apply_chebyshev_dbcsr(preconditioner_env, matrix_in, matrix_out)
160 : CASE DEFAULT
161 54745 : CPABORT("Wrong solver")
162 : END SELECT
163 : END IF
164 : CASE DEFAULT
165 80061 : CPABORT("Wrong preconditioner")
166 : END SELECT
167 :
168 80061 : CALL timestop(handle)
169 :
170 80061 : 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 53070 : 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 17690 : CALL timeset(routineN, handle)
190 :
191 17690 : 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 17690 : CPABORT("Unsupported complex K-point OT preconditioner")
197 : END SELECT
198 17690 : CPASSERT(ASSOCIATED(preconditioner_env%complex_fm))
199 17690 : CALL cp_cfm_get_info(matrix_in, nrow_global=n, ncol_global=k)
200 17690 : CALL cp_cfm_get_info(preconditioner_env%complex_fm, nrow_global=npre)
201 17690 : CPASSERT(n == npre)
202 :
203 29708 : SELECT CASE (preconditioner_env%in_use)
204 : CASE (ot_precond_full_all)
205 12018 : 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 17884 : matrix_in, z_zero, matrix_out)
224 : END SELECT
225 :
226 17690 : CALL timestop(handle)
227 :
228 17690 : 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 48848 : 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 12212 : POINTER :: local_data
244 : INTEGER :: i, j, k, n, ncol_local, nrow_local
245 12212 : 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 12212 : CPASSERT(ASSOCIATED(preconditioner_env%occ_evals))
251 12212 : CALL cp_cfm_get_info(matrix_in, nrow_global=n, ncol_global=k)
252 12212 : CPASSERT(n == SIZE(preconditioner_env%full_evals))
253 12212 : CPASSERT(k == SIZE(preconditioner_env%occ_evals))
254 : CALL cp_cfm_create(matrix_spectral, matrix_in%matrix_struct, &
255 12212 : name='complex FULL_ALL spectral input')
256 : CALL cp_cfm_gemm('C', 'N', n, k, n, z_one, preconditioner_env%complex_fm, &
257 12212 : matrix_in, z_zero, matrix_spectral)
258 : CALL cp_cfm_get_info(matrix_spectral, nrow_local=nrow_local, ncol_local=ncol_local, &
259 12212 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
260 223149 : DO j = 1, ncol_local
261 7135564 : DO i = 1, nrow_local
262 : scale = 1.0_dp/MAX(preconditioner_env%energy_gap, &
263 : preconditioner_env%full_evals(row_indices(i)) - &
264 6912415 : preconditioner_env%occ_evals(col_indices(j)))
265 7123352 : 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 12212 : matrix_spectral, z_zero, matrix_out)
270 12212 : CALL cp_cfm_release(matrix_spectral)
271 :
272 12212 : 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 41656 : 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 5207 : CALL timeset(routineN, handle)
348 :
349 5207 : CALL dbcsr_get_info(matrix_in_re, nfullrows_total=n, nfullcols_total=k)
350 5207 : CPASSERT(ASSOCIATED(preconditioner_env%complex_fm))
351 :
352 : CALL cp_fm_create(matrix_in_re_fm, preconditioner_env%complex_fm%matrix_struct, &
353 5207 : 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 5207 : 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 5207 : 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 5207 : nrow=n, ncol=k, name='complex preconditioner output imaginary')
360 5207 : CALL copy_dbcsr_to_fm(matrix_in_re, matrix_in_re_fm)
361 5207 : 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 5207 : name='complex preconditioner input')
365 : CALL cp_cfm_create(matrix_out, matrix_in_re_fm%matrix_struct, &
366 5207 : name='complex preconditioner output')
367 5207 : CALL cp_fm_to_cfm(matrix_in_re_fm, matrix_in_im_fm, matrix_in)
368 5207 : CALL apply_preconditioner_cfm_complex(preconditioner_env, matrix_in, matrix_out)
369 5207 : CALL cp_cfm_to_fm(matrix_out, matrix_out_re_fm, matrix_out_im_fm)
370 5207 : CALL dbcsr_set(matrix_out_re, 0.0_dp)
371 5207 : CALL dbcsr_set(matrix_out_im, 0.0_dp)
372 5207 : CALL copy_fm_to_dbcsr(matrix_out_re_fm, matrix_out_re)
373 5207 : CALL copy_fm_to_dbcsr(matrix_out_im_fm, matrix_out_im)
374 :
375 5207 : CALL cp_cfm_release(matrix_out)
376 5207 : CALL cp_cfm_release(matrix_in)
377 5207 : CALL cp_fm_release(matrix_out_im_fm)
378 5207 : CALL cp_fm_release(matrix_out_re_fm)
379 5207 : CALL cp_fm_release(matrix_in_im_fm)
380 5207 : CALL cp_fm_release(matrix_in_re_fm)
381 :
382 5207 : CALL timestop(handle)
383 :
384 5207 : 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 43000 : 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 21500 : CALL timeset(routineN, handle)
402 :
403 21500 : 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 21500 : matrix_in, 0.0_dp, matrix_out)
406 21500 : CALL timestop(handle)
407 :
408 21500 : 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 54947 : 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 54947 : CALL timeset(routineN, handle)
426 :
427 54947 : 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 54947 : 0.0_dp, matrix_out)
432 :
433 54947 : CALL timestop(handle)
434 :
435 54947 : 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 0 : 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 0 : degree = preconditioner_env%polynomial_degree
528 0 : CPASSERT(degree >= 1)
529 0 : theta = 0.5_dp*(preconditioner_env%polynomial_max + preconditioner_env%polynomial_min)
530 0 : delta = 0.5_dp*(preconditioner_env%polynomial_max - preconditioner_env%polynomial_min)
531 0 : sigma = theta/delta
532 0 : rho_previous = 1.0_dp/sigma
533 0 : CALL dbcsr_copy(direction, matrix_in, name='Chebyshev direction')
534 0 : CALL dbcsr_copy(residual, matrix_in, name='Chebyshev residual')
535 0 : CALL dbcsr_set(matrix_out, 0.0_dp)
536 0 : CALL dbcsr_add(matrix_out, direction, 1.0_dp, 1.0_dp/theta)
537 0 : CALL dbcsr_set(direction, 0.0_dp)
538 0 : CALL dbcsr_add(direction, matrix_in, 1.0_dp, 1.0_dp/theta)
539 0 : DO iteration = 2, degree
540 0 : CALL dbcsr_set(residual, 0.0_dp)
541 0 : 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 0 : matrix_out, 1.0_dp, residual)
544 0 : rho = 1.0_dp/(2.0_dp*sigma - rho_previous)
545 0 : CALL dbcsr_add(direction, residual, rho*rho_previous, 2.0_dp*rho/delta)
546 0 : CALL dbcsr_add(matrix_out, direction, 1.0_dp, 1.0_dp)
547 0 : rho_previous = rho
548 : END DO
549 0 : CALL dbcsr_release(residual)
550 0 : CALL dbcsr_release(direction)
551 :
552 0 : 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 0 : 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 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
570 : REAL(KIND=dp) :: correction
571 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
572 0 : POINTER :: local_data
573 : TYPE(cp_fm_type) :: matrix_tmp
574 :
575 0 : CALL timeset(routineN, handle)
576 :
577 0 : CPASSERT(ASSOCIATED(preconditioner_env%fm))
578 0 : CPASSERT(ASSOCIATED(preconditioner_env%full_evals))
579 0 : CPASSERT(ASSOCIATED(preconditioner_env%occ_evals))
580 :
581 0 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
582 0 : nocc = SIZE(preconditioner_env%occ_evals)
583 0 : CPASSERT(n == SIZE(preconditioner_env%full_evals))
584 0 : 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 0 : 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 0 : preconditioner_env%spectral_base_scale)
599 : CASE DEFAULT
600 0 : CPABORT("Unknown FERMI_LOW_RANK base operator during application")
601 : END SELECT
602 :
603 0 : IF (rank_used > 0) THEN
604 : CALL cp_fm_create(matrix_tmp, matrix_in%matrix_struct, name=routineN, &
605 0 : nrow=rank_used, ncol=k)
606 : CALL cp_fm_get_info(matrix_tmp, nrow_local=nrow_local, ncol_local=ncol_local, &
607 0 : 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 0 : matrix_in, 0.0_dp, matrix_tmp, a_first_col=nocc + 1)
611 :
612 0 : DO j = 1, ncol_local
613 0 : 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 0 : - preconditioner_env%spectral_base_scale
619 0 : 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 0 : matrix_tmp, 1.0_dp, matrix_out, a_first_col=nocc + 1)
625 0 : CALL cp_fm_release(matrix_tmp)
626 : END IF
627 :
628 0 : CALL timestop(handle)
629 :
630 0 : 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 156 : 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 156 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: base_input, coarse_projected, &
651 156 : coarse_residual, coarse_rhs, &
652 156 : local_input, local_output
653 :
654 156 : CALL timeset(routineN, handle)
655 156 : CPASSERT(ALLOCATED(preconditioner_env%lattice_inverse_k))
656 156 : scale = 1.0_dp
657 156 : IF (PRESENT(scale_factor)) scale = scale_factor
658 : CALL distribute_lattice_rhs(preconditioner_env, matrix_in, &
659 156 : preconditioner_env%lattice_rhs_distribution, local_input)
660 624 : ALLOCATE (local_output(SIZE(local_input, 1), SIZE(local_input, 2)))
661 156 : IF (ALLOCATED(preconditioner_env%lattice_local_indices)) THEN
662 0 : CPASSERT(ALLOCATED(preconditioner_env%lattice_local_operator))
663 0 : 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 0 : coarse_rhs(SIZE(preconditioner_env%lattice_local_indices), SIZE(local_input, 2)))
667 0 : DO i = 1, SIZE(preconditioner_env%lattice_local_indices)
668 0 : coarse_projected(i, :) = local_input(preconditioner_env%lattice_local_indices(i), :)
669 : END DO
670 0 : coarse_rhs(:, :) = MATMUL(preconditioner_env%lattice_local_coarse_inverse, &
671 0 : coarse_projected)
672 : base_input(:, :) = local_input - &
673 0 : MATMUL(preconditioner_env%lattice_local_operator, coarse_rhs)
674 : ELSE
675 13728 : ALLOCATE (base_input(SIZE(local_input, 1), SIZE(local_input, 2)), source=local_input)
676 : END IF
677 156 : used_fft = .TRUE.
678 156 : imaginary_residual = 0.0_dp
679 156 : 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 156 : used_fft, imaginary_residual)
683 : END IF
684 156 : IF (ALLOCATED(preconditioner_env%lattice_local_indices)) THEN
685 0 : ALLOCATE (coarse_residual(SIZE(coarse_rhs, 1), SIZE(coarse_rhs, 2)))
686 0 : coarse_projected(:, :) = MATMUL(TRANSPOSE(preconditioner_env%lattice_local_operator), &
687 0 : local_output)
688 0 : coarse_residual(:, :) = MATMUL(preconditioner_env%lattice_local_coarse_inverse, &
689 0 : coarse_projected)
690 : local_output(preconditioner_env%lattice_local_indices, :) = &
691 0 : local_output(preconditioner_env%lattice_local_indices, :) + coarse_rhs - coarse_residual
692 0 : DEALLOCATE (coarse_projected, coarse_residual, coarse_rhs)
693 : END IF
694 156 : fft_ok = MERGE(1, 0, used_fft)
695 156 : output_max = 0.0_dp
696 13728 : IF (SIZE(local_output) > 0) output_max = MAXVAL(ABS(local_output))
697 156 : CALL preconditioner_env%para_env%min(fft_ok)
698 156 : CALL preconditioner_env%para_env%max(imaginary_residual)
699 156 : CALL preconditioner_env%para_env%max(output_max)
700 156 : IF (fft_ok == 0) CPABORT("Lattice preconditioner inverse failed to use FFTs")
701 156 : 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 156 : local_output, matrix_out, scale)
706 156 : DEALLOCATE (base_input, local_input, local_output)
707 156 : CALL timestop(handle)
708 :
709 312 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: base_input, coarse_projected, &
729 0 : coarse_residual, coarse_rhs, complement_input, local_input, local_output, &
730 0 : occupied_component, occupied_update
731 :
732 0 : CALL timeset(routineN, handle)
733 0 : CPASSERT(ALLOCATED(preconditioner_env%lattice_state_inverse_k))
734 0 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
735 0 : CPASSERT(SIZE(preconditioner_env%lattice_state_inverse_k, 4) == k)
736 : CALL distribute_lattice_rhs(preconditioner_env, matrix_in, &
737 0 : preconditioner_env%lattice_rhs_distribution, local_input)
738 0 : CPASSERT(ALLOCATED(preconditioner_env%lattice_occ_vectors))
739 0 : 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 0 : occupied_update(n, SIZE(local_input, 2)))
744 0 : occupied_component(:, :) = MATMUL(TRANSPOSE(preconditioner_env%lattice_occ_vectors), &
745 0 : local_input)
746 : complement_input(:, :) = local_input - &
747 0 : MATMUL(preconditioner_env%lattice_occ_dual, occupied_component)
748 0 : ALLOCATE (base_input(n, SIZE(local_input, 2)), source=complement_input)
749 0 : IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
750 0 : CPASSERT(ALLOCATED(preconditioner_env%lattice_state_local_operator))
751 0 : 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 0 : SIZE(local_input, 2)))
756 0 : DO j = 1, SIZE(local_input, 2)
757 0 : state = preconditioner_env%lattice_rhs_distribution%global_columns(j)
758 0 : DO i = 1, SIZE(preconditioner_env%lattice_state_local_indices, 1)
759 : coarse_projected(i, j) = &
760 0 : 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 0 : coarse_projected(:, j))
765 : base_input(:, j) = complement_input(:, j) - &
766 : MATMUL(preconditioner_env%lattice_state_local_operator(:, :, state), &
767 0 : coarse_rhs(:, j))
768 : END DO
769 : END IF
770 0 : used_fft = .TRUE.
771 0 : imaginary_residual = 0.0_dp
772 0 : 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 0 : preconditioner_env%lattice_rhs_distribution%global_columns)
777 : END IF
778 0 : IF (ALLOCATED(preconditioner_env%lattice_state_local_indices)) THEN
779 0 : ALLOCATE (coarse_residual(SIZE(coarse_rhs, 1), SIZE(coarse_rhs, 2)))
780 0 : DO j = 1, SIZE(local_input, 2)
781 0 : state = preconditioner_env%lattice_rhs_distribution%global_columns(j)
782 : coarse_projected(:, j) = &
783 : MATMUL(TRANSPOSE(preconditioner_env%lattice_state_local_operator(:, :, state)), &
784 0 : local_output(:, j))
785 : coarse_residual(:, j) = MATMUL( &
786 : preconditioner_env%lattice_state_local_coarse_inverse(:, :, state), &
787 0 : 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 0 : coarse_rhs(:, j) - coarse_residual(:, j)
791 : END DO
792 0 : DEALLOCATE (coarse_projected, coarse_residual, coarse_rhs)
793 : END IF
794 0 : fft_ok = MERGE(1, 0, used_fft)
795 0 : output_max = 0.0_dp
796 0 : IF (SIZE(local_output) > 0) output_max = MAXVAL(ABS(local_output))
797 0 : CALL preconditioner_env%para_env%min(fft_ok)
798 0 : CALL preconditioner_env%para_env%max(imaginary_residual)
799 0 : CALL preconditioner_env%para_env%max(output_max)
800 0 : IF (fft_ok == 0) THEN
801 0 : CPABORT("State-dependent lattice preconditioner inverse failed to use FFTs")
802 : END IF
803 0 : 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 0 : occupied_component(:, :) = MATMUL(TRANSPOSE(preconditioner_env%lattice_occ_dual), &
807 0 : local_output)
808 0 : occupied_update(:, :) = MATMUL(preconditioner_env%lattice_occ_vectors, occupied_component)
809 0 : local_output(:, :) = local_output - occupied_update
810 0 : occupied_component(:, :) = MATMUL(TRANSPOSE(preconditioner_env%lattice_occ_vectors), &
811 0 : local_input)
812 0 : DO j = 1, SIZE(local_input, 2)
813 0 : global_column = preconditioner_env%lattice_rhs_distribution%global_columns(j)
814 0 : 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 0 : preconditioner_env%occ_evals(global_column))
819 : END DO
820 : END DO
821 0 : occupied_update(:, :) = MATMUL(preconditioner_env%lattice_occ_vectors, occupied_component)
822 0 : local_output(:, :) = local_output + occupied_update
823 : CALL collect_lattice_rhs(preconditioner_env, preconditioner_env%lattice_rhs_distribution, &
824 0 : local_output, matrix_out)
825 0 : DEALLOCATE (base_input, complement_input, local_input, local_output, occupied_component, occupied_update)
826 0 : CALL timestop(handle)
827 :
828 0 : 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 156 : 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 156 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cursor, send_global_columns, &
849 156 : send_global_rows
850 156 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
851 : LOGICAL :: rebuild
852 156 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: recv_values, send_values
853 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
854 156 : POINTER :: input_data
855 :
856 156 : 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 156 : 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 156 : distribution%num_pe /= num_pe
863 : IF (.NOT. rebuild) THEN
864 : rebuild = .NOT. ALLOCATED(distribution%local_rows) .OR. &
865 154 : .NOT. ALLOCATED(distribution%local_columns)
866 : END IF
867 : IF (.NOT. rebuild) THEN
868 : rebuild = ANY(distribution%local_rows /= row_indices) .OR. &
869 2926 : ANY(distribution%local_columns /= col_indices)
870 : END IF
871 154 : IF (rebuild) THEN
872 2 : CALL release_lattice_rhs_distribution(distribution)
873 2 : distribution%nrow_global = n
874 2 : distribution%ncol_global = k
875 2 : distribution%nrow_local = nrow_local
876 2 : distribution%ncol_local = ncol_local
877 2 : distribution%num_pe = num_pe
878 10 : ALLOCATE (distribution%local_rows(nrow_local), distribution%local_columns(ncol_local))
879 18 : distribution%local_rows(:) = row_indices
880 22 : 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 14 : cursor(num_pe))
884 6 : distribution%send_counts = 0
885 22 : DO j = 1, ncol_local
886 20 : destination = MODULO(col_indices(j) - 1, num_pe) + 1
887 : distribution%send_counts(destination) = &
888 22 : distribution%send_counts(destination) + nrow_local
889 : END DO
890 : CALL preconditioner_env%para_env%alltoall(distribution%send_counts, &
891 2 : distribution%recv_counts, 1)
892 2 : distribution%send_displacements(1) = 0
893 2 : distribution%recv_displacements(1) = 0
894 4 : DO i = 2, num_pe
895 : distribution%send_displacements(i) = distribution%send_displacements(i - 1) + &
896 2 : distribution%send_counts(i - 1)
897 : distribution%recv_displacements(i) = distribution%recv_displacements(i - 1) + &
898 4 : 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 50 : distribution%recv_global_columns(SUM(distribution%recv_counts)))
907 6 : cursor(:) = distribution%send_displacements
908 22 : DO j = 1, ncol_local
909 20 : destination = MODULO(col_indices(j) - 1, num_pe) + 1
910 182 : DO i = 1, nrow_local
911 160 : cursor(destination) = cursor(destination) + 1
912 160 : ientry = cursor(destination)
913 160 : send_global_rows(ientry) = row_indices(i)
914 160 : send_global_columns(ientry) = col_indices(j)
915 160 : distribution%send_local_rows(ientry) = i
916 180 : 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 2 : 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 2 : distribution%recv_global_columns, distribution%recv_counts, distribution%recv_displacements)
925 :
926 2 : nowned = 0
927 2 : DO j = preconditioner_env%para_env%mepos + 1, k, num_pe
928 10 : nowned = nowned + 1
929 : END DO
930 6 : ALLOCATE (distribution%global_columns(nowned))
931 2 : nowned = 0
932 2 : DO j = preconditioner_env%para_env%mepos + 1, k, num_pe
933 10 : nowned = nowned + 1
934 10 : distribution%global_columns(nowned) = j
935 : END DO
936 2 : DEALLOCATE (cursor, send_global_columns, send_global_rows)
937 : END IF
938 :
939 156 : nowned = SIZE(distribution%global_columns)
940 624 : ALLOCATE (local_input(n, nowned), source=0.0_dp)
941 : ALLOCATE (send_values(SUM(distribution%send_counts)), &
942 1404 : recv_values(SUM(distribution%recv_counts)))
943 12636 : DO ientry = 1, SIZE(send_values)
944 : send_values(ientry) = input_data(distribution%send_local_rows(ientry), &
945 12636 : 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 156 : recv_values, distribution%recv_counts, distribution%recv_displacements)
950 12636 : DO ientry = 1, SIZE(recv_values)
951 12480 : slot = (distribution%recv_global_columns(ientry) - 1)/num_pe + 1
952 12636 : local_input(distribution%recv_global_rows(ientry), slot) = recv_values(ientry)
953 : END DO
954 :
955 156 : DEALLOCATE (recv_values, send_values)
956 :
957 312 : 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 156 : 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 156 : POINTER :: output_data
980 :
981 156 : num_pe = preconditioner_env%para_env%num_pe
982 0 : CPASSERT(SIZE(local_output, 1) == distribution%nrow_global)
983 156 : CPASSERT(SIZE(local_output, 2) == SIZE(distribution%global_columns))
984 156 : scale = 1.0_dp
985 156 : IF (PRESENT(scale_factor)) scale = scale_factor
986 : ALLOCATE (send_values(SUM(distribution%recv_counts)), &
987 1404 : recv_values(SUM(distribution%send_counts)))
988 12636 : DO ientry = 1, SIZE(send_values)
989 12480 : slot = (distribution%recv_global_columns(ientry) - 1)/num_pe + 1
990 12636 : 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 156 : recv_values, distribution%send_counts, distribution%send_displacements)
995 :
996 156 : CALL cp_fm_get_info(matrix_out, local_data=output_data)
997 156 : output_data = 0.0_dp
998 12636 : DO ientry = 1, SIZE(recv_values)
999 : output_data(distribution%send_local_rows(ientry), &
1000 12636 : distribution%send_local_columns(ientry)) = scale*recv_values(ientry)
1001 : END DO
1002 156 : DEALLOCATE (recv_values, send_values)
1003 :
1004 156 : 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 468 : 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 156 : NULLIFY (fm_struct)
1022 156 : 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 156 : para_env=preconditioner_env%para_env)
1026 156 : CALL cp_fm_create(input_fm, fm_struct, name="lattice preconditioner input")
1027 156 : CALL cp_fm_create(output_fm, fm_struct, name="lattice preconditioner output")
1028 156 : CALL cp_fm_struct_release(fm_struct)
1029 156 : CALL copy_dbcsr_to_fm(matrix_in, input_fm)
1030 156 : CALL apply_lattice_operator_fm(preconditioner_env, input_fm, output_fm)
1031 156 : CALL copy_fm_to_dbcsr(output_fm, matrix_out)
1032 156 : CALL cp_fm_release(input_fm)
1033 156 : CALL cp_fm_release(output_fm)
1034 :
1035 156 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: DATA
1055 : TYPE(dbcsr_iterator_type) :: iter
1056 : TYPE(dbcsr_type) :: matrix_tmp
1057 :
1058 0 : CALL timeset(routineN, handle)
1059 :
1060 0 : IF (preconditioner_env%low_rank_base == ot_low_rank_base_lattice_fft) THEN
1061 0 : CALL apply_low_rank_dbcsr_via_fm(preconditioner_env, matrix_in, matrix_out)
1062 0 : CALL timestop(handle)
1063 0 : RETURN
1064 : END IF
1065 :
1066 0 : CPASSERT(ASSOCIATED(preconditioner_env%base_dbcsr_matrix))
1067 0 : CPASSERT(ASSOCIATED(preconditioner_env%dbcsr_matrix))
1068 0 : CPASSERT(ASSOCIATED(preconditioner_env%full_evals))
1069 0 : CPASSERT(ASSOCIATED(preconditioner_env%occ_evals))
1070 :
1071 0 : CALL dbcsr_get_info(matrix_in, nfullrows_total=n, nfullcols_total=k)
1072 0 : nocc = SIZE(preconditioner_env%occ_evals)
1073 0 : CPASSERT(n == SIZE(preconditioner_env%full_evals))
1074 0 : 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 0 : preconditioner_env%base_dbcsr_matrix, matrix_in, 0.0_dp, matrix_out)
1079 :
1080 0 : IF (rank_used > 0) THEN
1081 0 : CALL dbcsr_copy(matrix_tmp, matrix_in, name=routineN)
1082 0 : 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 0 : first_row=nocc + 1, last_row=nocc + rank_used)
1086 :
1087 0 : CALL dbcsr_iterator_start(iter, matrix_tmp)
1088 0 : 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 0 : row_offset=row_offset, col_offset=col_offset)
1092 0 : DO j = 1, col_size
1093 0 : DO i = 1, row_size
1094 0 : IF (row_offset + i - 1 >= nocc + 1 .AND. &
1095 0 : 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 0 : - preconditioner_env%spectral_base_scale
1101 0 : DATA(i, j) = DATA(i, j)*correction
1102 : END IF
1103 : END DO
1104 : END DO
1105 : END DO
1106 0 : 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 0 : first_k=nocc + 1, last_k=nocc + rank_used)
1111 0 : CALL dbcsr_release(matrix_tmp)
1112 : END IF
1113 :
1114 0 : CALL timestop(handle)
1115 :
1116 0 : 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 0 : 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 0 : NULLIFY (fm_struct)
1134 0 : 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 0 : para_env=preconditioner_env%para_env)
1138 0 : CALL cp_fm_create(input_fm, fm_struct, name="lattice low-rank input")
1139 0 : CALL cp_fm_create(output_fm, fm_struct, name="lattice low-rank output")
1140 0 : CALL cp_fm_struct_release(fm_struct)
1141 0 : CALL copy_dbcsr_to_fm(matrix_in, input_fm)
1142 0 : CALL apply_low_rank_fm(preconditioner_env, input_fm, output_fm)
1143 0 : CALL copy_fm_to_dbcsr(output_fm, matrix_out)
1144 0 : CALL cp_fm_release(input_fm)
1145 0 : CALL cp_fm_release(output_fm)
1146 :
1147 0 : 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 109880 : 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 27470 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1165 : REAL(KIND=dp) :: dum
1166 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1167 27470 : POINTER :: local_data
1168 : TYPE(cp_fm_type) :: matrix_tmp
1169 :
1170 27470 : CALL timeset(routineN, handle)
1171 :
1172 27470 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
1173 :
1174 27470 : 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 27470 : 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 27470 : matrix_in, 0.0_dp, matrix_tmp)
1181 :
1182 : ! do the right scaling
1183 218410 : DO j = 1, ncol_local
1184 2750241 : DO i = 1, nrow_local
1185 : dum = 1.0_dp/MAX(preconditioner_env%energy_gap, &
1186 2531831 : preconditioner_env%full_evals(row_indices(i)) - preconditioner_env%occ_evals(col_indices(j)))
1187 2722771 : 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 27470 : matrix_tmp, 0.0_dp, matrix_out)
1194 :
1195 27470 : CALL cp_fm_release(matrix_tmp)
1196 :
1197 27470 : CALL timestop(handle)
1198 :
1199 27470 : 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 49916 : 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 24958 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: DATA
1218 : TYPE(dbcsr_iterator_type) :: iter
1219 : TYPE(dbcsr_type) :: matrix_tmp
1220 :
1221 24958 : CALL timeset(routineN, handle)
1222 :
1223 24958 : 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 24958 : matrix_in, 0.0_dp, matrix_tmp)
1226 : ! do the right scaling
1227 24958 : CALL dbcsr_iterator_start(iter, matrix_tmp)
1228 69724 : 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 44766 : row_offset=row_offset, col_offset=col_offset)
1232 432560 : DO j = 1, col_size
1233 3263201 : 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 2855599 : - preconditioner_env%occ_evals(col_offset + j - 1))
1237 3218435 : DATA(i, j) = DATA(i, j)*dum
1238 : END DO
1239 : END DO
1240 : END DO
1241 24958 : CALL dbcsr_iterator_stop(iter)
1242 :
1243 : ! mult back
1244 : CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
1245 24958 : matrix_tmp, 0.0_dp, matrix_out)
1246 24958 : CALL dbcsr_release(matrix_tmp)
1247 24958 : CALL timestop(handle)
1248 :
1249 24958 : 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 48 : 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 12 : CALL timeset(routineN, handle)
1268 12 : CPASSERT(ASSOCIATED(preconditioner_env%occ_rotation))
1269 12 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
1270 :
1271 12 : CALL cp_fm_create(canonical_in, matrix_in%matrix_struct, name=routineN)
1272 12 : 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 12 : preconditioner_env%occ_rotation, 0.0_dp, canonical_in)
1275 12 : IF (preconditioner_env%lattice_fft_active) THEN
1276 0 : 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 12 : preconditioner_env%occ_rotation, 0.0_dp, matrix_out)
1282 12 : CALL cp_fm_release(canonical_in)
1283 12 : CALL cp_fm_release(canonical_out)
1284 12 : CALL timestop(handle)
1285 :
1286 12 : 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 0 : 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 0 : CALL timeset(routineN, handle)
1307 0 : CPASSERT(ASSOCIATED(preconditioner_env%occ_rotation))
1308 0 : CALL dbcsr_get_info(matrix_in, nfullrows_total=n, nfullcols_total=k)
1309 :
1310 0 : IF (preconditioner_env%lattice_fft_active) THEN
1311 0 : NULLIFY (fm_struct)
1312 : CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=k, &
1313 : context=preconditioner_env%ctxt, &
1314 0 : para_env=preconditioner_env%para_env)
1315 0 : CALL cp_fm_create(input_fm, fm_struct, name=routineN)
1316 0 : CALL cp_fm_create(output_fm, fm_struct, name=routineN)
1317 0 : CALL cp_fm_struct_release(fm_struct)
1318 0 : CALL copy_dbcsr_to_fm(matrix_in, input_fm)
1319 0 : CALL apply_full_all_covariant(preconditioner_env, input_fm, output_fm)
1320 0 : CALL copy_fm_to_dbcsr(output_fm, matrix_out)
1321 0 : CALL cp_fm_release(input_fm)
1322 0 : CALL cp_fm_release(output_fm)
1323 0 : CALL timestop(handle)
1324 0 : RETURN
1325 : END IF
1326 :
1327 0 : IF (.NOT. ASSOCIATED(preconditioner_env%occ_rotation_dbcsr)) THEN
1328 0 : 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 0 : sym=dbcsr_type_no_symmetry)
1332 : CALL copy_fm_to_dbcsr(preconditioner_env%occ_rotation, &
1333 0 : preconditioner_env%occ_rotation_dbcsr)
1334 : END IF
1335 :
1336 0 : CALL dbcsr_copy(canonical_in, matrix_in, name=routineN)
1337 0 : CALL dbcsr_set(canonical_in, 0.0_dp)
1338 0 : CALL dbcsr_copy(canonical_out, matrix_out, name=routineN)
1339 0 : CALL dbcsr_set(canonical_out, 0.0_dp)
1340 : CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_in, &
1341 0 : preconditioner_env%occ_rotation_dbcsr, 0.0_dp, canonical_in)
1342 0 : CALL apply_all(preconditioner_env, canonical_in, canonical_out)
1343 : CALL dbcsr_multiply('N', 'T', 1.0_dp, canonical_out, &
1344 0 : preconditioner_env%occ_rotation_dbcsr, 0.0_dp, matrix_out)
1345 0 : CALL dbcsr_release(canonical_in)
1346 0 : CALL dbcsr_release(canonical_out)
1347 0 : CALL timestop(handle)
1348 :
1349 0 : END SUBROUTINE apply_all_covariant
1350 :
1351 : END MODULE preconditioner_apply
|
