(git:34ef472)
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 ! **************************************************************************************************
15 MODULE 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 
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  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 
313 END MODULE preconditioner_apply
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