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