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