(git:374b731)
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-2024 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_fm_cholesky, ONLY: cp_fm_cholesky_restore
17 USE cp_fm_types, ONLY: cp_fm_create,&
18 cp_fm_get_info,&
19 cp_fm_release,&
20 cp_fm_type
21 USE dbcsr_api, ONLY: &
22 dbcsr_copy, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
23 dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_release, dbcsr_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 use_sp=matrix_in%use_sp)
201 CALL cp_fm_cholesky_restore(matrix_in, k, preconditioner_env%fm, work,&
202 & "SOLVE", transa="T")
203 CALL cp_fm_cholesky_restore(work, k, preconditioner_env%fm, matrix_out,&
204 & "SOLVE", transa="N")
205 CALL cp_fm_release(work)
206
207 CALL timestop(handle)
208
209 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
227 REAL(kind=dp) :: dum
228 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
229 POINTER :: local_data
230 TYPE(cp_fm_type) :: matrix_tmp
231
232 CALL timeset(routinen, handle)
233
234 CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
235
236 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 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 matrix_in, 0.0_dp, matrix_tmp)
243
244 ! do the right scaling
245 DO j = 1, ncol_local
246 DO i = 1, nrow_local
247 dum = 1.0_dp/max(preconditioner_env%energy_gap, &
248 preconditioner_env%full_evals(row_indices(i)) - preconditioner_env%occ_evals(col_indices(j)))
249 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 matrix_tmp, 0.0_dp, matrix_out)
256
257 CALL cp_fm_release(matrix_tmp)
258
259 CALL timestop(handle)
260
261 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: DATA
280 TYPE(dbcsr_iterator_type) :: iter
281 TYPE(dbcsr_type) :: matrix_tmp
282
283 CALL timeset(routinen, handle)
284
285 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 matrix_in, 0.0_dp, matrix_tmp)
288 ! do the right scaling
289 CALL dbcsr_iterator_start(iter, matrix_tmp)
290 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 row_offset=row_offset, col_offset=col_offset)
294 DO j = 1, col_size
295 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 - preconditioner_env%occ_evals(col_offset + j - 1))
299 DATA(i, j) = DATA(i, j)*dum
300 END DO
301 END DO
302 END DO
303 CALL dbcsr_iterator_stop(iter)
304
305 ! mult back
306 CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
307 matrix_tmp, 0.0_dp, matrix_out)
308 CALL dbcsr_release(matrix_tmp)
309 CALL timestop(handle)
310
311 END SUBROUTINE apply_all
312
313END MODULE preconditioner_apply
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
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(matrix, neig, matrixb, matrixout, op, pos, transa)
...
Definition cp_fm_cholesky.F:254
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:1016
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
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:513
input_constants::ot_precond_full_single
integer, parameter, public ot_precond_full_single
Definition input_constants.F:513
input_constants::ot_precond_solver_inv_chol
integer, parameter, public ot_precond_solver_inv_chol
Definition input_constants.F:520
input_constants::ot_precond_full_single_inverse
integer, parameter, public ot_precond_full_single_inverse
Definition input_constants.F:513
input_constants::ot_precond_s_inverse
integer, parameter, public ot_precond_s_inverse
Definition input_constants.F:513
input_constants::ot_precond_solver_update
integer, parameter, public ot_precond_solver_update
Definition input_constants.F:520
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition input_constants.F:513
input_constants::ot_precond_solver_direct
integer, parameter, public ot_precond_solver_direct
Definition input_constants.F:520
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_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
preconditioner_types::preconditioner_type
Definition preconditioner_types.F:42