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 DBCSR collocation, layouts and panel contractions for real-space resolution of the identity.
10 : ! **************************************************************************************************
11 : MODULE ri_rs_matrices
12 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
13 : gto_basis_set_type
14 : USE cp_dbcsr_api, ONLY: &
15 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_get, dbcsr_distribution_new, &
16 : dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_finalize, dbcsr_get_data_size, &
17 : dbcsr_get_info, dbcsr_get_occupation, dbcsr_iterator_blocks_left, &
18 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
19 : dbcsr_multiply, dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_set, dbcsr_type, &
20 : dbcsr_type_no_symmetry
21 : USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,&
22 : max_elements_per_block
23 : USE cp_fm_types, ONLY: cp_fm_create,&
24 : cp_fm_release,&
25 : cp_fm_to_fm,&
26 : cp_fm_type
27 : USE gw_auto_ri_types, ONLY: auto_ri_type
28 : USE gw_utils_dbcsr, ONLY: hadamard_product_inplace
29 : USE gw_utils_fm, ONLY: fm_contract_ABA
30 : USE kinds, ONLY: dp,&
31 : int_8,&
32 : max_line_length
33 : USE libint_2c_3c, ONLY: libint_potential_type
34 : USE machine, ONLY: m_flush
35 : USE message_passing, ONLY: mp_mem_avail_per_rank_GB,&
36 : mp_para_env_type
37 : USE mp2_ri_2c, ONLY: RI_2c_integral_mat
38 : USE qs_environment_types, ONLY: qs_environment_type
39 : USE ri_rs_types, ONLY: ri_rs_env
40 : USE ri_rs_utils, ONLY: evaluate_ao_basis_on_points
41 : #include "./base/base_uses.f90"
42 :
43 : IMPLICIT NONE
44 : PRIVATE
45 : PUBLIC :: atomic_basis_at_grid_point, setup_square_topology, &
46 : release_dbcsr_topology_and_matrices, reserve_blocks_within_radius, &
47 : print_matrix_occupation, &
48 : contract_grid_panels, mask_grid_blocks_near_panel, &
49 : panel_template_elems, panel_mem_estimate_GB, resolve_grid_panels, &
50 : build_ri_rs_bare_aux_interaction
51 :
52 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ri_rs_matrices'
53 :
54 : CONTAINS
55 :
56 : ! **************************************************************************************************
57 : !> \brief Builds the grid-by-AO collocation matrix with bounded dense evaluation panels.
58 : !> The molecular path uses atom-aligned row blocks and minimum-image AO values.
59 : !> The Gamma path uses uniform row blocks and sums neighbouring periodic images.
60 : !> \param ri_rs grid, borrowed geometry references and output collocation matrix
61 : !> \param basis_set_AO AO basis sets, one per atomic kind
62 : !> \param matrix_template template for the process grid and AO column block sizes
63 : !> \param para_env parallel environment
64 : !> \param eps_filter block screening threshold
65 : !> \param msg_elem_limit maximum number of elements per DBCSR message
66 : !> \param unit_nr output unit
67 : !> \param periodic_images select neighbouring-image summation for Gamma
68 : ! **************************************************************************************************
69 48 : SUBROUTINE atomic_basis_at_grid_point(ri_rs, basis_set_AO, matrix_template, &
70 : para_env, eps_filter, msg_elem_limit, unit_nr, periodic_images)
71 : TYPE(ri_rs_env), INTENT(INOUT) :: ri_rs
72 : TYPE(gto_basis_set_p_type), INTENT(IN) :: basis_set_AO(:)
73 : TYPE(dbcsr_type), INTENT(IN) :: matrix_template
74 : TYPE(mp_para_env_type), POINTER :: para_env
75 : REAL(KIND=dp), INTENT(IN) :: eps_filter
76 : INTEGER(KIND=int_8), INTENT(IN) :: msg_elem_limit
77 : INTEGER, INTENT(IN) :: unit_nr
78 : LOGICAL, INTENT(IN) :: periodic_images
79 :
80 : CHARACTER(LEN=*), PARAMETER :: routineN = 'atomic_basis_at_grid_point'
81 : INTEGER, PARAMETER :: grid_panel_size = 4096
82 :
83 : INTEGER :: blocks_per_panel, bs_eff, c_size, first_block, handle, i, i_blk, ia, iatom, &
84 : ikind, last_block, n_grid, natom, npcol, nprow, num_grid_chunks, panel_size, panel_start, &
85 : r_end, r_start, remaining, row_first, row_last, run, safe_max
86 48 : INTEGER, ALLOCATABLE :: blk_row_start(:)
87 48 : INTEGER, POINTER :: c_blk_sizes(:), col_dist(:), &
88 48 : r_blk_sizes(:), row_dist(:)
89 : REAL(KIND=dp) :: cutoff_squared
90 48 : REAL(KIND=dp), ALLOCATABLE :: atom_col_buffer(:, :)
91 : TYPE(dbcsr_distribution_type) :: dist, template_dist
92 : TYPE(gto_basis_set_type), POINTER :: ao_basis
93 :
94 48 : CALL timeset(routineN, handle)
95 48 : CPASSERT(ASSOCIATED(ri_rs%cell))
96 48 : CPASSERT(ASSOCIATED(para_env))
97 48 : CPASSERT(ASSOCIATED(ri_rs%particle_set))
98 48 : CALL dbcsr_get_info(matrix_template, distribution=template_dist, col_blk_size=c_blk_sizes)
99 48 : natom = SIZE(c_blk_sizes)
100 48 : n_grid = SIZE(ri_rs%grid_points, 2)
101 48 : CPASSERT(SIZE(ri_rs%particle_set) == natom)
102 48 : CPASSERT(n_grid == ri_rs%n_grid_points)
103 48 : CALL dbcsr_distribution_get(template_dist, nprows=nprow, npcols=npcol)
104 :
105 48 : bs_eff = max_elements_per_block
106 48 : IF (.NOT. periodic_images) THEN
107 48 : safe_max = INT(msg_elem_limit*MAX(MIN(nprow, npcol), 1)/n_grid/2)
108 48 : bs_eff = MAX(1, MIN(bs_eff, MAX(1, safe_max)))
109 : END IF
110 48 : IF (periodic_images) THEN
111 0 : num_grid_chunks = n_grid/bs_eff
112 0 : IF (MOD(n_grid, bs_eff) /= 0) num_grid_chunks = num_grid_chunks + 1
113 : ELSE
114 48 : num_grid_chunks = 0
115 166 : DO ia = 1, natom
116 118 : run = ri_rs%grid_atom_boundaries(ia + 1) - ri_rs%grid_atom_boundaries(ia)
117 166 : IF (run > 0) num_grid_chunks = num_grid_chunks + (run + bs_eff - 1)/bs_eff
118 : END DO
119 : END IF
120 192 : ALLOCATE (r_blk_sizes(num_grid_chunks), blk_row_start(num_grid_chunks))
121 48 : i_blk = 0
122 48 : r_start = 1
123 214 : DO ia = 1, MERGE(1, natom, periodic_images)
124 118 : IF (periodic_images) THEN
125 : remaining = n_grid
126 : ELSE
127 118 : remaining = ri_rs%grid_atom_boundaries(ia + 1) - ri_rs%grid_atom_boundaries(ia)
128 : END IF
129 780 : DO WHILE (remaining > 0)
130 614 : i_blk = i_blk + 1
131 614 : r_blk_sizes(i_blk) = MIN(bs_eff, remaining)
132 614 : blk_row_start(i_blk) = r_start
133 614 : r_start = r_start + r_blk_sizes(i_blk)
134 614 : remaining = remaining - r_blk_sizes(i_blk)
135 : END DO
136 : END DO
137 :
138 48 : IF (.NOT. periodic_images) THEN
139 48 : IF (unit_nr > 0) THEN
140 24 : WRITE (unit_nr, '(T2,A,T71,I12)') 'RI-RS grid row-blocks of ϕ_μ(r_l)', num_grid_chunks
141 24 : WRITE (unit_nr, '(T2,A,T69,I12)') 'RI-RS grid points per block (max)', bs_eff
142 : END IF
143 48 : IF (ri_rs%cutoff_radius_v_w > 0.0_dp .OR. ri_rs%cutoff_radius_w0 > 0.0_dp) THEN
144 8 : IF (ALLOCATED(ri_rs%chunk_centroids)) DEALLOCATE (ri_rs%chunk_centroids)
145 24 : ALLOCATE (ri_rs%chunk_centroids(3, num_grid_chunks))
146 110 : DO i_blk = 1, num_grid_chunks
147 102 : r_start = blk_row_start(i_blk)
148 102 : r_end = r_start + r_blk_sizes(i_blk) - 1
149 456 : DO i = 1, 3
150 : ri_rs%chunk_centroids(i, i_blk) = &
151 8952 : SUM(ri_rs%grid_points(i, r_start:r_end))/REAL(r_blk_sizes(i_blk), dp)
152 : END DO
153 : END DO
154 : END IF
155 : END IF
156 :
157 192 : ALLOCATE (row_dist(num_grid_chunks), col_dist(natom))
158 662 : DO i = 1, num_grid_chunks
159 662 : row_dist(i) = MOD(i - 1, nprow)
160 : END DO
161 166 : DO i = 1, natom
162 166 : col_dist(i) = MOD(i - 1, npcol)
163 : END DO
164 48 : CALL dbcsr_distribution_new(dist, template=template_dist, row_dist=row_dist, col_dist=col_dist)
165 : CALL dbcsr_create(ri_rs%mat_phi_mu_l, name='phi_val_sparse', dist=dist, &
166 48 : matrix_type=dbcsr_type_no_symmetry, row_blk_size=r_blk_sizes, col_blk_size=c_blk_sizes)
167 :
168 662 : blocks_per_panel = MAX(1, grid_panel_size/MAXVAL(r_blk_sizes))
169 48 : DO iatom = para_env%mepos + 1, natom, para_env%num_pe
170 59 : c_size = c_blk_sizes(iatom)
171 999 : ALLOCATE (atom_col_buffer(MIN(n_grid, blocks_per_panel*MAXVAL(r_blk_sizes)), c_size))
172 59 : cutoff_squared = ri_rs%radius_ao_per_atom(iatom)**2
173 59 : IF (.NOT. periodic_images .AND. ri_rs%cutoff_radius_ri_ao > 0.0_dp) THEN
174 3 : cutoff_squared = ri_rs%cutoff_radius_ri_ao**2
175 : END IF
176 59 : ikind = ri_rs%particle_set(iatom)%atomic_kind%kind_number
177 59 : ao_basis => basis_set_AO(ikind)%gto_basis_set
178 118 : DO first_block = 1, num_grid_chunks, blocks_per_panel
179 59 : last_block = MIN(first_block + blocks_per_panel - 1, num_grid_chunks)
180 59 : panel_start = blk_row_start(first_block)
181 59 : panel_size = blk_row_start(last_block) + r_blk_sizes(last_block) - panel_start
182 79323 : atom_col_buffer(:panel_size, :) = 0.0_dp
183 : CALL evaluate_ao_basis_on_points(atom_col_buffer(:panel_size, :), &
184 : ri_rs%grid_points(:, panel_start:panel_start + panel_size - 1), &
185 : ao_basis, ri_rs%particle_set(iatom)%r, ri_rs%cell, &
186 59 : cutoff_squared=cutoff_squared, periodic_images=periodic_images)
187 881 : DO i_blk = first_block, last_block
188 763 : row_first = blk_row_start(i_blk) - panel_start + 1
189 763 : row_last = row_first + r_blk_sizes(i_blk) - 1
190 82708 : IF (MAXVAL(ABS(atom_col_buffer(row_first:row_last, :))) > eps_filter) THEN
191 : CALL dbcsr_put_block(ri_rs%mat_phi_mu_l, row=i_blk, col=iatom, &
192 751 : block=atom_col_buffer(row_first:row_last, :))
193 : END IF
194 : END DO
195 : END DO
196 59 : DEALLOCATE (atom_col_buffer)
197 : END DO
198 48 : CALL dbcsr_finalize(ri_rs%mat_phi_mu_l)
199 48 : DEALLOCATE (r_blk_sizes, row_dist, col_dist)
200 48 : CALL dbcsr_distribution_release(dist)
201 48 : CALL timestop(handle)
202 240 : END SUBROUTINE atomic_basis_at_grid_point
203 :
204 : ! **************************************************************************************************
205 : !> \brief Creates a square topology from a selected row or column layout.
206 : !> \param matrix_template source matrix
207 : !> \param dim_type ROW or COL
208 : !> \param square_dist output distribution
209 : !> \param blk_sizes borrowed block sizes
210 : !> \param mapped_dist owned mapping in the other process-grid direction
211 : ! **************************************************************************************************
212 36060 : SUBROUTINE setup_square_topology(matrix_template, dim_type, square_dist, blk_sizes, mapped_dist)
213 :
214 : TYPE(dbcsr_type), INTENT(IN) :: matrix_template
215 : CHARACTER(LEN=*), INTENT(IN) :: dim_type
216 : TYPE(dbcsr_distribution_type), INTENT(OUT) :: square_dist
217 : INTEGER, DIMENSION(:), INTENT(OUT), POINTER :: blk_sizes, mapped_dist
218 :
219 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_square_topology'
220 :
221 : INTEGER :: handle, i, np, npcols, nprows
222 12020 : INTEGER, DIMENSION(:), POINTER :: col_blk, col_dist, row_blk, row_dist
223 : TYPE(dbcsr_distribution_type) :: dist_template
224 :
225 12020 : CALL timeset(routineN, handle)
226 :
227 : CALL dbcsr_get_info(matrix_template, distribution=dist_template, &
228 12020 : row_blk_size=row_blk, col_blk_size=col_blk)
229 : CALL dbcsr_distribution_get(dist_template, row_dist=row_dist, col_dist=col_dist, &
230 12020 : nprows=nprows, npcols=npcols)
231 :
232 12020 : IF (TRIM(dim_type) == 'ROW') THEN
233 : ! Creates ROW x ROW (e.g., Grid x Grid from mat_phi_mu_l)
234 600 : blk_sizes => row_blk
235 600 : np = npcols
236 1800 : ALLOCATE (mapped_dist(SIZE(blk_sizes)))
237 6960 : DO i = 1, SIZE(blk_sizes)
238 6960 : mapped_dist(i) = MOD(i - 1, np)
239 : END DO
240 : CALL dbcsr_distribution_new(square_dist, template=dist_template, &
241 600 : row_dist=row_dist, col_dist=mapped_dist)
242 :
243 11420 : ELSE IF (TRIM(dim_type) == 'COL') THEN
244 : ! Creates COL x COL (e.g., Aux x Aux from mat_Z_lP)
245 11420 : blk_sizes => col_blk
246 11420 : np = nprows
247 34260 : ALLOCATE (mapped_dist(SIZE(blk_sizes)))
248 35414 : DO i = 1, SIZE(blk_sizes)
249 35414 : mapped_dist(i) = MOD(i - 1, np)
250 : END DO
251 : CALL dbcsr_distribution_new(square_dist, template=dist_template, &
252 11420 : row_dist=mapped_dist, col_dist=col_dist)
253 : END IF
254 :
255 12020 : CALL timestop(handle)
256 :
257 12020 : END SUBROUTINE setup_square_topology
258 :
259 : ! **************************************************************************************************
260 : !> \brief Releases RI-RS topology mappings and temporary DBCSR matrices.
261 : !> \param dist distribution to release
262 : !> \param mapped_dist owned mapping to deallocate
263 : !> \param m1 optional matrix
264 : !> \param m2 optional matrix
265 : !> \param m3 optional matrix
266 : !> \param m4 optional matrix
267 : ! **************************************************************************************************
268 16620 : SUBROUTINE release_dbcsr_topology_and_matrices(dist, mapped_dist, m1, m2, m3, m4)
269 :
270 : TYPE(dbcsr_distribution_type), INTENT(INOUT), &
271 : OPTIONAL :: dist
272 : INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL, &
273 : POINTER :: mapped_dist
274 : TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: m1, m2, m3, m4
275 :
276 : CHARACTER(LEN=*), PARAMETER :: routineN = 'release_dbcsr_topology_and_matrices'
277 :
278 : INTEGER :: handle
279 :
280 16620 : CALL timeset(routineN, handle)
281 :
282 16620 : IF (PRESENT(dist)) CALL dbcsr_distribution_release(dist)
283 16620 : IF (PRESENT(mapped_dist)) THEN
284 12020 : IF (ASSOCIATED(mapped_dist)) THEN
285 12020 : DEALLOCATE (mapped_dist)
286 : NULLIFY (mapped_dist)
287 : END IF
288 : END IF
289 16620 : IF (PRESENT(m1)) CALL dbcsr_release(m1)
290 16620 : IF (PRESENT(m2)) CALL dbcsr_release(m2)
291 16620 : IF (PRESENT(m3)) CALL dbcsr_release(m3)
292 16620 : IF (PRESENT(m4)) CALL dbcsr_release(m4)
293 :
294 16620 : CALL timestop(handle)
295 :
296 16620 : END SUBROUTINE release_dbcsr_topology_and_matrices
297 :
298 : ! **************************************************************************************************
299 : !> \brief Reserves grid blocks with block-centre distances within a fixed radius.
300 : !> \param matrix target matrix
301 : !> \param centers block-centre coordinates
302 : !> \param radius truncation radius
303 : ! **************************************************************************************************
304 6 : SUBROUTINE reserve_blocks_within_radius(matrix, centers, radius)
305 :
306 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
307 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: centers
308 : REAL(KIND=dp), INTENT(IN) :: radius
309 :
310 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reserve_blocks_within_radius'
311 :
312 : INTEGER :: handle, i, j, my_pcol, my_prow, &
313 : nblkcols, nblkrows
314 6 : INTEGER, DIMENSION(:), POINTER :: col_blk, col_dist, row_blk, row_dist
315 : REAL(KIND=dp) :: radius2
316 6 : REAL(KIND=dp), ALLOCATABLE :: zero_blk(:, :)
317 : TYPE(dbcsr_distribution_type) :: dist
318 :
319 6 : CALL timeset(routineN, handle)
320 :
321 : CALL dbcsr_get_info(matrix, nblkrows_total=nblkrows, nblkcols_total=nblkcols, &
322 6 : row_blk_size=row_blk, col_blk_size=col_blk, distribution=dist)
323 : CALL dbcsr_distribution_get(dist, row_dist=row_dist, col_dist=col_dist, &
324 6 : myprow=my_prow, mypcol=my_pcol)
325 6 : CPASSERT(nblkrows == SIZE(centers, 2))
326 6 : CPASSERT(nblkcols == SIZE(centers, 2))
327 :
328 6 : radius2 = radius**2
329 144 : ALLOCATE (zero_blk(MAXVAL(row_blk(1:nblkrows)), MAXVAL(col_blk(1:nblkcols))))
330 6 : zero_blk(:, :) = 0.0_dp
331 :
332 66 : DO i = 1, nblkrows
333 60 : IF (row_dist(i) /= my_prow) CYCLE
334 336 : DO j = 1, nblkcols
335 300 : IF (col_dist(j) /= my_pcol) CYCLE
336 1260 : IF (DOT_PRODUCT(centers(:, i) - centers(:, j), centers(:, i) - centers(:, j)) <= radius2) THEN
337 168 : CALL dbcsr_put_block(matrix, i, j, zero_blk(1:row_blk(i), 1:col_blk(j)))
338 : END IF
339 : END DO
340 : END DO
341 6 : CALL dbcsr_finalize(matrix)
342 :
343 6 : DEALLOCATE (zero_blk)
344 6 : CALL timestop(handle)
345 :
346 6 : END SUBROUTINE reserve_blocks_within_radius
347 :
348 : ! **************************************************************************************************
349 : !> \brief Prints the percentage of non-zero elements in a distributed RI-RS matrix.
350 : !> \param matrix Distributed matrix whose occupation is reported
351 : !> \param label Mathematical matrix label used in the output
352 : !> \param para_env communicator shared by all matrix ranks
353 : !> \param unit_nr output unit; nonpositive values suppress printing
354 : !> \param msg_elem_limit maximum number of elements per DBCSR message
355 : !> \param suffix Optional text appended to the matrix label
356 : ! **************************************************************************************************
357 260 : SUBROUTINE print_matrix_occupation(matrix, label, para_env, unit_nr, msg_elem_limit, suffix)
358 :
359 : TYPE(dbcsr_type), INTENT(IN) :: matrix
360 : CHARACTER(LEN=*), INTENT(IN) :: label
361 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
362 : INTEGER, INTENT(IN) :: unit_nr
363 : INTEGER(KIND=int_8), INTENT(IN) :: msg_elem_limit
364 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: suffix
365 :
366 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_matrix_occupation'
367 :
368 : CHARACTER(LEN=32) :: output_format
369 : CHARACTER(LEN=max_line_length) :: msg, output_label
370 : INTEGER :: handle, i, unicode_shift
371 : REAL(KIND=dp) :: frac_2p31, max_loc, occ
372 :
373 260 : CALL timeset(routineN, handle)
374 :
375 260 : occ = dbcsr_get_occupation(matrix)
376 260 : max_loc = REAL(dbcsr_get_data_size(matrix), dp)
377 260 : CALL para_env%max(max_loc)
378 :
379 260 : IF (unit_nr > 0) THEN
380 130 : frac_2p31 = max_loc/REAL(msg_elem_limit, dp)
381 130 : output_label = 'Percentage of non-zero matrix elements in '//TRIM(label)
382 130 : IF (PRESENT(suffix)) output_label = TRIM(output_label)//TRIM(suffix)
383 : ! Fortran counts UTF-8 bytes, whereas the terminal displays each Greek letter in one
384 : ! column. Shift the absolute output tab once for every continuation byte.
385 130 : unicode_shift = 0
386 8171 : DO i = 1, LEN_TRIM(output_label)
387 8171 : IF (IAND(IACHAR(output_label(i:i)), 192) == 128) unicode_shift = unicode_shift + 1
388 : END DO
389 130 : WRITE (output_format, '(A,I0,A)') '(T2,A,T', 72 + unicode_shift, ',F7.2,A)'
390 130 : WRITE (unit_nr, output_format) TRIM(output_label), occ*100.0_dp, ' %'
391 130 : IF (frac_2p31 > 0.5_dp) THEN
392 : WRITE (msg, '(3A,F0.2,A)') &
393 0 : "The largest per-rank message of ", TRIM(label), " reaches ", frac_2p31, &
394 : " of the 32-bit limit that DBCSR uses for its message length. Beyond it the "// &
395 : "length overflows and multiply_cannon fails. Reduce the per-rank block size, "// &
396 0 : "for instance with more MPI ranks."
397 0 : CPWARN(TRIM(msg))
398 : END IF
399 130 : CALL m_flush(unit_nr)
400 : END IF
401 :
402 260 : CALL timestop(handle)
403 :
404 260 : END SUBROUTINE print_matrix_occupation
405 :
406 : ! **************************************************************************************************
407 : !> \brief Marks the grid blocks whose centroid lies within cutoff of the bounding box of the
408 : !> panel [blk0, blk1]'s chunk centroids.
409 : !> \param centroids ...
410 : !> \param blk0 ...
411 : !> \param blk1 ...
412 : !> \param cutoff ...
413 : !> \param used ...
414 : ! **************************************************************************************************
415 46 : SUBROUTINE mask_grid_blocks_near_panel(centroids, blk0, blk1, cutoff, used)
416 :
417 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: centroids
418 : INTEGER, INTENT(IN) :: blk0, blk1
419 : REAL(KIND=dp), INTENT(IN) :: cutoff
420 : LOGICAL, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: used
421 :
422 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mask_grid_blocks_near_panel'
423 :
424 : INTEGER :: c, handle, k
425 : REAL(KIND=dp) :: cutoff2, d2, dx
426 : REAL(KIND=dp), DIMENSION(3) :: hi, lo
427 :
428 46 : CALL timeset(routineN, handle)
429 :
430 46 : cutoff2 = cutoff**2
431 46 : lo(:) = MINVAL(centroids(:, blk0:blk1), DIM=2)
432 46 : hi(:) = MAXVAL(centroids(:, blk0:blk1), DIM=2)
433 :
434 138 : ALLOCATE (used(SIZE(centroids, 2)))
435 1012 : DO c = 1, SIZE(centroids, 2)
436 : d2 = 0.0_dp
437 3864 : DO k = 1, 3
438 2898 : dx = MAX(0.0_dp, lo(k) - centroids(k, c), centroids(k, c) - hi(k))
439 3864 : d2 = d2 + dx*dx
440 : END DO
441 1012 : used(c) = (d2 <= cutoff2)
442 : END DO
443 :
444 46 : CALL timestop(handle)
445 :
446 46 : END SUBROUTINE mask_grid_blocks_near_panel
447 :
448 : ! **************************************************************************************************
449 : !> \brief Exact allocated-element count of the geo template of panel [blk0, blk1]: the very same
450 : !> per-block-pair centroid test as build_geo_template_panel, so this is the true DBCSR
451 : !> data size of A_pan/B_pan (DBCSR stores whole blocks).
452 : !> \param r_blk_sizes ...
453 : !> \param centroids ...
454 : !> \param used ...
455 : !> \param blk0 ...
456 : !> \param blk1 ...
457 : !> \param cutoff ...
458 : !> \param nze_tmpl ...
459 : ! **************************************************************************************************
460 4 : SUBROUTINE panel_template_elems(r_blk_sizes, centroids, used, blk0, blk1, cutoff, nze_tmpl)
461 :
462 : INTEGER, DIMENSION(:), INTENT(IN) :: r_blk_sizes
463 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: centroids
464 : LOGICAL, DIMENSION(:), INTENT(IN) :: used
465 : INTEGER, INTENT(IN) :: blk0, blk1
466 : REAL(KIND=dp), INTENT(IN) :: cutoff
467 : INTEGER(KIND=int_8), INTENT(OUT) :: nze_tmpl
468 :
469 : CHARACTER(LEN=*), PARAMETER :: routineN = 'panel_template_elems'
470 :
471 : INTEGER :: c, handle, ib, n_used
472 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: used_idx
473 : REAL(KIND=dp) :: cutoff2
474 :
475 4 : CALL timeset(routineN, handle)
476 :
477 : ! Compress the near mask once so the pair loop only visits candidate columns.
478 88 : n_used = COUNT(used)
479 12 : ALLOCATE (used_idx(n_used))
480 4 : n_used = 0
481 88 : DO c = 1, SIZE(used)
482 88 : IF (used(c)) THEN
483 84 : n_used = n_used + 1
484 84 : used_idx(n_used) = c
485 : END IF
486 : END DO
487 :
488 4 : cutoff2 = cutoff**2
489 4 : nze_tmpl = 0_int_8
490 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(blk0, blk1, n_used, used_idx, centroids, cutoff2, &
491 4 : !$OMP r_blk_sizes) PRIVATE(ib, c) REDUCTION(+:nze_tmpl)
492 : DO ib = blk0, blk1
493 : DO c = 1, n_used
494 : IF (SUM((centroids(:, ib) - centroids(:, used_idx(c)))**2) <= cutoff2) THEN
495 : nze_tmpl = nze_tmpl + INT(r_blk_sizes(ib), int_8)*INT(r_blk_sizes(used_idx(c)), int_8)
496 : END IF
497 : END DO
498 : END DO
499 : !$OMP END PARALLEL DO
500 :
501 4 : CALL timestop(handle)
502 :
503 8 : END SUBROUTINE panel_template_elems
504 :
505 : ! **************************************************************************************************
506 : !> \brief Per-rank peak memory (GB) of one panel step of the neighborhood-restricted
507 : !> contractions: up to three grid x grid panels (two A operands and the shared B operand)
508 : !> plus the grid x RI intermediates (tmp2 and the accumulation operand) and the
509 : !> grid x AO intermediate (tmpA), whose column support is the panel's geometric
510 : !> neighborhood fraction f_near = width/n_grid. Used by planning and memory reporting.
511 : !> \param nze_tmpl ...
512 : !> \param pan_rows ...
513 : !> \param width ...
514 : !> \param n_grid_total ...
515 : !> \param n_RI ...
516 : !> \param n_ao ...
517 : !> \param n_procs ...
518 : !> \param mem_GB ...
519 : ! **************************************************************************************************
520 42 : SUBROUTINE panel_mem_estimate_GB(nze_tmpl, pan_rows, width, n_grid_total, n_RI, n_ao, &
521 : n_procs, mem_GB)
522 :
523 : INTEGER(KIND=int_8), INTENT(IN) :: nze_tmpl
524 : INTEGER, INTENT(IN) :: pan_rows, width, n_grid_total, n_RI, &
525 : n_ao, n_procs
526 : REAL(KIND=dp), INTENT(OUT) :: mem_GB
527 :
528 : REAL(KIND=dp) :: f_near
529 :
530 42 : f_near = REAL(width, dp)/REAL(MAX(n_grid_total, 1), dp)
531 : mem_GB = (3.0_dp*REAL(nze_tmpl, dp) + &
532 : REAL(pan_rows, dp)*f_near*(2.0_dp*REAL(n_RI, dp) + REAL(n_ao, dp)))* &
533 42 : 8.0_dp/REAL(MAX(n_procs, 1), dp)*1.0E-9_dp
534 :
535 42 : END SUBROUTINE panel_mem_estimate_GB
536 :
537 : ! **************************************************************************************************
538 : !> \brief Plans whole-grid-block panels using a message-size target without a cutoff.
539 : !> A positive cutoff additionally checks estimated message and memory costs.
540 : !> A single indivisible cutoff block exceeding an estimated limit is reported
541 : !> and accepted, since the panel cannot be split further.
542 : !> \param msg_elem_limit maximum number of elements per DBCSR message
543 : !> \param r_blk_sizes grid row-block sizes
544 : !> \param min_dim smaller dimension of the DBCSR process grid
545 : !> \param pan_first first grid block of each panel
546 : !> \param pan_last last grid block of each panel
547 : !> \param n_RI auxiliary basis dimension
548 : !> \param n_ao AO basis dimension
549 : !> \param n_procs number of MPI processes
550 : !> \param mem_budget_GB per-rank panel memory budget; nonpositive disables the memory check
551 : !> \param centroids optional grid-block centres
552 : !> \param cutoff optional geometric truncation radius
553 : ! **************************************************************************************************
554 40 : SUBROUTINE plan_grid_panels(msg_elem_limit, r_blk_sizes, min_dim, pan_first, pan_last, &
555 40 : n_RI, n_ao, n_procs, mem_budget_GB, centroids, cutoff)
556 :
557 : INTEGER(KIND=int_8), INTENT(IN) :: msg_elem_limit
558 : INTEGER, DIMENSION(:), INTENT(IN) :: r_blk_sizes
559 : INTEGER, INTENT(IN) :: min_dim
560 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: pan_first, pan_last
561 : INTEGER, INTENT(IN) :: n_RI, n_ao, n_procs
562 : REAL(KIND=dp), INTENT(IN) :: mem_budget_GB
563 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
564 : OPTIONAL :: centroids
565 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff
566 :
567 : CHARACTER(LEN=*), PARAMETER :: routineN = 'plan_grid_panels'
568 :
569 : INTEGER :: blk0, blk1, handle, ib, n_grid_blocks, &
570 : n_grid_total, n_panels, panel_size, &
571 : rows_acc, TARGET, width
572 : INTEGER(KIND=int_8) :: msg, msg_limit, nze_tmpl, side
573 40 : INTEGER, ALLOCATABLE :: tmp_first(:), tmp_last(:)
574 : LOGICAL :: fits, use_cutoff, warned_memory, &
575 : warned_message
576 40 : LOGICAL, ALLOCATABLE :: used(:)
577 : REAL(KIND=dp) :: f_near, mem_GB
578 :
579 40 : CALL timeset(routineN, handle)
580 :
581 40 : use_cutoff = PRESENT(centroids) .AND. PRESENT(cutoff)
582 40 : IF (use_cutoff) use_cutoff = cutoff > 0.0_dp
583 40 : msg_limit = msg_elem_limit/4
584 40 : warned_memory = .FALSE.
585 40 : warned_message = .FALSE.
586 40 : n_grid_blocks = SIZE(r_blk_sizes)
587 570 : n_grid_total = SUM(r_blk_sizes)
588 40 : panel_size = n_grid_total
589 40 : IF (.NOT. use_cutoff) THEN
590 : ! Preserve the original message-derived target and whole-block growth.
591 38 : panel_size = INT(0.5_dp*REAL(msg_elem_limit, dp)*REAL(min_dim, dp)/REAL(n_grid_total, dp))
592 38 : panel_size = MAX(1, MIN(panel_size, n_grid_total))
593 : END IF
594 160 : ALLOCATE (tmp_first(n_grid_blocks), tmp_last(n_grid_blocks))
595 :
596 40 : n_panels = 0
597 40 : blk0 = 1
598 80 : DO WHILE (blk0 <= n_grid_blocks)
599 : TARGET = panel_size
600 2 : DO
601 40 : rows_acc = 0
602 530 : DO ib = blk0, n_grid_blocks
603 530 : rows_acc = rows_acc + r_blk_sizes(ib)
604 530 : blk1 = ib
605 530 : IF (rows_acc >= TARGET) EXIT
606 : END DO
607 40 : IF (.NOT. use_cutoff) EXIT
608 2 : CALL mask_grid_blocks_near_panel(centroids, blk0, blk1, cutoff, used)
609 44 : width = SUM(r_blk_sizes, MASK=used)
610 2 : CALL panel_template_elems(r_blk_sizes, centroids, used, blk0, blk1, cutoff, nze_tmpl)
611 2 : f_near = REAL(width, dp)/REAL(MAX(n_grid_total, 1), dp)
612 2 : side = INT(REAL(rows_acc, dp)*f_near*REAL(MAX(n_RI, n_ao), dp), int_8)
613 2 : msg = MAX(nze_tmpl, side)/INT(MAX(min_dim, 1), int_8)
614 : CALL panel_mem_estimate_GB(nze_tmpl, rows_acc, width, n_grid_total, &
615 2 : n_RI, n_ao, n_procs, mem_GB)
616 2 : fits = msg <= msg_limit
617 2 : IF (mem_budget_GB > 0.0_dp) fits = fits .AND. mem_GB <= mem_budget_GB
618 0 : IF (fits) EXIT
619 0 : IF (blk1 == blk0) THEN
620 0 : IF (msg > msg_limit .AND. .NOT. warned_message) THEN
621 0 : CPWARN("A single RI-RS grid block exceeds the estimated panel message limit.")
622 0 : warned_message = .TRUE.
623 : END IF
624 0 : IF (mem_budget_GB > 0.0_dp .AND. mem_GB > mem_budget_GB .AND. .NOT. warned_memory) THEN
625 0 : CPWARN("A single RI-RS grid block exceeds the estimated panel memory budget.")
626 0 : warned_memory = .TRUE.
627 : END IF
628 : EXIT
629 : END IF
630 : ! Preserve whole blocks and prefer large panels to reduce the number of multiplies.
631 38 : TARGET = MAX(1, MIN(TARGET, rows_acc)/2)
632 : END DO
633 40 : n_panels = n_panels + 1
634 40 : tmp_first(n_panels) = blk0
635 40 : tmp_last(n_panels) = blk1
636 40 : blk0 = blk1 + 1
637 : END DO
638 :
639 160 : ALLOCATE (pan_first, SOURCE=tmp_first(:n_panels))
640 120 : ALLOCATE (pan_last, SOURCE=tmp_last(:n_panels))
641 40 : CALL timestop(handle)
642 :
643 80 : END SUBROUTINE plan_grid_panels
644 :
645 : ! **************************************************************************************************
646 : !> \brief Creates an empty (panel_chunks x neighborhood_chunks) DBCSR matrix with zero blocks
647 : !> pre-allocated only where |centroid(panel_row r) - centroid(column c)| <= cutoff.
648 : !> Used with retain_sparsity=.TRUE. in the subsequent dbcsr_multiply so distant blocks
649 : !> of the grid-basis panels (ϕ G ϕ^T, Z W Z^T, ...) are never computed at all.
650 : !> \param L_pan ...
651 : !> \param L_full ...
652 : !> \param centroids ...
653 : !> \param cutoff ...
654 : !> \param blk0 ...
655 : !> \param A_template ...
656 : !> \param col_map ...
657 : ! **************************************************************************************************
658 62 : SUBROUTINE build_geo_template_panel(L_pan, L_full, centroids, cutoff, blk0, A_template, col_map)
659 : TYPE(dbcsr_type), INTENT(IN) :: L_pan, L_full
660 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: centroids
661 : REAL(KIND=dp), INTENT(IN) :: cutoff
662 : INTEGER, INTENT(IN) :: blk0
663 : TYPE(dbcsr_type), INTENT(OUT) :: A_template
664 : INTEGER, DIMENSION(:), INTENT(IN) :: col_map
665 :
666 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_geo_template_panel'
667 :
668 : INTEGER :: c, cg, cs, handle, my_pcol, my_prow, &
669 : n_grid_blks, n_pan_blks, npcols, &
670 : nprows, r, rs
671 62 : INTEGER, DIMENSION(:), POINTER :: grid_blk_sizes, pan_blk_sizes
672 : REAL(KIND=dp) :: cutoff2
673 62 : REAL(KIND=dp), ALLOCATABLE :: zero_blk(:, :)
674 : TYPE(dbcsr_distribution_type) :: dist
675 :
676 62 : CALL timeset(routineN, handle)
677 :
678 62 : cutoff2 = cutoff**2
679 62 : CALL dbcsr_get_info(L_pan, nblkrows_total=n_pan_blks, row_blk_size=pan_blk_sizes)
680 62 : CALL dbcsr_get_info(L_full, nblkrows_total=n_grid_blks, row_blk_size=grid_blk_sizes)
681 :
682 : ! create_product_matrix assigns row r to process MOD(r-1,nprows) and
683 : ! col c to MOD(c-1,npcols), so we can determine local ownership analytically.
684 62 : CALL create_product_matrix(L_pan, L_full, 'T', A_template)
685 62 : CALL dbcsr_get_info(A_template, distribution=dist)
686 : CALL dbcsr_distribution_get(dist, nprows=nprows, npcols=npcols, &
687 62 : myprow=my_prow, mypcol=my_pcol)
688 :
689 : ALLOCATE (zero_blk(MAXVAL(pan_blk_sizes(1:n_pan_blks)), &
690 2852 : MAXVAL(grid_blk_sizes(1:n_grid_blks))))
691 62 : zero_blk(:, :) = 0.0_dp
692 :
693 1364 : DO r = 1, n_pan_blks
694 1302 : IF (MOD(r - 1, nprows) /= my_prow) CYCLE
695 651 : rs = pan_blk_sizes(r)
696 14384 : DO c = 1, n_grid_blks
697 13671 : IF (MOD(c - 1, npcols) /= my_pcol) CYCLE
698 13671 : cg = col_map(c)
699 13671 : IF ((centroids(1, blk0 + r - 1) - centroids(1, cg))**2 + &
700 : (centroids(2, blk0 + r - 1) - centroids(2, cg))**2 + &
701 1302 : (centroids(3, blk0 + r - 1) - centroids(3, cg))**2 <= cutoff2) THEN
702 12741 : cs = grid_blk_sizes(c)
703 12741 : CALL dbcsr_put_block(A_template, r, c, zero_blk(1:rs, 1:cs))
704 : END IF
705 : END DO
706 : END DO
707 62 : CALL dbcsr_finalize(A_template)
708 :
709 62 : DEALLOCATE (zero_blk)
710 62 : CALL timestop(handle)
711 :
712 62 : END SUBROUTINE build_geo_template_panel
713 :
714 : ! **************************************************************************************************
715 : !> \brief Slices a contiguous range of grid row-blocks [blk0, blk1] out of a (grid x n) DBCSR
716 : !> matrix into a new (P x n) panel matrix: iterate the source's local blocks, put the
717 : !> in-range ones into the panel with a remapped row-block index, then finalize. Row-block
718 : !> index i of the panel corresponds to source row-block blk0+i-1.
719 : !> \param mat_full ...
720 : !> \param blk0 ...
721 : !> \param blk1 ...
722 : !> \param mat_panel ...
723 : ! **************************************************************************************************
724 3012 : SUBROUTINE extract_grid_panel(mat_full, blk0, blk1, mat_panel)
725 :
726 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_full
727 : INTEGER, INTENT(IN) :: blk0, blk1
728 : TYPE(dbcsr_type), INTENT(OUT) :: mat_panel
729 :
730 : CHARACTER(LEN=*), PARAMETER :: routineN = 'extract_grid_panel'
731 :
732 : INTEGER :: handle, ib, jb, npb
733 3012 : INTEGER, DIMENSION(:), POINTER :: col_blk_full, col_dist_full, &
734 3012 : row_blk_full, row_blk_pan, &
735 3012 : row_dist_full, row_dist_pan
736 3012 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk
737 : TYPE(dbcsr_distribution_type) :: dist_full, dist_pan
738 : TYPE(dbcsr_iterator_type) :: iter
739 :
740 3012 : CALL timeset(routineN, handle)
741 :
742 : CALL dbcsr_get_info(mat_full, distribution=dist_full, &
743 3012 : row_blk_size=row_blk_full, col_blk_size=col_blk_full)
744 3012 : CALL dbcsr_distribution_get(dist_full, row_dist=row_dist_full, col_dist=col_dist_full)
745 :
746 3012 : npb = blk1 - blk0 + 1
747 12048 : ALLOCATE (row_dist_pan(npb), row_blk_pan(npb))
748 92352 : row_dist_pan(:) = row_dist_full(blk0:blk1)
749 92352 : row_blk_pan(:) = row_blk_full(blk0:blk1)
750 :
751 : CALL dbcsr_distribution_new(dist_pan, template=dist_full, &
752 3012 : row_dist=row_dist_pan, col_dist=col_dist_full)
753 : CALL dbcsr_create(mat_panel, name="grid_panel", dist=dist_pan, &
754 : matrix_type=dbcsr_type_no_symmetry, &
755 3012 : row_blk_size=row_blk_pan, col_blk_size=col_blk_full)
756 :
757 3012 : CALL dbcsr_iterator_start(iter, mat_full)
758 54093 : DO WHILE (dbcsr_iterator_blocks_left(iter))
759 51081 : CALL dbcsr_iterator_next_block(iter, ib, jb, blk)
760 51081 : IF (ib < blk0 .OR. ib > blk1) CYCLE
761 54093 : CALL dbcsr_put_block(mat_panel, ib - blk0 + 1, jb, blk)
762 : END DO
763 3012 : CALL dbcsr_iterator_stop(iter)
764 3012 : CALL dbcsr_finalize(mat_panel)
765 :
766 3012 : CALL dbcsr_distribution_release(dist_pan)
767 3012 : DEALLOCATE (row_dist_pan, row_blk_pan)
768 :
769 3012 : CALL timestop(handle)
770 :
771 9036 : END SUBROUTINE extract_grid_panel
772 :
773 : ! **************************************************************************************************
774 : !> \brief Marks which column blocks of a DBCSR matrix carry at least one non-zero block anywhere
775 : !> (global union). Used to restrict the inner index of the panel multiplies to the
776 : !> AO/RI atoms that actually touch the panel (exact: dropped rows only meet zeros).
777 : !> \param matrix ...
778 : !> \param para_env ...
779 : !> \param used ...
780 : ! **************************************************************************************************
781 2382 : SUBROUTINE collect_used_col_blocks(matrix, para_env, used)
782 :
783 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
784 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
785 : LOGICAL, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: used
786 :
787 : CHARACTER(LEN=*), PARAMETER :: routineN = 'collect_used_col_blocks'
788 :
789 : INTEGER :: handle, ib, jb, nblkcols
790 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iused
791 2382 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk
792 : TYPE(dbcsr_iterator_type) :: iter
793 :
794 2382 : CALL timeset(routineN, handle)
795 :
796 2382 : CALL dbcsr_get_info(matrix, nblkcols_total=nblkcols)
797 7146 : ALLOCATE (iused(nblkcols))
798 2382 : iused(:) = 0
799 :
800 2382 : CALL dbcsr_iterator_start(iter, matrix)
801 43648 : DO WHILE (dbcsr_iterator_blocks_left(iter))
802 41266 : CALL dbcsr_iterator_next_block(iter, ib, jb, blk)
803 41266 : iused(jb) = 1
804 : END DO
805 2382 : CALL dbcsr_iterator_stop(iter)
806 :
807 2382 : CALL para_env%sum(iused)
808 :
809 4764 : ALLOCATE (used(nblkcols))
810 8040 : used(:) = (iused(:) > 0)
811 2382 : DEALLOCATE (iused)
812 :
813 2382 : CALL timestop(handle)
814 :
815 4764 : END SUBROUTINE collect_used_col_blocks
816 :
817 : ! **************************************************************************************************
818 : !> \brief Copies the flagged block rows (compress_rows=.TRUE.) or block columns (.FALSE.) of a
819 : !> DBCSR matrix into a compressed matrix. The subset keeps the parent's process assignment
820 : !> along the compressed dimension, so every block stays on its owning rank: the extraction
821 : !> is purely local (zero communication), like extract_grid_panel.
822 : !> \param mat_full ...
823 : !> \param used ...
824 : !> \param mat_out ...
825 : !> \param compress_rows ...
826 : !> \param blk_map ...
827 : ! **************************************************************************************************
828 6308 : SUBROUTINE extract_masked_blocks(mat_full, used, mat_out, compress_rows, blk_map)
829 :
830 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_full
831 : LOGICAL, DIMENSION(:), INTENT(IN) :: used
832 : TYPE(dbcsr_type), INTENT(OUT) :: mat_out
833 : LOGICAL, INTENT(IN) :: compress_rows
834 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT), &
835 : OPTIONAL :: blk_map
836 :
837 : CHARACTER(LEN=*), PARAMETER :: routineN = 'extract_masked_blocks'
838 :
839 : INTEGER :: handle, ib, jb, n_blk, n_sub, r
840 6308 : INTEGER, ALLOCATABLE, DIMENSION(:) :: inv_map
841 6308 : INTEGER, DIMENSION(:), POINTER :: blk_full, blk_sub, col_blk_full, &
842 6308 : col_dist_full, dist_full_1d, &
843 6308 : dist_sub_1d, row_blk_full, &
844 6308 : row_dist_full
845 6308 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk
846 : TYPE(dbcsr_distribution_type) :: dist_full, dist_sub
847 : TYPE(dbcsr_iterator_type) :: iter
848 :
849 6308 : CALL timeset(routineN, handle)
850 :
851 : CALL dbcsr_get_info(mat_full, distribution=dist_full, &
852 6308 : row_blk_size=row_blk_full, col_blk_size=col_blk_full)
853 6308 : CALL dbcsr_distribution_get(dist_full, row_dist=row_dist_full, col_dist=col_dist_full)
854 :
855 6308 : IF (compress_rows) THEN
856 3926 : blk_full => row_blk_full
857 3926 : dist_full_1d => row_dist_full
858 : ELSE
859 2382 : blk_full => col_blk_full
860 2382 : dist_full_1d => col_dist_full
861 : END IF
862 6308 : n_blk = SIZE(blk_full)
863 :
864 22848 : n_sub = COUNT(used)
865 6308 : CPASSERT(n_sub > 0)
866 37848 : ALLOCATE (inv_map(n_blk), blk_sub(n_sub), dist_sub_1d(n_sub))
867 6350 : IF (PRESENT(blk_map)) ALLOCATE (blk_map(n_sub))
868 6308 : inv_map(:) = 0
869 6308 : r = 0
870 22848 : DO ib = 1, n_blk
871 22848 : IF (used(ib)) THEN
872 16540 : r = r + 1
873 16540 : inv_map(ib) = r
874 16540 : blk_sub(r) = blk_full(ib)
875 16540 : dist_sub_1d(r) = dist_full_1d(ib)
876 16540 : IF (PRESENT(blk_map)) blk_map(r) = ib
877 : END IF
878 : END DO
879 :
880 6308 : IF (compress_rows) THEN
881 : CALL dbcsr_distribution_new(dist_sub, template=dist_full, &
882 3926 : row_dist=dist_sub_1d, col_dist=col_dist_full)
883 : CALL dbcsr_create(mat_out, name="row_subset", dist=dist_sub, &
884 : matrix_type=dbcsr_type_no_symmetry, &
885 3926 : row_blk_size=blk_sub, col_blk_size=col_blk_full)
886 : ELSE
887 : CALL dbcsr_distribution_new(dist_sub, template=dist_full, &
888 2382 : row_dist=row_dist_full, col_dist=dist_sub_1d)
889 : CALL dbcsr_create(mat_out, name="col_subset", dist=dist_sub, &
890 : matrix_type=dbcsr_type_no_symmetry, &
891 2382 : row_blk_size=row_blk_full, col_blk_size=blk_sub)
892 : END IF
893 :
894 6308 : CALL dbcsr_iterator_start(iter, mat_full)
895 61048 : DO WHILE (dbcsr_iterator_blocks_left(iter))
896 54740 : CALL dbcsr_iterator_next_block(iter, ib, jb, blk)
897 61048 : IF (compress_rows) THEN
898 13474 : IF (inv_map(ib) > 0) CALL dbcsr_put_block(mat_out, inv_map(ib), jb, blk)
899 : ELSE
900 41266 : IF (inv_map(jb) > 0) CALL dbcsr_put_block(mat_out, ib, inv_map(jb), blk)
901 : END IF
902 : END DO
903 6308 : CALL dbcsr_iterator_stop(iter)
904 6308 : CALL dbcsr_finalize(mat_out)
905 :
906 6308 : CALL dbcsr_distribution_release(dist_sub)
907 6308 : DEALLOCATE (inv_map, blk_sub, dist_sub_1d)
908 :
909 6308 : CALL timestop(handle)
910 :
911 12616 : END SUBROUTINE extract_masked_blocks
912 :
913 : ! **************************************************************************************************
914 : !> \brief Creates the (empty) result matrix of op(mat_left) * op(mat_right) with the correct block
915 : !> structure and a distribution on the shared process grid, ready to be filled by
916 : !> dbcsr_multiply. Row structure comes from op(left), column structure from op(right).
917 : !> \param mat_left ...
918 : !> \param mat_right ...
919 : !> \param transb 'N' or 'T' applied to mat_right
920 : !> \param mat_out ...
921 : ! **************************************************************************************************
922 8514 : SUBROUTINE create_product_matrix(mat_left, mat_right, transb, mat_out)
923 :
924 : TYPE(dbcsr_type), INTENT(IN) :: mat_left, mat_right
925 : CHARACTER(LEN=1), INTENT(IN) :: transb
926 : TYPE(dbcsr_type), INTENT(OUT) :: mat_out
927 :
928 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_product_matrix'
929 :
930 : INTEGER :: handle, i, npcols, nprows
931 8514 : INTEGER, DIMENSION(:), POINTER :: col_blk_r, out_col_blk, out_col_dist, &
932 8514 : out_row_blk, out_row_dist, row_blk_l, &
933 8514 : row_blk_r
934 : TYPE(dbcsr_distribution_type) :: dist_l, dist_out
935 :
936 8514 : CALL timeset(routineN, handle)
937 :
938 8514 : CALL dbcsr_get_info(mat_left, distribution=dist_l, row_blk_size=row_blk_l)
939 8514 : CALL dbcsr_get_info(mat_right, row_blk_size=row_blk_r, col_blk_size=col_blk_r)
940 8514 : CALL dbcsr_distribution_get(dist_l, nprows=nprows, npcols=npcols)
941 :
942 : ! block SIZES follow op(left)/op(right); DISTRIBUTIONS are freshly round-robined onto the
943 : ! shared process grid (a transposed operand's row-dist is NOT a valid col-dist on a
944 : ! non-square grid). dbcsr_multiply redistributes internally, so any valid mapping works.
945 8514 : out_row_blk => row_blk_l
946 8514 : IF (transb == 'N') THEN
947 6178 : out_col_blk => col_blk_r
948 : ELSE
949 2336 : out_col_blk => row_blk_r
950 : END IF
951 :
952 42570 : ALLOCATE (out_row_dist(SIZE(out_row_blk)), out_col_dist(SIZE(out_col_blk)))
953 132732 : DO i = 1, SIZE(out_row_blk)
954 132732 : out_row_dist(i) = MOD(i - 1, nprows)
955 : END DO
956 57388 : DO i = 1, SIZE(out_col_blk)
957 57388 : out_col_dist(i) = MOD(i - 1, npcols)
958 : END DO
959 :
960 : CALL dbcsr_distribution_new(dist_out, template=dist_l, &
961 8514 : row_dist=out_row_dist, col_dist=out_col_dist)
962 : CALL dbcsr_create(mat_out, name="panel_product", dist=dist_out, &
963 : matrix_type=dbcsr_type_no_symmetry, &
964 8514 : row_blk_size=out_row_blk, col_blk_size=out_col_blk)
965 8514 : CALL dbcsr_distribution_release(dist_out)
966 8514 : DEALLOCATE (out_row_dist, out_col_dist)
967 :
968 8514 : CALL timestop(handle)
969 :
970 25542 : END SUBROUTINE create_product_matrix
971 :
972 : ! **************************************************************************************************
973 : !> \brief Projects two operators onto grid panels, takes their Hadamard product, and
974 : !> back-projects it: out += scale * L_out^T (A_grid ∘ B_grid) L_out,
975 : !> where A_grid = L_A M_A L_A^T and B_grid = L_B M_B L_B^T.
976 : !> Each L has shape (grid, basis), and its M has shape (basis, basis).
977 : !> The three L matrices must share the ordered grid rows, row blocking and process grid.
978 : !> Only panel x grid slabs are formed; contributions are accumulated over all panels.
979 : !> The sequence of EPS_FILTER screening and floating-point reductions may differ from
980 : !> a full-grid evaluation.
981 : !>
982 : !> The physical contractions and their signs/spin factors are specified in gw_ri_rs_kernels:
983 : !> χ: L_A=L_B=ϕ, M_A=G_occ, M_B=G_vir, L_out=Z -> RI x RI
984 : !> Σx: L_A=ϕ, M_A=D, L_B=Z, M_B=M^-1 V^tr M^-1, L_out=ϕ -> AO x AO
985 : !> Σc: L_A=ϕ, M_A=[G_occ,G_vir], L_B=Z, M_B=W, L_out=ϕ -> two AO x AO outputs.
986 : !> The second Σc output reuses the same W grid panel; its projection and back-projection
987 : !> otherwise follow the same operations as the first output.
988 : !> \param L_A grid-by-basis collocation or fitting matrix for the first operand
989 : !> \param M_A one or two square operators, with row/column blocking matching the columns of L_A
990 : !> \param L_B grid-by-basis collocation or fitting matrix for the second operand
991 : !> \param M_B second square operator, with row/column blocking matching the columns of L_B
992 : !> \param L_out grid-by-output-basis matrix for back-projection
993 : !> \param mat_out borrowed pointers to pre-created outputs, one per M_A, with the column blocking of L_out
994 : !> \param scale prefactor of the Hadamard product, including the caller's sign/spin convention
995 : !> \param eps DBCSR multiplication screening threshold
996 : !> \param para_env communicator containing all matrix ranks
997 : !> \param pan_first first global grid row-block of each ordered, nonoverlapping panel
998 : !> \param pan_last last global grid row-block of each panel; panels together cover the grid
999 : !> \param lb_eq_la L_B equals L_A in values and layout, allowing their grid/AO slices to be shared
1000 : !> \param lout_eq_la L_out equals L_A in values and layout, allowing back-projection slices to be shared
1001 : !> \param zero_out clear the output matrices before accumulating, otherwise add to their existing values
1002 : !> \param centroids optional global grid-block centres, paired with a positive cutoff for truncation
1003 : !> \param cutoff retain grid-block pairs within this centroid distance; nonpositive disables truncation
1004 : !> \param grid_occupation allocated-element fraction of the first grid operand before the Hadamard product
1005 : ! **************************************************************************************************
1006 1506 : SUBROUTINE contract_grid_panels(L_A, M_A, L_B, M_B, L_out, mat_out, scale, eps, para_env, &
1007 1506 : pan_first, pan_last, lb_eq_la, lout_eq_la, zero_out, &
1008 1506 : centroids, cutoff, grid_occupation)
1009 :
1010 : TYPE(dbcsr_type), INTENT(INOUT), TARGET :: L_A
1011 : TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: M_A
1012 : TYPE(dbcsr_type), INTENT(INOUT), TARGET :: L_B
1013 : TYPE(dbcsr_type), INTENT(INOUT) :: M_B
1014 : TYPE(dbcsr_type), INTENT(INOUT), TARGET :: L_out
1015 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: mat_out
1016 : REAL(KIND=dp), INTENT(IN) :: scale, eps
1017 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
1018 : INTEGER, DIMENSION(:), INTENT(IN) :: pan_first, pan_last
1019 : LOGICAL, INTENT(IN) :: lb_eq_la, lout_eq_la, zero_out
1020 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
1021 : OPTIONAL :: centroids
1022 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff
1023 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: grid_occupation
1024 :
1025 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_grid_panels'
1026 :
1027 : INTEGER :: blk0, blk1, handle, i_A, ipan, &
1028 : n_grid_total, ncols_pan, nrows_pan
1029 1506 : INTEGER, ALLOCATABLE, DIMENSION(:) :: gmap
1030 : LOGICAL :: use_cutoff
1031 1506 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: grid_used, usedA, usedB
1032 4518 : TYPE(dbcsr_type) :: A_pan(2), B_pan, LB_pan, LB_panC, &
1033 : MA_sub, MB_sub, tmp2, tmpA, tmpB
1034 : TYPE(dbcsr_type), POINTER :: panel_out, RB_A, RB_B, RB_out
1035 : TYPE(dbcsr_type), TARGET :: LA_near, LA_pan, LA_panC, LB_near, &
1036 : Lout_near, Lout_pan
1037 :
1038 1506 : CALL timeset(routineN, handle)
1039 :
1040 1506 : CPASSERT(SIZE(M_A) >= 1 .AND. SIZE(M_A) <= 2)
1041 1506 : CPASSERT(SIZE(mat_out) == SIZE(M_A))
1042 :
1043 1506 : use_cutoff = PRESENT(centroids) .AND. PRESENT(cutoff)
1044 1506 : IF (use_cutoff) use_cutoff = cutoff > 0.0_dp
1045 1506 : IF (PRESENT(grid_occupation)) grid_occupation = 0.0_dp
1046 :
1047 1506 : CALL dbcsr_get_info(L_A, nfullrows_total=n_grid_total)
1048 :
1049 1506 : IF (zero_out) THEN
1050 3682 : DO i_A = 1, SIZE(mat_out)
1051 3682 : CALL dbcsr_set(mat_out(i_A)%matrix, 0.0_dp)
1052 : END DO
1053 : END IF
1054 :
1055 3012 : DO ipan = 1, SIZE(pan_first)
1056 1506 : blk0 = pan_first(ipan)
1057 1506 : blk1 = pan_last(ipan)
1058 :
1059 : ! Slice only the current grid rows; compact away AO/RI columns that are exactly absent.
1060 1506 : CALL extract_grid_panel(L_A, blk0, blk1, LA_pan)
1061 1506 : CALL collect_used_col_blocks(LA_pan, para_env, usedA)
1062 1506 : IF (.NOT. lb_eq_la) THEN
1063 876 : CALL extract_grid_panel(L_B, blk0, blk1, LB_pan)
1064 876 : CALL collect_used_col_blocks(LB_pan, para_env, usedB)
1065 : ELSE
1066 630 : IF (ALLOCATED(usedB)) DEALLOCATE (usedB)
1067 3380 : ALLOCATE (usedB, SOURCE=usedA)
1068 : END IF
1069 1506 : IF (.NOT. ANY(usedA) .OR. .NOT. ANY(usedB)) THEN
1070 0 : CALL dbcsr_release(LA_pan)
1071 0 : IF (.NOT. lb_eq_la) CALL dbcsr_release(LB_pan)
1072 : CYCLE
1073 : END IF
1074 1506 : IF (lout_eq_la) THEN
1075 : panel_out => LA_pan
1076 : ELSE
1077 630 : CALL extract_grid_panel(L_out, blk0, blk1, Lout_pan)
1078 630 : panel_out => Lout_pan
1079 : END IF
1080 1506 : CALL extract_masked_blocks(LA_pan, usedA, LA_panC, compress_rows=.FALSE.)
1081 :
1082 : ! Grid rows within reach of the panel: with the CUTOFF_RADIUS_RL_W truncation only
1083 : ! they can appear as columns of the panel products / inner rows of the L_out multiply,
1084 : ! so the system-wide phi/Z right operands are cut down to this neighborhood slice
1085 : ! (local extraction, zero communication; exact w.r.t. the geo template).
1086 1506 : IF (use_cutoff) THEN
1087 42 : CALL mask_grid_blocks_near_panel(centroids, blk0, blk1, cutoff, grid_used)
1088 42 : CALL extract_masked_blocks(L_A, grid_used, LA_near, compress_rows=.TRUE., blk_map=gmap)
1089 42 : IF (.NOT. lb_eq_la) CALL extract_masked_blocks(L_B, grid_used, LB_near, compress_rows=.TRUE.)
1090 42 : IF (.NOT. lout_eq_la) CALL extract_masked_blocks(L_out, grid_used, Lout_near, compress_rows=.TRUE.)
1091 : RB_A => LA_near
1092 : ELSE
1093 : RB_A => L_A
1094 : END IF
1095 1506 : IF (lb_eq_la) THEN
1096 : RB_B => RB_A
1097 876 : ELSE IF (use_cutoff) THEN
1098 : RB_B => LB_near
1099 : ELSE
1100 854 : RB_B => L_B
1101 : END IF
1102 1506 : IF (lout_eq_la) THEN
1103 : RB_out => RB_A
1104 630 : ELSE IF (use_cutoff) THEN
1105 : RB_out => Lout_near
1106 : ELSE
1107 610 : RB_out => L_out
1108 : END IF
1109 :
1110 3842 : DO i_A = 1, SIZE(M_A)
1111 : ! A_pan(i_A) = LA_pan * M_A(i_A) * L_A^T (P x grid_near).
1112 : ! When cutoff is active, A_pan(i_A) is pre-seeded with only nearby blocks via
1113 : ! build_geo_template_panel, and the multiply uses retain_sparsity to skip
1114 : ! computing distant blocks entirely within the chosen grid cutoff.
1115 2336 : CALL extract_masked_blocks(M_A(i_A), usedA, MA_sub, compress_rows=.TRUE.)
1116 2336 : CALL create_product_matrix(LA_panC, MA_sub, 'N', tmpA)
1117 2336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, LA_panC, MA_sub, 0.0_dp, tmpA, filter_eps=eps)
1118 2336 : CALL dbcsr_release(MA_sub)
1119 2336 : IF (use_cutoff) THEN
1120 : CALL build_geo_template_panel(LA_pan, LA_near, centroids, cutoff, blk0, A_pan(i_A), &
1121 62 : col_map=gmap)
1122 : ELSE
1123 2274 : CALL create_product_matrix(tmpA, RB_A, 'T', A_pan(i_A))
1124 : END IF
1125 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpA, RB_A, 0.0_dp, A_pan(i_A), &
1126 2336 : filter_eps=eps, retain_sparsity=use_cutoff)
1127 2336 : CALL dbcsr_release(tmpA)
1128 :
1129 : ! Grid-basis occupation of A_pan(i_A) = ϕ G ϕ^T, accumulated over ALL panels into the
1130 : ! occupation of the (never formed) full grid x grid object:
1131 : ! sum_panels nnz(A_pan(i_A)) / n_grid^2, with nnz = occ * pan_rows * pan_cols.
1132 : ! Panel-independent by construction -- a single-panel sample would instead report the
1133 : ! local neighbor count of whichever region happens to land in that panel.
1134 6178 : IF (PRESENT(grid_occupation) .AND. i_A == 1) THEN
1135 46 : CALL dbcsr_get_info(A_pan(i_A), nfullrows_total=nrows_pan, nfullcols_total=ncols_pan)
1136 : grid_occupation = grid_occupation + dbcsr_get_occupation(A_pan(i_A))* &
1137 : REAL(ncols_pan, dp)*REAL(nrows_pan, dp)/ &
1138 46 : (REAL(n_grid_total, dp)*REAL(n_grid_total, dp))
1139 : END IF
1140 : END DO
1141 :
1142 : ! B_pan = LB_pan * M_B * L_B^T (P x grid_near); reuse the L_A slices when L_B == L_A.
1143 : ! B_pan retains the union of the A panels' block structures so that
1144 : ! retain_sparsity forces the final multiply to fill only those blocks (exact for ∘).
1145 1506 : IF (lb_eq_la) THEN
1146 630 : CALL extract_masked_blocks(M_B, usedA, MB_sub, compress_rows=.TRUE.)
1147 630 : CALL create_product_matrix(LA_panC, MB_sub, 'N', tmpB)
1148 630 : CALL dbcsr_multiply('N', 'N', 1.0_dp, LA_panC, MB_sub, 0.0_dp, tmpB, filter_eps=eps)
1149 : ELSE
1150 876 : CALL extract_masked_blocks(LB_pan, usedB, LB_panC, compress_rows=.FALSE.)
1151 876 : CALL extract_masked_blocks(M_B, usedB, MB_sub, compress_rows=.TRUE.)
1152 876 : CALL create_product_matrix(LB_panC, MB_sub, 'N', tmpB)
1153 876 : CALL dbcsr_multiply('N', 'N', 1.0_dp, LB_panC, MB_sub, 0.0_dp, tmpB, filter_eps=eps)
1154 876 : CALL dbcsr_release(LB_panC)
1155 : END IF
1156 1506 : CALL dbcsr_release(MB_sub)
1157 1506 : CALL dbcsr_create(B_pan, template=A_pan(1))
1158 1506 : CALL dbcsr_copy(B_pan, A_pan(1))
1159 1506 : IF (SIZE(M_A) == 2) CALL dbcsr_add(B_pan, A_pan(2), 1.0_dp, 1.0_dp)
1160 1506 : CALL dbcsr_set(B_pan, 0.0_dp)
1161 : ! The output pattern is already fixed. Omitting redundant on-the-fly filtering also
1162 : ! avoids overflowing DBCSR's single-precision screening norms for conditioned Z fits.
1163 : CALL dbcsr_multiply('N', 'T', 1.0_dp, tmpB, RB_B, 0.0_dp, B_pan, &
1164 1506 : retain_sparsity=.TRUE.)
1165 1506 : CALL dbcsr_release(tmpB)
1166 1506 : CALL dbcsr_release(LA_panC)
1167 :
1168 3842 : DO i_A = 1, SIZE(M_A)
1169 : ! Reuse A_pan storage for its Hadamard product with the shared B panel.
1170 2336 : CALL hadamard_product_inplace(A_pan(i_A), B_pan, scale)
1171 2336 : CALL create_product_matrix(A_pan(i_A), RB_out, 'N', tmp2)
1172 2336 : CALL dbcsr_multiply('N', 'N', 1.0_dp, A_pan(i_A), RB_out, 0.0_dp, tmp2, filter_eps=eps)
1173 2336 : CALL dbcsr_release(A_pan(i_A))
1174 :
1175 : ! mat_out += L_out_pan^T * tmp2.
1176 2336 : CALL dbcsr_multiply('T', 'N', 1.0_dp, panel_out, tmp2, 1.0_dp, mat_out(i_A)%matrix, filter_eps=eps)
1177 6178 : CALL dbcsr_release(tmp2)
1178 : END DO
1179 1506 : CALL dbcsr_release(B_pan)
1180 1506 : IF (.NOT. lout_eq_la) CALL dbcsr_release(Lout_pan)
1181 1506 : IF (.NOT. lb_eq_la) CALL dbcsr_release(LB_pan)
1182 1506 : CALL dbcsr_release(LA_pan)
1183 4518 : IF (use_cutoff) THEN
1184 42 : CALL dbcsr_release(LA_near)
1185 42 : IF (.NOT. lb_eq_la) CALL dbcsr_release(LB_near)
1186 42 : IF (.NOT. lout_eq_la) CALL dbcsr_release(Lout_near)
1187 : END IF
1188 :
1189 : END DO
1190 :
1191 1506 : CALL timestop(handle)
1192 :
1193 4518 : END SUBROUTINE contract_grid_panels
1194 :
1195 : ! **************************************************************************************************
1196 : !> \brief Build the unfiltered auxiliary interaction M^-1 V^tr M^-1 in the Z-column layout.
1197 : !> AUTO_RI uses its already transformed Coulomb matrix in the optimized auxiliary basis.
1198 : !> Filtering and cache lifetime are owned by the caller.
1199 : !> \param qs_env electronic-structure environment for ordinary RI integrals
1200 : !> \param metric_inverse inverse auxiliary metric
1201 : !> \param trunc_coulomb Coulomb operator for ordinary RI integrals
1202 : !> \param auto_ri optimized auxiliary basis and its Coulomb matrix
1203 : !> \param mat_Z_lP template defining the auxiliary blocking and process grid
1204 : !> \param matrix_V_aux output bare interaction in the auxiliary basis
1205 : ! **************************************************************************************************
1206 40 : SUBROUTINE build_ri_rs_bare_aux_interaction(qs_env, metric_inverse, trunc_coulomb, auto_ri, &
1207 : mat_Z_lP, matrix_V_aux)
1208 : TYPE(qs_environment_type), POINTER :: qs_env
1209 : TYPE(cp_fm_type), INTENT(IN) :: metric_inverse
1210 : TYPE(libint_potential_type), INTENT(IN) :: trunc_coulomb
1211 : TYPE(auto_ri_type), INTENT(IN) :: auto_ri
1212 : TYPE(dbcsr_type), INTENT(IN) :: mat_Z_lP
1213 : TYPE(dbcsr_type), INTENT(OUT) :: matrix_V_aux
1214 :
1215 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_ri_rs_bare_aux_interaction'
1216 :
1217 : INTEGER :: handle
1218 40 : INTEGER, POINTER :: blk_aux(:), dist_row_aux(:)
1219 40 : TYPE(cp_fm_type), ALLOCATABLE :: fm_Vtr_Gamma(:, :)
1220 : TYPE(dbcsr_distribution_type) :: dist_aux_aux
1221 :
1222 40 : CALL timeset(routineN, handle)
1223 40 : IF (auto_ri%enabled) THEN
1224 4 : CPASSERT(auto_ri%ready)
1225 12 : ALLOCATE (fm_Vtr_Gamma(1, 1))
1226 4 : CALL cp_fm_create(fm_Vtr_Gamma(1, 1), metric_inverse%matrix_struct)
1227 4 : CALL cp_fm_to_fm(auto_ri%V_pq, fm_Vtr_Gamma(1, 1))
1228 : ELSE
1229 : CALL RI_2c_integral_mat(qs_env, fm_Vtr_Gamma, metric_inverse%matrix_struct, &
1230 36 : metric_inverse%matrix_struct%nrow_global, trunc_coulomb)
1231 : END IF
1232 40 : CALL fm_contract_ABA(metric_inverse, fm_Vtr_Gamma(:, 1))
1233 40 : CALL setup_square_topology(mat_Z_lP, 'COL', dist_aux_aux, blk_aux, dist_row_aux)
1234 40 : CALL dbcsr_create(matrix_V_aux, "V_aux", dist_aux_aux, dbcsr_type_no_symmetry, blk_aux, blk_aux)
1235 40 : CALL copy_fm_to_dbcsr(fm_Vtr_Gamma(1, 1), matrix_V_aux, keep_sparsity=.FALSE.)
1236 40 : CALL cp_fm_release(fm_Vtr_Gamma)
1237 40 : CALL release_dbcsr_topology_and_matrices(dist=dist_aux_aux, mapped_dist=dist_row_aux)
1238 40 : CALL timestop(handle)
1239 80 : END SUBROUTINE build_ri_rs_bare_aux_interaction
1240 :
1241 : ! **************************************************************************************************
1242 : !> \brief Resolve and cache grid-panel boundaries for RI-RS streaming contractions.
1243 : !> \param ri_rs collocation and fitting matrices, grid geometry and owned panel cache
1244 : !> \param para_env parallel environment
1245 : !> \param msg_elem_limit maximum number of elements per DBCSR message
1246 : !> \param pan_first ...
1247 : !> \param pan_last ...
1248 : !> \param cutoff grid-centroid distance cutoff; zero disables geometric truncation
1249 : ! **************************************************************************************************
1250 174 : SUBROUTINE resolve_grid_panels(ri_rs, para_env, msg_elem_limit, pan_first, pan_last, cutoff)
1251 :
1252 : TYPE(ri_rs_env), INTENT(INOUT) :: ri_rs
1253 : TYPE(mp_para_env_type), POINTER :: para_env
1254 : INTEGER(KIND=int_8), INTENT(IN) :: msg_elem_limit
1255 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: pan_first, pan_last
1256 : REAL(KIND=dp), INTENT(IN) :: cutoff
1257 :
1258 : CHARACTER(LEN=*), PARAMETER :: routineN = 'resolve_grid_panels'
1259 :
1260 : INTEGER :: handle, min_dim, n_ao, n_grid_total, &
1261 : n_RI, npcols, nprows
1262 174 : INTEGER, DIMENSION(:), POINTER :: r_blk_sizes
1263 : LOGICAL :: use_cutoff
1264 : REAL(KIND=dp) :: effective_cutoff, mem_avail_GB, &
1265 : mem_budget_GB
1266 : TYPE(dbcsr_distribution_type) :: dist
1267 :
1268 174 : CALL timeset(routineN, handle)
1269 :
1270 : use_cutoff = cutoff > 0.0_dp .AND. &
1271 174 : ALLOCATED(ri_rs%chunk_centroids)
1272 166 : effective_cutoff = 0.0_dp
1273 8 : IF (use_cutoff) effective_cutoff = cutoff
1274 174 : IF (ALLOCATED(ri_rs%pan_first)) THEN
1275 : ! This is an exact cache key, not a tolerance-based geometric comparison.
1276 134 : IF (ri_rs%panel_plan_cutoff == effective_cutoff) THEN
1277 536 : ALLOCATE (pan_first, SOURCE=ri_rs%pan_first)
1278 536 : ALLOCATE (pan_last, SOURCE=ri_rs%pan_last)
1279 134 : CALL timestop(handle)
1280 : RETURN
1281 : END IF
1282 0 : DEALLOCATE (ri_rs%pan_first, ri_rs%pan_last)
1283 : END IF
1284 :
1285 : ! MIN(nprows, npcols) is the divisor that bounds the worst-rank Cannon message: a
1286 : ! P x n_grid panel is replicated into block row strips (P/nprows x n_grid) or column
1287 : ! strips (P x n_grid/npcols) during multiply_cannon, so the largest single-rank
1288 : ! message is ~ P*n_grid / MIN(nprows,npcols) elements.
1289 : CALL dbcsr_get_info(ri_rs%mat_phi_mu_l, nfullrows_total=n_grid_total, nfullcols_total=n_ao, &
1290 40 : row_blk_size=r_blk_sizes, distribution=dist)
1291 40 : CALL dbcsr_distribution_get(dist, nprows=nprows, npcols=npcols)
1292 40 : min_dim = MAX(MIN(nprows, npcols), 1)
1293 :
1294 40 : CALL dbcsr_get_info(ri_rs%mat_Z_lP, nfullcols_total=n_RI)
1295 : ! Only the cutoff path uses the measured free memory to limit panel sizes.
1296 40 : mem_budget_GB = 0.0_dp
1297 40 : IF (use_cutoff) THEN
1298 2 : CALL mp_mem_avail_per_rank_GB(para_env, mem_avail_GB)
1299 2 : mem_budget_GB = 0.5_dp*mem_avail_GB
1300 : END IF
1301 : CALL plan_grid_panels(msg_elem_limit, r_blk_sizes, min_dim, pan_first, pan_last, &
1302 : n_RI, n_ao, para_env%num_pe, mem_budget_GB, &
1303 40 : centroids=ri_rs%chunk_centroids, cutoff=effective_cutoff)
1304 :
1305 160 : ALLOCATE (ri_rs%pan_first, SOURCE=pan_first)
1306 160 : ALLOCATE (ri_rs%pan_last, SOURCE=pan_last)
1307 40 : ri_rs%panel_plan_cutoff = effective_cutoff
1308 :
1309 40 : CALL timestop(handle)
1310 :
1311 214 : END SUBROUTINE resolve_grid_panels
1312 :
1313 : END MODULE ri_rs_matrices
|
