Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief module that contains the algorithms to perform an iterative
10 : !> diagonalization by the block-Lanczos approach
11 : !> \par History
12 : !> 05.2009 created [MI]
13 : !> \author fawzi
14 : ! **************************************************************************************************
15 : MODULE qs_scf_lanczos
16 :
17 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
18 : cp_fm_qr_factorization,&
19 : cp_fm_scale_and_add,&
20 : cp_fm_transpose,&
21 : cp_fm_triangular_multiply
22 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
23 : USE cp_fm_diag, ONLY: choose_eigv_solver
24 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
25 : cp_fm_struct_release,&
26 : cp_fm_struct_type
27 : USE cp_fm_types, ONLY: cp_fm_create,&
28 : cp_fm_get_submatrix,&
29 : cp_fm_release,&
30 : cp_fm_set_all,&
31 : cp_fm_set_submatrix,&
32 : cp_fm_to_fm,&
33 : cp_fm_type,&
34 : cp_fm_vectorsnorm
35 : USE cp_log_handling, ONLY: cp_to_string
36 : USE kinds, ONLY: dp
37 : USE parallel_gemm_api, ONLY: parallel_gemm
38 : USE qs_mo_types, ONLY: get_mo_set,&
39 : mo_set_type
40 : USE qs_scf_types, ONLY: krylov_space_type
41 : USE scf_control_types, ONLY: scf_control_type
42 : #include "./base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 : PRIVATE
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_lanczos'
47 :
48 : PUBLIC :: krylov_space_allocate, lanczos_refinement, lanczos_refinement_2v
49 :
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 :
54 : ! **************************************************************************************************
55 : !> \brief allocates matrices and vectors used in the construction of
56 : !> the krylov space and for the lanczos refinement
57 : !> \param krylov_space ...
58 : !> \param scf_control ...
59 : !> \param mos ...
60 : !> \param
61 : !> \par History
62 : !> 05.2009 created [MI]
63 : ! **************************************************************************************************
64 :
65 10 : SUBROUTINE krylov_space_allocate(krylov_space, scf_control, mos)
66 :
67 : TYPE(krylov_space_type), INTENT(INOUT) :: krylov_space
68 : TYPE(scf_control_type), POINTER :: scf_control
69 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
70 :
71 : CHARACTER(LEN=*), PARAMETER :: routineN = 'krylov_space_allocate'
72 :
73 : INTEGER :: handle, ik, ispin, max_nmo, nao, nblock, &
74 : ndim, nk, nmo, nspin
75 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
76 : TYPE(cp_fm_type), POINTER :: mo_coeff
77 :
78 10 : CALL timeset(routineN, handle)
79 :
80 10 : IF (.NOT. ASSOCIATED(krylov_space%mo_conv)) THEN
81 6 : NULLIFY (fm_struct_tmp, mo_coeff)
82 :
83 6 : krylov_space%nkrylov = scf_control%diagonalization%nkrylov
84 6 : krylov_space%nblock = scf_control%diagonalization%nblock_krylov
85 :
86 6 : nk = krylov_space%nkrylov
87 6 : nblock = krylov_space%nblock
88 6 : nspin = SIZE(mos, 1)
89 :
90 24 : ALLOCATE (krylov_space%mo_conv(nspin))
91 18 : ALLOCATE (krylov_space%mo_refine(nspin))
92 18 : ALLOCATE (krylov_space%chc_mat(nspin))
93 18 : ALLOCATE (krylov_space%c_vec(nspin))
94 6 : max_nmo = 0
95 12 : DO ispin = 1, nspin
96 6 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
97 6 : CALL cp_fm_create(krylov_space%mo_conv(ispin), mo_coeff%matrix_struct)
98 6 : CALL cp_fm_create(krylov_space%mo_refine(ispin), mo_coeff%matrix_struct)
99 6 : NULLIFY (fm_struct_tmp)
100 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=nmo, &
101 : para_env=mo_coeff%matrix_struct%para_env, &
102 6 : context=mo_coeff%matrix_struct%context)
103 6 : CALL cp_fm_create(krylov_space%chc_mat(ispin), fm_struct_tmp, "chc")
104 6 : CALL cp_fm_create(krylov_space%c_vec(ispin), fm_struct_tmp, "vec")
105 6 : CALL cp_fm_struct_release(fm_struct_tmp)
106 18 : max_nmo = MAX(max_nmo, nmo)
107 : END DO
108 :
109 : !the use of max_nmo might not be ok, in this case allocate nspin matrices
110 18 : ALLOCATE (krylov_space%c_eval(max_nmo))
111 :
112 66 : ALLOCATE (krylov_space%v_mat(nk))
113 :
114 6 : NULLIFY (fm_struct_tmp)
115 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nblock, &
116 : para_env=mo_coeff%matrix_struct%para_env, &
117 6 : context=mo_coeff%matrix_struct%context)
118 54 : DO ik = 1, nk
119 : CALL cp_fm_create(krylov_space%v_mat(ik), matrix_struct=fm_struct_tmp, &
120 54 : name="v_mat_"//TRIM(ADJUSTL(cp_to_string(ik))))
121 : END DO
122 6 : ALLOCATE (krylov_space%tmp_mat)
123 : CALL cp_fm_create(krylov_space%tmp_mat, matrix_struct=fm_struct_tmp, &
124 6 : name="tmp_mat")
125 6 : CALL cp_fm_struct_release(fm_struct_tmp)
126 :
127 : ! NOTE: the following matrices are small and could be defined
128 : ! as standard array rather than distributed fm
129 6 : NULLIFY (fm_struct_tmp)
130 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nblock, ncol_global=nblock, &
131 : para_env=mo_coeff%matrix_struct%para_env, &
132 6 : context=mo_coeff%matrix_struct%context)
133 6 : ALLOCATE (krylov_space%block1_mat)
134 : CALL cp_fm_create(krylov_space%block1_mat, matrix_struct=fm_struct_tmp, &
135 6 : name="a_mat_"//TRIM(ADJUSTL(cp_to_string(ik))))
136 6 : ALLOCATE (krylov_space%block2_mat)
137 : CALL cp_fm_create(krylov_space%block2_mat, matrix_struct=fm_struct_tmp, &
138 6 : name="b_mat_"//TRIM(ADJUSTL(cp_to_string(ik))))
139 6 : ALLOCATE (krylov_space%block3_mat)
140 : CALL cp_fm_create(krylov_space%block3_mat, matrix_struct=fm_struct_tmp, &
141 6 : name="b2_mat_"//TRIM(ADJUSTL(cp_to_string(ik))))
142 6 : CALL cp_fm_struct_release(fm_struct_tmp)
143 :
144 6 : ndim = nblock*nk
145 6 : NULLIFY (fm_struct_tmp)
146 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ndim, ncol_global=ndim, &
147 : para_env=mo_coeff%matrix_struct%para_env, &
148 6 : context=mo_coeff%matrix_struct%context)
149 6 : ALLOCATE (krylov_space%block4_mat)
150 : CALL cp_fm_create(krylov_space%block4_mat, matrix_struct=fm_struct_tmp, &
151 6 : name="t_mat")
152 6 : ALLOCATE (krylov_space%block5_mat)
153 : CALL cp_fm_create(krylov_space%block5_mat, matrix_struct=fm_struct_tmp, &
154 6 : name="t_vec")
155 6 : CALL cp_fm_struct_release(fm_struct_tmp)
156 18 : ALLOCATE (krylov_space%t_eval(ndim))
157 :
158 : ELSE
159 : !Nothing should be done
160 : END IF
161 :
162 10 : CALL timestop(handle)
163 :
164 10 : END SUBROUTINE krylov_space_allocate
165 :
166 : ! **************************************************************************************************
167 : !> \brief lanczos refinement by blocks of non-converged MOs
168 : !> \param krylov_space ...
169 : !> \param ks ...
170 : !> \param c0 ...
171 : !> \param c1 ...
172 : !> \param eval ...
173 : !> \param nao ...
174 : !> \param eps_iter ...
175 : !> \param ispin ...
176 : !> \param check_moconv_only ...
177 : !> \param
178 : !> \par History
179 : !> 05.2009 created [MI]
180 : ! **************************************************************************************************
181 :
182 0 : SUBROUTINE lanczos_refinement(krylov_space, ks, c0, c1, eval, nao, &
183 : eps_iter, ispin, check_moconv_only)
184 :
185 : TYPE(krylov_space_type), POINTER :: krylov_space
186 : TYPE(cp_fm_type), INTENT(IN) :: ks, c0, c1
187 : REAL(dp), DIMENSION(:), POINTER :: eval
188 : INTEGER, INTENT(IN) :: nao
189 : REAL(dp), INTENT(IN) :: eps_iter
190 : INTEGER, INTENT(IN) :: ispin
191 : LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
192 :
193 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lanczos_refinement'
194 : REAL(KIND=dp), PARAMETER :: rmone = -1.0_dp, rone = 1.0_dp, &
195 : rzero = 0.0_dp
196 :
197 : INTEGER :: hand1, hand2, hand3, hand4, hand5, handle, ib, ik, imo, imo_low, imo_up, it, jt, &
198 : nblock, ndim, nmo, nmo_converged, nmo_nc, nmob, num_blocks
199 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: itaken
200 : LOGICAL :: my_check_moconv_only
201 : REAL(dp) :: max_norm, min_norm, vmax
202 0 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: q_mat, tblock, tvblock
203 : REAL(dp), DIMENSION(:), POINTER :: c_res, t_eval
204 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
205 : TYPE(cp_fm_type) :: c2_tmp, c3_tmp, c_tmp, hc
206 0 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: v_mat
207 : TYPE(cp_fm_type), POINTER :: a_mat, b2_mat, b_mat, chc, evec, t_mat, &
208 : t_vec
209 :
210 0 : CALL timeset(routineN, handle)
211 :
212 0 : NULLIFY (fm_struct_tmp)
213 0 : NULLIFY (chc, evec)
214 : NULLIFY (c_res, t_eval)
215 0 : NULLIFY (t_mat, t_vec)
216 0 : NULLIFY (a_mat, b_mat, b2_mat, v_mat)
217 :
218 0 : nmo = SIZE(eval, 1)
219 0 : my_check_moconv_only = .FALSE.
220 0 : IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
221 :
222 0 : chc => krylov_space%chc_mat(ispin)
223 0 : evec => krylov_space%c_vec(ispin)
224 0 : c_res => krylov_space%c_eval
225 0 : t_eval => krylov_space%t_eval
226 :
227 : NULLIFY (fm_struct_tmp)
228 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo, &
229 : para_env=c0%matrix_struct%para_env, &
230 0 : context=c0%matrix_struct%context)
231 : CALL cp_fm_create(c_tmp, matrix_struct=fm_struct_tmp, &
232 0 : name="c_tmp")
233 : CALL cp_fm_create(hc, matrix_struct=fm_struct_tmp, &
234 0 : name="hc")
235 0 : CALL cp_fm_struct_release(fm_struct_tmp)
236 :
237 : !Compute (C^t)HC
238 0 : CALL parallel_gemm('N', 'N', nao, nmo, nao, rone, ks, c0, rzero, hc)
239 0 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, hc, rzero, chc)
240 :
241 : !Diagonalize (C^t)HC
242 0 : CALL timeset(routineN//"diag_chc", hand1)
243 0 : CALL choose_eigv_solver(chc, evec, eval)
244 0 : CALL timestop(hand1)
245 :
246 : !Rotate the C vectors
247 0 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, c1)
248 :
249 : !Check for converged states
250 0 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, hc, evec, rzero, c_tmp)
251 0 : CALL cp_fm_to_fm(c1, c0, nmo, 1, 1)
252 0 : CALL cp_fm_column_scale(c1, eval)
253 0 : CALL cp_fm_scale_and_add(1.0_dp, c_tmp, rmone, c1)
254 0 : CALL cp_fm_vectorsnorm(c_tmp, c_res)
255 :
256 0 : nmo_converged = 0
257 0 : nmo_nc = 0
258 0 : max_norm = 0.0_dp
259 0 : min_norm = 1.e10_dp
260 0 : CALL cp_fm_set_all(c1, rzero)
261 0 : DO imo = 1, nmo
262 0 : max_norm = MAX(max_norm, c_res(imo))
263 0 : min_norm = MIN(min_norm, c_res(imo))
264 : END DO
265 0 : DO imo = 1, nmo
266 0 : IF (c_res(imo) <= eps_iter) THEN
267 0 : nmo_converged = nmo_converged + 1
268 : ELSE
269 0 : nmo_nc = nmo - nmo_converged
270 0 : EXIT
271 : END IF
272 : END DO
273 :
274 0 : nblock = krylov_space%nblock
275 0 : num_blocks = nmo_nc/nblock
276 :
277 0 : krylov_space%nmo_nc = nmo_nc
278 0 : krylov_space%nmo_conv = nmo_converged
279 0 : krylov_space%max_res_norm = max_norm
280 0 : krylov_space%min_res_norm = min_norm
281 :
282 0 : IF (my_check_moconv_only) THEN
283 0 : CALL cp_fm_release(c_tmp)
284 0 : CALL cp_fm_release(hc)
285 0 : CALL timestop(handle)
286 0 : RETURN
287 0 : ELSE IF (krylov_space%nmo_nc > 0) THEN
288 :
289 0 : CALL cp_fm_to_fm(c0, c1, nmo_nc, nmo_converged + 1, 1)
290 :
291 0 : nblock = krylov_space%nblock
292 0 : IF (MODULO(nmo_nc, nblock) > 0.0_dp) THEN
293 0 : num_blocks = nmo_nc/nblock + 1
294 : ELSE
295 0 : num_blocks = nmo_nc/nblock
296 : END IF
297 :
298 0 : DO ib = 1, num_blocks
299 :
300 0 : imo_low = (ib - 1)*nblock + 1
301 0 : imo_up = MIN(ib*nblock, nmo_nc)
302 0 : nmob = imo_up - imo_low + 1
303 0 : ndim = krylov_space%nkrylov*nmob
304 :
305 0 : NULLIFY (fm_struct_tmp)
306 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=ndim, &
307 : para_env=c0%matrix_struct%para_env, &
308 0 : context=c0%matrix_struct%context)
309 : CALL cp_fm_create(c2_tmp, matrix_struct=fm_struct_tmp, &
310 0 : name="c2_tmp")
311 0 : CALL cp_fm_struct_release(fm_struct_tmp)
312 0 : NULLIFY (fm_struct_tmp)
313 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmob, ncol_global=ndim, &
314 : para_env=c0%matrix_struct%para_env, &
315 0 : context=c0%matrix_struct%context)
316 : CALL cp_fm_create(c3_tmp, matrix_struct=fm_struct_tmp, &
317 0 : name="c3_tmp")
318 0 : CALL cp_fm_struct_release(fm_struct_tmp)
319 :
320 : ! Create local matrix of right size
321 0 : IF (nmob /= nblock) THEN
322 0 : NULLIFY (a_mat, b_mat, b2_mat, t_mat, t_vec, v_mat)
323 0 : ALLOCATE (a_mat, b_mat, b2_mat, t_mat, t_vec)
324 0 : NULLIFY (fm_struct_tmp)
325 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmob, ncol_global=nmob, &
326 : para_env=chc%matrix_struct%para_env, &
327 0 : context=chc%matrix_struct%context)
328 : CALL cp_fm_create(a_mat, matrix_struct=fm_struct_tmp, &
329 0 : name="a_mat")
330 : CALL cp_fm_create(b_mat, matrix_struct=fm_struct_tmp, &
331 0 : name="b_mat")
332 : CALL cp_fm_create(b2_mat, matrix_struct=fm_struct_tmp, &
333 0 : name="b2_mat")
334 0 : CALL cp_fm_struct_release(fm_struct_tmp)
335 0 : NULLIFY (fm_struct_tmp)
336 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ndim, ncol_global=ndim, &
337 : para_env=chc%matrix_struct%para_env, &
338 0 : context=chc%matrix_struct%context)
339 : CALL cp_fm_create(t_mat, matrix_struct=fm_struct_tmp, &
340 0 : name="t_mat")
341 : CALL cp_fm_create(t_vec, matrix_struct=fm_struct_tmp, &
342 0 : name="t_vec")
343 0 : CALL cp_fm_struct_release(fm_struct_tmp)
344 0 : NULLIFY (fm_struct_tmp)
345 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmob, &
346 : para_env=c0%matrix_struct%para_env, &
347 0 : context=c0%matrix_struct%context)
348 0 : ALLOCATE (v_mat(krylov_space%nkrylov))
349 0 : DO ik = 1, krylov_space%nkrylov
350 : CALL cp_fm_create(v_mat(ik), matrix_struct=fm_struct_tmp, &
351 0 : name="v_mat")
352 : END DO
353 0 : CALL cp_fm_struct_release(fm_struct_tmp)
354 : ELSE
355 0 : a_mat => krylov_space%block1_mat
356 0 : b_mat => krylov_space%block2_mat
357 0 : b2_mat => krylov_space%block3_mat
358 0 : t_mat => krylov_space%block4_mat
359 0 : t_vec => krylov_space%block5_mat
360 0 : v_mat => krylov_space%v_mat
361 : END IF
362 :
363 0 : ALLOCATE (tblock(nmob, nmob))
364 0 : ALLOCATE (tvblock(nmob, ndim))
365 :
366 0 : CALL timeset(routineN//"_kry_loop", hand2)
367 0 : CALL cp_fm_set_all(b_mat, rzero)
368 0 : CALL cp_fm_set_all(t_mat, rzero)
369 0 : CALL cp_fm_to_fm(c1, v_mat(1), nmob, imo_low, 1)
370 :
371 : !Compute A =(V^t)HV
372 0 : CALL parallel_gemm('N', 'N', nao, nmob, nao, rone, ks, v_mat(1), rzero, hc)
373 : CALL parallel_gemm('T', 'N', nmob, nmob, nao, rone, v_mat(1), hc, &
374 0 : rzero, a_mat)
375 :
376 : !Compute the residual matrix R for next
377 : !factorisation
378 : CALL parallel_gemm('N', 'N', nao, nmob, nmob, rone, v_mat(1), a_mat, &
379 0 : rzero, c_tmp)
380 0 : CALL cp_fm_scale_and_add(rmone, c_tmp, rone, hc)
381 :
382 : ! Build the block tridiagonal matrix
383 0 : CALL cp_fm_get_submatrix(a_mat, tblock, 1, 1, nmob, nmob)
384 0 : CALL cp_fm_set_submatrix(t_mat, tblock, 1, 1, nmob, nmob)
385 :
386 0 : DO ik = 2, krylov_space%nkrylov
387 :
388 : ! Call lapack for QR factorization
389 0 : CALL cp_fm_set_all(b_mat, rzero)
390 0 : CALL cp_fm_to_fm(c_tmp, v_mat(ik), nmob, 1, 1)
391 0 : CALL cp_fm_qr_factorization(c_tmp, b_mat, nao, nmob, 1, 1)
392 :
393 : CALL cp_fm_triangular_multiply(b_mat, v_mat(ik), side="R", invert_tr=.TRUE., &
394 0 : n_rows=nao, n_cols=nmob)
395 :
396 : !Compute A =(V^t)HV
397 0 : CALL parallel_gemm('N', 'N', nao, nmob, nao, rone, ks, v_mat(ik), rzero, hc)
398 0 : CALL parallel_gemm('T', 'N', nmob, nmob, nao, rone, v_mat(ik), hc, rzero, a_mat)
399 :
400 : !Compute the !residual matrix R !for next !factorisation
401 : CALL parallel_gemm('N', 'N', nao, nmob, nmob, rone, v_mat(ik), a_mat, &
402 0 : rzero, c_tmp)
403 0 : CALL cp_fm_scale_and_add(rmone, c_tmp, rone, hc)
404 0 : CALL cp_fm_to_fm(v_mat(ik - 1), hc, nmob, 1, 1)
405 : CALL cp_fm_triangular_multiply(b_mat, hc, side='R', transpose_tr=.TRUE., &
406 0 : n_rows=nao, n_cols=nmob, alpha=rmone)
407 0 : CALL cp_fm_scale_and_add(rone, c_tmp, rone, hc)
408 :
409 : ! Build the block tridiagonal matrix
410 0 : it = (ik - 2)*nmob + 1
411 0 : jt = (ik - 1)*nmob + 1
412 :
413 0 : CALL cp_fm_get_submatrix(a_mat, tblock, 1, 1, nmob, nmob)
414 0 : CALL cp_fm_set_submatrix(t_mat, tblock, jt, jt, nmob, nmob)
415 0 : CALL cp_fm_transpose(b_mat, a_mat)
416 0 : CALL cp_fm_get_submatrix(a_mat, tblock, 1, 1, nmob, nmob)
417 0 : CALL cp_fm_set_submatrix(t_mat, tblock, it, jt, nmob, nmob)
418 :
419 : END DO ! ik
420 0 : CALL timestop(hand2)
421 :
422 0 : DEALLOCATE (tblock)
423 :
424 0 : CALL timeset(routineN//"_diag_tri", hand3)
425 :
426 0 : CALL choose_eigv_solver(t_mat, t_vec, t_eval)
427 : ! Diagonalize the block-tridiagonal matrix
428 0 : CALL timestop(hand3)
429 :
430 0 : CALL timeset(routineN//"_build_cnew", hand4)
431 : ! !Compute the refined vectors
432 0 : CALL cp_fm_set_all(c2_tmp, rzero)
433 0 : DO ik = 1, krylov_space%nkrylov
434 0 : jt = (ik - 1)*nmob
435 : CALL parallel_gemm('N', 'N', nao, ndim, nmob, rone, v_mat(ik), t_vec, rone, c2_tmp, &
436 0 : b_first_row=(jt + 1))
437 : END DO
438 0 : DEALLOCATE (tvblock)
439 :
440 0 : CALL cp_fm_set_all(c3_tmp, rzero)
441 0 : CALL parallel_gemm('T', 'N', nmob, ndim, nao, rone, v_mat(1), c2_tmp, rzero, c3_tmp)
442 :
443 : !Try to avoid linear dependencies
444 0 : ALLOCATE (q_mat(nmob, ndim))
445 : !get max
446 0 : CALL cp_fm_get_submatrix(c3_tmp, q_mat, 1, 1, nmob, ndim)
447 :
448 0 : ALLOCATE (itaken(ndim))
449 0 : itaken = 0
450 0 : DO it = 1, nmob
451 0 : vmax = 0.0_dp
452 0 : ik = 0
453 : !select index ik
454 0 : DO jt = 1, ndim
455 0 : IF (itaken(jt) == 0 .AND. ABS(q_mat(it, jt)) > vmax) THEN
456 0 : vmax = ABS(q_mat(it, jt))
457 0 : ik = jt
458 : END IF
459 : END DO
460 0 : IF (ik == 0) CYCLE ! all candidates NaN or zero: keep the entry vector
461 0 : itaken(ik) = 1
462 :
463 0 : CALL cp_fm_to_fm(c2_tmp, v_mat(1), 1, ik, it)
464 : END DO
465 0 : DEALLOCATE (itaken)
466 0 : DEALLOCATE (q_mat)
467 :
468 : !Copy in the converged set to enlarge the converged subspace
469 0 : CALL cp_fm_to_fm(v_mat(1), c0, nmob, 1, (nmo_converged + imo_low))
470 0 : CALL timestop(hand4)
471 :
472 0 : IF (nmob < nblock) THEN
473 0 : CALL cp_fm_release(a_mat)
474 0 : CALL cp_fm_release(b_mat)
475 0 : CALL cp_fm_release(b2_mat)
476 0 : CALL cp_fm_release(t_mat)
477 0 : CALL cp_fm_release(t_vec)
478 0 : DEALLOCATE (a_mat, b_mat, b2_mat, t_mat, t_vec)
479 0 : CALL cp_fm_release(v_mat)
480 : END IF
481 0 : CALL cp_fm_release(c2_tmp)
482 0 : CALL cp_fm_release(c3_tmp)
483 : END DO ! ib
484 :
485 0 : CALL timeset(routineN//"_ortho", hand5)
486 0 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, c0, rzero, chc)
487 :
488 0 : CALL cp_fm_cholesky_decompose(chc, nmo)
489 0 : CALL cp_fm_triangular_multiply(chc, c0, 'R', invert_tr=.TRUE.)
490 0 : CALL timestop(hand5)
491 :
492 0 : CALL cp_fm_release(c_tmp)
493 0 : CALL cp_fm_release(hc)
494 : ELSE
495 0 : CALL cp_fm_release(c_tmp)
496 0 : CALL cp_fm_release(hc)
497 0 : CALL timestop(handle)
498 0 : RETURN
499 : END IF
500 :
501 0 : CALL timestop(handle)
502 :
503 0 : END SUBROUTINE lanczos_refinement
504 :
505 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
506 :
507 : ! **************************************************************************************************
508 : !> \brief ...
509 : !> \param krylov_space ...
510 : !> \param ks ...
511 : !> \param c0 ...
512 : !> \param c1 ...
513 : !> \param eval ...
514 : !> \param nao ...
515 : !> \param eps_iter ...
516 : !> \param ispin ...
517 : !> \param check_moconv_only ...
518 : ! **************************************************************************************************
519 744 : SUBROUTINE lanczos_refinement_2v(krylov_space, ks, c0, c1, eval, nao, &
520 : eps_iter, ispin, check_moconv_only)
521 :
522 : TYPE(krylov_space_type), POINTER :: krylov_space
523 : TYPE(cp_fm_type), INTENT(IN) :: ks, c0, c1
524 : REAL(dp), DIMENSION(:), POINTER :: eval
525 : INTEGER, INTENT(IN) :: nao
526 : REAL(dp), INTENT(IN) :: eps_iter
527 : INTEGER, INTENT(IN) :: ispin
528 : LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
529 :
530 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lanczos_refinement_2v'
531 : REAL(KIND=dp), PARAMETER :: rmone = -1.0_dp, rone = 1.0_dp, &
532 : rzero = 0.0_dp
533 :
534 : INTEGER :: hand1, hand2, hand3, hand4, hand5, hand6, handle, i, ia, ib, ik, imo, imo_low, &
535 : imo_up, info, it, j, jt, liwork, lwork, nblock, ndim, nmo, nmo_converged, nmo_nc, nmob, &
536 : num_blocks
537 744 : INTEGER, ALLOCATABLE, DIMENSION(:) :: itaken
538 744 : INTEGER, DIMENSION(:), POINTER :: iwork
539 : LOGICAL :: my_check_moconv_only
540 : REAL(dp) :: max_norm, min_norm, vmax
541 744 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a_block, b_block, q_mat, t_mat
542 744 : REAL(dp), DIMENSION(:), POINTER :: c_res, t_eval
543 744 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
544 744 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a_loc, b_loc
545 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
546 : TYPE(cp_fm_type) :: b_mat, c2_tmp, c_tmp, hc, v_tmp
547 744 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: v_mat
548 : TYPE(cp_fm_type), POINTER :: chc, evec
549 :
550 744 : CALL timeset(routineN, handle)
551 :
552 744 : NULLIFY (fm_struct_tmp)
553 744 : NULLIFY (chc, evec)
554 744 : NULLIFY (c_res, t_eval)
555 744 : NULLIFY (b_loc, a_loc)
556 :
557 744 : nmo = SIZE(eval, 1)
558 744 : my_check_moconv_only = .FALSE.
559 744 : IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
560 :
561 744 : chc => krylov_space%chc_mat(ispin)
562 744 : evec => krylov_space%c_vec(ispin)
563 744 : c_res => krylov_space%c_eval
564 744 : t_eval => krylov_space%t_eval
565 :
566 : NULLIFY (fm_struct_tmp)
567 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmo, &
568 : para_env=c0%matrix_struct%para_env, &
569 744 : context=c0%matrix_struct%context)
570 : CALL cp_fm_create(c_tmp, matrix_struct=fm_struct_tmp, &
571 744 : name="c_tmp")
572 : CALL cp_fm_create(hc, matrix_struct=fm_struct_tmp, &
573 744 : name="hc")
574 744 : CALL cp_fm_struct_release(fm_struct_tmp)
575 :
576 : !Compute (C^t)HC
577 744 : CALL parallel_gemm('N', 'N', nao, nmo, nao, rone, ks, c0, rzero, hc)
578 744 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, hc, rzero, chc)
579 :
580 : !Diagonalize (C^t)HC
581 744 : CALL timeset(routineN//"diag_chc", hand1)
582 744 : CALL choose_eigv_solver(chc, evec, eval)
583 :
584 744 : CALL timestop(hand1)
585 :
586 744 : CALL timeset(routineN//"check_conv", hand6)
587 : !Rotate the C vectors
588 744 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, c1)
589 :
590 : !Check for converged states
591 744 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, hc, evec, rzero, c_tmp)
592 744 : CALL cp_fm_to_fm(c1, c0, nmo, 1, 1)
593 744 : CALL cp_fm_column_scale(c1, eval)
594 744 : CALL cp_fm_scale_and_add(1.0_dp, c_tmp, rmone, c1)
595 744 : CALL cp_fm_vectorsnorm(c_tmp, c_res)
596 :
597 744 : nmo_converged = 0
598 744 : nmo_nc = 0
599 744 : max_norm = 0.0_dp
600 744 : min_norm = 1.e10_dp
601 744 : CALL cp_fm_set_all(c1, rzero)
602 24792 : DO imo = 1, nmo
603 24048 : max_norm = MAX(max_norm, c_res(imo))
604 24792 : min_norm = MIN(min_norm, c_res(imo))
605 : END DO
606 19970 : DO imo = 1, nmo
607 19970 : IF (c_res(imo) <= eps_iter) THEN
608 19226 : nmo_converged = nmo_converged + 1
609 : ELSE
610 734 : nmo_nc = nmo - nmo_converged
611 734 : EXIT
612 : END IF
613 : END DO
614 744 : CALL timestop(hand6)
615 :
616 744 : CALL cp_fm_release(c_tmp)
617 744 : CALL cp_fm_release(hc)
618 :
619 744 : krylov_space%nmo_nc = nmo_nc
620 744 : krylov_space%nmo_conv = nmo_converged
621 744 : krylov_space%max_res_norm = max_norm
622 744 : krylov_space%min_res_norm = min_norm
623 :
624 744 : IF (my_check_moconv_only) THEN
625 : ! Do nothing
626 744 : ELSE IF (krylov_space%nmo_nc > 0) THEN
627 :
628 734 : CALL cp_fm_to_fm(c0, c1, nmo_nc, nmo_converged + 1, 1)
629 :
630 734 : nblock = krylov_space%nblock
631 734 : IF (MODULO(nmo_nc, nblock) > 0.0_dp) THEN
632 100 : num_blocks = nmo_nc/nblock + 1
633 : ELSE
634 634 : num_blocks = nmo_nc/nblock
635 : END IF
636 :
637 4062 : DO ib = 1, num_blocks
638 :
639 3328 : imo_low = (ib - 1)*nblock + 1
640 3328 : imo_up = MIN(ib*nblock, nmo_nc)
641 3328 : nmob = imo_up - imo_low + 1
642 3328 : ndim = krylov_space%nkrylov*nmob
643 :
644 : ! Allocation
645 3328 : CALL timeset(routineN//"alloc", hand6)
646 13312 : ALLOCATE (a_block(nmob, nmob))
647 9984 : ALLOCATE (b_block(nmob, nmob))
648 13312 : ALLOCATE (t_mat(ndim, ndim))
649 :
650 3328 : NULLIFY (fm_struct_tmp)
651 : ! by forcing ncol_block=nmo, the needed part of the matrix is distributed on a subset of processes
652 : ! this is due to the use of two-dimensional grids of processes
653 : ! nrow_global is distributed over num_pe(1)
654 : ! a local_data array is anyway allocated for the processes non included
655 : ! this should have a minimum size
656 : ! with ncol_local=1.
657 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmob, &
658 : ncol_block=nmob, &
659 : para_env=c0%matrix_struct%para_env, &
660 : context=c0%matrix_struct%context, &
661 3328 : force_block=.TRUE.)
662 : CALL cp_fm_create(c_tmp, matrix_struct=fm_struct_tmp, &
663 3328 : name="c_tmp")
664 3328 : CALL cp_fm_set_all(c_tmp, rzero)
665 : CALL cp_fm_create(v_tmp, matrix_struct=fm_struct_tmp, &
666 3328 : name="v_tmp")
667 3328 : CALL cp_fm_set_all(v_tmp, rzero)
668 3328 : CALL cp_fm_struct_release(fm_struct_tmp)
669 3328 : NULLIFY (fm_struct_tmp)
670 36608 : ALLOCATE (v_mat(krylov_space%nkrylov))
671 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nmob, &
672 : ncol_block=nmob, &
673 : para_env=c0%matrix_struct%para_env, &
674 : context=c0%matrix_struct%context, &
675 3328 : force_block=.TRUE.)
676 29952 : DO ik = 1, krylov_space%nkrylov
677 : CALL cp_fm_create(v_mat(ik), matrix_struct=fm_struct_tmp, &
678 29952 : name="v_mat", set_zero=.TRUE.)
679 : END DO
680 : CALL cp_fm_create(hc, matrix_struct=fm_struct_tmp, &
681 3328 : name="hc", set_zero=.TRUE.)
682 3328 : CALL cp_fm_struct_release(fm_struct_tmp)
683 3328 : NULLIFY (fm_struct_tmp)
684 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=ndim, &
685 : ncol_block=ndim, &
686 : para_env=c0%matrix_struct%para_env, &
687 : context=c0%matrix_struct%context, &
688 3328 : force_block=.TRUE.)
689 : CALL cp_fm_create(c2_tmp, matrix_struct=fm_struct_tmp, &
690 3328 : name="c2_tmp")
691 3328 : CALL cp_fm_struct_release(fm_struct_tmp)
692 :
693 3328 : NULLIFY (fm_struct_tmp)
694 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmob, ncol_global=nmob, &
695 : para_env=c0%matrix_struct%para_env, &
696 3328 : context=c0%matrix_struct%context)
697 : CALL cp_fm_create(b_mat, matrix_struct=fm_struct_tmp, &
698 3328 : name="b_mat")
699 3328 : CALL cp_fm_struct_release(fm_struct_tmp)
700 3328 : CALL timestop(hand6)
701 : !End allocation
702 :
703 3328 : CALL cp_fm_set_all(b_mat, rzero)
704 3328 : CALL cp_fm_to_fm(c1, v_mat(1), nmob, imo_low, 1)
705 :
706 : ! Here starts the construction of krylov space
707 3328 : CALL timeset(routineN//"_kry_loop", hand2)
708 : !Compute A =(V^t)HV
709 3328 : CALL parallel_gemm('N', 'N', nao, nmob, nao, rone, ks, v_mat(1), rzero, hc)
710 :
711 3328 : a_block = 0.0_dp
712 3328 : a_loc => v_mat(1)%local_data
713 3328 : b_loc => hc%local_data
714 :
715 3328 : IF (SIZE(hc%local_data, 2) == nmob) THEN
716 : ! this is a work around to avoid problems due to the two dimensional grid of processes
717 : CALL dgemm('T', 'N', nmob, nmob, SIZE(hc%local_data, 1), 1.0_dp, a_loc(1, 1), &
718 3328 : SIZE(hc%local_data, 1), b_loc(1, 1), SIZE(hc%local_data, 1), 0.0_dp, a_block(1, 1), nmob)
719 : END IF
720 3328 : CALL hc%matrix_struct%para_env%sum(a_block)
721 :
722 : !Compute the residual matrix R for next
723 : !factorisation
724 206313 : c_tmp%local_data = 0.0_dp
725 3328 : IF (SIZE(c_tmp%local_data, 2) == nmob) THEN
726 3328 : b_loc => c_tmp%local_data
727 : CALL dgemm('N', 'N', SIZE(c_tmp%local_data, 1), nmob, nmob, 1.0_dp, a_loc(1, 1), &
728 : SIZE(c_tmp%local_data, 1), a_block(1, 1), nmob, 0.0_dp, &
729 3328 : b_loc(1, 1), SIZE(c_tmp%local_data, 1))
730 : END IF
731 3328 : CALL cp_fm_scale_and_add(rmone, c_tmp, rone, hc)
732 :
733 : ! Build the block tridiagonal matrix
734 3328 : t_mat = 0.0_dp
735 8150 : DO i = 1, nmob
736 18772 : t_mat(1:nmob, i) = a_block(1:nmob, i)
737 : END DO
738 :
739 26624 : DO ik = 2, krylov_space%nkrylov
740 : ! Call lapack for QR factorization
741 23296 : CALL cp_fm_set_all(b_mat, rzero)
742 23296 : CALL cp_fm_to_fm(c_tmp, v_mat(ik), nmob, 1, 1)
743 23296 : CALL cp_fm_qr_factorization(c_tmp, b_mat, nao, nmob, 1, 1)
744 :
745 : CALL cp_fm_triangular_multiply(b_mat, v_mat(ik), side="R", invert_tr=.TRUE., &
746 23296 : n_rows=nao, n_cols=nmob)
747 :
748 23296 : CALL cp_fm_get_submatrix(b_mat, b_block, 1, 1, nmob, nmob)
749 :
750 : !Compute A =(V^t)HV
751 23296 : CALL parallel_gemm('N', 'N', nao, nmob, nao, rone, ks, v_mat(ik), rzero, hc)
752 :
753 23296 : a_block = 0.0_dp
754 23296 : IF (SIZE(hc%local_data, 2) == nmob) THEN
755 23296 : a_loc => v_mat(ik)%local_data
756 23296 : b_loc => hc%local_data
757 : CALL dgemm('T', 'N', nmob, nmob, SIZE(hc%local_data, 1), 1.0_dp, a_loc(1, 1), &
758 23296 : SIZE(hc%local_data, 1), b_loc(1, 1), SIZE(hc%local_data, 1), 0.0_dp, a_block(1, 1), nmob)
759 : END IF
760 23296 : CALL hc%matrix_struct%para_env%sum(a_block)
761 :
762 : !Compute the residual matrix R for next
763 : !factorisation
764 1444191 : c_tmp%local_data = 0.0_dp
765 23296 : IF (SIZE(c_tmp%local_data, 2) == nmob) THEN
766 23296 : a_loc => v_mat(ik)%local_data
767 23296 : b_loc => c_tmp%local_data
768 : CALL dgemm('N', 'N', SIZE(c_tmp%local_data, 1), nmob, nmob, 1.0_dp, a_loc(1, 1), &
769 : SIZE(c_tmp%local_data, 1), a_block(1, 1), nmob, 0.0_dp, &
770 23296 : b_loc(1, 1), SIZE(c_tmp%local_data, 1))
771 : END IF
772 23296 : CALL cp_fm_scale_and_add(rmone, c_tmp, rone, hc)
773 :
774 23296 : IF (SIZE(c_tmp%local_data, 2) == nmob) THEN
775 23296 : a_loc => v_mat(ik - 1)%local_data
776 57050 : DO j = 1, nmob
777 111104 : DO i = 1, j
778 2012899 : DO ia = 1, SIZE(c_tmp%local_data, 1)
779 1979145 : b_loc(ia, i) = b_loc(ia, i) - a_loc(ia, j)*b_block(i, j)
780 : END DO
781 : END DO
782 : END DO
783 : END IF
784 :
785 : ! Build the block tridiagonal matrix
786 23296 : it = (ik - 2)*nmob
787 23296 : jt = (ik - 1)*nmob
788 60378 : DO j = 1, nmob
789 108108 : t_mat(jt + 1:jt + nmob, jt + j) = a_block(1:nmob, j)
790 131404 : DO i = 1, nmob
791 74354 : t_mat(it + i, jt + j) = b_block(j, i)
792 108108 : t_mat(jt + j, it + i) = b_block(j, i)
793 : END DO
794 : END DO
795 : END DO ! ik
796 3328 : CALL timestop(hand2)
797 :
798 3328 : CALL timeset(routineN//"_diag_tri", hand3)
799 3328 : lwork = 1 + 6*ndim + 2*ndim**2 + 5000
800 3328 : liwork = 5*ndim + 3
801 9984 : ALLOCATE (work(lwork))
802 9984 : ALLOCATE (iwork(liwork))
803 :
804 : ! Diagonalize the block-tridiagonal matrix
805 : CALL dsyevd('V', 'U', ndim, t_mat(1, 1), ndim, t_eval(1), &
806 3328 : work(1), lwork, iwork(1), liwork, info)
807 3328 : DEALLOCATE (work)
808 3328 : DEALLOCATE (iwork)
809 3328 : CALL timestop(hand3)
810 :
811 3328 : CALL timeset(routineN//"_build_cnew", hand4)
812 : ! !Compute the refined vectors
813 :
814 1627208 : c2_tmp%local_data = 0.0_dp
815 13312 : ALLOCATE (q_mat(nmob, ndim))
816 3328 : q_mat = 0.0_dp
817 3328 : b_loc => c2_tmp%local_data
818 29952 : DO ik = 1, krylov_space%nkrylov
819 26624 : CALL cp_fm_to_fm(v_mat(ik), v_tmp, nmob, 1, 1)
820 29952 : IF (SIZE(c2_tmp%local_data, 2) == ndim) THEN
821 : ! a_loc => v_mat(ik)%local_data
822 26624 : a_loc => v_tmp%local_data
823 26624 : it = (ik - 1)*nmob
824 : CALL dgemm('N', 'N', SIZE(c2_tmp%local_data, 1), ndim, nmob, 1.0_dp, a_loc(1, 1), &
825 : SIZE(c2_tmp%local_data, 1), t_mat(it + 1, 1), ndim, 1.0_dp, &
826 26624 : b_loc(1, 1), SIZE(c2_tmp%local_data, 1))
827 : END IF
828 : END DO !ik
829 :
830 : !Try to avoid linear dependencies
831 3328 : CALL cp_fm_to_fm(v_mat(1), v_tmp, nmob, 1, 1)
832 3328 : IF (SIZE(c2_tmp%local_data, 2) == ndim) THEN
833 : ! a_loc => v_mat(1)%local_data
834 3328 : a_loc => v_tmp%local_data
835 3328 : b_loc => c2_tmp%local_data
836 : CALL dgemm('T', 'N', nmob, ndim, SIZE(v_tmp%local_data, 1), 1.0_dp, a_loc(1, 1), &
837 : SIZE(v_tmp%local_data, 1), b_loc(1, 1), SIZE(v_tmp%local_data, 1), &
838 3328 : 0.0_dp, q_mat(1, 1), nmob)
839 : END IF
840 3328 : CALL hc%matrix_struct%para_env%sum(q_mat)
841 :
842 9984 : ALLOCATE (itaken(ndim))
843 3328 : itaken = 0
844 8150 : DO it = 1, nmob
845 4822 : vmax = 0.0_dp
846 4822 : ik = 0
847 : !select index ik
848 89798 : DO jt = 1, ndim
849 89798 : IF (itaken(jt) == 0 .AND. ABS(q_mat(it, jt)) > vmax) THEN
850 11064 : vmax = ABS(q_mat(it, jt))
851 11064 : ik = jt
852 : END IF
853 : END DO
854 4822 : IF (ik == 0) CYCLE ! all candidates NaN or zero: keep the entry vector
855 4822 : itaken(ik) = 1
856 :
857 8150 : CALL cp_fm_to_fm(c2_tmp, v_mat(1), 1, ik, it)
858 : END DO
859 3328 : DEALLOCATE (itaken)
860 3328 : DEALLOCATE (q_mat)
861 :
862 : !Copy in the converged set to enlarge the converged subspace
863 3328 : CALL cp_fm_to_fm(v_mat(1), c0, nmob, 1, (nmo_converged + imo_low))
864 3328 : CALL timestop(hand4)
865 :
866 3328 : CALL cp_fm_release(c2_tmp)
867 3328 : CALL cp_fm_release(c_tmp)
868 3328 : CALL cp_fm_release(hc)
869 3328 : CALL cp_fm_release(v_tmp)
870 3328 : CALL cp_fm_release(b_mat)
871 :
872 3328 : DEALLOCATE (t_mat)
873 3328 : DEALLOCATE (a_block)
874 3328 : DEALLOCATE (b_block)
875 :
876 30686 : CALL cp_fm_release(v_mat)
877 :
878 : END DO ! ib
879 :
880 734 : CALL timeset(routineN//"_ortho", hand5)
881 734 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, c0, rzero, chc)
882 :
883 734 : CALL cp_fm_cholesky_decompose(chc, nmo)
884 734 : CALL cp_fm_triangular_multiply(chc, c0, 'R', invert_tr=.TRUE.)
885 734 : CALL timestop(hand5)
886 : ELSE
887 : ! Do nothing
888 : END IF
889 :
890 744 : CALL timestop(handle)
891 1488 : END SUBROUTINE lanczos_refinement_2v
892 :
893 : END MODULE qs_scf_lanczos
|
