Line data Source code
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 solves the preconditioner, contains to utility function for
10 : !> fm<->dbcsr transfers, should be moved soon
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_solvers
16 : USE arnoldi_api, ONLY: arnoldi_data_type,&
17 : arnoldi_ev,&
18 : deallocate_arnoldi_data,&
19 : get_selected_ritz_val,&
20 : setup_arnoldi_data
21 : USE bibliography, ONLY: Schiffmann2015,&
22 : cite_reference
23 : USE cp_blacs_env, ONLY: cp_blacs_env_type
24 : USE cp_dbcsr_api, ONLY: &
25 : dbcsr_create, dbcsr_filter, dbcsr_get_info, dbcsr_get_occupation, dbcsr_init_p, &
26 : dbcsr_p_type, dbcsr_release, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_real_8
27 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 : copy_fm_to_dbcsr
29 : USE cp_fm_basic_linalg, ONLY: cp_fm_upper_to_full
30 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
31 : cp_fm_cholesky_invert
32 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
33 : cp_fm_struct_release,&
34 : cp_fm_struct_type
35 : USE cp_fm_types, ONLY: cp_fm_create,&
36 : cp_fm_release,&
37 : cp_fm_set_all,&
38 : cp_fm_type
39 : USE input_constants, ONLY: ot_precond_solver_default,&
40 : ot_precond_solver_direct,&
41 : ot_precond_solver_inv_chol,&
42 : ot_precond_solver_update
43 : USE iterate_matrix, ONLY: invert_Hotelling
44 : USE kinds, ONLY: dp
45 : USE message_passing, ONLY: mp_para_env_type
46 : USE preconditioner_types, ONLY: preconditioner_type
47 : #include "./base/base_uses.f90"
48 :
49 : IMPLICIT NONE
50 :
51 : PRIVATE
52 :
53 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner_solvers'
54 :
55 : PUBLIC :: solve_preconditioner, transfer_fm_to_dbcsr, transfer_dbcsr_to_fm
56 :
57 : CONTAINS
58 :
59 : ! **************************************************************************************************
60 : !> \brief ...
61 : !> \param my_solver_type ...
62 : !> \param preconditioner_env ...
63 : !> \param matrix_s ...
64 : !> \param matrix_h ...
65 : ! **************************************************************************************************
66 8381 : SUBROUTINE solve_preconditioner(my_solver_type, preconditioner_env, matrix_s, &
67 : matrix_h)
68 : INTEGER :: my_solver_type
69 : TYPE(preconditioner_type) :: preconditioner_env
70 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
71 : TYPE(dbcsr_type), POINTER :: matrix_h
72 :
73 : REAL(dp) :: occ_matrix
74 :
75 : ! here comes the solver
76 :
77 13696 : SELECT CASE (my_solver_type)
78 : CASE (ot_precond_solver_inv_chol)
79 : !
80 : ! compute the full inverse
81 5315 : preconditioner_env%solver = ot_precond_solver_inv_chol
82 5315 : CALL make_full_inverse_cholesky(preconditioner_env, matrix_s)
83 : CASE (ot_precond_solver_direct)
84 : !
85 : ! prepare for the direct solver
86 0 : preconditioner_env%solver = ot_precond_solver_direct
87 0 : CALL make_full_fact_cholesky(preconditioner_env, matrix_s)
88 : CASE (ot_precond_solver_update)
89 : !
90 : ! uses an update of the full inverse (needs to be computed the first time)
91 : ! make sure preconditioner_env is not destroyed in between
92 6 : occ_matrix = 1.0_dp
93 6 : IF (ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
94 6 : IF (preconditioner_env%condition_num < 0.0_dp) &
95 2 : CALL estimate_cond_num(preconditioner_env%sparse_matrix, preconditioner_env%condition_num)
96 : CALL dbcsr_filter(preconditioner_env%sparse_matrix, &
97 6 : 1.0_dp/preconditioner_env%condition_num*0.01_dp)
98 6 : occ_matrix = dbcsr_get_occupation(preconditioner_env%sparse_matrix)
99 : END IF
100 : ! check whether we are in the first step and if it is a good idea to use cholesky (matrix sparsity)
101 6 : IF (preconditioner_env%solver .NE. ot_precond_solver_update .AND. occ_matrix > 0.5_dp) THEN
102 2 : preconditioner_env%solver = ot_precond_solver_update
103 2 : CALL make_full_inverse_cholesky(preconditioner_env, matrix_s)
104 : ELSE
105 4 : preconditioner_env%solver = ot_precond_solver_update
106 4 : CALL make_inverse_update(preconditioner_env, matrix_h)
107 : END IF
108 : CASE (ot_precond_solver_default)
109 3060 : preconditioner_env%solver = ot_precond_solver_default
110 : CASE DEFAULT
111 : !
112 8381 : CPABORT("Doesn't know this type of solver")
113 : END SELECT
114 :
115 8381 : END SUBROUTINE solve_preconditioner
116 :
117 : ! **************************************************************************************************
118 : !> \brief Compute the inverse using cholseky factorization
119 : !> \param preconditioner_env ...
120 : !> \param matrix_s ...
121 : ! **************************************************************************************************
122 15951 : SUBROUTINE make_full_inverse_cholesky(preconditioner_env, matrix_s)
123 :
124 : TYPE(preconditioner_type) :: preconditioner_env
125 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
126 :
127 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_inverse_cholesky'
128 :
129 : INTEGER :: handle, info
130 : TYPE(cp_fm_type) :: fm_work
131 : TYPE(cp_fm_type), POINTER :: fm
132 :
133 5317 : CALL timeset(routineN, handle)
134 :
135 : ! Maybe we will get a sparse Cholesky at a given point then this can go,
136 : ! if stuff was stored in fm anyway this simple returns
137 : CALL transfer_dbcsr_to_fm(preconditioner_env%sparse_matrix, preconditioner_env%fm, &
138 5317 : preconditioner_env%para_env, preconditioner_env%ctxt)
139 5317 : fm => preconditioner_env%fm
140 :
141 5317 : CALL cp_fm_create(fm_work, fm%matrix_struct, name="fm_work")
142 : !
143 : ! compute the inverse of SPD matrix fm using the Cholesky factorization
144 5317 : CALL cp_fm_cholesky_decompose(fm, info_out=info)
145 :
146 : !
147 : ! if fm not SPD we go with the overlap matrix
148 5317 : IF (info /= 0) THEN
149 : !
150 : ! just the overlap matrix
151 0 : IF (PRESENT(matrix_s)) THEN
152 0 : CALL copy_dbcsr_to_fm(matrix_s, fm)
153 0 : CALL cp_fm_cholesky_decompose(fm)
154 : ELSE
155 0 : CALL cp_fm_set_all(fm, alpha=0._dp, beta=1._dp)
156 : END IF
157 : END IF
158 5317 : CALL cp_fm_cholesky_invert(fm)
159 :
160 5317 : CALL cp_fm_upper_to_full(fm, fm_work)
161 5317 : CALL cp_fm_release(fm_work)
162 :
163 5317 : CALL timestop(handle)
164 :
165 5317 : END SUBROUTINE make_full_inverse_cholesky
166 :
167 : ! **************************************************************************************************
168 : !> \brief Only perform the factorization, can be used later to solve the linear
169 : !> system on the fly
170 : !> \param preconditioner_env ...
171 : !> \param matrix_s ...
172 : ! **************************************************************************************************
173 0 : SUBROUTINE make_full_fact_cholesky(preconditioner_env, matrix_s)
174 :
175 : TYPE(preconditioner_type) :: preconditioner_env
176 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
177 :
178 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_fact_cholesky'
179 :
180 : INTEGER :: handle, info_out
181 : TYPE(cp_fm_type), POINTER :: fm
182 :
183 0 : CALL timeset(routineN, handle)
184 :
185 : ! Maybe we will get a sparse Cholesky at a given point then this can go,
186 : ! if stuff was stored in fm anyway this simple returns
187 : CALL transfer_dbcsr_to_fm(preconditioner_env%sparse_matrix, preconditioner_env%fm, &
188 0 : preconditioner_env%para_env, preconditioner_env%ctxt)
189 :
190 0 : fm => preconditioner_env%fm
191 : !
192 : ! compute the inverse of SPD matrix fm using the Cholesky factorization
193 0 : CALL cp_fm_cholesky_decompose(fm, info_out=info_out)
194 : !
195 : ! if fm not SPD we go with the overlap matrix
196 0 : IF (info_out .NE. 0) THEN
197 : !
198 : ! just the overlap matrix
199 0 : IF (PRESENT(matrix_s)) THEN
200 0 : CALL copy_dbcsr_to_fm(matrix_s, fm)
201 0 : CALL cp_fm_cholesky_decompose(fm)
202 : ELSE
203 0 : CALL cp_fm_set_all(fm, alpha=0._dp, beta=1._dp)
204 : END IF
205 : END IF
206 :
207 0 : CALL timestop(handle)
208 :
209 0 : END SUBROUTINE make_full_fact_cholesky
210 :
211 : ! **************************************************************************************************
212 : !> \brief computes an approximate inverse using Hotelling iterations
213 : !> \param preconditioner_env ...
214 : !> \param matrix_h as S is not always present this is a safe template for the transfer
215 : ! **************************************************************************************************
216 4 : SUBROUTINE make_inverse_update(preconditioner_env, matrix_h)
217 : TYPE(preconditioner_type) :: preconditioner_env
218 : TYPE(dbcsr_type), POINTER :: matrix_h
219 :
220 : CHARACTER(len=*), PARAMETER :: routineN = 'make_inverse_update'
221 :
222 : INTEGER :: handle
223 : LOGICAL :: use_guess
224 : REAL(KIND=dp) :: filter_eps
225 :
226 4 : CALL timeset(routineN, handle)
227 4 : use_guess = .TRUE.
228 : !
229 : ! uses an update of the full inverse (needs to be computed the first time)
230 : ! make sure preconditioner_env is not destroyed in between
231 :
232 4 : CALL cite_reference(Schiffmann2015)
233 :
234 : ! Maybe I gonna add a fm Hotelling, ... for now the same as above make sure we are dbcsr
235 4 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%sparse_matrix, matrix_h)
236 4 : IF (.NOT. ASSOCIATED(preconditioner_env%dbcsr_matrix)) THEN
237 0 : use_guess = .FALSE.
238 0 : CALL dbcsr_init_p(preconditioner_env%dbcsr_matrix)
239 : CALL dbcsr_create(preconditioner_env%dbcsr_matrix, "prec_dbcsr", &
240 0 : template=matrix_h, matrix_type=dbcsr_type_no_symmetry)
241 : END IF
242 :
243 : ! Try to get a reasonbale guess for the filtering threshold
244 4 : filter_eps = 1.0_dp/preconditioner_env%condition_num*0.1_dp
245 :
246 : ! Aggressive filtering on the initial guess is needed to avoid fill ins and retain sparsity
247 4 : CALL dbcsr_filter(preconditioner_env%dbcsr_matrix, filter_eps*100.0_dp)
248 : ! We don't need a high accuracy for the inverse so 0.4 is reasonable for convergence
249 : CALL invert_Hotelling(preconditioner_env%dbcsr_matrix, preconditioner_env%sparse_matrix, filter_eps*10.0_dp, &
250 4 : use_inv_as_guess=use_guess, norm_convergence=0.4_dp, filter_eps=filter_eps)
251 :
252 4 : CALL timestop(handle)
253 :
254 4 : END SUBROUTINE make_inverse_update
255 :
256 : ! **************************************************************************************************
257 : !> \brief computes an approximation to the condition number of a matrix using
258 : !> arnoldi iterations
259 : !> \param matrix ...
260 : !> \param cond_num ...
261 : ! **************************************************************************************************
262 2 : SUBROUTINE estimate_cond_num(matrix, cond_num)
263 : TYPE(dbcsr_type), POINTER :: matrix
264 : REAL(KIND=dp) :: cond_num
265 :
266 : CHARACTER(len=*), PARAMETER :: routineN = 'estimate_cond_num'
267 :
268 : INTEGER :: handle
269 : REAL(KIND=dp) :: max_ev, min_ev
270 : TYPE(arnoldi_data_type) :: my_arnoldi
271 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices
272 :
273 2 : CALL timeset(routineN, handle)
274 :
275 : ! its better to do 2 calculations as the maximum should quickly converge and the minimum won't need iram
276 4 : ALLOCATE (matrices(1))
277 2 : matrices(1)%matrix => matrix
278 : ! compute the minimum ev
279 : CALL setup_arnoldi_data(my_arnoldi, matrices, max_iter=20, threshold=5.0E-4_dp, selection_crit=2, &
280 2 : nval_request=1, nrestarts=15, generalized_ev=.FALSE., iram=.FALSE.)
281 2 : CALL arnoldi_ev(matrices, my_arnoldi)
282 2 : max_ev = REAL(get_selected_ritz_val(my_arnoldi, 1), dp)
283 2 : CALL deallocate_arnoldi_data(my_arnoldi)
284 :
285 : CALL setup_arnoldi_data(my_arnoldi, matrices, max_iter=20, threshold=5.0E-4_dp, selection_crit=3, &
286 2 : nval_request=1, nrestarts=15, generalized_ev=.FALSE., iram=.FALSE.)
287 2 : CALL arnoldi_ev(matrices, my_arnoldi)
288 2 : min_ev = REAL(get_selected_ritz_val(my_arnoldi, 1), dp)
289 2 : CALL deallocate_arnoldi_data(my_arnoldi)
290 :
291 2 : cond_num = max_ev/min_ev
292 2 : DEALLOCATE (matrices)
293 :
294 2 : CALL timestop(handle)
295 2 : END SUBROUTINE estimate_cond_num
296 :
297 : ! **************************************************************************************************
298 : !> \brief transfers a full matrix to a dbcsr
299 : !> \param fm_matrix a full matrix gets deallocated in the end
300 : !> \param dbcsr_matrix a dbcsr matrix, gets create from a template
301 : !> \param template_mat the template which is used for the structure
302 : ! **************************************************************************************************
303 6951 : SUBROUTINE transfer_fm_to_dbcsr(fm_matrix, dbcsr_matrix, template_mat)
304 :
305 : TYPE(cp_fm_type), POINTER :: fm_matrix
306 : TYPE(dbcsr_type), POINTER :: dbcsr_matrix, template_mat
307 :
308 : CHARACTER(len=*), PARAMETER :: routineN = 'transfer_fm_to_dbcsr'
309 :
310 : INTEGER :: handle
311 :
312 6951 : CALL timeset(routineN, handle)
313 6951 : IF (ASSOCIATED(fm_matrix)) THEN
314 6939 : IF (.NOT. ASSOCIATED(dbcsr_matrix)) THEN
315 3964 : CALL dbcsr_init_p(dbcsr_matrix)
316 : CALL dbcsr_create(dbcsr_matrix, template=template_mat, &
317 : name="preconditioner_env%dbcsr_matrix", &
318 : matrix_type=dbcsr_type_no_symmetry, &
319 3964 : nze=0, data_type=dbcsr_type_real_8)
320 : END IF
321 6939 : CALL copy_fm_to_dbcsr(fm_matrix, dbcsr_matrix)
322 6939 : CALL cp_fm_release(fm_matrix)
323 6939 : DEALLOCATE (fm_matrix)
324 : NULLIFY (fm_matrix)
325 : END IF
326 :
327 6951 : CALL timestop(handle)
328 :
329 6951 : END SUBROUTINE transfer_fm_to_dbcsr
330 :
331 : ! **************************************************************************************************
332 : !> \brief transfers a dbcsr to a full matrix
333 : !> \param dbcsr_matrix a dbcsr matrix, gets deallocated at the end
334 : !> \param fm_matrix a full matrix gets created if not yet done
335 : !> \param para_env the para_env
336 : !> \param context the blacs context
337 : ! **************************************************************************************************
338 6755 : SUBROUTINE transfer_dbcsr_to_fm(dbcsr_matrix, fm_matrix, para_env, context)
339 :
340 : TYPE(dbcsr_type), POINTER :: dbcsr_matrix
341 : TYPE(cp_fm_type), POINTER :: fm_matrix
342 : TYPE(mp_para_env_type), POINTER :: para_env
343 : TYPE(cp_blacs_env_type), POINTER :: context
344 :
345 : CHARACTER(len=*), PARAMETER :: routineN = 'transfer_dbcsr_to_fm'
346 :
347 : INTEGER :: handle, n
348 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
349 :
350 6755 : CALL timeset(routineN, handle)
351 6755 : IF (ASSOCIATED(dbcsr_matrix)) THEN
352 5317 : NULLIFY (fm_struct_tmp)
353 :
354 5317 : IF (ASSOCIATED(fm_matrix)) THEN
355 0 : CALL cp_fm_release(fm_matrix)
356 0 : DEALLOCATE (fm_matrix)
357 : END IF
358 :
359 5317 : CALL dbcsr_get_info(dbcsr_matrix, nfullrows_total=n)
360 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
361 5317 : context=context, para_env=para_env)
362 5317 : ALLOCATE (fm_matrix)
363 5317 : CALL cp_fm_create(fm_matrix, fm_struct_tmp)
364 5317 : CALL cp_fm_struct_release(fm_struct_tmp)
365 :
366 5317 : CALL copy_dbcsr_to_fm(dbcsr_matrix, fm_matrix)
367 5317 : CALL dbcsr_release(dbcsr_matrix)
368 5317 : DEALLOCATE (dbcsr_matrix)
369 : END IF
370 :
371 6755 : CALL timestop(handle)
372 :
373 6755 : END SUBROUTINE transfer_dbcsr_to_fm
374 :
375 : END MODULE preconditioner_solvers
|
