Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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_dbcsr_api, ONLY: &
17 : dbcsr_copy, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
18 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_release, dbcsr_type
19 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_restore
20 : USE cp_fm_types, ONLY: cp_fm_create,&
21 : cp_fm_get_info,&
22 : cp_fm_release,&
23 : cp_fm_type
24 : USE input_constants, ONLY: ot_precond_full_all,&
25 : ot_precond_full_kinetic,&
26 : ot_precond_full_single,&
27 : ot_precond_full_single_inverse,&
28 : ot_precond_s_inverse,&
29 : ot_precond_solver_direct,&
30 : ot_precond_solver_inv_chol,&
31 : ot_precond_solver_update
32 : USE kinds, ONLY: dp
33 : USE parallel_gemm_api, ONLY: parallel_gemm
34 : USE preconditioner_types, ONLY: preconditioner_type
35 : #include "./base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 : PRIVATE
39 :
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner_apply'
41 :
42 : PUBLIC :: apply_preconditioner_fm, apply_preconditioner_dbcsr
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief applies a previously created preconditioner to a full matrix
48 : !> \param preconditioner_env ...
49 : !> \param matrix_in ...
50 : !> \param matrix_out ...
51 : ! **************************************************************************************************
52 47856 : SUBROUTINE apply_preconditioner_fm(preconditioner_env, matrix_in, matrix_out)
53 :
54 : TYPE(preconditioner_type) :: preconditioner_env
55 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
56 :
57 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_fm'
58 :
59 : INTEGER :: handle
60 :
61 47856 : CALL timeset(routineN, handle)
62 :
63 47856 : SELECT CASE (preconditioner_env%in_use)
64 : CASE (0)
65 0 : CPABORT("No preconditioner in use")
66 : CASE (ot_precond_full_single)
67 1320 : CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
68 : CASE (ot_precond_full_all)
69 27280 : CALL apply_full_all(preconditioner_env, matrix_in, matrix_out)
70 : CASE (ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
71 38512 : SELECT CASE (preconditioner_env%solver)
72 : CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
73 19256 : CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
74 : CASE (ot_precond_solver_direct)
75 0 : CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
76 : CASE DEFAULT
77 19256 : CPABORT("Solver not implemented")
78 : END SELECT
79 : CASE DEFAULT
80 47856 : CPABORT("Unknown preconditioner")
81 : END SELECT
82 :
83 47856 : CALL timestop(handle)
84 :
85 47856 : END SUBROUTINE apply_preconditioner_fm
86 :
87 : ! **************************************************************************************************
88 : !> \brief ...
89 : !> \param preconditioner_env ...
90 : !> \param matrix_in ...
91 : !> \param matrix_out ...
92 : ! **************************************************************************************************
93 65489 : SUBROUTINE apply_preconditioner_dbcsr(preconditioner_env, matrix_in, matrix_out)
94 :
95 : TYPE(preconditioner_type) :: preconditioner_env
96 : TYPE(dbcsr_type) :: matrix_in, matrix_out
97 :
98 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_dbcsr'
99 :
100 : INTEGER :: handle
101 :
102 65489 : CALL timeset(routineN, handle)
103 :
104 65489 : SELECT CASE (preconditioner_env%in_use)
105 : CASE (0)
106 0 : CPABORT("No preconditioner in use")
107 : CASE (ot_precond_full_single)
108 208 : CALL apply_single(preconditioner_env, matrix_in, matrix_out)
109 : CASE (ot_precond_full_all)
110 17330 : CALL apply_all(preconditioner_env, matrix_in, matrix_out)
111 : CASE (ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
112 95902 : SELECT CASE (preconditioner_env%solver)
113 : CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
114 47951 : CALL apply_single(preconditioner_env, matrix_in, matrix_out)
115 : CASE (ot_precond_solver_direct)
116 0 : CPABORT("Apply_full_direct not supported with ot")
117 : !CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
118 : CASE DEFAULT
119 47951 : CPABORT("Wrong solver")
120 : END SELECT
121 : CASE DEFAULT
122 65489 : CPABORT("Wrong preconditioner")
123 : END SELECT
124 :
125 65489 : CALL timestop(handle)
126 :
127 65489 : END SUBROUTINE apply_preconditioner_dbcsr
128 :
129 : ! **************************************************************************************************
130 : !> \brief apply to full matrix, complete inversion has already been done
131 : !> \param preconditioner_env ...
132 : !> \param matrix_in ...
133 : !> \param matrix_out ...
134 : ! **************************************************************************************************
135 41152 : SUBROUTINE apply_full_single(preconditioner_env, matrix_in, matrix_out)
136 :
137 : TYPE(preconditioner_type) :: preconditioner_env
138 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
139 :
140 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_single'
141 :
142 : INTEGER :: handle, k, n
143 :
144 20576 : CALL timeset(routineN, handle)
145 :
146 20576 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
147 : CALL parallel_gemm('N', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
148 20576 : matrix_in, 0.0_dp, matrix_out)
149 20576 : CALL timestop(handle)
150 :
151 20576 : END SUBROUTINE apply_full_single
152 :
153 : ! **************************************************************************************************
154 : !> \brief apply to dbcsr matrix, complete inversion has already been done
155 : !> \param preconditioner_env ...
156 : !> \param matrix_in ...
157 : !> \param matrix_out ...
158 : ! **************************************************************************************************
159 48159 : SUBROUTINE apply_single(preconditioner_env, matrix_in, matrix_out)
160 :
161 : TYPE(preconditioner_type) :: preconditioner_env
162 : TYPE(dbcsr_type) :: matrix_in, matrix_out
163 :
164 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_single'
165 :
166 : INTEGER :: handle
167 :
168 48159 : CALL timeset(routineN, handle)
169 :
170 48159 : IF (.NOT. ASSOCIATED(preconditioner_env%dbcsr_matrix)) &
171 0 : CPABORT("NOT ASSOCIATED preconditioner_env%dbcsr_matrix")
172 : CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, matrix_in, &
173 48159 : 0.0_dp, matrix_out)
174 :
175 48159 : CALL timestop(handle)
176 :
177 48159 : END SUBROUTINE apply_single
178 :
179 : ! **************************************************************************************************
180 : !> \brief preconditioner contains the factorization, application done by
181 : !> solving the linear system
182 : !> \param preconditioner_env ...
183 : !> \param matrix_in ...
184 : !> \param matrix_out ...
185 : ! **************************************************************************************************
186 0 : SUBROUTINE apply_full_direct(preconditioner_env, matrix_in, matrix_out)
187 :
188 : TYPE(preconditioner_type) :: preconditioner_env
189 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
190 :
191 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_direct'
192 :
193 : INTEGER :: handle, k, n
194 : TYPE(cp_fm_type) :: work
195 :
196 0 : CALL timeset(routineN, handle)
197 :
198 0 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
199 : CALL cp_fm_create(work, matrix_in%matrix_struct, name="apply_full_single", &
200 0 : use_sp=matrix_in%use_sp)
201 : CALL cp_fm_cholesky_restore(matrix_in, k, preconditioner_env%fm, work,&
202 0 : & "SOLVE", transa="T")
203 : CALL cp_fm_cholesky_restore(work, k, preconditioner_env%fm, matrix_out,&
204 0 : & "SOLVE", transa="N")
205 0 : CALL cp_fm_release(work)
206 :
207 0 : CALL timestop(handle)
208 :
209 0 : END SUBROUTINE apply_full_direct
210 :
211 : ! **************************************************************************************************
212 : !> \brief full all to a full matrix
213 : !> \param preconditioner_env ...
214 : !> \param matrix_in ...
215 : !> \param matrix_out ...
216 : ! **************************************************************************************************
217 109120 : SUBROUTINE apply_full_all(preconditioner_env, matrix_in, matrix_out)
218 :
219 : TYPE(preconditioner_type) :: preconditioner_env
220 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
221 :
222 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_all'
223 :
224 : INTEGER :: handle, i, j, k, n, ncol_local, &
225 : nrow_local
226 27280 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
227 : REAL(KIND=dp) :: dum
228 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
229 27280 : POINTER :: local_data
230 : TYPE(cp_fm_type) :: matrix_tmp
231 :
232 27280 : CALL timeset(routineN, handle)
233 :
234 27280 : CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
235 :
236 27280 : CALL cp_fm_create(matrix_tmp, matrix_in%matrix_struct, name="apply_full_all")
237 : CALL cp_fm_get_info(matrix_tmp, nrow_local=nrow_local, ncol_local=ncol_local, &
238 27280 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
239 :
240 : !
241 : CALL parallel_gemm('T', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
242 27280 : matrix_in, 0.0_dp, matrix_tmp)
243 :
244 : ! do the right scaling
245 217266 : DO j = 1, ncol_local
246 2708530 : DO i = 1, nrow_local
247 : dum = 1.0_dp/MAX(preconditioner_env%energy_gap, &
248 2491264 : preconditioner_env%full_evals(row_indices(i)) - preconditioner_env%occ_evals(col_indices(j)))
249 2681250 : local_data(i, j) = local_data(i, j)*dum
250 : END DO
251 : END DO
252 :
253 : ! mult back
254 : CALL parallel_gemm('N', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
255 27280 : matrix_tmp, 0.0_dp, matrix_out)
256 :
257 27280 : CALL cp_fm_release(matrix_tmp)
258 :
259 27280 : CALL timestop(handle)
260 :
261 27280 : END SUBROUTINE apply_full_all
262 :
263 : ! **************************************************************************************************
264 : !> \brief full all to a dbcsr matrix
265 : !> \param preconditioner_env ...
266 : !> \param matrix_in ...
267 : !> \param matrix_out ...
268 : ! **************************************************************************************************
269 34660 : SUBROUTINE apply_all(preconditioner_env, matrix_in, matrix_out)
270 :
271 : TYPE(preconditioner_type) :: preconditioner_env
272 : TYPE(dbcsr_type) :: matrix_in, matrix_out
273 :
274 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_all'
275 :
276 : INTEGER :: col, col_offset, col_size, handle, i, j, &
277 : row, row_offset, row_size
278 : REAL(KIND=dp) :: dum
279 17330 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: DATA
280 : TYPE(dbcsr_iterator_type) :: iter
281 : TYPE(dbcsr_type) :: matrix_tmp
282 :
283 17330 : CALL timeset(routineN, handle)
284 :
285 17330 : CALL dbcsr_copy(matrix_tmp, matrix_in, name="apply_full_all")
286 : CALL dbcsr_multiply('T', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
287 17330 : matrix_in, 0.0_dp, matrix_tmp)
288 : ! do the right scaling
289 17330 : CALL dbcsr_iterator_start(iter, matrix_tmp)
290 49511 : DO WHILE (dbcsr_iterator_blocks_left(iter))
291 : CALL dbcsr_iterator_next_block(iter, row, col, DATA, &
292 : row_size=row_size, col_size=col_size, &
293 32181 : row_offset=row_offset, col_offset=col_offset)
294 346719 : DO j = 1, col_size
295 2838974 : DO i = 1, row_size
296 : dum = 1.0_dp/MAX(preconditioner_env%energy_gap, &
297 : preconditioner_env%full_evals(row_offset + i - 1) &
298 2509585 : - preconditioner_env%occ_evals(col_offset + j - 1))
299 2806793 : DATA(i, j) = DATA(i, j)*dum
300 : END DO
301 : END DO
302 : END DO
303 17330 : CALL dbcsr_iterator_stop(iter)
304 :
305 : ! mult back
306 : CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
307 17330 : matrix_tmp, 0.0_dp, matrix_out)
308 17330 : CALL dbcsr_release(matrix_tmp)
309 17330 : CALL timestop(handle)
310 :
311 17330 : END SUBROUTINE apply_all
312 :
313 : END MODULE preconditioner_apply
|
