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 Matrix-free application of the BSE matrices A and B to trial vectors from RI slabs that
10 : !> are sliced along the RI index over all MPI ranks
11 : !> \par History
12 : !> 09.2026 created [Maximilian Graml]
13 : ! **************************************************************************************************
14 : MODULE bse_matvec
15 :
16 : USE bse_util, ONLY: ia_of_occ_virt,&
17 : occ_of_ia,&
18 : virt_of_ia
19 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
20 : cp_blacs_env_release,&
21 : cp_blacs_env_type
22 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
23 : cp_fm_struct_release,&
24 : cp_fm_struct_type
25 : USE cp_fm_types, ONLY: cp_fm_create,&
26 : cp_fm_get_info,&
27 : cp_fm_get_submatrix,&
28 : cp_fm_release,&
29 : cp_fm_set_all,&
30 : cp_fm_to_fm_submat_general,&
31 : cp_fm_type
32 : USE input_constants, ONLY: bse_precond_full_diag
33 : USE kinds, ONLY: dp
34 : USE message_passing, ONLY: mp_mem_avail_per_rank_GB,&
35 : mp_mem_used_per_rank_GB,&
36 : mp_para_env_type
37 : #include "./base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 :
41 : PRIVATE
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_matvec'
44 :
45 : ! the detected MEMORY_PER_PROC claims this share of the free memory on top of the occupied memory,
46 : ! the share resolve_memory_per_proc of the GW claims; the probe counts reclaimable cache
47 : REAL(KIND=dp), PARAMETER :: detected_memory_fraction = 0.8_dp, &
48 : fallback_memory_per_proc_GB = 2.0_dp
49 :
50 : PUBLIC :: bse_matvec_env_type, bse_matvec_create, bse_matvec_release, bse_matvec_apply, &
51 : bse_memory_left_GB, bse_resolve_memory_per_proc, &
52 : bse_matvec_diagonal, bse_matvec_subblock, bse_matvec_vector_struct, bse_matvec_selfcheck
53 :
54 : ! **************************************************************************************************
55 : !> \brief RI slabs of one spin channel, every rank holding n_ri_loc whole RI slices
56 : !> \param alpha prefactor of the exchange term (2 singlet, 0 triplet)
57 : !> \param w_fac prefactor of the screened term (0 switches it off)
58 : !> \param eps_diff ε_a - ε_i at the compound index ia = (i-1)*virt + a of bse_util's ia_of_occ_virt
59 : !> \param B_ia B^P_ia at (a, i, P)
60 : !> \param B_bar_ia \bar{B}^P_ia at (a, i, P), only for do_abba
61 : !> \param B_bar_ij \bar{B}^P_ij at (j, i, P)
62 : !> \param B_ab B^P_ab at (b, a, P)
63 : !> \param row_count rows of a trial vector owned by each rank, in rank order
64 : !> \param row_displ first row of each rank minus one, so that row_count and row_displ describe
65 : !> the contiguous ia range per rank that allgatherv and sum_scatter need
66 : !> \param blacs_env npe x 1 process grid of the slabs and of all trial vectors
67 : !> \param block_cols trial vectors per pass in bse_matvec_apply, -1 to size it from the memory left
68 : !> \param memory_per_proc_gb the resolved MEMORY_PER_PROC of the BSE, which block_cols -1 sizes the pass from
69 : ! **************************************************************************************************
70 : TYPE bse_matvec_env_type
71 : INTEGER :: homo = 0, virt = 0, n_ov = 0, &
72 : n_ri_loc = 0, block_cols = -1
73 : REAL(KIND=dp) :: alpha = 0.0_dp, w_fac = 0.0_dp, &
74 : memory_per_proc_gb = -1.0_dp
75 : LOGICAL :: do_abba = .FALSE.
76 : INTEGER, ALLOCATABLE, DIMENSION(:) :: row_count, row_displ
77 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eps_diff
78 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: B_ia, B_bar_ia, B_bar_ij, B_ab
79 : TYPE(cp_blacs_env_type), POINTER :: blacs_env => NULL()
80 : TYPE(mp_para_env_type), POINTER :: para_env => NULL()
81 : END TYPE bse_matvec_env_type
82 :
83 : CONTAINS
84 :
85 : ! **************************************************************************************************
86 : !> \brief MEMORY_PER_PROC of the BSE as a number: the keyword's value, or for -1 the memory the rank
87 : !> holds plus the share detected_memory_fraction of the free memory, read once at the first
88 : !> consumer; when the free memory cannot be read, the GW's fallback of 2 GB with a warning
89 : !> \param memory_per_proc_gb the keyword's value on entry, the memory per MPI rank for the BSE on exit
90 : !> \param para_env the ranks that share the node's memory
91 : !> \param unit_nr output unit, positive on the writing rank only
92 : ! **************************************************************************************************
93 116 : SUBROUTINE bse_resolve_memory_per_proc(memory_per_proc_gb, para_env, unit_nr)
94 : REAL(KIND=dp), INTENT(INOUT) :: memory_per_proc_gb
95 : TYPE(mp_para_env_type), POINTER :: para_env
96 : INTEGER, INTENT(IN) :: unit_nr
97 :
98 : REAL(KIND=dp) :: mem_free_GB, mem_occ_GB
99 :
100 116 : IF (memory_per_proc_gb >= 0.0_dp) RETURN
101 116 : CALL mp_mem_avail_per_rank_GB(para_env, mem_free_GB)
102 116 : CALL mp_mem_used_per_rank_GB(para_env, mem_occ_GB)
103 116 : IF (mem_free_GB > 0.0_dp) THEN
104 116 : memory_per_proc_gb = mem_occ_GB + detected_memory_fraction*mem_free_GB
105 116 : IF (unit_nr > 0) THEN
106 58 : WRITE (unit_nr, '(T2,A4,T7,A,T67,F14.3)') 'BSE|', 'Memory per MPI rank for the BSE, detected (GB)', &
107 116 : memory_per_proc_gb
108 : END IF
109 : ELSE
110 0 : memory_per_proc_gb = fallback_memory_per_proc_GB
111 : CALL cp_warn(__LOCATION__, &
112 : "BSE: the free memory per MPI rank could not be detected; MEMORY_PER_PROC falls back "// &
113 0 : "to 2 GB. Set MEMORY_PER_PROC in the BSE section for a run that fits the node.")
114 : END IF
115 :
116 : END SUBROUTINE bse_resolve_memory_per_proc
117 :
118 : ! **************************************************************************************************
119 : !> \brief What MEMORY_PER_PROC leaves to a step: the resolved memory per MPI rank minus the memory the
120 : !> rank holds now, never negative; both probes reduce over the ranks, so every rank gets the same
121 : !> \param memory_per_proc_gb the resolved MEMORY_PER_PROC
122 : !> \param para_env the ranks that share the node's memory
123 : !> \return memory left per MPI rank in GB
124 : ! **************************************************************************************************
125 70 : FUNCTION bse_memory_left_GB(memory_per_proc_gb, para_env) RESULT(mem_left_GB)
126 : REAL(KIND=dp), INTENT(IN) :: memory_per_proc_gb
127 : TYPE(mp_para_env_type), POINTER :: para_env
128 : REAL(KIND=dp) :: mem_left_GB
129 :
130 : REAL(KIND=dp) :: mem_occ_GB
131 :
132 : ! a consumer before the resolution would size itself from the keyword's -1
133 70 : CPASSERT(memory_per_proc_gb >= 0.0_dp)
134 70 : CALL mp_mem_used_per_rank_GB(para_env, mem_occ_GB)
135 70 : mem_left_GB = MAX(memory_per_proc_gb - mem_occ_GB, 0.0_dp)
136 :
137 70 : END FUNCTION bse_memory_left_GB
138 :
139 : ! **************************************************************************************************
140 : !> \brief Moves the RI slabs onto an npe x 1 process grid, such that every rank owns whole RI
141 : !> slices, and stores them as local 3-index arrays: the contraction over each P in
142 : !> bse_matvec_apply is then a local DGEMM and only the final sum over P crosses the ranks
143 : !> \param mv_env the environment created here
144 : !> \param fm_S_ia B^P_ia, N_RI x homo*virt
145 : !> \param fm_S_bar_ij \bar{B}^P_ij, the slab that enters W_ij,ab, N_RI x homo*homo
146 : !> \param fm_S_ab B^P_ab, N_RI x virt*virt
147 : !> \param fm_S_bar_ia \bar{B}^P_ia, the slab that enters W_ib,aj, N_RI x homo*virt
148 : !> \param eps_reduced quasiparticle energies of the homo+virt active levels
149 : !> \param homo occupied levels of the active window
150 : !> \param virt virtual levels of the active window
151 : !> \param alpha prefactor of the exchange term, 2 singlet, 0 triplet
152 : !> \param w_fac prefactor of the screened term, 0 switches it off
153 : !> \param do_abba slice \bar{B}^P_ia as well, for the application of B
154 : !> \param unit_nr output unit, positive on the writing rank only
155 : !> \param block_cols trial vectors per pass of bse_matvec_apply; absent or -1 sizes it from the memory left
156 : !> \param memory_per_proc_gb the resolved MEMORY_PER_PROC of the BSE, what block_cols -1 sizes the pass from
157 : ! **************************************************************************************************
158 8 : SUBROUTINE bse_matvec_create(mv_env, fm_S_ia, fm_S_bar_ij, fm_S_ab, fm_S_bar_ia, eps_reduced, &
159 : homo, virt, alpha, w_fac, do_abba, unit_nr, block_cols, memory_per_proc_gb)
160 :
161 : TYPE(bse_matvec_env_type), INTENT(OUT) :: mv_env
162 : TYPE(cp_fm_type), INTENT(IN) :: fm_S_ia, fm_S_bar_ij, fm_S_ab, &
163 : fm_S_bar_ia
164 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: eps_reduced
165 : INTEGER, INTENT(IN) :: homo, virt
166 : REAL(KIND=dp), INTENT(IN) :: alpha, w_fac
167 : LOGICAL, INTENT(IN) :: do_abba
168 : INTEGER, INTENT(IN) :: unit_nr
169 : INTEGER, INTENT(IN), OPTIONAL :: block_cols
170 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: memory_per_proc_gb
171 :
172 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_matvec_create'
173 :
174 : CHARACTER(LEN=12) :: n_idle_str, n_ri_str, npe_str
175 : INTEGER :: a, handle, i, n_ri, n_ri_slab, &
176 : n_row_block
177 : REAL(KIND=dp) :: mem_slabs_GB, n_slab_entries
178 :
179 8 : CALL timeset(routineN, handle)
180 :
181 8 : mv_env%homo = homo
182 8 : mv_env%virt = virt
183 8 : IF (PRESENT(block_cols)) mv_env%block_cols = block_cols
184 8 : IF (PRESENT(memory_per_proc_gb)) mv_env%memory_per_proc_gb = memory_per_proc_gb
185 8 : mv_env%n_ov = homo*virt
186 8 : mv_env%alpha = alpha
187 8 : mv_env%w_fac = w_fac
188 8 : mv_env%do_abba = do_abba
189 8 : mv_env%para_env => fm_S_ia%matrix_struct%para_env
190 :
191 : CALL cp_blacs_env_create(mv_env%blacs_env, mv_env%para_env, &
192 24 : grid_2d=[mv_env%para_env%num_pe, 1])
193 :
194 : ! one contiguous ia range per rank, so that a trial vector is gathered with allgatherv and
195 : ! its image scattered back with sum_scatter, rather than replicated by two allreduces
196 8 : n_row_block = contiguous_row_block(mv_env%n_ov, mv_env%para_env%num_pe)
197 40 : ALLOCATE (mv_env%row_count(mv_env%para_env%num_pe), mv_env%row_displ(mv_env%para_env%num_pe))
198 24 : DO i = 1, mv_env%para_env%num_pe
199 16 : mv_env%row_displ(i) = MIN(mv_env%n_ov, (i - 1)*n_row_block)
200 24 : mv_env%row_count(i) = MIN(mv_env%n_ov, i*n_row_block) - mv_env%row_displ(i)
201 : END DO
202 :
203 : ! eps_diff_ia = ε_a - ε_i at the compound index ia = (i-1)*virt + a, the row index of every trial vector
204 24 : ALLOCATE (mv_env%eps_diff(mv_env%n_ov))
205 40 : DO i = 1, homo
206 424 : DO a = 1, virt
207 416 : mv_env%eps_diff(ia_of_occ_virt(i, a, virt)) = eps_reduced(homo + a) - eps_reduced(i)
208 : END DO
209 : END DO
210 :
211 : ! every slab is sliced on the same grid, so every rank owns the same RI indices of each slab
212 8 : CALL slice_slab(mv_env, fm_S_ia, virt, homo, mv_env%B_ia, mv_env%n_ri_loc)
213 8 : IF (w_fac /= 0.0_dp) THEN
214 8 : CALL slice_slab(mv_env, fm_S_bar_ij, homo, homo, mv_env%B_bar_ij, n_ri_slab)
215 8 : CPASSERT(n_ri_slab == mv_env%n_ri_loc)
216 8 : CALL slice_slab(mv_env, fm_S_ab, virt, virt, mv_env%B_ab, n_ri_slab)
217 8 : CPASSERT(n_ri_slab == mv_env%n_ri_loc)
218 8 : IF (do_abba) THEN
219 4 : CALL slice_slab(mv_env, fm_S_bar_ia, virt, homo, mv_env%B_bar_ia, n_ri_slab)
220 4 : CPASSERT(n_ri_slab == mv_env%n_ri_loc)
221 : END IF
222 : END IF
223 :
224 8 : CALL cp_fm_get_info(fm_S_ia, nrow_global=n_ri)
225 8 : IF (mv_env%para_env%num_pe > n_ri .AND. unit_nr > 0) THEN
226 0 : WRITE (npe_str, '(I0)') mv_env%para_env%num_pe
227 0 : WRITE (n_ri_str, '(I0)') n_ri
228 0 : WRITE (n_idle_str, '(I0)') mv_env%para_env%num_pe - n_ri
229 : CALL cp_warn(__LOCATION__, &
230 : "BSE iterative solver: more MPI ranks ("//TRIM(npe_str)// &
231 : ") than RI basis functions ("//TRIM(n_ri_str)//"); "// &
232 0 : TRIM(n_idle_str)//" ranks hold no RI slice and idle in the kernel application.")
233 : END IF
234 :
235 8 : IF (unit_nr > 0) THEN
236 : ! pair entries of the sliced slabs per RI index
237 4 : n_slab_entries = REAL(homo*virt, dp)
238 4 : IF (w_fac /= 0.0_dp) THEN
239 4 : n_slab_entries = n_slab_entries + REAL(homo, dp)**2 + REAL(virt, dp)**2
240 4 : IF (do_abba) n_slab_entries = n_slab_entries + REAL(homo*virt, dp)
241 : END IF
242 4 : mem_slabs_GB = n_slab_entries*REAL(mv_env%n_ri_loc, dp)*8.0E-9_dp
243 4 : WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', &
244 8 : 'Max. number of RI functions per MPI rank', mv_env%n_ri_loc
245 4 : WRITE (unit_nr, '(T2,A4,T7,A,T67,F14.3)') 'BSE|', &
246 8 : 'Memory of the RI slabs per MPI rank (GB)', mem_slabs_GB
247 : END IF
248 :
249 8 : CALL timestop(handle)
250 :
251 16 : END SUBROUTINE bse_matvec_create
252 :
253 : ! **************************************************************************************************
254 : !> \brief Rows of a trial vector per rank, the smallest block that leaves no rank beyond the last
255 : !> \param n_ov number of transitions, the global row count
256 : !> \param npe ranks of the npe x 1 grid
257 : !> \return rows of every full block; the last rank takes what remains, possibly none
258 : ! **************************************************************************************************
259 36 : PURE FUNCTION contiguous_row_block(n_ov, npe) RESULT(n_row_block)
260 :
261 : INTEGER, INTENT(IN) :: n_ov, npe
262 : INTEGER :: n_row_block
263 :
264 36 : n_row_block = (n_ov + npe - 1)/npe
265 :
266 36 : END FUNCTION contiguous_row_block
267 :
268 : ! **************************************************************************************************
269 : !> \brief Redistributes one N_RI x n_fast*n_slow slab to whole RI rows per rank and copies the owned
270 : !> rows into a contiguous (n_fast, n_slow, n_ri_loc) array, slab_loc(x, y, p) = slab(P_p, (y-1)*n_fast + x)
271 : !> for the p-th owned RI index P_p, so that every slice slab_loc(:, :, p) is a n_fast x n_slow DGEMM operand
272 : !> \param mv_env process grid and communicator of the slices
273 : !> \param fm_slab the slab on the grid of the caller, N_RI x n_fast*n_slow
274 : !> \param n_fast fast index of the pair index of the slab
275 : !> \param n_slow slow index of the pair index of the slab
276 : !> \param slab_loc the owned slices, allocated here
277 : !> \param n_ri_loc RI indices owned by this rank, the third extent of slab_loc
278 : ! **************************************************************************************************
279 28 : SUBROUTINE slice_slab(mv_env, fm_slab, n_fast, n_slow, slab_loc, n_ri_loc)
280 :
281 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
282 : TYPE(cp_fm_type), INTENT(IN) :: fm_slab
283 : INTEGER, INTENT(IN) :: n_fast, n_slow
284 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
285 : INTENT(OUT) :: slab_loc
286 : INTEGER, INTENT(OUT) :: n_ri_loc
287 :
288 : CHARACTER(LEN=*), PARAMETER :: routineN = 'slice_slab'
289 :
290 : INTEGER :: handle, n_ri, npair, p
291 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
292 : TYPE(cp_fm_type) :: fm_sliced
293 :
294 28 : CALL timeset(routineN, handle)
295 :
296 : ! row count from the struct: the slabs may carry fewer rows than the full RI basis
297 28 : CALL cp_fm_get_info(fm_slab, nrow_global=n_ri, ncol_global=npair)
298 28 : CPASSERT(npair == n_fast*n_slow)
299 :
300 : ! one-row blocks on the npe x 1 grid spread the RI rows cyclically over the ranks; the copy
301 : ! below turns each strided local_data row into a contiguous (n_fast, n_slow) slice
302 28 : NULLIFY (fm_struct)
303 : CALL cp_fm_struct_create(fm_struct, para_env=mv_env%para_env, context=mv_env%blacs_env, &
304 : nrow_global=n_ri, ncol_global=npair, nrow_block=1, &
305 28 : force_block=.TRUE.)
306 28 : CALL cp_fm_create(fm_sliced, fm_struct, name="fm_slab_ri_sliced")
307 28 : CALL cp_fm_struct_release(fm_struct)
308 : CALL cp_fm_to_fm_submat_general(fm_slab, fm_sliced, n_ri, npair, 1, 1, 1, 1, &
309 28 : fm_slab%matrix_struct%context)
310 :
311 28 : CALL cp_fm_get_info(fm_sliced, nrow_local=n_ri_loc)
312 140 : ALLOCATE (slab_loc(n_fast, n_slow, n_ri_loc))
313 1190 : DO p = 1, n_ri_loc
314 3514 : slab_loc(:, :, p) = RESHAPE(fm_sliced%local_data(p, 1:npair), [n_fast, n_slow])
315 : END DO
316 28 : CALL cp_fm_release(fm_sliced)
317 :
318 28 : CALL timestop(handle)
319 :
320 84 : END SUBROUTINE slice_slab
321 :
322 : ! **************************************************************************************************
323 : !> \brief Frees the sliced slabs, the transition energies and the process grid of mv_env
324 : !> \param mv_env the environment released
325 : ! **************************************************************************************************
326 8 : SUBROUTINE bse_matvec_release(mv_env)
327 :
328 : TYPE(bse_matvec_env_type), INTENT(INOUT) :: mv_env
329 :
330 8 : IF (ALLOCATED(mv_env%row_count)) DEALLOCATE (mv_env%row_count)
331 8 : IF (ALLOCATED(mv_env%row_displ)) DEALLOCATE (mv_env%row_displ)
332 8 : IF (ALLOCATED(mv_env%eps_diff)) DEALLOCATE (mv_env%eps_diff)
333 8 : IF (ALLOCATED(mv_env%B_ia)) DEALLOCATE (mv_env%B_ia)
334 8 : IF (ALLOCATED(mv_env%B_bar_ia)) DEALLOCATE (mv_env%B_bar_ia)
335 8 : IF (ALLOCATED(mv_env%B_bar_ij)) DEALLOCATE (mv_env%B_bar_ij)
336 8 : IF (ALLOCATED(mv_env%B_ab)) DEALLOCATE (mv_env%B_ab)
337 8 : IF (ASSOCIATED(mv_env%blacs_env)) CALL cp_blacs_env_release(mv_env%blacs_env)
338 8 : NULLIFY (mv_env%para_env)
339 :
340 8 : END SUBROUTINE bse_matvec_release
341 :
342 : ! **************************************************************************************************
343 : !> \brief Matrix structure of a block of trial vectors: rows ia distributed, all columns local
344 : !> \param mv_env process grid and communicator of the trial vectors
345 : !> \param ncol_global number of trial vectors
346 : !> \param fm_struct created here, released by the caller
347 : ! **************************************************************************************************
348 28 : SUBROUTINE bse_matvec_vector_struct(mv_env, ncol_global, fm_struct)
349 :
350 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
351 : INTEGER, INTENT(IN) :: ncol_global
352 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
353 :
354 28 : NULLIFY (fm_struct)
355 : ! the row block matches mv_env%row_count, which bse_matvec_apply's collectives rely on
356 : CALL cp_fm_struct_create(fm_struct, para_env=mv_env%para_env, context=mv_env%blacs_env, &
357 : nrow_global=mv_env%n_ov, ncol_global=ncol_global, &
358 : nrow_block=contiguous_row_block(mv_env%n_ov, mv_env%para_env%num_pe), &
359 28 : force_block=.TRUE.)
360 :
361 28 : END SUBROUTINE bse_matvec_vector_struct
362 :
363 : ! **************************************************************************************************
364 : !> \brief Applies A (and B) to ncol trial vectors without forming an N_ov x N_ov object,
365 : !> (A Z)_ia = (ε_a-ε_i) Z_ia + α sum_P B^P_ia sum_jb B^P_jb Z_jb - w sum_Pjb \bar{B}^P_ij B^P_ab Z_jb
366 : !> (B Z)_ia = α sum_P B^P_ia sum_jb B^P_jb Z_jb - w sum_Pjb \bar{B}^P_ib B^P_ja Z_jb
367 : !> α: exchange prefactor (2 singlet, 0 triplet), w: prefactor of the screened term,
368 : !> \bar{B}: slab bound by the caller of bse_matvec_create (B itself for TDHF and ALPHA
369 : !> screening, sum_Q [1+Q(0)]^-1_PQ B^Q otherwise); the sums over P run over the local RI
370 : !> slices and are completed by a sum over the ranks
371 : !> \param mv_env slabs, prefactors and transition energies
372 : !> \param fm_Z trial vectors, columns first_col .. first_col+ncol-1 are read
373 : !> \param first_col first trial vector read, and first column written
374 : !> \param ncol number of trial vectors
375 : !> \param fm_AZ receives A Z in the same columns
376 : !> \param fm_BZ receives B Z in the same columns, or from first_col_BZ on
377 : !> \param first_col_BZ first column of fm_BZ written, first_col by default
378 : ! **************************************************************************************************
379 84 : SUBROUTINE bse_matvec_apply(mv_env, fm_Z, first_col, ncol, fm_AZ, fm_BZ, first_col_BZ)
380 :
381 : TYPE(bse_matvec_env_type), INTENT(IN), TARGET :: mv_env
382 : TYPE(cp_fm_type), INTENT(IN) :: fm_Z
383 : INTEGER, INTENT(IN) :: first_col, ncol
384 : TYPE(cp_fm_type), INTENT(IN) :: fm_AZ
385 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_BZ
386 : INTEGER, INTENT(IN), OPTIONAL :: first_col_BZ
387 :
388 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_matvec_apply'
389 :
390 : INTEGER :: a, col, col_shift_B, first_row, handle, handle2, homo, i, ia, iloc, k, last_row, &
391 : me, n_outputs, n_small, n_work_bufs, nb, nb_max, nrow_local, nvec, offset, p, peer, rows, &
392 : virt
393 84 : INTEGER, DIMENSION(:), POINTER :: row_indices
394 : LOGICAL :: do_B, do_screened, has_exchange
395 : REAL(KIND=dp) :: mem_avail_GB
396 84 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: result_buf, work_buf
397 84 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: T_Pk
398 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
399 84 : POINTER :: B_ia_2d, local_result, packed, RA, &
400 84 : RA_virt, RB_occ, RB_virt, T_occ, &
401 84 : T_small_cols, T_small_rows, T_virt, Z, &
402 84 : Z_occ, Z_virt
403 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
404 84 : POINTER :: RA_screen, RB_screen, Z_screen
405 168 : INTEGER :: displs(SIZE(mv_env%row_displ)), &
406 168 : counts(SIZE(mv_env%row_count))
407 :
408 84 : CALL timeset(routineN, handle)
409 :
410 84 : homo = mv_env%homo
411 84 : virt = mv_env%virt
412 84 : n_small = MIN(homo, virt)
413 84 : do_B = PRESENT(fm_BZ)
414 84 : do_screened = mv_env%w_fac /= 0.0_dp .AND. mv_env%n_ri_loc > 0
415 84 : has_exchange = mv_env%alpha /= 0.0_dp .AND. mv_env%n_ri_loc > 0
416 84 : IF (do_B) THEN
417 : ! B Z needs the \bar{B}^P_ia slab, which is sliced for do_abba only
418 66 : CPASSERT(mv_env%do_abba)
419 : END IF
420 84 : col_shift_B = 0
421 84 : IF (PRESENT(first_col_BZ)) col_shift_B = first_col_BZ - first_col
422 :
423 84 : CALL cp_fm_get_info(fm_Z, nrow_local=nrow_local, row_indices=row_indices)
424 : ! the collectives below read the rows of one rank as one range, as bse_matvec_vector_struct
425 : ! lays them out; a vector built elsewhere would scatter into the wrong rows
426 84 : me = mv_env%para_env%mepos + 1
427 84 : CPASSERT(nrow_local == mv_env%row_count(me))
428 84 : CPASSERT(nrow_local == 0 .OR. row_indices(1) == mv_env%row_displ(me) + 1)
429 :
430 : ! columns per pass from three/four replicated blocks: trial vectors, scratch and A (and B)
431 84 : n_outputs = 1
432 84 : IF (do_B) n_outputs = 2
433 84 : n_work_bufs = 2 + n_outputs
434 84 : IF (mv_env%block_cols > 0) THEN
435 : ! pinned by input: the replicated work arrays are then a known n_work_bufs*n_ov*nb doubles
436 84 : nb_max = MIN(ncol, mv_env%block_cols)
437 : ELSE
438 : ! the pass takes what MEMORY_PER_PROC leaves after the solver's arrays
439 0 : mem_avail_GB = bse_memory_left_GB(mv_env%memory_per_proc_gb, mv_env%para_env)
440 : nb_max = INT(MIN(REAL(ncol, dp), &
441 0 : mem_avail_GB*1.0E9_dp/(8.0_dp*REAL(n_work_bufs, dp)*REAL(mv_env%n_ov, dp))))
442 : END IF
443 84 : nb_max = MAX(nb_max, 1)
444 : ! The paired work buffer and all collective counts/displacements use default INTEGER.
445 84 : nb_max = MIN(nb_max, (HUGE(0)/2)/mv_env%n_ov)
446 84 : CPASSERT(nb_max >= 1)
447 :
448 : ! Keep the same three/four full buffers, reused by communication, exchange and screening.
449 : ! Canonical Z and the exchange image start in the first slots of work_buf and result_buf.
450 420 : ALLOCATE (work_buf(2*mv_env%n_ov*nb_max), result_buf(n_outputs*mv_env%n_ov*nb_max))
451 336 : ALLOCATE (T_Pk(mv_env%n_ri_loc, nb_max))
452 84 : B_ia_2d(1:mv_env%n_ov, 1:mv_env%n_ri_loc) => mv_env%B_ia
453 :
454 168 : DO col = first_col, first_col + ncol - 1, nb_max
455 84 : nb = MIN(nb_max, first_col + ncol - col)
456 84 : nvec = mv_env%n_ov*nb
457 : ! Canonical ia-by-k views for the exchange block GEMMs.
458 84 : Z(1:mv_env%n_ov, 1:nb) => work_buf(1:nvec)
459 84 : RA(1:mv_env%n_ov, 1:nb) => result_buf(1:nvec)
460 :
461 : ! Z_ia,k on every rank, by design: each rank contracts its whole RI slices against all of Z (header of
462 : ! bse_matvec_create), so the trial vectors are the one replicated buffer, n_ov x nb per matvec at any rank
463 : ! count. Gather the entire block in rank-major order, then unpack the contiguous row ranges.
464 : ! result_buf is free until the exchange; the second work slot receives the gathered data.
465 84 : CALL timeset(routineN//"_gather", handle2)
466 348 : DO k = 1, nb
467 264 : offset = (k - 1)*nrow_local
468 6684 : result_buf(offset + 1:offset + nrow_local) = fm_Z%local_data(1:nrow_local, col + k - 1)
469 : END DO
470 252 : counts = mv_env%row_count*nb
471 252 : displs = mv_env%row_displ*nb
472 84 : CALL mv_env%para_env%allgatherv(result_buf(1:nrow_local*nb), work_buf(nvec + 1:2*nvec), counts, displs)
473 252 : DO peer = 1, mv_env%para_env%num_pe
474 168 : rows = mv_env%row_count(peer)
475 168 : first_row = mv_env%row_displ(peer) + 1
476 168 : last_row = mv_env%row_displ(peer) + rows
477 780 : DO k = 1, nb
478 528 : offset = nvec + displs(peer) + (k - 1)*rows
479 26040 : Z(first_row:last_row, k) = work_buf(offset + 1:offset + rows)
480 : END DO
481 : END DO
482 84 : CALL timestop(handle2)
483 :
484 : ! exchange, identical in A Z and B Z: R_ia,k = α sum_P B^P_ia t_Pk, t_Pk = sum_jb B^P_jb Z_jb,k
485 84 : CALL timeset(routineN//"_exchange", handle2)
486 84 : IF (has_exchange) THEN
487 : CALL DGEMM('T', 'N', mv_env%n_ri_loc, nb, mv_env%n_ov, 1.0_dp, B_ia_2d, mv_env%n_ov, Z, mv_env%n_ov, &
488 84 : 0.0_dp, T_Pk, mv_env%n_ri_loc)
489 : END IF
490 : ! Without screening, Z is no longer needed after forming T_Pk. Reuse its slot
491 : ! for the canonical exchange image, leaving result_buf free for the reduction.
492 84 : IF (.NOT. do_screened) RA(1:mv_env%n_ov, 1:nb) => work_buf(1:nvec)
493 84 : IF (has_exchange) THEN
494 : CALL DGEMM('N', 'N', mv_env%n_ov, nb, mv_env%n_ri_loc, mv_env%alpha, B_ia_2d, mv_env%n_ov, T_Pk, mv_env%n_ri_loc, &
495 84 : 0.0_dp, RA, mv_env%n_ov)
496 0 : ELSE IF (.NOT. do_screened) THEN
497 0 : RA(:, :) = 0.0_dp
498 : END IF
499 84 : CALL timestop(handle2)
500 :
501 84 : CALL timeset(routineN//"_pack", handle2)
502 84 : IF (do_screened) THEN
503 : ! Put the trial-vector index between the two orbital indices. One flat block then has
504 : ! both (i, (k a)) and ((i k), a) views, so neither contraction needs a per-vector GEMM.
505 : ! TDA keeps this Z in work slot 2; ABBA uses result slot 2, leaving both work slots for A/B.
506 84 : IF (do_B) THEN
507 66 : Z_screen(1:homo, 1:nb, 1:virt) => result_buf(nvec + 1:2*nvec)
508 : ELSE
509 18 : Z_screen(1:homo, 1:nb, 1:virt) => work_buf(nvec + 1:2*nvec)
510 : END IF
511 84 : Z_occ(1:homo, 1:nb*virt) => Z_screen
512 84 : Z_virt(1:homo*nb, 1:virt) => Z_screen
513 1092 : DO a = 1, virt
514 4260 : DO k = 1, nb
515 16848 : DO i = 1, homo
516 15840 : Z_screen(i, k, a) = Z(ia_of_occ_virt(i, a, virt), k)
517 : END DO
518 : END DO
519 : END DO
520 :
521 : ! Copy the shared exchange image into the two screened output layouts, then its original
522 : ! result slot is free for all RI-slice intermediates. These permutations happen once per pass.
523 84 : RA_screen(1:homo, 1:nb, 1:virt) => work_buf(1:nvec)
524 84 : RA_virt(1:homo*nb, 1:virt) => work_buf(1:nvec)
525 84 : IF (do_B) THEN
526 66 : RB_screen(1:virt, 1:nb, 1:homo) => work_buf(nvec + 1:2*nvec)
527 66 : RB_occ(1:virt, 1:nb*homo) => work_buf(nvec + 1:2*nvec)
528 66 : RB_virt(1:virt*nb, 1:homo) => work_buf(nvec + 1:2*nvec)
529 : END IF
530 84 : IF (has_exchange) THEN
531 1092 : DO a = 1, virt
532 4260 : DO k = 1, nb
533 16848 : DO i = 1, homo
534 12672 : ia = ia_of_occ_virt(i, a, virt)
535 12672 : RA_screen(i, k, a) = RA(ia, k)
536 15840 : IF (do_B) RB_screen(a, k, i) = RA(ia, k)
537 : END DO
538 : END DO
539 : END DO
540 : ELSE
541 0 : RA_virt(:, :) = 0.0_dp
542 0 : IF (do_B) RB_occ(:, :) = 0.0_dp
543 : END IF
544 84 : T_occ(1:homo, 1:nb*virt) => result_buf(1:nvec)
545 84 : T_virt(1:homo*nb, 1:virt) => result_buf(1:nvec)
546 : ! The B association chooses the smaller square intermediate: n_small^2*nb <= nvec.
547 84 : T_small_cols(1:n_small, 1:nb*n_small) => result_buf(1:n_small*n_small*nb)
548 84 : T_small_rows(1:n_small*nb, 1:n_small) => result_buf(1:n_small*n_small*nb)
549 : END IF
550 84 : CALL timestop(handle2)
551 :
552 84 : CALL timeset(routineN//"_W", handle2)
553 84 : IF (do_screened) THEN
554 3570 : DO p = 1, mv_env%n_ri_loc
555 : ! A: T_ikb = sum_j \bar{B}^P_ji Z_bjk, R^A_ika -= w sum_b T_ikb B^P_ba.
556 : ! Only the view of T changes between the occupied and virtual block GEMMs.
557 : CALL DGEMM('T', 'N', homo, nb*virt, homo, 1.0_dp, mv_env%B_bar_ij(:, :, p), homo, &
558 3486 : Z_occ, homo, 0.0_dp, T_occ, homo)
559 : CALL DGEMM('N', 'N', homo*nb, virt, virt, -mv_env%w_fac, T_virt, homo*nb, &
560 3486 : mv_env%B_ab(:, :, p), virt, 1.0_dp, RA_virt, homo*nb)
561 3570 : IF (do_B) THEN
562 : ! B: R^B_aik -= w sum_jb B^P_aj Z_bjk \bar{B}^P_bi.
563 : ! Associate the products as B (Z^T \bar{B}) for homo <= virt, otherwise (B Z^T) \bar{B}.
564 2739 : IF (homo <= virt) THEN
565 : CALL DGEMM('N', 'N', homo*nb, homo, virt, 1.0_dp, Z_virt, homo*nb, &
566 2739 : mv_env%B_bar_ia(:, :, p), virt, 0.0_dp, T_small_rows, homo*nb)
567 : CALL DGEMM('N', 'N', virt, nb*homo, homo, -mv_env%w_fac, mv_env%B_ia(:, :, p), virt, &
568 2739 : T_small_cols, homo, 1.0_dp, RB_occ, virt)
569 : ELSE
570 : CALL DGEMM('N', 'N', virt, nb*virt, homo, 1.0_dp, mv_env%B_ia(:, :, p), virt, &
571 0 : Z_occ, homo, 0.0_dp, T_small_cols, virt)
572 : CALL DGEMM('N', 'N', virt*nb, homo, virt, -mv_env%w_fac, T_small_rows, virt*nb, &
573 0 : mv_env%B_bar_ia(:, :, p), virt, 1.0_dp, RB_virt, virt*nb)
574 : END IF
575 : END IF
576 : END DO
577 : END IF
578 84 : CALL timestop(handle2)
579 :
580 : ! Pack by destination rank. Screening has distinct A/B layouts; without it,
581 : ! both outputs share the canonical exchange image. result_buf is free in either case.
582 84 : CALL timeset(routineN//"_reduce", handle2)
583 252 : DO peer = 1, mv_env%para_env%num_pe
584 168 : rows = mv_env%row_count(peer)
585 168 : first_row = mv_env%row_displ(peer) + 1
586 780 : DO k = 1, nb
587 528 : offset = n_outputs*displs(peer) + (k - 1)*rows
588 696 : IF (do_screened) THEN
589 13200 : DO iloc = 1, rows
590 12672 : ia = first_row + iloc - 1
591 12672 : i = occ_of_ia(ia, virt)
592 12672 : a = virt_of_ia(ia, virt)
593 12672 : result_buf(offset + iloc) = RA_screen(i, k, a)
594 13200 : IF (do_B) result_buf(offset + rows*nb + iloc) = RB_screen(a, k, i)
595 : END DO
596 : ELSE
597 0 : last_row = first_row + rows - 1
598 0 : result_buf(offset + 1:offset + rows) = RA(first_row:last_row, k)
599 0 : IF (do_B) result_buf(offset + rows*nb + 1:offset + rows*nb + rows) = RA(first_row:last_row, k)
600 : END IF
601 : END DO
602 : END DO
603 252 : counts = counts*n_outputs
604 : ! The wrapper's serial branch also expects a single-column flattened send buffer.
605 84 : packed(1:n_outputs*nvec, 1:1) => result_buf(1:n_outputs*nvec)
606 84 : CALL mv_env%para_env%sum_scatter(packed, work_buf(1:n_outputs*nrow_local*nb), counts)
607 84 : local_result(1:nrow_local, 1:n_outputs*nb) => work_buf(1:n_outputs*nrow_local*nb)
608 348 : DO k = 1, nb
609 6600 : DO iloc = 1, nrow_local
610 6336 : ia = row_indices(iloc)
611 : fm_AZ%local_data(iloc, col + k - 1) = local_result(iloc, k) + &
612 6600 : mv_env%eps_diff(ia)*fm_Z%local_data(iloc, col + k - 1)
613 : END DO
614 5292 : IF (do_B) fm_BZ%local_data(1:nrow_local, col + col_shift_B + k - 1) = local_result(:, nb + k)
615 : END DO
616 588 : CALL timestop(handle2)
617 : END DO
618 :
619 84 : DEALLOCATE (work_buf, result_buf, T_Pk)
620 :
621 84 : CALL timestop(handle)
622 :
623 168 : END SUBROUTINE bse_matvec_apply
624 :
625 : ! **************************************************************************************************
626 : !> \brief Diagonal used by the Davidson correction, either ε_a-ε_i or the full diagonal
627 : !> A_ia,ia = ε_a-ε_i + α sum_P (B^P_ia)^2 - w sum_P \bar{B}^P_ii B^P_aa
628 : !> with α, w and \bar{B} as in bse_matvec_apply
629 : !> \param mv_env slabs, prefactors and transition energies
630 : !> \param precond_kind bse_precond_full_diag or bse_precond_qp_diff
631 : !> \param diag replicated, size N_ov
632 : ! **************************************************************************************************
633 8 : SUBROUTINE bse_matvec_diagonal(mv_env, precond_kind, diag)
634 :
635 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
636 : INTEGER, INTENT(IN) :: precond_kind
637 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: diag
638 :
639 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_matvec_diagonal'
640 :
641 : INTEGER :: a, handle, i, ia, p
642 :
643 8 : CALL timeset(routineN, handle)
644 :
645 392 : diag(:) = 0.0_dp
646 8 : IF (precond_kind == bse_precond_full_diag) THEN
647 340 : DO p = 1, mv_env%n_ri_loc
648 1668 : DO i = 1, mv_env%homo
649 17596 : DO a = 1, mv_env%virt
650 15936 : ia = ia_of_occ_virt(i, a, mv_env%virt)
651 15936 : diag(ia) = diag(ia) + mv_env%alpha*mv_env%B_ia(a, i, p)**2
652 17264 : IF (mv_env%w_fac /= 0.0_dp) THEN
653 15936 : diag(ia) = diag(ia) - mv_env%w_fac*mv_env%B_bar_ij(i, i, p)*mv_env%B_ab(a, a, p)
654 : END IF
655 : END DO
656 : END DO
657 : END DO
658 : ! one n_ov vector, summed once per solve; replicated on every rank on purpose
659 776 : CALL mv_env%para_env%sum(diag)
660 : END IF
661 392 : diag(:) = diag(:) + mv_env%eps_diff(:)
662 :
663 8 : CALL timestop(handle)
664 :
665 8 : END SUBROUTINE bse_matvec_diagonal
666 :
667 : ! **************************************************************************************************
668 : !> \brief Exact A (and B) on a list of transitions, replicated on every rank,
669 : !> A_kl = δ_kl (ε_a-ε_i) + α sum_P B^P_ia B^P_jb - w sum_P \bar{B}^P_ij B^P_ab
670 : !> B_kl = α sum_P B^P_ia B^P_jb - w sum_P \bar{B}^P_ib B^P_ja
671 : !> with k = (i,a), l = (j,b) and α, w, \bar{B} as in bse_matvec_apply
672 : !> \param mv_env slabs, prefactors and transition energies
673 : !> \param ia_list global transition indices ia = (i-1)*virt + a of the block
674 : !> \param A_sub SIZE(ia_list) x SIZE(ia_list)
675 : !> \param B_sub same, only formed when present
676 : ! **************************************************************************************************
677 2 : SUBROUTINE bse_matvec_subblock(mv_env, ia_list, A_sub, B_sub)
678 :
679 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
680 : INTEGER, DIMENSION(:), INTENT(IN) :: ia_list
681 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: A_sub
682 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
683 : OPTIONAL :: B_sub
684 :
685 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_matvec_subblock'
686 :
687 : INTEGER :: handle, k, l, n_sub, p
688 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: a_virt, i_occ
689 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: B_Pk
690 :
691 2 : CALL timeset(routineN, handle)
692 :
693 2 : n_sub = SIZE(ia_list)
694 14 : ALLOCATE (i_occ(n_sub), a_virt(n_sub), B_Pk(MAX(mv_env%n_ri_loc, 1), n_sub))
695 98 : DO k = 1, n_sub
696 96 : i_occ(k) = occ_of_ia(ia_list(k), mv_env%virt)
697 98 : a_virt(k) = virt_of_ia(ia_list(k), mv_env%virt)
698 : END DO
699 :
700 : ! exchange, identical in A and B: α sum_P B_Pk B_Pl with B_Pk = B^P_ia at ia = ia_list(k)
701 4706 : A_sub(:, :) = 0.0_dp
702 2 : IF (mv_env%alpha /= 0.0_dp .AND. mv_env%n_ri_loc > 0) THEN
703 : ! P outermost: the gathered entries of one slice lie within one contiguous (a, i) block
704 85 : DO p = 1, mv_env%n_ri_loc
705 4069 : DO k = 1, n_sub
706 4067 : B_Pk(p, k) = mv_env%B_ia(a_virt(k), i_occ(k), p)
707 : END DO
708 : END DO
709 : CALL DGEMM('T', 'N', n_sub, n_sub, mv_env%n_ri_loc, mv_env%alpha, B_Pk, SIZE(B_Pk, 1), B_Pk, SIZE(B_Pk, 1), &
710 2 : 0.0_dp, A_sub, n_sub)
711 : END IF
712 2 : IF (PRESENT(B_sub)) B_sub(:, :) = A_sub(:, :)
713 :
714 2 : IF (mv_env%w_fac /= 0.0_dp) THEN
715 85 : DO p = 1, mv_env%n_ri_loc
716 4067 : DO l = 1, n_sub
717 195299 : DO k = 1, n_sub
718 195216 : A_sub(k, l) = A_sub(k, l) - mv_env%w_fac*mv_env%B_bar_ij(i_occ(l), i_occ(k), p)*mv_env%B_ab(a_virt(l), a_virt(k), p)
719 : END DO
720 : END DO
721 85 : IF (PRESENT(B_sub)) THEN
722 0 : DO l = 1, n_sub
723 0 : DO k = 1, n_sub
724 0 : B_sub(k, l) = B_sub(k, l) - mv_env%w_fac*mv_env%B_ia(a_virt(k), i_occ(l), p)*mv_env%B_bar_ia(a_virt(l), i_occ(k), p)
725 : END DO
726 : END DO
727 : END IF
728 : END DO
729 : END IF
730 :
731 9410 : CALL mv_env%para_env%sum(A_sub)
732 2 : IF (PRESENT(B_sub)) CALL mv_env%para_env%sum(B_sub)
733 98 : DO k = 1, n_sub
734 98 : A_sub(k, k) = A_sub(k, k) + mv_env%eps_diff(ia_list(k))
735 : END DO
736 :
737 2 : DEALLOCATE (i_occ, a_virt, B_Pk)
738 :
739 2 : CALL timestop(handle)
740 :
741 2 : END SUBROUTINE bse_matvec_subblock
742 :
743 : ! **************************************************************************************************
744 : !> \brief Debug check of the matrix-free application against the explicit matrices A (and B),
745 : !> dev = max_ia,k |(A Z)_ia,k - sum_jb A_ia,jb Z_jb,k| over up to eight unit vectors Z_ia,k = δ_ia,k
746 : !> and one dense vector Z_ia = sin(ia), the same for B, and max_ia |d_ia - A_ia,ia| for the diagonal
747 : !> \param mv_env the environment under test
748 : !> \param fm_A_explicit A from create_A_and_B, N_ov x N_ov on the grid of the caller
749 : !> \param unit_nr output unit, positive on the writing rank only
750 : !> \param fm_B_explicit B, present for an ABBA run
751 : ! **************************************************************************************************
752 0 : SUBROUTINE bse_matvec_selfcheck(mv_env, fm_A_explicit, unit_nr, fm_B_explicit)
753 :
754 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
755 : TYPE(cp_fm_type), INTENT(IN) :: fm_A_explicit
756 : INTEGER, INTENT(IN) :: unit_nr
757 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_B_explicit
758 :
759 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_matvec_selfcheck'
760 :
761 : INTEGER :: handle, ia, iloc, k, n_ov, n_unit, ncol, &
762 : nrow_local
763 0 : INTEGER, DIMENSION(:), POINTER :: row_indices
764 : REAL(KIND=dp) :: dev_A, dev_B, dev_diag
765 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag
766 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mv_result, ref_matrix, ref_result, Z
767 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
768 : TYPE(cp_fm_type) :: fm_AZ, fm_BZ, fm_Z
769 :
770 0 : CALL timeset(routineN, handle)
771 :
772 0 : n_ov = mv_env%n_ov
773 : ! eight unit vectors read eight columns of the explicit matrix exactly; the dense column covers the whole contraction
774 0 : n_unit = MIN(n_ov, 8)
775 0 : dev_B = 0.0_dp
776 0 : ncol = n_unit + 1
777 :
778 : ! trial vectors: e_k for k <= n_unit, then one dense column Z_ia = sin(ia)
779 0 : ALLOCATE (Z(n_ov, ncol))
780 0 : Z(:, :) = 0.0_dp
781 0 : DO k = 1, n_unit
782 0 : Z(k, k) = 1.0_dp
783 : END DO
784 0 : DO ia = 1, n_ov
785 0 : Z(ia, ncol) = SIN(REAL(ia, dp))
786 : END DO
787 :
788 0 : CALL bse_matvec_vector_struct(mv_env, ncol, fm_struct)
789 0 : CALL cp_fm_create(fm_Z, fm_struct, name="fm_Z_selfcheck")
790 0 : CALL cp_fm_create(fm_AZ, fm_struct, name="fm_AZ_selfcheck")
791 0 : CALL cp_fm_set_all(fm_AZ, 0.0_dp)
792 0 : IF (PRESENT(fm_B_explicit)) THEN
793 0 : CALL cp_fm_create(fm_BZ, fm_struct, name="fm_BZ_selfcheck")
794 0 : CALL cp_fm_set_all(fm_BZ, 0.0_dp)
795 : END IF
796 0 : CALL cp_fm_struct_release(fm_struct)
797 :
798 0 : CALL cp_fm_get_info(fm_Z, nrow_local=nrow_local, row_indices=row_indices)
799 0 : DO k = 1, ncol
800 0 : DO iloc = 1, nrow_local
801 0 : fm_Z%local_data(iloc, k) = Z(row_indices(iloc), k)
802 : END DO
803 : END DO
804 :
805 0 : IF (PRESENT(fm_B_explicit)) THEN
806 0 : CALL bse_matvec_apply(mv_env, fm_Z, 1, ncol, fm_AZ, fm_BZ)
807 : ELSE
808 0 : CALL bse_matvec_apply(mv_env, fm_Z, 1, ncol, fm_AZ)
809 : END IF
810 :
811 0 : ALLOCATE (ref_matrix(n_ov, n_ov), mv_result(n_ov, ncol), ref_result(n_ov, ncol), diag(n_ov))
812 :
813 0 : CALL cp_fm_get_submatrix(fm_A_explicit, ref_matrix)
814 0 : CALL cp_fm_get_submatrix(fm_AZ, mv_result)
815 0 : ref_result(:, :) = MATMUL(ref_matrix, Z)
816 0 : dev_A = MAXVAL(ABS(mv_result - ref_result))
817 0 : CALL bse_matvec_diagonal(mv_env, bse_precond_full_diag, diag)
818 0 : dev_diag = 0.0_dp
819 0 : DO ia = 1, n_ov
820 0 : dev_diag = MAX(dev_diag, ABS(diag(ia) - ref_matrix(ia, ia)))
821 : END DO
822 :
823 0 : IF (PRESENT(fm_B_explicit)) THEN
824 0 : CALL cp_fm_get_submatrix(fm_B_explicit, ref_matrix)
825 0 : CALL cp_fm_get_submatrix(fm_BZ, mv_result)
826 0 : ref_result(:, :) = MATMUL(ref_matrix, Z)
827 0 : dev_B = MAXVAL(ABS(mv_result - ref_result))
828 : END IF
829 :
830 0 : IF (unit_nr > 0) THEN
831 0 : WRITE (unit_nr, '(T2,A10,T13,A,T59,ES22.6)') 'BSE|DEBUG|', &
832 0 : 'Max deviation matvec vs explicit A', dev_A
833 0 : IF (PRESENT(fm_B_explicit)) THEN
834 0 : WRITE (unit_nr, '(T2,A10,T13,A,T59,ES22.6)') 'BSE|DEBUG|', &
835 0 : 'Max deviation matvec vs explicit B', dev_B
836 : END IF
837 0 : WRITE (unit_nr, '(T2,A10,T13,A,T59,ES22.6)') 'BSE|DEBUG|', &
838 0 : 'Max deviation diagonal vs explicit A', dev_diag
839 : END IF
840 :
841 0 : DEALLOCATE (Z, ref_matrix, mv_result, ref_result, diag)
842 0 : CALL cp_fm_release(fm_Z)
843 0 : CALL cp_fm_release(fm_AZ)
844 0 : IF (PRESENT(fm_B_explicit)) CALL cp_fm_release(fm_BZ)
845 :
846 0 : CALL timestop(handle)
847 :
848 0 : END SUBROUTINE bse_matvec_selfcheck
849 :
850 0 : END MODULE bse_matvec
|
