(git:e7e05ae)
qs_scf_block_davidson.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 module that contains the algorithms to perform an itrative
10 !> diagonalization by the block-Davidson approach
11 !> P. Blaha, et al J. Comp. Physics, 229, (2010), 453-460
12 !> Iterative diagonalization in augmented plane wave based
13 !> methods in electronic structure calculations
14 !> \par History
15 !> 05.2011 created [MI]
16 !> \author MI
17 ! **************************************************************************************************
18 MODULE qs_scf_block_davidson
19 
20  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
21  copy_fm_to_dbcsr,&
22  cp_dbcsr_m_by_n_from_row_template,&
23  cp_dbcsr_m_by_n_from_template,&
24  cp_dbcsr_sm_fm_multiply
25  USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
26  cp_fm_scale_and_add,&
27  cp_fm_symm,&
28  cp_fm_transpose,&
29  cp_fm_triangular_invert,&
30  cp_fm_upper_to_full
31  USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
32  cp_fm_cholesky_restore
33  USE cp_fm_diag, ONLY: choose_eigv_solver,&
34  cp_fm_power
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_release,&
41  cp_fm_set_all,&
42  cp_fm_to_fm,&
43  cp_fm_to_fm_submat,&
44  cp_fm_type,&
45  cp_fm_vectorsnorm
46  USE dbcsr_api, ONLY: &
47  dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_get_diag, dbcsr_get_info, dbcsr_init_p, &
48  dbcsr_multiply, dbcsr_norm, dbcsr_norm_column, dbcsr_release_p, dbcsr_scale_by_vector, &
49  dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_real_default, dbcsr_type_symmetric
50  USE kinds, ONLY: dp
51  USE machine, ONLY: m_walltime
52  USE parallel_gemm_api, ONLY: parallel_gemm
53  USE preconditioner, ONLY: apply_preconditioner
54  USE preconditioner_types, ONLY: preconditioner_type
55  USE qs_block_davidson_types, ONLY: davidson_type
56  USE qs_mo_types, ONLY: get_mo_set,&
57  mo_set_type
58 #include "./base/base_uses.f90"
59 
60  IMPLICIT NONE
61  PRIVATE
62  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_block_davidson'
63 
64  PUBLIC :: generate_extended_space, generate_extended_space_sparse
65 
66 CONTAINS
67 
68 ! **************************************************************************************************
69 !> \brief ...
70 !> \param bdav_env ...
71 !> \param mo_set ...
72 !> \param matrix_h ...
73 !> \param matrix_s ...
74 !> \param output_unit ...
75 !> \param preconditioner ...
76 ! **************************************************************************************************
77  SUBROUTINE generate_extended_space(bdav_env, mo_set, matrix_h, matrix_s, output_unit, &
78  preconditioner)
79 
80  TYPE(davidson_type) :: bdav_env
81  TYPE(mo_set_type), INTENT(IN) :: mo_set
82  TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
83  INTEGER, INTENT(IN) :: output_unit
84  TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner
85 
86  CHARACTER(len=*), PARAMETER :: routinen = 'generate_extended_space'
87 
88  INTEGER :: handle, homo, i_first, i_last, imo, iter, j, jj, max_iter, n, nao, nmat, nmat2, &
89  nmo, nmo_converged, nmo_not_converged, nset, nset_conv, nset_not_conv
90  INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv
91  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iconv_set, inotconv_set
92  LOGICAL :: converged, do_apply_preconditioner
93  REAL(dp) :: lambda, max_norm, min_norm, t1, t2
94  REAL(dp), ALLOCATABLE, DIMENSION(:) :: ritz_coeff, vnorm
95  REAL(dp), DIMENSION(:), POINTER :: eig_not_conv, eigenvalues, evals
96  TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
97  TYPE(cp_fm_type) :: c_conv, c_notconv, c_out, h_block, h_fm, &
98  m_hc, m_sc, m_tmp, mt_tmp, s_block, &
99  s_fm, v_block, w_block
100  TYPE(cp_fm_type), POINTER :: c_pz, c_z, mo_coeff
101  TYPE(dbcsr_type), POINTER :: mo_coeff_b
102 
103  CALL timeset(routinen, handle)
104 
105  NULLIFY (mo_coeff, mo_coeff_b, eigenvalues)
106 
107  do_apply_preconditioner = .false.
108  IF (PRESENT(preconditioner)) do_apply_preconditioner = .true.
109  CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, eigenvalues=eigenvalues, &
110  nao=nao, nmo=nmo, homo=homo)
111  IF (do_apply_preconditioner) THEN
112  max_iter = bdav_env%max_iter
113  ELSE
114  max_iter = 1
115  END IF
116 
117  NULLIFY (c_z, c_pz)
118  NULLIFY (evals, eig_not_conv)
119  t1 = m_walltime()
120  IF (output_unit > 0) THEN
121  WRITE (output_unit, "(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
122  " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", repeat("-", 60)
123  END IF
124 
125  ALLOCATE (iconv(nmo))
126  ALLOCATE (inotconv(nmo))
127  ALLOCATE (ritz_coeff(nmo))
128  ALLOCATE (vnorm(nmo))
129 
130  converged = .false.
131  DO iter = 1, max_iter
132 
133  ! compute Ritz values
134  ritz_coeff = 0.0_dp
135  CALL cp_fm_create(m_hc, mo_coeff%matrix_struct, name="hc")
136  CALL cp_dbcsr_sm_fm_multiply(matrix_h, mo_coeff, m_hc, nmo)
137  CALL cp_fm_create(m_sc, mo_coeff%matrix_struct, name="sc")
138  CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_coeff, m_sc, nmo)
139 
140  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=nmo, &
141  context=mo_coeff%matrix_struct%context, &
142  para_env=mo_coeff%matrix_struct%para_env)
143  CALL cp_fm_create(m_tmp, fm_struct_tmp, name="matrix_tmp")
144  CALL cp_fm_struct_release(fm_struct_tmp)
145 
146  CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, mo_coeff, m_hc, 0.0_dp, m_tmp)
147  CALL cp_fm_get_diag(m_tmp, ritz_coeff)
148  CALL cp_fm_release(m_tmp)
149 
150  ! Check for converged eigenvectors
151  c_z => bdav_env%matrix_z
152  c_pz => bdav_env%matrix_pz
153  CALL cp_fm_to_fm(m_sc, c_z)
154  CALL cp_fm_column_scale(c_z, ritz_coeff)
155  CALL cp_fm_scale_and_add(-1.0_dp, c_z, 1.0_dp, m_hc)
156  CALL cp_fm_vectorsnorm(c_z, vnorm)
157 
158  nmo_converged = 0
159  nmo_not_converged = 0
160  max_norm = 0.0_dp
161  min_norm = 1.e10_dp
162  DO imo = 1, nmo
163  max_norm = max(max_norm, vnorm(imo))
164  min_norm = min(min_norm, vnorm(imo))
165  END DO
166  iconv = 0
167  inotconv = 0
168  DO imo = 1, nmo
169  IF (vnorm(imo) <= bdav_env%eps_iter) THEN
170  nmo_converged = nmo_converged + 1
171  iconv(nmo_converged) = imo
172  ELSE
173  nmo_not_converged = nmo_not_converged + 1
174  inotconv(nmo_not_converged) = imo
175  END IF
176  END DO
177 
178  IF (nmo_converged > 0) THEN
179  ALLOCATE (iconv_set(nmo_converged, 2))
180  ALLOCATE (inotconv_set(nmo_not_converged, 2))
181  i_last = iconv(1)
182  nset = 0
183  DO j = 1, nmo_converged
184  imo = iconv(j)
185 
186  IF (imo == i_last + 1) THEN
187  i_last = imo
188  iconv_set(nset, 2) = imo
189  ELSE
190  i_last = imo
191  nset = nset + 1
192  iconv_set(nset, 1) = imo
193  iconv_set(nset, 2) = imo
194  END IF
195  END DO
196  nset_conv = nset
197 
198  i_last = inotconv(1)
199  nset = 0
200  DO j = 1, nmo_not_converged
201  imo = inotconv(j)
202 
203  IF (imo == i_last + 1) THEN
204  i_last = imo
205  inotconv_set(nset, 2) = imo
206  ELSE
207  i_last = imo
208  nset = nset + 1
209  inotconv_set(nset, 1) = imo
210  inotconv_set(nset, 2) = imo
211  END IF
212  END DO
213  nset_not_conv = nset
214  CALL cp_fm_release(m_sc)
215  CALL cp_fm_release(m_hc)
216  NULLIFY (c_z, c_pz)
217  END IF
218 
219  IF (real(nmo_converged, dp)/real(nmo, dp) > bdav_env%conv_percent) THEN
220  converged = .true.
221  DEALLOCATE (iconv_set)
222  DEALLOCATE (inotconv_set)
223  t2 = m_walltime()
224  IF (output_unit > 0) THEN
225  WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
226  iter, nmo_converged, max_norm, min_norm, t2 - t1
227 
228  WRITE (output_unit, *) " Reached convergence in ", iter, &
229  " Davidson iterations"
230  END IF
231 
232  EXIT
233  END IF
234 
235  IF (nmo_converged > 0) THEN
236  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nao, &
237  context=mo_coeff%matrix_struct%context, &
238  para_env=mo_coeff%matrix_struct%para_env)
239  !allocate h_fm
240  CALL cp_fm_create(h_fm, fm_struct_tmp, name="matrix_tmp")
241  !allocate s_fm
242  CALL cp_fm_create(s_fm, fm_struct_tmp, name="matrix_tmp")
243  !copy matrix_h in h_fm
244  CALL copy_dbcsr_to_fm(matrix_h, h_fm)
245  CALL cp_fm_upper_to_full(h_fm, s_fm)
246 
247  !copy matrix_s in s_fm
248 ! CALL cp_fm_set_all(s_fm,0.0_dp)
249  CALL copy_dbcsr_to_fm(matrix_s, s_fm)
250 
251  !allocate c_out
252  CALL cp_fm_create(c_out, fm_struct_tmp, name="matrix_tmp")
253  ! set c_out to zero
254  CALL cp_fm_set_all(c_out, 0.0_dp)
255  CALL cp_fm_struct_release(fm_struct_tmp)
256 
257  !allocate c_conv
258  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_converged, &
259  context=mo_coeff%matrix_struct%context, &
260  para_env=mo_coeff%matrix_struct%para_env)
261  CALL cp_fm_create(c_conv, fm_struct_tmp, name="c_conv")
262  CALL cp_fm_set_all(c_conv, 0.0_dp)
263  !allocate m_tmp
264  CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_nxmc")
265  CALL cp_fm_struct_release(fm_struct_tmp)
266  END IF
267 
268  !allocate c_notconv
269  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_not_converged, &
270  context=mo_coeff%matrix_struct%context, &
271  para_env=mo_coeff%matrix_struct%para_env)
272  CALL cp_fm_create(c_notconv, fm_struct_tmp, name="c_notconv")
273  CALL cp_fm_set_all(c_notconv, 0.0_dp)
274  IF (nmo_converged > 0) THEN
275  CALL cp_fm_create(m_hc, fm_struct_tmp, name="m_hc")
276  CALL cp_fm_create(m_sc, fm_struct_tmp, name="m_sc")
277  !allocate c_z
278  ALLOCATE (c_z, c_pz)
279  CALL cp_fm_create(c_z, fm_struct_tmp, name="c_z")
280  CALL cp_fm_create(c_pz, fm_struct_tmp, name="c_pz")
281  CALL cp_fm_set_all(c_z, 0.0_dp)
282 
283  ! sum contributions to c_out
284  jj = 1
285  DO j = 1, nset_conv
286  i_first = iconv_set(j, 1)
287  i_last = iconv_set(j, 2)
288  n = i_last - i_first + 1
289  CALL cp_fm_to_fm_submat(mo_coeff, c_conv, nao, n, 1, i_first, 1, jj)
290  jj = jj + n
291  END DO
292  CALL cp_fm_symm('L', 'U', nao, nmo_converged, 1.0_dp, s_fm, c_conv, 0.0_dp, m_tmp)
293  CALL parallel_gemm('N', 'T', nao, nao, nmo_converged, 1.0_dp, m_tmp, m_tmp, 0.0_dp, c_out)
294 
295  ! project c_out out of H
296  lambda = 100.0_dp*abs(eigenvalues(homo))
297  CALL cp_fm_scale_and_add(lambda, c_out, 1.0_dp, h_fm)
298  CALL cp_fm_release(m_tmp)
299  CALL cp_fm_release(h_fm)
300 
301  END IF
302 
303  !allocate m_tmp
304  CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_nxm")
305  CALL cp_fm_struct_release(fm_struct_tmp)
306  IF (nmo_converged > 0) THEN
307  ALLOCATE (eig_not_conv(nmo_not_converged))
308  jj = 1
309  DO j = 1, nset_not_conv
310  i_first = inotconv_set(j, 1)
311  i_last = inotconv_set(j, 2)
312  n = i_last - i_first + 1
313  CALL cp_fm_to_fm_submat(mo_coeff, c_notconv, nao, n, 1, i_first, 1, jj)
314  eig_not_conv(jj:jj + n - 1) = ritz_coeff(i_first:i_last)
315  jj = jj + n
316  END DO
317  CALL parallel_gemm('N', 'N', nao, nmo_not_converged, nao, 1.0_dp, c_out, c_notconv, 0.0_dp, m_hc)
318  CALL cp_fm_symm('L', 'U', nao, nmo_not_converged, 1.0_dp, s_fm, c_notconv, 0.0_dp, m_sc)
319  ! extend suspace using only the not converged vectors
320  CALL cp_fm_to_fm(m_sc, m_tmp)
321  CALL cp_fm_column_scale(m_tmp, eig_not_conv)
322  CALL cp_fm_scale_and_add(-1.0_dp, m_tmp, 1.0_dp, m_hc)
323  DEALLOCATE (eig_not_conv)
324  CALL cp_fm_to_fm(m_tmp, c_z)
325  ELSE
326  CALL cp_fm_to_fm(mo_coeff, c_notconv)
327  END IF
328 
329  !preconditioner
330  IF (do_apply_preconditioner) THEN
331  IF (preconditioner%in_use /= 0) THEN
332  CALL apply_preconditioner(preconditioner, c_z, c_pz)
333  ELSE
334  CALL cp_fm_to_fm(c_z, c_pz)
335  END IF
336  ELSE
337  CALL cp_fm_to_fm(c_z, c_pz)
338  END IF
339  CALL cp_fm_release(m_tmp)
340 
341  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo_not_converged, ncol_global=nmo_not_converged, &
342  context=mo_coeff%matrix_struct%context, &
343  para_env=mo_coeff%matrix_struct%para_env)
344 
345  CALL cp_fm_create(m_tmp, fm_struct_tmp, name="m_tmp_mxm")
346  CALL cp_fm_create(mt_tmp, fm_struct_tmp, name="mt_tmp_mxm")
347  CALL cp_fm_struct_release(fm_struct_tmp)
348 
349  nmat = nmo_not_converged
350  nmat2 = 2*nmo_not_converged
351  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmat2, ncol_global=nmat2, &
352  context=mo_coeff%matrix_struct%context, &
353  para_env=mo_coeff%matrix_struct%para_env)
354 
355  CALL cp_fm_create(s_block, fm_struct_tmp, name="sb")
356  CALL cp_fm_create(h_block, fm_struct_tmp, name="hb")
357  CALL cp_fm_create(v_block, fm_struct_tmp, name="vb")
358  CALL cp_fm_create(w_block, fm_struct_tmp, name="wb")
359  ALLOCATE (evals(nmat2))
360 
361  CALL cp_fm_struct_release(fm_struct_tmp)
362 
363  ! compute CSC
364  CALL cp_fm_set_all(s_block, 0.0_dp, 1.0_dp)
365 
366  ! compute CHC
367  CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_notconv, m_hc, 0.0_dp, m_tmp)
368  CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1, 1)
369 
370  ! compute ZSC
371  CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_sc, 0.0_dp, m_tmp)
372  CALL cp_fm_to_fm_submat(m_tmp, s_block, nmat, nmat, 1, 1, 1 + nmat, 1)
373  CALL cp_fm_transpose(m_tmp, mt_tmp)
374  CALL cp_fm_to_fm_submat(mt_tmp, s_block, nmat, nmat, 1, 1, 1, 1 + nmat)
375  ! compute ZHC
376  CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_hc, 0.0_dp, m_tmp)
377  CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1 + nmat, 1)
378  CALL cp_fm_transpose(m_tmp, mt_tmp)
379  CALL cp_fm_to_fm_submat(mt_tmp, h_block, nmat, nmat, 1, 1, 1, 1 + nmat)
380 
381  CALL cp_fm_release(mt_tmp)
382 
383  ! reuse m_sc and m_hc to computr HZ and SZ
384  IF (nmo_converged > 0) THEN
385  CALL parallel_gemm('N', 'N', nao, nmat, nao, 1.0_dp, c_out, c_pz, 0.0_dp, m_hc)
386  CALL cp_fm_symm('L', 'U', nao, nmo_not_converged, 1.0_dp, s_fm, c_pz, 0.0_dp, m_sc)
387 
388  CALL cp_fm_release(c_out)
389  CALL cp_fm_release(c_conv)
390  CALL cp_fm_release(s_fm)
391  ELSE
392  CALL cp_dbcsr_sm_fm_multiply(matrix_h, c_pz, m_hc, nmo)
393  CALL cp_dbcsr_sm_fm_multiply(matrix_s, c_pz, m_sc, nmo)
394  END IF
395 
396  ! compute ZSZ
397  CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_sc, 0.0_dp, m_tmp)
398  CALL cp_fm_to_fm_submat(m_tmp, s_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
399  ! compute ZHZ
400  CALL parallel_gemm('T', 'N', nmat, nmat, nao, 1.0_dp, c_pz, m_hc, 0.0_dp, m_tmp)
401  CALL cp_fm_to_fm_submat(m_tmp, h_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
402 
403  CALL cp_fm_release(m_sc)
404 
405  ! solution of the reduced eigenvalues problem
406  CALL reduce_extended_space(s_block, h_block, v_block, w_block, evals, nmat2)
407 
408  ! extract egenvectors
409  CALL cp_fm_to_fm_submat(v_block, m_tmp, nmat, nmat, 1, 1, 1, 1)
410  CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, c_notconv, m_tmp, 0.0_dp, m_hc)
411  CALL cp_fm_to_fm_submat(v_block, m_tmp, nmat, nmat, 1 + nmat, 1, 1, 1)
412  CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, c_pz, m_tmp, 1.0_dp, m_hc)
413 
414  CALL cp_fm_release(m_tmp)
415 
416  CALL cp_fm_release(c_notconv)
417  CALL cp_fm_release(s_block)
418  CALL cp_fm_release(h_block)
419  CALL cp_fm_release(w_block)
420  CALL cp_fm_release(v_block)
421 
422  IF (nmo_converged > 0) THEN
423  CALL cp_fm_release(c_z)
424  CALL cp_fm_release(c_pz)
425  DEALLOCATE (c_z, c_pz)
426  jj = 1
427  DO j = 1, nset_not_conv
428  i_first = inotconv_set(j, 1)
429  i_last = inotconv_set(j, 2)
430  n = i_last - i_first + 1
431  CALL cp_fm_to_fm_submat(m_hc, mo_coeff, nao, n, 1, jj, 1, i_first)
432  eigenvalues(i_first:i_last) = evals(jj:jj + n - 1)
433  jj = jj + n
434  END DO
435  DEALLOCATE (iconv_set)
436  DEALLOCATE (inotconv_set)
437  ELSE
438  CALL cp_fm_to_fm(m_hc, mo_coeff)
439  eigenvalues(1:nmo) = evals(1:nmo)
440  END IF
441  DEALLOCATE (evals)
442 
443  CALL cp_fm_release(m_hc)
444 
445  CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b) !fm->dbcsr
446 
447  t2 = m_walltime()
448  IF (output_unit > 0) THEN
449  WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
450  iter, nmo_converged, max_norm, min_norm, t2 - t1
451  END IF
452  t1 = m_walltime()
453 
454  END DO ! iter
455 
456  DEALLOCATE (iconv)
457  DEALLOCATE (inotconv)
458  DEALLOCATE (ritz_coeff)
459  DEALLOCATE (vnorm)
460 
461  CALL timestop(handle)
462  END SUBROUTINE generate_extended_space
463 
464 ! **************************************************************************************************
465 !> \brief ...
466 !> \param bdav_env ...
467 !> \param mo_set ...
468 !> \param matrix_h ...
469 !> \param matrix_s ...
470 !> \param output_unit ...
471 !> \param preconditioner ...
472 ! **************************************************************************************************
473  SUBROUTINE generate_extended_space_sparse(bdav_env, mo_set, matrix_h, matrix_s, output_unit, &
474  preconditioner)
475 
476  TYPE(davidson_type) :: bdav_env
477  TYPE(mo_set_type), INTENT(IN) :: mo_set
478  TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
479  INTEGER, INTENT(IN) :: output_unit
480  TYPE(preconditioner_type), OPTIONAL, POINTER :: preconditioner
481 
482  CHARACTER(len=*), PARAMETER :: routinen = 'generate_extended_space_sparse'
483 
484  INTEGER :: handle, homo, i_first, i_last, imo, iter, j, jj, k, max_iter, n, nao, nmat, &
485  nmat2, nmo, nmo_converged, nmo_not_converged, nset, nset_conv, nset_not_conv
486  INTEGER, ALLOCATABLE, DIMENSION(:) :: iconv, inotconv
487  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iconv_set, inotconv_set
488  LOGICAL :: converged, do_apply_preconditioner
489  REAL(dp) :: lambda, max_norm, min_norm, t1, t2
490  REAL(dp), ALLOCATABLE, DIMENSION(:) :: eig_not_conv, evals, ritz_coeff, vnorm
491  REAL(dp), DIMENSION(:), POINTER :: eigenvalues
492  TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
493  TYPE(cp_fm_type) :: h_block, matrix_mm_fm, matrix_mmt_fm, &
494  matrix_nm_fm, matrix_z_fm, mo_conv_fm, &
495  s_block, v_block, w_block
496  TYPE(cp_fm_type), POINTER :: mo_coeff, mo_notconv_fm
497  TYPE(dbcsr_type), POINTER :: c_out, matrix_hc, matrix_mm, matrix_pz, &
498  matrix_sc, matrix_z, mo_coeff_b, &
499  mo_conv, mo_notconv, smo_conv
500 
501  CALL timeset(routinen, handle)
502 
503  do_apply_preconditioner = .false.
504  IF (PRESENT(preconditioner)) do_apply_preconditioner = .true.
505 
506  NULLIFY (mo_coeff, mo_coeff_b, matrix_hc, matrix_sc, matrix_z, matrix_pz, matrix_mm)
507  NULLIFY (mo_notconv_fm, mo_conv, mo_notconv, smo_conv, c_out)
508  NULLIFY (fm_struct_tmp)
509  CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, &
510  eigenvalues=eigenvalues, homo=homo, nao=nao, nmo=nmo)
511  IF (do_apply_preconditioner) THEN
512  max_iter = bdav_env%max_iter
513  ELSE
514  max_iter = 1
515  END IF
516 
517  t1 = m_walltime()
518  IF (output_unit > 0) THEN
519  WRITE (output_unit, "(T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
520  " Cycle ", " conv. MOS ", " B2MAX ", " B2MIN ", " Time", repeat("-", 60)
521  END IF
522 
523  ! Allocate array for Ritz values
524  ALLOCATE (ritz_coeff(nmo))
525  ALLOCATE (iconv(nmo))
526  ALLOCATE (inotconv(nmo))
527  ALLOCATE (vnorm(nmo))
528 
529  converged = .false.
530  DO iter = 1, max_iter
531  NULLIFY (c_out, mo_conv, mo_notconv_fm, mo_notconv)
532  ! Prepare HC and SC, using mo_coeff_b (sparse), these are still sparse
533  CALL dbcsr_init_p(matrix_hc)
534  CALL dbcsr_create(matrix_hc, template=mo_coeff_b, &
535  name="matrix_hc", &
536  matrix_type=dbcsr_type_no_symmetry, &
537  nze=0, data_type=dbcsr_type_real_default)
538  CALL dbcsr_init_p(matrix_sc)
539  CALL dbcsr_create(matrix_sc, template=mo_coeff_b, &
540  name="matrix_sc", &
541  matrix_type=dbcsr_type_no_symmetry, &
542  nze=0, data_type=dbcsr_type_real_default)
543 
544  CALL dbcsr_get_info(mo_coeff_b, nfullrows_total=n, nfullcols_total=k)
545  CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_h, mo_coeff_b, 0.0_dp, matrix_hc, last_column=k)
546  CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_coeff_b, 0.0_dp, matrix_sc, last_column=k)
547 
548  ! compute Ritz values
549  ritz_coeff = 0.0_dp
550  ! Allocate Sparse matrices: nmoxnmo
551  ! matrix_mm
552 
553  CALL dbcsr_init_p(matrix_mm)
554  CALL cp_dbcsr_m_by_n_from_template(matrix_mm, template=matrix_s, m=nmo, n=nmo, &
555  sym=dbcsr_type_no_symmetry)
556 
557  CALL dbcsr_multiply('t', 'n', 1.0_dp, mo_coeff_b, matrix_hc, 0.0_dp, matrix_mm, last_column=k)
558  CALL dbcsr_get_diag(matrix_mm, ritz_coeff)
559  CALL mo_coeff%matrix_struct%para_env%sum(ritz_coeff)
560 
561  ! extended subspace P Z = P [H - theta S]C this ia another matrix of type and size as mo_coeff_b
562  CALL dbcsr_init_p(matrix_z)
563  CALL dbcsr_create(matrix_z, template=mo_coeff_b, &
564  name="matrix_z", &
565  matrix_type=dbcsr_type_no_symmetry, &
566  nze=0, data_type=dbcsr_type_real_default)
567  CALL dbcsr_copy(matrix_z, matrix_sc)
568  CALL dbcsr_scale_by_vector(matrix_z, ritz_coeff, side='right')
569  CALL dbcsr_add(matrix_z, matrix_hc, -1.0_dp, 1.0_dp)
570 
571  ! Check for converged eigenvectors
572  vnorm = 0.0_dp
573  CALL dbcsr_norm(matrix_z, which_norm=dbcsr_norm_column, norm_vector=vnorm)
574  nmo_converged = 0
575  nmo_not_converged = 0
576  max_norm = 0.0_dp
577  min_norm = 1.e10_dp
578  DO imo = 1, nmo
579  max_norm = max(max_norm, vnorm(imo))
580  min_norm = min(min_norm, vnorm(imo))
581  END DO
582  iconv = 0
583  inotconv = 0
584 
585  DO imo = 1, nmo
586  IF (vnorm(imo) <= bdav_env%eps_iter) THEN
587  nmo_converged = nmo_converged + 1
588  iconv(nmo_converged) = imo
589  ELSE
590  nmo_not_converged = nmo_not_converged + 1
591  inotconv(nmo_not_converged) = imo
592  END IF
593  END DO
594 
595  IF (nmo_converged > 0) THEN
596  ALLOCATE (iconv_set(nmo_converged, 2))
597  ALLOCATE (inotconv_set(nmo_not_converged, 2))
598  i_last = iconv(1)
599  nset = 0
600  DO j = 1, nmo_converged
601  imo = iconv(j)
602 
603  IF (imo == i_last + 1) THEN
604  i_last = imo
605  iconv_set(nset, 2) = imo
606  ELSE
607  i_last = imo
608  nset = nset + 1
609  iconv_set(nset, 1) = imo
610  iconv_set(nset, 2) = imo
611  END IF
612  END DO
613  nset_conv = nset
614 
615  i_last = inotconv(1)
616  nset = 0
617  DO j = 1, nmo_not_converged
618  imo = inotconv(j)
619 
620  IF (imo == i_last + 1) THEN
621  i_last = imo
622  inotconv_set(nset, 2) = imo
623  ELSE
624  i_last = imo
625  nset = nset + 1
626  inotconv_set(nset, 1) = imo
627  inotconv_set(nset, 2) = imo
628  END IF
629  END DO
630  nset_not_conv = nset
631 
632  CALL dbcsr_release_p(matrix_hc)
633  CALL dbcsr_release_p(matrix_sc)
634  CALL dbcsr_release_p(matrix_z)
635  CALL dbcsr_release_p(matrix_mm)
636  END IF
637 
638  IF (real(nmo_converged, dp)/real(nmo, dp) > bdav_env%conv_percent) THEN
639  DEALLOCATE (iconv_set)
640 
641  DEALLOCATE (inotconv_set)
642 
643  converged = .true.
644  t2 = m_walltime()
645  IF (output_unit > 0) THEN
646  WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
647  iter, nmo_converged, max_norm, min_norm, t2 - t1
648 
649  WRITE (output_unit, *) " Reached convergence in ", iter, &
650  " Davidson iterations"
651  END IF
652 
653  EXIT
654  END IF
655 
656  IF (nmo_converged > 0) THEN
657 
658  !allocate mo_conv_fm
659  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_converged, &
660  context=mo_coeff%matrix_struct%context, &
661  para_env=mo_coeff%matrix_struct%para_env)
662  CALL cp_fm_create(mo_conv_fm, fm_struct_tmp, name="mo_conv_fm")
663 
664  CALL cp_fm_struct_release(fm_struct_tmp)
665 
666  ! extract mo_conv from mo_coeff full matrix
667  jj = 1
668  DO j = 1, nset_conv
669  i_first = iconv_set(j, 1)
670  i_last = iconv_set(j, 2)
671  n = i_last - i_first + 1
672  CALL cp_fm_to_fm_submat(mo_coeff, mo_conv_fm, nao, n, 1, i_first, 1, jj)
673  jj = jj + n
674  END DO
675 
676  ! allocate c_out sparse matrix, to project out the converged MOS
677  CALL dbcsr_init_p(c_out)
678  CALL dbcsr_create(c_out, template=matrix_s, &
679  name="c_out", &
680  matrix_type=dbcsr_type_symmetric, &
681  nze=0, data_type=dbcsr_type_real_default)
682 
683  ! allocate mo_conv sparse
684  CALL dbcsr_init_p(mo_conv)
685  CALL cp_dbcsr_m_by_n_from_row_template(mo_conv, template=matrix_s, n=nmo_converged, &
686  sym=dbcsr_type_no_symmetry)
687 
688  CALL dbcsr_init_p(smo_conv)
689  CALL cp_dbcsr_m_by_n_from_row_template(smo_conv, template=matrix_s, n=nmo_converged, &
690  sym=dbcsr_type_no_symmetry)
691 
692  CALL copy_fm_to_dbcsr(mo_conv_fm, mo_conv) !fm->dbcsr
693 
694  CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_conv, 0.0_dp, smo_conv, last_column=nmo_converged)
695  CALL dbcsr_multiply('n', 't', 1.0_dp, smo_conv, smo_conv, 0.0_dp, c_out, last_column=nao)
696  ! project c_out out of H
697  lambda = 100.0_dp*abs(eigenvalues(homo))
698  CALL dbcsr_add(c_out, matrix_h, lambda, 1.0_dp)
699 
700  CALL dbcsr_release_p(mo_conv)
701  CALL dbcsr_release_p(smo_conv)
702  CALL cp_fm_release(mo_conv_fm)
703 
704  !allocate c_notconv_fm
705  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo_not_converged, &
706  context=mo_coeff%matrix_struct%context, &
707  para_env=mo_coeff%matrix_struct%para_env)
708  ALLOCATE (mo_notconv_fm)
709  CALL cp_fm_create(mo_notconv_fm, fm_struct_tmp, name="mo_notconv_fm")
710  CALL cp_fm_struct_release(fm_struct_tmp)
711 
712  ! extract mo_notconv from mo_coeff full matrix
713  ALLOCATE (eig_not_conv(nmo_not_converged))
714 
715  jj = 1
716  DO j = 1, nset_not_conv
717  i_first = inotconv_set(j, 1)
718  i_last = inotconv_set(j, 2)
719  n = i_last - i_first + 1
720  CALL cp_fm_to_fm_submat(mo_coeff, mo_notconv_fm, nao, n, 1, i_first, 1, jj)
721  eig_not_conv(jj:jj + n - 1) = ritz_coeff(i_first:i_last)
722  jj = jj + n
723  END DO
724 
725  ! allocate mo_conv sparse
726  CALL dbcsr_init_p(mo_notconv)
727  CALL cp_dbcsr_m_by_n_from_row_template(mo_notconv, template=matrix_s, n=nmo_not_converged, &
728  sym=dbcsr_type_no_symmetry)
729 
730  CALL dbcsr_init_p(matrix_hc)
731  CALL cp_dbcsr_m_by_n_from_row_template(matrix_hc, template=matrix_s, n=nmo_not_converged, &
732  sym=dbcsr_type_no_symmetry)
733 
734  CALL dbcsr_init_p(matrix_sc)
735  CALL cp_dbcsr_m_by_n_from_row_template(matrix_sc, template=matrix_s, n=nmo_not_converged, &
736  sym=dbcsr_type_no_symmetry)
737 
738  CALL dbcsr_init_p(matrix_z)
739  CALL cp_dbcsr_m_by_n_from_row_template(matrix_z, template=matrix_s, n=nmo_not_converged, &
740  sym=dbcsr_type_no_symmetry)
741 
742  CALL copy_fm_to_dbcsr(mo_notconv_fm, mo_notconv) !fm->dbcsr
743 
744  CALL dbcsr_multiply('n', 'n', 1.0_dp, c_out, mo_notconv, 0.0_dp, matrix_hc, &
745  last_column=nmo_not_converged)
746  CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, mo_notconv, 0.0_dp, matrix_sc, &
747  last_column=nmo_not_converged)
748 
749  CALL dbcsr_copy(matrix_z, matrix_sc)
750  CALL dbcsr_scale_by_vector(matrix_z, eig_not_conv, side='right')
751  CALL dbcsr_add(matrix_z, matrix_hc, -1.0_dp, 1.0_dp)
752 
753  DEALLOCATE (eig_not_conv)
754 
755  ! matrix_mm
756  CALL dbcsr_init_p(matrix_mm)
757  CALL cp_dbcsr_m_by_n_from_template(matrix_mm, template=matrix_s, m=nmo_not_converged, n=nmo_not_converged, &
758  sym=dbcsr_type_no_symmetry)
759 
760  CALL dbcsr_multiply('t', 'n', 1.0_dp, mo_notconv, matrix_hc, 0.0_dp, matrix_mm, &
761  last_column=nmo_not_converged)
762 
763  ELSE
764  mo_notconv => mo_coeff_b
765  mo_notconv_fm => mo_coeff
766  c_out => matrix_h
767  END IF
768 
769  ! allocate matrix_pz using as template matrix_z
770  CALL dbcsr_init_p(matrix_pz)
771  CALL dbcsr_create(matrix_pz, template=matrix_z, &
772  name="matrix_pz", &
773  matrix_type=dbcsr_type_no_symmetry, &
774  nze=0, data_type=dbcsr_type_real_default)
775 
776  IF (do_apply_preconditioner) THEN
777  IF (preconditioner%in_use /= 0) THEN
778  CALL apply_preconditioner(preconditioner, matrix_z, matrix_pz)
779  ELSE
780  CALL dbcsr_copy(matrix_pz, matrix_z)
781  END IF
782  ELSE
783  CALL dbcsr_copy(matrix_pz, matrix_z)
784  END IF
785 
786  !allocate NMOxNMO full matrices
787  nmat = nmo_not_converged
788  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmat, ncol_global=nmat, &
789  context=mo_coeff%matrix_struct%context, &
790  para_env=mo_coeff%matrix_struct%para_env)
791  CALL cp_fm_create(matrix_mm_fm, fm_struct_tmp, name="m_tmp_mxm")
792  CALL cp_fm_create(matrix_mmt_fm, fm_struct_tmp, name="mt_tmp_mxm")
793  CALL cp_fm_struct_release(fm_struct_tmp)
794 
795  !allocate 2NMOx2NMO full matrices
796  nmat2 = 2*nmo_not_converged
797  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmat2, ncol_global=nmat2, &
798  context=mo_coeff%matrix_struct%context, &
799  para_env=mo_coeff%matrix_struct%para_env)
800 
801  CALL cp_fm_create(s_block, fm_struct_tmp, name="sb")
802  CALL cp_fm_create(h_block, fm_struct_tmp, name="hb")
803  CALL cp_fm_create(v_block, fm_struct_tmp, name="vb")
804  CALL cp_fm_create(w_block, fm_struct_tmp, name="wb")
805  ALLOCATE (evals(nmat2))
806  CALL cp_fm_struct_release(fm_struct_tmp)
807 
808  ! compute CSC
809  CALL cp_fm_set_all(s_block, 0.0_dp, 1.0_dp)
810  ! compute CHC
811  CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
812  CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1, 1)
813 
814  ! compute the bottom left ZSC (top right is transpose)
815  CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_sc, 0.0_dp, matrix_mm, last_column=nmat)
816  ! set the bottom left part of S[C,Z] block matrix ZSC
817  !copy sparse to full
818  CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
819  CALL cp_fm_to_fm_submat(matrix_mm_fm, s_block, nmat, nmat, 1, 1, 1 + nmat, 1)
820  CALL cp_fm_transpose(matrix_mm_fm, matrix_mmt_fm)
821  CALL cp_fm_to_fm_submat(matrix_mmt_fm, s_block, nmat, nmat, 1, 1, 1, 1 + nmat)
822 
823  ! compute the bottom left ZHC (top right is transpose)
824  CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_hc, 0.0_dp, matrix_mm, last_column=nmat)
825  ! set the bottom left part of S[C,Z] block matrix ZHC
826  CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
827  CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1 + nmat, 1)
828  CALL cp_fm_transpose(matrix_mm_fm, matrix_mmt_fm)
829  CALL cp_fm_to_fm_submat(matrix_mmt_fm, h_block, nmat, nmat, 1, 1, 1, 1 + nmat)
830 
831  CALL cp_fm_release(matrix_mmt_fm)
832 
833  ! (reuse matrix_sc and matrix_hc to computr HZ and SZ)
834  CALL dbcsr_get_info(matrix_pz, nfullrows_total=n, nfullcols_total=k)
835  CALL dbcsr_multiply('n', 'n', 1.0_dp, c_out, matrix_pz, 0.0_dp, matrix_hc, last_column=k)
836  CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_s, matrix_pz, 0.0_dp, matrix_sc, last_column=k)
837 
838  ! compute the bottom right ZSZ
839  CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_sc, 0.0_dp, matrix_mm, last_column=k)
840  ! set the bottom right part of S[C,Z] block matrix ZSZ
841  CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
842  CALL cp_fm_to_fm_submat(matrix_mm_fm, s_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
843 
844  ! compute the bottom right ZHZ
845  CALL dbcsr_multiply('t', 'n', 1.0_dp, matrix_pz, matrix_hc, 0.0_dp, matrix_mm, last_column=k)
846  ! set the bottom right part of H[C,Z] block matrix ZHZ
847  CALL copy_dbcsr_to_fm(matrix_mm, matrix_mm_fm)
848  CALL cp_fm_to_fm_submat(matrix_mm_fm, h_block, nmat, nmat, 1, 1, 1 + nmat, 1 + nmat)
849 
850  CALL dbcsr_release_p(matrix_mm)
851  CALL dbcsr_release_p(matrix_sc)
852  CALL dbcsr_release_p(matrix_hc)
853 
854  CALL reduce_extended_space(s_block, h_block, v_block, w_block, evals, nmat2)
855 
856  ! allocate two (nao x nmat) full matrix
857  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmat, &
858  context=mo_coeff%matrix_struct%context, &
859  para_env=mo_coeff%matrix_struct%para_env)
860  CALL cp_fm_create(matrix_nm_fm, fm_struct_tmp, name="m_nxm")
861  CALL cp_fm_create(matrix_z_fm, fm_struct_tmp, name="m_nxm")
862  CALL cp_fm_struct_release(fm_struct_tmp)
863 
864  CALL copy_dbcsr_to_fm(matrix_pz, matrix_z_fm)
865  ! extract egenvectors
866  CALL cp_fm_to_fm_submat(v_block, matrix_mm_fm, nmat, nmat, 1, 1, 1, 1)
867  CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, mo_notconv_fm, matrix_mm_fm, 0.0_dp, matrix_nm_fm)
868  CALL cp_fm_to_fm_submat(v_block, matrix_mm_fm, nmat, nmat, 1 + nmat, 1, 1, 1)
869  CALL parallel_gemm('N', 'N', nao, nmat, nmat, 1.0_dp, matrix_z_fm, matrix_mm_fm, 1.0_dp, matrix_nm_fm)
870 
871  CALL dbcsr_release_p(matrix_z)
872  CALL dbcsr_release_p(matrix_pz)
873  CALL cp_fm_release(matrix_z_fm)
874  CALL cp_fm_release(s_block)
875  CALL cp_fm_release(h_block)
876  CALL cp_fm_release(w_block)
877  CALL cp_fm_release(v_block)
878  CALL cp_fm_release(matrix_mm_fm)
879 
880  ! in case some vector are already converged only a subset of vectors are copied in the MOS
881  IF (nmo_converged > 0) THEN
882  jj = 1
883  DO j = 1, nset_not_conv
884  i_first = inotconv_set(j, 1)
885  i_last = inotconv_set(j, 2)
886  n = i_last - i_first + 1
887  CALL cp_fm_to_fm_submat(matrix_nm_fm, mo_coeff, nao, n, 1, jj, 1, i_first)
888  eigenvalues(i_first:i_last) = evals(jj:jj + n - 1)
889  jj = jj + n
890  END DO
891  DEALLOCATE (iconv_set)
892  DEALLOCATE (inotconv_set)
893 
894  CALL dbcsr_release_p(mo_notconv)
895  CALL dbcsr_release_p(c_out)
896  CALL cp_fm_release(mo_notconv_fm)
897  DEALLOCATE (mo_notconv_fm)
898  ELSE
899  CALL cp_fm_to_fm(matrix_nm_fm, mo_coeff)
900  eigenvalues(1:nmo) = evals(1:nmo)
901  END IF
902  DEALLOCATE (evals)
903 
904  CALL cp_fm_release(matrix_nm_fm)
905  CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b) !fm->dbcsr
906 
907  t2 = m_walltime()
908  IF (output_unit > 0) THEN
909  WRITE (output_unit, '(T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
910  iter, nmo_converged, max_norm, min_norm, t2 - t1
911  END IF
912  t1 = m_walltime()
913 
914  END DO ! iter
915 
916  DEALLOCATE (ritz_coeff)
917  DEALLOCATE (iconv)
918  DEALLOCATE (inotconv)
919  DEALLOCATE (vnorm)
920 
921  CALL timestop(handle)
922 
923  END SUBROUTINE generate_extended_space_sparse
924 
925 ! **************************************************************************************************
926 
927 ! **************************************************************************************************
928 !> \brief ...
929 !> \param s_block ...
930 !> \param h_block ...
931 !> \param v_block ...
932 !> \param w_block ...
933 !> \param evals ...
934 !> \param ndim ...
935 ! **************************************************************************************************
936  SUBROUTINE reduce_extended_space(s_block, h_block, v_block, w_block, evals, ndim)
937 
938  TYPE(cp_fm_type), INTENT(IN) :: s_block, h_block, v_block, w_block
939  REAL(dp), DIMENSION(:) :: evals
940  INTEGER :: ndim
941 
942  CHARACTER(len=*), PARAMETER :: routinen = 'reduce_extended_space'
943 
944  INTEGER :: handle, info
945 
946  CALL timeset(routinen, handle)
947 
948  CALL cp_fm_to_fm(s_block, w_block)
949  CALL cp_fm_cholesky_decompose(s_block, info_out=info)
950  IF (info == 0) THEN
951  CALL cp_fm_triangular_invert(s_block)
952  CALL cp_fm_cholesky_restore(h_block, ndim, s_block, w_block, "MULTIPLY", pos="RIGHT")
953  CALL cp_fm_cholesky_restore(w_block, ndim, s_block, h_block, "MULTIPLY", pos="LEFT", transa="T")
954  CALL choose_eigv_solver(h_block, w_block, evals)
955  CALL cp_fm_cholesky_restore(w_block, ndim, s_block, v_block, "MULTIPLY")
956  ELSE
957 ! S^(-1/2)
958  CALL cp_fm_power(w_block, s_block, -0.5_dp, 1.0e-5_dp, info)
959  CALL cp_fm_to_fm(w_block, s_block)
960  CALL parallel_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, h_block, s_block, 0.0_dp, w_block)
961  CALL parallel_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, s_block, w_block, 0.0_dp, h_block)
962  CALL choose_eigv_solver(h_block, w_block, evals)
963  CALL parallel_gemm('N', 'N', ndim, ndim, ndim, 1.0_dp, s_block, w_block, 0.0_dp, v_block)
964  END IF
965 
966  CALL timestop(handle)
967 
968  END SUBROUTINE reduce_extended_space
969 
970 END MODULE qs_scf_block_davidson
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_dbcsr_m_by_n_from_row_template
subroutine, public cp_dbcsr_m_by_n_from_row_template(matrix, template, n, sym, data_type)
Utility function to create dbcsr matrix, m x n matrix (n arbitrary) with the same processor grid and ...
Definition: cp_dbcsr_operations.F:1156
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition: cp_dbcsr_operations.F:118
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_transpose
subroutine, public cp_fm_transpose(matrix, matrixt)
transposes a matrix matrixt = matrix ^ T
Definition: cp_fm_basic_linalg.F:1644
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_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition: cp_fm_basic_linalg.F:101
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_symm
subroutine, public cp_fm_symm(side, uplo, m, n, alpha, matrix_a, matrix_b, beta, matrix_c)
computes matrix_c = beta * matrix_c + alpha * matrix_a * matrix_b computes matrix_c = beta * matrix_c...
Definition: cp_fm_basic_linalg.F:510
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_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::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition: cp_fm_diag.F:896
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_vectorsnorm
subroutine, public cp_fm_vectorsnorm(matrix, norm_array)
find the inorm of each column norm_{j}= sqrt( \sum_{i} A_{ij}*A_{ij} )
Definition: cp_fm_types.F:1207
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition: cp_fm_types.F:1473
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
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
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition: machine.F:123
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
preconditioner_types
types of preconditioners
Definition: preconditioner_types.F:14
preconditioner
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition: preconditioner.F:15
qs_block_davidson_types
module that contains the algorithms to perform an itrative diagonalization by the block-Davidson appr...
Definition: qs_block_davidson_types.F:18
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_scf_block_davidson
module that contains the algorithms to perform an itrative diagonalization by the block-Davidson appr...
Definition: qs_scf_block_davidson.F:18
qs_scf_block_davidson::generate_extended_space_sparse
subroutine, public generate_extended_space_sparse(bdav_env, mo_set, matrix_h, matrix_s, output_unit, preconditioner)
...
Definition: qs_scf_block_davidson.F:475
qs_scf_block_davidson::generate_extended_space
subroutine, public generate_extended_space(bdav_env, mo_set, matrix_h, matrix_s, output_unit, preconditioner)
...
Definition: qs_scf_block_davidson.F:79