Line data Source code
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 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_env_type,&
17 : arnoldi_ev,&
18 : deallocate_arnoldi_env,&
19 : get_selected_ritz_val,&
20 : setup_arnoldi_env
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
27 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 : copy_fm_to_dbcsr
29 : USE cp_fm_basic_linalg, ONLY: cp_fm_uplo_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 9466 : 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 15240 : SELECT CASE (my_solver_type)
78 : CASE (ot_precond_solver_inv_chol)
79 : !
80 : ! compute the full inverse
81 5774 : preconditioner_env%solver = ot_precond_solver_inv_chol
82 5774 : 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) THEN
95 2 : CALL estimate_cond_num(preconditioner_env%sparse_matrix, preconditioner_env%condition_num)
96 : END IF
97 : CALL dbcsr_filter(preconditioner_env%sparse_matrix, &
98 6 : 1.0_dp/preconditioner_env%condition_num*0.01_dp)
99 6 : occ_matrix = dbcsr_get_occupation(preconditioner_env%sparse_matrix)
100 : END IF
101 : ! check whether we are in the first step and if it is a good idea to use cholesky (matrix sparsity)
102 6 : IF (preconditioner_env%solver /= ot_precond_solver_update .AND. occ_matrix > 0.5_dp) THEN
103 2 : preconditioner_env%solver = ot_precond_solver_update
104 2 : CALL make_full_inverse_cholesky(preconditioner_env, matrix_s)
105 : ELSE
106 4 : preconditioner_env%solver = ot_precond_solver_update
107 4 : CALL make_inverse_update(preconditioner_env, matrix_h)
108 : END IF
109 : CASE (ot_precond_solver_default)
110 3686 : preconditioner_env%solver = ot_precond_solver_default
111 : CASE DEFAULT
112 : !
113 9466 : CPABORT("Doesn't know this type of solver")
114 : END SELECT
115 :
116 9466 : END SUBROUTINE solve_preconditioner
117 :
118 : ! **************************************************************************************************
119 : !> \brief Compute the inverse using cholseky factorization
120 : !> \param preconditioner_env ...
121 : !> \param matrix_s ...
122 : ! **************************************************************************************************
123 17328 : SUBROUTINE make_full_inverse_cholesky(preconditioner_env, matrix_s)
124 :
125 : TYPE(preconditioner_type) :: preconditioner_env
126 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
127 :
128 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_inverse_cholesky'
129 :
130 : INTEGER :: handle, info
131 : TYPE(cp_fm_type) :: fm_work
132 : TYPE(cp_fm_type), POINTER :: fm
133 :
134 5776 : CALL timeset(routineN, handle)
135 :
136 : ! Maybe we will get a sparse Cholesky at a given point then this can go,
137 : ! if stuff was stored in fm anyway this simple returns
138 : CALL transfer_dbcsr_to_fm(preconditioner_env%sparse_matrix, preconditioner_env%fm, &
139 5776 : preconditioner_env%para_env, preconditioner_env%ctxt)
140 5776 : fm => preconditioner_env%fm
141 :
142 5776 : CALL cp_fm_create(fm_work, fm%matrix_struct, name="fm_work")
143 : !
144 : ! compute the inverse of SPD matrix fm using the Cholesky factorization
145 5776 : CALL cp_fm_cholesky_decompose(fm, info_out=info)
146 :
147 : !
148 : ! if fm not SPD we go with the overlap matrix
149 5776 : IF (info /= 0) THEN
150 : !
151 : ! just the overlap matrix
152 0 : IF (PRESENT(matrix_s)) THEN
153 0 : CALL copy_dbcsr_to_fm(matrix_s, fm)
154 0 : CALL cp_fm_cholesky_decompose(fm)
155 : ELSE
156 0 : CALL cp_fm_set_all(fm, alpha=0._dp, beta=1._dp)
157 : END IF
158 : END IF
159 5776 : CALL cp_fm_cholesky_invert(fm)
160 :
161 5776 : CALL cp_fm_uplo_to_full(fm, fm_work)
162 5776 : CALL cp_fm_release(fm_work)
163 :
164 5776 : CALL timestop(handle)
165 :
166 5776 : END SUBROUTINE make_full_inverse_cholesky
167 :
168 : ! **************************************************************************************************
169 : !> \brief Only perform the factorization, can be used later to solve the linear
170 : !> system on the fly
171 : !> \param preconditioner_env ...
172 : !> \param matrix_s ...
173 : ! **************************************************************************************************
174 0 : SUBROUTINE make_full_fact_cholesky(preconditioner_env, matrix_s)
175 :
176 : TYPE(preconditioner_type) :: preconditioner_env
177 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
178 :
179 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_fact_cholesky'
180 :
181 : INTEGER :: handle, info_out
182 : TYPE(cp_fm_type), POINTER :: fm
183 :
184 0 : CALL timeset(routineN, handle)
185 :
186 : ! Maybe we will get a sparse Cholesky at a given point then this can go,
187 : ! if stuff was stored in fm anyway this simple returns
188 : CALL transfer_dbcsr_to_fm(preconditioner_env%sparse_matrix, preconditioner_env%fm, &
189 0 : preconditioner_env%para_env, preconditioner_env%ctxt)
190 :
191 0 : fm => preconditioner_env%fm
192 : !
193 : ! compute the inverse of SPD matrix fm using the Cholesky factorization
194 0 : CALL cp_fm_cholesky_decompose(fm, info_out=info_out)
195 : !
196 : ! if fm not SPD we go with the overlap matrix
197 0 : IF (info_out /= 0) THEN
198 : !
199 : ! just the overlap matrix
200 0 : IF (PRESENT(matrix_s)) THEN
201 0 : CALL copy_dbcsr_to_fm(matrix_s, fm)
202 0 : CALL cp_fm_cholesky_decompose(fm)
203 : ELSE
204 0 : CALL cp_fm_set_all(fm, alpha=0._dp, beta=1._dp)
205 : END IF
206 : END IF
207 :
208 0 : CALL timestop(handle)
209 :
210 0 : END SUBROUTINE make_full_fact_cholesky
211 :
212 : ! **************************************************************************************************
213 : !> \brief computes an approximate inverse using Hotelling iterations
214 : !> \param preconditioner_env ...
215 : !> \param matrix_h as S is not always present this is a safe template for the transfer
216 : ! **************************************************************************************************
217 4 : SUBROUTINE make_inverse_update(preconditioner_env, matrix_h)
218 : TYPE(preconditioner_type) :: preconditioner_env
219 : TYPE(dbcsr_type), POINTER :: matrix_h
220 :
221 : CHARACTER(len=*), PARAMETER :: routineN = 'make_inverse_update'
222 :
223 : INTEGER :: handle
224 : LOGICAL :: use_guess
225 : REAL(KIND=dp) :: filter_eps
226 :
227 4 : CALL timeset(routineN, handle)
228 4 : use_guess = .TRUE.
229 : !
230 : ! uses an update of the full inverse (needs to be computed the first time)
231 : ! make sure preconditioner_env is not destroyed in between
232 :
233 4 : CALL cite_reference(Schiffmann2015)
234 :
235 : ! Maybe I gonna add a fm Hotelling, ... for now the same as above make sure we are dbcsr
236 4 : CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%sparse_matrix, matrix_h)
237 4 : IF (.NOT. ASSOCIATED(preconditioner_env%dbcsr_matrix)) THEN
238 0 : use_guess = .FALSE.
239 0 : CALL dbcsr_init_p(preconditioner_env%dbcsr_matrix)
240 : CALL dbcsr_create(preconditioner_env%dbcsr_matrix, "prec_dbcsr", &
241 0 : template=matrix_h, matrix_type=dbcsr_type_no_symmetry)
242 : END IF
243 :
244 : ! Try to get a reasonbale guess for the filtering threshold
245 4 : filter_eps = 1.0_dp/preconditioner_env%condition_num*0.1_dp
246 :
247 : ! Aggressive filtering on the initial guess is needed to avoid fill ins and retain sparsity
248 4 : CALL dbcsr_filter(preconditioner_env%dbcsr_matrix, filter_eps*100.0_dp)
249 : ! We don't need a high accuracy for the inverse so 0.4 is reasonable for convergence
250 : CALL invert_Hotelling(preconditioner_env%dbcsr_matrix, preconditioner_env%sparse_matrix, filter_eps*10.0_dp, &
251 4 : use_inv_as_guess=use_guess, norm_convergence=0.4_dp, filter_eps=filter_eps)
252 :
253 4 : CALL timestop(handle)
254 :
255 4 : END SUBROUTINE make_inverse_update
256 :
257 : ! **************************************************************************************************
258 : !> \brief computes an approximation to the condition number of a matrix using
259 : !> arnoldi iterations
260 : !> \param matrix ...
261 : !> \param cond_num ...
262 : ! **************************************************************************************************
263 2 : SUBROUTINE estimate_cond_num(matrix, cond_num)
264 : TYPE(dbcsr_type), POINTER :: matrix
265 : REAL(KIND=dp) :: cond_num
266 :
267 : CHARACTER(len=*), PARAMETER :: routineN = 'estimate_cond_num'
268 :
269 : INTEGER :: handle
270 : REAL(KIND=dp) :: max_ev, min_ev
271 : TYPE(arnoldi_env_type) :: arnoldi_env
272 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices
273 :
274 2 : CALL timeset(routineN, handle)
275 :
276 : ! its better to do 2 calculations as the maximum should quickly converge and the minimum won't need iram
277 4 : ALLOCATE (matrices(1))
278 2 : matrices(1)%matrix => matrix
279 : ! compute the minimum ev
280 : CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=5.0E-4_dp, selection_crit=2, &
281 2 : nval_request=1, nrestarts=15, generalized_ev=.FALSE., iram=.FALSE.)
282 2 : CALL arnoldi_ev(matrices, arnoldi_env)
283 2 : max_ev = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
284 2 : CALL deallocate_arnoldi_env(arnoldi_env)
285 :
286 : CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=5.0E-4_dp, selection_crit=3, &
287 2 : nval_request=1, nrestarts=15, generalized_ev=.FALSE., iram=.FALSE.)
288 2 : CALL arnoldi_ev(matrices, arnoldi_env)
289 2 : min_ev = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
290 2 : CALL deallocate_arnoldi_env(arnoldi_env)
291 :
292 2 : cond_num = max_ev/min_ev
293 2 : DEALLOCATE (matrices)
294 :
295 2 : CALL timestop(handle)
296 2 : END SUBROUTINE estimate_cond_num
297 :
298 : ! **************************************************************************************************
299 : !> \brief transfers a full matrix to a dbcsr
300 : !> \param fm_matrix a full matrix gets deallocated in the end
301 : !> \param dbcsr_matrix a dbcsr matrix, gets create from a template
302 : !> \param template_mat the template which is used for the structure
303 : ! **************************************************************************************************
304 7816 : SUBROUTINE transfer_fm_to_dbcsr(fm_matrix, dbcsr_matrix, template_mat)
305 :
306 : TYPE(cp_fm_type), POINTER :: fm_matrix
307 : TYPE(dbcsr_type), POINTER :: dbcsr_matrix, template_mat
308 :
309 : CHARACTER(len=*), PARAMETER :: routineN = 'transfer_fm_to_dbcsr'
310 :
311 : INTEGER :: handle
312 :
313 7816 : CALL timeset(routineN, handle)
314 7816 : IF (ASSOCIATED(fm_matrix)) THEN
315 7804 : IF (.NOT. ASSOCIATED(dbcsr_matrix)) THEN
316 4642 : CALL dbcsr_init_p(dbcsr_matrix)
317 : CALL dbcsr_create(dbcsr_matrix, template=template_mat, &
318 : name="preconditioner_env%dbcsr_matrix", &
319 4642 : matrix_type=dbcsr_type_no_symmetry)
320 : END IF
321 7804 : CALL copy_fm_to_dbcsr(fm_matrix, dbcsr_matrix)
322 7804 : CALL cp_fm_release(fm_matrix)
323 7804 : DEALLOCATE (fm_matrix)
324 : NULLIFY (fm_matrix)
325 : END IF
326 :
327 7816 : CALL timestop(handle)
328 :
329 7816 : 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 7434 : 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 7434 : CALL timeset(routineN, handle)
351 7434 : IF (ASSOCIATED(dbcsr_matrix)) THEN
352 5776 : NULLIFY (fm_struct_tmp)
353 :
354 5776 : IF (ASSOCIATED(fm_matrix)) THEN
355 0 : CALL cp_fm_release(fm_matrix)
356 0 : DEALLOCATE (fm_matrix)
357 : END IF
358 :
359 5776 : 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 5776 : context=context, para_env=para_env)
362 5776 : ALLOCATE (fm_matrix)
363 5776 : CALL cp_fm_create(fm_matrix, fm_struct_tmp)
364 5776 : CALL cp_fm_struct_release(fm_struct_tmp)
365 :
366 5776 : CALL copy_dbcsr_to_fm(dbcsr_matrix, fm_matrix)
367 5776 : CALL dbcsr_release(dbcsr_matrix)
368 5776 : DEALLOCATE (dbcsr_matrix)
369 : END IF
370 :
371 7434 : CALL timestop(handle)
372 :
373 7434 : END SUBROUTINE transfer_dbcsr_to_fm
374 :
375 : END MODULE preconditioner_solvers
|
