(git:98357aa)
Loading...
Searching...
No Matches
preconditioner_makes.F
Go to the documentation of this file.
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! **************************************************************************************************
15MODULE 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
65CONTAINS
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 SUBROUTINE make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
80 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 precon_type = preconditioner_env%in_use
94 SELECT CASE (precon_type)
95 CASE (ot_precond_full_single)
96 IF (my_solver_type /= ot_precond_solver_default) THEN
97 cpabort("Only PRECOND_SOLVER DEFAULT for the moment")
98 END IF
99 IF (PRESENT(matrix_s)) THEN
100 CALL make_full_single(preconditioner_env, preconditioner_env%fm, &
101 matrix_h, matrix_s, energy_homo, energy_gap)
102 ELSE
103 CALL make_full_single_ortho(preconditioner_env, preconditioner_env%fm, &
104 matrix_h, energy_homo, energy_gap)
105 END IF
106
107 CASE (ot_precond_s_inverse)
108 IF (my_solver_type == ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
109 IF (.NOT. PRESENT(matrix_s)) THEN
110 cpabort("Type for S=1 not implemented")
111 END IF
112 CALL make_full_s_inverse(preconditioner_env, matrix_s)
113
114 CASE (ot_precond_full_kinetic)
115 IF (my_solver_type == ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
116 IF (.NOT. (PRESENT(matrix_s) .AND. PRESENT(matrix_t))) THEN
117 cpabort("Type for S=1 not implemented")
118 END IF
119 CALL make_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
120 CASE (ot_precond_full_single_inverse)
121 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 matrix_s=matrix_s)
124 CASE (ot_precond_full_all)
125 IF (my_solver_type /= ot_precond_solver_default) THEN
126 cpabort("Only PRECOND_SOLVER DEFAULT for the moment")
127 END IF
128 IF (PRESENT(matrix_s)) THEN
129 CALL make_full_all(preconditioner_env, mo_coeff, matrix_h, matrix_s, &
130 eigenvalues_ot, energy_gap)
131 ELSE
132 CALL make_full_all_ortho(preconditioner_env, mo_coeff, matrix_h, &
133 eigenvalues_ot, energy_gap)
134 END IF
135
136 CASE DEFAULT
137 cpabort("Type not implemented")
138 END SELECT
139
140 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 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 CALL timeset(routinen, handle)
156
157 cpassert(ASSOCIATED(matrix_s))
158
159 IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
160 ALLOCATE (preconditioner_env%sparse_matrix)
161 END IF
162 CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_s, name="full_kinetic")
163
164 CALL timestop(handle)
165
166 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 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 CALL timeset(routinen, handle)
188
189 cpassert(ASSOCIATED(matrix_t))
190 cpassert(ASSOCIATED(matrix_s))
191
192 IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
193 ALLOCATE (preconditioner_env%sparse_matrix)
194 END IF
195 CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_t, name="full_kinetic")
196
197 shift = max(0.0_dp, energy_gap)
198
199 CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, &
200 alpha_scalar=1.0_dp, beta_scalar=shift)
201
202 CALL timestop(handle)
203
204 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 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 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 CALL timeset(routinen, handle)
230
231 NULLIFY (fm_struct_tmp, evals)
232
233 IF (ASSOCIATED(fm)) THEN
234 CALL cp_fm_release(fm)
235 DEALLOCATE (fm)
236 NULLIFY (fm)
237 END IF
238 CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
239 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 para_env=preconditioner_env%para_env)
244 ALLOCATE (fm)
245 CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
246 CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
247 CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
248 CALL cp_fm_struct_release(fm_struct_tmp)
249
250 CALL copy_dbcsr_to_fm(matrix_h, fm_h)
251 CALL copy_dbcsr_to_fm(matrix_s, fm_s)
252 CALL cp_fm_cholesky_decompose(fm_s)
253
254 SELECT CASE (preconditioner_env%cholesky_use)
255 CASE (cholesky_inverse)
256! if cho inverse
257 CALL cp_fm_triangular_invert(fm_s)
258 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 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 invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
264 CASE (cholesky_reduce)
265 CALL cp_fm_cholesky_reduce(fm_h, fm_s)
266 CASE DEFAULT
267 cpabort("cholesky type not implemented")
268 END SELECT
269
270 CALL choose_eigv_solver(fm_h, fm, evals)
271
272 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 invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
276 DO i = 1, n
277 evals(i) = 1.0_dp/max(evals(i) - energy_homo, energy_gap)
278 END DO
279 CALL cp_fm_to_fm(fm, fm_h)
280 CASE (cholesky_reduce)
281 CALL cp_fm_cholesky_restore(fm, n, fm_s, fm_h, "SOLVE")
282 DO i = 1, n
283 evals(i) = 1.0_dp/max(evals(i) - energy_homo, energy_gap)
284 END DO
285 CALL cp_fm_to_fm(fm_h, fm)
286 END SELECT
287
288 CALL cp_fm_column_scale(fm, evals)
289 CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
290 CALL cp_fm_to_fm(fm_s, fm)
291
292 DEALLOCATE (evals)
293 CALL cp_fm_release(fm_h)
294 CALL cp_fm_release(fm_s)
295
296 CALL timestop(handle)
297
298 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 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 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 CALL timeset(routinen, handle)
323 NULLIFY (fm_struct_tmp, evals)
324
325 IF (ASSOCIATED(fm)) THEN
326 CALL cp_fm_release(fm)
327 DEALLOCATE (fm)
328 NULLIFY (fm)
329 END IF
330 CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
331 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 para_env=preconditioner_env%para_env)
336 ALLOCATE (fm)
337 CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
338 CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
339 CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
340 CALL cp_fm_struct_release(fm_struct_tmp)
341
342 CALL copy_dbcsr_to_fm(matrix_h, fm_h)
343
344 CALL choose_eigv_solver(fm_h, fm, evals)
345 DO i = 1, n
346 evals(i) = 1.0_dp/max(evals(i) - energy_homo, energy_gap)
347 END DO
348 CALL cp_fm_to_fm(fm, fm_h)
349 CALL cp_fm_column_scale(fm, evals)
350 CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
351 CALL cp_fm_to_fm(fm_s, fm)
352
353 DEALLOCATE (evals)
354 CALL cp_fm_release(fm_h)
355 CALL cp_fm_release(fm_s)
356
357 CALL timestop(handle)
358
359 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 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 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 CALL timeset(routinen, handle)
401
402 IF (ASSOCIATED(preconditioner_env%fm)) THEN
403 CALL cp_fm_release(preconditioner_env%fm)
404 DEALLOCATE (preconditioner_env%fm)
405 NULLIFY (preconditioner_env%fm)
406 END IF
407 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 para_env=preconditioner_env%para_env)
411 ALLOCATE (preconditioner_env%fm)
412 CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
413 CALL cp_fm_create(ortho, fm_struct_tmp, name="ortho")
414 CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
415 CALL cp_fm_struct_release(fm_struct_tmp)
416 ALLOCATE (preconditioner_env%full_evals(n))
417 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 CALL copy_dbcsr_to_fm(matrix_s, ortho)
422 CALL cp_fm_cholesky_decompose(ortho)
423! if cho inverse
424 IF (preconditioner_env%cholesky_use == cholesky_inverse) THEN
425 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 CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
434 CALL cp_dbcsr_sm_fm_multiply(matrix_s, matrix_c0, matrix_sc0, k)
435 CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
436 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 CALL cp_fm_create(matrix_shc0, matrix_c0%matrix_struct, name="shc0")
440
441 SELECT CASE (preconditioner_env%cholesky_use)
442 CASE (cholesky_inverse)
443! if cho inverse
444 CALL cp_fm_to_fm(matrix_hc0, matrix_shc0)
445 CALL cp_fm_triangular_multiply(ortho, matrix_shc0, side="L", transpose_tr=.true., &
446 invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=k, alpha=1.0_dp)
447 CASE (cholesky_reduce)
448 CALL cp_fm_cholesky_restore(matrix_hc0, k, ortho, matrix_shc0, "SOLVE", transa="T")
449 CASE DEFAULT
450 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 para_env=preconditioner_env%para_env)
455 CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
456 CALL cp_fm_struct_release(fm_struct_tmp)
457 ! since we only use diagonal elements this is a bit of a waste
458 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_shc0, matrix_shc0, 0.0_dp, matrix_s1)
459 ALLOCATE (diag(k))
460 CALL cp_fm_get_diag(matrix_s1, diag)
461 error_estimate = maxval(sqrt(abs(diag - c0_evals**2)))
462 DEALLOCATE (diag)
463 CALL cp_fm_release(matrix_s1)
464 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 preconditioner_env%energy_gap = max(energy_gap, error_estimate*fudge_factor)
470 CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
471 matrix_pre => preconditioner_env%fm
472 CALL cp_fm_uplo_to_full(matrix_tmp, matrix_pre)
473 ! tmp = H ( 1 - PS )
474 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 para_env=preconditioner_env%para_env)
479 CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
480 CALL cp_fm_struct_release(fm_struct_tmp)
481 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 CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
484 CALL cp_fm_release(matrix_left)
485
486 ALLOCATE (shifted_evals(k))
487 lambda = lambda_base + error_estimate
488 shifted_evals = c0_evals - lambda
489 CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
490 CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
491 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 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 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 invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
500 CASE (cholesky_reduce)
501 CALL cp_fm_cholesky_reduce(matrix_tmp, ortho)
502 END SELECT
503 CALL choose_eigv_solver(matrix_tmp, matrix_pre, preconditioner_env%full_evals)
504 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 invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
508 CALL cp_fm_to_fm(matrix_pre, matrix_tmp)
509 CASE (cholesky_reduce)
510 CALL cp_fm_cholesky_restore(matrix_pre, n, ortho, matrix_tmp, "SOLVE")
511 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 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 preconditioner_env%full_evals(1:k) = c0_evals
535 CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
536
537 CALL cp_fm_release(matrix_sc0)
538 CALL cp_fm_release(matrix_hc0)
539 CALL cp_fm_release(ortho)
540 CALL cp_fm_release(matrix_tmp)
541 DEALLOCATE (shifted_evals)
542 CALL timestop(handle)
543
544 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 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 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 CALL timeset(routinen, handle)
575
576 IF (ASSOCIATED(preconditioner_env%fm)) THEN
577 CALL cp_fm_release(preconditioner_env%fm)
578 DEALLOCATE (preconditioner_env%fm)
579 NULLIFY (preconditioner_env%fm)
580 END IF
581 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 para_env=preconditioner_env%para_env)
585 ALLOCATE (preconditioner_env%fm)
586 CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
587 CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
588 CALL cp_fm_struct_release(fm_struct_tmp)
589 ALLOCATE (preconditioner_env%full_evals(n))
590 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 CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
599 CALL cp_fm_to_fm(matrix_c0, matrix_sc0)
600 CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
601 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 para_env=preconditioner_env%para_env)
607 CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
608 CALL cp_fm_struct_release(fm_struct_tmp)
609 ! since we only use diagonal elements this is a bit of a waste
610 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_hc0, matrix_hc0, 0.0_dp, matrix_s1)
611 ALLOCATE (diag(k))
612 CALL cp_fm_get_diag(matrix_s1, diag)
613 error_estimate = maxval(sqrt(abs(diag - c0_evals**2)))
614 DEALLOCATE (diag)
615 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 preconditioner_env%energy_gap = max(energy_gap, error_estimate*fudge_factor)
621
622 matrix_pre => preconditioner_env%fm
623 CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
624 CALL cp_fm_uplo_to_full(matrix_tmp, matrix_pre)
625 ! tmp = H ( 1 - PS )
626 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 para_env=preconditioner_env%para_env)
631 CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
632 CALL cp_fm_struct_release(fm_struct_tmp)
633 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 CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
636 CALL cp_fm_release(matrix_left)
637
638 ALLOCATE (shifted_evals(k))
639 lambda = lambda_base + error_estimate
640 shifted_evals = c0_evals - lambda
641 CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
642 CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
643 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 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 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 preconditioner_env%full_evals(1:k) = c0_evals
671 CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
672
673 CALL cp_fm_release(matrix_sc0)
674 CALL cp_fm_release(matrix_hc0)
675 CALL cp_fm_release(matrix_tmp)
676 DEALLOCATE (shifted_evals)
677
678 CALL timestop(handle)
679
680 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 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices
712 TYPE(dbcsr_type), TARGET :: dbcsr_cthc, dbcsr_hc, dbcsr_sc, mo_dbcsr
713
714 CALL timeset(routinen, handle)
715
716 ! Allocate all working matrices needed
717 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 CALL cp_fm_to_dbcsr_row_template(mo_dbcsr, matrix_c0, matrix_h)
721 CALL dbcsr_create(dbcsr_sc, template=mo_dbcsr)
722 CALL dbcsr_create(dbcsr_hc, template=mo_dbcsr)
723 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 IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
727 ALLOCATE (preconditioner_env%sparse_matrix)
728 END IF
729
730 ! Put the first term of the preconditioner (H) into the output matrix
731 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 IF (PRESENT(matrix_s)) THEN
736 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, mo_dbcsr, 0.0_dp, dbcsr_sc)
737 ELSE
738 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 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_h, mo_dbcsr, 0.0_dp, dbcsr_hc)
744 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 ALLOCATE (matrices(1))
748 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 nval_request=1, nrestarts=8, generalized_ev=.false., iram=.false.)
751 IF (ASSOCIATED(preconditioner_env%max_ev_vector)) THEN
752 CALL set_arnoldi_initial_vector(arnoldi_env, preconditioner_env%max_ev_vector)
753 END IF
754 CALL arnoldi_ev(matrices, arnoldi_env)
755 max_ev = real(get_selected_ritz_val(arnoldi_env, 1), dp)
756
757 ! save the ev as guess for the next time
758 IF (.NOT. ASSOCIATED(preconditioner_env%max_ev_vector)) ALLOCATE (preconditioner_env%max_ev_vector)
759 CALL get_selected_ritz_vector(arnoldi_env, 1, matrices(1)%matrix, preconditioner_env%max_ev_vector)
760 CALL deallocate_arnoldi_env(arnoldi_env)
761 DEALLOCATE (matrices)
762
763 ! Lets shift the occupied states a bit further up, -1.0 because we gonna subtract it from H
764 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 CALL dbcsr_multiply("N", "N", 2.0_dp, dbcsr_sc, dbcsr_cthc, 0.0_dp, dbcsr_hc)
767 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 IF (PRESENT(matrix_s)) THEN
772 ALLOCATE (matrices(2))
773 matrices(1)%matrix => preconditioner_env%sparse_matrix
774 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 nval_request=1, nrestarts=21, generalized_ev=.true., iram=.false.)
777 ELSE
778 ALLOCATE (matrices(1))
779 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 nval_request=1, nrestarts=8, generalized_ev=.false., iram=.false.)
782 END IF
783 IF (ASSOCIATED(preconditioner_env%min_ev_vector)) THEN
784 CALL set_arnoldi_initial_vector(arnoldi_env, preconditioner_env%min_ev_vector)
785 END IF
786
787 ! compute the LUMO energy
788 CALL arnoldi_ev(matrices, arnoldi_env)
789 min_ev = real(get_selected_ritz_val(arnoldi_env, 1), dp)
790
791 ! save the lumo vector for restarting in the next step
792 IF (.NOT. ASSOCIATED(preconditioner_env%min_ev_vector)) ALLOCATE (preconditioner_env%min_ev_vector)
793 CALL get_selected_ritz_vector(arnoldi_env, 1, matrices(1)%matrix, preconditioner_env%min_ev_vector)
794 CALL deallocate_arnoldi_env(arnoldi_env)
795 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 pre_shift = max(1.5_dp*(min_ev - max_ev), energy_gap)
801 IF (min_ev < pre_shift) THEN
802 pre_shift = pre_shift - min_ev
803 ELSE
804 pre_shift = 0.0_dp
805 END IF
806 IF (PRESENT(matrix_s)) THEN
807 CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, 1.0_dp, pre_shift)
808 ELSE
809 CALL dbcsr_add_on_diag(preconditioner_env%sparse_matrix, pre_shift)
810 END IF
811
812 CALL dbcsr_release(mo_dbcsr)
813 CALL dbcsr_release(dbcsr_hc)
814 CALL dbcsr_release(dbcsr_sc)
815 CALL dbcsr_release(dbcsr_cthc)
816
817 CALL timestop(handle)
818
819 END SUBROUTINE make_full_single_inverse
820
821END MODULE preconditioner_makes
822
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:16
arnoldi_api::arnoldi_ev
subroutine, public arnoldi_ev(matrix, arnoldi_env)
Driver routine for different arnoldi eigenvalue methods the selection which one is to be taken is mad...
Definition arnoldi_api.F:69
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_add_on_diag
subroutine, public dbcsr_add_on_diag(matrix, alpha)
Adds the given scalar to the diagonal of the matrix. Reserves any missing diagonal blocks.
Definition cp_dbcsr_contrib.F:331
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:548
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:212
cp_dbcsr_operations::cp_fm_to_dbcsr_row_template
subroutine, public cp_fm_to_dbcsr_row_template(matrix, fm_in, template)
Utility function to copy a specially shaped fm to dbcsr_matrix The result matrix will be the matrix i...
Definition cp_dbcsr_operations.F:858
cp_dbcsr_operations::cp_dbcsr_m_by_n_from_template
subroutine, public cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym)
Utility function to create an arbitrary shaped dbcsr matrix with the same processor grid as the templ...
Definition cp_dbcsr_operations.F:892
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1690
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1584
cp_fm_basic_linalg::cp_fm_triangular_invert
subroutine, public cp_fm_triangular_invert(matrix_a, uplo_tr)
inverts a triangular matrix
Definition cp_fm_basic_linalg.F:1968
cp_fm_basic_linalg::cp_fm_triangular_multiply
subroutine, public cp_fm_triangular_multiply(triangular_matrix, matrix_b, side, transpose_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
Definition cp_fm_basic_linalg.F:1413
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa)
apply Cholesky decomposition op can be "SOLVE" (out = U^-1 * in) or "MULTIPLY" (out = U * in) pos can...
Definition cp_fm_cholesky.F:313
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_fm_cholesky.F:61
cp_fm_cholesky::cp_fm_cholesky_reduce
subroutine, public cp_fm_cholesky_reduce(matrix, matrixb, itype)
reduce a matrix pencil A,B to normal form B has to be cholesky decomposed with cp_fm_cholesky_decompo...
Definition cp_fm_cholesky.F:255
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:245
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition cp_fm_types.F:564
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::ot_precond_full_kinetic
integer, parameter, public ot_precond_full_kinetic
Definition input_constants.F:531
input_constants::cholesky_reduce
integer, parameter, public cholesky_reduce
Definition input_constants.F:206
input_constants::cholesky_inverse
integer, parameter, public cholesky_inverse
Definition input_constants.F:206
input_constants::ot_precond_solver_default
integer, parameter, public ot_precond_solver_default
Definition input_constants.F:538
input_constants::ot_precond_full_single
integer, parameter, public ot_precond_full_single
Definition input_constants.F:531
input_constants::ot_precond_solver_inv_chol
integer, parameter, public ot_precond_solver_inv_chol
Definition input_constants.F:538
input_constants::ot_precond_full_single_inverse
integer, parameter, public ot_precond_full_single_inverse
Definition input_constants.F:531
input_constants::ot_precond_s_inverse
integer, parameter, public ot_precond_s_inverse
Definition input_constants.F:531
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition input_constants.F:531
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
preconditioner_makes
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition preconditioner_makes.F:15
preconditioner_makes::make_preconditioner_matrix
subroutine, public make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, energy_homo, eigenvalues_ot, energy_gap, my_solver_type)
...
Definition preconditioner_makes.F:82
preconditioner_types
types of preconditioners
Definition preconditioner_types.F:14
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
preconditioner_types::preconditioner_type
Definition preconditioner_types.F:42