(git:e8f5963)
Loading...
Searching...
No Matches
preconditioner_apply.F
Go to the documentation of this file.
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! **************************************************************************************************
15MODULE 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
44CONTAINS
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 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 CALL timeset(routinen, handle)
62
63 SELECT CASE (preconditioner_env%in_use)
64 CASE (0)
65 cpabort("No preconditioner in use")
66 CASE (ot_precond_full_single)
67 CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
68 CASE (ot_precond_full_all)
69 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 SELECT CASE (preconditioner_env%solver)
72 CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
73 CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
74 CASE (ot_precond_solver_direct)
75 CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
76 CASE DEFAULT
77 cpabort("Solver not implemented")
78 END SELECT
79 CASE DEFAULT
80 cpabort("Unknown preconditioner")
81 END SELECT
82
83 CALL timestop(handle)
84
85 END SUBROUTINE apply_preconditioner_fm
86
87! **************************************************************************************************
88!> \brief ...
89!> \param preconditioner_env ...
90!> \param matrix_in ...
91!> \param matrix_out ...
92! **************************************************************************************************
93 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 CALL timeset(routinen, handle)
103
104 SELECT CASE (preconditioner_env%in_use)
105 CASE (0)
106 cpabort("No preconditioner in use")
107 CASE (ot_precond_full_single)
108 CALL apply_single(preconditioner_env, matrix_in, matrix_out)
109 CASE (ot_precond_full_all)
110 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 SELECT CASE (preconditioner_env%solver)
113 CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
114 CALL apply_single(preconditioner_env, matrix_in, matrix_out)
115 CASE (ot_precond_solver_direct)
116 cpabort("Apply_full_direct not supported with ot")
117 !CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
118 CASE DEFAULT
119 cpabort("Wrong solver")
120 END SELECT
121 CASE DEFAULT
122 cpabort("Wrong preconditioner")
123 END SELECT
124
125 CALL timestop(handle)
126
127 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 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 CALL timeset(routinen, handle)
145
146 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 matrix_in, 0.0_dp, matrix_out)
149 CALL timestop(handle)
150
151 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 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 CALL timeset(routinen, handle)
169
170 IF (.NOT. ASSOCIATED(preconditioner_env%dbcsr_matrix)) &
171 cpabort("NOT ASSOCIATED preconditioner_env%dbcsr_matrix")
172 CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, matrix_in, &
173 0.0_dp, matrix_out)
174
175 CALL timestop(handle)
176
177 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 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 CALL timeset(routinen, handle)
197
198 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 CALL cp_fm_cholesky_restore(matrix_in, k, preconditioner_env%fm, work,&
201 & "SOLVE", transa="T")
202 CALL cp_fm_cholesky_restore(work, k, preconditioner_env%fm, matrix_out,&
203 & "SOLVE", transa="N")
204 CALL cp_fm_release(work)
205
206 CALL timestop(handle)
207
208 END SUBROUTINE apply_full_direct
209
210! **************************************************************************************************
211!> \brief full all to a full matrix
212!> \param preconditioner_env ...
213!> \param matrix_in ...
214!> \param matrix_out ...
215! **************************************************************************************************
216 SUBROUTINE apply_full_all(preconditioner_env, matrix_in, matrix_out)
217
218 TYPE(preconditioner_type) :: preconditioner_env
219 TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
220
221 CHARACTER(len=*), PARAMETER :: routinen = 'apply_full_all'
222
223 INTEGER :: handle, i, j, k, n, ncol_local, &
224 nrow_local
225 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
226 REAL(kind=dp) :: dum
227 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
228 POINTER :: local_data
229 TYPE(cp_fm_type) :: matrix_tmp
230
231 CALL timeset(routinen, handle)
232
233 CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
234
235 CALL cp_fm_create(matrix_tmp, matrix_in%matrix_struct, name="apply_full_all")
236 CALL cp_fm_get_info(matrix_tmp, nrow_local=nrow_local, ncol_local=ncol_local, &
237 row_indices=row_indices, col_indices=col_indices, local_data=local_data)
238
239 !
240 CALL parallel_gemm('T', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
241 matrix_in, 0.0_dp, matrix_tmp)
242
243 ! do the right scaling
244 DO j = 1, ncol_local
245 DO i = 1, nrow_local
246 dum = 1.0_dp/max(preconditioner_env%energy_gap, &
247 preconditioner_env%full_evals(row_indices(i)) - preconditioner_env%occ_evals(col_indices(j)))
248 local_data(i, j) = local_data(i, j)*dum
249 END DO
250 END DO
251
252 ! mult back
253 CALL parallel_gemm('N', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
254 matrix_tmp, 0.0_dp, matrix_out)
255
256 CALL cp_fm_release(matrix_tmp)
257
258 CALL timestop(handle)
259
260 END SUBROUTINE apply_full_all
261
262! **************************************************************************************************
263!> \brief full all to a dbcsr matrix
264!> \param preconditioner_env ...
265!> \param matrix_in ...
266!> \param matrix_out ...
267! **************************************************************************************************
268 SUBROUTINE apply_all(preconditioner_env, matrix_in, matrix_out)
269
270 TYPE(preconditioner_type) :: preconditioner_env
271 TYPE(dbcsr_type) :: matrix_in, matrix_out
272
273 CHARACTER(len=*), PARAMETER :: routinen = 'apply_all'
274
275 INTEGER :: col, col_offset, col_size, handle, i, j, &
276 row, row_offset, row_size
277 REAL(kind=dp) :: dum
278 REAL(kind=dp), DIMENSION(:, :), POINTER :: DATA
279 TYPE(dbcsr_iterator_type) :: iter
280 TYPE(dbcsr_type) :: matrix_tmp
281
282 CALL timeset(routinen, handle)
283
284 CALL dbcsr_copy(matrix_tmp, matrix_in, name="apply_full_all")
285 CALL dbcsr_multiply('T', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
286 matrix_in, 0.0_dp, matrix_tmp)
287 ! do the right scaling
288 CALL dbcsr_iterator_start(iter, matrix_tmp)
289 DO WHILE (dbcsr_iterator_blocks_left(iter))
290 CALL dbcsr_iterator_next_block(iter, row, col, DATA, &
291 row_size=row_size, col_size=col_size, &
292 row_offset=row_offset, col_offset=col_offset)
293 DO j = 1, col_size
294 DO i = 1, row_size
295 dum = 1.0_dp/max(preconditioner_env%energy_gap, &
296 preconditioner_env%full_evals(row_offset + i - 1) &
297 - preconditioner_env%occ_evals(col_offset + j - 1))
298 DATA(i, j) = DATA(i, j)*dum
299 END DO
300 END DO
301 END DO
302 CALL dbcsr_iterator_stop(iter)
303
304 ! mult back
305 CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
306 matrix_tmp, 0.0_dp, matrix_out)
307 CALL dbcsr_release(matrix_tmp)
308 CALL timestop(handle)
309
310 END SUBROUTINE apply_all
311
312END MODULE preconditioner_apply
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:943
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1049
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1095
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:970
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1011
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1141
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa)
apply Cholesky decomposition op can be "SOLVE" (out = U^-1 * in) or "MULTIPLY" (out = U * in) pos can...
Definition cp_fm_cholesky.F:313
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1037
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::ot_precond_full_kinetic
integer, parameter, public ot_precond_full_kinetic
Definition input_constants.F:527
input_constants::ot_precond_full_single
integer, parameter, public ot_precond_full_single
Definition input_constants.F:527
input_constants::ot_precond_solver_inv_chol
integer, parameter, public ot_precond_solver_inv_chol
Definition input_constants.F:534
input_constants::ot_precond_full_single_inverse
integer, parameter, public ot_precond_full_single_inverse
Definition input_constants.F:527
input_constants::ot_precond_s_inverse
integer, parameter, public ot_precond_s_inverse
Definition input_constants.F:527
input_constants::ot_precond_solver_update
integer, parameter, public ot_precond_solver_update
Definition input_constants.F:534
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition input_constants.F:527
input_constants::ot_precond_solver_direct
integer, parameter, public ot_precond_solver_direct
Definition input_constants.F:534
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
preconditioner_apply
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition preconditioner_apply.F:15
preconditioner_apply::apply_preconditioner_fm
subroutine, public apply_preconditioner_fm(preconditioner_env, matrix_in, matrix_out)
applies a previously created preconditioner to a full matrix
Definition preconditioner_apply.F:53
preconditioner_apply::apply_preconditioner_dbcsr
subroutine, public apply_preconditioner_dbcsr(preconditioner_env, matrix_in, matrix_out)
...
Definition preconditioner_apply.F:94
preconditioner_types
types of preconditioners
Definition preconditioner_types.F:14
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
preconditioner_types::preconditioner_type
Definition preconditioner_types.F:42