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_types, ONLY: cp_fm_create,&
33 : cp_fm_get_info,&
34 : cp_fm_release,&
35 : cp_fm_to_fm,&
36 : cp_fm_type
37 : USE input_constants, ONLY: &
38 : ot_precond_full_all, ot_precond_full_kinetic, ot_precond_full_single, &
39 : ot_precond_full_single_inverse, ot_precond_s_inverse, ot_precond_solver_chebyshev, &
40 : ot_precond_solver_direct, ot_precond_solver_inv_chol, ot_precond_solver_update
41 : USE kinds, ONLY: dp
42 : USE mathconstants, ONLY: z_one,&
43 : z_zero
44 : USE parallel_gemm_api, ONLY: parallel_gemm
45 : USE preconditioner_types, ONLY: preconditioner_type
46 : #include "./base/base_uses.f90"
47 :
48 : IMPLICIT NONE
49 : PRIVATE
50 :
51 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner_apply'
52 :
53 : PUBLIC :: apply_preconditioner_cfm_complex, apply_preconditioner_dbcsr_complex, &
54 : apply_preconditioner_fm, apply_preconditioner_dbcsr
55 :
56 : CONTAINS
57 :
58 : ! **************************************************************************************************
59 : !> \brief applies a previously created preconditioner to a full matrix
60 : !> \param preconditioner_env ...
61 : !> \param matrix_in ...
62 : !> \param matrix_out ...
63 : ! **************************************************************************************************
64 50362 : SUBROUTINE apply_preconditioner_fm(preconditioner_env, matrix_in, matrix_out)
65 :
66 : TYPE(preconditioner_type) :: preconditioner_env
67 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
68 :
69 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_fm'
70 :
71 : INTEGER :: handle
72 :
73 50362 : CALL timeset(routineN, handle)
74 :
75 50362 : SELECT CASE (preconditioner_env%in_use)
76 : CASE (0)
77 0 : CPABORT("No preconditioner in use")
78 : CASE (ot_precond_full_single)
79 1320 : CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
80 : CASE (ot_precond_full_all)
81 28380 : CALL apply_full_all(preconditioner_env, matrix_in, matrix_out)
82 : CASE (ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
83 41324 : SELECT CASE (preconditioner_env%solver)
84 : CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
85 20662 : CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
86 : CASE (ot_precond_solver_direct)
87 0 : CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
88 : CASE (ot_precond_solver_chebyshev)
89 0 : CALL apply_chebyshev_fm(preconditioner_env, matrix_in, matrix_out)
90 : CASE DEFAULT
91 20662 : CPABORT("Solver not implemented")
92 : END SELECT
93 : CASE DEFAULT
94 50362 : CPABORT("Unknown preconditioner")
95 : END SELECT
96 :
97 50362 : CALL timestop(handle)
98 :
99 50362 : END SUBROUTINE apply_preconditioner_fm
100 :
101 : ! **************************************************************************************************
102 : !> \brief ...
103 : !> \param preconditioner_env ...
104 : !> \param matrix_in ...
105 : !> \param matrix_out ...
106 : ! **************************************************************************************************
107 76441 : SUBROUTINE apply_preconditioner_dbcsr(preconditioner_env, matrix_in, matrix_out)
108 :
109 : TYPE(preconditioner_type) :: preconditioner_env
110 : TYPE(dbcsr_type) :: matrix_in, matrix_out
111 :
112 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_dbcsr'
113 :
114 : INTEGER :: handle
115 :
116 76441 : CALL timeset(routineN, handle)
117 :
118 76441 : SELECT CASE (preconditioner_env%in_use)
119 : CASE (0)
120 0 : CPABORT("No preconditioner in use")
121 : CASE (ot_precond_full_single)
122 202 : CALL apply_single(preconditioner_env, matrix_in, matrix_out)
123 : CASE (ot_precond_full_all)
124 23636 : CALL apply_all(preconditioner_env, matrix_in, matrix_out)
125 : CASE (ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
126 105168 : SELECT CASE (preconditioner_env%solver)
127 : CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
128 52565 : CALL apply_single(preconditioner_env, matrix_in, matrix_out)
129 : CASE (ot_precond_solver_direct)
130 0 : CPABORT("Apply_full_direct not supported with ot")
131 : !CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
132 : CASE (ot_precond_solver_chebyshev)
133 38 : CALL apply_chebyshev_dbcsr(preconditioner_env, matrix_in, matrix_out)
134 : CASE DEFAULT
135 52603 : CPABORT("Wrong solver")
136 : END SELECT
137 : CASE DEFAULT
138 76441 : CPABORT("Wrong preconditioner")
139 : END SELECT
140 :
141 76441 : CALL timestop(handle)
142 :
143 76441 : END SUBROUTINE apply_preconditioner_dbcsr
144 :
145 : ! **************************************************************************************************
146 : !> \brief Apply a complex k-point orbital preconditioner.
147 : !> \param preconditioner_env complex preconditioner storage
148 : !> \param matrix_in complex input channel
149 : !> \param matrix_out complex output channel
150 : ! **************************************************************************************************
151 57864 : SUBROUTINE apply_preconditioner_cfm_complex(preconditioner_env, matrix_in, matrix_out)
152 :
153 : TYPE(preconditioner_type) :: preconditioner_env
154 : TYPE(cp_cfm_type), INTENT(IN) :: matrix_in, matrix_out
155 :
156 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_cfm_complex'
157 :
158 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
159 19288 : POINTER :: local_data
160 : INTEGER :: handle, i, j, k, n, ncol_local, npre, &
161 : nrow_local
162 19288 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
163 : REAL(KIND=dp) :: scale
164 : TYPE(cp_cfm_type) :: matrix_spectral
165 :
166 19288 : CALL timeset(routineN, handle)
167 :
168 19288 : SELECT CASE (preconditioner_env%in_use)
169 : CASE (ot_precond_full_all, ot_precond_full_single, ot_precond_full_single_inverse, &
170 : ot_precond_full_kinetic, ot_precond_s_inverse)
171 : CASE DEFAULT
172 19288 : CPABORT("Unsupported complex K-point OT preconditioner")
173 : END SELECT
174 19288 : CPASSERT(ASSOCIATED(preconditioner_env%complex_fm))
175 19288 : CALL cp_cfm_get_info(matrix_in, nrow_global=n, ncol_global=k)
176 19288 : CALL cp_cfm_get_info(preconditioner_env%complex_fm, nrow_global=npre)
177 19288 : CPASSERT(n == npre)
178 :
179 19288 : IF (preconditioner_env%in_use == ot_precond_full_all) THEN
180 15321 : CPASSERT(ASSOCIATED(preconditioner_env%full_evals))
181 15321 : CPASSERT(ASSOCIATED(preconditioner_env%occ_evals))
182 15321 : CPASSERT(n == SIZE(preconditioner_env%full_evals))
183 15321 : CPASSERT(k == SIZE(preconditioner_env%occ_evals))
184 : CALL cp_cfm_create(matrix_spectral, matrix_in%matrix_struct, &
185 15321 : name='complex FULL_ALL spectral input')
186 : CALL cp_cfm_gemm('C', 'N', n, k, n, z_one, preconditioner_env%complex_fm, &
187 15321 : matrix_in, z_zero, matrix_spectral)
188 : CALL cp_cfm_get_info(matrix_spectral, nrow_local=nrow_local, ncol_local=ncol_local, &
189 15321 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
190 276066 : DO j = 1, ncol_local
191 7408673 : DO i = 1, nrow_local
192 : scale = 1.0_dp/MAX(preconditioner_env%energy_gap, &
193 : preconditioner_env%full_evals(row_indices(i)) - &
194 7132607 : preconditioner_env%occ_evals(col_indices(j)))
195 7393352 : local_data(i, j) = scale*local_data(i, j)
196 : END DO
197 : END DO
198 : CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, preconditioner_env%complex_fm, &
199 15321 : matrix_spectral, z_zero, matrix_out)
200 15321 : CALL cp_cfm_release(matrix_spectral)
201 : ELSE
202 : CALL cp_cfm_gemm('N', 'N', n, k, n, z_one, preconditioner_env%complex_fm, &
203 3967 : matrix_in, z_zero, matrix_out)
204 : END IF
205 :
206 19288 : CALL timestop(handle)
207 :
208 19288 : END SUBROUTINE apply_preconditioner_cfm_complex
209 :
210 : ! **************************************************************************************************
211 : !> \brief Apply a complex orbital preconditioner to paired real/imaginary DBCSR matrices.
212 : !> \param preconditioner_env complex k-point preconditioner
213 : !> \param matrix_in_re real input channel
214 : !> \param matrix_in_im imaginary input channel
215 : !> \param matrix_out_re real output channel
216 : !> \param matrix_out_im imaginary output channel
217 : ! **************************************************************************************************
218 29096 : SUBROUTINE apply_preconditioner_dbcsr_complex(preconditioner_env, matrix_in_re, matrix_in_im, &
219 : matrix_out_re, matrix_out_im)
220 :
221 : TYPE(preconditioner_type) :: preconditioner_env
222 : TYPE(dbcsr_type) :: matrix_in_re, matrix_in_im, &
223 : matrix_out_re, matrix_out_im
224 :
225 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_dbcsr_complex'
226 :
227 : INTEGER :: handle, k, n
228 : TYPE(cp_cfm_type) :: matrix_in, matrix_out
229 : TYPE(cp_fm_type) :: matrix_in_im_fm, matrix_in_re_fm, &
230 : matrix_out_im_fm, matrix_out_re_fm
231 :
232 3637 : CALL timeset(routineN, handle)
233 :
234 3637 : CALL dbcsr_get_info(matrix_in_re, nfullrows_total=n, nfullcols_total=k)
235 3637 : CPASSERT(ASSOCIATED(preconditioner_env%complex_fm))
236 :
237 : CALL cp_fm_create(matrix_in_re_fm, preconditioner_env%complex_fm%matrix_struct, &
238 3637 : nrow=n, ncol=k, name='complex preconditioner input real')
239 : CALL cp_fm_create(matrix_in_im_fm, preconditioner_env%complex_fm%matrix_struct, &
240 3637 : nrow=n, ncol=k, name='complex preconditioner input imaginary')
241 : CALL cp_fm_create(matrix_out_re_fm, preconditioner_env%complex_fm%matrix_struct, &
242 3637 : nrow=n, ncol=k, name='complex preconditioner output real')
243 : CALL cp_fm_create(matrix_out_im_fm, preconditioner_env%complex_fm%matrix_struct, &
244 3637 : nrow=n, ncol=k, name='complex preconditioner output imaginary')
245 3637 : CALL copy_dbcsr_to_fm(matrix_in_re, matrix_in_re_fm)
246 3637 : CALL copy_dbcsr_to_fm(matrix_in_im, matrix_in_im_fm)
247 :
248 : CALL cp_cfm_create(matrix_in, matrix_in_re_fm%matrix_struct, &
249 3637 : name='complex preconditioner input')
250 : CALL cp_cfm_create(matrix_out, matrix_in_re_fm%matrix_struct, &
251 3637 : name='complex preconditioner output')
252 3637 : CALL cp_fm_to_cfm(matrix_in_re_fm, matrix_in_im_fm, matrix_in)
253 3637 : CALL apply_preconditioner_cfm_complex(preconditioner_env, matrix_in, matrix_out)
254 3637 : CALL cp_cfm_to_fm(matrix_out, matrix_out_re_fm, matrix_out_im_fm)
255 3637 : CALL dbcsr_set(matrix_out_re, 0.0_dp)
256 3637 : CALL dbcsr_set(matrix_out_im, 0.0_dp)
257 3637 : CALL copy_fm_to_dbcsr(matrix_out_re_fm, matrix_out_re)
258 3637 : CALL copy_fm_to_dbcsr(matrix_out_im_fm, matrix_out_im)
259 :
260 3637 : CALL cp_cfm_release(matrix_out)
261 3637 : CALL cp_cfm_release(matrix_in)
262 3637 : CALL cp_fm_release(matrix_out_im_fm)
263 3637 : CALL cp_fm_release(matrix_out_re_fm)
264 3637 : CALL cp_fm_release(matrix_in_im_fm)
265 3637 : CALL cp_fm_release(matrix_in_re_fm)
266 :
267 3637 : CALL timestop(handle)
268 :
269 3637 : END SUBROUTINE apply_preconditioner_dbcsr_complex
270 :
271 : ! **************************************************************************************************
272 : !> \brief apply to full matrix, complete inversion has already been done
273 : !> \param preconditioner_env ...
274 : !> \param matrix_in ...
275 : !> \param matrix_out ...
276 : ! **************************************************************************************************
277 43964 : SUBROUTINE apply_full_single(preconditioner_env, matrix_in, matrix_out)
278 :
279 : TYPE(preconditioner_type) :: preconditioner_env
280 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
281 :
282 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_single'
283 :
284 : INTEGER :: handle, k, n
285 :
286 21982 : CALL timeset(routineN, handle)
287 :
288 21982 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
289 : CALL parallel_gemm('N', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
290 21982 : matrix_in, 0.0_dp, matrix_out)
291 21982 : CALL timestop(handle)
292 :
293 21982 : END SUBROUTINE apply_full_single
294 :
295 : ! **************************************************************************************************
296 : !> \brief apply to dbcsr matrix, complete inversion has already been done
297 : !> \param preconditioner_env ...
298 : !> \param matrix_in ...
299 : !> \param matrix_out ...
300 : ! **************************************************************************************************
301 52767 : SUBROUTINE apply_single(preconditioner_env, matrix_in, matrix_out)
302 :
303 : TYPE(preconditioner_type) :: preconditioner_env
304 : TYPE(dbcsr_type) :: matrix_in, matrix_out
305 :
306 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_single'
307 :
308 : INTEGER :: handle
309 :
310 52767 : CALL timeset(routineN, handle)
311 :
312 52767 : IF (.NOT. ASSOCIATED(preconditioner_env%dbcsr_matrix)) THEN
313 0 : CPABORT("NOT ASSOCIATED preconditioner_env%dbcsr_matrix")
314 : END IF
315 : CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, matrix_in, &
316 52767 : 0.0_dp, matrix_out)
317 :
318 52767 : CALL timestop(handle)
319 :
320 52767 : END SUBROUTINE apply_single
321 :
322 : ! **************************************************************************************************
323 : !> \brief preconditioner contains the factorization, application done by
324 : !> solving the linear system
325 : !> \param preconditioner_env ...
326 : !> \param matrix_in ...
327 : !> \param matrix_out ...
328 : ! **************************************************************************************************
329 0 : SUBROUTINE apply_full_direct(preconditioner_env, matrix_in, matrix_out)
330 :
331 : TYPE(preconditioner_type) :: preconditioner_env
332 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
333 :
334 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_direct'
335 :
336 : INTEGER :: handle, k, n
337 : TYPE(cp_fm_type) :: work
338 :
339 0 : CALL timeset(routineN, handle)
340 :
341 0 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
342 0 : CALL cp_fm_create(work, matrix_in%matrix_struct, name="apply_full_single")
343 : CALL cp_fm_cholesky_restore(matrix_in, k, preconditioner_env%fm, work,&
344 0 : & "SOLVE", transa="T")
345 : CALL cp_fm_cholesky_restore(work, k, preconditioner_env%fm, matrix_out,&
346 0 : & "SOLVE", transa="N")
347 0 : CALL cp_fm_release(work)
348 :
349 0 : CALL timestop(handle)
350 :
351 0 : END SUBROUTINE apply_full_direct
352 :
353 : ! **************************************************************************************************
354 : !> \brief Apply a Chebyshev approximation to the inverse of the stored dense SPD operator.
355 : !> \param preconditioner_env ...
356 : !> \param matrix_in ...
357 : !> \param matrix_out ...
358 : ! **************************************************************************************************
359 0 : SUBROUTINE apply_chebyshev_fm(preconditioner_env, matrix_in, matrix_out)
360 :
361 : TYPE(preconditioner_type) :: preconditioner_env
362 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
363 :
364 : INTEGER :: degree, iteration, k, n
365 : REAL(KIND=dp) :: delta, rho, rho_previous, sigma, theta
366 : TYPE(cp_fm_type) :: direction, residual
367 :
368 0 : CPASSERT(ASSOCIATED(preconditioner_env%fm))
369 0 : degree = preconditioner_env%polynomial_degree
370 0 : CPASSERT(degree >= 1)
371 0 : theta = 0.5_dp*(preconditioner_env%polynomial_max + preconditioner_env%polynomial_min)
372 0 : delta = 0.5_dp*(preconditioner_env%polynomial_max - preconditioner_env%polynomial_min)
373 0 : sigma = theta/delta
374 0 : rho_previous = 1.0_dp/sigma
375 0 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
376 0 : CALL cp_fm_create(direction, matrix_in%matrix_struct, name='Chebyshev direction')
377 0 : CALL cp_fm_create(residual, matrix_in%matrix_struct, name='Chebyshev residual')
378 0 : CALL cp_fm_to_fm(matrix_in, direction)
379 0 : CALL cp_fm_scale(1.0_dp/theta, direction)
380 0 : CALL cp_fm_to_fm(direction, matrix_out)
381 0 : DO iteration = 2, degree
382 0 : CALL cp_fm_to_fm(matrix_in, residual)
383 : CALL parallel_gemm('N', 'N', n, k, n, -1.0_dp, preconditioner_env%fm, &
384 0 : matrix_out, 1.0_dp, residual)
385 0 : rho = 1.0_dp/(2.0_dp*sigma - rho_previous)
386 0 : CALL cp_fm_scale(rho*rho_previous, direction)
387 0 : CALL cp_fm_scale_and_add(1.0_dp, direction, 2.0_dp*rho/delta, residual)
388 0 : CALL cp_fm_scale_and_add(1.0_dp, matrix_out, 1.0_dp, direction)
389 0 : rho_previous = rho
390 : END DO
391 0 : CALL cp_fm_release(residual)
392 0 : CALL cp_fm_release(direction)
393 :
394 0 : END SUBROUTINE apply_chebyshev_fm
395 :
396 : ! **************************************************************************************************
397 : !> \brief Apply a Chebyshev approximation to the inverse of the stored sparse SPD operator.
398 : !> \param preconditioner_env ...
399 : !> \param matrix_in ...
400 : !> \param matrix_out ...
401 : ! **************************************************************************************************
402 38 : SUBROUTINE apply_chebyshev_dbcsr(preconditioner_env, matrix_in, matrix_out)
403 :
404 : TYPE(preconditioner_type) :: preconditioner_env
405 : TYPE(dbcsr_type) :: matrix_in, matrix_out
406 :
407 : INTEGER :: degree, iteration
408 : REAL(KIND=dp) :: delta, rho, rho_previous, sigma, theta
409 : TYPE(dbcsr_type) :: direction, residual
410 :
411 0 : CPASSERT(ASSOCIATED(preconditioner_env%dbcsr_matrix))
412 38 : degree = preconditioner_env%polynomial_degree
413 38 : CPASSERT(degree >= 1)
414 38 : theta = 0.5_dp*(preconditioner_env%polynomial_max + preconditioner_env%polynomial_min)
415 38 : delta = 0.5_dp*(preconditioner_env%polynomial_max - preconditioner_env%polynomial_min)
416 38 : sigma = theta/delta
417 38 : rho_previous = 1.0_dp/sigma
418 38 : CALL dbcsr_copy(direction, matrix_in, name='Chebyshev direction')
419 38 : CALL dbcsr_copy(residual, matrix_in, name='Chebyshev residual')
420 38 : CALL dbcsr_set(matrix_out, 0.0_dp)
421 38 : CALL dbcsr_add(matrix_out, direction, 1.0_dp, 1.0_dp/theta)
422 38 : CALL dbcsr_set(direction, 0.0_dp)
423 38 : CALL dbcsr_add(direction, matrix_in, 1.0_dp, 1.0_dp/theta)
424 304 : DO iteration = 2, degree
425 266 : CALL dbcsr_set(residual, 0.0_dp)
426 266 : CALL dbcsr_add(residual, matrix_in, 1.0_dp, 1.0_dp)
427 : CALL dbcsr_multiply('N', 'N', -1.0_dp, preconditioner_env%dbcsr_matrix, &
428 266 : matrix_out, 1.0_dp, residual)
429 266 : rho = 1.0_dp/(2.0_dp*sigma - rho_previous)
430 266 : CALL dbcsr_add(direction, residual, rho*rho_previous, 2.0_dp*rho/delta)
431 266 : CALL dbcsr_add(matrix_out, direction, 1.0_dp, 1.0_dp)
432 304 : rho_previous = rho
433 : END DO
434 38 : CALL dbcsr_release(residual)
435 38 : CALL dbcsr_release(direction)
436 :
437 38 : END SUBROUTINE apply_chebyshev_dbcsr
438 :
439 : ! **************************************************************************************************
440 : !> \brief full all to a full matrix
441 : !> \param preconditioner_env ...
442 : !> \param matrix_in ...
443 : !> \param matrix_out ...
444 : ! **************************************************************************************************
445 113520 : SUBROUTINE apply_full_all(preconditioner_env, matrix_in, matrix_out)
446 :
447 : TYPE(preconditioner_type) :: preconditioner_env
448 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
449 :
450 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_all'
451 :
452 : INTEGER :: handle, i, j, k, n, ncol_local, &
453 : nrow_local
454 28380 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
455 : REAL(KIND=dp) :: dum
456 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
457 28380 : POINTER :: local_data
458 : TYPE(cp_fm_type) :: matrix_tmp
459 :
460 28380 : CALL timeset(routineN, handle)
461 :
462 28380 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
463 :
464 28380 : CALL cp_fm_create(matrix_tmp, matrix_in%matrix_struct, name="apply_full_all")
465 : CALL cp_fm_get_info(matrix_tmp, nrow_local=nrow_local, ncol_local=ncol_local, &
466 28380 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
467 :
468 : !
469 : CALL parallel_gemm('T', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
470 28380 : matrix_in, 0.0_dp, matrix_tmp)
471 :
472 : ! do the right scaling
473 224454 : DO j = 1, ncol_local
474 2822205 : DO i = 1, nrow_local
475 : dum = 1.0_dp/MAX(preconditioner_env%energy_gap, &
476 2597751 : preconditioner_env%full_evals(row_indices(i)) - preconditioner_env%occ_evals(col_indices(j)))
477 2793825 : local_data(i, j) = local_data(i, j)*dum
478 : END DO
479 : END DO
480 :
481 : ! mult back
482 : CALL parallel_gemm('N', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
483 28380 : matrix_tmp, 0.0_dp, matrix_out)
484 :
485 28380 : CALL cp_fm_release(matrix_tmp)
486 :
487 28380 : CALL timestop(handle)
488 :
489 28380 : END SUBROUTINE apply_full_all
490 :
491 : ! **************************************************************************************************
492 : !> \brief full all to a dbcsr matrix
493 : !> \param preconditioner_env ...
494 : !> \param matrix_in ...
495 : !> \param matrix_out ...
496 : ! **************************************************************************************************
497 47272 : SUBROUTINE apply_all(preconditioner_env, matrix_in, matrix_out)
498 :
499 : TYPE(preconditioner_type) :: preconditioner_env
500 : TYPE(dbcsr_type) :: matrix_in, matrix_out
501 :
502 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_all'
503 :
504 : INTEGER :: col, col_offset, col_size, handle, i, j, &
505 : row, row_offset, row_size
506 : REAL(KIND=dp) :: dum
507 23636 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: DATA
508 : TYPE(dbcsr_iterator_type) :: iter
509 : TYPE(dbcsr_type) :: matrix_tmp
510 :
511 23636 : CALL timeset(routineN, handle)
512 :
513 23636 : CALL dbcsr_copy(matrix_tmp, matrix_in, name="apply_full_all")
514 : CALL dbcsr_multiply('T', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
515 23636 : matrix_in, 0.0_dp, matrix_tmp)
516 : ! do the right scaling
517 23636 : CALL dbcsr_iterator_start(iter, matrix_tmp)
518 66616 : DO WHILE (dbcsr_iterator_blocks_left(iter))
519 : CALL dbcsr_iterator_next_block(iter, row, col, DATA, &
520 : row_size=row_size, col_size=col_size, &
521 42980 : row_offset=row_offset, col_offset=col_offset)
522 423944 : DO j = 1, col_size
523 3228737 : DO i = 1, row_size
524 : dum = 1.0_dp/MAX(preconditioner_env%energy_gap, &
525 : preconditioner_env%full_evals(row_offset + i - 1) &
526 2828429 : - preconditioner_env%occ_evals(col_offset + j - 1))
527 3185757 : DATA(i, j) = DATA(i, j)*dum
528 : END DO
529 : END DO
530 : END DO
531 23636 : CALL dbcsr_iterator_stop(iter)
532 :
533 : ! mult back
534 : CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
535 23636 : matrix_tmp, 0.0_dp, matrix_out)
536 23636 : CALL dbcsr_release(matrix_tmp)
537 23636 : CALL timestop(handle)
538 :
539 23636 : END SUBROUTINE apply_all
540 :
541 : END MODULE preconditioner_apply
|
