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_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
27 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 : copy_fm_to_dbcsr
29 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
30 : cp_fm_scale_and_add
31 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_restore
32 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
33 : cp_fm_struct_release,&
34 : cp_fm_struct_type
35 : USE cp_fm_types, ONLY: cp_fm_create,&
36 : cp_fm_get_info,&
37 : cp_fm_release,&
38 : cp_fm_to_fm,&
39 : cp_fm_type
40 : USE input_constants, ONLY: &
41 : ot_low_rank_base_lattice_fft, ot_low_rank_base_overlap, ot_precond_fermi_low_rank, &
42 : ot_precond_full_all, ot_precond_full_kinetic, ot_precond_full_single, &
43 : ot_precond_full_single_inverse, ot_precond_s_inverse, ot_precond_solver_chebyshev, &
44 : ot_precond_solver_direct, ot_precond_solver_inv_chol, ot_precond_solver_update
45 : USE kinds, ONLY: dp
46 : USE lattice_low_rank_preconditioner, ONLY: apply_lattice_inverse_dense
47 : USE low_rank_preconditioner_model, ONLY: low_rank_inverse_weight
48 : USE mathconstants, ONLY: z_one,&
49 : z_zero
50 : USE parallel_gemm_api, ONLY: parallel_gemm
51 : USE preconditioner_types, ONLY: preconditioner_type
52 : #include "./base/base_uses.f90"
53 :
54 : IMPLICIT NONE
55 : PRIVATE
56 :
57 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner_apply'
58 :
59 : PUBLIC :: apply_preconditioner_cfm_complex, apply_preconditioner_dbcsr_complex, &
60 : apply_preconditioner_fm, apply_preconditioner_dbcsr
61 :
62 : CONTAINS
63 :
64 : ! **************************************************************************************************
65 : !> \brief applies a previously created preconditioner to a full matrix
66 : !> \param preconditioner_env ...
67 : !> \param matrix_in ...
68 : !> \param matrix_out ...
69 : ! **************************************************************************************************
70 50898 : SUBROUTINE apply_preconditioner_fm(preconditioner_env, matrix_in, matrix_out)
71 :
72 : TYPE(preconditioner_type) :: preconditioner_env
73 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
74 :
75 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_fm'
76 :
77 : INTEGER :: handle
78 :
79 50898 : CALL timeset(routineN, handle)
80 :
81 50898 : SELECT CASE (preconditioner_env%in_use)
82 : CASE (0)
83 0 : CPABORT("No preconditioner in use")
84 : CASE (ot_precond_full_single)
85 1320 : CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
86 : CASE (ot_precond_full_all)
87 28380 : CALL apply_full_all(preconditioner_env, matrix_in, matrix_out)
88 : CASE (ot_precond_fermi_low_rank)
89 0 : CALL apply_low_rank_fm(preconditioner_env, matrix_in, matrix_out)
90 : CASE (ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
91 42396 : SELECT CASE (preconditioner_env%solver)
92 : CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
93 21198 : CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
94 : CASE (ot_precond_solver_direct)
95 0 : CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
96 : CASE (ot_precond_solver_chebyshev)
97 0 : CALL apply_chebyshev_fm(preconditioner_env, matrix_in, matrix_out)
98 : CASE DEFAULT
99 21198 : CPABORT("Solver not implemented")
100 : END SELECT
101 : CASE DEFAULT
102 50898 : CPABORT("Unknown preconditioner")
103 : END SELECT
104 :
105 50898 : CALL timestop(handle)
106 :
107 50898 : END SUBROUTINE apply_preconditioner_fm
108 :
109 : ! **************************************************************************************************
110 : !> \brief ...
111 : !> \param preconditioner_env ...
112 : !> \param matrix_in ...
113 : !> \param matrix_out ...
114 : ! **************************************************************************************************
115 77787 : SUBROUTINE apply_preconditioner_dbcsr(preconditioner_env, matrix_in, matrix_out)
116 :
117 : TYPE(preconditioner_type) :: preconditioner_env
118 : TYPE(dbcsr_type) :: matrix_in, matrix_out
119 :
120 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_dbcsr'
121 :
122 : INTEGER :: handle
123 :
124 77787 : CALL timeset(routineN, handle)
125 :
126 77787 : SELECT CASE (preconditioner_env%in_use)
127 : CASE (0)
128 0 : CPABORT("No preconditioner in use")
129 : CASE (ot_precond_full_single)
130 202 : CALL apply_single(preconditioner_env, matrix_in, matrix_out)
131 : CASE (ot_precond_full_all)
132 24638 : CALL apply_all(preconditioner_env, matrix_in, matrix_out)
133 : CASE (ot_precond_fermi_low_rank)
134 96 : CALL apply_low_rank_dbcsr(preconditioner_env, matrix_in, matrix_out)
135 : CASE (ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
136 105664 : SELECT CASE (preconditioner_env%solver)
137 : CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
138 52813 : CALL apply_single(preconditioner_env, matrix_in, matrix_out)
139 : CASE (ot_precond_solver_direct)
140 0 : CPABORT("Apply_full_direct not supported with ot")
141 : !CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
142 : CASE (ot_precond_solver_chebyshev)
143 38 : CALL apply_chebyshev_dbcsr(preconditioner_env, matrix_in, matrix_out)
144 : CASE DEFAULT
145 52851 : CPABORT("Wrong solver")
146 : END SELECT
147 : CASE DEFAULT
148 77787 : CPABORT("Wrong preconditioner")
149 : END SELECT
150 :
151 77787 : CALL timestop(handle)
152 :
153 77787 : END SUBROUTINE apply_preconditioner_dbcsr
154 :
155 : ! **************************************************************************************************
156 : !> \brief Apply a complex k-point orbital preconditioner.
157 : !> \param preconditioner_env complex preconditioner storage
158 : !> \param matrix_in complex input channel
159 : !> \param matrix_out complex output channel
160 : ! **************************************************************************************************
161 57372 : SUBROUTINE apply_preconditioner_cfm_complex(preconditioner_env, matrix_in, matrix_out)
162 :
163 : TYPE(preconditioner_type) :: preconditioner_env
164 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_in, matrix_out
165 :
166 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_cfm_complex'
167 :
168 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
169 19124 : POINTER :: local_data
170 : INTEGER :: handle, i, j, k, n, ncol_local, npre, &
171 : nrow_local
172 19124 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
173 : REAL(KIND=dp) :: scale
174 : TYPE(cp_cfm_type) :: matrix_spectral
175 :
176 19124 : CALL timeset(routineN, handle)
177 :
178 19124 : SELECT CASE (preconditioner_env%in_use)
179 : CASE (ot_precond_full_all, ot_precond_full_single, ot_precond_full_single_inverse, &
180 : ot_precond_full_kinetic, ot_precond_s_inverse)
181 : CASE DEFAULT
182 19124 : CPABORT("Unsupported complex K-point OT preconditioner")
183 : END SELECT
184 19124 : CPASSERT(ASSOCIATED(preconditioner_env%complex_fm))
185 19124 : CALL cp_cfm_get_info(matrix_in, nrow_global=n, ncol_global=k)
186 19124 : CALL cp_cfm_get_info(preconditioner_env%complex_fm, nrow_global=npre)
187 19124 : CPASSERT(n == npre)
188 :
189 19124 : IF (preconditioner_env%in_use == ot_precond_full_all) THEN
190 15177 : CPASSERT(ASSOCIATED(preconditioner_env%full_evals))
191 15177 : CPASSERT(ASSOCIATED(preconditioner_env%occ_evals))
192 15177 : CPASSERT(n == SIZE(preconditioner_env%full_evals))
193 15177 : CPASSERT(k == SIZE(preconditioner_env%occ_evals))
194 : CALL cp_cfm_create(matrix_spectral, matrix_in%matrix_struct, &
195 15177 : name='complex FULL_ALL spectral input')
196 : CALL cp_cfm_gemm('C', 'N', n, k, n, z_one, preconditioner_env%complex_fm, &
197 15177 : matrix_in, z_zero, matrix_spectral)
198 : CALL cp_cfm_get_info(matrix_spectral, nrow_local=nrow_local, ncol_local=ncol_local, &
199 15177 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
200 274434 : DO j = 1, ncol_local
201 7318945 : DO i = 1, nrow_local
202 : scale = 1.0_dp/MAX(preconditioner_env%energy_gap, &
203 : preconditioner_env%full_evals(row_indices(i)) - &
204 7044511 : preconditioner_env%occ_evals(col_indices(j)))
205 7303768 : local_data(i, j) = scale*local_data(i, j)
206 : END DO
207 : END DO
208 : CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, preconditioner_env%complex_fm, &
209 15177 : matrix_spectral, z_zero, matrix_out)
210 15177 : CALL cp_cfm_release(matrix_spectral)
211 : ELSE
212 : CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, preconditioner_env%complex_fm, &
213 3947 : matrix_in, z_zero, matrix_out)
214 : END IF
215 :
216 19124 : CALL timestop(handle)
217 :
218 19124 : END SUBROUTINE apply_preconditioner_cfm_complex
219 :
220 : ! **************************************************************************************************
221 : !> \brief Apply a complex orbital preconditioner to paired real/imaginary DBCSR matrices.
222 : !> \param preconditioner_env complex k-point preconditioner
223 : !> \param matrix_in_re real input channel
224 : !> \param matrix_in_im imaginary input channel
225 : !> \param matrix_out_re real output channel
226 : !> \param matrix_out_im imaginary output channel
227 : ! **************************************************************************************************
228 27784 : SUBROUTINE apply_preconditioner_dbcsr_complex(preconditioner_env, matrix_in_re, matrix_in_im, &
229 : matrix_out_re, matrix_out_im)
230 :
231 : TYPE(preconditioner_type) :: preconditioner_env
232 : TYPE(dbcsr_type) :: matrix_in_re, matrix_in_im, &
233 : matrix_out_re, matrix_out_im
234 :
235 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_dbcsr_complex'
236 :
237 : INTEGER :: handle, k, n
238 : TYPE(cp_cfm_type) :: matrix_in, matrix_out
239 : TYPE(cp_fm_type) :: matrix_in_im_fm, matrix_in_re_fm, &
240 : matrix_out_im_fm, matrix_out_re_fm
241 :
242 3473 : CALL timeset(routineN, handle)
243 :
244 3473 : CALL dbcsr_get_info(matrix_in_re, nfullrows_total=n, nfullcols_total=k)
245 3473 : CPASSERT(ASSOCIATED(preconditioner_env%complex_fm))
246 :
247 : CALL cp_fm_create(matrix_in_re_fm, preconditioner_env%complex_fm%matrix_struct, &
248 3473 : nrow=n, ncol=k, name='complex preconditioner input real')
249 : CALL cp_fm_create(matrix_in_im_fm, preconditioner_env%complex_fm%matrix_struct, &
250 3473 : nrow=n, ncol=k, name='complex preconditioner input imaginary')
251 : CALL cp_fm_create(matrix_out_re_fm, preconditioner_env%complex_fm%matrix_struct, &
252 3473 : nrow=n, ncol=k, name='complex preconditioner output real')
253 : CALL cp_fm_create(matrix_out_im_fm, preconditioner_env%complex_fm%matrix_struct, &
254 3473 : nrow=n, ncol=k, name='complex preconditioner output imaginary')
255 3473 : CALL copy_dbcsr_to_fm(matrix_in_re, matrix_in_re_fm)
256 3473 : CALL copy_dbcsr_to_fm(matrix_in_im, matrix_in_im_fm)
257 :
258 : CALL cp_cfm_create(matrix_in, matrix_in_re_fm%matrix_struct, &
259 3473 : name='complex preconditioner input')
260 : CALL cp_cfm_create(matrix_out, matrix_in_re_fm%matrix_struct, &
261 3473 : name='complex preconditioner output')
262 3473 : CALL cp_fm_to_cfm(matrix_in_re_fm, matrix_in_im_fm, matrix_in)
263 3473 : CALL apply_preconditioner_cfm_complex(preconditioner_env, matrix_in, matrix_out)
264 3473 : CALL cp_cfm_to_fm(matrix_out, matrix_out_re_fm, matrix_out_im_fm)
265 3473 : CALL dbcsr_set(matrix_out_re, 0.0_dp)
266 3473 : CALL dbcsr_set(matrix_out_im, 0.0_dp)
267 3473 : CALL copy_fm_to_dbcsr(matrix_out_re_fm, matrix_out_re)
268 3473 : CALL copy_fm_to_dbcsr(matrix_out_im_fm, matrix_out_im)
269 :
270 3473 : CALL cp_cfm_release(matrix_out)
271 3473 : CALL cp_cfm_release(matrix_in)
272 3473 : CALL cp_fm_release(matrix_out_im_fm)
273 3473 : CALL cp_fm_release(matrix_out_re_fm)
274 3473 : CALL cp_fm_release(matrix_in_im_fm)
275 3473 : CALL cp_fm_release(matrix_in_re_fm)
276 :
277 3473 : CALL timestop(handle)
278 :
279 3473 : END SUBROUTINE apply_preconditioner_dbcsr_complex
280 :
281 : ! **************************************************************************************************
282 : !> \brief apply to full matrix, complete inversion has already been done
283 : !> \param preconditioner_env ...
284 : !> \param matrix_in ...
285 : !> \param matrix_out ...
286 : ! **************************************************************************************************
287 45036 : SUBROUTINE apply_full_single(preconditioner_env, matrix_in, matrix_out)
288 :
289 : TYPE(preconditioner_type) :: preconditioner_env
290 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
291 :
292 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_single'
293 :
294 : INTEGER :: handle, k, n
295 :
296 22518 : CALL timeset(routineN, handle)
297 :
298 22518 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
299 : CALL parallel_gemm('N', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
300 22518 : matrix_in, 0.0_dp, matrix_out)
301 22518 : CALL timestop(handle)
302 :
303 22518 : END SUBROUTINE apply_full_single
304 :
305 : ! **************************************************************************************************
306 : !> \brief apply to dbcsr matrix, complete inversion has already been done
307 : !> \param preconditioner_env ...
308 : !> \param matrix_in ...
309 : !> \param matrix_out ...
310 : ! **************************************************************************************************
311 53015 : SUBROUTINE apply_single(preconditioner_env, matrix_in, matrix_out)
312 :
313 : TYPE(preconditioner_type) :: preconditioner_env
314 : TYPE(dbcsr_type) :: matrix_in, matrix_out
315 :
316 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_single'
317 :
318 : INTEGER :: handle
319 :
320 53015 : CALL timeset(routineN, handle)
321 :
322 53015 : IF (.NOT. ASSOCIATED(preconditioner_env%dbcsr_matrix)) THEN
323 0 : CPABORT("NOT ASSOCIATED preconditioner_env%dbcsr_matrix")
324 : END IF
325 : CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, matrix_in, &
326 53015 : 0.0_dp, matrix_out)
327 :
328 53015 : CALL timestop(handle)
329 :
330 53015 : END SUBROUTINE apply_single
331 :
332 : ! **************************************************************************************************
333 : !> \brief preconditioner contains the factorization, application done by
334 : !> solving the linear system
335 : !> \param preconditioner_env ...
336 : !> \param matrix_in ...
337 : !> \param matrix_out ...
338 : ! **************************************************************************************************
339 0 : SUBROUTINE apply_full_direct(preconditioner_env, matrix_in, matrix_out)
340 :
341 : TYPE(preconditioner_type) :: preconditioner_env
342 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
343 :
344 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_direct'
345 :
346 : INTEGER :: handle, k, n
347 : TYPE(cp_fm_type) :: work
348 :
349 0 : CALL timeset(routineN, handle)
350 :
351 0 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
352 0 : CALL cp_fm_create(work, matrix_in%matrix_struct, name="apply_full_single")
353 : CALL cp_fm_cholesky_restore(matrix_in, k, preconditioner_env%fm, work,&
354 0 : & "SOLVE", transa="T")
355 : CALL cp_fm_cholesky_restore(work, k, preconditioner_env%fm, matrix_out,&
356 0 : & "SOLVE", transa="N")
357 0 : CALL cp_fm_release(work)
358 :
359 0 : CALL timestop(handle)
360 :
361 0 : END SUBROUTINE apply_full_direct
362 :
363 : ! **************************************************************************************************
364 : !> \brief Apply a Chebyshev approximation to the inverse of the stored dense SPD operator.
365 : !> \param preconditioner_env ...
366 : !> \param matrix_in ...
367 : !> \param matrix_out ...
368 : ! **************************************************************************************************
369 0 : SUBROUTINE apply_chebyshev_fm(preconditioner_env, matrix_in, matrix_out)
370 :
371 : TYPE(preconditioner_type) :: preconditioner_env
372 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
373 :
374 : INTEGER :: degree, iteration, k, n
375 : REAL(KIND=dp) :: delta, rho, rho_previous, sigma, theta
376 : TYPE(cp_fm_type) :: direction, residual
377 :
378 0 : CPASSERT(ASSOCIATED(preconditioner_env%fm))
379 0 : degree = preconditioner_env%polynomial_degree
380 0 : CPASSERT(degree >= 1)
381 0 : theta = 0.5_dp*(preconditioner_env%polynomial_max + preconditioner_env%polynomial_min)
382 0 : delta = 0.5_dp*(preconditioner_env%polynomial_max - preconditioner_env%polynomial_min)
383 0 : sigma = theta/delta
384 0 : rho_previous = 1.0_dp/sigma
385 0 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
386 0 : CALL cp_fm_create(direction, matrix_in%matrix_struct, name='Chebyshev direction')
387 0 : CALL cp_fm_create(residual, matrix_in%matrix_struct, name='Chebyshev residual')
388 0 : CALL cp_fm_to_fm(matrix_in, direction)
389 0 : CALL cp_fm_scale(1.0_dp/theta, direction)
390 0 : CALL cp_fm_to_fm(direction, matrix_out)
391 0 : DO iteration = 2, degree
392 0 : CALL cp_fm_to_fm(matrix_in, residual)
393 : CALL parallel_gemm('N', 'N', n, k, n, -1.0_dp, preconditioner_env%fm, &
394 0 : matrix_out, 1.0_dp, residual)
395 0 : rho = 1.0_dp/(2.0_dp*sigma - rho_previous)
396 0 : CALL cp_fm_scale(rho*rho_previous, direction)
397 0 : CALL cp_fm_scale_and_add(1.0_dp, direction, 2.0_dp*rho/delta, residual)
398 0 : CALL cp_fm_scale_and_add(1.0_dp, matrix_out, 1.0_dp, direction)
399 0 : rho_previous = rho
400 : END DO
401 0 : CALL cp_fm_release(residual)
402 0 : CALL cp_fm_release(direction)
403 :
404 0 : END SUBROUTINE apply_chebyshev_fm
405 :
406 : ! **************************************************************************************************
407 : !> \brief Apply a Chebyshev approximation to the inverse of the stored sparse SPD operator.
408 : !> \param preconditioner_env ...
409 : !> \param matrix_in ...
410 : !> \param matrix_out ...
411 : ! **************************************************************************************************
412 38 : SUBROUTINE apply_chebyshev_dbcsr(preconditioner_env, matrix_in, matrix_out)
413 :
414 : TYPE(preconditioner_type) :: preconditioner_env
415 : TYPE(dbcsr_type) :: matrix_in, matrix_out
416 :
417 : INTEGER :: degree, iteration
418 : REAL(KIND=dp) :: delta, rho, rho_previous, sigma, theta
419 : TYPE(dbcsr_type) :: direction, residual
420 :
421 0 : CPASSERT(ASSOCIATED(preconditioner_env%dbcsr_matrix))
422 38 : degree = preconditioner_env%polynomial_degree
423 38 : CPASSERT(degree >= 1)
424 38 : theta = 0.5_dp*(preconditioner_env%polynomial_max + preconditioner_env%polynomial_min)
425 38 : delta = 0.5_dp*(preconditioner_env%polynomial_max - preconditioner_env%polynomial_min)
426 38 : sigma = theta/delta
427 38 : rho_previous = 1.0_dp/sigma
428 38 : CALL dbcsr_copy(direction, matrix_in, name='Chebyshev direction')
429 38 : CALL dbcsr_copy(residual, matrix_in, name='Chebyshev residual')
430 38 : CALL dbcsr_set(matrix_out, 0.0_dp)
431 38 : CALL dbcsr_add(matrix_out, direction, 1.0_dp, 1.0_dp/theta)
432 38 : CALL dbcsr_set(direction, 0.0_dp)
433 38 : CALL dbcsr_add(direction, matrix_in, 1.0_dp, 1.0_dp/theta)
434 304 : DO iteration = 2, degree
435 266 : CALL dbcsr_set(residual, 0.0_dp)
436 266 : CALL dbcsr_add(residual, matrix_in, 1.0_dp, 1.0_dp)
437 : CALL dbcsr_multiply('N', 'N', -1.0_dp, preconditioner_env%dbcsr_matrix, &
438 266 : matrix_out, 1.0_dp, residual)
439 266 : rho = 1.0_dp/(2.0_dp*sigma - rho_previous)
440 266 : CALL dbcsr_add(direction, residual, rho*rho_previous, 2.0_dp*rho/delta)
441 266 : CALL dbcsr_add(matrix_out, direction, 1.0_dp, 1.0_dp)
442 304 : rho_previous = rho
443 : END DO
444 38 : CALL dbcsr_release(residual)
445 38 : CALL dbcsr_release(direction)
446 :
447 38 : END SUBROUTINE apply_chebyshev_dbcsr
448 :
449 : ! **************************************************************************************************
450 : !> \brief Applies the overlap-inverse base and a bounded low-rank spectral correction.
451 : !> \param preconditioner_env ...
452 : !> \param matrix_in ...
453 : !> \param matrix_out ...
454 : ! **************************************************************************************************
455 96 : SUBROUTINE apply_low_rank_fm(preconditioner_env, matrix_in, matrix_out)
456 :
457 : TYPE(preconditioner_type) :: preconditioner_env
458 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
459 :
460 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_low_rank_fm'
461 :
462 : INTEGER :: handle, i, j, k, n, ncol_local, nocc, &
463 : nrow_local, rank_used
464 48 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
465 : REAL(KIND=dp) :: correction
466 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
467 48 : POINTER :: local_data
468 : TYPE(cp_fm_type) :: matrix_tmp
469 :
470 48 : CALL timeset(routineN, handle)
471 :
472 48 : CPASSERT(ASSOCIATED(preconditioner_env%fm))
473 48 : CPASSERT(ASSOCIATED(preconditioner_env%full_evals))
474 48 : CPASSERT(ASSOCIATED(preconditioner_env%occ_evals))
475 :
476 48 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
477 48 : nocc = SIZE(preconditioner_env%occ_evals)
478 48 : CPASSERT(n == SIZE(preconditioner_env%full_evals))
479 48 : rank_used = MIN(preconditioner_env%spectral_rank, n - nocc)
480 :
481 : ! With the exact base V V^T = S^-1. The optional lattice base replaces
482 : ! S^-1 by the inverse of its block-circulant translational projection:
483 : ! P = B^-1/window + V_r diag(1/max(gap, epsilon_a-mu) - 1/window) V_r^T.
484 : ! The common reference mu makes the operator right-covariant for any number
485 : ! of orbital columns. The lattice path retains only non-negative corrections.
486 48 : SELECT CASE (preconditioner_env%low_rank_base)
487 : CASE (ot_low_rank_base_overlap)
488 0 : CPASSERT(ASSOCIATED(preconditioner_env%base_fm))
489 : CALL parallel_gemm('N', 'N', n, k, n, preconditioner_env%spectral_base_scale, &
490 0 : preconditioner_env%base_fm, matrix_in, 0.0_dp, matrix_out)
491 : CASE (ot_low_rank_base_lattice_fft)
492 48 : CALL apply_lattice_base_fm(preconditioner_env, matrix_in, matrix_out)
493 : CASE DEFAULT
494 48 : CPABORT("Unknown FERMI_LOW_RANK base operator during application")
495 : END SELECT
496 :
497 48 : IF (rank_used > 0) THEN
498 : CALL cp_fm_create(matrix_tmp, matrix_in%matrix_struct, name=routineN, &
499 48 : nrow=rank_used, ncol=k)
500 : CALL cp_fm_get_info(matrix_tmp, nrow_local=nrow_local, ncol_local=ncol_local, &
501 48 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
502 :
503 : CALL parallel_gemm('T', 'N', rank_used, k, n, 1.0_dp, preconditioner_env%fm, &
504 48 : matrix_in, 0.0_dp, matrix_tmp, a_first_col=nocc + 1)
505 :
506 432 : DO j = 1, ncol_local
507 7344 : DO i = 1, nrow_local
508 : correction = low_rank_inverse_weight( &
509 : preconditioner_env%full_evals(nocc + row_indices(i)), &
510 : preconditioner_env%spectral_reference, &
511 : preconditioner_env%energy_gap) &
512 6912 : - preconditioner_env%spectral_base_scale
513 7296 : local_data(i, j) = local_data(i, j)*correction
514 : END DO
515 : END DO
516 :
517 : CALL parallel_gemm('N', 'N', n, k, rank_used, 1.0_dp, preconditioner_env%fm, &
518 48 : matrix_tmp, 1.0_dp, matrix_out, a_first_col=nocc + 1)
519 48 : CALL cp_fm_release(matrix_tmp)
520 : END IF
521 :
522 48 : CALL timestop(handle)
523 :
524 48 : END SUBROUTINE apply_low_rank_fm
525 :
526 : ! **************************************************************************************************
527 : !> \brief Apply the replicated-cell FFT base to a distributed full matrix.
528 : !> \param preconditioner_env ...
529 : !> \param matrix_in ...
530 : !> \param matrix_out ...
531 : ! **************************************************************************************************
532 48 : SUBROUTINE apply_lattice_base_fm(preconditioner_env, matrix_in, matrix_out)
533 :
534 : TYPE(preconditioner_type) :: preconditioner_env
535 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
536 :
537 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_lattice_base_fm'
538 :
539 : INTEGER :: handle, i, j, k, n, ncol_local, &
540 : nrow_local
541 48 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
542 : LOGICAL :: used_fft
543 : REAL(KIND=dp) :: imaginary_residual, scale
544 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: global_input, global_output
545 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
546 48 : POINTER :: input_data, output_data
547 :
548 48 : CALL timeset(routineN, handle)
549 48 : CPASSERT(ALLOCATED(preconditioner_env%lattice_inverse_k))
550 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k, nrow_local=nrow_local, &
551 : ncol_local=ncol_local, row_indices=row_indices, &
552 48 : col_indices=col_indices, local_data=input_data)
553 48 : CALL cp_fm_get_info(matrix_out, local_data=output_data)
554 192 : ALLOCATE (global_input(n, k), source=0.0_dp)
555 144 : ALLOCATE (global_output(n, k))
556 432 : DO j = 1, ncol_local
557 13104 : DO i = 1, nrow_local
558 13056 : global_input(row_indices(i), col_indices(j)) = input_data(i, j)
559 : END DO
560 : END DO
561 48 : CALL preconditioner_env%para_env%sum(global_input)
562 : CALL apply_lattice_inverse_dense(preconditioner_env%lattice_inverse_k, &
563 : preconditioner_env%lattice_dims, global_input, global_output, &
564 48 : used_fft, imaginary_residual)
565 48 : IF (.NOT. used_fft) CPABORT("FERMI_LOW_RANK lattice inverse failed to use FFTs")
566 25776 : scale = MAX(1.0_dp, MAXVAL(ABS(global_output)))
567 48 : IF (imaginary_residual > 1.0E-10_dp*scale) THEN
568 0 : CPABORT("FERMI_LOW_RANK lattice inverse produced a non-real Gamma-point result")
569 : END IF
570 432 : DO j = 1, ncol_local
571 13104 : DO i = 1, nrow_local
572 : output_data(i, j) = preconditioner_env%spectral_base_scale* &
573 13056 : global_output(row_indices(i), col_indices(j))
574 : END DO
575 : END DO
576 48 : DEALLOCATE (global_input, global_output)
577 48 : CALL timestop(handle)
578 :
579 96 : END SUBROUTINE apply_lattice_base_fm
580 :
581 : ! **************************************************************************************************
582 : !> \brief DBCSR variant of the bounded low-rank spectral correction.
583 : !> \param preconditioner_env ...
584 : !> \param matrix_in ...
585 : !> \param matrix_out ...
586 : ! **************************************************************************************************
587 144 : SUBROUTINE apply_low_rank_dbcsr(preconditioner_env, matrix_in, matrix_out)
588 :
589 : TYPE(preconditioner_type) :: preconditioner_env
590 : TYPE(dbcsr_type) :: matrix_in, matrix_out
591 :
592 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_low_rank_dbcsr'
593 :
594 : INTEGER :: col, col_offset, col_size, handle, i, j, &
595 : k, n, nocc, rank_used, row, &
596 : row_offset, row_size
597 : REAL(KIND=dp) :: correction
598 96 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: DATA
599 : TYPE(dbcsr_iterator_type) :: iter
600 : TYPE(dbcsr_type) :: matrix_tmp
601 :
602 96 : CALL timeset(routineN, handle)
603 :
604 96 : IF (preconditioner_env%low_rank_base == ot_low_rank_base_lattice_fft) THEN
605 48 : CALL apply_low_rank_dbcsr_via_fm(preconditioner_env, matrix_in, matrix_out)
606 48 : CALL timestop(handle)
607 48 : RETURN
608 : END IF
609 :
610 48 : CPASSERT(ASSOCIATED(preconditioner_env%base_dbcsr_matrix))
611 48 : CPASSERT(ASSOCIATED(preconditioner_env%dbcsr_matrix))
612 48 : CPASSERT(ASSOCIATED(preconditioner_env%full_evals))
613 48 : CPASSERT(ASSOCIATED(preconditioner_env%occ_evals))
614 :
615 48 : CALL dbcsr_get_info(matrix_in, nfullrows_total=n, nfullcols_total=k)
616 48 : nocc = SIZE(preconditioner_env%occ_evals)
617 48 : CPASSERT(n == SIZE(preconditioner_env%full_evals))
618 48 : rank_used = MIN(preconditioner_env%spectral_rank, n - nocc)
619 :
620 : ! DBCSR form of the same overlap-inverse base plus positive low-rank update.
621 : CALL dbcsr_multiply('N', 'N', preconditioner_env%spectral_base_scale, &
622 48 : preconditioner_env%base_dbcsr_matrix, matrix_in, 0.0_dp, matrix_out)
623 :
624 48 : IF (rank_used > 0) THEN
625 48 : CALL dbcsr_copy(matrix_tmp, matrix_in, name=routineN)
626 48 : CALL dbcsr_set(matrix_tmp, 0.0_dp)
627 : CALL dbcsr_multiply('T', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
628 : matrix_in, 0.0_dp, matrix_tmp, &
629 48 : first_row=nocc + 1, last_row=nocc + rank_used)
630 :
631 48 : CALL dbcsr_iterator_start(iter, matrix_tmp)
632 288 : DO WHILE (dbcsr_iterator_blocks_left(iter))
633 : CALL dbcsr_iterator_next_block(iter, row, col, DATA, &
634 : row_size=row_size, col_size=col_size, &
635 240 : row_offset=row_offset, col_offset=col_offset)
636 2208 : DO j = 1, col_size
637 14832 : DO i = 1, row_size
638 12672 : IF (row_offset + i - 1 >= nocc + 1 .AND. &
639 1920 : row_offset + i - 1 <= nocc + rank_used) THEN
640 : correction = low_rank_inverse_weight( &
641 : preconditioner_env%full_evals(row_offset + i - 1), &
642 : preconditioner_env%spectral_reference, &
643 : preconditioner_env%energy_gap) &
644 9216 : - preconditioner_env%spectral_base_scale
645 9216 : DATA(i, j) = DATA(i, j)*correction
646 : END IF
647 : END DO
648 : END DO
649 : END DO
650 48 : CALL dbcsr_iterator_stop(iter)
651 :
652 : CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
653 : matrix_tmp, 1.0_dp, matrix_out, &
654 48 : first_k=nocc + 1, last_k=nocc + rank_used)
655 48 : CALL dbcsr_release(matrix_tmp)
656 : END IF
657 :
658 48 : CALL timestop(handle)
659 :
660 96 : END SUBROUTINE apply_low_rank_dbcsr
661 :
662 : ! **************************************************************************************************
663 : !> \brief Bridge DBCSR OT gradients to the full-matrix lattice-FFT implementation.
664 : !> \param preconditioner_env ...
665 : !> \param matrix_in ...
666 : !> \param matrix_out ...
667 : ! **************************************************************************************************
668 144 : SUBROUTINE apply_low_rank_dbcsr_via_fm(preconditioner_env, matrix_in, matrix_out)
669 :
670 : TYPE(preconditioner_type) :: preconditioner_env
671 : TYPE(dbcsr_type) :: matrix_in, matrix_out
672 :
673 : INTEGER :: k, n
674 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
675 : TYPE(cp_fm_type) :: input_fm, output_fm
676 :
677 48 : NULLIFY (fm_struct)
678 48 : CALL dbcsr_get_info(matrix_in, nfullrows_total=n, nfullcols_total=k)
679 : CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=k, &
680 : context=preconditioner_env%ctxt, &
681 48 : para_env=preconditioner_env%para_env)
682 48 : CALL cp_fm_create(input_fm, fm_struct, name="lattice low-rank input")
683 48 : CALL cp_fm_create(output_fm, fm_struct, name="lattice low-rank output")
684 48 : CALL cp_fm_struct_release(fm_struct)
685 48 : CALL copy_dbcsr_to_fm(matrix_in, input_fm)
686 48 : CALL apply_low_rank_fm(preconditioner_env, input_fm, output_fm)
687 48 : CALL copy_fm_to_dbcsr(output_fm, matrix_out)
688 48 : CALL cp_fm_release(input_fm)
689 48 : CALL cp_fm_release(output_fm)
690 :
691 48 : END SUBROUTINE apply_low_rank_dbcsr_via_fm
692 :
693 : ! **************************************************************************************************
694 : !> \brief full all to a full matrix
695 : !> \param preconditioner_env ...
696 : !> \param matrix_in ...
697 : !> \param matrix_out ...
698 : ! **************************************************************************************************
699 113520 : SUBROUTINE apply_full_all(preconditioner_env, matrix_in, matrix_out)
700 :
701 : TYPE(preconditioner_type) :: preconditioner_env
702 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
703 :
704 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_all'
705 :
706 : INTEGER :: handle, i, j, k, n, ncol_local, &
707 : nrow_local
708 28380 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
709 : REAL(KIND=dp) :: dum
710 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
711 28380 : POINTER :: local_data
712 : TYPE(cp_fm_type) :: matrix_tmp
713 :
714 28380 : CALL timeset(routineN, handle)
715 :
716 28380 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
717 :
718 28380 : CALL cp_fm_create(matrix_tmp, matrix_in%matrix_struct, name="apply_full_all")
719 : CALL cp_fm_get_info(matrix_tmp, nrow_local=nrow_local, ncol_local=ncol_local, &
720 28380 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
721 :
722 : !
723 : CALL parallel_gemm('T', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
724 28380 : matrix_in, 0.0_dp, matrix_tmp)
725 :
726 : ! do the right scaling
727 224454 : DO j = 1, ncol_local
728 2822205 : DO i = 1, nrow_local
729 : dum = 1.0_dp/MAX(preconditioner_env%energy_gap, &
730 2597751 : preconditioner_env%full_evals(row_indices(i)) - preconditioner_env%occ_evals(col_indices(j)))
731 2793825 : local_data(i, j) = local_data(i, j)*dum
732 : END DO
733 : END DO
734 :
735 : ! mult back
736 : CALL parallel_gemm('N', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
737 28380 : matrix_tmp, 0.0_dp, matrix_out)
738 :
739 28380 : CALL cp_fm_release(matrix_tmp)
740 :
741 28380 : CALL timestop(handle)
742 :
743 28380 : END SUBROUTINE apply_full_all
744 :
745 : ! **************************************************************************************************
746 : !> \brief full all to a dbcsr matrix
747 : !> \param preconditioner_env ...
748 : !> \param matrix_in ...
749 : !> \param matrix_out ...
750 : ! **************************************************************************************************
751 49276 : SUBROUTINE apply_all(preconditioner_env, matrix_in, matrix_out)
752 :
753 : TYPE(preconditioner_type) :: preconditioner_env
754 : TYPE(dbcsr_type) :: matrix_in, matrix_out
755 :
756 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_all'
757 :
758 : INTEGER :: col, col_offset, col_size, handle, i, j, &
759 : row, row_offset, row_size
760 : REAL(KIND=dp) :: dum
761 24638 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: DATA
762 : TYPE(dbcsr_iterator_type) :: iter
763 : TYPE(dbcsr_type) :: matrix_tmp
764 :
765 24638 : CALL timeset(routineN, handle)
766 :
767 24638 : CALL dbcsr_copy(matrix_tmp, matrix_in, name="apply_full_all")
768 : CALL dbcsr_multiply('T', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
769 24638 : matrix_in, 0.0_dp, matrix_tmp)
770 : ! do the right scaling
771 24638 : CALL dbcsr_iterator_start(iter, matrix_tmp)
772 69121 : DO WHILE (dbcsr_iterator_blocks_left(iter))
773 : CALL dbcsr_iterator_next_block(iter, row, col, DATA, &
774 : row_size=row_size, col_size=col_size, &
775 44483 : row_offset=row_offset, col_offset=col_offset)
776 431918 : DO j = 1, col_size
777 3266916 : DO i = 1, row_size
778 : dum = 1.0_dp/MAX(preconditioner_env%energy_gap, &
779 : preconditioner_env%full_evals(row_offset + i - 1) &
780 2859636 : - preconditioner_env%occ_evals(col_offset + j - 1))
781 3222433 : DATA(i, j) = DATA(i, j)*dum
782 : END DO
783 : END DO
784 : END DO
785 24638 : CALL dbcsr_iterator_stop(iter)
786 :
787 : ! mult back
788 : CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
789 24638 : matrix_tmp, 0.0_dp, matrix_out)
790 24638 : CALL dbcsr_release(matrix_tmp)
791 24638 : CALL timestop(handle)
792 :
793 24638 : END SUBROUTINE apply_all
794 :
795 : END MODULE preconditioner_apply
|
