(git:f2099e5)
Loading...
Searching...
No Matches
bse_davidson.F
Go to the documentation of this file.
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 Block Davidson solvers for the lowest excitations of the Bethe-Salpeter equation on top
10!> of the matrix-free application of A and B: the TDA problem sum_jb A_ia,jb X_jb = Ω X_ia
11!> (bse_davidson_tda) and the full problem, either as sum_jb [(A+B)(A-B)]_ia,jb x_jb = Ω^2 x_ia
12!> (bse_davidson_abba_mk) or as the pair sum_jb (A+B)_ia,jb y_jb = Ω x_ia, sum_jb (A-B)_ia,jb x_jb = Ω y_ia
13!> (bse_davidson_abba_os). Every solver projects its problem onto an orthonormal basis, takes the
14!> lowest Ritz pairs, extends the basis by the preconditioned residuals of the roots still in need
15!> and restarts thick at the subspace ceiling; the helpers below are shared by the three drivers
16!> \par History
17!> 09.2026 created [Maximilian Graml]
18! **************************************************************************************************
19MODULE bse_davidson
20
21 USE bibliography, ONLY: bai2012,&
22 davidson1975,&
23 fukaya2014,&
24 olsen1988,&
25 stathopoulos1998,&
26 stratmann1998,&
27 vecharynski2017,&
28 cite_reference
29 USE bse_full_diag, ONLY: create_hermitian_form_of_abba
30 USE bse_matvec, ONLY: bse_matvec_apply,&
31 bse_matvec_diagonal,&
32 bse_matvec_env_type,&
33 bse_matvec_subblock,&
34 bse_matvec_vector_struct,&
35 mem_fraction
36 USE bse_types, ONLY: bse_env_type
37 USE cp_blacs_env, ONLY: cp_blacs_env_create,&
38 cp_blacs_env_release,&
39 cp_blacs_env_type
40 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
41 cp_fm_scale_and_add,&
42 cp_fm_transpose
43 USE cp_fm_diag, ONLY: choose_eigv_solver
44 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
45 cp_fm_struct_release,&
46 cp_fm_struct_type
47 USE cp_fm_types, ONLY: &
48 cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_release, cp_fm_set_all, &
49 cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_to_fm_submat, cp_fm_type
50 USE input_constants, ONLY: bse_iter_and_cond,&
51 bse_iter_en_cond,&
52 bse_iter_or_cond,&
53 bse_iter_res_cond,&
54 bse_memcheck_abort,&
55 bse_memcheck_clamp,&
56 bse_memcheck_off,&
57 bse_memcheck_warn
58 USE kinds, ONLY: dp
59 USE machine, ONLY: m_walltime
60 USE mathlib, ONLY: diamat_all
61 USE message_passing, ONLY: mp_mem_avail_per_rank_gb,&
62 mp_para_env_type
63 USE parallel_gemm_api, ONLY: parallel_gemm
64 USE physcon, ONLY: evolt
65 USE util, ONLY: sort
66#include "./base/base_uses.f90"
67
68 IMPLICIT NONE
69
70 PRIVATE
71
72 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_davidson'
73
74 ! energies closer than this belong to one multiplet (Hartree)
75 REAL(KIND=dp), PARAMETER, PRIVATE :: deg_thresh = 1.0e-6_dp
76 ! a basis that is not orthonormal to this accuracy makes the residuals meaningless
77 REAL(KIND=dp), PARAMETER, PRIVATE :: ortho_tol = 1.0e-8_dp
78
79 ! outcome of an ABBA solver: solved, or A-B not positive definite on the trial space
80 INTEGER, PARAMETER, PUBLIC :: abba_ok = 0, abba_indefinite = 1
81
82 ! the three drivers, for the memory estimate
83 INTEGER, PARAMETER, PRIVATE :: driver_tda = 1, driver_mk = 2, driver_os = 3
84
85 PUBLIC :: bse_davidson_tda, bse_davidson_abba_mk, bse_davidson_abba_os, &
86 bse_davidson_refcheck
87
88CONTAINS
89
90! **************************************************************************************************
91!> \brief Block Davidson with thick restart for the NUM_EXC_EN lowest solutions of
92!> sum_jb A_ia,jb X_jb^n = Ω^n X_ia^n. With the orthonormal basis Z_ia,m the reduced matrix is
93!> Ar_mn = sum_ia Z_ia,m (A Z)_ia,n with eigenpairs (θ_k, c_mk), X_ia^k = sum_m Z_ia,m c_mk,
94!> and the basis grows by the corrections t_ia,k = r_ia,k/(d_ia - θ_k) of the residuals
95!> r_ia,k = sum_m (A Z)_ia,m c_mk - θ_k X_ia^k, d_ia being the diagonal chosen by PRECONDITIONER.
96!> Arrays: Ar_mn (the reduced A) in fm_red_A on the process grid, c_mk in coef_ritz, θ_k c_mk in
97!> coef_ritz_theta
98!> \param mv_env slabs, prefactors and transition energies of the matrix-free A
99!> \param bse_env bse_env carries the BSE_ITERAT settings and BSE_DEBUG_PRINT
100!> \param unit_nr output unit, positive on the writing rank only
101!> \param exc_ens Ω^n in Hartree
102!> \param fm_X X_ia^n on the process grid of mv_env, column n
103! **************************************************************************************************
104 SUBROUTINE bse_davidson_tda(mv_env, bse_env, unit_nr, exc_ens, fm_X)
105
106 TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
107 TYPE(bse_env_type), INTENT(IN) :: bse_env
108 INTEGER, INTENT(IN) :: unit_nr
109 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
110 INTENT(OUT) :: exc_ens
111 TYPE(cp_fm_type), INTENT(OUT) :: fm_x
112
113 CHARACTER(LEN=*), PARAMETER :: routinen = 'bse_davidson_tda'
114
115 INTEGER :: block_size, handle, iter, k, m, m_max, &
116 n_act, n_cand, n_dependent, n_kernel, &
117 n_ov, n_req, n_restart, n_want, nt
118 INTEGER, ALLOCATABLE, DIMENSION(:) :: selected
119 LOGICAL :: all_conv, has_prev, stalled
120 LOGICAL, ALLOCATABLE, DIMENSION(:) :: conv, need
121 REAL(kind=dp) :: dev, t_iter
122 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: de, diag, res, theta, theta_prev, &
123 theta_sub
124 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coef_restart, coef_restart_cand, &
125 coef_ritz, coef_ritz_theta, red_block
126 TYPE(cp_blacs_env_type), POINTER :: diag_blacs_env
127 TYPE(cp_fm_struct_type), POINTER :: fm_struct
128 TYPE(cp_fm_type) :: fm_az, fm_red_a, fm_work, fm_z
129 TYPE(mp_para_env_type), POINTER :: para_env
130
131 CALL timeset(routinen, handle)
132 CALL cite_reference(davidson1975)
133 CALL cite_reference(stathopoulos1998)
134 CALL cite_reference(fukaya2014)
135
136 para_env => mv_env%para_env
137 ! every local_data use below takes whole rows per rank: the npe x 1 grid of bse_matvec_create
138 cpassert(mv_env%blacs_env%num_pe(2) == 1)
139 ! the trial vectors live on an npe x 1 grid, which is a bad one to diagonalize on
140 NULLIFY (diag_blacs_env)
141 CALL cp_blacs_env_create(diag_blacs_env, para_env)
142 n_ov = mv_env%n_ov
143 CALL davidson_sizes(bse_env, n_ov, unit_nr, n_want, n_act, block_size)
144
145 ALLOCATE (diag(n_ov))
146 CALL bse_matvec_diagonal(mv_env, bse_env%preconditioner, diag)
147
148 CALL subspace_ceiling(bse_env, mv_env, driver_tda, n_act, block_size, unit_nr, m_max)
149
150 CALL bse_matvec_vector_struct(mv_env, m_max, fm_struct)
151 CALL cp_fm_create(fm_z, fm_struct, name="fm_Z_davidson", set_zero=.true.)
152 CALL cp_fm_create(fm_az, fm_struct, name="fm_AZ_davidson", set_zero=.true.)
153 CALL cp_fm_struct_release(fm_struct)
154 CALL bse_matvec_vector_struct(mv_env, min(2*n_act, m_max), fm_struct)
155 CALL cp_fm_create(fm_work, fm_struct, name="fm_work_davidson", set_zero=.true.)
156 CALL cp_fm_struct_release(fm_struct)
157 NULLIFY (fm_struct)
158 CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=diag_blacs_env, &
159 nrow_global=m_max, ncol_global=m_max)
160 CALL cp_fm_create(fm_red_a, fm_struct, name="fm_red_A_davidson", set_zero=.true.)
161 CALL cp_fm_struct_release(fm_struct)
162
163 ! c_mk and θ_k c_mk, the restart candidates (current and previous Ritz vectors) and the restart basis Q_pk
164 ALLOCATE (coef_ritz(m_max, n_act), coef_restart_cand(m_max, 2*n_act), coef_ritz_theta(m_max, n_act))
165 ALLOCATE (coef_restart(m_max, 2*n_act))
166 ALLOCATE (theta(m_max), theta_prev(n_act), res(n_act), de(n_act), conv(n_act), need(n_act), &
167 selected(n_act))
168 coef_restart_cand(:, :) = 0.0_dp
169 theta_prev(:) = 0.0_dp
170 has_prev = .false.
171 stalled = .false.
172 n_restart = 0
173 n_dependent = 0
174
175 IF (bse_env%num_guess_transitions > 0) THEN
176 CALL initial_guess_subblock(mv_env, bse_env, diag, n_act, m_max, deg_thresh, fm_z, m, theta_sub, &
177 unit_nr)
178 ELSE
179 CALL initial_guess(diag, n_act, m_max, deg_thresh, fm_z, m)
180 END IF
181 ! both guess routines return orthonormal columns by construction; one Gram matrix checks it
182 IF (extension_deviation(fm_z, fm_z, 0, m, para_env) > ortho_tol) THEN
183 cpabort("BSE Davidson: the initial guess is not orthonormal")
184 END IF
185 ! (A Z)_ia,m of the guess and Ar_mn = sum_ia Z_ia,m (A Z)_ia,n
186 CALL bse_matvec_apply(mv_env, fm_z, 1, m, fm_az)
187 n_kernel = m
188 CALL reduced_gram_blocks(fm_z, m, fm_az, 1, m, block_size, para_env, fm_red_a)
189
190 CALL print_iteration_header('Block Davidson iterations within the TDA:', unit_nr)
191
192 all_conv = .false.
193 DO iter = 1, bse_env%max_iter
194 t_iter = m_walltime()
195
196 ! sum_n Ar_mn c_nk = θ_k c_mk, θ ascending
197 CALL solve_reduced(fm_red_a, m, n_act, para_env, diag_blacs_env, theta, coef_ritz)
198
199 ! thick restart onto the current and the previous Ritz vectors, no kernel application; with the
200 ! restart basis Q_pk in coef_restart, Z_ia,k <- sum_p Z_ia,p Q_pk, A Z alike, Ar_kl <- sum_pq Q_pk Ar_pq Q_ql
201 IF (m + block_size > m_max .AND. m > n_act .AND. m < n_ov) THEN
202 coef_restart_cand(1:m, 1:n_act) = coef_ritz(1:m, 1:n_act)
203 n_cand = n_act
204 IF (has_prev) n_cand = 2*n_act
205 CALL restart_basis(coef_restart_cand, m, n_cand, 2*n_act, coef_restart, k)
206 CALL rotate_in_place(fm_z, m, coef_restart, k, fm_work, fm_az)
207 CALL reduced_rotate(fm_red_a, m, coef_restart, k, para_env, diag_blacs_env)
208 m = k
209 n_restart = n_restart + 1
210 CALL solve_reduced(fm_red_a, m, n_act, para_env, diag_blacs_env, theta, coef_ritz)
211 END IF
212
213 ! r_ia,k = sum_m (A Z)_ia,m c_mk - θ_k sum_m Z_ia,m c_mk for all tracked roots k
214 DO k = 1, n_act
215 coef_ritz_theta(1:m, k) = coef_ritz(1:m, k)*theta(k)
216 END DO
217 CALL subspace_rotate(fm_az, m, coef_ritz(1:m, 1:n_act), n_act, fm_work, 1.0_dp, 0.0_dp)
218 CALL subspace_rotate(fm_z, m, coef_ritz_theta(1:m, 1:n_act), n_act, fm_work, -1.0_dp, 1.0_dp)
219 CALL column_norms(fm_work, 1, n_act, para_env, res)
220
221 CALL convergence_step(bse_env, theta, theta_prev, res, m, n_ov, n_want, iter, t_iter, unit_nr, &
222 de, conv, n_req, need)
223 all_conv = .NOT. any(need)
224 IF (all_conv) EXIT
225
226 IF (iter == bse_env%max_iter) EXIT
227
228 ! corrections for the lowest unconverged roots, t_ia,k = r_ia,k/(d_ia - θ_k)
229 CALL select_roots(bse_env, need, res, min(block_size, m_max - m), selected, nt)
230 CALL davidson_corrections(fm_work, selected(1:nt), theta, diag, fm_z, m + 1)
231 CALL extend_orthonormal(fm_z, m, nt, para_env, n_dependent)
232 ! no new vector: with an unchanged basis the next iteration finds dE = 0, which lets the roots that
233 ! only wait for the energy criterion converge; a second empty pass means the residuals cannot improve
234 IF (nt == 0 .AND. stalled) THEN
235 cpabort("BSE Davidson: correction vectors vanished before convergence")
236 END IF
237 stalled = nt == 0
238 ! post-condition of the extension, one Gram matrix; the residual norms below rely on it
239 IF (extension_deviation(fm_z, fm_z, m, nt, para_env) > ortho_tol) THEN
240 cpabort("BSE Davidson: the basis lost orthonormality")
241 END IF
242
243 IF (nt > 0) THEN
244 CALL bse_matvec_apply(mv_env, fm_z, m + 1, nt, fm_az)
245 n_kernel = n_kernel + nt
246
247 ! new columns Ar_mn = sum_ia Z_ia,m (A Z)_ia,n of the symmetric reduced matrix
248 ALLOCATE (red_block(m + nt, nt))
249 CALL subspace_gram(fm_z, 1, m + nt, fm_az, m + 1, nt, para_env, red_block)
250 CALL reduced_extend(fm_red_a, red_block, m, nt)
251 DEALLOCATE (red_block)
252 END IF
253
254 ! the Ritz vectors of this iteration are the second block of restart candidates
255 coef_restart_cand(:, n_act + 1:2*n_act) = 0.0_dp
256 coef_restart_cand(1:m, n_act + 1:2*n_act) = coef_ritz(1:m, 1:n_act)
257 has_prev = .true.
258 theta_prev(:) = theta(1:n_act)
259 m = m + nt
260 END DO
261
262 IF (.NOT. all_conv) CALL abort_unconverged(theta, res, conv, need, n_req, unit_nr)
263 IF (bse_env%bse_debug_print) CALL print_tracked_roots(theta, res, n_want, n_act, unit_nr)
264
265 IF (bse_env%bse_debug_print) THEN
266 dev = orthonormality_deviation(fm_z, fm_z, m, block_size, para_env)
267 IF (unit_nr > 0) THEN
268 WRITE (unit_nr, '(T2,A10,T13,A,T63,ES18.6)') 'BSE|DEBUG|', &
269 'Max deviation of the basis from orthonormality', dev
270 END IF
271 END IF
272
273 IF (ALLOCATED(theta_sub)) CALL check_guess_bound(theta_sub, theta, n_req, unit_nr)
274 CALL print_summary(iter, n_kernel, n_restart, n_want, n_req, n_act, n_ov, unit_nr, n_dependent=n_dependent)
275
276 ALLOCATE (exc_ens(n_want))
277 exc_ens(:) = theta(1:n_want)
278 CALL bse_matvec_vector_struct(mv_env, n_want, fm_struct)
279 CALL cp_fm_create(fm_x, fm_struct, name="fm_X_davidson", set_zero=.true.)
280 CALL cp_fm_struct_release(fm_struct)
281 ! X_ia^n = sum_m Z_ia,m c_mn
282 CALL subspace_rotate(fm_z, m, coef_ritz(1:m, 1:n_want), n_want, fm_x, 1.0_dp, 0.0_dp)
283
284 CALL cp_fm_release(fm_z)
285 CALL cp_fm_release(fm_az)
286 CALL cp_fm_release(fm_work)
287 CALL cp_fm_release(fm_red_a)
288 CALL cp_blacs_env_release(diag_blacs_env)
289 DEALLOCATE (coef_ritz, coef_restart_cand, coef_ritz_theta, coef_restart, theta, theta_prev, res, de, conv, &
290 need, selected, diag)
291 IF (ALLOCATED(theta_sub)) DEALLOCATE (theta_sub)
292
293 CALL timestop(handle)
294
295 END SUBROUTINE bse_davidson_tda
296
297! **************************************************************************************************
298!> \brief Block Davidson with thick restart for the NUM_EXC_EN lowest solutions of the full problem,
299!> written as sum_jb [(A+B)(A-B)]_ia,jb x_jb^n = (Ω^n)^2 x_ia^n with x^n = X^n - Y^n, which is
300!> symmetric in the inner product <u,v> = sum_ia,jb u_ia (A-B)_ia,jb v_jb. The basis V_ia,m is
301!> orthonormal in that inner product, W_ia,m = sum_jb (A-B)_ia,jb V_jb,m, and the reduced matrix
302!> Hr_mn = sum_ia,jb W_ia,m (A+B)_ia,jb W_jb,n has the eigenpairs ((Ω^k)^2, c_mk). Then
303!> x_ia^k = (Ω^k)^1/2 sum_m V_ia,m c_mk and y_ia^k = (X^k+Y^k)_ia = (Ω^k)^-1/2 sum_m W_ia,m c_mk
304!> obey sum_ia x_ia^k y_ia^k = 1, the residual is r_ia,k = sum_jb (A+B)_ia,jb y_jb^k - Ω^k x_ia^k
305!> and the correction t_ia,k = r_ia,k/(d_ia^2 - (Ω^k)^2), d_ia chosen by PRECONDITIONER.
306!> Arrays: Hr_mn in fm_red_M on the process grid, c_mk in coef_ritz, (Ω^k)^1/2 c_mk in coef_x,
307!> (Ω^k)^-1/2 c_mk in coef_y, Ω^k (Ω^k)^1/2 c_mk in coef_x_omega
308!> \param mv_env slabs, prefactors and transition energies of the matrix-free A and B
309!> \param bse_env bse_env carries the BSE_ITERAT settings and BSE_DEBUG_PRINT
310!> \param unit_nr output unit, positive on the writing rank only
311!> \param exc_ens Ω^n in Hartree
312!> \param fm_X X_ia^n on the process grid of mv_env, column n; created only for abba_ok
313!> \param fm_Y Y_ia^n, as fm_X
314!> \param abba_status abba_indefinite if A-B is not positive definite on the trial space
315!> \param ab_margin smallest eigenvalue of A-B on the trial space in Hartree
316!> \param fm_X_tda converged TDA vectors as initial guess
317! **************************************************************************************************
318 SUBROUTINE bse_davidson_abba_mk(mv_env, bse_env, unit_nr, exc_ens, fm_X, fm_Y, abba_status, &
319 ab_margin, fm_X_tda)
320
321 TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
322 TYPE(bse_env_type), INTENT(IN) :: bse_env
323 INTEGER, INTENT(IN) :: unit_nr
324 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
325 INTENT(OUT) :: exc_ens
326 TYPE(cp_fm_type), INTENT(OUT) :: fm_x, fm_y
327 INTEGER, INTENT(OUT) :: abba_status
328 REAL(kind=dp), INTENT(OUT) :: ab_margin
329 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_x_tda
330
331 CHARACTER(LEN=*), PARAMETER :: routinen = 'bse_davidson_abba_mk'
332
333 INTEGER :: block_size, guess_col, handle, iter, k, m, m_max, n_act, n_cand, n_dependent, &
334 n_guess_cols, n_kernel, n_ov, n_req, n_restart, n_seed, n_want, nrow_local, nt, nt_guess
335 INTEGER, ALLOCATABLE, DIMENSION(:) :: selected
336 LOGICAL :: all_conv, has_prev, stalled
337 LOGICAL, ALLOCATABLE, DIMENSION(:) :: conv, need
338 REAL(kind=dp) :: dev, t_iter
339 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: de, diag, diag_sq, omega, omega_prev, &
340 omega_sq, res, theta, theta_sub
341 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coef_restart, coef_restart_cand, &
342 coef_ritz, coef_x, coef_x_omega, &
343 coef_y, red_block
344 TYPE(cp_blacs_env_type), POINTER :: diag_blacs_env
345 TYPE(cp_fm_struct_type), POINTER :: fm_struct
346 TYPE(cp_fm_type) :: fm_metric, fm_metric_vec, fm_mw, &
347 fm_red_m, fm_scratch, fm_v, fm_w, &
348 fm_work
349 TYPE(mp_para_env_type), POINTER :: para_env
350
351 CALL timeset(routinen, handle)
352 CALL cite_reference(davidson1975)
353 CALL cite_reference(vecharynski2017)
354 CALL cite_reference(bai2012)
355 CALL cite_reference(stathopoulos1998)
356 CALL cite_reference(fukaya2014)
357
358 abba_status = abba_ok
359 ab_margin = 0.0_dp
360 para_env => mv_env%para_env
361 ! every local_data use below takes whole rows per rank: the npe x 1 grid of bse_matvec_create
362 cpassert(mv_env%blacs_env%num_pe(2) == 1)
363 ! the trial vectors live on an npe x 1 grid, which is a bad one to diagonalize on
364 NULLIFY (diag_blacs_env)
365 CALL cp_blacs_env_create(diag_blacs_env, para_env)
366 n_ov = mv_env%n_ov
367 CALL davidson_sizes(bse_env, n_ov, unit_nr, n_want, n_act, block_size)
368
369 ALLOCATE (diag(n_ov))
370 CALL bse_matvec_diagonal(mv_env, bse_env%preconditioner, diag)
371
372 CALL subspace_ceiling(bse_env, mv_env, driver_mk, n_act, block_size, unit_nr, m_max)
373
374 CALL bse_matvec_vector_struct(mv_env, m_max, fm_struct)
375 CALL cp_fm_create(fm_v, fm_struct, name="fm_V_davidson", set_zero=.true.)
376 CALL cp_fm_create(fm_w, fm_struct, name="fm_W_davidson", set_zero=.true.)
377 CALL cp_fm_create(fm_mw, fm_struct, name="fm_MW_davidson", set_zero=.true.)
378 CALL cp_fm_struct_release(fm_struct)
379 CALL bse_matvec_vector_struct(mv_env, min(2*n_act, m_max), fm_struct)
380 CALL cp_fm_create(fm_work, fm_struct, name="fm_work_davidson", set_zero=.true.)
381 CALL cp_fm_struct_release(fm_struct)
382 CALL bse_matvec_vector_struct(mv_env, block_size, fm_struct)
383 CALL cp_fm_create(fm_scratch, fm_struct, name="fm_scratch_davidson", set_zero=.true.)
384 CALL cp_fm_struct_release(fm_struct)
385 NULLIFY (fm_struct)
386 CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=diag_blacs_env, &
387 nrow_global=m_max, ncol_global=m_max)
388 CALL cp_fm_create(fm_red_m, fm_struct, name="fm_red_M_davidson", set_zero=.true.)
389 CALL cp_fm_struct_release(fm_struct)
390 CALL cp_fm_get_info(fm_v, nrow_local=nrow_local)
391
392 ! c_mk, the restart candidates and basis, then (Ω^k)^1/2 c_mk, (Ω^k)^-1/2 c_mk and Ω^k (Ω^k)^1/2 c_mk
393 ALLOCATE (coef_ritz(m_max, n_act), coef_restart_cand(m_max, 2*n_act), coef_restart(m_max, 2*n_act))
394 ALLOCATE (coef_x(m_max, n_act), coef_y(m_max, n_act), coef_x_omega(m_max, n_act))
395 ALLOCATE (theta(m_max), omega(n_act), omega_prev(n_act), res(n_act), de(n_act), conv(n_act), &
396 need(n_act), selected(n_act))
397 coef_restart_cand(:, :) = 0.0_dp
398 omega_prev(:) = 0.0_dp
399 has_prev = .false.
400 stalled = .false.
401 all_conv = .false.
402 n_restart = 0
403 n_kernel = 0
404 n_dependent = 0
405
406 ! guess: the TDA vectors if given, then the eigenvectors of the exact guess block or unit
407 ! vectors on the lowest entries of the diagonal
408 n_seed = 0
409 IF (PRESENT(fm_x_tda)) THEN
410 CALL cp_fm_get_info(fm_x_tda, ncol_global=n_seed)
411 n_seed = min(n_seed, n_act)
412 CALL cp_fm_to_fm(fm_x_tda, fm_v, n_seed)
413 END IF
414 IF (bse_env%num_guess_transitions > 0) THEN
415 CALL initial_guess_subblock(mv_env, bse_env, diag, n_act, m_max, deg_thresh, fm_v, n_guess_cols, theta_sub, &
416 unit_nr, n_seed, abba_status, ab_margin)
417 ELSE
418 CALL initial_guess(diag, n_act, m_max, deg_thresh, fm_v, n_guess_cols, n_seed)
419 END IF
420 ! the corrections divide by d_ia^2 - (Ω^k)^2
421 ALLOCATE (diag_sq(n_ov))
422 diag_sq(:) = diag(:)**2
423
424 ! the guess enters the basis in blocks of the scratch width
425 m = 0
426 guess_col = 1
427 DO WHILE (guess_col <= n_guess_cols)
428 nt_guess = min(block_size, n_guess_cols - guess_col + 1)
429 nt = nt_guess
430 IF (guess_col > m + 1) THEN
431 CALL cp_fm_to_fm(fm_v, fm_v, nt, guess_col, m + 1)
432 END IF
433 CALL extend_basis_mk(mv_env, fm_v, fm_w, fm_mw, fm_scratch, m, nt, n_kernel, n_dependent, abba_status, &
434 ab_margin)
435 IF (abba_status /= abba_ok) EXIT
436 guess_col = guess_col + nt_guess
437 m = m + nt
438 END DO
439
440 IF (abba_status == abba_ok) THEN
441 fm_v%local_data(1:nrow_local, m + 1:n_guess_cols) = 0.0_dp
442 ! Hr_mn = sum_ia W_ia,m (M W)_ia,n on the guess basis
443 CALL reduced_gram_blocks(fm_w, m, fm_mw, 1, m, block_size, para_env, fm_red_m)
444
445 CALL print_iteration_header('Block Davidson iterations for ABBA, (A+B)(A-B) x = E^2 x:', unit_nr)
446
447 DO iter = 1, bse_env%max_iter
448 t_iter = m_walltime()
449
450 ! sum_n Hr_mn c_nk = θ_k c_mk with θ_k = (Ω^k)^2, ascending
451 CALL solve_reduced(fm_red_m, m, n_act, para_env, diag_blacs_env, theta, coef_ritz)
452
453 ! thick restart: V, W and MW are rotated alike by Q_pk in coef_restart, which keeps
454 ! sum_ia V_ia,m W_ia,n = δ_mn, and Hr_kl <- sum_pq Q_pk Hr_pq Q_ql
455 IF (m + block_size > m_max .AND. m > n_act .AND. m < n_ov) THEN
456 coef_restart_cand(1:m, 1:n_act) = coef_ritz(1:m, 1:n_act)
457 n_cand = n_act
458 IF (has_prev) n_cand = 2*n_act
459 CALL restart_basis(coef_restart_cand, m, n_cand, 2*n_act, coef_restart, k)
460 CALL rotate_in_place(fm_v, m, coef_restart, k, fm_work, fm_w, fm_mw)
461 CALL reduced_rotate(fm_red_m, m, coef_restart, k, para_env, diag_blacs_env)
462 m = k
463 n_restart = n_restart + 1
464 CALL solve_reduced(fm_red_m, m, n_act, para_env, diag_blacs_env, theta, coef_ritz)
465 END IF
466
467 ! a Ritz value bounds an eigenvalue of (A+B)(A-B) from above
468 IF (theta(1) <= 0.0_dp) THEN
469 CALL cp_abort(__location__, &
470 "BSE Davidson (MK_DAVIDSON): negative squared excitation energy. Matrix "// &
471 "(A+B) is not positive definite, or the basis lost its orthonormality.")
472 END IF
473 omega(:) = sqrt(theta(1:n_act))
474
475 ! r_ia,k = sum_m (M W)_ia,m c_mk (Ω^k)^-1/2 - Ω^k (Ω^k)^1/2 sum_m V_ia,m c_mk
476 DO k = 1, n_act
477 coef_y(1:m, k) = coef_ritz(1:m, k)/sqrt(omega(k))
478 coef_x_omega(1:m, k) = coef_ritz(1:m, k)*omega(k)*sqrt(omega(k))
479 END DO
480 CALL subspace_rotate(fm_mw, m, coef_y(1:m, 1:n_act), n_act, fm_work, 1.0_dp, 0.0_dp)
481 CALL subspace_rotate(fm_v, m, coef_x_omega(1:m, 1:n_act), n_act, fm_work, -1.0_dp, 1.0_dp)
482 CALL column_norms(fm_work, 1, n_act, para_env, res)
483
484 CALL convergence_step(bse_env, omega, omega_prev, res, m, n_ov, n_want, iter, t_iter, unit_nr, &
485 de, conv, n_req, need)
486 all_conv = .NOT. any(need)
487 IF (all_conv) EXIT
488
489 IF (iter == bse_env%max_iter) EXIT
490
491 ! corrections t_ia,k = r_ia,k/(d_ia^2 - (Ω^k)^2) for the lowest unconverged roots
492 CALL select_roots(bse_env, need, res, min(block_size, m_max - m), selected, nt)
493 CALL davidson_corrections(fm_work, selected(1:nt), theta, diag_sq, fm_v, m + 1)
494 CALL extend_basis_mk(mv_env, fm_v, fm_w, fm_mw, fm_scratch, m, nt, n_kernel, n_dependent, abba_status, &
495 ab_margin)
496 IF (abba_status /= abba_ok) EXIT
497 ! no new vector: with an unchanged basis the next iteration finds dE = 0, which lets the roots that
498 ! only wait for the energy criterion converge; a second empty pass means the residuals cannot improve
499 IF (nt == 0 .AND. stalled) THEN
500 cpabort("BSE Davidson: correction vectors vanished before convergence")
501 END IF
502 stalled = nt == 0
503
504 IF (nt > 0) THEN
505 ! new columns Hr_mn = sum_ia W_ia,m (M W)_ia,n of the symmetric reduced matrix
506 ALLOCATE (red_block(m + nt, nt))
507 CALL subspace_gram(fm_w, 1, m + nt, fm_mw, m + 1, nt, para_env, red_block)
508 CALL reduced_extend(fm_red_m, red_block, m, nt)
509 DEALLOCATE (red_block)
510 END IF
511
512 ! the Ritz vectors of this iteration are the second block of restart candidates
513 coef_restart_cand(:, n_act + 1:2*n_act) = 0.0_dp
514 coef_restart_cand(1:m, n_act + 1:2*n_act) = coef_ritz(1:m, 1:n_act)
515 has_prev = .true.
516 omega_prev(:) = omega(:)
517 m = m + nt
518 END DO
519 END IF
520
521 IF (abba_status == abba_ok) THEN
522 IF (.NOT. all_conv) CALL abort_unconverged(omega, res, conv, need, n_req, unit_nr)
523 IF (bse_env%bse_debug_print) CALL print_tracked_roots(omega, res, n_want, n_act, unit_nr)
524
525 IF (bse_env%bse_debug_print) THEN
526 dev = orthonormality_deviation(fm_v, fm_w, m, block_size, para_env)
527 IF (unit_nr > 0) THEN
528 WRITE (unit_nr, '(T2,A10,T13,A,T63,ES18.6)') 'BSE|DEBUG|', &
529 'Max deviation of the basis from orthonormality', dev
530 END IF
531 END IF
532
533 ! with sum_ia V_ia,m W_ia,n = δ_mn the smallest eigenvalue of A-B on the span of V is the
534 ! inverse of the largest eigenvalue of the metric sum_ia V_ia,m V_ia,n
535 NULLIFY (fm_struct)
536 CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=diag_blacs_env, &
537 nrow_global=m, ncol_global=m)
538 CALL cp_fm_create(fm_metric, fm_struct, name="bse_basis_metric")
539 CALL cp_fm_create(fm_metric_vec, fm_struct, name="bse_basis_metric_vectors")
540 CALL cp_fm_struct_release(fm_struct)
541 CALL reduced_gram_blocks(fm_v, m, fm_v, 1, m, block_size, para_env, fm_metric)
542 ALLOCATE (omega_sq(m))
543 omega_sq(:) = 0.0_dp
544 CALL choose_eigv_solver(fm_metric, fm_metric_vec, omega_sq)
545 ab_margin = 1.0_dp/omega_sq(m)
546 DEALLOCATE (omega_sq)
547 CALL cp_fm_release(fm_metric)
548 CALL cp_fm_release(fm_metric_vec)
549
550 IF (ALLOCATED(theta_sub)) CALL check_guess_bound(theta_sub, omega, n_req, unit_nr)
551 CALL print_summary(iter, n_kernel, n_restart, n_want, n_req, n_act, n_ov, unit_nr, ab_margin, n_dependent)
552
553 ! x_ia^n = (Ω^n)^1/2 sum_m V_ia,m c_mn, y_ia^n = (Ω^n)^-1/2 sum_m W_ia,m c_mn,
554 ! X_ia^n = (y_ia^n + x_ia^n)/2, Y_ia^n = (y_ia^n - x_ia^n)/2
555 ALLOCATE (exc_ens(n_want))
556 exc_ens(:) = omega(1:n_want)
557 DO k = 1, n_want
558 coef_x(1:m, k) = coef_ritz(1:m, k)*sqrt(omega(k))
559 coef_y(1:m, k) = coef_ritz(1:m, k)/sqrt(omega(k))
560 END DO
561 CALL bse_matvec_vector_struct(mv_env, n_want, fm_struct)
562 CALL cp_fm_create(fm_x, fm_struct, name="fm_X_davidson", set_zero=.true.)
563 CALL cp_fm_create(fm_y, fm_struct, name="fm_Y_davidson", set_zero=.true.)
564 CALL cp_fm_struct_release(fm_struct)
565 CALL subspace_rotate(fm_w, m, coef_y(1:m, 1:n_want), n_want, fm_x, 0.5_dp, 0.0_dp)
566 CALL subspace_rotate(fm_v, m, coef_x(1:m, 1:n_want), n_want, fm_x, 0.5_dp, 1.0_dp)
567 CALL subspace_rotate(fm_w, m, coef_y(1:m, 1:n_want), n_want, fm_y, 0.5_dp, 0.0_dp)
568 CALL subspace_rotate(fm_v, m, coef_x(1:m, 1:n_want), n_want, fm_y, -0.5_dp, 1.0_dp)
569 END IF
570
571 CALL cp_fm_release(fm_v)
572 CALL cp_fm_release(fm_w)
573 CALL cp_fm_release(fm_mw)
574 CALL cp_fm_release(fm_work)
575 CALL cp_fm_release(fm_scratch)
576 CALL cp_fm_release(fm_red_m)
577 CALL cp_blacs_env_release(diag_blacs_env)
578 DEALLOCATE (coef_ritz, coef_restart_cand, coef_x, coef_y, coef_x_omega, coef_restart, theta, omega, omega_prev, &
579 res, de, conv, need, selected, diag, diag_sq)
580 IF (ALLOCATED(theta_sub)) DEALLOCATE (theta_sub)
581
582 CALL timestop(handle)
583
584 END SUBROUTINE bse_davidson_abba_mk
585
586! **************************************************************************************************
587!> \brief Block Davidson with thick restart for the NUM_EXC_EN lowest solutions of the full problem
588!> in the paired form of Olsen and Stratmann: with x^n = X^n - Y^n, y^n = X^n + Y^n, M = A+B and
589!> K = A-B, sum_jb M_ia,jb y_jb^n = Ω^n x_ia^n and sum_jb K_ia,jb x_jb^n = Ω^n y_ia^n. One
590!> orthonormal basis b_ia,m carries both vectors: Mr_mn = sum_ia b_ia,m (M b)_ia,n and
591!> Kr_mn = sum_ia b_ia,m (K b)_ia,n are solved by solve_reduced_paired for Ω^k, R_mk and L_mk with
592!> y_ia^k = sum_m b_ia,m R_mk, x_ia^k = sum_m b_ia,m L_mk and sum_ia x_ia^k y_ia^k = 1. The basis
593!> grows by the corrections t_ia,k = r_ia,k/(d_ia - Ω^k) of both residuals,
594!> rR_ia,k = sum_m (M b)_ia,m R_mk - Ω^k x_ia^k and rL_ia,k = sum_m (K b)_ia,m L_mk - Ω^k y_ia^k,
595!> d_ia chosen by PRECONDITIONER; the rR block enters before the rL block, which is nearly
596!> parallel to it for a small B. Arrays: Mr_mn and Kr_mn in fm_red_M and fm_red_K on the process
597!> grid, R_mk in coef_right, L_mk in coef_left
598!> \param mv_env slabs, prefactors and transition energies of the matrix-free A and B
599!> \param bse_env bse_env carries the BSE_ITERAT settings and BSE_DEBUG_PRINT
600!> \param unit_nr output unit, positive on the writing rank only
601!> \param exc_ens Ω^n in Hartree
602!> \param fm_X X_ia^n on the process grid of mv_env, column n; created only for abba_ok
603!> \param fm_Y Y_ia^n, as fm_X
604!> \param abba_status abba_indefinite if A-B is not positive definite on the trial space
605!> \param ab_margin smallest eigenvalue of A-B on the trial space in Hartree
606!> \param fm_X_tda converged TDA vectors as initial guess, x = y = X
607! **************************************************************************************************
608 SUBROUTINE bse_davidson_abba_os(mv_env, bse_env, unit_nr, exc_ens, fm_X, fm_Y, abba_status, &
609 ab_margin, fm_X_tda)
610
611 TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
612 TYPE(bse_env_type), INTENT(IN) :: bse_env
613 INTEGER, INTENT(IN) :: unit_nr
614 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
615 INTENT(OUT) :: exc_ens
616 TYPE(cp_fm_type), INTENT(OUT) :: fm_x, fm_y
617 INTEGER, INTENT(OUT) :: abba_status
618 REAL(kind=dp), INTENT(OUT) :: ab_margin
619 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_x_tda
620
621 CHARACTER(LEN=*), PARAMETER :: routinen = 'bse_davidson_abba_os'
622
623 INTEGER :: block_size, guess_col, handle, iter, k, m, m_max, n_act, n_cand, n_dependent, &
624 n_guess_cols, n_kernel, n_ov, n_pair, n_req, n_restart, n_seed, n_want, nrow_local, nt, &
625 nt_guess, nt_l, nt_r, nt_sel
626 INTEGER, ALLOCATABLE, DIMENSION(:) :: selected
627 LOGICAL :: all_conv, has_prev, stalled
628 LOGICAL, ALLOCATABLE, DIMENSION(:) :: conv, need
629 REAL(kind=dp) :: dev, kr_min, t_iter
630 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: de, diag, omega, omega_prev, res, res_l, &
631 theta_sub
632 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coef_left, coef_restart, &
633 coef_restart_cand, coef_right, &
634 coef_work, red_block_k, red_block_m
635 TYPE(cp_blacs_env_type), POINTER :: diag_blacs_env
636 TYPE(cp_fm_struct_type), POINTER :: fm_struct
637 TYPE(cp_fm_type) :: fm_b, fm_kb, fm_mb, fm_red_k, fm_red_m, &
638 fm_scratch, fm_work
639 TYPE(mp_para_env_type), POINTER :: para_env
640
641 CALL timeset(routinen, handle)
642 CALL cite_reference(davidson1975)
643 CALL cite_reference(olsen1988)
644 CALL cite_reference(stratmann1998)
645 CALL cite_reference(bai2012)
646 CALL cite_reference(stathopoulos1998)
647 CALL cite_reference(fukaya2014)
648
649 abba_status = abba_ok
650 ab_margin = 0.0_dp
651 para_env => mv_env%para_env
652 ! every local_data use below takes whole rows per rank: the npe x 1 grid of bse_matvec_create
653 cpassert(mv_env%blacs_env%num_pe(2) == 1)
654 ! the trial vectors live on an npe x 1 grid, which is a bad one to diagonalize on
655 NULLIFY (diag_blacs_env)
656 CALL cp_blacs_env_create(diag_blacs_env, para_env)
657 n_ov = mv_env%n_ov
658 CALL davidson_sizes(bse_env, n_ov, unit_nr, n_want, n_act, block_size)
659
660 ALLOCATE (diag(n_ov))
661 CALL bse_matvec_diagonal(mv_env, bse_env%preconditioner, diag)
662
663 CALL subspace_ceiling(bse_env, mv_env, driver_os, n_act, block_size, unit_nr, m_max)
664
665 CALL bse_matvec_vector_struct(mv_env, m_max, fm_struct)
666 CALL cp_fm_create(fm_b, fm_struct, name="fm_b_davidson", set_zero=.true.)
667 CALL cp_fm_create(fm_mb, fm_struct, name="fm_Mb_davidson", set_zero=.true.)
668 CALL cp_fm_create(fm_kb, fm_struct, name="fm_Kb_davidson", set_zero=.true.)
669 CALL cp_fm_struct_release(fm_struct)
670 ! both residual blocks, or the restart basis
671 CALL bse_matvec_vector_struct(mv_env, max(2*n_act, min(4*n_act, m_max)), fm_struct)
672 CALL cp_fm_create(fm_work, fm_struct, name="fm_work_davidson", set_zero=.true.)
673 CALL cp_fm_struct_release(fm_struct)
674 CALL bse_matvec_vector_struct(mv_env, 2*block_size, fm_struct)
675 CALL cp_fm_create(fm_scratch, fm_struct, name="fm_scratch_davidson", set_zero=.true.)
676 CALL cp_fm_struct_release(fm_struct)
677 NULLIFY (fm_struct)
678 CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=diag_blacs_env, &
679 nrow_global=m_max, ncol_global=m_max)
680 CALL cp_fm_create(fm_red_m, fm_struct, name="fm_red_M_davidson", set_zero=.true.)
681 CALL cp_fm_create(fm_red_k, fm_struct, name="fm_red_K_davidson", set_zero=.true.)
682 CALL cp_fm_struct_release(fm_struct)
683 CALL cp_fm_get_info(fm_b, nrow_local=nrow_local)
684
685 ! R_mk, L_mk and a work block, then the restart candidates (current and previous pairs) and the basis Q_pk
686 ALLOCATE (coef_right(m_max, n_act), coef_left(m_max, n_act), coef_work(m_max, n_act))
687 ALLOCATE (coef_restart_cand(m_max, 4*n_act), coef_restart(m_max, 4*n_act))
688 ALLOCATE (omega(n_act), omega_prev(n_act), res(n_act), res_l(n_act), de(n_act), conv(n_act), &
689 need(n_act), selected(n_act))
690 coef_restart_cand(:, :) = 0.0_dp
691 omega_prev(:) = 0.0_dp
692 kr_min = 0.0_dp
693 has_prev = .false.
694 stalled = .false.
695 all_conv = .false.
696 n_restart = 0
697 n_kernel = 0
698 n_dependent = 0
699
700 ! guess: the TDA vectors if given (each a pair by itself), then the pairs x_k, y_k of the exact
701 ! guess block or unit vectors on the lowest entries of the diagonal
702 n_seed = 0
703 IF (PRESENT(fm_x_tda)) THEN
704 CALL cp_fm_get_info(fm_x_tda, ncol_global=n_seed)
705 n_seed = min(n_seed, n_act)
706 CALL cp_fm_to_fm(fm_x_tda, fm_b, n_seed)
707 END IF
708 n_pair = 0
709 IF (bse_env%num_guess_transitions > 0) THEN
710 CALL initial_guess_subblock(mv_env, bse_env, diag, n_act, m_max, deg_thresh, fm_b, n_guess_cols, theta_sub, &
711 unit_nr, n_seed, abba_status, ab_margin, n_pair)
712 ELSE
713 CALL initial_guess(diag, n_act, m_max, deg_thresh, fm_b, n_guess_cols, n_seed)
714 END IF
715 n_guess_cols = n_guess_cols + n_pair
716
717 ! the guess enters the basis in blocks of the scratch width
718 m = 0
719 guess_col = 1
720 DO WHILE (guess_col <= n_guess_cols)
721 nt_guess = min(2*block_size, n_guess_cols - guess_col + 1)
722 nt = nt_guess
723 IF (guess_col > m + 1) THEN
724 CALL cp_fm_to_fm(fm_b, fm_b, nt, guess_col, m + 1)
725 END IF
726 CALL extend_basis_os(mv_env, fm_b, fm_mb, fm_kb, fm_scratch, m, nt, n_kernel, n_dependent)
727 guess_col = guess_col + nt_guess
728 m = m + nt
729 END DO
730
731 IF (abba_status == abba_ok) THEN
732 fm_b%local_data(1:nrow_local, m + 1:n_guess_cols) = 0.0_dp
733 ! Mr_mn = sum_ia b_ia,m (M b)_ia,n and Kr_mn = sum_ia b_ia,m (K b)_ia,n on the guess basis
734 CALL reduced_gram_blocks(fm_b, m, fm_mb, 1, m, block_size, para_env, fm_red_m)
735 CALL reduced_gram_blocks(fm_b, m, fm_kb, 1, m, block_size, para_env, fm_red_k)
736
737 CALL print_iteration_header('Block Davidson iterations for ABBA, Olsen-Stratmann paired subspace:', &
738 unit_nr)
739
740 DO iter = 1, bse_env%max_iter
741 t_iter = m_walltime()
742
743 ! Ω^k ascending with the coefficients R_mk of y^k and L_mk of x^k, sum_m R_mk L_mk = 1
744 CALL solve_reduced_paired(fm_red_m, fm_red_k, m, n_act, para_env, omega, coef_right, coef_left, kr_min, &
745 abba_status, diag_blacs_env)
746 ! the smallest eigenvalue of A-B on the span of b is that of the last Kr
747 ab_margin = kr_min
748 IF (abba_status /= abba_ok) EXIT
749
750 ! thick restart: b, Mb and Kb are rotated alike; the current pairs come first, so that
751 ! they survive the cap that leaves room for one block of corrections; with Q_pk in
752 ! coef_restart, b_ia,k <- sum_p b_ia,p Q_pk and Mr_kl <- sum_pq Q_pk Mr_pq Q_ql, Kr alike
753 IF (m + 2*block_size > m_max .AND. m > 2*n_act .AND. m < n_ov) THEN
754 coef_restart_cand(1:m, 1:n_act) = coef_right(1:m, :)
755 coef_restart_cand(1:m, n_act + 1:2*n_act) = coef_left(1:m, :)
756 n_cand = 2*n_act
757 IF (has_prev) n_cand = 4*n_act
758 CALL restart_basis(coef_restart_cand, m, n_cand, m_max - 2*block_size, coef_restart, k)
759 CALL rotate_in_place(fm_b, m, coef_restart, k, fm_work, fm_mb, fm_kb)
760 CALL reduced_rotate(fm_red_m, m, coef_restart, k, para_env, diag_blacs_env)
761 CALL reduced_rotate(fm_red_k, m, coef_restart, k, para_env, diag_blacs_env)
762 m = k
763 n_restart = n_restart + 1
764 CALL solve_reduced_paired(fm_red_m, fm_red_k, m, n_act, para_env, omega, coef_right, coef_left, kr_min, &
765 abba_status, diag_blacs_env)
766 ab_margin = kr_min
767 IF (abba_status /= abba_ok) EXIT
768 END IF
769
770 ! rR_ia,k = sum_m (M b)_ia,m R_mk - Ω^k sum_m b_ia,m L_mk into the columns 1..n_act and
771 ! rL_ia,k = sum_m (K b)_ia,m L_mk - Ω^k sum_m b_ia,m R_mk into the columns n_act+1..2 n_act
772 DO k = 1, n_act
773 coef_work(1:m, k) = coef_left(1:m, k)*omega(k)
774 END DO
775 CALL subspace_rotate(fm_mb, m, coef_right(1:m, 1:n_act), n_act, fm_work, 1.0_dp, 0.0_dp)
776 CALL subspace_rotate(fm_b, m, coef_work(1:m, 1:n_act), n_act, fm_work, -1.0_dp, 1.0_dp)
777 DO k = 1, n_act
778 coef_work(1:m, k) = coef_right(1:m, k)*omega(k)
779 END DO
780 CALL subspace_rotate(fm_kb, m, coef_left(1:m, 1:n_act), n_act, fm_work, 1.0_dp, 0.0_dp, n_act + 1)
781 CALL subspace_rotate(fm_b, m, coef_work(1:m, 1:n_act), n_act, fm_work, -1.0_dp, 1.0_dp, n_act + 1)
782 CALL column_norms(fm_work, 1, n_act, para_env, res)
783 CALL column_norms(fm_work, n_act + 1, n_act, para_env, res_l)
784 res(:) = max(res(:), res_l(:))
785
786 CALL convergence_step(bse_env, omega, omega_prev, res, m, n_ov, n_want, iter, t_iter, unit_nr, &
787 de, conv, n_req, need)
788 all_conv = .NOT. any(need)
789 IF (all_conv) EXIT
790
791 IF (iter == bse_env%max_iter) EXIT
792
793 ! corrections t_ia,k = r_ia,k/(d_ia - Ω^k) of the lowest unconverged roots: the rR block
794 ! against the basis, then the rL block of the same roots against the extended basis;
795 ! both fit the room, which is odd only when the ceiling is the number of transitions
796 CALL select_roots(bse_env, need, res, min(block_size, (m_max - m + 1)/2), selected, nt_sel)
797 nt_r = nt_sel
798 CALL davidson_corrections(fm_work, selected(1:nt_r), omega, diag, fm_b, m + 1)
799 CALL extend_basis_os(mv_env, fm_b, fm_mb, fm_kb, fm_scratch, m, nt_r, n_kernel, n_dependent)
800 nt_l = min(nt_sel, m_max - m - nt_r)
801 CALL davidson_corrections(fm_work, selected(1:nt_l), omega, diag, fm_b, m + nt_r + 1, n_act)
802 CALL extend_basis_os(mv_env, fm_b, fm_mb, fm_kb, fm_scratch, m + nt_r, nt_l, n_kernel, n_dependent)
803 nt = nt_r + nt_l
804 ! no new vector: with an unchanged basis the next iteration finds dE = 0, which lets the roots that
805 ! only wait for the energy criterion converge; a second empty pass means the residuals cannot improve
806 IF (nt == 0 .AND. stalled) THEN
807 cpabort("BSE Davidson: correction vectors vanished before convergence")
808 END IF
809 stalled = nt == 0
810
811 IF (nt > 0) THEN
812 ! new columns n of Mr_mn = sum_ia b_ia,m (M b)_ia,n and Kr_mn = sum_ia b_ia,m (K b)_ia,n
813 ALLOCATE (red_block_m(m + nt, nt), red_block_k(m + nt, nt))
814 CALL subspace_gram(fm_b, 1, m + nt, fm_mb, m + 1, nt, para_env, red_block_m)
815 CALL subspace_gram(fm_b, 1, m + nt, fm_kb, m + 1, nt, para_env, red_block_k)
816 CALL reduced_extend(fm_red_m, red_block_m, m, nt)
817 CALL reduced_extend(fm_red_k, red_block_k, m, nt)
818 DEALLOCATE (red_block_m, red_block_k)
819 END IF
820
821 ! the pairs of this iteration are the second block of restart candidates
822 coef_restart_cand(:, 2*n_act + 1:4*n_act) = 0.0_dp
823 coef_restart_cand(1:m, 2*n_act + 1:3*n_act) = coef_right(1:m, :)
824 coef_restart_cand(1:m, 3*n_act + 1:4*n_act) = coef_left(1:m, :)
825 has_prev = .true.
826 omega_prev(:) = omega(:)
827 m = m + nt
828 END DO
829 END IF
830
831 IF (abba_status == abba_ok) THEN
832 IF (.NOT. all_conv) CALL abort_unconverged(omega, res, conv, need, n_req, unit_nr)
833 IF (bse_env%bse_debug_print) CALL print_tracked_roots(omega, res, n_want, n_act, unit_nr)
834
835 IF (bse_env%bse_debug_print) THEN
836 dev = orthonormality_deviation(fm_b, fm_b, m, block_size, para_env)
837 IF (unit_nr > 0) THEN
838 WRITE (unit_nr, '(T2,A10,T13,A,T63,ES18.6)') 'BSE|DEBUG|', &
839 'Max deviation of the basis from orthonormality', dev
840 END IF
841 END IF
842
843 IF (ALLOCATED(theta_sub)) CALL check_guess_bound(theta_sub, omega, n_req, unit_nr)
844 CALL print_summary(iter, n_kernel, n_restart, n_want, n_req, n_act, n_ov, unit_nr, ab_margin, n_dependent)
845
846 ! X_ia^n = sum_m b_ia,m (R_mn + L_mn)/2, Y_ia^n = sum_m b_ia,m (R_mn - L_mn)/2
847 ALLOCATE (exc_ens(n_want))
848 exc_ens(:) = omega(1:n_want)
849 CALL bse_matvec_vector_struct(mv_env, n_want, fm_struct)
850 CALL cp_fm_create(fm_x, fm_struct, name="fm_X_davidson", set_zero=.true.)
851 CALL cp_fm_create(fm_y, fm_struct, name="fm_Y_davidson", set_zero=.true.)
852 CALL cp_fm_struct_release(fm_struct)
853 coef_work(1:m, 1:n_want) = 0.5_dp*(coef_right(1:m, 1:n_want) + coef_left(1:m, 1:n_want))
854 CALL subspace_rotate(fm_b, m, coef_work(1:m, 1:n_want), n_want, fm_x, 1.0_dp, 0.0_dp)
855 coef_work(1:m, 1:n_want) = 0.5_dp*(coef_right(1:m, 1:n_want) - coef_left(1:m, 1:n_want))
856 CALL subspace_rotate(fm_b, m, coef_work(1:m, 1:n_want), n_want, fm_y, 1.0_dp, 0.0_dp)
857 END IF
858
859 CALL cp_fm_release(fm_b)
860 CALL cp_fm_release(fm_mb)
861 CALL cp_fm_release(fm_kb)
862 CALL cp_fm_release(fm_work)
863 CALL cp_fm_release(fm_scratch)
864 CALL cp_fm_release(fm_red_m)
865 CALL cp_fm_release(fm_red_k)
866 CALL cp_blacs_env_release(diag_blacs_env)
867 DEALLOCATE (coef_right, coef_left, coef_work, coef_restart_cand, coef_restart, omega, omega_prev, res, res_l, &
868 de, conv, need, selected, diag)
869 IF (ALLOCATED(theta_sub)) DEALLOCATE (theta_sub)
870
871 CALL timestop(handle)
872
873 END SUBROUTINE bse_davidson_abba_os
874
875! **************************************************************************************************
876!> \brief Makes the nt columns T that follow column m of V orthonormal in the (A-B) inner product and
877!> stores their products: T_ia,k <- T_ia,k - sum_m V_ia,m sum_jb W_jb,m T_jb,k, Euclidean
878!> orthonormalisation of T, W_ia,k = sum_jb (A-B)_ia,jb T_jb,k, T_ia,k <- sum_l T_ia,l (R^-1)_lk and
879!> W alike with sum_p R_pk R_pl = sum_ia T_ia,k W_ia,l, and (M W)_ia,k = sum_jb (A+B)_ia,jb W_jb,k
880!> \param mv_env the matrix-free A and B
881!> \param fm_V basis; columns 1..m are kept, m+1..m+nt hold T on entry and the kept columns on exit
882!> \param fm_W (A-B) V, extended in the same columns
883!> \param fm_MW (A+B) W, extended in the same columns
884!> \param fm_scratch scratch of at least nt columns
885!> \param m columns of the basis before the extension
886!> \param nt on exit the number of columns kept
887!> \param n_kernel counter of the kernel applications
888!> \param n_dependent counter of the correction vectors dropped as dependent
889!> \param abba_status abba_indefinite if the Gram matrix of T in the (A-B) inner product is not positive definite
890!> \param ab_margin smallest eigenvalue of T^T (A-B) T if that matrix is not positive definite
891! **************************************************************************************************
892 SUBROUTINE extend_basis_mk(mv_env, fm_V, fm_W, fm_MW, fm_scratch, m, nt, n_kernel, n_dependent, abba_status, &
893 ab_margin)
894
895 TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
896 TYPE(cp_fm_type), INTENT(IN) :: fm_v, fm_w, fm_mw, fm_scratch
897 INTEGER, INTENT(IN) :: m
898 INTEGER, INTENT(INOUT) :: nt, n_kernel, n_dependent, abba_status
899 REAL(kind=dp), INTENT(INOUT) :: ab_margin
900
901 CHARACTER(LEN=*), PARAMETER :: routinen = 'extend_basis_mk'
902
903 INTEGER :: handle
904
905 CALL timeset(routinen, handle)
906
907 ! T_ia,k <- T_ia,k - sum_m V_ia,m sum_jb W_jb,m T_jb,k, then Euclidean orthonormalisation
908 CALL extend_orthonormal(fm_v, m, nt, mv_env%para_env, n_dependent, fm_w)
909
910 IF (nt > 0) THEN
911 ! W_ia,k = sum_jb (A-B)_ia,jb T_jb,k: A T into W, B T into the scratch
912 CALL bse_matvec_apply(mv_env, fm_v, m + 1, nt, fm_w, fm_scratch, 1)
913 CALL columns_axpy(-1.0_dp, fm_scratch, 1, fm_w, m + 1, nt)
914 n_kernel = n_kernel + nt
915
916 ! T_ia,k <- sum_l T_ia,l (R^-1)_lk and W alike, with sum_p R_pk R_pl = sum_ia T_ia,k W_ia,l
917 CALL cholqr2_metric(fm_v, fm_w, m + 1, nt, mv_env%para_env, abba_status, ab_margin)
918
919 IF (abba_status == abba_ok) THEN
920 ! post-condition of the extension, one Gram matrix; the residual norms below rely on it
921 IF (extension_deviation(fm_v, fm_w, m, nt, mv_env%para_env) > ortho_tol) THEN
922 cpabort("BSE Davidson: the basis lost orthonormality in the (A-B) inner product")
923 END IF
924 ! (M W)_ia,k = sum_jb (A+B)_ia,jb W_jb,k: A W into MW, B W into the scratch
925 CALL bse_matvec_apply(mv_env, fm_w, m + 1, nt, fm_mw, fm_scratch, 1)
926 CALL columns_axpy(1.0_dp, fm_scratch, 1, fm_mw, m + 1, nt)
927 n_kernel = n_kernel + nt
928 END IF
929 END IF
930
931 CALL timestop(handle)
932
933 END SUBROUTINE extend_basis_mk
934
935! **************************************************************************************************
936!> \brief Two Cholesky QR passes in the inner product of a symmetric K with the products KT at hand:
937!> G_kl = sum_ia T_ia,k (KT)_ia,l = sum_p R_pk R_pl, T_ia,k <- sum_l T_ia,l (R^-1)_lk, KT alike. A failed
938!> factorization of a G with a negative eigenvalue reports K as indefinite, any other failure
939!> is repeated once with a diagonal shift and a third pass. The Gram is the only communication:
940!> R is factorised on one rank and broadcast, the triangular solve runs on the local rows
941!> \param fm_T the block to orthonormalise, columns first_col .. first_col+nt-1
942!> \param fm_KT K T in the same columns, transformed alike
943!> \param first_col first column of the block
944!> \param nt number of columns in the block
945!> \param para_env communicator of the row-distributed matrices
946!> \param abba_status abba_indefinite if gram has a negative eigenvalue, else untouched
947!> \param eig_min smallest eigenvalue of gram, set if it is negative
948! **************************************************************************************************
949 SUBROUTINE cholqr2_metric(fm_T, fm_KT, first_col, nt, para_env, abba_status, eig_min)
950
951 TYPE(cp_fm_type), INTENT(IN) :: fm_t, fm_kt
952 INTEGER, INTENT(IN) :: first_col, nt
953 TYPE(mp_para_env_type), INTENT(IN) :: para_env
954 INTEGER, INTENT(INOUT) :: abba_status
955 REAL(kind=dp), INTENT(INOUT) :: eig_min
956
957 CHARACTER(LEN=*), PARAMETER :: routinen = 'cholqr2_metric'
958
959 INTEGER :: handle, info, ipass, k, n_ov, n_pass, &
960 nrow_local
961 REAL(kind=dp) :: shift
962 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigval
963 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: chol_factor, gram, gram_raw
964
965 CALL timeset(routinen, handle)
966 CALL cp_fm_get_info(fm_t, nrow_local=nrow_local, nrow_global=n_ov)
967 ALLOCATE (gram_raw(nt, nt), gram(nt, nt), chol_factor(nt, nt), eigval(nt))
968
969 n_pass = 2
970 ipass = 0
971 DO WHILE (ipass < n_pass)
972 ipass = ipass + 1
973 ! G_kl = (sum_ia T_ia,k (KT)_ia,l + sum_ia T_ia,l (KT)_ia,k)/2 = sum_p R_pk R_pl
974 CALL subspace_gram(fm_t, first_col, nt, fm_kt, first_col, nt, para_env, gram_raw)
975 gram(:, :) = 0.5_dp*(gram_raw(:, :) + transpose(gram_raw(:, :)))
976 chol_factor(:, :) = gram(:, :)
977 info = 0
978 IF (para_env%is_source()) CALL dpotrf('U', nt, chol_factor, nt, info)
979 CALL para_env%bcast(info)
980 IF (info /= 0) THEN
981 IF (n_pass == 3) THEN
982 CALL cp_abort(__location__, &
983 "BSE Davidson: orthonormalisation in the (A-B) inner product broke down")
984 END IF
985 chol_factor(:, :) = gram(:, :)
986 eigval(:) = 0.0_dp
987 IF (para_env%is_source()) CALL diamat_all(chol_factor, eigval)
988 CALL para_env%bcast(eigval)
989 IF (eigval(1) < 0.0_dp) THEN
990 abba_status = abba_indefinite
991 eig_min = eigval(1)
992 EXIT
993 END IF
994 ! a positive G that failed numerically: factorise G_kl + shift δ_kl and pass once more, the
995 ! shift of shifted Cholesky QR (Fukaya et al. 2020), 11 (N_ov nt + nt (nt+1)) EPSILON |G|_F
996 n_pass = 3
997 shift = 11.0_dp*(real(n_ov, dp)*real(nt, dp) + real(nt, dp)*real(nt + 1, dp))* &
998 epsilon(1.0_dp)*norm2(gram)
999 chol_factor(:, :) = gram(:, :)
1000 DO k = 1, nt
1001 chol_factor(k, k) = chol_factor(k, k) + shift
1002 END DO
1003 IF (para_env%is_source()) CALL dpotrf('U', nt, chol_factor, nt, info)
1004 CALL para_env%bcast(info)
1005 IF (info /= 0) THEN
1006 CALL cp_abort(__location__, &
1007 "BSE Davidson: orthonormalisation in the (A-B) inner product broke down")
1008 END IF
1009 END IF
1010 ! T_ia,k <- sum_l T_ia,l (R^-1)_lk and KT alike, on the local rows
1011 CALL para_env%bcast(chol_factor)
1012 IF (nrow_local > 0) THEN
1013 CALL dtrsm('R', 'U', 'N', 'N', nrow_local, nt, 1.0_dp, chol_factor, nt, &
1014 fm_t%local_data(:, first_col:first_col + nt - 1), SIZE(fm_t%local_data, 1))
1015 CALL dtrsm('R', 'U', 'N', 'N', nrow_local, nt, 1.0_dp, chol_factor, nt, &
1016 fm_kt%local_data(:, first_col:first_col + nt - 1), SIZE(fm_kt%local_data, 1))
1017 END IF
1018 END DO
1019
1020 DEALLOCATE (gram_raw, gram, chol_factor, eigval)
1021 CALL timestop(handle)
1022
1023 END SUBROUTINE cholqr2_metric
1024
1025! **************************************************************************************************
1026!> \brief Paired reduced problem on the orthonormal basis b_ia,m,
1027!> Mr_mn = sum_ia b_ia,m (M b)_ia,n and Kr_mn = sum_ia b_ia,m (K b)_ia,n,
1028!> sum_n Mr_mn R_nk = Ω_k L_mk and sum_n Kr_mn L_nk = Ω_k R_mk,
1029!> solved as in initial_guess_subblock through the symmetric form:
1030!> Kr_mn = sum_p U_mp κ_p U_np and S_mn = sum_p U_mp √κ_p U_np, the square root of Kr,
1031!> sum_nop S_mn Mr_no S_op T_pk = λ_k T_mk with Ω_k = √λ_k,
1032!> R_mk = sum_n S_mn T_nk/√Ω_k and L_mk = sum_n Mr_mn R_nk/Ω_k.
1033!> κ_1 <= 0 reports K as indefinite on the span of b.
1034!> Every m x m step runs on the grid of blacs_env: the two eigensolves through
1035!> choose_eigv_solver and the four products through parallel_gemm; only the pair
1036!> coefficients come back replicated, R_mk in coef_right and L_mk in coef_left
1037!> \param fm_red_M Mr in its live block, on the grid of blacs_env
1038!> \param fm_red_K Kr in its live block, on the grid of blacs_env
1039!> \param m dimension of the live block
1040!> \param n_act number of pairs formed
1041!> \param para_env communicator of the grid
1042!> \param omega Ω_k in Hartree
1043!> \param coef_right coefficients of y_k in column k, zero below row m
1044!> \param coef_left coefficients of x_k in column k, zero below row m
1045!> \param kr_min κ_1, the smallest eigenvalue of Kr
1046!> \param abba_status abba_indefinite when κ_1 <= 0; omega, coef_right and coef_left are then untouched
1047!> \param blacs_env grid the reduced problem is solved on
1048! **************************************************************************************************
1049 SUBROUTINE solve_reduced_paired(fm_red_M, fm_red_K, m, n_act, para_env, omega, coef_right, coef_left, kr_min, &
1050 abba_status, blacs_env)
1051
1052 TYPE(cp_fm_type), INTENT(IN) :: fm_red_m, fm_red_k
1053 INTEGER, INTENT(IN) :: m, n_act
1054 TYPE(mp_para_env_type), INTENT(IN), TARGET :: para_env
1055 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: omega
1056 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: coef_right, coef_left
1057 REAL(kind=dp), INTENT(OUT) :: kr_min
1058 INTEGER, INTENT(INOUT) :: abba_status
1059 TYPE(cp_blacs_env_type), INTENT(IN), TARGET :: blacs_env
1060
1061 CHARACTER(LEN=*), PARAMETER :: routinen = 'solve_reduced_paired'
1062
1063 INTEGER :: handle
1064 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: kappa, lambda
1065 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1066 TYPE(cp_fm_type) :: fm_h, fm_h_eigvec, fm_k_eigvec, &
1067 fm_k_sqrt, fm_l, fm_mr, fm_r, fm_work
1068
1069 CALL timeset(routinen, handle)
1070 ! n_act pairs are extracted from an m-dimensional basis
1071 cpassert(n_act <= m)
1072
1073 ! Kr_mn = sum_p U_mp κ_p U_np, S_mn = sum_p U_mp √κ_p U_np, H_mp = sum_no S_mn Mr_no S_op, T its eigenvectors:
1074 ! κ, U in fm_K_eigvec, S in fm_K_sqrt, the symmetrised Mr, H, T in fm_H_eigvec and a work matrix,
1075 ! all m x m on the grid of blacs_env
1076 ALLOCATE (kappa(m), lambda(m))
1077 NULLIFY (fm_struct)
1078 CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
1079 nrow_global=m, ncol_global=m)
1080 CALL cp_fm_create(fm_k_eigvec, fm_struct, name="bse_paired_U")
1081 CALL cp_fm_create(fm_k_sqrt, fm_struct, name="bse_paired_S")
1082 CALL cp_fm_create(fm_mr, fm_struct, name="bse_paired_Mr")
1083 CALL cp_fm_create(fm_h, fm_struct, name="bse_paired_H")
1084 CALL cp_fm_create(fm_h_eigvec, fm_struct, name="bse_paired_T")
1085 CALL cp_fm_create(fm_work, fm_struct, name="bse_paired_work")
1086 CALL cp_fm_struct_release(fm_struct)
1087
1088 ! Kr_mn = sum_p U_mp κ_p U_np on the symmetrised live block, (Kr_mn + Kr_nm)/2
1089 CALL cp_fm_to_fm_submat(fm_red_k, fm_work, m, m, 1, 1, 1, 1)
1090 CALL symmetrise_in_place(fm_work, fm_k_eigvec)
1091 kappa(:) = 0.0_dp
1092 CALL choose_eigv_solver(fm_work, fm_k_eigvec, kappa)
1093 kr_min = kappa(1)
1094
1095 IF (kappa(1) <= 0.0_dp) THEN
1096 abba_status = abba_indefinite
1097 ELSE
1098 ! S_mn = sum_p U_mp √κ_p U_np
1099 CALL cp_fm_to_fm(fm_k_eigvec, fm_work)
1100 CALL cp_fm_column_scale(fm_work, sqrt(kappa))
1101 CALL parallel_gemm('N', 'T', m, m, m, 1.0_dp, fm_work, fm_k_eigvec, 0.0_dp, fm_k_sqrt)
1102
1103 ! H_mp = sum_no S_mn Mr_no S_op with Mr_no <- (Mr_no + Mr_on)/2, then T holds its eigenvectors
1104 CALL cp_fm_to_fm_submat(fm_red_m, fm_mr, m, m, 1, 1, 1, 1)
1105 CALL symmetrise_in_place(fm_mr, fm_h_eigvec)
1106 CALL parallel_gemm('N', 'N', m, m, m, 1.0_dp, fm_mr, fm_k_sqrt, 0.0_dp, fm_work)
1107 CALL parallel_gemm('N', 'N', m, m, m, 1.0_dp, fm_k_sqrt, fm_work, 0.0_dp, fm_h)
1108 lambda(:) = 0.0_dp
1109 CALL choose_eigv_solver(fm_h, fm_h_eigvec, lambda)
1110 ! H is congruent to Mr = sum_ia b (M b), so λ_1 <= 0 means (A+B) indefinite on the span of b, or
1111 ! round-off in a nearly dependent basis, which extension_deviation excludes at every extension
1112 IF (lambda(1) <= 0.0_dp) THEN
1113 CALL cp_abort(__location__, &
1114 "BSE Davidson (OLSEN_STRATMANN): negative squared excitation energy. Matrix "// &
1115 "(A+B) is not positive definite on the trial space, or the basis lost its "// &
1116 "orthonormality.")
1117 END IF
1118 omega(1:n_act) = sqrt(lambda(1:n_act))
1119
1120 ! R_mk = sum_n S_mn T_nk/√Ω_k, L_mk = sum_n Mr_mn R_nk/Ω_k
1121 NULLIFY (fm_struct)
1122 CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
1123 nrow_global=m, ncol_global=n_act)
1124 CALL cp_fm_create(fm_r, fm_struct, name="bse_paired_R")
1125 CALL cp_fm_create(fm_l, fm_struct, name="bse_paired_L")
1126 CALL cp_fm_struct_release(fm_struct)
1127 CALL parallel_gemm('N', 'N', m, n_act, m, 1.0_dp, fm_k_sqrt, fm_h_eigvec, 0.0_dp, fm_r)
1128 CALL cp_fm_column_scale(fm_r, 1.0_dp/sqrt(omega(1:n_act)))
1129 CALL parallel_gemm('N', 'N', m, n_act, m, 1.0_dp, fm_mr, fm_r, 0.0_dp, fm_l)
1130 CALL cp_fm_column_scale(fm_l, 1.0_dp/omega(1:n_act))
1131
1132 coef_right(:, :) = 0.0_dp
1133 coef_left(:, :) = 0.0_dp
1134 CALL cp_fm_get_submatrix(fm_r, coef_right(1:m, 1:n_act))
1135 CALL cp_fm_get_submatrix(fm_l, coef_left(1:m, 1:n_act))
1136 CALL cp_fm_release(fm_r)
1137 CALL cp_fm_release(fm_l)
1138 END IF
1139
1140 CALL cp_fm_release(fm_k_eigvec)
1141 CALL cp_fm_release(fm_k_sqrt)
1142 CALL cp_fm_release(fm_mr)
1143 CALL cp_fm_release(fm_h)
1144 CALL cp_fm_release(fm_h_eigvec)
1145 CALL cp_fm_release(fm_work)
1146 DEALLOCATE (kappa, lambda)
1147 CALL timestop(handle)
1148
1149 END SUBROUTINE solve_reduced_paired
1150
1151! **************************************************************************************************
1152!> \brief Makes the nt columns T that follow column m of b orthonormal to b and to each other,
1153!> T_ia,k <- T_ia,k - sum_m b_ia,m sum_jb b_jb,m T_jb,k, then T_ia,k <- sum_l T_ia,l (R^-1)_lk with
1154!> sum_p R_pk R_pl = sum_ia T_ia,k T_ia,l (cholqr2), and stores the products
1155!> (M b)_ia,k = sum_jb (A+B)_ia,jb b_jb,k, (K b)_ia,k = sum_jb (A-B)_ia,jb b_jb,k of the kept
1156!> columns from one application of A and B
1157!> \param mv_env the matrix-free A and B
1158!> \param fm_b basis; columns 1..m are kept, m+1..m+nt hold T on entry and the kept columns on exit
1159!> \param fm_Mb (A+B) b, extended in the same columns
1160!> \param fm_Kb (A-B) b, extended in the same columns
1161!> \param fm_scratch scratch of at least nt columns
1162!> \param m columns of the basis before the extension
1163!> \param nt on exit the number of columns kept
1164!> \param n_kernel counter of the kernel applications
1165!> \param n_dependent counter of the correction vectors dropped as dependent
1166! **************************************************************************************************
1167 SUBROUTINE extend_basis_os(mv_env, fm_b, fm_Mb, fm_Kb, fm_scratch, m, nt, n_kernel, n_dependent)
1168
1169 TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
1170 TYPE(cp_fm_type), INTENT(IN) :: fm_b, fm_mb, fm_kb, fm_scratch
1171 INTEGER, INTENT(IN) :: m
1172 INTEGER, INTENT(INOUT) :: nt, n_kernel, n_dependent
1173
1174 CHARACTER(LEN=*), PARAMETER :: routinen = 'extend_basis_os'
1175
1176 INTEGER :: handle
1177
1178 CALL timeset(routinen, handle)
1179
1180 CALL extend_orthonormal(fm_b, m, nt, mv_env%para_env, n_dependent)
1181
1182 IF (nt > 0) THEN
1183 ! post-condition of the extension, one Gram matrix; the residual norms below rely on it
1184 IF (extension_deviation(fm_b, fm_b, m, nt, mv_env%para_env) > ortho_tol) THEN
1185 cpabort("BSE Davidson: the basis lost orthonormality")
1186 END IF
1187 ! A b into Mb and B b into the scratch, then K b = A b - B b and M b = A b + B b
1188 CALL bse_matvec_apply(mv_env, fm_b, m + 1, nt, fm_mb, fm_scratch, 1)
1189 CALL cp_fm_to_fm(fm_mb, fm_kb, nt, m + 1, m + 1)
1190 CALL columns_axpy(-1.0_dp, fm_scratch, 1, fm_kb, m + 1, nt)
1191 CALL columns_axpy(1.0_dp, fm_scratch, 1, fm_mb, m + 1, nt)
1192 n_kernel = n_kernel + nt
1193 END IF
1194
1195 CALL timestop(handle)
1196
1197 END SUBROUTINE extend_basis_os
1198
1199! **************************************************************************************************
1200!> \brief One convergence step of a driver: dE_k = |E_k - E_k^prev|, the flags from CONVERGENCE_CRITERION,
1201!> the last required root, every root converged once the basis spans all transitions, the roots
1202!> in need, and the line of the iteration table
1203!> \param bse_env CONVERGENCE_CRITERION, EPS_ENERGY and EPS_RES
1204!> \param energies tracked energies in Hartree, ascending; the first SIZE(res) are read
1205!> \param energies_prev the same of the previous iteration
1206!> \param res residual norms of the tracked roots
1207!> \param m subspace dimension
1208!> \param n_ov number of transitions
1209!> \param n_want requested states
1210!> \param iter iteration count
1211!> \param t_iter wall time at the start of the iteration
1212!> \param unit_nr output unit, positive on the writing rank only
1213!> \param dE energy changes in Hartree
1214!> \param conv convergence flags
1215!> \param n_req last required root
1216!> \param need roots in need
1217! **************************************************************************************************
1218 SUBROUTINE convergence_step(bse_env, energies, energies_prev, res, m, n_ov, n_want, iter, t_iter, unit_nr, &
1219 dE, conv, n_req, need)
1220
1221 TYPE(bse_env_type), INTENT(IN) :: bse_env
1222 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: energies, energies_prev, res
1223 INTEGER, INTENT(IN) :: m, n_ov, n_want, iter
1224 REAL(kind=dp), INTENT(IN) :: t_iter
1225 INTEGER, INTENT(IN) :: unit_nr
1226 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: de
1227 LOGICAL, DIMENSION(:), INTENT(OUT) :: conv
1228 INTEGER, INTENT(OUT) :: n_req
1229 LOGICAL, DIMENSION(:), INTENT(OUT) :: need
1230
1231 INTEGER :: n_act
1232
1233 n_act = SIZE(res)
1234 de(:) = abs(energies(1:n_act) - energies_prev(1:n_act))
1235 CALL convergence_flags(bse_env, de, res, conv)
1236 n_req = multiplet_end(energies, n_want, n_act, deg_thresh)
1237 ! a subspace that spans all transitions is exact
1238 IF (m == n_ov) conv(:) = .true.
1239 CALL roots_in_need(conv, energies, res, n_req, need)
1240 ! the masked MAXVAL is -HUGE when no root is in need; the unmasked term covers that case
1241 CALL print_iteration(iter, m, count(conv(1:n_want)), &
1242 max(maxval(res(1:n_req)), maxval(res, mask=need)), &
1243 maxval(de(1:n_req)), energies(1), m_walltime() - t_iter, unit_nr)
1244
1245 END SUBROUTINE convergence_step
1246
1247! **************************************************************************************************
1248!> \brief Convergence flag of every tracked root k from CONVERGENCE_CRITERION: EN dE_k < EPS_ENERGY,
1249!> RES |r_k| < EPS_RES, OR either of the two, AND both
1250!> \param bse_env CONVERGENCE_CRITERION, EPS_ENERGY and EPS_RES
1251!> \param dE energy changes in Hartree
1252!> \param res residual norms
1253!> \param conv one flag per tracked root
1254! **************************************************************************************************
1255 SUBROUTINE convergence_flags(bse_env, dE, res, conv)
1256
1257 TYPE(bse_env_type), INTENT(IN) :: bse_env
1258 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: de, res
1259 LOGICAL, DIMENSION(:), INTENT(OUT) :: conv
1260
1261 INTEGER :: k
1262
1263 DO k = 1, SIZE(conv)
1264 SELECT CASE (bse_env%convergence_criterion)
1265 CASE (bse_iter_en_cond)
1266 conv(k) = de(k) < bse_env%eps_energy
1267 CASE (bse_iter_res_cond)
1268 conv(k) = res(k) < bse_env%eps_res
1269 CASE (bse_iter_or_cond)
1270 conv(k) = de(k) < bse_env%eps_energy .OR. res(k) < bse_env%eps_res
1271 CASE (bse_iter_and_cond)
1272 conv(k) = de(k) < bse_env%eps_energy .AND. res(k) < bse_env%eps_res
1273 END SELECT
1274 END DO
1275
1276 END SUBROUTINE convergence_flags
1277
1278! **************************************************************************************************
1279!> \brief Last member of the multiplet that contains state n_first: the scan stops at the first gap of
1280!> at least thresh, or at n_limit. A multiplet cut by NUM_EXC_EN is thereby taken as a whole
1281!> \param energies ascending energies in Hartree
1282!> \param n_first first state of the scan
1283!> \param n_limit last state the scan may reach
1284!> \param thresh gap that ends a multiplet, in Hartree
1285!> \return the last member, n_first <= n_last <= n_limit
1286! **************************************************************************************************
1287 PURE FUNCTION multiplet_end(energies, n_first, n_limit, thresh) RESULT(n_last)
1288
1289 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: energies
1290 INTEGER, INTENT(IN) :: n_first, n_limit
1291 REAL(kind=dp), INTENT(IN) :: thresh
1292 INTEGER :: n_last
1293
1294 n_last = n_first
1295 DO WHILE (n_last < n_limit)
1296 IF (energies(n_last + 1) - energies(n_last) >= thresh) EXIT
1297 n_last = n_last + 1
1298 END DO
1299
1300 END FUNCTION multiplet_end
1301! **************************************************************************************************
1302!> \brief Roots that still need correction vectors: the unconverged ones up to n_req, and a tracked
1303!> root k above them as long as E_k - |r_k| < E_n_req. An eigenvalue of a symmetric matrix
1304!> lies within |r_k| of E_k, so such a root can still belong to the requested states; the ABBA
1305!> solvers, whose residuals are not those of a symmetric problem, use the same rule as a heuristic
1306!> \param conv convergence flags of the tracked roots
1307!> \param energies tracked energies in Hartree, ascending
1308!> \param res residual norms
1309!> \param n_req last required root
1310!> \param need one flag per tracked root, set here
1311! **************************************************************************************************
1312 SUBROUTINE roots_in_need(conv, energies, res, n_req, need)
1313
1314 LOGICAL, DIMENSION(:), INTENT(IN) :: conv
1315 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: energies, res
1316 INTEGER, INTENT(IN) :: n_req
1317 LOGICAL, DIMENSION(:), INTENT(OUT) :: need
1318
1319 INTEGER :: k
1320
1321 DO k = 1, SIZE(need)
1322 need(k) = .NOT. conv(k)
1323 IF (k > n_req) need(k) = need(k) .AND. energies(k) - res(k) < energies(n_req)
1324 END DO
1325
1326 END SUBROUTINE roots_in_need
1327
1328! **************************************************************************************************
1329!> \brief The lowest roots in need, which receive a correction vector
1330!> \param bse_env CONVERGENCE_CRITERION and EPS_RES
1331!> \param need roots in need, from roots_in_need
1332!> \param res residual norms
1333!> \param n_max largest number of roots
1334!> \param selected indices of the selected roots in ascending order, zero beyond nt
1335!> \param nt number of selected roots
1336! **************************************************************************************************
1337 SUBROUTINE select_roots(bse_env, need, res, n_max, selected, nt)
1338
1339 TYPE(bse_env_type), INTENT(IN) :: bse_env
1340 LOGICAL, DIMENSION(:), INTENT(IN) :: need
1341 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: res
1342 INTEGER, INTENT(IN) :: n_max
1343 INTEGER, DIMENSION(:), INTENT(OUT) :: selected
1344 INTEGER, INTENT(OUT) :: nt
1345
1346 INTEGER :: k
1347
1348 selected(:) = 0
1349 nt = 0
1350 DO k = 1, SIZE(need)
1351 IF (.NOT. need(k)) cycle
1352 ! a root below EPS_RES only waits for its energy criterion
1353 IF (bse_env%convergence_criterion == bse_iter_and_cond .AND. &
1354 res(k) < bse_env%eps_res) cycle
1355 IF (nt == n_max) EXIT
1356 nt = nt + 1
1357 selected(nt) = k
1358 END DO
1359
1360 END SUBROUTINE select_roots
1361
1362! **************************************************************************************************
1363!> \brief Title and column heads of the iteration table
1364!> \param title line printed above the column heads
1365!> \param unit_nr output unit, positive on the writing rank only
1366! **************************************************************************************************
1367 SUBROUTINE print_iteration_header(title, unit_nr)
1368
1369 CHARACTER(LEN=*), INTENT(IN) :: title
1370 INTEGER, INTENT(IN) :: unit_nr
1371
1372 IF (unit_nr > 0) THEN
1373 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1374 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', title
1375 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1376 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
1377 '|r| in a.u., energies in eV, t the wall time of the iteration'
1378 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1379 WRITE (unit_nr, '(T2,A4,T7,A5,T14,A8,T24,A6,T33,A12,T47,A12,T61,A12,T74,A7)') 'BSE|', &
1380 'Iter.', 'Z-space', 'Conv.', 'Max |r|', 'Max dE', 'Lowest E', 't (s)'
1381 END IF
1382
1383 END SUBROUTINE print_iteration_header
1384
1385! **************************************************************************************************
1386!> \brief One line of the iteration table
1387!> \param iter iteration count
1388!> \param m subspace dimension
1389!> \param n_conv converged roots among the requested ones
1390!> \param max_res largest residual norm over the required roots and the roots in need
1391!> \param max_dE in Hartree
1392!> \param lowest_E in Hartree
1393!> \param t_iter wall time of the iteration in seconds
1394!> \param unit_nr output unit, positive on the writing rank only
1395! **************************************************************************************************
1396 SUBROUTINE print_iteration(iter, m, n_conv, max_res, max_dE, lowest_E, t_iter, unit_nr)
1397
1398 INTEGER, INTENT(IN) :: iter, m, n_conv
1399 REAL(kind=dp), INTENT(IN) :: max_res, max_de, lowest_e, t_iter
1400 INTEGER, INTENT(IN) :: unit_nr
1401
1402 IF (unit_nr > 0) THEN
1403 WRITE (unit_nr, '(T2,A4,T7,I5,T14,I8,T24,I6,T33,ES12.4,T47,ES12.4,T61,F12.6,T74,F7.2)') &
1404 'BSE|', iter, m, n_conv, max_res, max_de*evolt, lowest_e*evolt, t_iter
1405 END IF
1406
1407 END SUBROUTINE print_iteration
1408
1409! **************************************************************************************************
1410!> \brief Table of the requested roots and of the tracked roots in need, then the abort for MAX_ITER
1411!> \param energies in Hartree
1412!> \param res residual norms
1413!> \param conv convergence flags
1414!> \param need roots in need
1415!> \param n_req last required root
1416!> \param unit_nr output unit, positive on the writing rank only
1417! **************************************************************************************************
1418 SUBROUTINE abort_unconverged(energies, res, conv, need, n_req, unit_nr)
1419
1420 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: energies, res
1421 LOGICAL, DIMENSION(:), INTENT(IN) :: conv, need
1422 INTEGER, INTENT(IN) :: n_req, unit_nr
1423
1424 INTEGER :: k
1425
1426 IF (unit_nr > 0) THEN
1427 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1428 WRITE (unit_nr, '(T2,A4,T7,A12,T30,A11,T50,A10,T67,A14)') 'BSE|', &
1429 'Excitation n', 'Energy (eV)', '|r| (a.u.)', 'Converged'
1430 DO k = 1, SIZE(need)
1431 IF (k > n_req .AND. .NOT. need(k)) cycle
1432 WRITE (unit_nr, '(T2,A4,T7,I12,T27,F14.6,T46,ES14.4,T67,L14)') 'BSE|', &
1433 k, energies(k)*evolt, res(k), conv(k)
1434 END DO
1435 END IF
1436 ! unconverged roots would enter postprocess_bse as converged states; the table above is the diagnostic
1437 cpabort("BSE Davidson: MAX_ITER reached before convergence")
1438
1439 END SUBROUTINE abort_unconverged
1440
1441! **************************************************************************************************
1442!> \brief Debug table of all tracked roots of a converged solve
1443!> \param energies in Hartree
1444!> \param res residual norms
1445!> \param n_want requested states
1446!> \param n_act tracked states
1447!> \param unit_nr output unit, positive on the writing rank only
1448! **************************************************************************************************
1449 SUBROUTINE print_tracked_roots(energies, res, n_want, n_act, unit_nr)
1450
1451 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: energies, res
1452 INTEGER, INTENT(IN) :: n_want, n_act, unit_nr
1453
1454 INTEGER :: k
1455
1456 IF (unit_nr > 0) THEN
1457 WRITE (unit_nr, '(T2,A10,T13,A13,T27,A9,T45,A18,T67,A14)') 'BSE|DEBUG|', &
1458 'Tracked state', 'Requested', 'Energy (eV)', '|r| (a.u.)'
1459 DO k = 1, n_act
1460 WRITE (unit_nr, '(T2,A10,T13,I13,T27,L9,T45,F18.10,T67,ES14.4)') 'BSE|DEBUG|', &
1461 k, k <= n_want, energies(k)*evolt, res(k)
1462 END DO
1463 END IF
1464
1465 END SUBROUTINE print_tracked_roots
1466
1467! **************************************************************************************************
1468!> \brief Counters of a converged solve, the multiplet note and the warning for a degenerate group
1469!> that reaches the last buffer state
1470!> \param iter iterations run
1471!> \param n_kernel kernel applications
1472!> \param n_restart thick restarts
1473!> \param n_want requested states
1474!> \param n_req last required root
1475!> \param n_act tracked states
1476!> \param n_ov number of transitions
1477!> \param unit_nr output unit, positive on the writing rank only
1478!> \param ab_margin smallest eigenvalue of A-B on the trial space in Hartree
1479!> \param n_dependent correction vectors dropped as dependent, printed when not zero
1480! **************************************************************************************************
1481 SUBROUTINE print_summary(iter, n_kernel, n_restart, n_want, n_req, n_act, n_ov, unit_nr, ab_margin, &
1482 n_dependent)
1483
1484 INTEGER, INTENT(IN) :: iter, n_kernel, n_restart, n_want, &
1485 n_req, n_act, n_ov, unit_nr
1486 REAL(kind=dp), INTENT(IN), OPTIONAL :: ab_margin
1487 INTEGER, INTENT(IN), OPTIONAL :: n_dependent
1488
1489 IF (unit_nr > 0) THEN
1490 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1491 WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Davidson iterations', iter
1492 WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Kernel applications', n_kernel
1493 WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Thick restarts', n_restart
1494 IF (PRESENT(n_dependent)) THEN
1495 IF (n_dependent > 0) THEN
1496 WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Correction vectors dropped as dependent', n_dependent
1497 END IF
1498 END IF
1499 IF (PRESENT(ab_margin)) THEN
1500 WRITE (unit_nr, '(T2,A4,T7,A,T65,F16.6)') 'BSE|', &
1501 'A-B margin on the trial space (eV)', ab_margin*evolt
1502 END IF
1503 IF (n_req > n_want) THEN
1504 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
1505 'Note: the requested states end inside a degenerate multiplet.'
1506 END IF
1507 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1508 END IF
1509
1510 IF (n_req > n_want .AND. n_req == n_act .AND. n_act < n_ov .AND. unit_nr > 0) THEN
1511 CALL cp_warn(__location__, &
1512 "BSE Davidson: the degenerate group at NUM_EXC_EN reaches the last buffer state, "// &
1513 "further members may lie above it. Increase BSE_ITERAT%NUM_BUFFER_STATES.")
1514 END IF
1515
1516 END SUBROUTINE print_summary
1517
1518! **************************************************************************************************
1519!> \brief Number of wanted and tracked states and the block size, clamped to the problem size
1520!> \param bse_env NUM_EXC_EN, NUM_BUFFER_STATES, BLOCK_SIZE and MAX_ITER
1521!> \param n_ov number of transitions
1522!> \param unit_nr output unit, positive on the writing rank only
1523!> \param n_want NUM_EXC_EN, at most n_ov
1524!> \param n_act wanted plus buffer states
1525!> \param block_size BLOCK_SIZE, or MIN(32, n_act) for -1; at most n_act
1526! **************************************************************************************************
1527 SUBROUTINE davidson_sizes(bse_env, n_ov, unit_nr, n_want, n_act, block_size)
1528
1529 TYPE(bse_env_type), INTENT(IN) :: bse_env
1530 INTEGER, INTENT(IN) :: n_ov, unit_nr
1531 INTEGER, INTENT(OUT) :: n_want, n_act, block_size
1532
1533 IF (bse_env%num_exc_en < 1) cpabort("BSE_ITERAT%NUM_EXC_EN must be at least 1")
1534 IF (bse_env%num_buffer_states < 0) cpabort("BSE_ITERAT%NUM_BUFFER_STATES must not be negative")
1535 IF (bse_env%block_size < 1 .AND. bse_env%block_size /= -1) THEN
1536 cpabort("BSE_ITERAT%BLOCK_SIZE must be at least 1, or -1 for the default")
1537 END IF
1538 IF (bse_env%max_iter < 1) cpabort("BSE_ITERAT%MAX_ITER must be at least 1")
1539
1540 n_want = bse_env%num_exc_en
1541 IF (n_want > n_ov) THEN
1542 IF (unit_nr > 0) CALL cp_warn(__location__, &
1543 "BSE_ITERAT%NUM_EXC_EN exceeds the number of transitions and is reduced to it.")
1544 n_want = n_ov
1545 END IF
1546 n_act = min(n_want + bse_env%num_buffer_states, n_ov)
1547 block_size = bse_env%block_size
1548 ! 32 was within 5 % of the best width on the L8 deck at 200 states; wider passes cost restarts
1549 IF (block_size == -1) block_size = min(32, n_act)
1550 IF (block_size > n_act) THEN
1551 IF (unit_nr > 0) CALL cp_warn(__location__, &
1552 "BSE_ITERAT%BLOCK_SIZE exceeds the number of tracked states and is reduced to it.")
1553 block_size = n_act
1554 END IF
1555
1556 END SUBROUTINE davidson_sizes
1557
1558! **************************************************************************************************
1559!> \brief Largest subspace dimension m_max: MAX_SUBSPACE_FACTOR x n_act, or 20 x n_act for -1, with
1560!> the solver's arrays at m_max checked against the memory budget per rank as MEMORY_CHECK says;
1561!> CLAMP lowers m_max to the largest that fits, never below 3 x n_act (4 x n_act for OLSEN_STRATMANN)
1562!> \param bse_env MAX_SUBSPACE_FACTOR, MEMORY_CHECK and MEMORY_BUDGET_GB
1563!> \param mv_env n_ov, the RI share per rank and block_cols
1564!> \param driver driver_tda, driver_mk or driver_os
1565!> \param n_act tracked states
1566!> \param block_size correction vectors per iteration
1567!> \param unit_nr output unit, positive on the writing rank only
1568!> \param m_max the ceiling, 3 or 4 x n_act <= m_max <= n_ov
1569! **************************************************************************************************
1570 SUBROUTINE subspace_ceiling(bse_env, mv_env, driver, n_act, block_size, unit_nr, m_max)
1571
1572 TYPE(bse_env_type), INTENT(IN) :: bse_env
1573 TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
1574 INTEGER, INTENT(IN) :: driver, n_act, block_size, unit_nr
1575 INTEGER, INTENT(OUT) :: m_max
1576
1577 CHARACTER(LEN=16) :: bud_str, fac_str, mem_str
1578 INTEGER :: m_floor, m_hi, m_lo, m_mid, m_start, &
1579 min_fac, n_ri_max, nb
1580 LOGICAL :: from_free, over, skipped
1581 REAL(kind=dp) :: budget_gb, dist_gb, mem_avail_gb, repl_gb
1582
1583 ! the smallest ceiling that leaves room for the restart and one block: TDA and MK restart onto at
1584 ! most 2 n_act vectors and add a block of at most n_act; OS restarts onto right and left pairs
1585 ! under the cap m_max - 2 block_size, which has to hold the current 2 n_act
1586 IF (driver == driver_os) THEN
1587 min_fac = 4
1588 ELSE
1589 min_fac = 3
1590 END IF
1591 WRITE (fac_str, '(I16)') min_fac
1592
1593 ! the start value: the keyword, or for -1 the generous 20 x n_act that the memory check below trims to what fits
1594 IF (bse_env%max_subspace_factor == -1) THEN
1595 m_start = 20*n_act
1596 ELSE IF (bse_env%max_subspace_factor >= min_fac) THEN
1597 m_start = bse_env%max_subspace_factor*n_act
1598 ELSE
1599 CALL cp_abort(__location__, &
1600 "BSE_ITERAT%MAX_SUBSPACE_FACTOR must be at least "//trim(adjustl(fac_str))// &
1601 " (or -1) for the chosen solver")
1602 END IF
1603 m_floor = min(min_fac*n_act, mv_env%n_ov)
1604 m_start = min(max(m_start, m_floor), mv_env%n_ov)
1605 m_max = m_start
1606
1607 ! the widest kernel call of the driver caps the pass width; the largest RI share of any rank
1608 ! keeps the estimate, and so m_max, identical on all ranks
1609 SELECT CASE (driver)
1610 CASE (driver_tda)
1611 nb = n_act
1612 CASE (driver_mk)
1613 nb = block_size
1614 CASE DEFAULT
1615 nb = 2*block_size
1616 END SELECT
1617 IF (mv_env%block_cols > 0) nb = min(nb, mv_env%block_cols)
1618 n_ri_max = mv_env%n_ri_loc
1619 CALL mv_env%para_env%max(n_ri_max)
1620
1621 ! the budget: as given, or the share mem_fraction of the free memory per rank
1622 from_free = bse_env%memory_budget_gb < 0.0_dp
1623 skipped = .false.
1624 budget_gb = bse_env%memory_budget_gb
1625 IF (from_free) THEN
1626 CALL mp_mem_avail_per_rank_gb(mv_env%para_env, mem_avail_gb)
1627 budget_gb = mem_fraction*mem_avail_gb
1628 skipped = mem_avail_gb <= 0.0_dp
1629 IF (skipped .AND. bse_env%memory_check /= bse_memcheck_off .AND. unit_nr > 0) THEN
1630 CALL cp_warn(__location__, "BSE Davidson: free memory not detectable, memory check skipped")
1631 END IF
1632 END IF
1633
1634 CALL davidson_footprint(driver, mv_env, m_max, n_act, block_size, nb, n_ri_max, dist_gb, repl_gb)
1635 over = .NOT. skipped .AND. bse_env%memory_check /= bse_memcheck_off .AND. &
1636 dist_gb + repl_gb > budget_gb
1637
1638 IF (over .AND. bse_env%memory_check == bse_memcheck_clamp) THEN
1639 CALL davidson_footprint(driver, mv_env, m_floor, n_act, block_size, nb, n_ri_max, dist_gb, repl_gb)
1640 IF (dist_gb + repl_gb > budget_gb) THEN
1641 WRITE (mem_str, '(F16.3)') dist_gb + repl_gb
1642 WRITE (bud_str, '(F16.3)') budget_gb
1643 CALL cp_abort(__location__, &
1644 "BSE Davidson: not enough memory for the smallest subspace of "//trim(adjustl(fac_str))// &
1645 " vectors per state: "//trim(adjustl(mem_str))//" GB per MPI rank, the budget is "// &
1646 trim(adjustl(bud_str))//" GB. Raise MEMORY_BUDGET_GB or use more MPI ranks.")
1647 END IF
1648 ! the footprint grows with m: bisect for the largest m within the budget
1649 m_lo = m_floor
1650 m_hi = m_start
1651 DO WHILE (m_hi - m_lo > 1)
1652 m_mid = (m_lo + m_hi)/2
1653 CALL davidson_footprint(driver, mv_env, m_mid, n_act, block_size, nb, n_ri_max, dist_gb, repl_gb)
1654 IF (dist_gb + repl_gb <= budget_gb) THEN
1655 m_lo = m_mid
1656 ELSE
1657 m_hi = m_mid
1658 END IF
1659 END DO
1660 m_max = m_lo
1661 CALL davidson_footprint(driver, mv_env, m_max, n_act, block_size, nb, n_ri_max, dist_gb, repl_gb)
1662 END IF
1663
1664 IF (unit_nr > 0) THEN
1665 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1666 IF (bse_env%max_subspace_factor == -1) THEN
1667 WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Maximum subspace dimension (20 x buffered states)', &
1668 m_start
1669 ELSE
1670 WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Maximum subspace dimension (MAX_SUBSPACE_FACTOR)', &
1671 m_start
1672 END IF
1673 IF (m_max < m_start) THEN
1674 WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Maximum subspace dimension after the memory check', &
1675 m_max
1676 END IF
1677 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory per MPI rank at the subspace ceiling (GB)'
1678 WRITE (unit_nr, '(T2,A4,T9,A,T67,F14.3)') 'BSE|', 'Distributed arrays', dist_gb
1679 WRITE (unit_nr, '(T2,A4,T9,A,T67,F14.3)') 'BSE|', 'Replicated arrays', repl_gb
1680 WRITE (unit_nr, '(T2,A4,T9,A,T67,F14.3)') 'BSE|', 'Total', dist_gb + repl_gb
1681 IF (.NOT. skipped) THEN
1682 IF (from_free) THEN
1683 WRITE (unit_nr, '(T2,A4,T7,A,T67,F14.3)') 'BSE|', 'Fraction of the free memory made available', &
1684 mem_fraction
1685 END IF
1686 WRITE (unit_nr, '(T2,A4,T7,A,T67,F14.3)') 'BSE|', 'Memory budget per MPI rank (GB)', budget_gb
1687 END IF
1688 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1689 IF (bse_env%memory_check == bse_memcheck_off) THEN
1690 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: OFF'
1691 ELSE IF (skipped) THEN
1692 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: skipped, free memory not detectable'
1693 ELSE
1694 SELECT CASE (bse_env%memory_check)
1695 CASE (bse_memcheck_warn)
1696 IF (over) THEN
1697 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: WARN, the total exceeds the budget'
1698 ELSE
1699 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: within the budget, no WARN raised'
1700 END IF
1701 CASE (bse_memcheck_clamp)
1702 IF (over) THEN
1703 WRITE (unit_nr, '(T2,A4,T7,A,I0,A,I0)') 'BSE|', &
1704 'Memory check: CLAMP applied, subspace ceiling reduced from ', m_start, ' to ', m_max
1705 ELSE
1706 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: within the budget, no CLAMP applied'
1707 END IF
1708 CASE (bse_memcheck_abort)
1709 IF (over) THEN
1710 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: ABORT, the total exceeds the budget'
1711 ELSE
1712 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'Memory check: within the budget, no ABORT raised'
1713 END IF
1714 END SELECT
1715 END IF
1716 END IF
1717
1718 IF (over) THEN
1719 WRITE (mem_str, '(F16.3)') dist_gb + repl_gb
1720 WRITE (bud_str, '(F16.3)') budget_gb
1721 SELECT CASE (bse_env%memory_check)
1722 CASE (bse_memcheck_warn)
1723 IF (unit_nr > 0) CALL cp_warn(__location__, &
1724 "BSE Davidson: the solver's arrays at the subspace ceiling need "// &
1725 trim(adjustl(mem_str))//" GB per MPI rank, above the budget of "// &
1726 trim(adjustl(bud_str))//" GB.")
1727 CASE (bse_memcheck_clamp)
1728 IF (unit_nr > 0) CALL cp_warn(__location__, &
1729 "BSE Davidson: subspace ceiling reduced to fit the memory budget; raise "// &
1730 "MEMORY_BUDGET_GB or lower MAX_SUBSPACE_FACTOR to silence this.")
1731 CASE (bse_memcheck_abort)
1732 CALL cp_abort(__location__, &
1733 "BSE Davidson: the solver's arrays at the subspace ceiling need "//trim(adjustl(mem_str))// &
1734 " GB per MPI rank, the budget is "//trim(adjustl(bud_str))//" GB. Raise MEMORY_BUDGET_GB, "// &
1735 "lower MAX_SUBSPACE_FACTOR, or set MEMORY_CHECK CLAMP.")
1736 END SELECT
1737 END IF
1738
1739 END SUBROUTINE subspace_ceiling
1740
1741! **************************************************************************************************
1742!> \brief Memory per MPI rank of the solver's arrays at subspace dimension m, split into the family
1743!> spread over the ranks and the family every rank holds whole. The counts are those of the
1744!> ALLOCATE and cp_fm_create statements of the drivers at a thick restart, their peak:
1745!> distributed 8 B (N_ov (n_basis m + n_work + n_scratch) + (n_red + n_solve) m^2)/P,
1746!> replicated 8 B (c_coef m n_act + m block_size + n_mv N_ov nb + n_ri nb + n_virt virt^2)
1747!> \param driver driver_tda, driver_mk or driver_os
1748!> \param mv_env n_ov, virt and the communicator
1749!> \param m subspace dimension
1750!> \param n_act tracked states
1751!> \param block_size correction vectors per iteration
1752!> \param nb kernel pass width
1753!> \param n_ri RI functions of a rank, the largest over the ranks
1754!> \param dist_GB distributed family, falls as 1/P
1755!> \param repl_GB replicated family, flat in P
1756! **************************************************************************************************
1757 SUBROUTINE davidson_footprint(driver, mv_env, m, n_act, block_size, nb, n_ri, dist_GB, repl_GB)
1758
1759 INTEGER, INTENT(IN) :: driver
1760 TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
1761 INTEGER, INTENT(IN) :: m, n_act, block_size, nb, n_ri
1762 REAL(kind=dp), INTENT(OUT) :: dist_gb, repl_gb
1763
1764 INTEGER :: c_coef, n_basis, n_mv, n_red, n_scratch, &
1765 n_solve, n_virt, n_work
1766 REAL(kind=dp) :: m_real, n_ov_real
1767
1768 ! distributed 8 B (N_ov (n_basis m + n_work + n_scratch) + (n_red + n_solve) m^2)/P: the trial vector
1769 ! arrays, the work and scratch columns, the reduced and the per-solve m x m matrices;
1770 ! replicated 8 B (c_coef m n_act + m block_size + n_mv N_ov nb + n_ri nb + n_virt virt^2): the m x n_act
1771 ! coefficient blocks with the restart work array, the pass buffers of the kernel, virt x virt for B
1772 SELECT CASE (driver)
1773 CASE (driver_tda)
1774 n_basis = 2
1775 n_work = min(2*n_act, m)
1776 n_scratch = 0
1777 n_red = 1
1778 n_solve = 2
1779 c_coef = 8
1780 n_mv = 3
1781 n_virt = 0
1782 CASE (driver_mk)
1783 n_basis = 3
1784 n_work = min(2*n_act, m)
1785 n_scratch = block_size
1786 n_red = 1
1787 n_solve = 2
1788 c_coef = 10
1789 n_mv = 4
1790 n_virt = 1
1791 CASE (driver_os)
1792 n_basis = 3
1793 n_work = max(2*n_act, min(4*n_act, m))
1794 n_scratch = 2*block_size
1795 n_red = 2
1796 n_solve = 6
1797 c_coef = 15
1798 n_mv = 4
1799 n_virt = 1
1800 CASE DEFAULT
1801 ! unreachable for the three driver constants; the DEFAULT keeps every counter defined
1802 cpabort("BSE Davidson: unknown driver in the memory estimate")
1803 END SELECT
1804
1805 ! in REAL: m^2 overflows INTEGER(4) above 46340
1806 m_real = real(m, dp)
1807 n_ov_real = real(mv_env%n_ov, dp)
1808 dist_gb = 8.0e-9_dp*(n_ov_real*real(n_basis*m + n_work + n_scratch, dp) + real(n_red + n_solve, dp)*m_real*m_real)/ &
1809 REAL(mv_env%para_env%num_pe, dp)
1810 repl_gb = 8.0e-9_dp*(real(c_coef, dp)*m_real*real(n_act, dp) + m_real*real(block_size, dp) + &
1811 REAL(n_mv, dp)*n_ov_real*REAL(nb, dp) + REAL(n_ri, dp)*REAL(nb, dp) + &
1812 REAL(n_virt, dp)*REAL(mv_env%virt, dp)**2)
1813
1814 END SUBROUTINE davidson_footprint
1815
1816! **************************************************************************************************
1817!> \brief Unit vectors Z_ia,k = δ_ia,diag_order(k) on the entries diag_order(1..n_used) with the lowest diagonal, never
1818!> cutting a degenerate group. The sorted diagonal is replicated; each rank writes the entries
1819!> that fall on its own rows
1820!> \param diag the diagonal chosen by PRECONDITIONER, replicated
1821!> \param n_guess guess vectors wanted
1822!> \param n_max largest number of guess vectors
1823!> \param deg_thresh_diag entries of the sorted diagonal closer than this form one degenerate group
1824!> \param fm_Z receives the unit vectors in its columns 1..n_used
1825!> \param n_used guess vectors written, n_guess <= n_used <= n_max
1826!> \param n_given leading columns that the caller has filled; they replace the lowest unit vectors
1827! **************************************************************************************************
1828 SUBROUTINE initial_guess(diag, n_guess, n_max, deg_thresh_diag, fm_Z, n_used, n_given)
1829
1830 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: diag
1831 INTEGER, INTENT(IN) :: n_guess, n_max
1832 REAL(kind=dp), INTENT(IN) :: deg_thresh_diag
1833 TYPE(cp_fm_type), INTENT(IN) :: fm_z
1834 INTEGER, INTENT(OUT) :: n_used
1835 INTEGER, INTENT(IN), OPTIONAL :: n_given
1836
1837 INTEGER :: iloc, k, k_first, n_ov, nrow_local
1838 INTEGER, ALLOCATABLE, DIMENSION(:) :: diag_order
1839 INTEGER, DIMENSION(:), POINTER :: row_indices
1840 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: diag_sorted
1841
1842 n_ov = SIZE(diag)
1843 ALLOCATE (diag_sorted(n_ov), diag_order(n_ov))
1844 diag_sorted(:) = diag(:)
1845 CALL sort(diag_sorted, n_ov, diag_order)
1846
1847 n_used = multiplet_end(diag_sorted, n_guess, n_max, deg_thresh_diag)
1848
1849 k_first = 1
1850 IF (PRESENT(n_given)) k_first = n_given + 1
1851
1852 CALL cp_fm_get_info(fm_z, nrow_local=nrow_local, row_indices=row_indices)
1853 DO k = k_first, n_used
1854 DO iloc = 1, nrow_local
1855 IF (row_indices(iloc) == diag_order(k)) fm_z%local_data(iloc, k) = 1.0_dp
1856 END DO
1857 END DO
1858
1859 DEALLOCATE (diag_sorted, diag_order)
1860
1861 END SUBROUTINE initial_guess
1862
1863! **************************************************************************************************
1864!> \brief Initial guess from the exact A (and B) on the n_sub transitions with the lowest diagonal:
1865!> TDA: sum_q A_sub,pq V_qk = θ_k V_pk; ABBA: x_pk = sum_q (K_sub^-1/2)_pq t_qk with
1866!> sum_qrs (K_sub^1/2)_pq M_sub,qr (K_sub^1/2)_rs t_sk = θ_k^2 t_pk, K = A-B, M = A+B. The Ritz
1867!> values θ_k are upper bounds of the true energies (Rayleigh-Ritz for TDA, Cauchy-type
1868!> interlacing of Bai and Li for the linear response problem), kept for check_guess_bound. A
1869!> degenerate group at the upper boundary is taken as a whole. The block and its eigenvectors
1870!> are replicated; each rank writes the entries that fall on its own rows
1871!> \param mv_env the matrix-free A and B, for the exact block
1872!> \param bse_env NUM_GUESS_TRANSITIONS and BSE_DEBUG_PRINT
1873!> \param diag the diagonal chosen by PRECONDITIONER, replicated; its lowest entries select the block
1874!> \param n_guess guess vectors wanted
1875!> \param n_max largest number of guess vectors
1876!> \param deg_thresh_diag Ritz values closer than this form one degenerate group
1877!> \param fm_Z receives the guess vectors in its columns 1..n_used
1878!> \param n_used guess vectors written, n_guess <= n_used <= n_max; 0 for abba_indefinite
1879!> \param theta_sub allocated here, the n_sub Ritz values in Hartree
1880!> \param unit_nr output unit, positive on the writing rank only
1881!> \param n_given leading columns that the caller has filled; they replace the lowest guess vectors
1882!> \param abba_status ABBA problem when present; abba_indefinite if K_sub is not positive definite
1883!> \param ab_margin smallest eigenvalue of K_sub in Hartree
1884!> \param n_pair ABBA: partners y_k = K_sub x_k/θ_k of the guess vectors x_k not given by the caller,
1885!> written to the n_pair columns after n_used; n_used + n_pair stays within n_max
1886! **************************************************************************************************
1887 SUBROUTINE initial_guess_subblock(mv_env, bse_env, diag, n_guess, n_max, deg_thresh_diag, fm_Z, n_used, &
1888 theta_sub, unit_nr, n_given, abba_status, ab_margin, n_pair)
1889
1890 TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
1891 TYPE(bse_env_type), INTENT(IN) :: bse_env
1892 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: diag
1893 INTEGER, INTENT(IN) :: n_guess, n_max
1894 REAL(kind=dp), INTENT(IN) :: deg_thresh_diag
1895 TYPE(cp_fm_type), INTENT(IN) :: fm_z
1896 INTEGER, INTENT(OUT) :: n_used
1897 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
1898 INTENT(OUT) :: theta_sub
1899 INTEGER, INTENT(IN) :: unit_nr
1900 INTEGER, INTENT(IN), OPTIONAL :: n_given
1901 INTEGER, INTENT(OUT), OPTIONAL :: abba_status
1902 REAL(kind=dp), INTENT(OUT), OPTIONAL :: ab_margin
1903 INTEGER, INTENT(OUT), OPTIONAL :: n_pair
1904
1905 CHARACTER(LEN=*), PARAMETER :: routinen = 'initial_guess_subblock'
1906
1907 INTEGER :: handle, iloc, k, k_first, n_ov, n_sub, &
1908 n_top, nrow_local
1909 INTEGER, ALLOCATABLE, DIMENSION(:) :: diag_order, rank_in_block
1910 INTEGER, DIMENSION(:), POINTER :: row_indices
1911 LOGICAL :: ok
1912 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: diag_sorted, eig_k
1913 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: a_sub, b_sub, eigvec_h, eigvec_k, &
1914 guess_x, guess_y, h_sub, k_mhalf, &
1915 k_phalf, k_sub, m_sub
1916 TYPE(mp_para_env_type), POINTER :: para_env
1917
1918 CALL timeset(routinen, handle)
1919
1920 IF (PRESENT(n_pair)) n_pair = 0
1921 para_env => mv_env%para_env
1922 n_ov = SIZE(diag)
1923 n_sub = min(n_ov, max(bse_env%num_guess_transitions, n_guess))
1924 ALLOCATE (diag_sorted(n_ov), diag_order(n_ov), rank_in_block(n_ov))
1925 diag_sorted(:) = diag(:)
1926 CALL sort(diag_sorted, n_ov, diag_order)
1927
1928 ALLOCATE (a_sub(n_sub, n_sub), guess_x(n_sub, n_sub), theta_sub(n_sub))
1929 ok = .true.
1930 IF (PRESENT(abba_status)) THEN
1931 abba_status = abba_ok
1932 ALLOCATE (b_sub(n_sub, n_sub), k_sub(n_sub, n_sub), m_sub(n_sub, n_sub), h_sub(n_sub, n_sub), &
1933 eigvec_k(n_sub, n_sub), eigvec_h(n_sub, n_sub), k_phalf(n_sub, n_sub), &
1934 k_mhalf(n_sub, n_sub), eig_k(n_sub))
1935 CALL bse_matvec_subblock(mv_env, diag_order(1:n_sub), a_sub, b_sub)
1936 ! (K_sub^±1/2)_pq = sum_r U_pr κ_r^±1/2 U_qr from K_sub,pq = sum_r U_pr κ_r U_qr
1937 k_sub(:, :) = a_sub(:, :) - b_sub(:, :)
1938 CALL solve_replicated(k_sub, n_sub, para_env, eig_k, eigvec_k)
1939 ab_margin = eig_k(1)
1940 IF (eig_k(1) <= 0.0_dp) THEN
1941 abba_status = abba_indefinite
1942 ok = .false.
1943 n_used = 0
1944 ELSE
1945 DO k = 1, n_sub
1946 k_phalf(:, k) = eigvec_k(:, k)*sqrt(eig_k(k))
1947 k_mhalf(:, k) = eigvec_k(:, k)/sqrt(eig_k(k))
1948 END DO
1949 k_phalf(:, :) = matmul(k_phalf, transpose(eigvec_k))
1950 k_mhalf(:, :) = matmul(k_mhalf, transpose(eigvec_k))
1951 ! M_sub = 2 A_sub - K_sub, H_ps = sum_qr (K^1/2)_pq M_qr (K^1/2)_rs, x_pk = sum_q (K^-1/2)_pq t_qk
1952 m_sub(:, :) = 2.0_dp*a_sub(:, :) - k_sub(:, :)
1953 h_sub(:, :) = matmul(k_phalf, matmul(m_sub, k_phalf))
1954 CALL solve_replicated(h_sub, n_sub, para_env, theta_sub, eigvec_h)
1955 ! the same failure as in solve_reduced_paired: with K_sub positive definite, M_sub is indefinite
1956 IF (theta_sub(1) <= 0.0_dp) THEN
1957 cpabort("BSE Davidson: the guess block gives a non-positive squared excitation energy")
1958 END IF
1959 theta_sub(:) = sqrt(theta_sub(:))
1960 guess_x(:, :) = matmul(k_mhalf, eigvec_h)
1961 IF (PRESENT(n_pair)) THEN
1962 ! partners guess_y_k = K_sub guess_x_k/θ_k
1963 ALLOCATE (guess_y(n_sub, n_sub))
1964 guess_y(:, :) = matmul(k_sub, guess_x)
1965 DO k = 1, n_sub
1966 guess_y(:, k) = guess_y(:, k)/theta_sub(k)
1967 END DO
1968 END IF
1969 END IF
1970 DEALLOCATE (b_sub, k_sub, m_sub, h_sub, eigvec_k, eigvec_h, k_phalf, k_mhalf, eig_k)
1971 ELSE
1972 CALL bse_matvec_subblock(mv_env, diag_order(1:n_sub), a_sub)
1973 CALL solve_replicated(a_sub, n_sub, para_env, theta_sub, guess_x)
1974 END IF
1975
1976 IF (ok) THEN
1977 k_first = 1
1978 IF (PRESENT(n_given)) k_first = n_given + 1
1979
1980 n_top = min(n_max, n_sub)
1981 ! the partners of the columns k_first..n_used share the n_max columns with the guess
1982 IF (PRESENT(n_pair)) n_top = min(n_top, (n_max + k_first - 1)/2)
1983 n_used = multiplet_end(theta_sub, n_guess, n_top, deg_thresh_diag)
1984
1985 rank_in_block(:) = 0
1986 DO k = 1, n_sub
1987 rank_in_block(diag_order(k)) = k
1988 END DO
1989 CALL cp_fm_get_info(fm_z, nrow_local=nrow_local, row_indices=row_indices)
1990 DO k = k_first, n_used
1991 DO iloc = 1, nrow_local
1992 IF (rank_in_block(row_indices(iloc)) > 0) fm_z%local_data(iloc, k) = guess_x(rank_in_block(row_indices(iloc)), k)
1993 END DO
1994 END DO
1995 IF (PRESENT(n_pair)) THEN
1996 n_pair = min(n_used - k_first + 1, n_max - n_used)
1997 DO k = 1, n_pair
1998 DO iloc = 1, nrow_local
1999 IF (rank_in_block(row_indices(iloc)) > 0) THEN
2000 fm_z%local_data(iloc, n_used + k) = guess_y(rank_in_block(row_indices(iloc)), k_first + k - 1)
2001 END IF
2002 END DO
2003 END DO
2004 END IF
2005
2006 IF (unit_nr > 0) THEN
2007 WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Transitions in the exact guess block', n_sub
2008 IF (bse_env%bse_debug_print) THEN
2009 WRITE (unit_nr, '(T2,A10,T13,A,T67,F14.6)') 'BSE|DEBUG|', &
2010 'Lowest Ritz value of the guess block (eV)', theta_sub(1)*evolt
2011 END IF
2012 END IF
2013 END IF
2014
2015 DEALLOCATE (diag_sorted, diag_order, rank_in_block, a_sub, guess_x)
2016 IF (ALLOCATED(guess_y)) DEALLOCATE (guess_y)
2017
2018 CALL timestop(handle)
2019
2020 END SUBROUTINE initial_guess_subblock
2021
2022! **************************************************************************************************
2023!> \brief A converged energy above the Ritz value of the guess block with the same index is a
2024!> provable miss: both are upper bounds of the true eigenvalue with that index
2025!> \param theta_sub Ritz values of the guess block in Hartree, ascending
2026!> \param energies converged energies in Hartree
2027!> \param n_req roots compared, the required ones
2028!> \param unit_nr output unit, positive on the writing rank only
2029! **************************************************************************************************
2030 SUBROUTINE check_guess_bound(theta_sub, energies, n_req, unit_nr)
2031
2032 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: theta_sub, energies
2033 INTEGER, INTENT(IN) :: n_req, unit_nr
2034
2035 REAL(kind=dp) :: excess
2036
2037 excess = maxval(energies(1:n_req) - theta_sub(1:n_req))
2038 IF (unit_nr > 0) THEN
2039 WRITE (unit_nr, '(T2,A4,T7,A,T65,F16.6)') 'BSE|', &
2040 'Largest excess over the guess-block bound (eV)', excess*evolt
2041 END IF
2042 IF (excess > deg_thresh .AND. unit_nr > 0) THEN
2043 CALL cp_warn(__location__, &
2044 "A converged BSE state lies above the upper bound given by the guess block, "// &
2045 "so a lower state was missed. Raise BSE_ITERAT%NUM_GUESS_TRANSITIONS or "// &
2046 "BSE_ITERAT%NUM_BUFFER_STATES.")
2047 END IF
2048
2049 END SUBROUTINE check_guess_bound
2050
2051! **************************************************************************************************
2052!> \brief sum_q (red_pq + red_qp)/2 c_qk = θ_k c_pk for p, q <= m, red_pq the live block of fm_red, solved
2053!> on the grid of blacs_env: θ_k in theta ascending, c_pk of the first n_vec k gathered into
2054!> coef_ritz
2055!> \param fm_red the reduced matrix; its live block is the first m rows and columns
2056!> \param m dimension of the live block
2057!> \param n_vec eigenvectors gathered; the drivers never read more than n_act
2058!> \param para_env communicator of the grid
2059!> \param blacs_env grid the reduced problem is solved on
2060!> \param theta θ_k in the first m entries, zero beyond
2061!> \param coef_ritz c_pk in the first m rows and n_vec columns, zero beyond
2062! **************************************************************************************************
2063 SUBROUTINE solve_reduced(fm_red, m, n_vec, para_env, blacs_env, theta, coef_ritz)
2064
2065 TYPE(cp_fm_type), INTENT(IN) :: fm_red
2066 INTEGER, INTENT(IN) :: m, n_vec
2067 TYPE(mp_para_env_type), INTENT(IN), TARGET :: para_env
2068 TYPE(cp_blacs_env_type), INTENT(IN), TARGET :: blacs_env
2069 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: theta
2070 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: coef_ritz
2071
2072 CHARACTER(LEN=*), PARAMETER :: routinen = 'solve_reduced'
2073
2074 INTEGER :: handle, n_get
2075 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigval
2076 TYPE(cp_fm_type) :: fm_reduced_eigvec, fm_reduced_sym
2077
2078 CALL timeset(routinen, handle)
2079
2080 CALL reduced_live_block(fm_red, m, para_env, blacs_env, "bse_reduced", fm_reduced_sym)
2081 CALL cp_fm_create(fm_reduced_eigvec, fm_reduced_sym%matrix_struct, name="bse_reduced_vectors")
2082 ! a_pq <- (red_pq + red_qp)/2
2083 CALL symmetrise_in_place(fm_reduced_sym, fm_reduced_eigvec)
2084 ALLOCATE (eigval(m))
2085 eigval(:) = 0.0_dp
2086 CALL choose_eigv_solver(fm_reduced_sym, fm_reduced_eigvec, eigval)
2087
2088 theta(:) = 0.0_dp
2089 theta(1:m) = eigval(:)
2090 coef_ritz(:, :) = 0.0_dp
2091 n_get = min(m, n_vec)
2092 CALL cp_fm_get_submatrix(fm_reduced_eigvec, coef_ritz(1:m, 1:n_get), 1, 1, m, n_get)
2093
2094 DEALLOCATE (eigval)
2095 CALL cp_fm_release(fm_reduced_sym)
2096 CALL cp_fm_release(fm_reduced_eigvec)
2097
2098 CALL timestop(handle)
2099
2100 END SUBROUTINE solve_reduced
2101
2102! **************************************************************************************************
2103!> \brief sum_q (mat_pq + mat_qp)/2 coef_qk = θ_k coef_pk for p, q <= m, a small replicated block solved
2104!> on one rank and broadcast: the fixed-size blocks of the guess and of cholqr2, never the
2105!> reduced matrix
2106!> \param mat the block; its first m rows and columns are read
2107!> \param m dimension of the block
2108!> \param para_env communicator; the solve runs on its source rank
2109!> \param theta ascending
2110!> \param coef all m eigenvectors in the columns
2111! **************************************************************************************************
2112 SUBROUTINE solve_replicated(mat, m, para_env, theta, coef)
2113
2114 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: mat
2115 INTEGER, INTENT(IN) :: m
2116 TYPE(mp_para_env_type), INTENT(IN) :: para_env
2117 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: theta
2118 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: coef
2119
2120 CHARACTER(LEN=*), PARAMETER :: routinen = 'solve_replicated'
2121
2122 INTEGER :: handle
2123 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigval
2124 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_sym
2125
2126 CALL timeset(routinen, handle)
2127
2128 ALLOCATE (mat_sym(m, m), eigval(m))
2129 mat_sym(:, :) = 0.5_dp*(mat(1:m, 1:m) + transpose(mat(1:m, 1:m)))
2130 eigval(:) = 0.0_dp
2131 IF (para_env%is_source()) CALL diamat_all(mat_sym, eigval)
2132 CALL para_env%bcast(mat_sym)
2133 CALL para_env%bcast(eigval)
2134
2135 theta(:) = 0.0_dp
2136 theta(1:m) = eigval(:)
2137 coef(:, :) = 0.0_dp
2138 coef(1:m, 1:m) = mat_sym(:, :)
2139 DEALLOCATE (mat_sym, eigval)
2140
2141 CALL timestop(handle)
2142
2143 END SUBROUTINE solve_replicated
2144
2145! **************************************************************************************************
2146!> \brief a_pq <- (a_pq + a_qp)/2 through a scratch matrix of the same shape: the eigensolvers read one triangle
2147!> \param fm_a symmetrised in place
2148!> \param fm_scratch overwritten with the transpose
2149! **************************************************************************************************
2150 SUBROUTINE symmetrise_in_place(fm_a, fm_scratch)
2151
2152 TYPE(cp_fm_type), INTENT(IN) :: fm_a, fm_scratch
2153
2154 CALL cp_fm_transpose(fm_a, fm_scratch)
2155 CALL cp_fm_scale_and_add(0.5_dp, fm_a, 0.5_dp, fm_scratch)
2156
2157 END SUBROUTINE symmetrise_in_place
2158
2159! **************************************************************************************************
2160!> \brief Orthonormal coefficient vectors of the restart basis from the candidate columns in their
2161!> order of priority, at most n_cap of them: per candidate k, u_p = cand_pk -
2162!> sum_j (sum_q coef_restart_qj cand_qk) coef_restart_pj over the vectors j kept so far, twice,
2163!> then coef_restart_p,new = u_p/|u| unless |u| is below coef_norm_drop
2164!> \param cand candidate coefficient columns, zero-padded to m rows: the current Ritz vectors, then
2165!> the previous ones
2166!> \param m rows of the coefficient vectors, the current basis size
2167!> \param n_cand candidate columns
2168!> \param n_cap largest number of vectors kept
2169!> \param coef_restart m x n_new restart basis
2170!> \param n_new vectors kept, at most MIN(n_cap, n_cand)
2171! **************************************************************************************************
2172 SUBROUTINE restart_basis(cand, m, n_cand, n_cap, coef_restart, n_new)
2173
2174 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: cand
2175 INTEGER, INTENT(IN) :: m, n_cand, n_cap
2176 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: coef_restart
2177 INTEGER, INTENT(OUT) :: n_new
2178
2179 CHARACTER(LEN=*), PARAMETER :: routinen = 'restart_basis'
2180 REAL(kind=dp), PARAMETER :: coef_norm_drop = 1.0e-8_dp
2181
2182 INTEGER :: handle, ipass, j, k
2183 REAL(kind=dp) :: norm
2184 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: basis
2185
2186 CALL timeset(routinen, handle)
2187
2188 ! never more columns than candidates: the paired solver's cap is close to m_max
2189 ALLOCATE (basis(m, min(n_cap, n_cand)))
2190 n_new = 0
2191 DO k = 1, n_cand
2192 IF (n_new == n_cap) EXIT
2193 basis(:, n_new + 1) = cand(1:m, k)
2194 ! a previous vector is nearly parallel to its successor, one pass leaves the basis
2195 ! non-orthonormal at the 1e-6 level and stalls the residuals there
2196 DO ipass = 1, 2
2197 DO j = 1, n_new
2198 basis(:, n_new + 1) = basis(:, n_new + 1) - &
2199 dot_product(basis(:, j), basis(:, n_new + 1))*basis(:, j)
2200 END DO
2201 norm = norm2(basis(:, n_new + 1))
2202 IF (norm < coef_norm_drop) EXIT
2203 basis(:, n_new + 1) = basis(:, n_new + 1)/norm
2204 END DO
2205 ! relative: a candidate that keeps less than this of its unit norm repeats the vectors kept
2206 IF (norm < coef_norm_drop) cycle
2207 n_new = n_new + 1
2208 END DO
2209 coef_restart(:, :) = 0.0_dp
2210 coef_restart(1:m, 1:n_new) = basis(:, 1:n_new)
2211 DEALLOCATE (basis)
2212
2213 CALL timestop(handle)
2214
2215 END SUBROUTINE restart_basis
2216
2217! **************************************************************************************************
2218!> \brief out_ia,k = alpha sum_m V_ia,m coef_mk + beta out_ia,k for m <= nv, k <= nc on the local rows ia
2219!> On the npe x 1 grid whole columns are local, so the rotation is one local DGEMM
2220!> \param fm_V basis, columns 1..nv read
2221!> \param nv columns of fm_V read
2222!> \param coef nv x nc coefficients, replicated
2223!> \param nc columns written
2224!> \param fm_out receives the rotated columns
2225!> \param alpha scale of the product
2226!> \param beta scale of the previous content of fm_out
2227!> \param out_col first column of fm_out written, 1 by default
2228! **************************************************************************************************
2229 SUBROUTINE subspace_rotate(fm_V, nv, coef, nc, fm_out, alpha, beta, out_col)
2230
2231 TYPE(cp_fm_type), INTENT(IN) :: fm_v
2232 INTEGER, INTENT(IN) :: nv
2233 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
2234 INTENT(IN) :: coef
2235 INTEGER, INTENT(IN) :: nc
2236 TYPE(cp_fm_type), INTENT(IN) :: fm_out
2237 REAL(kind=dp), INTENT(IN) :: alpha, beta
2238 INTEGER, INTENT(IN), OPTIONAL :: out_col
2239
2240 CHARACTER(LEN=*), PARAMETER :: routinen = 'subspace_rotate'
2241
2242 INTEGER :: handle, nrow_local, o_col
2243
2244 CALL timeset(routinen, handle)
2245
2246 o_col = 1
2247 IF (PRESENT(out_col)) o_col = out_col
2248 CALL cp_fm_get_info(fm_v, nrow_local=nrow_local)
2249 IF (nrow_local > 0) THEN
2250 CALL dgemm('N', 'N', nrow_local, nc, nv, alpha, fm_v%local_data, SIZE(fm_v%local_data, 1), &
2251 coef, SIZE(coef, 1), beta, fm_out%local_data(:, o_col:o_col + nc - 1), SIZE(fm_out%local_data, 1))
2252 END IF
2253
2254 CALL timestop(handle)
2255
2256 END SUBROUTINE subspace_rotate
2257
2258! **************************************************************************************************
2259!> \brief V_ia,k <- sum_p V_ia,p coef_pk for p <= m, k <= nq through the work matrix, and the same rotation
2260!> of up to two further matrices that share the basis, so that a thick restart is one call
2261!> \param fm_V basis; on exit columns 1..nq hold the rotation and nq+1..m are zero
2262!> \param m columns of fm_V read
2263!> \param coef restart basis, m x nq
2264!> \param nq columns of the rotated basis
2265!> \param fm_work at least nq columns
2266!> \param fm_V2 rotated alike
2267!> \param fm_V3 rotated alike
2268! **************************************************************************************************
2269 SUBROUTINE rotate_in_place(fm_V, m, coef, nq, fm_work, fm_V2, fm_V3)
2270
2271 TYPE(cp_fm_type), INTENT(IN) :: fm_v
2272 INTEGER, INTENT(IN) :: m
2273 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: coef
2274 INTEGER, INTENT(IN) :: nq
2275 TYPE(cp_fm_type), INTENT(IN) :: fm_work
2276 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_v2, fm_v3
2277
2278 CALL rotate_one(fm_v)
2279 IF (PRESENT(fm_v2)) CALL rotate_one(fm_v2)
2280 IF (PRESENT(fm_v3)) CALL rotate_one(fm_v3)
2281
2282 CONTAINS
2283
2284! **************************************************************************************************
2285!> \brief The rotation of one matrix
2286!> \param fm rotated in place
2287! **************************************************************************************************
2288 SUBROUTINE rotate_one(fm)
2289
2290 TYPE(cp_fm_type), INTENT(IN) :: fm
2291
2292 INTEGER :: nrow_local
2293
2294 CALL subspace_rotate(fm, m, coef(1:m, 1:nq), nq, fm_work, 1.0_dp, 0.0_dp)
2295 CALL cp_fm_to_fm(fm_work, fm, nq)
2296 CALL cp_fm_get_info(fm, nrow_local=nrow_local)
2297 fm%local_data(1:nrow_local, nq + 1:m) = 0.0_dp
2298
2299 END SUBROUTINE rotate_one
2300
2301 END SUBROUTINE rotate_in_place
2302! **************************************************************************************************
2303!> \brief T_ia,t0+k-1 <- T_ia,t0+k-1 + alpha U_ia,u0+k-1 for k <= n
2304!> On the npe x 1 grid whole columns are local, so the update needs no communication
2305!> \param alpha scale of the added columns
2306!> \param fm_U columns added
2307!> \param u0 first column read
2308!> \param fm_T columns updated
2309!> \param t0 first column written
2310!> \param n number of columns
2311! **************************************************************************************************
2312 SUBROUTINE columns_axpy(alpha, fm_U, u0, fm_T, t0, n)
2313
2314 REAL(kind=dp), INTENT(IN) :: alpha
2315 TYPE(cp_fm_type), INTENT(IN) :: fm_u
2316 INTEGER, INTENT(IN) :: u0
2317 TYPE(cp_fm_type), INTENT(IN) :: fm_t
2318 INTEGER, INTENT(IN) :: t0, n
2319
2320 INTEGER :: nrow_local
2321
2322 CALL cp_fm_get_info(fm_t, nrow_local=nrow_local)
2323 fm_t%local_data(1:nrow_local, t0:t0 + n - 1) = fm_t%local_data(1:nrow_local, t0:t0 + n - 1) + &
2324 alpha*fm_u%local_data(1:nrow_local, u0:u0 + n - 1)
2325
2326 END SUBROUTINE columns_axpy
2327
2328! **************************************************************************************************
2329!> \brief G_kl = sum_ia U_ia,u0+k-1 V_ia,v0+l-1 for k <= nu, l <= nv, local product and sum over the ranks
2330!> The npe x 1 grid splits the rows over the ranks, which is what the one allreduce completes
2331!> \param fm_U left factor
2332!> \param u0 first column of fm_U
2333!> \param nu columns of fm_U
2334!> \param fm_V right factor
2335!> \param v0 first column of fm_V
2336!> \param nv columns of fm_V
2337!> \param para_env communicator of the row-distributed matrices
2338!> \param gram nu x nv, replicated
2339! **************************************************************************************************
2340 SUBROUTINE subspace_gram(fm_U, u0, nu, fm_V, v0, nv, para_env, gram)
2341
2342 TYPE(cp_fm_type), INTENT(IN) :: fm_u
2343 INTEGER, INTENT(IN) :: u0, nu
2344 TYPE(cp_fm_type), INTENT(IN) :: fm_v
2345 INTEGER, INTENT(IN) :: v0, nv
2346 TYPE(mp_para_env_type), INTENT(IN) :: para_env
2347 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
2348 INTENT(OUT) :: gram
2349
2350 CHARACTER(LEN=*), PARAMETER :: routinen = 'subspace_gram'
2351
2352 INTEGER :: handle, nrow_local
2353
2354 CALL timeset(routinen, handle)
2355
2356 CALL cp_fm_get_info(fm_u, nrow_local=nrow_local)
2357 gram(:, :) = 0.0_dp
2358 IF (nrow_local > 0) THEN
2359 CALL dgemm('T', 'N', nu, nv, nrow_local, 1.0_dp, &
2360 fm_u%local_data(:, u0:u0 + nu - 1), SIZE(fm_u%local_data, 1), &
2361 fm_v%local_data(:, v0:v0 + nv - 1), SIZE(fm_v%local_data, 1), &
2362 0.0_dp, gram, SIZE(gram, 1))
2363 END IF
2364 CALL para_env%sum(gram)
2365
2366 CALL timestop(handle)
2367
2368 END SUBROUTINE subspace_gram
2369
2370! **************************************************************************************************
2371!> \brief Overlaps of the whole basis U with a block of nv columns of V, written into the distributed
2372!> reduced matrix: fm_red_k,v0+l-1 <- sum_ia U_ia,k V_ia,v0+l-1 for k <= nu, l <= nv, in column
2373!> blocks of nb so that the replicated temporary is nu x nb and never nu x nv
2374!> \param fm_U the basis, columns 1..nu
2375!> \param nu columns of fm_U, rows of fm_red written
2376!> \param fm_V the block of columns
2377!> \param v0 first column of fm_V read and first column of fm_red written
2378!> \param nv columns of fm_V read
2379!> \param nb width of one column block
2380!> \param para_env communicator of the row-distributed matrices
2381!> \param fm_red the distributed reduced matrix, on the diagonalisation grid
2382! **************************************************************************************************
2383 SUBROUTINE reduced_gram_blocks(fm_U, nu, fm_V, v0, nv, nb, para_env, fm_red)
2384
2385 TYPE(cp_fm_type), INTENT(IN) :: fm_u
2386 INTEGER, INTENT(IN) :: nu
2387 TYPE(cp_fm_type), INTENT(IN) :: fm_v
2388 INTEGER, INTENT(IN) :: v0, nv, nb
2389 TYPE(mp_para_env_type), INTENT(IN) :: para_env
2390 TYPE(cp_fm_type), INTENT(IN) :: fm_red
2391
2392 CHARACTER(LEN=*), PARAMETER :: routinen = 'reduced_gram_blocks'
2393
2394 INTEGER :: c0, handle, nc
2395 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gram_block
2396
2397 CALL timeset(routinen, handle)
2398
2399 ALLOCATE (gram_block(nu, nb))
2400 DO c0 = v0, v0 + nv - 1, nb
2401 nc = min(nb, v0 + nv - c0)
2402 CALL subspace_gram(fm_u, 1, nu, fm_v, c0, nc, para_env, gram_block(:, 1:nc))
2403 CALL cp_fm_set_submatrix(fm_red, gram_block(:, 1:nc), 1, c0, nu, nc)
2404 END DO
2405 DEALLOCATE (gram_block)
2406
2407 CALL timestop(handle)
2408
2409 END SUBROUTINE reduced_gram_blocks
2410
2411! **************************************************************************************************
2412!> \brief fm_red_p,m+l <- red_block_pl for p <= m+nt and fm_red_m+l,p <- red_block_pl for p <= m,
2413!> l <= nt: the two blocks a basis extension by nt columns adds to the symmetric reduced matrix
2414!> \param fm_red the reduced matrix, at least m+nt rows and columns
2415!> \param red_block (m+nt) x nt, the new columns' overlaps with the whole extended basis
2416!> \param m live dimension before the extension
2417!> \param nt columns added
2418! **************************************************************************************************
2419 SUBROUTINE reduced_extend(fm_red, red_block, m, nt)
2420
2421 TYPE(cp_fm_type), INTENT(IN) :: fm_red
2422 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: red_block
2423 INTEGER, INTENT(IN) :: m, nt
2424
2425 CALL cp_fm_set_submatrix(fm_red, red_block, 1, m + 1, m + nt, nt)
2426 IF (m > 0) CALL cp_fm_set_submatrix(fm_red, red_block(1:m, :), m + 1, 1, nt, m, transpose=.true.)
2427
2428 END SUBROUTINE reduced_extend
2429
2430! **************************************************************************************************
2431!> \brief A new m x m matrix on blacs_env holding the live block fm_red(1:m, 1:m)
2432!> \param fm_red the reduced matrix
2433!> \param m dimension of the live block
2434!> \param para_env communicator of the grid
2435!> \param blacs_env grid of the new matrix
2436!> \param name name of the new matrix
2437!> \param fm_live created here, released by the caller
2438! **************************************************************************************************
2439 SUBROUTINE reduced_live_block(fm_red, m, para_env, blacs_env, name, fm_live)
2440
2441 TYPE(cp_fm_type), INTENT(IN) :: fm_red
2442 INTEGER, INTENT(IN) :: m
2443 TYPE(mp_para_env_type), INTENT(IN), TARGET :: para_env
2444 TYPE(cp_blacs_env_type), INTENT(IN), TARGET :: blacs_env
2445 CHARACTER(LEN=*), INTENT(IN) :: name
2446 TYPE(cp_fm_type), INTENT(OUT) :: fm_live
2447
2448 TYPE(cp_fm_struct_type), POINTER :: fm_struct
2449
2450 NULLIFY (fm_struct)
2451 CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
2452 nrow_global=m, ncol_global=m)
2453 CALL cp_fm_create(fm_live, fm_struct, name=name)
2454 CALL cp_fm_struct_release(fm_struct)
2455 CALL cp_fm_to_fm_submat(fm_red, fm_live, m, m, 1, 1, 1, 1)
2456
2457 END SUBROUTINE reduced_live_block
2458
2459! **************************************************************************************************
2460!> \brief fm_red_ij <- sum_pq coef_restart_pi fm_red_pq coef_restart_qj for p, q <= m and i, j <= k,
2461!> the rest of the matrix zeroed: two distributed products, no replicated m x m
2462!> \param fm_red the reduced matrix, rotated in place
2463!> \param m live dimension before the restart
2464!> \param coef_restart replicated, zero-padded to m rows
2465!> \param k vectors of the restart basis, the live dimension after
2466!> \param para_env communicator of the grid
2467!> \param blacs_env grid of fm_red and of the products
2468! **************************************************************************************************
2469 SUBROUTINE reduced_rotate(fm_red, m, coef_restart, k, para_env, blacs_env)
2470
2471 TYPE(cp_fm_type), INTENT(IN) :: fm_red
2472 INTEGER, INTENT(IN) :: m
2473 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: coef_restart
2474 INTEGER, INTENT(IN) :: k
2475 TYPE(mp_para_env_type), INTENT(IN), TARGET :: para_env
2476 TYPE(cp_blacs_env_type), INTENT(IN), TARGET :: blacs_env
2477
2478 CHARACTER(LEN=*), PARAMETER :: routinen = 'reduced_rotate'
2479
2480 INTEGER :: handle
2481 TYPE(cp_fm_struct_type), POINTER :: fm_struct
2482 TYPE(cp_fm_type) :: fm_gq, fm_live, fm_q, fm_qgq
2483
2484 CALL timeset(routinen, handle)
2485
2486 CALL reduced_live_block(fm_red, m, para_env, blacs_env, "bse_reduced_live", fm_live)
2487 NULLIFY (fm_struct)
2488 CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
2489 nrow_global=m, ncol_global=k)
2490 CALL cp_fm_create(fm_q, fm_struct, name="bse_restart_basis")
2491 CALL cp_fm_create(fm_gq, fm_struct, name="bse_reduced_GQ")
2492 CALL cp_fm_struct_release(fm_struct)
2493 CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
2494 nrow_global=k, ncol_global=k)
2495 CALL cp_fm_create(fm_qgq, fm_struct, name="bse_reduced_QGQ")
2496 CALL cp_fm_struct_release(fm_struct)
2497
2498 CALL cp_fm_set_submatrix(fm_q, coef_restart(1:m, 1:k))
2499 ! GQ_pj = sum_q fm_red_pq coef_restart_qj, then QGQ_ij = sum_p coef_restart_pi GQ_pj
2500 CALL parallel_gemm('N', 'N', m, k, m, 1.0_dp, fm_live, fm_q, 0.0_dp, fm_gq)
2501 CALL parallel_gemm('T', 'N', k, k, m, 1.0_dp, fm_q, fm_gq, 0.0_dp, fm_qgq)
2502 CALL cp_fm_set_all(fm_red, 0.0_dp)
2503 CALL cp_fm_to_fm_submat(fm_qgq, fm_red, k, k, 1, 1, 1, 1)
2504
2505 CALL cp_fm_release(fm_live)
2506 CALL cp_fm_release(fm_q)
2507 CALL cp_fm_release(fm_gq)
2508 CALL cp_fm_release(fm_qgq)
2509
2510 CALL timestop(handle)
2511
2512 END SUBROUTINE reduced_rotate
2513
2514! **************************************************************************************************
2515!> \brief max_kl |sum_ia U_ia,k V_ia,l - δ_kl| over k, l <= m, in column blocks of nb
2516!> \param fm_U left basis
2517!> \param fm_V right basis, fm_U for the Euclidean case
2518!> \param m columns compared
2519!> \param nb width of one column block
2520!> \param para_env communicator of the row-distributed matrices
2521!> \return the largest deviation
2522! **************************************************************************************************
2523 FUNCTION orthonormality_deviation(fm_U, fm_V, m, nb, para_env) RESULT(dev)
2524
2525 TYPE(cp_fm_type), INTENT(IN) :: fm_u, fm_v
2526 INTEGER, INTENT(IN) :: m, nb
2527 TYPE(mp_para_env_type), INTENT(IN) :: para_env
2528 REAL(kind=dp) :: dev
2529
2530 INTEGER :: c0, l, nc
2531 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gram_block
2532
2533 ALLOCATE (gram_block(m, nb))
2534 dev = 0.0_dp
2535 DO c0 = 1, m, nb
2536 nc = min(nb, m - c0 + 1)
2537 CALL subspace_gram(fm_u, 1, m, fm_v, c0, nc, para_env, gram_block(:, 1:nc))
2538 DO l = 1, nc
2539 gram_block(c0 + l - 1, l) = gram_block(c0 + l - 1, l) - 1.0_dp
2540 END DO
2541 dev = max(dev, maxval(abs(gram_block(:, 1:nc))))
2542 END DO
2543 DEALLOCATE (gram_block)
2544
2545 END FUNCTION orthonormality_deviation
2546
2547! **************************************************************************************************
2548!> \brief Euclidean norms of nt columns of a row-distributed matrix
2549!> Each rank sums the squares of its own rows of the npe x 1 grid, then one allreduce
2550!> \param fm_T row-distributed matrix
2551!> \param first_col first column measured
2552!> \param nt number of columns
2553!> \param para_env communicator of the rows
2554!> \param norms nt Euclidean norms
2555! **************************************************************************************************
2556 SUBROUTINE column_norms(fm_T, first_col, nt, para_env, norms)
2557
2558 TYPE(cp_fm_type), INTENT(IN) :: fm_t
2559 INTEGER, INTENT(IN) :: first_col, nt
2560 TYPE(mp_para_env_type), INTENT(IN) :: para_env
2561 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: norms
2562
2563 INTEGER :: k, nrow_local
2564
2565 CALL cp_fm_get_info(fm_t, nrow_local=nrow_local)
2566 norms(:) = 0.0_dp
2567 DO k = 1, nt
2568 norms(k) = sum(fm_t%local_data(1:nrow_local, first_col + k - 1)**2)
2569 END DO
2570 CALL para_env%sum(norms)
2571 norms(:) = sqrt(norms(:))
2572
2573 END SUBROUTINE column_norms
2574
2575! **************************************************************************************************
2576!> \brief Davidson corrections t_ia,k = r_ia,k/(d_ia - θ_k) of the selected roots k, written to fm_Z
2577!> Each rank divides the entries of its own rows of the npe x 1 grid, no communication
2578!> \param fm_R residuals, column r_offset + k belongs to root k
2579!> \param selected indices of the roots corrected
2580!> \param theta Ritz values, indexed by root
2581!> \param diag d_ia, replicated
2582!> \param fm_Z receives the corrections
2583!> \param first_col first column of fm_Z written
2584!> \param r_offset columns of fm_R before the residuals, 0 by default
2585! **************************************************************************************************
2586 SUBROUTINE davidson_corrections(fm_R, selected, theta, diag, fm_Z, first_col, r_offset)
2587
2588 TYPE(cp_fm_type), INTENT(IN) :: fm_r
2589 INTEGER, DIMENSION(:), INTENT(IN) :: selected
2590 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: theta, diag
2591 TYPE(cp_fm_type), INTENT(IN) :: fm_z
2592 INTEGER, INTENT(IN) :: first_col
2593 INTEGER, INTENT(IN), OPTIONAL :: r_offset
2594
2595 CHARACTER(LEN=*), PARAMETER :: routinen = 'davidson_corrections'
2596 REAL(kind=dp), PARAMETER :: eref_scale = 0.99_dp, threshold = 16.0_dp*epsilon(1.0_dp)
2597
2598 INTEGER :: handle, iloc, it, k, nrow_local, r_off
2599 INTEGER, DIMENSION(:), POINTER :: row_indices
2600 REAL(kind=dp) :: denom
2601
2602 CALL timeset(routinen, handle)
2603
2604 r_off = 0
2605 IF (PRESENT(r_offset)) r_off = r_offset
2606 CALL cp_fm_get_info(fm_z, nrow_local=nrow_local, row_indices=row_indices)
2607 DO it = 1, SIZE(selected)
2608 k = selected(it)
2609 DO iloc = 1, nrow_local
2610 denom = diag(row_indices(iloc)) - theta(k)
2611 ! near-degeneracy guard as in tddfpt_compute_residual_vects: threshold detects a numerically
2612 ! zero d_ia - θ_k, which the shift (1 - eref_scale) θ_k replaces by the denominator at 0.99 θ_k
2613 IF (abs(denom) < threshold) denom = denom + (1.0_dp - eref_scale)*theta(k)
2614 fm_z%local_data(iloc, first_col + it - 1) = fm_r%local_data(iloc, r_off + k)/denom
2615 END DO
2616 END DO
2617
2618 CALL timestop(handle)
2619
2620 END SUBROUTINE davidson_corrections
2621
2622! **************************************************************************************************
2623!> \brief T_ia,k <- T_ia,k - sum_m V_ia,m sum_jb U_jb,m T_jb,k for m <= nv and nt columns k, twice;
2624!> U = V for orthonormal columns V, else the dual basis with sum_ia U_ia,m V_ia,n = δ_mn
2625!> One allreduce per pass inside the Gram, the subtraction stays on the rows of the npe x 1 grid
2626!> \param fm_T columns projected
2627!> \param first_col first column of fm_T projected
2628!> \param nt number of columns projected
2629!> \param fm_V the basis projected out
2630!> \param nv columns of fm_V
2631!> \param para_env communicator of the row-distributed matrices
2632!> \param fm_dual U
2633! **************************************************************************************************
2634 SUBROUTINE project_out(fm_T, first_col, nt, fm_V, nv, para_env, fm_dual)
2635
2636 TYPE(cp_fm_type), INTENT(IN) :: fm_t
2637 INTEGER, INTENT(IN) :: first_col, nt
2638 TYPE(cp_fm_type), INTENT(IN) :: fm_v
2639 INTEGER, INTENT(IN) :: nv
2640 TYPE(mp_para_env_type), INTENT(IN) :: para_env
2641 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_dual
2642
2643 CHARACTER(LEN=*), PARAMETER :: routinen = 'project_out'
2644
2645 INTEGER :: handle, ipass, nrow_local
2646 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff
2647
2648 IF (nt == 0 .OR. nv == 0) RETURN
2649 CALL timeset(routinen, handle)
2650 CALL cp_fm_get_info(fm_t, nrow_local=nrow_local)
2651 ALLOCATE (coeff(nv, nt))
2652 DO ipass = 1, 2
2653 IF (PRESENT(fm_dual)) THEN
2654 CALL subspace_gram(fm_dual, 1, nv, fm_t, first_col, nt, para_env, coeff)
2655 ELSE
2656 CALL subspace_gram(fm_v, 1, nv, fm_t, first_col, nt, para_env, coeff)
2657 END IF
2658 IF (nrow_local > 0) THEN
2659 CALL dgemm('N', 'N', nrow_local, nt, nv, -1.0_dp, &
2660 fm_v%local_data, SIZE(fm_v%local_data, 1), coeff, nv, 1.0_dp, &
2661 fm_t%local_data(:, first_col:first_col + nt - 1), SIZE(fm_t%local_data, 1))
2662 END IF
2663 END DO
2664 DEALLOCATE (coeff)
2665 CALL timestop(handle)
2666
2667 END SUBROUTINE project_out
2668
2669! **************************************************************************************************
2670!> \brief Normalises nt columns and removes those that vanished in the projection
2671!> Only the norms communicate; the scaling and the zero-fill stay on the local rows
2672!> \param fm_T row-distributed matrix; the kept columns move to the front of the block
2673!> \param first_col first column of the block
2674!> \param nt on exit the number of columns kept
2675!> \param para_env communicator of the rows
2676! **************************************************************************************************
2677 SUBROUTINE drop_small_columns(fm_T, first_col, nt, para_env)
2678
2679 TYPE(cp_fm_type), INTENT(IN) :: fm_t
2680 INTEGER, INTENT(IN) :: first_col
2681 INTEGER, INTENT(INOUT) :: nt
2682 TYPE(mp_para_env_type), INTENT(IN) :: para_env
2683
2684 REAL(kind=dp), PARAMETER :: col_norm_drop = 1.0e-10_dp
2685
2686 INTEGER :: k, n_kept, nrow_local
2687 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: norms
2688
2689 IF (nt == 0) RETURN
2690 CALL cp_fm_get_info(fm_t, nrow_local=nrow_local)
2691 ALLOCATE (norms(nt))
2692 CALL column_norms(fm_t, first_col, nt, para_env, norms)
2693 n_kept = 0
2694 DO k = 1, nt
2695 ! absolute: the corrections r/(d - θ) arrive unnormalised, a column below this is numerically zero
2696 IF (norms(k) < col_norm_drop) cycle
2697 n_kept = n_kept + 1
2698 fm_t%local_data(1:nrow_local, first_col + n_kept - 1) = &
2699 fm_t%local_data(1:nrow_local, first_col + k - 1)/norms(k)
2700 END DO
2701 fm_t%local_data(1:nrow_local, first_col + n_kept:first_col + nt - 1) = 0.0_dp
2702 nt = n_kept
2703 DEALLOCATE (norms)
2704
2705 END SUBROUTINE drop_small_columns
2706
2707! **************************************************************************************************
2708!> \brief Orthonormalises nt columns by two Cholesky QR passes on the replicated Gram matrix,
2709!> T_ia,k <- sum_l T_ia,l (R^-1)_lk with sum_p R_pk R_pl = G_kl = sum_ia T_ia,k T_ia,l. A
2710!> failed factorization means a numerically dependent block: the directions of the Gram matrix
2711!> below the drop tolerance are removed (canonical orthogonalisation), nt shrinks, and the
2712!> passes start again. The Gram is the only communication: R is factorised on one rank and
2713!> broadcast, the triangular solve runs on the rows of the npe x 1 grid
2714!> \param fm_T the block to orthonormalise, columns first_col .. first_col+nt-1
2715!> \param first_col first column of the block
2716!> \param nt columns of the block; on exit the number kept
2717!> \param para_env communicator of the row-distributed matrix
2718!> \param ok .FALSE. when two drop rounds still leave no factorisable Gram matrix, or no direction survives
2719!> \param n_dropped columns removed as dependent
2720! **************************************************************************************************
2721 SUBROUTINE cholqr2(fm_T, first_col, nt, para_env, ok, n_dropped)
2722
2723 TYPE(cp_fm_type), INTENT(IN) :: fm_t
2724 INTEGER, INTENT(IN) :: first_col
2725 INTEGER, INTENT(INOUT) :: nt
2726 TYPE(mp_para_env_type), INTENT(IN) :: para_env
2727 LOGICAL, INTENT(OUT) :: ok
2728 INTEGER, INTENT(OUT), OPTIONAL :: n_dropped
2729
2730 CHARACTER(LEN=*), PARAMETER :: routinen = 'cholqr2'
2731 REAL(kind=dp), PARAMETER :: eig_drop_rel = 1.0e-10_dp
2732
2733 INTEGER :: handle, info, ipass, k, n_drop_calls, &
2734 n_keep, nrow_local
2735 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eig
2736 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: chol_factor, gram, gram_eigvec, t_reduced
2737
2738 ok = .true.
2739 IF (PRESENT(n_dropped)) n_dropped = 0
2740 IF (nt == 0) RETURN
2741 CALL timeset(routinen, handle)
2742 CALL cp_fm_get_info(fm_t, nrow_local=nrow_local)
2743 ALLOCATE (gram(nt, nt), chol_factor(nt, nt))
2744
2745 n_drop_calls = 0
2746 ipass = 0
2747 DO WHILE (ipass < 2)
2748 ipass = ipass + 1
2749 ! G_kl = sum_ia T_ia,k T_ia,l = sum_p R_pk R_pl
2750 CALL subspace_gram(fm_t, first_col, nt, fm_t, first_col, nt, para_env, gram)
2751 chol_factor(:, :) = gram(:, :)
2752 info = 0
2753 IF (para_env%is_source()) CALL dpotrf('U', nt, chol_factor, nt, info)
2754 CALL para_env%bcast(info)
2755 IF (info /= 0) THEN
2756 n_drop_calls = n_drop_calls + 1
2757 IF (n_drop_calls > 2) THEN
2758 ok = .false.
2759 EXIT
2760 END IF
2761 ! T_ia,k <- sum_l T_ia,l U_lk Λ_k^-1/2 over the eigenpairs of G above the drop tolerance,
2762 ! eig_drop_rel of the largest eigenvalue: directions below it are numerically dependent
2763 ALLOCATE (eig(nt), gram_eigvec(nt, nt))
2764 CALL solve_replicated(gram, nt, para_env, eig, gram_eigvec)
2765 n_keep = count(eig > eig_drop_rel*eig(nt))
2766 IF (n_keep == 0) THEN
2767 ok = .false.
2768 DEALLOCATE (eig, gram_eigvec)
2769 EXIT
2770 END IF
2771 DO k = 1, n_keep
2772 gram_eigvec(:, nt - n_keep + k) = gram_eigvec(:, nt - n_keep + k)/sqrt(eig(nt - n_keep + k))
2773 END DO
2774 IF (nrow_local > 0) THEN
2775 ALLOCATE (t_reduced(nrow_local, n_keep))
2776 CALL dgemm('N', 'N', nrow_local, n_keep, nt, 1.0_dp, &
2777 fm_t%local_data(:, first_col:first_col + nt - 1), SIZE(fm_t%local_data, 1), &
2778 gram_eigvec(:, nt - n_keep + 1:nt), nt, 0.0_dp, t_reduced, nrow_local)
2779 fm_t%local_data(1:nrow_local, first_col:first_col + n_keep - 1) = t_reduced(:, :)
2780 fm_t%local_data(1:nrow_local, first_col + n_keep:first_col + nt - 1) = 0.0_dp
2781 DEALLOCATE (t_reduced)
2782 END IF
2783 IF (PRESENT(n_dropped)) n_dropped = n_dropped + nt - n_keep
2784 nt = n_keep
2785 DEALLOCATE (eig, gram_eigvec, gram, chol_factor)
2786 ALLOCATE (gram(nt, nt), chol_factor(nt, nt))
2787 ipass = 0
2788 cycle
2789 END IF
2790 ! T_ia,k <- sum_l T_ia,l (R^-1)_lk on the local rows
2791 CALL para_env%bcast(chol_factor)
2792 IF (nrow_local > 0) THEN
2793 CALL dtrsm('R', 'U', 'N', 'N', nrow_local, nt, 1.0_dp, chol_factor, nt, &
2794 fm_t%local_data(:, first_col:first_col + nt - 1), SIZE(fm_t%local_data, 1))
2795 END IF
2796 END DO
2797
2798 DEALLOCATE (gram, chol_factor)
2799 CALL timestop(handle)
2800
2801 END SUBROUTINE cholqr2
2802
2803! **************************************************************************************************
2804!> \brief Makes the nt columns T after column m of fm_V orthonormal to the basis and to each other:
2805!> T_ia,k <- T_ia,k - sum_p V_ia,p sum_jb U_jb,p T_jb,k twice, the columns that vanished dropped, then
2806!> the Euclidean cholqr2; nt shrinks by the vanished and the dependent columns
2807!> \param fm_V basis in columns 1..m, the block T in columns m+1..m+nt
2808!> \param m columns of the basis
2809!> \param nt columns of T; on exit the number kept
2810!> \param para_env communicator of the row-distributed matrix
2811!> \param n_dependent counter of the columns dropped as dependent, incremented
2812!> \param fm_dual U, the dual basis with sum_ia U_ia,p V_ia,q = δ_pq; fm_V itself when absent
2813! **************************************************************************************************
2814 SUBROUTINE extend_orthonormal(fm_V, m, nt, para_env, n_dependent, fm_dual)
2815
2816 TYPE(cp_fm_type), INTENT(IN) :: fm_v
2817 INTEGER, INTENT(IN) :: m
2818 INTEGER, INTENT(INOUT) :: nt
2819 TYPE(mp_para_env_type), INTENT(IN) :: para_env
2820 INTEGER, INTENT(INOUT) :: n_dependent
2821 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_dual
2822
2823 INTEGER :: n_dropped
2824 LOGICAL :: ok
2825
2826 CALL project_out(fm_v, m + 1, nt, fm_v, m, para_env, fm_dual)
2827 CALL drop_small_columns(fm_v, m + 1, nt, para_env)
2828 CALL cholqr2(fm_v, m + 1, nt, para_env, ok, n_dropped)
2829 ! cholqr2 gives up after two rounds of dropping dependent directions, or when none is left
2830 IF (.NOT. ok) cpabort("BSE Davidson: orthonormalisation broke down")
2831 n_dependent = n_dependent + n_dropped
2832
2833 END SUBROUTINE extend_orthonormal
2834
2835! **************************************************************************************************
2836!> \brief Largest deviation of the nt new columns from orthonormality against the whole basis,
2837!> sum_ia V_ia,m D_ia,n - δ_mn for m = 1..nv+nt, n = nv+1..nv+nt; D = V in the Euclidean case
2838!> \param fm_V the basis with its nt new columns
2839!> \param fm_dual D, fm_V itself in the Euclidean case
2840!> \param nv columns of the basis before the extension
2841!> \param nt new columns
2842!> \param para_env communicator of the row-distributed matrices
2843!> \return the largest deviation
2844! **************************************************************************************************
2845 FUNCTION extension_deviation(fm_V, fm_dual, nv, nt, para_env) RESULT(dev)
2846
2847 TYPE(cp_fm_type), INTENT(IN) :: fm_v, fm_dual
2848 INTEGER, INTENT(IN) :: nv, nt
2849 TYPE(mp_para_env_type), INTENT(IN) :: para_env
2850 REAL(kind=dp) :: dev
2851
2852 INTEGER :: k
2853 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gram
2854
2855 dev = 0.0_dp
2856 IF (nt == 0) RETURN
2857 ALLOCATE (gram(nv + nt, nt))
2858 CALL subspace_gram(fm_v, 1, nv + nt, fm_dual, nv + 1, nt, para_env, gram)
2859 DO k = 1, nt
2860 gram(nv + k, k) = gram(nv + k, k) - 1.0_dp
2861 END DO
2862 dev = maxval(abs(gram))
2863 DEALLOCATE (gram)
2864
2865 END FUNCTION extension_deviation
2866
2867! **************************************************************************************************
2868!> \brief Debug check of a Davidson result against the full diagonalization of the explicit
2869!> matrices: energies, and per multiplet the smallest singular value of X_ref^T X for the
2870!> TDA, or of X_ref^T X - Y_ref^T Y for ABBA. ABBA diagonalizes the Hermitian reduction
2871!> C = (A-B)^1/2 (A+B) (A-B)^1/2 with C T^n = (Ω^n)^2 T^n,
2872!> (X+Y)^n = (Ω^n)^-1/2 (A-B)^1/2 T^n, (X-Y)^n = (Ω^n)^1/2 (A-B)^-1/2 T^n.
2873!> A failed reference solve aborts: this is debug output with nothing to fall back on
2874!> \param fm_A_explicit A from create_A_and_B, N_ov x N_ov
2875!> \param exc_ens Davidson energies in Hartree
2876!> \param fm_X Davidson X_ia^n, column n
2877!> \param bse_env handed to create_hermitian_form_of_ABBA
2878!> \param unit_nr output unit, positive on the writing rank only
2879!> \param fm_B_explicit present for an ABBA result, together with fm_Y
2880!> \param fm_Y Davidson Y_ia^n, with fm_B_explicit
2881! **************************************************************************************************
2882 SUBROUTINE bse_davidson_refcheck(fm_A_explicit, exc_ens, fm_X, bse_env, unit_nr, fm_B_explicit, fm_Y)
2883
2884 TYPE(cp_fm_type), INTENT(IN) :: fm_a_explicit
2885 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: exc_ens
2886 TYPE(cp_fm_type), INTENT(IN) :: fm_x
2887 TYPE(bse_env_type), INTENT(IN) :: bse_env
2888 INTEGER, INTENT(IN) :: unit_nr
2889 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_b_explicit, fm_y
2890
2891 CHARACTER(LEN=*), PARAMETER :: routinen = 'bse_davidson_refcheck'
2892
2893 INTEGER :: diag_info, handle, k, n_ov, n_ref, n_want
2894 REAL(kind=dp) :: min_overlap
2895 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ref_ens
2896 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: diff_ref, sum_ref, x_dav, x_ref, y_dav, &
2897 y_ref
2898 TYPE(cp_fm_type) :: fm_c, fm_eigvec, fm_inv_sqrt_a_minus_b, &
2899 fm_sqrt_a_minus_b, fm_work
2900
2901 CALL timeset(routinen, handle)
2902
2903 ! B and Y arrive together or not at all
2904 cpassert(PRESENT(fm_b_explicit) .EQV. PRESENT(fm_y))
2905 IF (unit_nr > 0) THEN
2906 WRITE (unit_nr, '(T2,A4)') 'BSE|'
2907 IF (PRESENT(fm_b_explicit)) THEN
2908 WRITE (unit_nr, '(T2,A10,T13,A)') 'BSE|DEBUG|', 'Reference check against the explicit A and B'
2909 ELSE
2910 WRITE (unit_nr, '(T2,A10,T13,A)') 'BSE|DEBUG|', 'Reference check against the explicit A'
2911 END IF
2912 END IF
2913 CALL cp_fm_get_info(fm_a_explicit, nrow_global=n_ov)
2914 n_want = SIZE(exc_ens)
2915 ALLOCATE (ref_ens(n_ov))
2916
2917 IF (PRESENT(fm_b_explicit)) THEN
2918 CALL create_hermitian_form_of_abba(fm_a_explicit, fm_b_explicit, fm_c, fm_sqrt_a_minus_b, &
2919 fm_inv_sqrt_a_minus_b, unit_nr, bse_env, 0.0_dp)
2920 CALL cp_fm_create(fm_eigvec, fm_c%matrix_struct)
2921 CALL choose_eigv_solver(fm_c, fm_eigvec, ref_ens, diag_info)
2922 IF (diag_info /= 0) cpabort("Reference diagonalization of C failed in the BSE Davidson check")
2923 IF (ref_ens(1) <= 0.0_dp) THEN
2924 CALL cp_abort(__location__, &
2925 "Reference matrix C has a non-positive eigenvalue in the BSE Davidson check")
2926 END IF
2927 ref_ens(:) = sqrt(ref_ens(:))
2928
2929 n_ref = multiplet_end(ref_ens, n_want, n_ov, deg_thresh)
2930 ALLOCATE (sum_ref(n_ov, n_ref), diff_ref(n_ov, n_ref))
2931 CALL cp_fm_create(fm_work, fm_c%matrix_struct)
2932 CALL parallel_gemm("N", "N", n_ov, n_ref, n_ov, 1.0_dp, fm_sqrt_a_minus_b, fm_eigvec, 0.0_dp, &
2933 fm_work)
2934 CALL cp_fm_get_submatrix(fm_work, sum_ref, 1, 1, n_ov, n_ref)
2935 CALL parallel_gemm("N", "N", n_ov, n_ref, n_ov, 1.0_dp, fm_inv_sqrt_a_minus_b, fm_eigvec, &
2936 0.0_dp, fm_work)
2937 CALL cp_fm_get_submatrix(fm_work, diff_ref, 1, 1, n_ov, n_ref)
2938 CALL cp_fm_release(fm_work)
2939 CALL cp_fm_release(fm_c)
2940 CALL cp_fm_release(fm_sqrt_a_minus_b)
2941 CALL cp_fm_release(fm_inv_sqrt_a_minus_b)
2942
2943 ALLOCATE (x_ref(n_ov, n_ref), y_ref(n_ov, n_ref), y_dav(n_ov, n_want))
2944 DO k = 1, n_ref
2945 sum_ref(:, k) = sum_ref(:, k)/sqrt(ref_ens(k))
2946 diff_ref(:, k) = diff_ref(:, k)*sqrt(ref_ens(k))
2947 END DO
2948 x_ref(:, :) = 0.5_dp*(sum_ref(:, :) + diff_ref(:, :))
2949 y_ref(:, :) = 0.5_dp*(sum_ref(:, :) - diff_ref(:, :))
2950 DEALLOCATE (sum_ref, diff_ref)
2951 CALL cp_fm_get_submatrix(fm_y, y_dav)
2952 ELSE
2953 ! the eigensolver destroys its input
2954 CALL cp_fm_create(fm_work, fm_a_explicit%matrix_struct)
2955 CALL cp_fm_to_fm(fm_a_explicit, fm_work)
2956 CALL cp_fm_create(fm_eigvec, fm_a_explicit%matrix_struct)
2957 CALL choose_eigv_solver(fm_work, fm_eigvec, ref_ens, diag_info)
2958 IF (diag_info /= 0) cpabort("Reference diagonalization of A failed in the BSE Davidson check")
2959 CALL cp_fm_release(fm_work)
2960
2961 n_ref = multiplet_end(ref_ens, n_want, n_ov, deg_thresh)
2962 ALLOCATE (x_ref(n_ov, n_ref))
2963 CALL cp_fm_get_submatrix(fm_eigvec, x_ref, 1, 1, n_ov, n_ref)
2964 END IF
2965 CALL cp_fm_release(fm_eigvec)
2966
2967 ALLOCATE (x_dav(n_ov, n_want))
2968 CALL cp_fm_get_submatrix(fm_x, x_dav)
2969
2970 ! Y_ref and Y_dav stay unallocated for the TDA and then count as absent (Fortran 2008)
2971 CALL min_multiplet_overlap(ref_ens, n_want, n_ref, x_ref, x_dav, min_overlap, y_ref, y_dav)
2972 CALL print_refcheck(maxval(abs(exc_ens(:) - ref_ens(1:n_want))), min_overlap, unit_nr)
2973
2974 DEALLOCATE (ref_ens, x_ref, x_dav)
2975 IF (ALLOCATED(y_ref)) DEALLOCATE (y_ref, y_dav)
2976
2977 CALL timestop(handle)
2978
2979 END SUBROUTINE bse_davidson_refcheck
2980
2981! **************************************************************************************************
2982!> \brief Smallest singular value over the multiplets of O_kl = sum_ia (X_ref,ia^k X_ia^l - Y_ref,ia^k Y_ia^l),
2983!> k over a multiplet of the reference and l over the Davidson states inside it. Degenerate
2984!> eigenvectors are arbitrary within their multiplet, so the vectors are not compared one by one:
2985!> every singular value of O is 1 when the Davidson states lie in the span of the multiplet
2986!> \param ref_ens reference energies in Hartree, ascending
2987!> \param n_want Davidson states compared
2988!> \param n_ref reference states up to the end of the multiplet that contains state n_want
2989!> \param X_ref reference X, N_ov x n_ref
2990!> \param X_dav Davidson X, N_ov x n_want
2991!> \param min_overlap the smallest singular value found, 1 if every multiplet is reproduced
2992!> \param Y_ref reference Y for ABBA, with Y_dav
2993!> \param Y_dav Davidson Y for ABBA
2994! **************************************************************************************************
2995 SUBROUTINE min_multiplet_overlap(ref_ens, n_want, n_ref, X_ref, X_dav, min_overlap, Y_ref, Y_dav)
2996
2997 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: ref_ens
2998 INTEGER, INTENT(IN) :: n_want, n_ref
2999 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: x_ref, x_dav
3000 REAL(kind=dp), INTENT(OUT) :: min_overlap
3001 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
3002 OPTIONAL :: y_ref, y_dav
3003
3004 INTEGER :: info, lwork, mult_first, mult_last, &
3005 mult_last_dav
3006 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sing_vals, work
3007 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: overlap
3008 REAL(kind=dp), DIMENSION(1, 1) :: dummy
3009
3010 min_overlap = 1.0_dp
3011 mult_first = 1
3012 DO WHILE (mult_first <= n_want)
3013 mult_last = multiplet_end(ref_ens, mult_first, n_ref, deg_thresh)
3014 mult_last_dav = min(mult_last, n_want)
3015 ALLOCATE (overlap(mult_last - mult_first + 1, mult_last_dav - mult_first + 1), sing_vals(mult_last_dav - mult_first + 1))
3016 overlap(:, :) = matmul(transpose(x_ref(:, mult_first:mult_last)), x_dav(:, mult_first:mult_last_dav))
3017 IF (PRESENT(y_ref)) THEN
3018 overlap(:, :) = overlap(:, :) - matmul(transpose(y_ref(:, mult_first:mult_last)), y_dav(:, mult_first:mult_last_dav))
3019 END IF
3020 ! DGESVD wants lwork >= MAX(3 n + m, 5 n) for the m x n overlap with n <= m; 5 m + 10 covers both
3021 lwork = 5*(mult_last - mult_first + 1) + 10
3022 ALLOCATE (work(lwork))
3023 CALL dgesvd('N', 'N', mult_last - mult_first + 1, mult_last_dav - mult_first + 1, overlap, mult_last - mult_first + 1, sing_vals, &
3024 dummy, 1, dummy, 1, work, lwork, info)
3025 IF (info /= 0) cpabort("SVD failed in the BSE Davidson check")
3026 min_overlap = min(min_overlap, minval(sing_vals))
3027 DEALLOCATE (overlap, sing_vals, work)
3028 mult_first = mult_last + 1
3029 END DO
3030
3031 END SUBROUTINE min_multiplet_overlap
3032
3033! **************************************************************************************************
3034!> \brief The two lines of a reference check
3035!> \param dev_E largest energy deviation in Hartree
3036!> \param min_overlap smallest singular value from min_multiplet_overlap
3037!> \param unit_nr output unit, positive on the writing rank only
3038! **************************************************************************************************
3039 SUBROUTINE print_refcheck(dev_E, min_overlap, unit_nr)
3040
3041 REAL(kind=dp), INTENT(IN) :: dev_e, min_overlap
3042 INTEGER, INTENT(IN) :: unit_nr
3043
3044 IF (unit_nr > 0) THEN
3045 WRITE (unit_nr, '(T2,A10,T13,A,T63,ES18.6)') 'BSE|DEBUG|', &
3046 'Max energy deviation vs full diagonalization (eV)', dev_e*evolt
3047 WRITE (unit_nr, '(T2,A10,T13,A,T63,F18.12)') 'BSE|DEBUG|', &
3048 'Min subspace overlap vs full diagonalization', min_overlap
3049 END IF
3050
3051 END SUBROUTINE print_refcheck
3052
3053END MODULE bse_davidson
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::olsen1988
integer, save, public olsen1988
Definition bibliography.F:36
bibliography::fukaya2014
integer, save, public fukaya2014
Definition bibliography.F:36
bibliography::davidson1975
integer, save, public davidson1975
Definition bibliography.F:36
bibliography::bai2012
integer, save, public bai2012
Definition bibliography.F:36
bibliography::stathopoulos1998
integer, save, public stathopoulos1998
Definition bibliography.F:36
bibliography::stratmann1998
integer, save, public stratmann1998
Definition bibliography.F:36
bibliography::vecharynski2017
integer, save, public vecharynski2017
Definition bibliography.F:36
bse_davidson
Block Davidson solvers for the lowest excitations of the Bethe-Salpeter equation on top of the matrix...
Definition bse_davidson.F:19
bse_davidson::bse_davidson_abba_os
subroutine, public bse_davidson_abba_os(mv_env, bse_env, unit_nr, exc_ens, fm_x, fm_y, abba_status, ab_margin, fm_x_tda)
Block Davidson with thick restart for the NUM_EXC_EN lowest solutions of the full problem in the pair...
Definition bse_davidson.F:610
bse_davidson::abba_ok
integer, parameter, public abba_ok
Definition bse_davidson.F:80
bse_davidson::abba_indefinite
integer, parameter, public abba_indefinite
Definition bse_davidson.F:80
bse_davidson::bse_davidson_abba_mk
subroutine, public bse_davidson_abba_mk(mv_env, bse_env, unit_nr, exc_ens, fm_x, fm_y, abba_status, ab_margin, fm_x_tda)
Block Davidson with thick restart for the NUM_EXC_EN lowest solutions of the full problem,...
Definition bse_davidson.F:320
bse_davidson::bse_davidson_tda
subroutine, public bse_davidson_tda(mv_env, bse_env, unit_nr, exc_ens, fm_x)
Block Davidson with thick restart for the NUM_EXC_EN lowest solutions of sum_jb A_ia,...
Definition bse_davidson.F:105
bse_davidson::bse_davidson_refcheck
subroutine, public bse_davidson_refcheck(fm_a_explicit, exc_ens, fm_x, bse_env, unit_nr, fm_b_explicit, fm_y)
Debug check of a Davidson result against the full diagonalization of the explicit matrices: energies,...
Definition bse_davidson.F:2883
bse_full_diag
Routines for the full diagonalization of GW + Bethe-Salpeter for computing electronic excitations.
Definition bse_full_diag.F:14
bse_full_diag::create_hermitian_form_of_abba
subroutine, public create_hermitian_form_of_abba(fm_a, fm_b, fm_c, fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b, unit_nr, bse_env, diag_est)
Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem (cf....
Definition bse_full_diag.F:508
bse_matvec
Matrix-free application of the BSE matrices A and B to trial vectors from RI slabs that are sliced al...
Definition bse_matvec.F:14
bse_matvec::bse_matvec_subblock
subroutine, public bse_matvec_subblock(mv_env, ia_list, a_sub, b_sub)
Exact A (and B) on a list of transitions, replicated on every rank, A_kl = δ_kl (ε_a-ε_i) + α sum_P B...
Definition bse_matvec.F:524
bse_matvec::bse_matvec_apply
subroutine, public bse_matvec_apply(mv_env, fm_z, first_col, ncol, fm_az, fm_bz, first_col_bz)
Applies A (and B) to ncol trial vectors without forming an N_ov x N_ov object, (A Z)_ia = (ε_a-ε_i) Z...
Definition bse_matvec.F:319
bse_matvec::bse_matvec_diagonal
subroutine, public bse_matvec_diagonal(mv_env, precond_kind, diag)
Diagonal used by the Davidson correction, either ε_a-ε_i or the full diagonal A_ia,...
Definition bse_matvec.F:481
bse_matvec::mem_fraction
real(kind=dp), parameter, public mem_fraction
Definition bse_matvec.F:47
bse_matvec::bse_matvec_vector_struct
subroutine, public bse_matvec_vector_struct(mv_env, ncol_global, fm_struct)
Matrix structure of a block of trial vectors: rows ia distributed, all columns local.
Definition bse_matvec.F:288
bse_types
The BSE environment: the settings of the &BSE section and the state a GW path prepares for the solver...
Definition bse_types.F:14
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_blacs_env::cp_blacs_env_release
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
Definition cp_blacs_env.F:274
cp_blacs_env::cp_blacs_env_create
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Definition cp_blacs_env.F:123
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:1690
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:1532
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:170
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::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:262
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:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition cp_fm_types.F:750
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:1455
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:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:925
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::bse_memcheck_warn
integer, parameter, public bse_memcheck_warn
Definition input_constants.F:1422
input_constants::bse_memcheck_off
integer, parameter, public bse_memcheck_off
Definition input_constants.F:1422
input_constants::bse_iter_res_cond
integer, parameter, public bse_iter_res_cond
Definition input_constants.F:1408
input_constants::bse_memcheck_clamp
integer, parameter, public bse_memcheck_clamp
Definition input_constants.F:1422
input_constants::bse_iter_and_cond
integer, parameter, public bse_iter_and_cond
Definition input_constants.F:1408
input_constants::bse_iter_en_cond
integer, parameter, public bse_iter_en_cond
Definition input_constants.F:1408
input_constants::bse_memcheck_abort
integer, parameter, public bse_memcheck_abort
Definition input_constants.F:1422
input_constants::bse_iter_or_cond
integer, parameter, public bse_iter_or_cond
Definition input_constants.F:1408
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:141
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition mathlib.F:381
mathlib::diag
subroutine, public diag(n, a, d, v)
Diagonalize matrix a. The eigenvalues are returned in vector d and the eigenvectors are returned in m...
Definition mathlib.F:1633
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
message_passing::mp_mem_avail_per_rank_gb
subroutine, public mp_mem_avail_per_rank_gb(comm, mem_avail_gb)
Memory that is currently free on the node, per MPI rank of that node, in GB.
Definition message_passing.F:2734
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
util
All kind of helpful little routines.
Definition util.F:14
bse_matvec::bse_matvec_env_type
RI slabs of one spin channel, every rank holding n_ri_loc whole RI slices.
Definition bse_matvec.F:67
bse_types::bse_env_type
Settings of the &BSE section (read_bse_section, re-read by prepare_bse_env, normalised in place by ad...
Definition bse_types.F:41
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
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:114
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768