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 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_makes
16 : USE arnoldi_api, ONLY: arnoldi_env_type,&
17 : arnoldi_ev,&
18 : deallocate_arnoldi_env,&
19 : get_selected_ritz_val,&
20 : get_selected_ritz_vector,&
21 : set_arnoldi_initial_vector,&
22 : setup_arnoldi_env
23 : USE cp_dbcsr_api, ONLY: &
24 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_multiply, dbcsr_p_type, &
25 : dbcsr_release, dbcsr_type, dbcsr_type_symmetric
26 : USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag
27 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 : cp_dbcsr_m_by_n_from_template,&
29 : cp_dbcsr_sm_fm_multiply,&
30 : cp_fm_to_dbcsr_row_template
31 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
32 : cp_fm_triangular_invert,&
33 : cp_fm_triangular_multiply,&
34 : cp_fm_uplo_to_full
35 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
36 : cp_fm_cholesky_reduce,&
37 : cp_fm_cholesky_restore
38 : USE cp_fm_diag, ONLY: choose_eigv_solver
39 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
40 : cp_fm_struct_release,&
41 : cp_fm_struct_type
42 : USE cp_fm_types, ONLY: cp_fm_create,&
43 : cp_fm_get_diag,&
44 : cp_fm_get_info,&
45 : cp_fm_release,&
46 : cp_fm_to_fm,&
47 : cp_fm_type
48 : USE input_constants, ONLY: &
49 : cholesky_inverse, cholesky_reduce, ot_precond_full_all, ot_precond_full_kinetic, &
50 : ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_s_inverse, &
51 : ot_precond_solver_default, ot_precond_solver_inv_chol
52 : USE kinds, ONLY: dp
53 : USE parallel_gemm_api, ONLY: parallel_gemm
54 : USE preconditioner_types, ONLY: preconditioner_type
55 : #include "./base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 :
59 : PRIVATE
60 :
61 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner_makes'
62 :
63 : PUBLIC :: make_preconditioner_matrix
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief ...
69 : !> \param preconditioner_env ...
70 : !> \param matrix_h ...
71 : !> \param matrix_s ...
72 : !> \param matrix_t ...
73 : !> \param mo_coeff ...
74 : !> \param energy_homo ...
75 : !> \param eigenvalues_ot ...
76 : !> \param energy_gap ...
77 : !> \param my_solver_type ...
78 : ! **************************************************************************************************
79 9466 : SUBROUTINE make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
80 9466 : energy_homo, eigenvalues_ot, energy_gap, &
81 : my_solver_type)
82 : TYPE(preconditioner_type) :: preconditioner_env
83 : TYPE(dbcsr_type), POINTER :: matrix_h
84 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s, matrix_t
85 : TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
86 : REAL(KIND=dp) :: energy_homo
87 : REAL(KIND=dp), DIMENSION(:) :: eigenvalues_ot
88 : REAL(KIND=dp) :: energy_gap
89 : INTEGER :: my_solver_type
90 :
91 : INTEGER :: precon_type
92 :
93 9466 : precon_type = preconditioner_env%in_use
94 38 : SELECT CASE (precon_type)
95 : CASE (ot_precond_full_single)
96 38 : IF (my_solver_type /= ot_precond_solver_default) THEN
97 0 : CPABORT("Only PRECOND_SOLVER DEFAULT for the moment")
98 : END IF
99 38 : IF (PRESENT(matrix_s)) THEN
100 : CALL make_full_single(preconditioner_env, preconditioner_env%fm, &
101 32 : matrix_h, matrix_s, energy_homo, energy_gap)
102 : ELSE
103 : CALL make_full_single_ortho(preconditioner_env, preconditioner_env%fm, &
104 6 : matrix_h, energy_homo, energy_gap)
105 : END IF
106 :
107 : CASE (ot_precond_s_inverse)
108 72 : IF (my_solver_type == ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
109 72 : IF (.NOT. PRESENT(matrix_s)) THEN
110 0 : CPABORT("Type for S=1 not implemented")
111 : END IF
112 72 : CALL make_full_s_inverse(preconditioner_env, matrix_s)
113 :
114 : CASE (ot_precond_full_kinetic)
115 1309 : IF (my_solver_type == ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
116 1309 : IF (.NOT. (PRESENT(matrix_s) .AND. PRESENT(matrix_t))) THEN
117 0 : CPABORT("Type for S=1 not implemented")
118 : END IF
119 1309 : CALL make_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
120 : CASE (ot_precond_full_single_inverse)
121 4399 : IF (my_solver_type == ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
122 : CALL make_full_single_inverse(preconditioner_env, mo_coeff, matrix_h, energy_gap, &
123 4399 : matrix_s=matrix_s)
124 : CASE (ot_precond_full_all)
125 3648 : IF (my_solver_type /= ot_precond_solver_default) THEN
126 0 : CPABORT("Only PRECOND_SOLVER DEFAULT for the moment")
127 : END IF
128 3648 : IF (PRESENT(matrix_s)) THEN
129 : CALL make_full_all(preconditioner_env, mo_coeff, matrix_h, matrix_s, &
130 3568 : eigenvalues_ot, energy_gap)
131 : ELSE
132 : CALL make_full_all_ortho(preconditioner_env, mo_coeff, matrix_h, &
133 80 : eigenvalues_ot, energy_gap)
134 : END IF
135 :
136 : CASE DEFAULT
137 9466 : CPABORT("Type not implemented")
138 : END SELECT
139 :
140 9466 : END SUBROUTINE make_preconditioner_matrix
141 :
142 : ! **************************************************************************************************
143 : !> \brief Simply takes the overlap matrix as preconditioner
144 : !> \param preconditioner_env ...
145 : !> \param matrix_s ...
146 : ! **************************************************************************************************
147 72 : SUBROUTINE make_full_s_inverse(preconditioner_env, matrix_s)
148 : TYPE(preconditioner_type) :: preconditioner_env
149 : TYPE(dbcsr_type), POINTER :: matrix_s
150 :
151 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_s_inverse'
152 :
153 : INTEGER :: handle
154 :
155 72 : CALL timeset(routineN, handle)
156 :
157 72 : CPASSERT(ASSOCIATED(matrix_s))
158 :
159 72 : IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
160 72 : ALLOCATE (preconditioner_env%sparse_matrix)
161 : END IF
162 72 : CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_s, name="full_kinetic")
163 :
164 72 : CALL timestop(handle)
165 :
166 72 : END SUBROUTINE make_full_s_inverse
167 :
168 : ! **************************************************************************************************
169 : !> \brief kinetic matrix+shift*overlap as preconditioner. Cheap but could
170 : !> be better
171 : !> \param preconditioner_env ...
172 : !> \param matrix_t ...
173 : !> \param matrix_s ...
174 : !> \param energy_gap ...
175 : ! **************************************************************************************************
176 1309 : SUBROUTINE make_full_kinetic(preconditioner_env, matrix_t, matrix_s, &
177 : energy_gap)
178 : TYPE(preconditioner_type) :: preconditioner_env
179 : TYPE(dbcsr_type), POINTER :: matrix_t, matrix_s
180 : REAL(KIND=dp) :: energy_gap
181 :
182 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_kinetic'
183 :
184 : INTEGER :: handle
185 : REAL(KIND=dp) :: shift
186 :
187 1309 : CALL timeset(routineN, handle)
188 :
189 1309 : CPASSERT(ASSOCIATED(matrix_t))
190 1309 : CPASSERT(ASSOCIATED(matrix_s))
191 :
192 1309 : IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
193 1307 : ALLOCATE (preconditioner_env%sparse_matrix)
194 : END IF
195 1309 : CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_t, name="full_kinetic")
196 :
197 1309 : shift = MAX(0.0_dp, energy_gap)
198 :
199 : CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, &
200 1309 : alpha_scalar=1.0_dp, beta_scalar=shift)
201 :
202 1309 : CALL timestop(handle)
203 :
204 1309 : END SUBROUTINE make_full_kinetic
205 :
206 : ! **************************************************************************************************
207 : !> \brief full_single_preconditioner
208 : !> \param preconditioner_env ...
209 : !> \param fm ...
210 : !> \param matrix_h ...
211 : !> \param matrix_s ...
212 : !> \param energy_homo ...
213 : !> \param energy_gap ...
214 : ! **************************************************************************************************
215 32 : SUBROUTINE make_full_single(preconditioner_env, fm, matrix_h, matrix_s, &
216 : energy_homo, energy_gap)
217 : TYPE(preconditioner_type) :: preconditioner_env
218 : TYPE(cp_fm_type), POINTER :: fm
219 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
220 : REAL(KIND=dp) :: energy_homo, energy_gap
221 :
222 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single'
223 :
224 : INTEGER :: handle, i, n
225 32 : REAL(KIND=dp), DIMENSION(:), POINTER :: evals
226 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
227 : TYPE(cp_fm_type) :: fm_h, fm_s
228 :
229 32 : CALL timeset(routineN, handle)
230 :
231 32 : NULLIFY (fm_struct_tmp, evals)
232 :
233 32 : IF (ASSOCIATED(fm)) THEN
234 0 : CALL cp_fm_release(fm)
235 0 : DEALLOCATE (fm)
236 : NULLIFY (fm)
237 : END IF
238 32 : CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
239 96 : ALLOCATE (evals(n))
240 :
241 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
242 : context=preconditioner_env%ctxt, &
243 32 : para_env=preconditioner_env%para_env)
244 32 : ALLOCATE (fm)
245 32 : CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
246 32 : CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
247 32 : CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
248 32 : CALL cp_fm_struct_release(fm_struct_tmp)
249 :
250 32 : CALL copy_dbcsr_to_fm(matrix_h, fm_h)
251 32 : CALL copy_dbcsr_to_fm(matrix_s, fm_s)
252 32 : CALL cp_fm_cholesky_decompose(fm_s)
253 :
254 32 : SELECT CASE (preconditioner_env%cholesky_use)
255 : CASE (cholesky_inverse)
256 : ! if cho inverse
257 0 : CALL cp_fm_triangular_invert(fm_s)
258 0 : CALL cp_fm_uplo_to_full(fm_h, fm)
259 :
260 : CALL cp_fm_triangular_multiply(fm_s, fm_h, side="R", transpose_tr=.FALSE., &
261 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
262 : CALL cp_fm_triangular_multiply(fm_s, fm_h, side="L", transpose_tr=.TRUE., &
263 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
264 : CASE (cholesky_reduce)
265 32 : CALL cp_fm_cholesky_reduce(fm_h, fm_s)
266 : CASE DEFAULT
267 32 : CPABORT("cholesky type not implemented")
268 : END SELECT
269 :
270 32 : CALL choose_eigv_solver(fm_h, fm, evals)
271 :
272 32 : SELECT CASE (preconditioner_env%cholesky_use)
273 : CASE (cholesky_inverse)
274 : CALL cp_fm_triangular_multiply(fm_s, fm, side="L", transpose_tr=.FALSE., &
275 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
276 0 : DO i = 1, n
277 0 : evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap)
278 : END DO
279 0 : CALL cp_fm_to_fm(fm, fm_h)
280 : CASE (cholesky_reduce)
281 32 : CALL cp_fm_cholesky_restore(fm, n, fm_s, fm_h, "SOLVE")
282 568 : DO i = 1, n
283 568 : evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap)
284 : END DO
285 64 : CALL cp_fm_to_fm(fm_h, fm)
286 : END SELECT
287 :
288 32 : CALL cp_fm_column_scale(fm, evals)
289 32 : CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
290 32 : CALL cp_fm_to_fm(fm_s, fm)
291 :
292 32 : DEALLOCATE (evals)
293 32 : CALL cp_fm_release(fm_h)
294 32 : CALL cp_fm_release(fm_s)
295 :
296 32 : CALL timestop(handle)
297 :
298 96 : END SUBROUTINE make_full_single
299 :
300 : ! **************************************************************************************************
301 : !> \brief full single in the orthonormal basis
302 : !> \param preconditioner_env ...
303 : !> \param fm ...
304 : !> \param matrix_h ...
305 : !> \param energy_homo ...
306 : !> \param energy_gap ...
307 : ! **************************************************************************************************
308 6 : SUBROUTINE make_full_single_ortho(preconditioner_env, fm, matrix_h, &
309 : energy_homo, energy_gap)
310 : TYPE(preconditioner_type) :: preconditioner_env
311 : TYPE(cp_fm_type), POINTER :: fm
312 : TYPE(dbcsr_type), POINTER :: matrix_h
313 : REAL(KIND=dp) :: energy_homo, energy_gap
314 :
315 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single_ortho'
316 :
317 : INTEGER :: handle, i, n
318 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: evals
319 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
320 : TYPE(cp_fm_type) :: fm_h, fm_s
321 :
322 6 : CALL timeset(routineN, handle)
323 6 : NULLIFY (fm_struct_tmp, evals)
324 :
325 6 : IF (ASSOCIATED(fm)) THEN
326 0 : CALL cp_fm_release(fm)
327 0 : DEALLOCATE (fm)
328 : NULLIFY (fm)
329 : END IF
330 6 : CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
331 18 : ALLOCATE (evals(n))
332 :
333 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
334 : context=preconditioner_env%ctxt, &
335 6 : para_env=preconditioner_env%para_env)
336 6 : ALLOCATE (fm)
337 6 : CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
338 6 : CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
339 6 : CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
340 6 : CALL cp_fm_struct_release(fm_struct_tmp)
341 :
342 6 : CALL copy_dbcsr_to_fm(matrix_h, fm_h)
343 :
344 6 : CALL choose_eigv_solver(fm_h, fm, evals)
345 282 : DO i = 1, n
346 282 : evals(i) = 1.0_dp/MAX(evals(i) - energy_homo, energy_gap)
347 : END DO
348 6 : CALL cp_fm_to_fm(fm, fm_h)
349 6 : CALL cp_fm_column_scale(fm, evals)
350 6 : CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
351 6 : CALL cp_fm_to_fm(fm_s, fm)
352 :
353 6 : DEALLOCATE (evals)
354 6 : CALL cp_fm_release(fm_h)
355 6 : CALL cp_fm_release(fm_s)
356 :
357 6 : CALL timestop(handle)
358 :
359 18 : END SUBROUTINE make_full_single_ortho
360 :
361 : ! **************************************************************************************************
362 : !> \brief generates a state by state preconditioner based on the full hamiltonian matrix
363 : !> \param preconditioner_env ...
364 : !> \param matrix_c0 ...
365 : !> \param matrix_h ...
366 : !> \param matrix_s ...
367 : !> \param c0_evals ...
368 : !> \param energy_gap should be a slight underestimate of the physical energy gap for almost all systems
369 : !> the c0 are already ritz states of (h,s)
370 : !> \par History
371 : !> 10.2006 made more stable [Joost VandeVondele]
372 : !> \note
373 : !> includes error estimate on the hamiltonian matrix to result in a stable preconditioner
374 : !> a preconditioner for each eigenstate i is generated by keeping the factorized form
375 : !> U diag( something i ) U^T. It is important to only precondition in the subspace orthogonal to c0.
376 : !> not only is it the only part that matters, it also simplifies the computation of
377 : !> the lagrangian multipliers in the OT minimization (i.e. if the c0 here is different
378 : !> from the c0 used in the OT setup, there will be a bug).
379 : ! **************************************************************************************************
380 3568 : SUBROUTINE make_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, c0_evals, energy_gap)
381 : TYPE(preconditioner_type) :: preconditioner_env
382 : TYPE(cp_fm_type), INTENT(IN) :: matrix_c0
383 : TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
384 : REAL(KIND=dp), DIMENSION(:) :: c0_evals
385 : REAL(KIND=dp) :: energy_gap
386 :
387 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all'
388 : REAL(KIND=dp), PARAMETER :: fudge_factor = 0.25_dp, &
389 : lambda_base = 10.0_dp
390 :
391 : INTEGER :: handle, k, n
392 : REAL(KIND=dp) :: error_estimate, lambda
393 3568 : REAL(KIND=dp), DIMENSION(:), POINTER :: diag, norms, shifted_evals
394 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
395 : TYPE(cp_fm_type) :: matrix_hc0, matrix_left, matrix_s1, &
396 : matrix_s2, matrix_sc0, matrix_shc0, &
397 : matrix_tmp, ortho
398 : TYPE(cp_fm_type), POINTER :: matrix_pre
399 :
400 3568 : CALL timeset(routineN, handle)
401 :
402 3568 : IF (ASSOCIATED(preconditioner_env%fm)) THEN
403 0 : CALL cp_fm_release(preconditioner_env%fm)
404 0 : DEALLOCATE (preconditioner_env%fm)
405 : NULLIFY (preconditioner_env%fm)
406 : END IF
407 3568 : CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
408 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
409 : context=preconditioner_env%ctxt, &
410 3568 : para_env=preconditioner_env%para_env)
411 3568 : ALLOCATE (preconditioner_env%fm)
412 3568 : CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
413 3568 : CALL cp_fm_create(ortho, fm_struct_tmp, name="ortho")
414 3568 : CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
415 3568 : CALL cp_fm_struct_release(fm_struct_tmp)
416 10704 : ALLOCATE (preconditioner_env%full_evals(n))
417 10598 : ALLOCATE (preconditioner_env%occ_evals(k))
418 :
419 : ! 0) cholesky decompose the overlap matrix, if this fails the basis is singular,
420 : ! more than EPS_DEFAULT
421 3568 : CALL copy_dbcsr_to_fm(matrix_s, ortho)
422 3568 : CALL cp_fm_cholesky_decompose(ortho)
423 : ! if cho inverse
424 3568 : IF (preconditioner_env%cholesky_use == cholesky_inverse) THEN
425 0 : CALL cp_fm_triangular_invert(ortho)
426 : END IF
427 : ! 1) Construct a new H matrix, which has the current C0 as eigenvectors,
428 : ! possibly shifted by an amount lambda,
429 : ! and the same spectrum as the original H matrix in the space orthogonal to the C0
430 : ! with P=C0 C0 ^ T
431 : ! (1 - PS)^T H (1-PS) + (PS)^T (H - lambda S ) (PS)
432 : ! we exploit that the C0 are already the ritz states of H
433 3568 : CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
434 3568 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, matrix_c0, matrix_sc0, k)
435 3568 : CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
436 3568 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, matrix_c0, matrix_hc0, k)
437 :
438 : ! An aside, try to estimate the error on the ritz values, we'll need it later on
439 3568 : CALL cp_fm_create(matrix_shc0, matrix_c0%matrix_struct, name="shc0")
440 :
441 3568 : SELECT CASE (preconditioner_env%cholesky_use)
442 : CASE (cholesky_inverse)
443 : ! if cho inverse
444 0 : CALL cp_fm_to_fm(matrix_hc0, matrix_shc0)
445 : CALL cp_fm_triangular_multiply(ortho, matrix_shc0, side="L", transpose_tr=.TRUE., &
446 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=k, alpha=1.0_dp)
447 : CASE (cholesky_reduce)
448 3568 : CALL cp_fm_cholesky_restore(matrix_hc0, k, ortho, matrix_shc0, "SOLVE", transa="T")
449 : CASE DEFAULT
450 3568 : CPABORT("cholesky type not implemented")
451 : END SELECT
452 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
453 : context=preconditioner_env%ctxt, &
454 3568 : para_env=preconditioner_env%para_env)
455 3568 : CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
456 3568 : CALL cp_fm_struct_release(fm_struct_tmp)
457 : ! since we only use diagonal elements this is a bit of a waste
458 3568 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_shc0, matrix_shc0, 0.0_dp, matrix_s1)
459 7030 : ALLOCATE (diag(k))
460 3568 : CALL cp_fm_get_diag(matrix_s1, diag)
461 20646 : error_estimate = MAXVAL(SQRT(ABS(diag - c0_evals**2)))
462 3568 : DEALLOCATE (diag)
463 3568 : CALL cp_fm_release(matrix_s1)
464 3568 : CALL cp_fm_release(matrix_shc0)
465 : ! we'll only use the energy gap, if our estimate of the error on the eigenvalues
466 : ! is small enough. A large error combined with a small energy gap would otherwise lead to
467 : ! an aggressive but bad preconditioner. Only when the error is small (MD), we can precondition
468 : ! aggressively
469 3568 : preconditioner_env%energy_gap = MAX(energy_gap, error_estimate*fudge_factor)
470 3568 : CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
471 3568 : matrix_pre => preconditioner_env%fm
472 3568 : CALL cp_fm_uplo_to_full(matrix_tmp, matrix_pre)
473 : ! tmp = H ( 1 - PS )
474 3568 : CALL parallel_gemm('N', 'T', n, n, k, -1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
475 :
476 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=n, &
477 : context=preconditioner_env%ctxt, &
478 3568 : para_env=preconditioner_env%para_env)
479 3568 : CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
480 3568 : CALL cp_fm_struct_release(fm_struct_tmp)
481 3568 : CALL parallel_gemm('T', 'N', k, n, n, 1.0_dp, matrix_c0, matrix_tmp, 0.0_dp, matrix_left)
482 : ! tmp = (1 - PS)^T H (1-PS)
483 3568 : CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
484 3568 : CALL cp_fm_release(matrix_left)
485 :
486 7030 : ALLOCATE (shifted_evals(k))
487 3568 : lambda = lambda_base + error_estimate
488 20540 : shifted_evals = c0_evals - lambda
489 3568 : CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
490 3568 : CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
491 3568 : CALL parallel_gemm('N', 'T', n, n, k, 1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
492 :
493 : ! 2) diagonalize this operator
494 3568 : SELECT CASE (preconditioner_env%cholesky_use)
495 : CASE (cholesky_inverse)
496 : CALL cp_fm_triangular_multiply(ortho, matrix_tmp, side="R", transpose_tr=.FALSE., &
497 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
498 : CALL cp_fm_triangular_multiply(ortho, matrix_tmp, side="L", transpose_tr=.TRUE., &
499 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
500 : CASE (cholesky_reduce)
501 3568 : CALL cp_fm_cholesky_reduce(matrix_tmp, ortho)
502 : END SELECT
503 3568 : CALL choose_eigv_solver(matrix_tmp, matrix_pre, preconditioner_env%full_evals)
504 3568 : SELECT CASE (preconditioner_env%cholesky_use)
505 : CASE (cholesky_inverse)
506 : CALL cp_fm_triangular_multiply(ortho, matrix_pre, side="L", transpose_tr=.FALSE., &
507 0 : invert_tr=.FALSE., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
508 0 : CALL cp_fm_to_fm(matrix_pre, matrix_tmp)
509 : CASE (cholesky_reduce)
510 3568 : CALL cp_fm_cholesky_restore(matrix_pre, n, ortho, matrix_tmp, "SOLVE")
511 7136 : CALL cp_fm_to_fm(matrix_tmp, matrix_pre)
512 : END SELECT
513 :
514 : ! test that the subspace remained conserved
515 : IF (.FALSE.) THEN
516 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
517 : context=preconditioner_env%ctxt, &
518 : para_env=preconditioner_env%para_env)
519 : CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
520 : CALL cp_fm_create(matrix_s2, fm_struct_tmp, name="matrix_s2")
521 : CALL cp_fm_struct_release(fm_struct_tmp)
522 : ALLOCATE (norms(k))
523 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_sc0, matrix_tmp, 0.0_dp, matrix_s1)
524 : CALL choose_eigv_solver(matrix_s1, matrix_s2, norms)
525 : WRITE (*, *) "matrix norm deviation (should be close to zero): ", MAXVAL(ABS(ABS(norms) - 1.0_dp))
526 : DEALLOCATE (norms)
527 : CALL cp_fm_release(matrix_s1)
528 : CALL cp_fm_release(matrix_s2)
529 : END IF
530 :
531 : ! 3) replace the lowest k evals and evecs with what they should be
532 20540 : preconditioner_env%occ_evals = c0_evals
533 : ! notice, this choice causes the preconditioner to be constant when applied to sc0 (see apply_full_all)
534 20540 : preconditioner_env%full_evals(1:k) = c0_evals
535 3568 : CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
536 :
537 3568 : CALL cp_fm_release(matrix_sc0)
538 3568 : CALL cp_fm_release(matrix_hc0)
539 3568 : CALL cp_fm_release(ortho)
540 3568 : CALL cp_fm_release(matrix_tmp)
541 3568 : DEALLOCATE (shifted_evals)
542 3568 : CALL timestop(handle)
543 :
544 28544 : END SUBROUTINE make_full_all
545 :
546 : ! **************************************************************************************************
547 : !> \brief full all in the orthonormal basis
548 : !> \param preconditioner_env ...
549 : !> \param matrix_c0 ...
550 : !> \param matrix_h ...
551 : !> \param c0_evals ...
552 : !> \param energy_gap ...
553 : ! **************************************************************************************************
554 80 : SUBROUTINE make_full_all_ortho(preconditioner_env, matrix_c0, matrix_h, c0_evals, energy_gap)
555 :
556 : TYPE(preconditioner_type) :: preconditioner_env
557 : TYPE(cp_fm_type), INTENT(IN) :: matrix_c0
558 : TYPE(dbcsr_type), POINTER :: matrix_h
559 : REAL(KIND=dp), DIMENSION(:) :: c0_evals
560 : REAL(KIND=dp) :: energy_gap
561 :
562 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all_ortho'
563 : REAL(KIND=dp), PARAMETER :: fudge_factor = 0.25_dp, &
564 : lambda_base = 10.0_dp
565 :
566 : INTEGER :: handle, k, n
567 : REAL(KIND=dp) :: error_estimate, lambda
568 80 : REAL(KIND=dp), DIMENSION(:), POINTER :: diag, norms, shifted_evals
569 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
570 : TYPE(cp_fm_type) :: matrix_hc0, matrix_left, matrix_s1, &
571 : matrix_s2, matrix_sc0, matrix_tmp
572 : TYPE(cp_fm_type), POINTER :: matrix_pre
573 :
574 80 : CALL timeset(routineN, handle)
575 :
576 80 : IF (ASSOCIATED(preconditioner_env%fm)) THEN
577 0 : CALL cp_fm_release(preconditioner_env%fm)
578 0 : DEALLOCATE (preconditioner_env%fm)
579 : NULLIFY (preconditioner_env%fm)
580 : END IF
581 80 : CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
582 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
583 : context=preconditioner_env%ctxt, &
584 80 : para_env=preconditioner_env%para_env)
585 80 : ALLOCATE (preconditioner_env%fm)
586 80 : CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
587 80 : CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
588 80 : CALL cp_fm_struct_release(fm_struct_tmp)
589 240 : ALLOCATE (preconditioner_env%full_evals(n))
590 240 : ALLOCATE (preconditioner_env%occ_evals(k))
591 :
592 : ! 1) Construct a new H matrix, which has the current C0 as eigenvectors,
593 : ! possibly shifted by an amount lambda,
594 : ! and the same spectrum as the original H matrix in the space orthogonal to the C0
595 : ! with P=C0 C0 ^ T
596 : ! (1 - PS)^T H (1-PS) + (PS)^T (H - lambda S ) (PS)
597 : ! we exploit that the C0 are already the ritz states of H
598 80 : CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
599 80 : CALL cp_fm_to_fm(matrix_c0, matrix_sc0)
600 80 : CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
601 80 : CALL cp_dbcsr_sm_fm_multiply(matrix_h, matrix_c0, matrix_hc0, k)
602 :
603 : ! An aside, try to estimate the error on the ritz values, we'll need it later on
604 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
605 : context=preconditioner_env%ctxt, &
606 80 : para_env=preconditioner_env%para_env)
607 80 : CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
608 80 : CALL cp_fm_struct_release(fm_struct_tmp)
609 : ! since we only use diagonal elements this is a bit of a waste
610 80 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_hc0, matrix_hc0, 0.0_dp, matrix_s1)
611 160 : ALLOCATE (diag(k))
612 80 : CALL cp_fm_get_diag(matrix_s1, diag)
613 826 : error_estimate = MAXVAL(SQRT(ABS(diag - c0_evals**2)))
614 80 : DEALLOCATE (diag)
615 80 : CALL cp_fm_release(matrix_s1)
616 : ! we'll only use the energy gap, if our estimate of the error on the eigenvalues
617 : ! is small enough. A large error combined with a small energy gap would otherwise lead to
618 : ! an aggressive but bad preconditioner. Only when the error is small (MD), we can precondition
619 : ! aggressively
620 80 : preconditioner_env%energy_gap = MAX(energy_gap, error_estimate*fudge_factor)
621 :
622 80 : matrix_pre => preconditioner_env%fm
623 80 : CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
624 80 : CALL cp_fm_uplo_to_full(matrix_tmp, matrix_pre)
625 : ! tmp = H ( 1 - PS )
626 80 : CALL parallel_gemm('N', 'T', n, n, k, -1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
627 :
628 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=n, &
629 : context=preconditioner_env%ctxt, &
630 80 : para_env=preconditioner_env%para_env)
631 80 : CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
632 80 : CALL cp_fm_struct_release(fm_struct_tmp)
633 80 : CALL parallel_gemm('T', 'N', k, n, n, 1.0_dp, matrix_c0, matrix_tmp, 0.0_dp, matrix_left)
634 : ! tmp = (1 - PS)^T H (1-PS)
635 80 : CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
636 80 : CALL cp_fm_release(matrix_left)
637 :
638 160 : ALLOCATE (shifted_evals(k))
639 80 : lambda = lambda_base + error_estimate
640 826 : shifted_evals = c0_evals - lambda
641 80 : CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
642 80 : CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
643 80 : CALL parallel_gemm('N', 'T', n, n, k, 1.0_dp, matrix_hc0, matrix_sc0, 1.0_dp, matrix_tmp)
644 :
645 : ! 2) diagonalize this operator
646 80 : CALL choose_eigv_solver(matrix_tmp, matrix_pre, preconditioner_env%full_evals)
647 :
648 : ! test that the subspace remained conserved
649 : IF (.FALSE.) THEN
650 : CALL cp_fm_to_fm(matrix_pre, matrix_tmp)
651 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=k, ncol_global=k, &
652 : context=preconditioner_env%ctxt, &
653 : para_env=preconditioner_env%para_env)
654 : CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
655 : CALL cp_fm_create(matrix_s2, fm_struct_tmp, name="matrix_s2")
656 : CALL cp_fm_struct_release(fm_struct_tmp)
657 : ALLOCATE (norms(k))
658 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_sc0, matrix_tmp, 0.0_dp, matrix_s1)
659 : CALL choose_eigv_solver(matrix_s1, matrix_s2, norms)
660 :
661 : WRITE (*, *) "matrix norm deviation (should be close to zero): ", MAXVAL(ABS(ABS(norms) - 1.0_dp))
662 : DEALLOCATE (norms)
663 : CALL cp_fm_release(matrix_s1)
664 : CALL cp_fm_release(matrix_s2)
665 : END IF
666 :
667 : ! 3) replace the lowest k evals and evecs with what they should be
668 826 : preconditioner_env%occ_evals = c0_evals
669 : ! notice, this choice causes the preconditioner to be constant when applied to sc0 (see apply_full_all)
670 826 : preconditioner_env%full_evals(1:k) = c0_evals
671 80 : CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
672 :
673 80 : CALL cp_fm_release(matrix_sc0)
674 80 : CALL cp_fm_release(matrix_hc0)
675 80 : CALL cp_fm_release(matrix_tmp)
676 80 : DEALLOCATE (shifted_evals)
677 :
678 80 : CALL timestop(handle)
679 :
680 560 : END SUBROUTINE make_full_all_ortho
681 :
682 : ! **************************************************************************************************
683 : !> \brief generates a preconditioner matrix H-lambda S+(SC)(2.0*CT*H*C+delta)(SC)^T
684 : !> for later inversion.
685 : !> H is the Kohn Sham matrix
686 : !> lambda*S shifts the spectrum of the generalized form up by lambda
687 : !> the last term only shifts the occupied space (reversing them in energy order)
688 : !> This form is implicitly multiplied from both sides by S^0.5
689 : !> This ensures we precondition the correct quantity
690 : !> Before this reads S^-0.5 H S^-0.5 + lambda + (S^0.5 C)shifts(S^0.5 C)T
691 : !> which might be a bit more obvious
692 : !> Replaced the old full_single_inverse at revision 14616
693 : !> \param preconditioner_env the preconditioner env
694 : !> \param matrix_c0 the MO coefficient matrix (fm)
695 : !> \param matrix_h Kohn-Sham matrix (dbcsr)
696 : !> \param energy_gap an additional shift in lambda=-E_homo+energy_gap
697 : !> \param matrix_s the overlap matrix if not orthonormal (dbcsr, optional)
698 : ! **************************************************************************************************
699 4399 : SUBROUTINE make_full_single_inverse(preconditioner_env, matrix_c0, matrix_h, energy_gap, matrix_s)
700 : TYPE(preconditioner_type) :: preconditioner_env
701 : TYPE(cp_fm_type), INTENT(IN) :: matrix_c0
702 : TYPE(dbcsr_type), POINTER :: matrix_h
703 : REAL(KIND=dp) :: energy_gap
704 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s
705 :
706 : CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single_inverse'
707 :
708 : INTEGER :: handle, k, n
709 : REAL(KIND=dp) :: max_ev, min_ev, pre_shift
710 : TYPE(arnoldi_env_type) :: arnoldi_env
711 4399 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices
712 : TYPE(dbcsr_type), TARGET :: dbcsr_cThc, dbcsr_hc, dbcsr_sc, mo_dbcsr
713 :
714 4399 : CALL timeset(routineN, handle)
715 :
716 : ! Allocate all working matrices needed
717 4399 : CALL cp_fm_get_info(matrix_c0, nrow_global=n, ncol_global=k)
718 : ! copy the fm MO's to a sparse matrix, can be solved better if the sparse version is already present
719 : ! but for the time beeing this will do
720 4399 : CALL cp_fm_to_dbcsr_row_template(mo_dbcsr, matrix_c0, matrix_h)
721 4399 : CALL dbcsr_create(dbcsr_sc, template=mo_dbcsr)
722 4399 : CALL dbcsr_create(dbcsr_hc, template=mo_dbcsr)
723 4399 : CALL cp_dbcsr_m_by_n_from_template(dbcsr_cThc, matrix_h, k, k, sym=dbcsr_type_symmetric)
724 :
725 : ! Check whether the output matrix was already created, if not do it now
726 4399 : IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
727 4399 : ALLOCATE (preconditioner_env%sparse_matrix)
728 : END IF
729 :
730 : ! Put the first term of the preconditioner (H) into the output matrix
731 4399 : CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_h)
732 :
733 : ! Precompute some matrices
734 : ! S*C, if orthonormal this will be simply C so a copy will do
735 4399 : IF (PRESENT(matrix_s)) THEN
736 4011 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, mo_dbcsr, 0.0_dp, dbcsr_sc)
737 : ELSE
738 388 : CALL dbcsr_copy(dbcsr_sc, mo_dbcsr)
739 : END IF
740 :
741 : !----------------------------compute the occupied subspace and shift it ------------------------------------
742 : ! cT*H*C which will be used to shift the occupied states to 0
743 4399 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_h, mo_dbcsr, 0.0_dp, dbcsr_hc)
744 4399 : CALL dbcsr_multiply("T", "N", 1.0_dp, mo_dbcsr, dbcsr_hc, 0.0_dp, dbcsr_cThc)
745 :
746 : ! Compute the Energy of the HOMO. We will use this as a reference energy
747 8798 : ALLOCATE (matrices(1))
748 4399 : matrices(1)%matrix => dbcsr_cThc
749 : CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=1.0E-3_dp, selection_crit=2, &
750 4399 : nval_request=1, nrestarts=8, generalized_ev=.FALSE., iram=.FALSE.)
751 4399 : IF (ASSOCIATED(preconditioner_env%max_ev_vector)) THEN
752 2416 : CALL set_arnoldi_initial_vector(arnoldi_env, preconditioner_env%max_ev_vector)
753 : END IF
754 4399 : CALL arnoldi_ev(matrices, arnoldi_env)
755 4399 : max_ev = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
756 :
757 : ! save the ev as guess for the next time
758 4399 : IF (.NOT. ASSOCIATED(preconditioner_env%max_ev_vector)) ALLOCATE (preconditioner_env%max_ev_vector)
759 4399 : CALL get_selected_ritz_vector(arnoldi_env, 1, matrices(1)%matrix, preconditioner_env%max_ev_vector)
760 4399 : CALL deallocate_arnoldi_env(arnoldi_env)
761 4399 : DEALLOCATE (matrices)
762 :
763 : ! Lets shift the occupied states a bit further up, -1.0 because we gonna subtract it from H
764 4399 : CALL dbcsr_add_on_diag(dbcsr_cThc, -0.5_dp)
765 : ! Get the AO representation of the shift (see above why S is needed), W-matrix like object
766 4399 : CALL dbcsr_multiply("N", "N", 2.0_dp, dbcsr_sc, dbcsr_cThc, 0.0_dp, dbcsr_hc)
767 4399 : CALL dbcsr_multiply("N", "T", -1.0_dp, dbcsr_hc, dbcsr_sc, 1.0_dp, preconditioner_env%sparse_matrix)
768 :
769 : !-------------------------------------compute eigenvalues of H ----------------------------------------------
770 : ! Setup the arnoldi procedure to compute the lowest ev. if S is present this has to be the generalized ev
771 4399 : IF (PRESENT(matrix_s)) THEN
772 12033 : ALLOCATE (matrices(2))
773 4011 : matrices(1)%matrix => preconditioner_env%sparse_matrix
774 4011 : matrices(2)%matrix => matrix_s
775 : CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=2.0E-2_dp, selection_crit=3, &
776 4011 : nval_request=1, nrestarts=21, generalized_ev=.TRUE., iram=.FALSE.)
777 : ELSE
778 776 : ALLOCATE (matrices(1))
779 388 : matrices(1)%matrix => preconditioner_env%sparse_matrix
780 : CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=2.0E-2_dp, selection_crit=3, &
781 388 : nval_request=1, nrestarts=8, generalized_ev=.FALSE., iram=.FALSE.)
782 : END IF
783 4399 : IF (ASSOCIATED(preconditioner_env%min_ev_vector)) THEN
784 2416 : CALL set_arnoldi_initial_vector(arnoldi_env, preconditioner_env%min_ev_vector)
785 : END IF
786 :
787 : ! compute the LUMO energy
788 4399 : CALL arnoldi_ev(matrices, arnoldi_env)
789 4399 : min_eV = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
790 :
791 : ! save the lumo vector for restarting in the next step
792 4399 : IF (.NOT. ASSOCIATED(preconditioner_env%min_ev_vector)) ALLOCATE (preconditioner_env%min_ev_vector)
793 4399 : CALL get_selected_ritz_vector(arnoldi_env, 1, matrices(1)%matrix, preconditioner_env%min_ev_vector)
794 4399 : CALL deallocate_arnoldi_env(arnoldi_env)
795 4399 : DEALLOCATE (matrices)
796 :
797 : !-------------------------------------compute eigenvalues of H ----------------------------------------------
798 : ! Shift the Lumo to the 1.5*the computed energy_gap or the external energy gap value
799 : ! The factor 1.5 is determined by trying. If the LUMO is positive, enough, just leave it alone
800 4399 : pre_shift = MAX(1.5_dp*(min_ev - max_ev), energy_gap)
801 4399 : IF (min_ev < pre_shift) THEN
802 4381 : pre_shift = pre_shift - min_ev
803 : ELSE
804 18 : pre_shift = 0.0_dp
805 : END IF
806 4399 : IF (PRESENT(matrix_s)) THEN
807 4011 : CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, 1.0_dp, pre_shift)
808 : ELSE
809 388 : CALL dbcsr_add_on_diag(preconditioner_env%sparse_matrix, pre_shift)
810 : END IF
811 :
812 4399 : CALL dbcsr_release(mo_dbcsr)
813 4399 : CALL dbcsr_release(dbcsr_hc)
814 4399 : CALL dbcsr_release(dbcsr_sc)
815 4399 : CALL dbcsr_release(dbcsr_cThc)
816 :
817 4399 : CALL timestop(handle)
818 :
819 4399 : END SUBROUTINE make_full_single_inverse
820 :
821 : END MODULE preconditioner_makes
822 :
|
