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 Complex block-Lanczos refinement kernel for k-point SCF
10 : !> diagonalization. Processes unconverged columns in blocks of
11 : !> NBLOCK with two-pass DGKS reorthogonalization against the
12 : !> working basis and deflation prefix, BLAS-3 CholeskyQR2 batch
13 : !> orthonormalization (with CGS2 fallback on rank-deficient
14 : !> batches), and an explicit Hermitian Rayleigh-Ritz projection
15 : !> T = V^dag(A V). Scratch matrices and local LAPACK buffers live
16 : !> in a caller-owned workspace reused across sweeps.
17 : ! **************************************************************************************************
18 : MODULE qs_scf_lanczos_c
19 : USE cp_blacs_env, ONLY: cp_blacs_env_type
20 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
21 : cp_cfm_gemm,&
22 : cp_cfm_triangular_multiply,&
23 : cp_cfm_vectorsnorm
24 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
25 : USE cp_cfm_diag, ONLY: cp_cfm_heevd
26 : USE cp_cfm_types, ONLY: cp_cfm_create,&
27 : cp_cfm_get_info,&
28 : cp_cfm_release,&
29 : cp_cfm_to_cfm,&
30 : cp_cfm_type
31 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
32 : cp_fm_struct_release,&
33 : cp_fm_struct_type
34 : USE kinds, ONLY: dp
35 : USE mathconstants, ONLY: z_mone,&
36 : z_one,&
37 : z_zero
38 : USE message_passing, ONLY: mp_para_env_type
39 : #include "./base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 : PRIVATE
43 :
44 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_lanczos_c'
45 : INTEGER, PARAMETER, PRIVATE :: shared_budget_factor = 4
46 :
47 : PUBLIC :: lanczos_refinement_2v_c
48 : PUBLIC :: lanczos_channel_begin, lanczos_sweep_suffix, lanczos_channel_end
49 : PUBLIC :: lanczos_workspace_create, lanczos_workspace_release
50 :
51 : ! **************************************************************************************************
52 : !> \brief Caller-owned scratch of the refinement kernel: every matrix and
53 : !> array that the sweeps of one channel reuse. The driver creates
54 : !> the workspace once per diagonalization call and hands it to
55 : !> every kernel call. The widths follow the worst case of the run.
56 : !> The basis budget holds the larger of the NKRYLOV block budget
57 : !> and four times the widest single block, and every per-call
58 : !> clamp stays below it. The small Rayleigh-Ritz pair keeps its
59 : !> per-block right-sized allocation, because the solver stack
60 : !> diagonalizes whole matrices and knows no leading-block argument
61 : ! **************************************************************************************************
62 : TYPE, PUBLIC :: lanczos_workspace_type
63 : ! entry Rayleigh-Ritz, its rotation and the shared scratch
64 : TYPE(cp_cfm_type) :: hmo, mo_gram, evec, rot_scr
65 : ! working basis and its A-images
66 : TYPE(cp_cfm_type) :: vmat, hvmat
67 : ! candidate block, within-batch scratch column, batch Gram, deflation and
68 : ! projection-coefficient scratch
69 : TYPE(cp_cfm_type) :: cand_blk, one_col, batch_gram, proj_coeff, defl_coeff
70 : REAL(dp), ALLOCATABLE :: theta(:), nrm(:), rank_scale(:), pair_res(:), entry_res(:)
71 : LOGICAL, ALLOCATABLE :: ldone(:)
72 : COMPLEX(dp), ALLOCATABLE :: eval_scale(:)
73 : END TYPE lanczos_workspace_type
74 :
75 : CONTAINS
76 :
77 : ! **************************************************************************************************
78 : !> \brief Refine the non-converged columns of c0 on the standardized operator ks.
79 : !> The entry Rayleigh-Ritz pass updates the convergence statistics.
80 : !> Every non-converged column block then runs a single-pass
81 : !> Lanczos refinement. A full Cholesky sweep closes the call.
82 : !> \param ks standardized operator A = U^-dag H U^-1, nao x nao
83 : !> \param c0 current MOS in the standardized basis, refined in place.
84 : !> The matrix dimensions of the call follow c0
85 : !> \param eval eigenvalues of c0, updated in place by the entry Rayleigh-Ritz
86 : !> \param eps_col per-column convergence thresholds on the residual norms:
87 : !> the driver assigns the relaxed empty threshold to the unoccupied
88 : !> columns. The occupation decides, not the column index. This
89 : !> mirrors the Davidson k-point driver
90 : !> \param nkrylov number of block expansion steps in the block regime.
91 : !> The kernel uses it as a floor on the budget only when NBLOCK
92 : !> covers every pending column
93 : !> \param nblock number of columns refined jointly, one or more
94 : !> \param ws caller-owned scratch of the call, created by the driver at
95 : !> the worst-case widths of the run
96 : !> \param max_res_norm largest residual column norm on exit: the entry
97 : !> Rayleigh-Ritz values, or the wanted-pair residuals of the
98 : !> refinement blocks when a refinement ran
99 : !> \param min_res_norm smallest residual column norm on exit (same basis)
100 : !> \param nmo_conv number of leading converged columns on exit
101 : !> \param nmo_nc columns outside the leading converged prefix: every
102 : !> still-unconverged column, plus converged columns beyond a gap.
103 : !> The count stays conservative for the not-converged reporting
104 : !> \param check_moconv_only only run the entry Rayleigh-Ritz and the
105 : !> convergence count, skip the refinement
106 : ! **************************************************************************************************
107 0 : SUBROUTINE lanczos_refinement_2v_c(ks, c0, eval, eps_col, nkrylov, nblock, ws, &
108 : max_res_norm, min_res_norm, nmo_conv, nmo_nc, check_moconv_only)
109 :
110 : TYPE(cp_cfm_type), INTENT(IN) :: ks
111 : TYPE(cp_cfm_type), INTENT(INOUT) :: c0
112 : REAL(dp), DIMENSION(:), POINTER :: eval
113 : REAL(dp), DIMENSION(:), INTENT(IN) :: eps_col
114 : INTEGER, INTENT(IN) :: nkrylov, nblock
115 : TYPE(lanczos_workspace_type), INTENT(INOUT) :: ws
116 : REAL(dp), INTENT(OUT) :: max_res_norm, min_res_norm
117 : INTEGER, INTENT(OUT) :: nmo_conv, nmo_nc
118 : LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
119 :
120 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lanczos_refinement_2v_c'
121 :
122 : INTEGER :: handle
123 : LOGICAL :: my_check_moconv_only, next_suffix_rr, &
124 : refined
125 :
126 0 : CALL timeset(routineN, handle)
127 :
128 0 : my_check_moconv_only = .FALSE.
129 0 : IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
130 :
131 : CALL lanczos_channel_begin(ks, c0, eval, eps_col, ws, max_res_norm, min_res_norm, &
132 0 : nmo_conv, nmo_nc, check_moconv_only=my_check_moconv_only)
133 :
134 0 : IF (.NOT. my_check_moconv_only) THEN
135 0 : refined = .FALSE.
136 0 : IF (nmo_nc > 0) THEN
137 : CALL lanczos_sweep_suffix(ks, c0, eval, eps_col, nkrylov, nblock, .TRUE., .FALSE., &
138 : ws, max_res_norm, min_res_norm, nmo_conv, nmo_nc, &
139 0 : refined, next_suffix_rr)
140 : END IF
141 0 : IF (refined) THEN
142 0 : CALL lanczos_channel_end(c0, ws)
143 : END IF
144 : END IF
145 :
146 0 : CALL timestop(handle)
147 0 : END SUBROUTINE lanczos_refinement_2v_c
148 :
149 : ! **************************************************************************************************
150 : !> \brief Entry preparation of one (k-point, spin) channel before its inner
151 : !> Lanczos sweeps: orthonormalizes the input MOS (skipped for a
152 : !> read-only convergence probe), runs the full-width entry
153 : !> Rayleigh-Ritz pass, and computes the initial leading converged
154 : !> prefix and residual bounds.
155 : !> \param ks standardized operator A = U^-dag H U^-1, nao x nao
156 : !> \param c0 current MOS in the standardized basis, rotated in place
157 : !> \param eval eigenvalues of c0, updated in place
158 : !> \param eps_col per-column convergence thresholds
159 : !> \param ws caller-owned channel workspace
160 : !> \param max_res_norm largest entry residual norm
161 : !> \param min_res_norm smallest entry residual norm
162 : !> \param nmo_conv number of leading converged columns on entry
163 : !> \param nmo_nc number of columns outside the leading converged prefix
164 : !> \param check_moconv_only skip the entry Cholesky orthonormalization for
165 : !> a read-only probe
166 : ! **************************************************************************************************
167 764 : SUBROUTINE lanczos_channel_begin(ks, c0, eval, eps_col, ws, max_res_norm, min_res_norm, &
168 : nmo_conv, nmo_nc, check_moconv_only)
169 : TYPE(cp_cfm_type), INTENT(IN) :: ks
170 : TYPE(cp_cfm_type), INTENT(INOUT) :: c0
171 : REAL(dp), DIMENSION(:), POINTER :: eval
172 : REAL(dp), DIMENSION(:), INTENT(IN) :: eps_col
173 : TYPE(lanczos_workspace_type), INTENT(INOUT) :: ws
174 : REAL(dp), INTENT(OUT) :: max_res_norm, min_res_norm
175 : INTEGER, INTENT(OUT) :: nmo_conv, nmo_nc
176 : LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
177 :
178 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lanczos_channel_begin'
179 :
180 : INTEGER :: handle, info, iwork, nao, nmo
181 : LOGICAL :: my_check_moconv_only
182 :
183 764 : CALL timeset(routineN, handle)
184 :
185 764 : CALL cp_cfm_get_info(c0, nrow_global=nao, ncol_global=nmo)
186 :
187 764 : my_check_moconv_only = .FALSE.
188 764 : IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
189 :
190 764 : IF (.NOT. my_check_moconv_only) THEN
191 764 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, c0, c0, z_zero, ws%mo_gram)
192 764 : CALL cp_cfm_cholesky_decompose(ws%mo_gram, nmo, info)
193 764 : IF (info /= 0) THEN
194 : CALL cp_abort(__LOCATION__, &
195 0 : "The Gram matrix of the entry MOS lost positive definiteness.")
196 : END IF
197 764 : CALL cp_cfm_triangular_multiply(ws%mo_gram, c0, 'R', invert_tr=.TRUE.)
198 : END IF
199 :
200 764 : CALL cp_cfm_gemm('N', 'N', nao, nmo, nao, z_one, ks, c0, z_zero, ws%hmo)
201 : CALL rayleigh_ritz(ws%hmo, c0, nmo, nmo, eval, ws%mo_gram, ws%evec, &
202 764 : ws%eval_scale, ws%rot_scr, 1, ws%entry_res)
203 :
204 764 : nmo_conv = 0
205 764 : nmo_nc = 0
206 7916 : max_res_norm = MAXVAL(ws%entry_res(1:nmo))
207 7916 : min_res_norm = MINVAL(ws%entry_res(1:nmo))
208 2509 : DO iwork = 1, nmo
209 2509 : IF (ws%entry_res(iwork) <= eps_col(iwork)) THEN
210 1745 : nmo_conv = nmo_conv + 1
211 : ELSE
212 594 : nmo_nc = nmo - nmo_conv
213 594 : EXIT
214 : END IF
215 : END DO
216 :
217 764 : CALL timestop(handle)
218 764 : END SUBROUTINE lanczos_channel_begin
219 :
220 : ! **************************************************************************************************
221 : !> \brief Run one inner Lanczos sweep on the unconverged suffix of c0.
222 : !> When the previous sweep ran in multi-block mode or had a
223 : !> truncated block, a local suffix Rayleigh-Ritz pass first
224 : !> decouples the unconverged tail of width nmo_nc without touching
225 : !> the locked prefix 1..nmo_conv.
226 : !> \param ks standardized operator, nao x nao
227 : !> \param c0 MOS, refined in place on columns nmo_conv+1..nmo
228 : !> \param eval eigenvalues, updated in place on the unconverged suffix
229 : !> \param eps_col per-column convergence thresholds
230 : !> \param nkrylov number of block expansion steps in the block regime
231 : !> \param nblock number of columns refined jointly
232 : !> \param first_sweep true on the first sweep right after lanczos_channel_begin
233 : !> \param need_suffix_rr whether to run a local suffix Rayleigh-Ritz before
234 : !> refining
235 : !> \param ws caller-owned channel workspace
236 : !> \param max_res_norm largest residual norm of the active suffix on exit
237 : !> \param min_res_norm smallest residual norm of the active suffix on exit
238 : !> \param nmo_conv leading converged column count, advanced in place
239 : !> \param nmo_nc unconverged suffix width, updated in place
240 : !> \param step_refined whether this sweep modified c0
241 : !> \param next_need_suffix_rr whether the next sweep requires a suffix RR
242 : ! **************************************************************************************************
243 839 : SUBROUTINE lanczos_sweep_suffix(ks, c0, eval, eps_col, nkrylov, nblock, first_sweep, &
244 : need_suffix_rr, ws, max_res_norm, min_res_norm, &
245 : nmo_conv, nmo_nc, step_refined, next_need_suffix_rr)
246 : TYPE(cp_cfm_type), INTENT(IN) :: ks
247 : TYPE(cp_cfm_type), INTENT(INOUT) :: c0
248 : REAL(dp), DIMENSION(:), POINTER :: eval
249 : REAL(dp), DIMENSION(:), INTENT(IN) :: eps_col
250 : INTEGER, INTENT(IN) :: nkrylov, nblock
251 : LOGICAL, INTENT(IN) :: first_sweep, need_suffix_rr
252 : TYPE(lanczos_workspace_type), INTENT(INOUT) :: ws
253 : REAL(dp), INTENT(INOUT) :: max_res_norm, min_res_norm
254 : INTEGER, INTENT(INOUT) :: nmo_conv, nmo_nc
255 : LOGICAL, INTENT(OUT) :: step_refined, next_need_suffix_rr
256 :
257 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lanczos_sweep_suffix'
258 :
259 : INTEGER :: handle, n_conv_cols, nao, nmo_converged, &
260 : nmo_nc_loc
261 : LOGICAL :: refined, skip_first_defl
262 : REAL(dp) :: prev_max, prev_min
263 :
264 839 : CALL timeset(routineN, handle)
265 :
266 839 : CALL cp_cfm_get_info(c0, nrow_global=nao)
267 :
268 839 : step_refined = .FALSE.
269 839 : next_need_suffix_rr = .FALSE.
270 :
271 839 : IF (need_suffix_rr .AND. nmo_nc > 0) THEN
272 : CALL suffix_rayleigh_ritz(ks, c0, eval, eps_col, ws, nmo_conv, nmo_nc, &
273 193 : max_res_norm, min_res_norm)
274 193 : step_refined = .TRUE.
275 193 : IF (nmo_nc == 0) THEN
276 0 : CALL timestop(handle)
277 0 : RETURN
278 : END IF
279 : END IF
280 :
281 839 : IF (nmo_nc > 0) THEN
282 839 : prev_max = max_res_norm
283 839 : prev_min = min_res_norm
284 839 : nmo_converged = nmo_conv
285 839 : nmo_nc_loc = nmo_nc
286 839 : skip_first_defl = first_sweep .OR. need_suffix_rr
287 : CALL lanczos_refine(ks, c0, eval, nao, eps_col, nkrylov, nblock, &
288 : nmo_converged, nmo_nc_loc, skip_first_defl, n_conv_cols, &
289 839 : max_res_norm, min_res_norm, refined, next_need_suffix_rr, ws)
290 839 : IF (.NOT. refined) THEN
291 0 : max_res_norm = prev_max
292 0 : min_res_norm = prev_min
293 : ELSE
294 839 : step_refined = .TRUE.
295 : END IF
296 839 : nmo_conv = nmo_converged + n_conv_cols
297 839 : nmo_nc = nmo_nc_loc - n_conv_cols
298 : END IF
299 :
300 839 : CALL timestop(handle)
301 : END SUBROUTINE lanczos_sweep_suffix
302 :
303 : ! **************************************************************************************************
304 : !> \brief Defensive full-width Cholesky orthonormalization at the end of a
305 : !> channel's refinement sweeps, executed at most once per channel.
306 : !> \param c0 refined MOS in the standardized basis
307 : !> \param ws caller-owned channel workspace
308 : ! **************************************************************************************************
309 1782 : SUBROUTINE lanczos_channel_end(c0, ws)
310 : TYPE(cp_cfm_type), INTENT(INOUT) :: c0
311 : TYPE(lanczos_workspace_type), INTENT(INOUT) :: ws
312 :
313 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lanczos_channel_end'
314 :
315 : INTEGER :: handle, info, nao, nmo
316 :
317 594 : CALL timeset(routineN, handle)
318 :
319 594 : CALL cp_cfm_get_info(c0, nrow_global=nao, ncol_global=nmo)
320 594 : CALL cp_cfm_gemm('C', 'N', nmo, nmo, nao, z_one, c0, c0, z_zero, ws%mo_gram)
321 594 : CALL cp_cfm_cholesky_decompose(ws%mo_gram, nmo, info)
322 594 : IF (info /= 0) THEN
323 : CALL cp_abort(__LOCATION__, &
324 0 : "The Gram matrix of the refined set lost positive definiteness.")
325 : END IF
326 594 : CALL cp_cfm_triangular_multiply(ws%mo_gram, c0, 'R', invert_tr=.TRUE.)
327 :
328 594 : CALL timestop(handle)
329 594 : END SUBROUTINE lanczos_channel_end
330 :
331 : ! **************************************************************************************************
332 : !> \brief Orthonormalize and run a local Rayleigh-Ritz pass on the
333 : !> unconverged suffix c0(:, nmo_conv+1:nmo) after a multi-block or
334 : !> truncated sweep, advancing nmo_conv if newly decoupled pairs meet
335 : !> eps_col.
336 : !> \param ks standardized operator, nao x nao
337 : !> \param c0 MOS, updated in place on the suffix nmo_conv+1..nmo
338 : !> \param eval eigenvalues, updated in place on the suffix
339 : !> \param eps_col per-column convergence thresholds
340 : !> \param ws caller-owned channel workspace
341 : !> \param nmo_conv leading converged count, advanced in place
342 : !> \param nmo_nc unconverged suffix width, updated in place
343 : !> \param max_res_norm largest residual norm of the suffix
344 : !> \param min_res_norm smallest residual norm of the suffix
345 : ! **************************************************************************************************
346 193 : SUBROUTINE suffix_rayleigh_ritz(ks, c0, eval, eps_col, ws, nmo_conv, nmo_nc, &
347 : max_res_norm, min_res_norm)
348 : TYPE(cp_cfm_type), INTENT(IN) :: ks
349 : TYPE(cp_cfm_type), INTENT(INOUT) :: c0
350 : REAL(dp), DIMENSION(:), POINTER :: eval
351 : REAL(dp), DIMENSION(:), INTENT(IN) :: eps_col
352 : TYPE(lanczos_workspace_type), INTENT(INOUT) :: ws
353 : INTEGER, INTENT(INOUT) :: nmo_conv, nmo_nc
354 : REAL(dp), INTENT(OUT) :: max_res_norm, min_res_norm
355 :
356 : INTEGER :: col_lo, info, ipass, iwork, m, nao, nmo
357 : TYPE(cp_cfm_type) :: esolve, hsolve
358 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
359 :
360 193 : NULLIFY (fm_struct_tmp)
361 193 : CALL cp_cfm_get_info(c0, nrow_global=nao, ncol_global=nmo)
362 :
363 193 : col_lo = nmo_conv + 1
364 193 : m = nmo_nc
365 :
366 193 : CALL cp_cfm_to_cfm(c0, ws%rot_scr, m, col_lo, 1)
367 193 : IF (nmo_conv > 0) THEN
368 291 : DO ipass = 1, 2
369 194 : CALL cp_cfm_gemm('C', 'N', nmo_conv, m, nao, z_one, c0, ws%rot_scr, z_zero, ws%mo_gram)
370 291 : CALL cp_cfm_gemm('N', 'N', nao, m, nmo_conv, z_mone, c0, ws%mo_gram, z_one, ws%rot_scr)
371 : END DO
372 : END IF
373 193 : CALL cp_cfm_gemm('C', 'N', m, m, nao, z_one, ws%rot_scr, ws%rot_scr, z_zero, ws%mo_gram)
374 193 : CALL cp_cfm_cholesky_decompose(ws%mo_gram, m, info)
375 193 : IF (info /= 0) THEN
376 : CALL cp_abort(__LOCATION__, &
377 0 : "The Gram matrix of the suffix MOS lost positive definiteness.")
378 : END IF
379 : CALL cp_cfm_triangular_multiply(ws%mo_gram, ws%rot_scr, side='R', invert_tr=.TRUE., &
380 193 : n_rows=nao, n_cols=m)
381 193 : CALL cp_cfm_to_cfm(ws%rot_scr, c0, m, 1, col_lo)
382 :
383 193 : CALL cp_cfm_gemm('N', 'N', nao, m, nao, z_one, ks, c0, z_zero, ws%hmo, b_first_col=col_lo)
384 193 : IF (m == nmo) THEN
385 : CALL rayleigh_ritz(ws%hmo, c0, m, m, eval(col_lo:nmo), ws%mo_gram, ws%evec, &
386 96 : ws%eval_scale, ws%rot_scr, col_lo, ws%entry_res, bas_col=col_lo)
387 : ELSE
388 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=m, ncol_global=m, &
389 : para_env=c0%matrix_struct%para_env, &
390 97 : context=c0%matrix_struct%context)
391 97 : CALL cp_cfm_create(hsolve, matrix_struct=fm_struct_tmp, name="lan_sh", set_zero=.TRUE.)
392 97 : CALL cp_cfm_create(esolve, matrix_struct=fm_struct_tmp, name="lan_se", set_zero=.TRUE.)
393 97 : CALL cp_fm_struct_release(fm_struct_tmp)
394 : CALL rayleigh_ritz(ws%hmo, c0, m, m, eval(col_lo:nmo), hsolve, esolve, &
395 97 : ws%eval_scale, ws%rot_scr, col_lo, ws%entry_res, bas_col=col_lo)
396 97 : CALL cp_cfm_release(hsolve)
397 97 : CALL cp_cfm_release(esolve)
398 : END IF
399 :
400 2556 : max_res_norm = MAXVAL(ws%entry_res(1:m))
401 2556 : min_res_norm = MINVAL(ws%entry_res(1:m))
402 195 : DO iwork = 1, m
403 195 : IF (ws%entry_res(iwork) <= eps_col(col_lo + iwork - 1)) THEN
404 2 : nmo_conv = nmo_conv + 1
405 : ELSE
406 : EXIT
407 : END IF
408 : END DO
409 193 : nmo_nc = nmo - nmo_conv
410 193 : END SUBROUTINE suffix_rayleigh_ritz
411 :
412 : ! **************************************************************************************************
413 : !> \brief Per-block single-pass Lanczos refinement on the unconverged tail.
414 : !> Processes non-converged columns in blocks of nblock, expanding
415 : !> each seed block by A-images with two-pass DGKS reorthogonalization
416 : !> against the working basis and deflation prefix, and solving the
417 : !> projected Hermitian Rayleigh-Ritz problem T = V^dag(A V).
418 : !> \param ks standardized operator, nao x nao
419 : !> \param c0 MOS, refined in place block by block
420 : !> \param eval eigenvalues, updated per converged column
421 : !> \param nao number of atomic orbitals
422 : !> \param eps_col per-column convergence thresholds
423 : !> \param nkrylov number of block expansion steps in the block regime
424 : !> \param nblock number of columns refined jointly
425 : !> \param nmo_converged number of leading locked columns
426 : !> \param nmo_nc_loc number of non-converged tail columns
427 : !> \param skip_first_defl skip deflation on the first seed block when the
428 : !> suffix just underwent a Rayleigh-Ritz rotation
429 : !> \param n_conv_cols length of the leading converged run across blocks
430 : !> \param max_res_norm largest residual norm among the wanted pairs
431 : !> \param min_res_norm smallest residual norm among the wanted pairs
432 : !> \param refined whether any block completed a Rayleigh-Ritz pass
433 : !> \param need_suffix_rr whether a subsequent sweep must run a suffix RR
434 : !> (multi-block run or any block truncated below its target width)
435 : !> \param ws caller-owned scratch of the call. Only the small
436 : !> Rayleigh-Ritz pair stays local to a block, at its right size,
437 : !> because the solver stack diagonalizes whole matrices
438 : ! **************************************************************************************************
439 839 : SUBROUTINE lanczos_refine(ks, c0, eval, nao, eps_col, nkrylov, nblock, &
440 : nmo_converged, nmo_nc_loc, skip_first_defl, n_conv_cols, &
441 : max_res_norm, min_res_norm, refined, need_suffix_rr, ws)
442 : TYPE(cp_cfm_type), INTENT(IN) :: ks
443 : TYPE(cp_cfm_type), INTENT(INOUT) :: c0
444 : REAL(dp), DIMENSION(:), POINTER :: eval
445 : INTEGER, INTENT(IN) :: nao
446 : REAL(dp), DIMENSION(:), INTENT(IN) :: eps_col
447 : INTEGER, INTENT(IN) :: nkrylov, nblock, nmo_converged, &
448 : nmo_nc_loc
449 : LOGICAL, INTENT(IN) :: skip_first_defl
450 : INTEGER, INTENT(OUT) :: n_conv_cols
451 : REAL(dp), INTENT(OUT) :: max_res_norm, min_res_norm
452 : LOGICAL, INTENT(OUT) :: refined, need_suffix_rr
453 : TYPE(lanczos_workspace_type), INTENT(INOUT) :: ws
454 :
455 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lanczos_refine'
456 :
457 : INTEGER :: col_hi, col_lo, handle, ib, j, nbasis, &
458 : nblock_e, nbmax, nmob_c, nmobt, nres, &
459 : nseed, nstep
460 : LOGICAL :: any_truncated, prefix_open, single_block
461 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
462 : TYPE(cp_cfm_type) :: esolve, hsolve
463 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
464 : TYPE(mp_para_env_type), POINTER :: para_env
465 :
466 839 : CALL timeset(routineN, handle)
467 :
468 839 : NULLIFY (fm_struct_tmp, blacs_env, para_env)
469 839 : para_env => c0%matrix_struct%para_env
470 839 : blacs_env => c0%matrix_struct%context
471 :
472 839 : nblock_e = MIN(MAX(nblock, 1), MAX(nmo_nc_loc, 1))
473 839 : single_block = nblock_e >= nmo_nc_loc
474 839 : IF (single_block) THEN
475 : ! one block covers every pending column: the Krylov budget holds
476 : ! about four times the pending width, floored at NKRYLOV and at
477 : ! the block width. The shared subspace holds every target pair
478 : ! plus a threefold expansion room. NKRYLOV plays no role here
479 370 : nbmax = MIN(nao - 2, MAX(shared_budget_factor*nmo_nc_loc, nkrylov, nblock_e + 2))
480 : ELSE
481 469 : nbmax = nkrylov*nblock_e
482 : END IF
483 839 : IF (nbmax < 2) THEN
484 : ! no working basis beyond the seed exists (nao below four, or a
485 : ! degenerate NKRYLOV/NBLOCK combination): leave the columns at
486 : ! their entry state for the caller to report, instead of aborting
487 0 : n_conv_cols = 0
488 0 : max_res_norm = 0.0_dp
489 0 : min_res_norm = HUGE(1.0_dp)
490 0 : refined = .FALSE.
491 0 : need_suffix_rr = .FALSE.
492 0 : CPWARN("Krylov budget below two columns, Lanczos refinement skipped")
493 0 : CALL timestop(handle)
494 0 : RETURN
495 : END IF
496 : ! the budget bounds the block width: a wider block would append
497 : ! past the end of vmat at its first seed call (nao - 2 caps nbmax
498 : ! below the pending width on minimal bases). The pending set then
499 : ! splits into narrower blocks, and no column is dropped
500 839 : IF (nblock_e > nbmax - 2) THEN
501 110 : nblock_e = MAX(nbmax - 2, 1)
502 110 : single_block = nblock_e >= nmo_nc_loc
503 : END IF
504 :
505 839 : n_conv_cols = 0
506 839 : max_res_norm = 0.0_dp
507 839 : min_res_norm = HUGE(1.0_dp)
508 839 : prefix_open = .TRUE.
509 839 : refined = .FALSE.
510 839 : any_truncated = .FALSE.
511 :
512 7236 : block_loop: DO ib = 1, nmo_nc_loc, nblock_e
513 6397 : col_lo = nmo_converged + ib
514 6397 : col_hi = MIN(col_lo + nblock_e - 1, nmo_converged + nmo_nc_loc)
515 6397 : nmobt = col_hi - col_lo + 1
516 14217 : ws%ldone(1:nmobt) = .FALSE.
517 : nres = 0
518 :
519 : ! ---------- seed block: the pending columns, purified and
520 : ! orthonormalized into the working basis
521 6397 : CALL cp_cfm_to_cfm(c0, ws%cand_blk, nmobt, col_lo, 1)
522 6397 : CALL cp_cfm_vectorsnorm(ws%cand_blk, ws%nrm)
523 14217 : ws%rank_scale(1:nmobt) = ws%nrm(1:nmobt)
524 6397 : nbasis = 0
525 : CALL purify_and_append(c0, col_lo - 1, ib == 1 .AND. skip_first_defl, ws, &
526 12007 : MIN(nmobt, nbmax), nbasis, nseed)
527 6397 : IF (nseed == 0) THEN
528 : ! no seed survived the purification: no column of this block
529 : ! converged, and the leading prefix cannot extend past it
530 : prefix_open = .FALSE.
531 : any_truncated = .TRUE.
532 : CYCLE block_loop
533 : END IF
534 6397 : nstep = nseed
535 :
536 : ! ---------- expansion: extend the basis with the A-images of the
537 : ! last batch until the budget is exhausted
538 50658 : DO WHILE (nbasis < nbmax .AND. nstep > 0)
539 : CALL cp_cfm_gemm('N', 'N', nao, nstep, nao, z_one, ks, ws%vmat, z_zero, ws%hvmat, &
540 44261 : b_first_col=nbasis - nstep + 1, c_first_col=nbasis - nstep + 1)
541 44261 : CALL cp_cfm_to_cfm(ws%hvmat, ws%cand_blk, nstep, nbasis - nstep + 1, 1)
542 44261 : CALL cp_cfm_vectorsnorm(ws%cand_blk, ws%nrm)
543 97087 : ws%rank_scale(1:nstep) = ws%nrm(1:nstep)
544 44261 : nmob_c = MIN(nstep, nbmax - nbasis)
545 44261 : CALL purify_and_append(c0, col_lo - 1, .FALSE., ws, nmob_c, nbasis, nseed)
546 44365 : nstep = nseed
547 : END DO
548 : ! the last batch has no successor: its A-images must still be
549 : ! formed explicitly, or the projection reads the stale tail of
550 : ! the reused product buffer and silently degrades the Ritz pairs
551 6397 : IF (nstep > 0) THEN
552 : CALL cp_cfm_gemm('N', 'N', nao, nstep, nao, z_one, ks, ws%vmat, z_zero, ws%hvmat, &
553 6293 : b_first_col=nbasis - nstep + 1, c_first_col=nbasis - nstep + 1)
554 : END IF
555 :
556 : ! ---------- Rayleigh-Ritz: T = V^dag (A V) as a distributed
557 : ! problem on a right-sized matrix, through the shared solver
558 : ! stack
559 6397 : NULLIFY (fm_struct_tmp)
560 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nbasis, ncol_global=nbasis, &
561 6397 : para_env=para_env, context=blacs_env)
562 6397 : CALL cp_cfm_create(hsolve, matrix_struct=fm_struct_tmp, name="lan_h", set_zero=.TRUE.)
563 6397 : CALL cp_cfm_create(esolve, matrix_struct=fm_struct_tmp, name="lan_e", set_zero=.TRUE.)
564 6397 : CALL cp_fm_struct_release(fm_struct_tmp)
565 6397 : nres = MIN(nmobt, nbasis)
566 6397 : IF (nres < nmobt) any_truncated = .TRUE.
567 : ! one shared Rayleigh-Ritz pass per block, on the stored
568 : ! A-images of the working basis: rotates the leading nres Ritz
569 : ! vectors into c0(:, col_lo:col_lo+nres-1) and forms their exact
570 : ! residual norms in ws%nrm
571 : CALL rayleigh_ritz(ws%hvmat, ws%vmat, nbasis, nres, ws%theta, hsolve, esolve, &
572 6397 : ws%eval_scale, ws%cand_blk, col_lo, ws%nrm, c_dst=c0)
573 14217 : ws%pair_res(1:nres) = ws%nrm(1:nres)
574 :
575 : ! ---------- write-back metadata: no pair-to-column identity ever
576 : ! latches across calls, heevd re-sorts ascending. Columns beyond
577 : ! nres keep their entry vectors and fall to the end-of-call
578 : ! fallback
579 14217 : DO j = 1, nres
580 7820 : eval(col_lo + j - 1) = ws%theta(j)
581 14217 : ws%ldone(j) = ws%pair_res(j) < eps_col(col_lo + j - 1)
582 : END DO
583 6397 : refined = .TRUE.
584 6397 : CALL cp_cfm_release(hsolve)
585 6397 : CALL cp_cfm_release(esolve)
586 :
587 : ! exit statistics: the residuals describe the state the columns
588 : ! are left in, not the worst value seen during the call. The
589 : ! converged count is the global leading converged run
590 14217 : DO j = 1, nres
591 7820 : max_res_norm = MAX(max_res_norm, ws%pair_res(j))
592 14217 : min_res_norm = MIN(min_res_norm, ws%pair_res(j))
593 : END DO
594 7236 : IF (prefix_open) THEN
595 8117 : DO j = 1, nmobt
596 4829 : IF (.NOT. ws%ldone(j)) EXIT
597 8117 : n_conv_cols = n_conv_cols + 1
598 : END DO
599 8117 : prefix_open = ALL(ws%ldone(1:nmobt))
600 : END IF
601 : END DO block_loop
602 :
603 839 : need_suffix_rr = (.NOT. single_block) .OR. any_truncated
604 :
605 839 : CALL timestop(handle)
606 839 : END SUBROUTINE lanczos_refine
607 :
608 : ! **************************************************************************************************
609 : !> \brief Purify a block of candidate directions and append the survivors to
610 : !> the working basis: projection passes against the working basis and
611 : !> the leading ndefl columns of c0, with the second pass decided by
612 : !> the DGKS norm-ratio test, modified Gram-Schmidt within the block,
613 : !> then a per-column rank safeguard relative to rank_scale,
614 : !> normalization and append. The count of appended columns is
615 : !> returned in napp. The kernel advances nbasis accordingly.
616 : !> \param c0 MOS holding the deflation columns
617 : !> \param ndefl number of leading c0 columns to deflate against
618 : !> \param skip_defl skip deflation when the seed block is already orthogonal
619 : !> to the deflation prefix (first block after a full-width RR)
620 : !> \param ws kernel workspace: cand_blk holds the candidate directions, nmob
621 : !> of them live. vmat is the working basis appended in place.
622 : !> rank_scale carries the per-column rank-safeguard scale, nrm the
623 : !> column norms, proj_coeff the projection coefficients, defl_coeff the
624 : !> deflation-projection coefficients and one_col the scratch column
625 : !> of the within-batch orthogonalization
626 : !> \param nmob number of live candidate columns
627 : !> \param nbasis number of basis columns on entry, advanced on exit
628 : !> \param napp number of appended columns
629 : ! **************************************************************************************************
630 50658 : SUBROUTINE purify_and_append(c0, ndefl, skip_defl, ws, nmob, nbasis, napp)
631 : TYPE(cp_cfm_type), INTENT(IN) :: c0
632 : INTEGER, INTENT(IN) :: ndefl
633 : LOGICAL, INTENT(IN) :: skip_defl
634 : TYPE(lanczos_workspace_type), INTENT(INOUT) :: ws
635 : INTEGER, INTENT(IN) :: nmob
636 : INTEGER, INTENT(INOUT) :: nbasis
637 : INTEGER, INTENT(OUT) :: napp
638 :
639 : REAL(dp), PARAMETER :: eta_dgks = 0.70710678118654752_dp
640 :
641 : COMPLEX(dp) :: col_scale1(1)
642 : INTEGER :: icol, info, ipass, iqr, kept, nao
643 : LOGICAL :: chol_ok, do_defl, do_purge
644 :
645 50658 : CALL cp_cfm_get_info(ws%cand_blk, nrow_global=nao)
646 :
647 : ! batch purge: two passes of block gemms project the candidate set
648 : ! against the working basis and the deflation prefix; the first seed
649 : ! block after a full-width rotation skips it, later seed blocks must
650 : ! deflate against the columns earlier blocks wrote
651 50658 : do_defl = (ndefl > 0) .AND. (.NOT. skip_defl)
652 50658 : do_purge = (nbasis > 0) .OR. do_defl
653 : IF (do_purge) THEN
654 138307 : DO ipass = 1, 2
655 94063 : IF (nbasis > 0) THEN
656 88505 : CALL cp_cfm_gemm('C', 'N', nbasis, nmob, nao, z_one, ws%vmat, ws%cand_blk, z_zero, ws%proj_coeff)
657 88505 : CALL cp_cfm_gemm('N', 'N', nao, nmob, nbasis, z_mone, ws%vmat, ws%proj_coeff, z_one, ws%cand_blk)
658 : END IF
659 94063 : IF (do_defl) THEN
660 87949 : CALL cp_cfm_gemm('C', 'N', ndefl, nmob, nao, z_one, c0, ws%cand_blk, z_zero, ws%defl_coeff)
661 87949 : CALL cp_cfm_gemm('N', 'N', nao, nmob, ndefl, z_mone, c0, ws%defl_coeff, z_one, ws%cand_blk)
662 : END IF
663 94063 : CALL cp_cfm_vectorsnorm(ws%cand_blk, ws%nrm)
664 143882 : IF (ipass == 1) THEN
665 : ! DGKS norm-ratio test (Daniel, Gragg, Kaufman and Stewart): a
666 : ! first pass that left every candidate norm within eta_dgks of
667 : ! its pre-purge norm already orthogonalized the batch to the
668 : ! projected subspace to working precision
669 56029 : IF (ALL(ws%nrm(1:nmob) >= eta_dgks*ws%rank_scale(1:nmob))) EXIT
670 : END IF
671 : END DO
672 : END IF
673 :
674 : ! rank safeguard: once the Krylov space saturates the available
675 : ! dimensions, the orthogonal complement is pure cancellation noise
676 : ! and a residual of 1e-16 passes a plain positivity test
677 50898 : IF (ALL(ws%nrm(1:nmob) <= 1.0E-12_dp*ws%rank_scale(1:nmob))) THEN
678 104 : napp = 0
679 50619 : RETURN
680 : END IF
681 :
682 : ! Path A: single-column fast path (norm guard already satisfied above)
683 50554 : IF (nmob == 1) THEN
684 46970 : ws%eval_scale(1) = CMPLX(1.0_dp/ws%nrm(1), 0.0_dp, KIND=dp)
685 858982 : IF (SIZE(ws%eval_scale) > 1) ws%eval_scale(2:) = z_zero
686 46970 : CALL cp_cfm_column_scale(ws%cand_blk, ws%eval_scale)
687 46970 : CALL cp_cfm_to_cfm(ws%cand_blk, ws%vmat, 1, 1, nbasis + 1)
688 46970 : napp = 1
689 46970 : nbasis = nbasis + 1
690 46970 : RETURN
691 : END IF
692 :
693 : ! Path B: healthy multi-column batch fast path via two-pass BLAS-3 CholeskyQR2
694 16975 : IF (ALL(ws%nrm(1:nmob) > 1.0E-12_dp*ws%rank_scale(1:nmob))) THEN
695 3579 : CALL cp_cfm_to_cfm(ws%cand_blk, ws%vmat, nmob, 1, nbasis + 1)
696 3579 : chol_ok = .TRUE.
697 10669 : DO iqr = 1, 2
698 : CALL cp_cfm_gemm('C', 'N', nmob, nmob, nao, z_one, ws%cand_blk, ws%cand_blk, &
699 7124 : z_zero, ws%batch_gram)
700 7124 : CALL cp_cfm_cholesky_decompose(ws%batch_gram, nmob, info)
701 7124 : IF (info /= 0) THEN
702 : chol_ok = .FALSE.
703 : EXIT
704 : END IF
705 : CALL cp_cfm_triangular_multiply(ws%batch_gram, ws%cand_blk, side='R', &
706 17759 : invert_tr=.TRUE., n_rows=nao, n_cols=nmob)
707 : END DO
708 3579 : IF (chol_ok) THEN
709 3545 : CALL cp_cfm_to_cfm(ws%cand_blk, ws%vmat, nmob, 1, nbasis + 1)
710 3545 : napp = nmob
711 3545 : nbasis = nbasis + nmob
712 3545 : RETURN
713 : END IF
714 34 : CALL cp_cfm_to_cfm(ws%vmat, ws%cand_blk, nmob, nbasis + 1, 1)
715 : END IF
716 :
717 : ! Path C: rank-deficient / near-dependent fallback by two-pass CGS2
718 : ! with per-column norm retirement
719 39 : kept = 0
720 602 : DO icol = 1, nmob
721 563 : CALL cp_cfm_to_cfm(ws%cand_blk, ws%one_col, 1, icol, 1)
722 1689 : DO ipass = 1, 2
723 1689 : IF (kept > 0) THEN
724 : CALL cp_cfm_gemm('C', 'N', kept, 1, nao, z_one, ws%cand_blk, &
725 1048 : ws%one_col, z_zero, ws%proj_coeff)
726 : CALL cp_cfm_gemm('N', 'N', nao, 1, kept, z_mone, ws%cand_blk, &
727 1048 : ws%proj_coeff, z_one, ws%one_col)
728 : END IF
729 : END DO
730 563 : CALL cp_cfm_vectorsnorm(ws%one_col, ws%nrm)
731 602 : IF (ws%nrm(1) > 1.0E-12_dp*ws%rank_scale(icol)) THEN
732 150 : col_scale1(1) = CMPLX(1.0_dp/ws%nrm(1), 0.0_dp, KIND=dp)
733 150 : CALL cp_cfm_column_scale(ws%one_col, col_scale1)
734 150 : kept = kept + 1
735 150 : CALL cp_cfm_to_cfm(ws%one_col, ws%cand_blk, 1, 1, kept)
736 : END IF
737 : END DO
738 39 : napp = kept
739 39 : IF (napp > 0) CALL cp_cfm_to_cfm(ws%cand_blk, ws%vmat, napp, 1, nbasis + 1)
740 39 : nbasis = nbasis + napp
741 : END SUBROUTINE purify_and_append
742 :
743 : ! **************************************************************************************************
744 : !> \brief One Rayleigh-Ritz pass on a precomputed basis pair: project
745 : !> the operator onto the basis, solve the small Hermitian
746 : !> problem, rotate the leading Ritz vectors into the target
747 : !> matrix window, and return their exact residual norms.
748 : !> \param abas A-images of the basis columns, A*B
749 : !> \param basis orthonormal basis B; updated in place when c_dst is absent
750 : !> \param nbas live basis width
751 : !> \param npairs number of leading Ritz pairs to write back and form residuals
752 : !> \param ritz_val eigenvalues of the projected problem, at least nbas
753 : !> \param small_t scratch holding the projected matrix on entry, right
754 : !> sized to nbas
755 : !> \param small_z scratch receiving the eigenvectors, right sized to nbas
756 : !> \param eval_scale complex scaling array spanning every global column
757 : !> of work_blk
758 : !> \param work_blk scratch receiving the rotated Ritz vectors and then
759 : !> the residual columns
760 : !> \param dst_col starting destination column in c_dst (or basis)
761 : !> \param res_norm residual norms, spans every global column of work_blk
762 : !> \param c_dst optional target matrix receiving the rotated Ritz vectors
763 : !> \param bas_col optional starting column of basis (defaults to 1)
764 : ! **************************************************************************************************
765 7354 : SUBROUTINE rayleigh_ritz(abas, basis, nbas, npairs, ritz_val, small_t, small_z, &
766 7354 : eval_scale, work_blk, dst_col, res_norm, c_dst, bas_col)
767 : TYPE(cp_cfm_type), INTENT(IN) :: abas
768 : TYPE(cp_cfm_type), INTENT(INOUT) :: basis
769 : INTEGER, INTENT(IN) :: nbas, npairs
770 : REAL(dp), DIMENSION(:) :: ritz_val
771 : TYPE(cp_cfm_type) :: small_t, small_z
772 : COMPLEX(dp), DIMENSION(:) :: eval_scale
773 : TYPE(cp_cfm_type) :: work_blk
774 : INTEGER, INTENT(IN) :: dst_col
775 : REAL(dp), DIMENSION(:) :: res_norm
776 : TYPE(cp_cfm_type), INTENT(INOUT), OPTIONAL :: c_dst
777 : INTEGER, INTENT(IN), OPTIONAL :: bas_col
778 :
779 : INTEGER :: my_bas_col, nao, ncol_work
780 :
781 7354 : my_bas_col = 1
782 7354 : IF (PRESENT(bas_col)) my_bas_col = bas_col
783 :
784 7354 : CALL cp_cfm_get_info(work_blk, nrow_global=nao, ncol_global=ncol_work)
785 :
786 : CALL cp_cfm_gemm('C', 'N', nbas, nbas, nao, z_one, basis, abas, z_zero, small_t, &
787 7354 : a_first_col=my_bas_col)
788 7354 : CALL cp_cfm_heevd(small_t, small_z, ritz_val(1:nbas))
789 :
790 : CALL cp_cfm_gemm('N', 'N', nao, npairs, nbas, z_one, basis, small_z, z_zero, work_blk, &
791 7354 : a_first_col=my_bas_col)
792 7354 : IF (PRESENT(c_dst)) THEN
793 6397 : CALL cp_cfm_to_cfm(work_blk, c_dst, npairs, 1, dst_col)
794 : ELSE
795 957 : CALL cp_cfm_to_cfm(work_blk, basis, npairs, 1, dst_col)
796 : END IF
797 :
798 24689 : eval_scale(1:npairs) = CMPLX(-ritz_val(1:npairs), 0.0_dp, KIND=dp)
799 9126 : IF (npairs < ncol_work) eval_scale(npairs + 1:ncol_work) = z_zero
800 7354 : CALL cp_cfm_column_scale(work_blk, eval_scale)
801 7354 : CALL cp_cfm_gemm('N', 'N', nao, npairs, nbas, z_one, abas, small_z, z_one, work_blk)
802 7354 : CALL cp_cfm_vectorsnorm(work_blk, res_norm)
803 7354 : END SUBROUTINE rayleigh_ritz
804 :
805 : ! **************************************************************************************************
806 : !> \brief Create the kernel workspace at the worst-case widths of the run.
807 : !> The basis budget holds the larger of the NKRYLOV block budget
808 : !> and four times the widest single block, and every per-call
809 : !> clamp of the kernel stays below it. A single block covers
810 : !> every pending column only when NBLOCK spans them, so four
811 : !> times MIN(NBLOCK, nmo) bounds that regime. The candidate and
812 : !> Ritz widths hold the widest block a run can form
813 : !> \param ws workspace to create
814 : !> \param mo_struct structure of the channel MOS, nao x nmo
815 : !> \param ao_struct structure of the channel operators, nao x nao
816 : !> \param nkrylov number of block expansion steps of the block regime
817 : !> \param nblock number of columns refined jointly
818 : ! **************************************************************************************************
819 586 : SUBROUTINE lanczos_workspace_create(ws, mo_struct, ao_struct, nkrylov, nblock)
820 : TYPE(lanczos_workspace_type), INTENT(OUT) :: ws
821 : TYPE(cp_fm_struct_type), POINTER :: mo_struct, ao_struct
822 : INTEGER, INTENT(IN) :: nkrylov, nblock
823 :
824 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lanczos_workspace_create'
825 :
826 : INTEGER :: handle, nao, nbe, nmo, wdth
827 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
828 :
829 586 : CALL timeset(routineN, handle)
830 586 : NULLIFY (fm_struct_tmp)
831 :
832 586 : nao = ao_struct%nrow_global
833 586 : nmo = mo_struct%ncol_global
834 586 : nbe = MIN(MAX(nblock, 1), nmo)
835 : ! a single block covering every pending column does not multiply
836 : ! the budget by the block count
837 586 : IF (MAX(nblock, 1) >= nmo) THEN
838 258 : wdth = MAX(shared_budget_factor*nbe, nkrylov, nbe + 2)
839 : ELSE
840 328 : wdth = MAX(nkrylov*nbe, shared_budget_factor*nbe, nbe + 2)
841 : END IF
842 :
843 586 : CALL cp_cfm_create(ws%hmo, mo_struct, name="lan_hmo")
844 586 : CALL cp_cfm_create(ws%rot_scr, mo_struct, name="lan_rot_scr")
845 :
846 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=nmo, &
847 586 : para_env=mo_struct%para_env, context=mo_struct%context)
848 586 : CALL cp_cfm_create(ws%mo_gram, matrix_struct=fm_struct_tmp, name="lan_mo_gram")
849 586 : CALL cp_cfm_create(ws%evec, matrix_struct=fm_struct_tmp, name="lan_evec")
850 586 : CALL cp_fm_struct_release(fm_struct_tmp)
851 :
852 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=wdth, &
853 586 : para_env=ao_struct%para_env, context=ao_struct%context)
854 586 : CALL cp_cfm_create(ws%vmat, matrix_struct=fm_struct_tmp, name="lan_v", set_zero=.TRUE.)
855 586 : CALL cp_cfm_create(ws%hvmat, matrix_struct=fm_struct_tmp, name="lan_hv", set_zero=.TRUE.)
856 586 : CALL cp_fm_struct_release(fm_struct_tmp)
857 :
858 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nbe, &
859 586 : para_env=ao_struct%para_env, context=ao_struct%context)
860 586 : CALL cp_cfm_create(ws%cand_blk, matrix_struct=fm_struct_tmp, name="lan_r", set_zero=.TRUE.)
861 586 : CALL cp_fm_struct_release(fm_struct_tmp)
862 :
863 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=1, &
864 586 : para_env=ao_struct%para_env, context=ao_struct%context)
865 586 : CALL cp_cfm_create(ws%one_col, matrix_struct=fm_struct_tmp, name="lan_w", set_zero=.TRUE.)
866 586 : CALL cp_fm_struct_release(fm_struct_tmp)
867 :
868 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nbe, ncol_global=nbe, &
869 586 : para_env=ao_struct%para_env, context=ao_struct%context)
870 586 : CALL cp_cfm_create(ws%batch_gram, matrix_struct=fm_struct_tmp, name="lan_bgram", set_zero=.TRUE.)
871 586 : CALL cp_fm_struct_release(fm_struct_tmp)
872 :
873 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=wdth, ncol_global=nbe, &
874 586 : para_env=ao_struct%para_env, context=ao_struct%context)
875 586 : CALL cp_cfm_create(ws%proj_coeff, matrix_struct=fm_struct_tmp, name="lan_pcoeff", set_zero=.TRUE.)
876 586 : CALL cp_fm_struct_release(fm_struct_tmp)
877 :
878 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=nbe, &
879 586 : para_env=mo_struct%para_env, context=mo_struct%context)
880 586 : CALL cp_cfm_create(ws%defl_coeff, matrix_struct=fm_struct_tmp, name="lan_pcb", set_zero=.TRUE.)
881 586 : CALL cp_fm_struct_release(fm_struct_tmp)
882 :
883 5274 : ALLOCATE (ws%theta(wdth), ws%nrm(nbe), ws%rank_scale(nbe), ws%pair_res(nbe), ws%entry_res(nmo))
884 1758 : ALLOCATE (ws%ldone(nbe))
885 1758 : ALLOCATE (ws%eval_scale(MAX(nmo, nbe)))
886 :
887 586 : CALL timestop(handle)
888 586 : END SUBROUTINE lanczos_workspace_create
889 :
890 : ! **************************************************************************************************
891 : !> \brief Release the kernel workspace. Safe on a workspace that was
892 : !> never created: the matrix components default to disassociated
893 : !> and the arrays to unallocated
894 : !> \param ws workspace to release
895 : ! **************************************************************************************************
896 36736 : SUBROUTINE lanczos_workspace_release(ws)
897 : TYPE(lanczos_workspace_type), INTENT(INOUT) :: ws
898 :
899 36736 : CALL cp_cfm_release(ws%hmo)
900 36736 : CALL cp_cfm_release(ws%mo_gram)
901 36736 : CALL cp_cfm_release(ws%evec)
902 36736 : CALL cp_cfm_release(ws%rot_scr)
903 36736 : CALL cp_cfm_release(ws%vmat)
904 36736 : CALL cp_cfm_release(ws%hvmat)
905 36736 : CALL cp_cfm_release(ws%cand_blk)
906 36736 : CALL cp_cfm_release(ws%one_col)
907 36736 : CALL cp_cfm_release(ws%batch_gram)
908 36736 : CALL cp_cfm_release(ws%proj_coeff)
909 36736 : CALL cp_cfm_release(ws%defl_coeff)
910 36736 : IF (ALLOCATED(ws%theta)) DEALLOCATE (ws%theta)
911 36736 : IF (ALLOCATED(ws%nrm)) DEALLOCATE (ws%nrm)
912 36736 : IF (ALLOCATED(ws%rank_scale)) DEALLOCATE (ws%rank_scale)
913 36736 : IF (ALLOCATED(ws%pair_res)) DEALLOCATE (ws%pair_res)
914 36736 : IF (ALLOCATED(ws%entry_res)) DEALLOCATE (ws%entry_res)
915 36736 : IF (ALLOCATED(ws%ldone)) DEALLOCATE (ws%ldone)
916 36736 : IF (ALLOCATED(ws%eval_scale)) DEALLOCATE (ws%eval_scale)
917 36736 : END SUBROUTINE lanczos_workspace_release
918 :
919 0 : END MODULE qs_scf_lanczos_c
|
