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