(git:e7e05ae)
preconditioner_makes.F
Go to the documentation of this file.
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_operations, ONLY: copy_dbcsr_to_fm,&
24  cp_dbcsr_m_by_n_from_template,&
25  cp_dbcsr_sm_fm_multiply,&
26  cp_fm_to_dbcsr_row_template
27  USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
28  cp_fm_triangular_invert,&
29  cp_fm_triangular_multiply,&
30  cp_fm_upper_to_full
31  USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
32  cp_fm_cholesky_reduce,&
33  cp_fm_cholesky_restore
34  USE cp_fm_diag, ONLY: choose_eigv_solver
35  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
36  cp_fm_struct_release,&
37  cp_fm_struct_type
38  USE cp_fm_types, ONLY: cp_fm_create,&
39  cp_fm_get_diag,&
40  cp_fm_get_info,&
41  cp_fm_release,&
42  cp_fm_to_fm,&
43  cp_fm_type
44  USE dbcsr_api, ONLY: &
45  dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_multiply, &
46  dbcsr_p_type, dbcsr_release, dbcsr_type, dbcsr_type_symmetric
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  SUBROUTINE make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
79  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  precon_type = preconditioner_env%in_use
93  SELECT CASE (precon_type)
94  CASE (ot_precond_full_single)
95  IF (my_solver_type .NE. ot_precond_solver_default) &
96  cpabort("Only PRECOND_SOLVER DEFAULT for the moment")
97  IF (PRESENT(matrix_s)) THEN
98  CALL make_full_single(preconditioner_env, preconditioner_env%fm, &
99  matrix_h, matrix_s, energy_homo, energy_gap)
100  ELSE
101  CALL make_full_single_ortho(preconditioner_env, preconditioner_env%fm, &
102  matrix_h, energy_homo, energy_gap)
103  END IF
104 
105  CASE (ot_precond_s_inverse)
106  IF (my_solver_type .EQ. ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
107  IF (.NOT. PRESENT(matrix_s)) &
108  cpabort("Type for S=1 not implemented")
109  CALL make_full_s_inverse(preconditioner_env, matrix_s)
110 
111  CASE (ot_precond_full_kinetic)
112  IF (my_solver_type .EQ. ot_precond_solver_default) my_solver_type = ot_precond_solver_inv_chol
113  IF (.NOT. (PRESENT(matrix_s) .AND. PRESENT(matrix_t))) &
114  cpabort("Type for S=1 not implemented")
115  CALL make_full_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap)
116  CASE (ot_precond_full_single_inverse)
117  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  matrix_s=matrix_s)
120  CASE (ot_precond_full_all)
121  IF (my_solver_type .NE. ot_precond_solver_default) THEN
122  cpabort("Only PRECOND_SOLVER DEFAULT for the moment")
123  END IF
124  IF (PRESENT(matrix_s)) THEN
125  CALL make_full_all(preconditioner_env, mo_coeff, matrix_h, matrix_s, &
126  eigenvalues_ot, energy_gap)
127  ELSE
128  CALL make_full_all_ortho(preconditioner_env, mo_coeff, matrix_h, &
129  eigenvalues_ot, energy_gap)
130  END IF
131 
132  CASE DEFAULT
133  cpabort("Type not implemented")
134  END SELECT
135 
136  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  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  CALL timeset(routinen, handle)
152 
153  cpassert(ASSOCIATED(matrix_s))
154 
155  IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
156  ALLOCATE (preconditioner_env%sparse_matrix)
157  END IF
158  CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_s, name="full_kinetic")
159 
160  CALL timestop(handle)
161 
162  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  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  CALL timeset(routinen, handle)
184 
185  cpassert(ASSOCIATED(matrix_t))
186  cpassert(ASSOCIATED(matrix_s))
187 
188  IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
189  ALLOCATE (preconditioner_env%sparse_matrix)
190  END IF
191  CALL dbcsr_copy(preconditioner_env%sparse_matrix, matrix_t, name="full_kinetic")
192 
193  shift = max(0.0_dp, energy_gap)
194 
195  CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, &
196  alpha_scalar=1.0_dp, beta_scalar=shift)
197 
198  CALL timestop(handle)
199 
200  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  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  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  CALL timeset(routinen, handle)
226 
227  NULLIFY (fm_struct_tmp, evals)
228 
229  IF (ASSOCIATED(fm)) THEN
230  CALL cp_fm_release(fm)
231  DEALLOCATE (fm)
232  NULLIFY (fm)
233  END IF
234  CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
235  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  para_env=preconditioner_env%para_env)
240  ALLOCATE (fm)
241  CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
242  CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
243  CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
244  CALL cp_fm_struct_release(fm_struct_tmp)
245 
246  CALL copy_dbcsr_to_fm(matrix_h, fm_h)
247  CALL copy_dbcsr_to_fm(matrix_s, fm_s)
248  CALL cp_fm_cholesky_decompose(fm_s)
249 
250  SELECT CASE (preconditioner_env%cholesky_use)
251  CASE (cholesky_inverse)
252 ! if cho inverse
253  CALL cp_fm_triangular_invert(fm_s)
254  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  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  invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
260  CASE (cholesky_reduce)
261  CALL cp_fm_cholesky_reduce(fm_h, fm_s)
262  CASE DEFAULT
263  cpabort("cholesky type not implemented")
264  END SELECT
265 
266  CALL choose_eigv_solver(fm_h, fm, evals)
267 
268  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  invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
272  DO i = 1, n
273  evals(i) = 1.0_dp/max(evals(i) - energy_homo, energy_gap)
274  END DO
275  CALL cp_fm_to_fm(fm, fm_h)
276  CASE (cholesky_reduce)
277  CALL cp_fm_cholesky_restore(fm, n, fm_s, fm_h, "SOLVE")
278  DO i = 1, n
279  evals(i) = 1.0_dp/max(evals(i) - energy_homo, energy_gap)
280  END DO
281  CALL cp_fm_to_fm(fm_h, fm)
282  END SELECT
283 
284  CALL cp_fm_column_scale(fm, evals)
285  CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
286  CALL cp_fm_to_fm(fm_s, fm)
287 
288  DEALLOCATE (evals)
289  CALL cp_fm_release(fm_h)
290  CALL cp_fm_release(fm_s)
291 
292  CALL timestop(handle)
293 
294  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  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  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  CALL timeset(routinen, handle)
319  NULLIFY (fm_struct_tmp, evals)
320 
321  IF (ASSOCIATED(fm)) THEN
322  CALL cp_fm_release(fm)
323  DEALLOCATE (fm)
324  NULLIFY (fm)
325  END IF
326  CALL dbcsr_get_info(matrix_h, nfullrows_total=n)
327  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  para_env=preconditioner_env%para_env)
332  ALLOCATE (fm)
333  CALL cp_fm_create(fm, fm_struct_tmp, name="preconditioner")
334  CALL cp_fm_create(fm_h, fm_struct_tmp, name="fm_h")
335  CALL cp_fm_create(fm_s, fm_struct_tmp, name="fm_s")
336  CALL cp_fm_struct_release(fm_struct_tmp)
337 
338  CALL copy_dbcsr_to_fm(matrix_h, fm_h)
339 
340  CALL choose_eigv_solver(fm_h, fm, evals)
341  DO i = 1, n
342  evals(i) = 1.0_dp/max(evals(i) - energy_homo, energy_gap)
343  END DO
344  CALL cp_fm_to_fm(fm, fm_h)
345  CALL cp_fm_column_scale(fm, evals)
346  CALL parallel_gemm('N', 'T', n, n, n, 1.0_dp, fm, fm_h, 0.0_dp, fm_s)
347  CALL cp_fm_to_fm(fm_s, fm)
348 
349  DEALLOCATE (evals)
350  CALL cp_fm_release(fm_h)
351  CALL cp_fm_release(fm_s)
352 
353  CALL timestop(handle)
354 
355  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  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  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  CALL timeset(routinen, handle)
397 
398  IF (ASSOCIATED(preconditioner_env%fm)) THEN
399  CALL cp_fm_release(preconditioner_env%fm)
400  DEALLOCATE (preconditioner_env%fm)
401  NULLIFY (preconditioner_env%fm)
402  END IF
403  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  para_env=preconditioner_env%para_env)
407  ALLOCATE (preconditioner_env%fm)
408  CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
409  CALL cp_fm_create(ortho, fm_struct_tmp, name="ortho")
410  CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
411  CALL cp_fm_struct_release(fm_struct_tmp)
412  ALLOCATE (preconditioner_env%full_evals(n))
413  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  CALL copy_dbcsr_to_fm(matrix_s, ortho)
418  CALL cp_fm_cholesky_decompose(ortho)
419 ! if cho inverse
420  IF (preconditioner_env%cholesky_use == cholesky_inverse) THEN
421  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  CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
430  CALL cp_dbcsr_sm_fm_multiply(matrix_s, matrix_c0, matrix_sc0, k)
431  CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
432  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  CALL cp_fm_create(matrix_shc0, matrix_c0%matrix_struct, name="shc0")
436 
437  SELECT CASE (preconditioner_env%cholesky_use)
438  CASE (cholesky_inverse)
439 ! if cho inverse
440  CALL cp_fm_to_fm(matrix_hc0, matrix_shc0)
441  CALL cp_fm_triangular_multiply(ortho, matrix_shc0, side="L", transpose_tr=.true., &
442  invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=k, alpha=1.0_dp)
443  CASE (cholesky_reduce)
444  CALL cp_fm_cholesky_restore(matrix_hc0, k, ortho, matrix_shc0, "SOLVE", transa="T")
445  CASE DEFAULT
446  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  para_env=preconditioner_env%para_env)
451  CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
452  CALL cp_fm_struct_release(fm_struct_tmp)
453  ! since we only use diagonal elements this is a bit of a waste
454  CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_shc0, matrix_shc0, 0.0_dp, matrix_s1)
455  ALLOCATE (diag(k))
456  CALL cp_fm_get_diag(matrix_s1, diag)
457  error_estimate = maxval(sqrt(abs(diag - c0_evals**2)))
458  DEALLOCATE (diag)
459  CALL cp_fm_release(matrix_s1)
460  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  preconditioner_env%energy_gap = max(energy_gap, error_estimate*fudge_factor)
466  CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
467  matrix_pre => preconditioner_env%fm
468  CALL cp_fm_upper_to_full(matrix_tmp, matrix_pre)
469  ! tmp = H ( 1 - PS )
470  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  para_env=preconditioner_env%para_env)
475  CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
476  CALL cp_fm_struct_release(fm_struct_tmp)
477  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  CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
480  CALL cp_fm_release(matrix_left)
481 
482  ALLOCATE (shifted_evals(k))
483  lambda = lambda_base + error_estimate
484  shifted_evals = c0_evals - lambda
485  CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
486  CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
487  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  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  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  invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
496  CASE (cholesky_reduce)
497  CALL cp_fm_cholesky_reduce(matrix_tmp, ortho)
498  END SELECT
499  CALL choose_eigv_solver(matrix_tmp, matrix_pre, preconditioner_env%full_evals)
500  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  invert_tr=.false., uplo_tr="U", n_rows=n, n_cols=n, alpha=1.0_dp)
504  CALL cp_fm_to_fm(matrix_pre, matrix_tmp)
505  CASE (cholesky_reduce)
506  CALL cp_fm_cholesky_restore(matrix_pre, n, ortho, matrix_tmp, "SOLVE")
507  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  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  preconditioner_env%full_evals(1:k) = c0_evals
531  CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
532 
533  CALL cp_fm_release(matrix_sc0)
534  CALL cp_fm_release(matrix_hc0)
535  CALL cp_fm_release(ortho)
536  CALL cp_fm_release(matrix_tmp)
537  DEALLOCATE (shifted_evals)
538  CALL timestop(handle)
539 
540  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  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  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  CALL timeset(routinen, handle)
571 
572  IF (ASSOCIATED(preconditioner_env%fm)) THEN
573  CALL cp_fm_release(preconditioner_env%fm)
574  DEALLOCATE (preconditioner_env%fm)
575  NULLIFY (preconditioner_env%fm)
576  END IF
577  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  para_env=preconditioner_env%para_env)
581  ALLOCATE (preconditioner_env%fm)
582  CALL cp_fm_create(preconditioner_env%fm, fm_struct_tmp, name="preconditioner_env%fm")
583  CALL cp_fm_create(matrix_tmp, fm_struct_tmp, name="matrix_tmp")
584  CALL cp_fm_struct_release(fm_struct_tmp)
585  ALLOCATE (preconditioner_env%full_evals(n))
586  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  CALL cp_fm_create(matrix_sc0, matrix_c0%matrix_struct, name="sc0")
595  CALL cp_fm_to_fm(matrix_c0, matrix_sc0)
596  CALL cp_fm_create(matrix_hc0, matrix_c0%matrix_struct, name="hc0")
597  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  para_env=preconditioner_env%para_env)
603  CALL cp_fm_create(matrix_s1, fm_struct_tmp, name="matrix_s1")
604  CALL cp_fm_struct_release(fm_struct_tmp)
605  ! since we only use diagonal elements this is a bit of a waste
606  CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, matrix_hc0, matrix_hc0, 0.0_dp, matrix_s1)
607  ALLOCATE (diag(k))
608  CALL cp_fm_get_diag(matrix_s1, diag)
609  error_estimate = maxval(sqrt(abs(diag - c0_evals**2)))
610  DEALLOCATE (diag)
611  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  preconditioner_env%energy_gap = max(energy_gap, error_estimate*fudge_factor)
617 
618  matrix_pre => preconditioner_env%fm
619  CALL copy_dbcsr_to_fm(matrix_h, matrix_tmp)
620  CALL cp_fm_upper_to_full(matrix_tmp, matrix_pre)
621  ! tmp = H ( 1 - PS )
622  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  para_env=preconditioner_env%para_env)
627  CALL cp_fm_create(matrix_left, fm_struct_tmp, name="matrix_left")
628  CALL cp_fm_struct_release(fm_struct_tmp)
629  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  CALL parallel_gemm('N', 'N', n, n, k, -1.0_dp, matrix_sc0, matrix_left, 1.0_dp, matrix_tmp)
632  CALL cp_fm_release(matrix_left)
633 
634  ALLOCATE (shifted_evals(k))
635  lambda = lambda_base + error_estimate
636  shifted_evals = c0_evals - lambda
637  CALL cp_fm_to_fm(matrix_sc0, matrix_hc0)
638  CALL cp_fm_column_scale(matrix_hc0, shifted_evals)
639  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  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  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  preconditioner_env%full_evals(1:k) = c0_evals
667  CALL cp_fm_to_fm(matrix_c0, matrix_pre, k, 1, 1)
668 
669  CALL cp_fm_release(matrix_sc0)
670  CALL cp_fm_release(matrix_hc0)
671  CALL cp_fm_release(matrix_tmp)
672  DEALLOCATE (shifted_evals)
673 
674  CALL timestop(handle)
675 
676  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  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices
708  TYPE(dbcsr_type), TARGET :: dbcsr_cthc, dbcsr_hc, dbcsr_sc, mo_dbcsr
709 
710  CALL timeset(routinen, handle)
711 
712  ! Allocate all working matrices needed
713  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  CALL cp_fm_to_dbcsr_row_template(mo_dbcsr, matrix_c0, matrix_h)
717  CALL dbcsr_create(dbcsr_sc, template=mo_dbcsr)
718  CALL dbcsr_create(dbcsr_hc, template=mo_dbcsr)
719  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  IF (.NOT. ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
723  ALLOCATE (preconditioner_env%sparse_matrix)
724  END IF
725 
726  ! Put the first term of the preconditioner (H) into the output matrix
727  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  IF (PRESENT(matrix_s)) THEN
732  CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, mo_dbcsr, 0.0_dp, dbcsr_sc)
733  ELSE
734  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  CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_h, mo_dbcsr, 0.0_dp, dbcsr_hc)
740  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  ALLOCATE (matrices(1))
744  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  nval_request=1, nrestarts=8, generalized_ev=.false., iram=.false.)
747  IF (ASSOCIATED(preconditioner_env%max_ev_vector)) &
748  CALL set_arnoldi_initial_vector(my_arnoldi, preconditioner_env%max_ev_vector)
749  CALL arnoldi_ev(matrices, my_arnoldi)
750  max_ev = real(get_selected_ritz_val(my_arnoldi, 1), dp)
751 
752  ! save the ev as guess for the next time
753  IF (.NOT. ASSOCIATED(preconditioner_env%max_ev_vector)) ALLOCATE (preconditioner_env%max_ev_vector)
754  CALL get_selected_ritz_vector(my_arnoldi, 1, matrices(1)%matrix, preconditioner_env%max_ev_vector)
755  CALL deallocate_arnoldi_data(my_arnoldi)
756  DEALLOCATE (matrices)
757 
758  ! Lets shift the occupied states a bit further up, -1.0 because we gonna subtract it from H
759  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  CALL dbcsr_multiply("N", "N", 2.0_dp, dbcsr_sc, dbcsr_cthc, 0.0_dp, dbcsr_hc)
762  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  IF (PRESENT(matrix_s)) THEN
767  ALLOCATE (matrices(2))
768  matrices(1)%matrix => preconditioner_env%sparse_matrix
769  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  nval_request=1, nrestarts=21, generalized_ev=.true., iram=.false.)
772  ELSE
773  ALLOCATE (matrices(1))
774  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  nval_request=1, nrestarts=8, generalized_ev=.false., iram=.false.)
777  END IF
778  IF (ASSOCIATED(preconditioner_env%min_ev_vector)) &
779  CALL set_arnoldi_initial_vector(my_arnoldi, preconditioner_env%min_ev_vector)
780 
781  ! compute the LUMO energy
782  CALL arnoldi_ev(matrices, my_arnoldi)
783  min_ev = real(get_selected_ritz_val(my_arnoldi, 1), dp)
784 
785  ! save the lumo vector for restarting in the next step
786  IF (.NOT. ASSOCIATED(preconditioner_env%min_ev_vector)) ALLOCATE (preconditioner_env%min_ev_vector)
787  CALL get_selected_ritz_vector(my_arnoldi, 1, matrices(1)%matrix, preconditioner_env%min_ev_vector)
788  CALL deallocate_arnoldi_data(my_arnoldi)
789  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  pre_shift = max(1.5_dp*(min_ev - max_ev), energy_gap)
795  IF (min_ev .LT. pre_shift) THEN
796  pre_shift = pre_shift - min_ev
797  ELSE
798  pre_shift = 0.0_dp
799  END IF
800  IF (PRESENT(matrix_s)) THEN
801  CALL dbcsr_add(preconditioner_env%sparse_matrix, matrix_s, 1.0_dp, pre_shift)
802  ELSE
803  CALL dbcsr_add_on_diag(preconditioner_env%sparse_matrix, pre_shift)
804  END IF
805 
806  CALL dbcsr_release(mo_dbcsr)
807  CALL dbcsr_release(dbcsr_hc)
808  CALL dbcsr_release(dbcsr_sc)
809  CALL dbcsr_release(dbcsr_cthc)
810 
811  CALL timestop(handle)
812 
813  END SUBROUTINE make_full_single_inverse
814 
815 END MODULE preconditioner_makes
816 
arnoldi_api
arnoldi iteration using dbcsr
Definition: arnoldi_api.F:15
arnoldi_api::arnoldi_ev
subroutine, public arnoldi_ev(matrix, arnoldi_data)
Driver routine for different arnoldi eigenvalue methods the selection which one is to be taken is mad...
Definition: arnoldi_api.F:69
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:744
cp_dbcsr_operations::cp_dbcsr_m_by_n_from_template
subroutine, public cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym, data_type)
Utility function to create an arbitrary shaped dbcsr matrix with the same processor grid as the templ...
Definition: cp_dbcsr_operations.F:1105
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:208
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:1070
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:1816
cp_fm_basic_linalg::cp_fm_upper_to_full
subroutine, public cp_fm_upper_to_full(matrix, work)
given an upper triangular matrix computes the corresponding full matrix
Definition: cp_fm_basic_linalg.F:1693
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:2127
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:1526
cp_fm_cholesky
various cholesky decomposition related routines
Definition: cp_fm_cholesky.F:14
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:50
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:192
cp_fm_cholesky::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa)
...
Definition: cp_fm_cholesky.F:254
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:208
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:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
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:570
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:1016
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::ot_precond_full_kinetic
integer, parameter, public ot_precond_full_kinetic
Definition: input_constants.F:513
input_constants::cholesky_reduce
integer, parameter, public cholesky_reduce
Definition: input_constants.F:195
input_constants::cholesky_inverse
integer, parameter, public cholesky_inverse
Definition: input_constants.F:195
input_constants::ot_precond_solver_default
integer, parameter, public ot_precond_solver_default
Definition: input_constants.F:520
input_constants::ot_precond_full_single
integer, parameter, public ot_precond_full_single
Definition: input_constants.F:513
input_constants::ot_precond_solver_inv_chol
integer, parameter, public ot_precond_solver_inv_chol
Definition: input_constants.F:520
input_constants::ot_precond_full_single_inverse
integer, parameter, public ot_precond_full_single_inverse
Definition: input_constants.F:513
input_constants::ot_precond_s_inverse
integer, parameter, public ot_precond_s_inverse
Definition: input_constants.F:513
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition: input_constants.F:513
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:81
preconditioner_types
types of preconditioners
Definition: preconditioner_types.F:14