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