(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
18MODULE 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
66CONTAINS
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
970END MODULE qs_scf_block_davidson
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
preconditioner::apply_preconditioner
Definition preconditioner.F:68
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:1882
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:1710
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:1759
cp_fm_basic_linalg::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa)
...
Definition cp_fm_basic_linalg.F:3162
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:167
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:2193
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:576
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_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
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:116
preconditioner_types::preconditioner_type
Definition preconditioner_types.F:42
qs_block_davidson_types::davidson_type
Definition qs_block_davidson_types.F:38
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46