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 :: col, col_shift_B, handle, handle2, ia, &
391 : iloc, k, me, n_work_bufs, nb, nb_max, &
392 : nrow_local, p
393 84 : INTEGER, DIMENSION(:), POINTER :: row_indices
394 : LOGICAL :: do_B
395 : REAL(KIND=dp) :: mem_avail_GB
396 84 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: R_loc
397 84 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: BabZ_buf, RA_buf, RB_buf, Z_buf
398 84 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: BjaZ_ab, T_Pk
399 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
400 84 : POINTER :: B_ia_2d, BabZ_wide, RA, RB, Z, Z_wide
401 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
402 84 : POINTER :: BabZ_3d, RA_3d, RB_3d, Z_3d
403 :
404 84 : CALL timeset(routineN, handle)
405 :
406 84 : do_B = PRESENT(fm_BZ)
407 84 : IF (do_B) THEN
408 : ! B Z needs the \bar{B}^P_ia slab, which is sliced for do_abba only
409 66 : CPASSERT(mv_env%do_abba)
410 : END IF
411 84 : col_shift_B = 0
412 84 : IF (PRESENT(first_col_BZ)) col_shift_B = first_col_BZ - first_col
413 :
414 84 : CALL cp_fm_get_info(fm_Z, nrow_local=nrow_local, row_indices=row_indices)
415 : ! the collectives below read the rows of one rank as one range, as bse_matvec_vector_struct
416 : ! lays them out; a vector built elsewhere would scatter into the wrong rows
417 84 : me = mv_env%para_env%mepos + 1
418 84 : CPASSERT(nrow_local == mv_env%row_count(me))
419 84 : CPASSERT(nrow_local == 0 .OR. row_indices(1) == mv_env%row_displ(me) + 1)
420 252 : ALLOCATE (R_loc(nrow_local))
421 :
422 : ! columns per pass from the memory of the replicated work arrays Z, B_ab Z, R^A (and R^B)
423 84 : IF (do_B) THEN
424 : n_work_bufs = 4
425 : ELSE
426 18 : n_work_bufs = 3
427 : END IF
428 84 : IF (mv_env%block_cols > 0) THEN
429 : ! pinned by input: the replicated work arrays are then a known n_work_bufs*n_ov*nb doubles
430 84 : nb_max = MIN(ncol, mv_env%block_cols)
431 : ELSE
432 : ! the pass takes what MEMORY_PER_PROC leaves after the solver's arrays
433 0 : mem_avail_GB = bse_memory_left_GB(mv_env%memory_per_proc_gb, mv_env%para_env)
434 : nb_max = INT(MIN(REAL(ncol, dp), &
435 0 : mem_avail_GB*1.0E9_dp/(8.0_dp*REAL(n_work_bufs, dp)*REAL(mv_env%n_ov, dp))))
436 : END IF
437 84 : nb_max = MAX(nb_max, 1)
438 :
439 420 : ALLOCATE (Z_buf(mv_env%n_ov*nb_max), RA_buf(mv_env%n_ov*nb_max), BabZ_buf(mv_env%n_ov*nb_max))
440 84 : IF (do_B) THEN
441 330 : ALLOCATE (RB_buf(mv_env%n_ov*nb_max), BjaZ_ab(mv_env%virt, mv_env%virt))
442 : END IF
443 336 : ALLOCATE (T_Pk(mv_env%n_ri_loc, nb_max))
444 84 : B_ia_2d(1:mv_env%n_ov, 1:mv_env%n_ri_loc) => mv_env%B_ia
445 :
446 168 : DO col = first_col, first_col + ncol - 1, nb_max
447 84 : nb = MIN(nb_max, first_col + ncol - col)
448 : ! three views of one contiguous buffer, bounds-remapping pointer assignment (Fortran 2003):
449 : ! Z_ia,k for the exchange DGEMM, Z_3d(a, i, k) per column, Z_wide(a, (i k)) for the wide product
450 84 : Z(1:mv_env%n_ov, 1:nb) => Z_buf(1:mv_env%n_ov*nb)
451 84 : Z_3d(1:mv_env%virt, 1:mv_env%homo, 1:nb) => Z_buf(1:mv_env%n_ov*nb)
452 84 : Z_wide(1:mv_env%virt, 1:mv_env%homo*nb) => Z_buf(1:mv_env%n_ov*nb)
453 84 : RA(1:mv_env%n_ov, 1:nb) => RA_buf(1:mv_env%n_ov*nb)
454 84 : RA_3d(1:mv_env%virt, 1:mv_env%homo, 1:nb) => RA_buf(1:mv_env%n_ov*nb)
455 84 : BabZ_3d(1:mv_env%virt, 1:mv_env%homo, 1:nb) => BabZ_buf(1:mv_env%n_ov*nb)
456 84 : BabZ_wide(1:mv_env%virt, 1:mv_env%homo*nb) => BabZ_buf(1:mv_env%n_ov*nb)
457 84 : IF (do_B) THEN
458 66 : RB(1:mv_env%n_ov, 1:nb) => RB_buf(1:mv_env%n_ov*nb)
459 66 : RB_3d(1:mv_env%virt, 1:mv_env%homo, 1:nb) => RB_buf(1:mv_env%n_ov*nb)
460 : END IF
461 :
462 : ! Z_ia,k on every rank, by design: each rank contracts its whole RI slices against all of Z (header of
463 : ! bse_matvec_create), so the trial vectors are the one replicated buffer, n_ov x nb per matvec at any rank
464 : ! count; every rank owns a contiguous ia range, so the column arrives sorted
465 84 : CALL timeset(routineN//"_gather", handle2)
466 348 : DO k = 1, nb
467 : CALL mv_env%para_env%allgatherv(fm_Z%local_data(1:nrow_local, col + k - 1), Z(:, k), &
468 348 : mv_env%row_count, mv_env%row_displ)
469 : END DO
470 84 : CALL timestop(handle2)
471 :
472 : ! 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
473 84 : CALL timeset(routineN//"_exchange", handle2)
474 84 : RA(:, :) = 0.0_dp
475 84 : IF (mv_env%alpha /= 0.0_dp .AND. mv_env%n_ri_loc > 0) THEN
476 : 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, &
477 84 : 0.0_dp, T_Pk, mv_env%n_ri_loc)
478 : 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, &
479 84 : 0.0_dp, RA, mv_env%n_ov)
480 : END IF
481 20338 : IF (do_B) RB(:, :) = RA(:, :)
482 84 : CALL timestop(handle2)
483 :
484 84 : CALL timeset(routineN//"_W", handle2)
485 84 : IF (mv_env%w_fac /= 0.0_dp) THEN
486 3570 : DO p = 1, mv_env%n_ri_loc
487 : ! A: R_ia,k -= w sum_jb \bar{B}^P_ij B^P_ab Z_jb,k, BabZ_ib,k = sum_b B^P_ab Z_jb,k first
488 : CALL DGEMM('T', 'N', mv_env%virt, mv_env%homo*nb, mv_env%virt, 1.0_dp, mv_env%B_ab(:, :, p), mv_env%virt, &
489 3486 : Z_wide, mv_env%virt, 0.0_dp, BabZ_wide, mv_env%virt)
490 14442 : DO k = 1, nb
491 : CALL DGEMM('N', 'N', mv_env%virt, mv_env%homo, mv_env%homo, -mv_env%w_fac, BabZ_3d(:, :, k), mv_env%virt, &
492 14442 : mv_env%B_bar_ij(:, :, p), mv_env%homo, 1.0_dp, RA_3d(:, :, k), mv_env%virt)
493 : END DO
494 : ! B: R_ia,k -= w sum_jb \bar{B}^P_ib B^P_ja Z_jb,k, BjaZ_ab = sum_j B^P_ja Z_jb,k first
495 3570 : IF (do_B) THEN
496 11288 : DO k = 1, nb
497 : CALL DGEMM('N', 'T', mv_env%virt, mv_env%virt, mv_env%homo, 1.0_dp, mv_env%B_ia(:, :, p), mv_env%virt, &
498 8549 : Z_3d(:, :, k), mv_env%virt, 0.0_dp, BjaZ_ab, mv_env%virt)
499 : CALL DGEMM('N', 'N', mv_env%virt, mv_env%homo, mv_env%virt, -mv_env%w_fac, BjaZ_ab, mv_env%virt, &
500 11288 : mv_env%B_bar_ia(:, :, p), mv_env%virt, 1.0_dp, RB_3d(:, :, k), mv_env%virt)
501 : END DO
502 : END IF
503 : END DO
504 : END IF
505 84 : CALL timestop(handle2)
506 :
507 : ! the sum over the ranks lands on the owner of each row, never replicated
508 84 : CALL timeset(routineN//"_reduce", handle2)
509 : ! (A Z)_ia,k = (ε_a-ε_i) Z_ia,k + R^A_ia,k and (B Z)_ia,k = R^B_ia,k on the local rows
510 348 : DO k = 1, nb
511 264 : CALL mv_env%para_env%sum_scatter(RA(:, k:k), R_loc, mv_env%row_count)
512 6600 : DO iloc = 1, nrow_local
513 6336 : ia = row_indices(iloc)
514 6600 : fm_AZ%local_data(iloc, col + k - 1) = R_loc(iloc) + mv_env%eps_diff(ia)*Z(ia, k)
515 : END DO
516 348 : IF (do_B) THEN
517 206 : CALL mv_env%para_env%sum_scatter(RB(:, k:k), R_loc, mv_env%row_count)
518 5150 : fm_BZ%local_data(1:nrow_local, col + col_shift_B + k - 1) = R_loc(1:nrow_local)
519 : END IF
520 : END DO
521 504 : CALL timestop(handle2)
522 : END DO
523 :
524 84 : DEALLOCATE (Z_buf, RA_buf, BabZ_buf, T_Pk)
525 84 : IF (do_B) THEN
526 66 : DEALLOCATE (RB_buf, BjaZ_ab)
527 : END IF
528 :
529 84 : CALL timestop(handle)
530 :
531 252 : END SUBROUTINE bse_matvec_apply
532 :
533 : ! **************************************************************************************************
534 : !> \brief Diagonal used by the Davidson correction, either ε_a-ε_i or the full diagonal
535 : !> A_ia,ia = ε_a-ε_i + α sum_P (B^P_ia)^2 - w sum_P \bar{B}^P_ii B^P_aa
536 : !> with α, w and \bar{B} as in bse_matvec_apply
537 : !> \param mv_env slabs, prefactors and transition energies
538 : !> \param precond_kind bse_precond_full_diag or bse_precond_qp_diff
539 : !> \param diag replicated, size N_ov
540 : ! **************************************************************************************************
541 8 : SUBROUTINE bse_matvec_diagonal(mv_env, precond_kind, diag)
542 :
543 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
544 : INTEGER, INTENT(IN) :: precond_kind
545 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: diag
546 :
547 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_matvec_diagonal'
548 :
549 : INTEGER :: a, handle, i, ia, p
550 :
551 8 : CALL timeset(routineN, handle)
552 :
553 392 : diag(:) = 0.0_dp
554 8 : IF (precond_kind == bse_precond_full_diag) THEN
555 340 : DO p = 1, mv_env%n_ri_loc
556 1668 : DO i = 1, mv_env%homo
557 17596 : DO a = 1, mv_env%virt
558 15936 : ia = ia_of_occ_virt(i, a, mv_env%virt)
559 15936 : diag(ia) = diag(ia) + mv_env%alpha*mv_env%B_ia(a, i, p)**2
560 17264 : IF (mv_env%w_fac /= 0.0_dp) THEN
561 15936 : diag(ia) = diag(ia) - mv_env%w_fac*mv_env%B_bar_ij(i, i, p)*mv_env%B_ab(a, a, p)
562 : END IF
563 : END DO
564 : END DO
565 : END DO
566 : ! one n_ov vector, summed once per solve; replicated on every rank on purpose
567 776 : CALL mv_env%para_env%sum(diag)
568 : END IF
569 392 : diag(:) = diag(:) + mv_env%eps_diff(:)
570 :
571 8 : CALL timestop(handle)
572 :
573 8 : END SUBROUTINE bse_matvec_diagonal
574 :
575 : ! **************************************************************************************************
576 : !> \brief Exact A (and B) on a list of transitions, replicated on every rank,
577 : !> A_kl = δ_kl (ε_a-ε_i) + α sum_P B^P_ia B^P_jb - w sum_P \bar{B}^P_ij B^P_ab
578 : !> B_kl = α sum_P B^P_ia B^P_jb - w sum_P \bar{B}^P_ib B^P_ja
579 : !> with k = (i,a), l = (j,b) and α, w, \bar{B} as in bse_matvec_apply
580 : !> \param mv_env slabs, prefactors and transition energies
581 : !> \param ia_list global transition indices ia = (i-1)*virt + a of the block
582 : !> \param A_sub SIZE(ia_list) x SIZE(ia_list)
583 : !> \param B_sub same, only formed when present
584 : ! **************************************************************************************************
585 2 : SUBROUTINE bse_matvec_subblock(mv_env, ia_list, A_sub, B_sub)
586 :
587 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
588 : INTEGER, DIMENSION(:), INTENT(IN) :: ia_list
589 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: A_sub
590 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
591 : OPTIONAL :: B_sub
592 :
593 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_matvec_subblock'
594 :
595 : INTEGER :: handle, k, l, n_sub, p
596 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: a_virt, i_occ
597 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: B_Pk
598 :
599 2 : CALL timeset(routineN, handle)
600 :
601 2 : n_sub = SIZE(ia_list)
602 14 : ALLOCATE (i_occ(n_sub), a_virt(n_sub), B_Pk(MAX(mv_env%n_ri_loc, 1), n_sub))
603 98 : DO k = 1, n_sub
604 96 : i_occ(k) = occ_of_ia(ia_list(k), mv_env%virt)
605 98 : a_virt(k) = virt_of_ia(ia_list(k), mv_env%virt)
606 : END DO
607 :
608 : ! exchange, identical in A and B: α sum_P B_Pk B_Pl with B_Pk = B^P_ia at ia = ia_list(k)
609 4706 : A_sub(:, :) = 0.0_dp
610 2 : IF (mv_env%alpha /= 0.0_dp .AND. mv_env%n_ri_loc > 0) THEN
611 : ! P outermost: the gathered entries of one slice lie within one contiguous (a, i) block
612 85 : DO p = 1, mv_env%n_ri_loc
613 4069 : DO k = 1, n_sub
614 4067 : B_Pk(p, k) = mv_env%B_ia(a_virt(k), i_occ(k), p)
615 : END DO
616 : END DO
617 : 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), &
618 2 : 0.0_dp, A_sub, n_sub)
619 : END IF
620 2 : IF (PRESENT(B_sub)) B_sub(:, :) = A_sub(:, :)
621 :
622 2 : IF (mv_env%w_fac /= 0.0_dp) THEN
623 85 : DO p = 1, mv_env%n_ri_loc
624 4067 : DO l = 1, n_sub
625 195299 : DO k = 1, n_sub
626 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)
627 : END DO
628 : END DO
629 85 : IF (PRESENT(B_sub)) THEN
630 0 : DO l = 1, n_sub
631 0 : DO k = 1, n_sub
632 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)
633 : END DO
634 : END DO
635 : END IF
636 : END DO
637 : END IF
638 :
639 9410 : CALL mv_env%para_env%sum(A_sub)
640 2 : IF (PRESENT(B_sub)) CALL mv_env%para_env%sum(B_sub)
641 98 : DO k = 1, n_sub
642 98 : A_sub(k, k) = A_sub(k, k) + mv_env%eps_diff(ia_list(k))
643 : END DO
644 :
645 2 : DEALLOCATE (i_occ, a_virt, B_Pk)
646 :
647 2 : CALL timestop(handle)
648 :
649 2 : END SUBROUTINE bse_matvec_subblock
650 :
651 : ! **************************************************************************************************
652 : !> \brief Debug check of the matrix-free application against the explicit matrices A (and B),
653 : !> 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
654 : !> and one dense vector Z_ia = sin(ia), the same for B, and max_ia |d_ia - A_ia,ia| for the diagonal
655 : !> \param mv_env the environment under test
656 : !> \param fm_A_explicit A from create_A_and_B, N_ov x N_ov on the grid of the caller
657 : !> \param unit_nr output unit, positive on the writing rank only
658 : !> \param fm_B_explicit B, present for an ABBA run
659 : ! **************************************************************************************************
660 0 : SUBROUTINE bse_matvec_selfcheck(mv_env, fm_A_explicit, unit_nr, fm_B_explicit)
661 :
662 : TYPE(bse_matvec_env_type), INTENT(IN) :: mv_env
663 : TYPE(cp_fm_type), INTENT(IN) :: fm_A_explicit
664 : INTEGER, INTENT(IN) :: unit_nr
665 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_B_explicit
666 :
667 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_matvec_selfcheck'
668 :
669 : INTEGER :: handle, ia, iloc, k, n_ov, n_unit, ncol, &
670 : nrow_local
671 0 : INTEGER, DIMENSION(:), POINTER :: row_indices
672 : REAL(KIND=dp) :: dev_A, dev_B, dev_diag
673 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag
674 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mv_result, ref_matrix, ref_result, Z
675 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
676 : TYPE(cp_fm_type) :: fm_AZ, fm_BZ, fm_Z
677 :
678 0 : CALL timeset(routineN, handle)
679 :
680 0 : n_ov = mv_env%n_ov
681 : ! eight unit vectors read eight columns of the explicit matrix exactly; the dense column covers the whole contraction
682 0 : n_unit = MIN(n_ov, 8)
683 0 : dev_B = 0.0_dp
684 0 : ncol = n_unit + 1
685 :
686 : ! trial vectors: e_k for k <= n_unit, then one dense column Z_ia = sin(ia)
687 0 : ALLOCATE (Z(n_ov, ncol))
688 0 : Z(:, :) = 0.0_dp
689 0 : DO k = 1, n_unit
690 0 : Z(k, k) = 1.0_dp
691 : END DO
692 0 : DO ia = 1, n_ov
693 0 : Z(ia, ncol) = SIN(REAL(ia, dp))
694 : END DO
695 :
696 0 : CALL bse_matvec_vector_struct(mv_env, ncol, fm_struct)
697 0 : CALL cp_fm_create(fm_Z, fm_struct, name="fm_Z_selfcheck")
698 0 : CALL cp_fm_create(fm_AZ, fm_struct, name="fm_AZ_selfcheck")
699 0 : CALL cp_fm_set_all(fm_AZ, 0.0_dp)
700 0 : IF (PRESENT(fm_B_explicit)) THEN
701 0 : CALL cp_fm_create(fm_BZ, fm_struct, name="fm_BZ_selfcheck")
702 0 : CALL cp_fm_set_all(fm_BZ, 0.0_dp)
703 : END IF
704 0 : CALL cp_fm_struct_release(fm_struct)
705 :
706 0 : CALL cp_fm_get_info(fm_Z, nrow_local=nrow_local, row_indices=row_indices)
707 0 : DO k = 1, ncol
708 0 : DO iloc = 1, nrow_local
709 0 : fm_Z%local_data(iloc, k) = Z(row_indices(iloc), k)
710 : END DO
711 : END DO
712 :
713 0 : IF (PRESENT(fm_B_explicit)) THEN
714 0 : CALL bse_matvec_apply(mv_env, fm_Z, 1, ncol, fm_AZ, fm_BZ)
715 : ELSE
716 0 : CALL bse_matvec_apply(mv_env, fm_Z, 1, ncol, fm_AZ)
717 : END IF
718 :
719 0 : ALLOCATE (ref_matrix(n_ov, n_ov), mv_result(n_ov, ncol), ref_result(n_ov, ncol), diag(n_ov))
720 :
721 0 : CALL cp_fm_get_submatrix(fm_A_explicit, ref_matrix)
722 0 : CALL cp_fm_get_submatrix(fm_AZ, mv_result)
723 0 : ref_result(:, :) = MATMUL(ref_matrix, Z)
724 0 : dev_A = MAXVAL(ABS(mv_result - ref_result))
725 0 : CALL bse_matvec_diagonal(mv_env, bse_precond_full_diag, diag)
726 0 : dev_diag = 0.0_dp
727 0 : DO ia = 1, n_ov
728 0 : dev_diag = MAX(dev_diag, ABS(diag(ia) - ref_matrix(ia, ia)))
729 : END DO
730 :
731 0 : IF (PRESENT(fm_B_explicit)) THEN
732 0 : CALL cp_fm_get_submatrix(fm_B_explicit, ref_matrix)
733 0 : CALL cp_fm_get_submatrix(fm_BZ, mv_result)
734 0 : ref_result(:, :) = MATMUL(ref_matrix, Z)
735 0 : dev_B = MAXVAL(ABS(mv_result - ref_result))
736 : END IF
737 :
738 0 : IF (unit_nr > 0) THEN
739 0 : WRITE (unit_nr, '(T2,A10,T13,A,T59,ES22.6)') 'BSE|DEBUG|', &
740 0 : 'Max deviation matvec vs explicit A', dev_A
741 0 : IF (PRESENT(fm_B_explicit)) THEN
742 0 : WRITE (unit_nr, '(T2,A10,T13,A,T59,ES22.6)') 'BSE|DEBUG|', &
743 0 : 'Max deviation matvec vs explicit B', dev_B
744 : END IF
745 0 : WRITE (unit_nr, '(T2,A10,T13,A,T59,ES22.6)') 'BSE|DEBUG|', &
746 0 : 'Max deviation diagonal vs explicit A', dev_diag
747 : END IF
748 :
749 0 : DEALLOCATE (Z, ref_matrix, mv_result, ref_result, diag)
750 0 : CALL cp_fm_release(fm_Z)
751 0 : CALL cp_fm_release(fm_AZ)
752 0 : IF (PRESENT(fm_B_explicit)) CALL cp_fm_release(fm_BZ)
753 :
754 0 : CALL timestop(handle)
755 :
756 0 : END SUBROUTINE bse_matvec_selfcheck
757 :
758 0 : END MODULE bse_matvec
|
