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