(git:374b731)
Loading...
Searching...
No Matches
mp2_ri_grad_util.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines for calculating RI-MP2 gradients
10!> \par History
11!> 10.2013 created [Mauro Del Ben]
12! **************************************************************************************************
13MODULE mp2_ri_grad_util
14!
15 USE cp_blacs_env, ONLY: cp_blacs_env_create,&
16 cp_blacs_env_release,&
17 cp_blacs_env_type
18 USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,&
19 cp_dbcsr_m_by_n_from_template
20 USE cp_fm_basic_linalg, ONLY: cp_fm_geadd
21 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
22 cp_fm_struct_get,&
23 cp_fm_struct_release,&
24 cp_fm_struct_type
25 USE cp_fm_types, ONLY: &
26 cp_fm_create, cp_fm_get_info, cp_fm_indxg2l, cp_fm_indxg2p, cp_fm_indxl2g, cp_fm_release, &
27 cp_fm_set_all, cp_fm_to_fm, cp_fm_type
28 USE dbcsr_api, ONLY: dbcsr_p_type,&
29 dbcsr_type,&
30 dbcsr_type_no_symmetry
31 USE group_dist_types, ONLY: create_group_dist,&
32 get_group_dist,&
33 group_dist_d1_type,&
34 release_group_dist
35 USE kinds, ONLY: dp
36 USE libint_2c_3c, ONLY: compare_potential_types
37 USE message_passing, ONLY: mp_para_env_type,&
38 mp_request_null,&
39 mp_request_type,&
40 mp_waitall
41 USE mp2_types, ONLY: integ_mat_buffer_type,&
42 mp2_type
43 USE parallel_gemm_api, ONLY: parallel_gemm
44 USE util, ONLY: get_limit
45#include "./base/base_uses.f90"
46
47 IMPLICIT NONE
48
49 PRIVATE
50
51 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_ri_grad_util'
52
53 PUBLIC :: complete_gamma, array2fm, fm2array, prepare_redistribution, create_dbcsr_gamma
54
55 TYPE index_map
56 INTEGER, DIMENSION(:, :), ALLOCATABLE :: map
57 END TYPE
58
59CONTAINS
60
61! **************************************************************************************************
62!> \brief complete the calculation of the Gamma matrices
63!> \param mp2_env ...
64!> \param B_ia_Q ...
65!> \param dimen_RI ...
66!> \param homo ...
67!> \param virtual ...
68!> \param para_env ...
69!> \param para_env_sub ...
70!> \param ngroup ...
71!> \param my_group_L_size ...
72!> \param my_group_L_start ...
73!> \param my_group_L_end ...
74!> \param my_B_size ...
75!> \param my_B_virtual_start ...
76!> \param gd_array ...
77!> \param gd_B_virtual ...
78!> \param kspin ...
79!> \author Mauro Del Ben
80! **************************************************************************************************
81 SUBROUTINE complete_gamma(mp2_env, B_ia_Q, dimen_RI, homo, virtual, para_env, para_env_sub, ngroup, &
82 my_group_L_size, my_group_L_start, my_group_L_end, &
83 my_B_size, my_B_virtual_start, gd_array, gd_B_virtual, kspin)
84
85 TYPE(mp2_type) :: mp2_env
86 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
87 INTENT(INOUT) :: b_ia_q
88 INTEGER, INTENT(IN) :: dimen_ri, homo, virtual
89 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env, para_env_sub
90 INTEGER, INTENT(IN) :: ngroup, my_group_l_size, &
91 my_group_l_start, my_group_l_end, &
92 my_b_size, my_b_virtual_start
93 TYPE(group_dist_d1_type), INTENT(IN) :: gd_array, gd_b_virtual
94 INTEGER, INTENT(IN) :: kspin
95
96 CHARACTER(LEN=*), PARAMETER :: routinen = 'complete_gamma'
97
98 INTEGER :: dimen_ia, handle, i, ispin, kkb, my_ia_end, my_ia_size, my_ia_start, my_p_end, &
99 my_p_size, my_p_start, nspins, pos_group, pos_sub
100 INTEGER, ALLOCATABLE, DIMENSION(:) :: pos_info
101 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: group_grid_2_mepos, mepos_2_grid_group
102 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: bib_c_2d, gamma_2d, gamma_pq
103 TYPE(cp_blacs_env_type), POINTER :: blacs_env
104 TYPE(cp_fm_struct_type), POINTER :: fm_struct_ia, fm_struct_ri
105 TYPE(cp_fm_type) :: fm_gamma, fm_gamma_pq, fm_gamma_pq_2, fm_gamma_pq_temp, &
106 fm_gamma_pq_temp_2, fm_ia_p, fm_y, operator_half, pq_half
107 TYPE(group_dist_d1_type) :: gd_array_new, gd_ia, gd_ia_new, gd_p, &
108 gd_p_new
109
110 CALL timeset(routinen, handle)
111
112 ! reshape the data into a 2D array
113 ! reorder the local data in such a way to help the next stage of matrix creation
114 ! now the data inside the group are divided into a ia x K matrix
115 dimen_ia = homo*virtual
116 CALL create_group_dist(gd_ia, para_env_sub%num_pe, dimen_ia)
117 CALL get_group_dist(gd_ia, para_env_sub%mepos, my_ia_start, my_ia_end, my_ia_size)
118
119 CALL mat_3d_to_2d(b_ia_q, bib_c_2d, homo, my_b_size, virtual, my_b_virtual_start, &
120 my_ia_start, my_ia_end, my_ia_size, my_group_l_size, para_env_sub, gd_b_virtual)
121
122 CALL mat_3d_to_2d_gamma(mp2_env%ri_grad%Gamma_P_ia(kspin)%array, gamma_2d, homo, my_b_size, virtual, my_b_virtual_start, &
123 my_ia_start, my_ia_end, my_ia_size, my_group_l_size, para_env_sub, gd_b_virtual)
124
125 ! create the processor map and size arrays
126 CALL create_group_dist(gd_ia_new, para_env%num_pe)
127 CALL create_group_dist(gd_array_new, para_env%num_pe)
128 CALL create_group_dist(gd_p, para_env_sub%num_pe, dimen_ri)
129 CALL create_group_dist(gd_p_new, para_env%num_pe)
130
131 CALL get_group_dist(gd_p, para_env_sub%mepos, my_p_start, my_p_end, my_p_size)
132
133 CALL prepare_redistribution(para_env, para_env_sub, ngroup, &
134 group_grid_2_mepos, mepos_2_grid_group, pos_info=pos_info)
135
136 DO i = 0, para_env%num_pe - 1
137 ! calculate position of the group
138 pos_group = i/para_env_sub%num_pe
139 ! calculate position in the subgroup
140 pos_sub = pos_info(i)
141 ! 1 -> rows, 2 -> cols
142 CALL get_group_dist(gd_ia, pos_sub, gd_ia_new, i)
143 CALL get_group_dist(gd_array, pos_group, gd_array_new, i)
144 CALL get_group_dist(gd_p, pos_sub, gd_p_new, i)
145 END DO
146
147 ! create the blacs env
148 NULLIFY (blacs_env)
149 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
150
151 NULLIFY (fm_struct_ia)
152 CALL cp_fm_struct_create(fm_struct_ia, context=blacs_env, nrow_global=dimen_ia, &
153 ncol_global=dimen_ri, para_env=para_env)
154
155 ! create the fm matrix Gamma
156 CALL array2fm(gamma_2d, fm_struct_ia, my_ia_start, my_ia_end, &
157 my_group_l_start, my_group_l_end, &
158 gd_ia_new, gd_array_new, &
159 group_grid_2_mepos, para_env_sub%num_pe, ngroup, &
160 fm_y)
161 ! create the fm matrix B_ia_P
162 CALL array2fm(bib_c_2d, fm_struct_ia, my_ia_start, my_ia_end, &
163 my_group_l_start, my_group_l_end, &
164 gd_ia_new, gd_array_new, &
165 group_grid_2_mepos, para_env_sub%num_pe, ngroup, &
166 fm_ia_p)
167
168 NULLIFY (fm_struct_ri)
169 CALL cp_fm_struct_create(fm_struct_ri, context=blacs_env, nrow_global=dimen_ri, &
170 ncol_global=dimen_ri, para_env=para_env)
171
172 ! since we will need (P|Q)^(-1/2) in the future, make a copy
173 CALL array2fm(mp2_env%ri_grad%PQ_half, fm_struct_ri, my_p_start, my_p_end, &
174 my_group_l_start, my_group_l_end, &
175 gd_p_new, gd_array_new, &
176 group_grid_2_mepos, para_env_sub%num_pe, ngroup, &
177 pq_half, do_release_mat=.false.)
178
179 IF (.NOT. compare_potential_types(mp2_env%ri_metric, mp2_env%potential_parameter)) THEN
180 CALL array2fm(mp2_env%ri_grad%operator_half, fm_struct_ri, my_p_start, my_p_end, &
181 my_group_l_start, my_group_l_end, &
182 gd_p_new, gd_array_new, &
183 group_grid_2_mepos, para_env_sub%num_pe, ngroup, &
184 operator_half, do_release_mat=.false.)
185 END IF
186
187 CALL release_group_dist(gd_p_new)
188 CALL release_group_dist(gd_ia_new)
189 CALL release_group_dist(gd_array_new)
190
191 ! complete the gamma matrix Gamma_ia^P = sum_Q (Y_ia^Q * (Q|P)^(-1/2))
192 IF (compare_potential_types(mp2_env%ri_metric, mp2_env%potential_parameter)) THEN
193 CALL cp_fm_create(fm_gamma, fm_struct_ia, name="fm_Gamma")
194 CALL cp_fm_struct_release(fm_struct_ia)
195 ! perform the matrix multiplication
196 CALL parallel_gemm(transa="N", transb="T", m=dimen_ia, n=dimen_ri, k=dimen_ri, alpha=1.0_dp, &
197 matrix_a=fm_y, matrix_b=pq_half, beta=0.0_dp, &
198 matrix_c=fm_gamma)
199 ! release the Y matrix
200 CALL cp_fm_release(fm_y)
201
202 ! complete gamma small (fm_Gamma_PQ)
203 ! create temp matrix
204 CALL cp_fm_create(fm_gamma_pq_temp, fm_struct_ri, name="fm_Gamma_PQ_temp")
205 CALL parallel_gemm(transa="T", transb="N", m=dimen_ri, n=dimen_ri, k=dimen_ia, alpha=1.0_dp, &
206 matrix_a=fm_gamma, matrix_b=fm_ia_p, beta=0.0_dp, &
207 matrix_c=fm_gamma_pq_temp)
208 CALL cp_fm_release(fm_ia_p)
209 ! create fm_Gamma_PQ matrix
210 CALL cp_fm_create(fm_gamma_pq, fm_struct_ri, name="fm_Gamma_PQ")
211 CALL cp_fm_struct_release(fm_struct_ri)
212 ! perform matrix multiplication
213 CALL parallel_gemm(transa="N", transb="T", m=dimen_ri, n=dimen_ri, k=dimen_ri, alpha=1.0_dp, &
214 matrix_a=fm_gamma_pq_temp, matrix_b=pq_half, beta=0.0_dp, &
215 matrix_c=fm_gamma_pq)
216 CALL cp_fm_release(fm_gamma_pq_temp)
217 CALL cp_fm_release(pq_half)
218
219 ELSE
220 ! create temp matrix
221 CALL cp_fm_create(fm_gamma_pq_temp, fm_struct_ri, name="fm_Gamma_PQ_temp")
222 CALL parallel_gemm(transa="T", transb="N", m=dimen_ri, n=dimen_ri, k=dimen_ia, alpha=1.0_dp, &
223 matrix_a=fm_y, matrix_b=fm_ia_p, beta=0.0_dp, &
224 matrix_c=fm_gamma_pq_temp)
225 CALL cp_fm_release(fm_ia_p)
226
227 CALL cp_fm_create(fm_gamma, fm_struct_ia, name="fm_Gamma")
228 CALL cp_fm_struct_release(fm_struct_ia)
229 ! perform the matrix multiplication
230 CALL parallel_gemm(transa="N", transb="T", m=dimen_ia, n=dimen_ri, k=dimen_ri, alpha=1.0_dp, &
231 matrix_a=fm_y, matrix_b=pq_half, beta=0.0_dp, &
232 matrix_c=fm_gamma)
233 CALL cp_fm_release(fm_y)
234
235 CALL cp_fm_create(fm_gamma_pq_temp_2, fm_struct_ri, name="fm_Gamma_PQ_temp_2")
236 CALL parallel_gemm(transa="N", transb="T", m=dimen_ri, n=dimen_ri, k=dimen_ri, alpha=1.0_dp, &
237 matrix_a=fm_gamma_pq_temp, matrix_b=operator_half, beta=0.0_dp, &
238 matrix_c=fm_gamma_pq_temp_2)
239
240 CALL cp_fm_create(fm_gamma_pq_2, fm_struct_ri, name="fm_Gamma_PQ_2")
241 CALL parallel_gemm(transa="N", transb="N", m=dimen_ri, n=dimen_ri, k=dimen_ri, alpha=1.0_dp, &
242 matrix_a=pq_half, matrix_b=fm_gamma_pq_temp_2, beta=0.0_dp, &
243 matrix_c=fm_gamma_pq_temp)
244 CALL cp_fm_to_fm(fm_gamma_pq_temp, fm_gamma_pq_2)
245 CALL cp_fm_geadd(1.0_dp, "T", fm_gamma_pq_temp, 1.0_dp, fm_gamma_pq_2)
246 CALL cp_fm_release(fm_gamma_pq_temp)
247 CALL cp_fm_release(pq_half)
248
249 CALL cp_fm_create(fm_gamma_pq, fm_struct_ri)
250 CALL cp_fm_struct_release(fm_struct_ri)
251 CALL parallel_gemm(transa="N", transb="N", m=dimen_ri, n=dimen_ri, k=dimen_ri, alpha=-1.0_dp, &
252 matrix_a=operator_half, matrix_b=fm_gamma_pq_temp_2, beta=0.0_dp, &
253 matrix_c=fm_gamma_pq)
254 CALL cp_fm_release(operator_half)
255 CALL cp_fm_release(fm_gamma_pq_temp_2)
256 END IF
257
258 ! now redistribute the data, in this case we go back
259 ! to the original way the integrals were distributed
260 CALL fm2array(gamma_2d, my_ia_size, my_ia_start, my_ia_end, &
261 my_group_l_size, my_group_l_start, my_group_l_end, &
262 group_grid_2_mepos, mepos_2_grid_group, &
263 para_env_sub%num_pe, ngroup, &
264 fm_gamma)
265
266 ALLOCATE (gamma_pq(my_p_size, my_group_l_size))
267 CALL fm2array(gamma_pq, my_p_size, my_p_start, my_p_end, &
268 my_group_l_size, my_group_l_start, my_group_l_end, &
269 group_grid_2_mepos, mepos_2_grid_group, &
270 para_env_sub%num_pe, ngroup, &
271 fm_gamma_pq)
272 IF (.NOT. ALLOCATED(mp2_env%ri_grad%Gamma_PQ)) THEN
273 CALL move_alloc(gamma_pq, mp2_env%ri_grad%Gamma_PQ)
274 ELSE
275 mp2_env%ri_grad%Gamma_PQ(:, :) = mp2_env%ri_grad%Gamma_PQ + gamma_pq
276 DEALLOCATE (gamma_pq)
277 END IF
278
279 IF (.NOT. compare_potential_types(mp2_env%ri_metric, mp2_env%potential_parameter)) THEN
280 ALLOCATE (gamma_pq(my_p_size, my_group_l_size))
281 CALL fm2array(gamma_pq, my_p_size, my_p_start, my_p_end, &
282 my_group_l_size, my_group_l_start, my_group_l_end, &
283 group_grid_2_mepos, mepos_2_grid_group, &
284 para_env_sub%num_pe, ngroup, &
285 fm_gamma_pq_2)
286 IF (.NOT. ALLOCATED(mp2_env%ri_grad%Gamma_PQ_2)) THEN
287 CALL move_alloc(gamma_pq, mp2_env%ri_grad%Gamma_PQ_2)
288 ELSE
289 mp2_env%ri_grad%Gamma_PQ_2(:, :) = mp2_env%ri_grad%Gamma_PQ_2 + gamma_pq
290 DEALLOCATE (gamma_pq)
291 END IF
292 END IF
293
294 ! allocate G_P_ia (DBCSR)
295 IF (.NOT. ALLOCATED(mp2_env%ri_grad%G_P_ia)) THEN
296 nspins = SIZE(mp2_env%ri_grad%mo_coeff_o)
297 ALLOCATE (mp2_env%ri_grad%G_P_ia(my_group_l_size, nspins))
298 DO ispin = 1, nspins
299 DO kkb = 1, my_group_l_size
300 NULLIFY (mp2_env%ri_grad%G_P_ia(kkb, ispin)%matrix)
301 END DO
302 END DO
303 END IF
304
305 ! create the Gamma_ia_P in DBCSR style
306 CALL create_dbcsr_gamma(gamma_2d, homo, virtual, dimen_ia, para_env_sub, &
307 my_ia_start, my_ia_end, my_group_l_size, gd_ia, &
308 mp2_env%ri_grad%G_P_ia(:, kspin), mp2_env%ri_grad%mo_coeff_o(1)%matrix)
309
310 DEALLOCATE (pos_info)
311 DEALLOCATE (group_grid_2_mepos, mepos_2_grid_group)
312 CALL release_group_dist(gd_ia)
313 CALL release_group_dist(gd_p)
314
315 ! release blacs_env
316 CALL cp_blacs_env_release(blacs_env)
317
318 CALL timestop(handle)
319
320 END SUBROUTINE complete_gamma
321
322! **************************************************************************************************
323!> \brief Redistribute a 3d matrix to a 2d matrix
324!> \param B_ia_Q ...
325!> \param BIb_C_2D ...
326!> \param homo ...
327!> \param my_B_size ...
328!> \param virtual ...
329!> \param my_B_virtual_start ...
330!> \param my_ia_start ...
331!> \param my_ia_end ...
332!> \param my_ia_size ...
333!> \param my_group_L_size ...
334!> \param para_env_sub ...
335!> \param gd_B_virtual ...
336! **************************************************************************************************
337 SUBROUTINE mat_3d_to_2d(B_ia_Q, BIb_C_2D, homo, my_B_size, virtual, my_B_virtual_start, &
338 my_ia_start, my_ia_end, my_ia_size, my_group_L_size, para_env_sub, gd_B_virtual)
339 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
340 INTENT(INOUT) :: b_ia_q
341 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
342 INTENT(OUT) :: bib_c_2d
343 INTEGER, INTENT(IN) :: homo, my_b_size, virtual, &
344 my_b_virtual_start, my_ia_start, &
345 my_ia_end, my_ia_size, my_group_l_size
346 TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
347 TYPE(group_dist_d1_type), INTENT(IN) :: gd_b_virtual
348
349 CHARACTER(LEN=*), PARAMETER :: routinen = 'mat_3d_to_2d'
350
351 INTEGER :: handle, ia_global, iib, jjb, proc_receive, proc_send, proc_shift, rec_b_size, &
352 rec_b_virtual_end, rec_b_virtual_start
353 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: bib_c_rec
354
355 CALL timeset(routinen, handle)
356
357 ALLOCATE (bib_c_2d(my_ia_size, my_group_l_size))
358 bib_c_2d = 0.0_dp
359
360 DO iib = 1, homo
361 DO jjb = 1, my_b_size
362 ia_global = (iib - 1)*virtual + my_b_virtual_start + jjb - 1
363 IF (ia_global >= my_ia_start .AND. ia_global <= my_ia_end) THEN
364 bib_c_2d(ia_global - my_ia_start + 1, :) = b_ia_q(iib, jjb, 1:my_group_l_size)
365 END IF
366 END DO
367 END DO
368
369 DO proc_shift = 1, para_env_sub%num_pe - 1
370 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
371 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
372
373 CALL get_group_dist(gd_b_virtual, proc_receive, rec_b_virtual_start, rec_b_virtual_end, rec_b_size)
374
375 ALLOCATE (bib_c_rec(homo, rec_b_size, my_group_l_size))
376 bib_c_rec = 0.0_dp
377
378 CALL para_env_sub%sendrecv(b_ia_q, proc_send, &
379 bib_c_rec, proc_receive)
380
381 DO iib = 1, homo
382 DO jjb = 1, rec_b_size
383 ia_global = (iib - 1)*virtual + rec_b_virtual_start + jjb - 1
384 IF (ia_global >= my_ia_start .AND. ia_global <= my_ia_end) THEN
385 bib_c_2d(ia_global - my_ia_start + 1, :) = bib_c_rec(iib, jjb, 1:my_group_l_size)
386 END IF
387 END DO
388 END DO
389
390 DEALLOCATE (bib_c_rec)
391 END DO
392 DEALLOCATE (b_ia_q)
393
394 CALL timestop(handle)
395 END SUBROUTINE mat_3d_to_2d
396
397! **************************************************************************************************
398!> \brief Redistribute a 3d matrix to a 2d matrix, specialized for Gamma_P_ia to account for a different order of indices
399!> \param B_ia_Q ...
400!> \param BIb_C_2D ...
401!> \param homo ...
402!> \param my_B_size ...
403!> \param virtual ...
404!> \param my_B_virtual_start ...
405!> \param my_ia_start ...
406!> \param my_ia_end ...
407!> \param my_ia_size ...
408!> \param my_group_L_size ...
409!> \param para_env_sub ...
410!> \param gd_B_virtual ...
411! **************************************************************************************************
412 SUBROUTINE mat_3d_to_2d_gamma(B_ia_Q, BIb_C_2D, homo, my_B_size, virtual, my_B_virtual_start, &
413 my_ia_start, my_ia_end, my_ia_size, my_group_L_size, para_env_sub, gd_B_virtual)
414 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
415 INTENT(INOUT) :: b_ia_q
416 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
417 INTENT(OUT) :: bib_c_2d
418 INTEGER, INTENT(IN) :: homo, my_b_size, virtual, &
419 my_b_virtual_start, my_ia_start, &
420 my_ia_end, my_ia_size, my_group_l_size
421 TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
422 TYPE(group_dist_d1_type), INTENT(IN) :: gd_b_virtual
423
424 CHARACTER(LEN=*), PARAMETER :: routinen = 'mat_3d_to_2d_gamma'
425
426 INTEGER :: handle, ia_global, iib, jjb, proc_receive, proc_send, proc_shift, rec_b_size, &
427 rec_b_virtual_end, rec_b_virtual_start
428 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: bib_c_rec
429
430 CALL timeset(routinen, handle)
431
432 ALLOCATE (bib_c_2d(my_ia_size, my_group_l_size))
433 bib_c_2d = 0.0_dp
434
435 DO iib = 1, homo
436 DO jjb = 1, my_b_size
437 ia_global = (iib - 1)*virtual + my_b_virtual_start + jjb - 1
438 IF (ia_global >= my_ia_start .AND. ia_global <= my_ia_end) THEN
439 bib_c_2d(ia_global - my_ia_start + 1, :) = b_ia_q(jjb, iib, 1:my_group_l_size)
440 END IF
441 END DO
442 END DO
443
444 DO proc_shift = 1, para_env_sub%num_pe - 1
445 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
446 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
447
448 CALL get_group_dist(gd_b_virtual, proc_receive, rec_b_virtual_start, rec_b_virtual_end, rec_b_size)
449
450 ALLOCATE (bib_c_rec(rec_b_size, homo, my_group_l_size))
451 bib_c_rec = 0.0_dp
452
453 CALL para_env_sub%sendrecv(b_ia_q, proc_send, &
454 bib_c_rec, proc_receive)
455
456 DO iib = 1, homo
457 DO jjb = 1, rec_b_size
458 ia_global = (iib - 1)*virtual + rec_b_virtual_start + jjb - 1
459 IF (ia_global >= my_ia_start .AND. ia_global <= my_ia_end) THEN
460 bib_c_2d(ia_global - my_ia_start + 1, :) = bib_c_rec(jjb, iib, 1:my_group_l_size)
461 END IF
462 END DO
463 END DO
464
465 DEALLOCATE (bib_c_rec)
466 END DO
467 DEALLOCATE (b_ia_q)
468
469 CALL timestop(handle)
470 END SUBROUTINE mat_3d_to_2d_gamma
471
472! **************************************************************************************************
473!> \brief redistribute local part of array to fm
474!> \param mat2D ...
475!> \param fm_struct ...
476!> \param my_start_row ...
477!> \param my_end_row ...
478!> \param my_start_col ...
479!> \param my_end_col ...
480!> \param gd_row ...
481!> \param gd_col ...
482!> \param group_grid_2_mepos ...
483!> \param ngroup_row ...
484!> \param ngroup_col ...
485!> \param fm_mat ...
486!> \param integ_group_size ...
487!> \param color_group ...
488!> \param do_release_mat whether to release the array (default: yes)
489! **************************************************************************************************
490 SUBROUTINE array2fm(mat2D, fm_struct, my_start_row, my_end_row, &
491 my_start_col, my_end_col, &
492 gd_row, gd_col, &
493 group_grid_2_mepos, ngroup_row, ngroup_col, &
494 fm_mat, integ_group_size, color_group, do_release_mat)
495
496 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
497 INTENT(INOUT) :: mat2d
498 TYPE(cp_fm_struct_type), POINTER :: fm_struct
499 INTEGER, INTENT(IN) :: my_start_row, my_end_row, my_start_col, &
500 my_end_col
501 TYPE(group_dist_d1_type), INTENT(IN) :: gd_row, gd_col
502 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN) :: group_grid_2_mepos
503 INTEGER, INTENT(IN) :: ngroup_row, ngroup_col
504 TYPE(cp_fm_type), INTENT(OUT) :: fm_mat
505 INTEGER, INTENT(IN), OPTIONAL :: integ_group_size, color_group
506 LOGICAL, INTENT(IN), OPTIONAL :: do_release_mat
507
508 CHARACTER(LEN=*), PARAMETER :: routinen = 'array2fm'
509
510 INTEGER :: dummy_proc, end_col_block, end_row_block, handle, handle2, i_global, i_local, &
511 i_sub, iib, iii, itmp(2), j_global, j_local, j_sub, jjb, my_num_col_blocks, &
512 my_num_row_blocks, mypcol, myprow, ncol_block, ncol_local, npcol, nprow, nrow_block, &
513 nrow_local, num_cols, num_rec_cols, num_rows, number_of_rec, number_of_send, &
514 proc_receive, proc_send, proc_shift, rec_col_end, rec_col_size, rec_col_start, &
515 rec_counter, rec_pcol, rec_prow, rec_row_end, rec_row_size, rec_row_start, ref_send_pcol, &
516 ref_send_prow, send_counter, send_pcol, send_prow, size_rec_buffer, size_send_buffer
517 INTEGER :: start_col_block, start_row_block
518 INTEGER, ALLOCATABLE, DIMENSION(:) :: iii_vet, index_col_rec, map_rec_size, &
519 map_send_size
520 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: blocks_ranges_col, blocks_ranges_row, &
521 grid_2_mepos, grid_ref_2_send_pos, &
522 mepos_2_grid
523 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
524 LOGICAL :: convert_pos, my_do_release_mat
525 REAL(kind=dp) :: part_col, part_row
526 TYPE(integ_mat_buffer_type), ALLOCATABLE, &
527 DIMENSION(:) :: buffer_rec, buffer_send
528 TYPE(mp_para_env_type), POINTER :: para_env
529 TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: req_send
530
531 CALL timeset(routinen, handle)
532
533 my_do_release_mat = .true.
534 IF (PRESENT(do_release_mat)) my_do_release_mat = do_release_mat
535
536 CALL cp_fm_struct_get(fm_struct, para_env=para_env, nrow_global=num_rows, ncol_global=num_cols)
537
538 ! create the full matrix, (num_rows x num_cols)
539 CALL cp_fm_create(fm_mat, fm_struct, name="fm_mat")
540 CALL cp_fm_set_all(matrix=fm_mat, alpha=0.0_dp)
541
542 ! start filling procedure
543 ! fill the matrix
544 CALL cp_fm_get_info(matrix=fm_mat, &
545 nrow_local=nrow_local, &
546 ncol_local=ncol_local, &
547 row_indices=row_indices, &
548 col_indices=col_indices, &
549 nrow_block=nrow_block, &
550 ncol_block=ncol_block)
551 myprow = fm_mat%matrix_struct%context%mepos(1)
552 mypcol = fm_mat%matrix_struct%context%mepos(2)
553 nprow = fm_mat%matrix_struct%context%num_pe(1)
554 npcol = fm_mat%matrix_struct%context%num_pe(2)
555
556 ! start the communication part
557 ! 0) create array containing the processes position
558 ! and supporting infos
559 CALL timeset(routinen//"_info", handle2)
560 ALLOCATE (grid_2_mepos(0:nprow - 1, 0:npcol - 1))
561 grid_2_mepos = 0
562 ALLOCATE (mepos_2_grid(2, 0:para_env%num_pe - 1))
563 grid_2_mepos(myprow, mypcol) = para_env%mepos
564 ! sum infos
565 CALL para_env%sum(grid_2_mepos)
566 CALL para_env%allgather([myprow, mypcol], mepos_2_grid)
567
568 ! 1) loop over my local data and define a map for the proc to send data
569 ALLOCATE (map_send_size(0:para_env%num_pe - 1))
570 map_send_size = 0
571 dummy_proc = 0
572 DO jjb = my_start_col, my_end_col
573 send_pcol = cp_fm_indxg2p(jjb, ncol_block, dummy_proc, &
574 fm_mat%matrix_struct%first_p_pos(2), npcol)
575 DO iib = my_start_row, my_end_row
576 send_prow = cp_fm_indxg2p(iib, nrow_block, dummy_proc, &
577 fm_mat%matrix_struct%first_p_pos(1), nprow)
578 proc_send = grid_2_mepos(send_prow, send_pcol)
579 map_send_size(proc_send) = map_send_size(proc_send) + 1
580 END DO
581 END DO
582
583 ! 2) loop over my local data of fm_mat and define a map for the proc from which rec data
584 ALLOCATE (map_rec_size(0:para_env%num_pe - 1))
585 map_rec_size = 0
586 part_row = real(num_rows, kind=dp)/real(ngroup_row, kind=dp)
587 part_col = real(num_cols, kind=dp)/real(ngroup_col, kind=dp)
588
589 ! In some cases we have to convert global positions to positions in a subgroup (RPA)
590 convert_pos = .false.
591 IF (PRESENT(integ_group_size) .AND. PRESENT(color_group)) convert_pos = .true.
592
593 DO jjb = 1, nrow_local
594 j_global = row_indices(jjb)
595 ! check the group holding this element
596 ! dirty way, if someone has a better idea ...
597 rec_prow = int(real(j_global - 1, kind=dp)/part_row)
598 rec_prow = max(0, rec_prow)
599 rec_prow = min(rec_prow, ngroup_row - 1)
600 DO
601 itmp = get_limit(num_rows, ngroup_row, rec_prow)
602 IF (j_global >= itmp(1) .AND. j_global <= itmp(2)) EXIT
603 IF (j_global < itmp(1)) rec_prow = rec_prow - 1
604 IF (j_global > itmp(2)) rec_prow = rec_prow + 1
605 END DO
606
607 IF (convert_pos) THEN
608 ! if the group is not in the same RPA group cycle
609 IF ((rec_prow/integ_group_size) .NE. color_group) cycle
610 ! convert global position to position into the RPA group
611 rec_prow = mod(rec_prow, integ_group_size)
612 END IF
613
614 DO iib = 1, ncol_local
615 i_global = col_indices(iib)
616 ! check the process in the group holding this element
617 rec_pcol = int(real(i_global - 1, kind=dp)/part_col)
618 rec_pcol = max(0, rec_pcol)
619 rec_pcol = min(rec_pcol, ngroup_col - 1)
620 DO
621 itmp = get_limit(num_cols, ngroup_col, rec_pcol)
622 IF (i_global >= itmp(1) .AND. i_global <= itmp(2)) EXIT
623 IF (i_global < itmp(1)) rec_pcol = rec_pcol - 1
624 IF (i_global > itmp(2)) rec_pcol = rec_pcol + 1
625 END DO
626
627 proc_receive = group_grid_2_mepos(rec_prow, rec_pcol)
628
629 map_rec_size(proc_receive) = map_rec_size(proc_receive) + 1
630
631 END DO ! i_global
632 END DO ! j_global
633
634 ! 3) check if the local data has to be stored in the new fm_mat
635 IF (map_rec_size(para_env%mepos) > 0) THEN
636 DO jjb = 1, ncol_local
637 j_global = col_indices(jjb)
638 IF (j_global >= my_start_col .AND. j_global <= my_end_col) THEN
639 DO iib = 1, nrow_local
640 i_global = row_indices(iib)
641 IF (i_global >= my_start_row .AND. i_global <= my_end_row) THEN
642 fm_mat%local_data(iib, jjb) = mat2d(i_global - my_start_row + 1, j_global - my_start_col + 1)
643 END IF
644 END DO
645 END IF
646 END DO
647 END IF
648 CALL timestop(handle2)
649
650 ! 4) reorder data in the send_buffer
651 CALL timeset(routinen//"_buffer_send", handle2)
652 ! count the number of messages to send
653 number_of_send = 0
654 DO proc_shift = 1, para_env%num_pe - 1
655 proc_send = modulo(para_env%mepos + proc_shift, para_env%num_pe)
656 IF (map_send_size(proc_send) > 0) THEN
657 number_of_send = number_of_send + 1
658 END IF
659 END DO
660 ! allocate the structure that will hold the messages to be sent
661 ALLOCATE (buffer_send(number_of_send))
662 ! grid_ref_2_send_pos is an array (map) that given a pair
663 ! (ref_send_prow,ref_send_pcol) returns
664 ! the position in the buffer_send associated to that process
665 ALLOCATE (grid_ref_2_send_pos(0:nprow - 1, 0:npcol - 1))
666 grid_ref_2_send_pos = 0
667 ! finalize the allocation of buffer_send with the actual size
668 ! of each message (actual size is size_send_buffer)
669 send_counter = 0
670 DO proc_shift = 1, para_env%num_pe - 1
671 proc_send = modulo(para_env%mepos + proc_shift, para_env%num_pe)
672 size_send_buffer = map_send_size(proc_send)
673 IF (map_send_size(proc_send) > 0) THEN
674 send_counter = send_counter + 1
675 ! allocate the sending buffer (msg)
676 ALLOCATE (buffer_send(send_counter)%msg(size_send_buffer))
677 buffer_send(send_counter)%msg = 0.0_dp
678 buffer_send(send_counter)%proc = proc_send
679 ! get the pointer to prow, pcol of the process that has
680 ! to receive this message
681 ref_send_prow = mepos_2_grid(1, proc_send)
682 ref_send_pcol = mepos_2_grid(2, proc_send)
683 ! save the rank of the process that has to receive this message
684 grid_ref_2_send_pos(ref_send_prow, ref_send_pcol) = send_counter
685 END IF
686 END DO
687
688 ! loop over the locally held data and fill the buffer_send
689 ! for doing that we need an array that keep track if the
690 ! sequential increase of the index for each message
691 ALLOCATE (iii_vet(number_of_send))
692 iii_vet = 0
693 DO iib = my_start_row, my_end_row
694 send_prow = cp_fm_indxg2p(iib, nrow_block, dummy_proc, &
695 fm_mat%matrix_struct%first_p_pos(1), nprow)
696 DO jjb = my_start_col, my_end_col
697 send_pcol = cp_fm_indxg2p(jjb, ncol_block, dummy_proc, &
698 fm_mat%matrix_struct%first_p_pos(2), npcol)
699 ! we don't need to send to ourselves
700 IF (grid_2_mepos(send_prow, send_pcol) == para_env%mepos) cycle
701
702 send_counter = grid_ref_2_send_pos(send_prow, send_pcol)
703 iii_vet(send_counter) = iii_vet(send_counter) + 1
704 iii = iii_vet(send_counter)
705 buffer_send(send_counter)%msg(iii) = mat2d(iib - my_start_row + 1, jjb - my_start_col + 1)
706 END DO
707 END DO
708
709 DEALLOCATE (iii_vet)
710 DEALLOCATE (grid_ref_2_send_pos)
711 IF (my_do_release_mat) DEALLOCATE (mat2d)
712 CALL timestop(handle2)
713
714 ! 5) similarly to what done for the buffer_send
715 ! create the buffer for receive, post the message with irecv
716 ! and send the messages non-blocking
717 CALL timeset(routinen//"_isendrecv", handle2)
718 ! count the number of messages to be received
719 number_of_rec = 0
720 DO proc_shift = 1, para_env%num_pe - 1
721 proc_receive = modulo(para_env%mepos - proc_shift, para_env%num_pe)
722 IF (map_rec_size(proc_receive) > 0) THEN
723 number_of_rec = number_of_rec + 1
724 END IF
725 END DO
726
727 ALLOCATE (buffer_rec(number_of_rec))
728
729 rec_counter = 0
730 DO proc_shift = 1, para_env%num_pe - 1
731 proc_receive = modulo(para_env%mepos - proc_shift, para_env%num_pe)
732 size_rec_buffer = map_rec_size(proc_receive)
733 IF (map_rec_size(proc_receive) > 0) THEN
734 rec_counter = rec_counter + 1
735 ! prepare the buffer for receive
736 ALLOCATE (buffer_rec(rec_counter)%msg(size_rec_buffer))
737 buffer_rec(rec_counter)%msg = 0.0_dp
738 buffer_rec(rec_counter)%proc = proc_receive
739 ! post the message to be received
740 CALL para_env%irecv(buffer_rec(rec_counter)%msg, proc_receive, &
741 buffer_rec(rec_counter)%msg_req)
742 END IF
743 END DO
744
745 ! send messages
746 ALLOCATE (req_send(number_of_send))
747 send_counter = 0
748 DO proc_shift = 1, para_env%num_pe - 1
749 proc_send = modulo(para_env%mepos + proc_shift, para_env%num_pe)
750 IF (map_send_size(proc_send) > 0) THEN
751 send_counter = send_counter + 1
752 CALL para_env%isend(buffer_send(send_counter)%msg, proc_send, &
753 buffer_send(send_counter)%msg_req)
754 req_send(send_counter) = buffer_send(send_counter)%msg_req
755 END IF
756 END DO
757 CALL timestop(handle2)
758
759 ! 6) fill the fm_mat matrix with the messages received from the
760 ! other processes
761 CALL timeset(routinen//"_fill", handle2)
762 ! In order to perform this step efficiently first we have to know the
763 ! ranges of the blocks over which a given process hold its local data.
764 ! Start with the rows ...
765 my_num_row_blocks = 1
766 DO iib = 1, nrow_local - 1
767 IF (abs(row_indices(iib + 1) - row_indices(iib)) == 1) cycle
768 my_num_row_blocks = my_num_row_blocks + 1
769 END DO
770 ALLOCATE (blocks_ranges_row(2, my_num_row_blocks))
771 blocks_ranges_row = 0
772 blocks_ranges_row(1, 1) = row_indices(1)
773 iii = 1
774 DO iib = 1, nrow_local - 1
775 IF (abs(row_indices(iib + 1) - row_indices(iib)) == 1) cycle
776 iii = iii + 1
777 blocks_ranges_row(2, iii - 1) = row_indices(iib)
778 blocks_ranges_row(1, iii) = row_indices(iib + 1)
779 END DO
780 blocks_ranges_row(2, my_num_row_blocks) = row_indices(max(nrow_local, 1))
781 ! ... and columns
782 my_num_col_blocks = 1
783 DO jjb = 1, ncol_local - 1
784 IF (abs(col_indices(jjb + 1) - col_indices(jjb)) == 1) cycle
785 my_num_col_blocks = my_num_col_blocks + 1
786 END DO
787 ALLOCATE (blocks_ranges_col(2, my_num_col_blocks))
788 blocks_ranges_col = 0
789 blocks_ranges_col(1, 1) = col_indices(1)
790 iii = 1
791 DO jjb = 1, ncol_local - 1
792 IF (abs(col_indices(jjb + 1) - col_indices(jjb)) == 1) cycle
793 iii = iii + 1
794 blocks_ranges_col(2, iii - 1) = col_indices(jjb)
795 blocks_ranges_col(1, iii) = col_indices(jjb + 1)
796 END DO
797 blocks_ranges_col(2, my_num_col_blocks) = col_indices(max(ncol_local, 1))
798
799 rec_counter = 0
800 DO proc_shift = 1, para_env%num_pe - 1
801 proc_receive = modulo(para_env%mepos - proc_shift, para_env%num_pe)
802 size_rec_buffer = map_rec_size(proc_receive)
803
804 IF (map_rec_size(proc_receive) > 0) THEN
805 rec_counter = rec_counter + 1
806
807 CALL get_group_dist(gd_col, proc_receive, rec_col_start, rec_col_end, rec_col_size)
808
809 ! precompute the number of received columns and relative index
810 num_rec_cols = 0
811 DO jjb = 1, my_num_col_blocks
812 start_col_block = max(blocks_ranges_col(1, jjb), rec_col_start)
813 end_col_block = min(blocks_ranges_col(2, jjb), rec_col_end)
814 DO j_sub = start_col_block, end_col_block
815 num_rec_cols = num_rec_cols + 1
816 END DO
817 END DO
818 ALLOCATE (index_col_rec(num_rec_cols))
819 index_col_rec = 0
820 iii = 0
821 DO jjb = 1, my_num_col_blocks
822 start_col_block = max(blocks_ranges_col(1, jjb), rec_col_start)
823 end_col_block = min(blocks_ranges_col(2, jjb), rec_col_end)
824 DO j_sub = start_col_block, end_col_block
825 iii = iii + 1
826 j_local = cp_fm_indxg2l(j_sub, ncol_block, dummy_proc, &
827 fm_mat%matrix_struct%first_p_pos(2), npcol)
828 index_col_rec(iii) = j_local
829 END DO
830 END DO
831
832 CALL get_group_dist(gd_row, proc_receive, rec_row_start, rec_row_end, rec_row_size)
833
834 ! wait for the message
835 CALL buffer_rec(rec_counter)%msg_req%wait()
836
837 ! fill the local data in fm_mat
838 iii = 0
839 DO iib = 1, my_num_row_blocks
840 start_row_block = max(blocks_ranges_row(1, iib), rec_row_start)
841 end_row_block = min(blocks_ranges_row(2, iib), rec_row_end)
842 DO i_sub = start_row_block, end_row_block
843 i_local = cp_fm_indxg2l(i_sub, nrow_block, dummy_proc, &
844 fm_mat%matrix_struct%first_p_pos(1), nprow)
845 DO jjb = 1, num_rec_cols
846 iii = iii + 1
847 j_local = index_col_rec(jjb)
848 fm_mat%local_data(i_local, j_local) = buffer_rec(rec_counter)%msg(iii)
849 END DO
850 END DO
851 END DO
852
853 DEALLOCATE (buffer_rec(rec_counter)%msg)
854 DEALLOCATE (index_col_rec)
855 END IF
856 END DO
857 DEALLOCATE (buffer_rec)
858
859 DEALLOCATE (blocks_ranges_row)
860 DEALLOCATE (blocks_ranges_col)
861
862 CALL timestop(handle2)
863
864 ! 7) Finally wait for all messeges to be sent
865 CALL timeset(routinen//"_waitall", handle2)
866 CALL mp_waitall(req_send(:))
867 DO send_counter = 1, number_of_send
868 DEALLOCATE (buffer_send(send_counter)%msg)
869 END DO
870 DEALLOCATE (buffer_send)
871 CALL timestop(handle2)
872
873 DEALLOCATE (map_send_size)
874 DEALLOCATE (map_rec_size)
875 DEALLOCATE (grid_2_mepos)
876 DEALLOCATE (mepos_2_grid)
877
878 CALL timestop(handle)
879
880 END SUBROUTINE array2fm
881
882! **************************************************************************************************
883!> \brief redistribute fm to local part of array
884!> \param mat2D ...
885!> \param my_rows ...
886!> \param my_start_row ...
887!> \param my_end_row ...
888!> \param my_cols ...
889!> \param my_start_col ...
890!> \param my_end_col ...
891!> \param group_grid_2_mepos ...
892!> \param mepos_2_grid_group ...
893!> \param ngroup_row ...
894!> \param ngroup_col ...
895!> \param fm_mat ...
896! **************************************************************************************************
897 SUBROUTINE fm2array(mat2D, my_rows, my_start_row, my_end_row, &
898 my_cols, my_start_col, my_end_col, &
899 group_grid_2_mepos, mepos_2_grid_group, &
900 ngroup_row, ngroup_col, &
901 fm_mat)
902
903 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
904 INTENT(OUT) :: mat2d
905 INTEGER, INTENT(IN) :: my_rows, my_start_row, my_end_row, &
906 my_cols, my_start_col, my_end_col
907 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN) :: group_grid_2_mepos, mepos_2_grid_group
908 INTEGER, INTENT(IN) :: ngroup_row, ngroup_col
909 TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat
910
911 CHARACTER(LEN=*), PARAMETER :: routinen = 'fm2array'
912
913 INTEGER :: dummy_proc, handle, handle2, i_global, iib, iii, itmp(2), j_global, jjb, mypcol, &
914 myprow, ncol_block, ncol_local, npcol, nprow, nrow_block, nrow_local, num_cols, &
915 num_rec_rows, num_rows, number_of_rec, number_of_send, proc_receive, proc_send, &
916 proc_shift, rec_col_size, rec_counter, rec_pcol, rec_prow, rec_row_size, ref_send_pcol, &
917 ref_send_prow, send_counter, send_pcol, send_prow, size_rec_buffer, size_send_buffer
918 INTEGER, ALLOCATABLE, DIMENSION(:) :: iii_vet, index_row_rec, map_rec_size, &
919 map_send_size
920 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: grid_2_mepos, grid_ref_2_send_pos, &
921 mepos_2_grid, sizes
922 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
923 REAL(kind=dp) :: part_col, part_row
924 TYPE(integ_mat_buffer_type), ALLOCATABLE, &
925 DIMENSION(:) :: buffer_rec, buffer_send
926 TYPE(mp_para_env_type), POINTER :: para_env
927 TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: req_send
928
929 CALL timeset(routinen, handle)
930
931 ! allocate the array
932 ALLOCATE (mat2d(my_rows, my_cols))
933 mat2d = 0.0_dp
934
935 ! start procedure
936 ! get information of from the full matrix
937 CALL cp_fm_get_info(matrix=fm_mat, &
938 nrow_local=nrow_local, &
939 ncol_local=ncol_local, &
940 row_indices=row_indices, &
941 col_indices=col_indices, &
942 nrow_block=nrow_block, &
943 ncol_block=ncol_block, &
944 nrow_global=num_rows, &
945 ncol_global=num_cols, &
946 para_env=para_env)
947 myprow = fm_mat%matrix_struct%context%mepos(1)
948 mypcol = fm_mat%matrix_struct%context%mepos(2)
949 nprow = fm_mat%matrix_struct%context%num_pe(1)
950 npcol = fm_mat%matrix_struct%context%num_pe(2)
951
952 ! start the communication part
953 ! 0) create array containing the processes position
954 ! and supporting infos
955 CALL timeset(routinen//"_info", handle2)
956 ALLOCATE (grid_2_mepos(0:nprow - 1, 0:npcol - 1))
957 grid_2_mepos = 0
958 ALLOCATE (mepos_2_grid(2, 0:para_env%num_pe - 1))
959
960 ! fill the info array
961 grid_2_mepos(myprow, mypcol) = para_env%mepos
962
963 ! sum infos
964 CALL para_env%sum(grid_2_mepos)
965 CALL para_env%allgather([myprow, mypcol], mepos_2_grid)
966
967 ! 1) loop over my local data and define a map for the proc to send data
968 ALLOCATE (map_send_size(0:para_env%num_pe - 1))
969 map_send_size = 0
970 part_row = real(num_rows, kind=dp)/real(ngroup_row, kind=dp)
971 part_col = real(num_cols, kind=dp)/real(ngroup_col, kind=dp)
972
973 DO jjb = 1, ncol_local
974 j_global = col_indices(jjb)
975 ! check the group holding this element
976 ! dirty way, if someone has a better idea ...
977 send_pcol = int(real(j_global - 1, kind=dp)/part_col)
978 send_pcol = max(0, send_pcol)
979 send_pcol = min(send_pcol, ngroup_col - 1)
980 DO
981 itmp = get_limit(num_cols, ngroup_col, send_pcol)
982 IF (j_global >= itmp(1) .AND. j_global <= itmp(2)) EXIT
983 IF (j_global < itmp(1)) send_pcol = send_pcol - 1
984 IF (j_global > itmp(2)) send_pcol = send_pcol + 1
985 END DO
986
987 DO iib = 1, nrow_local
988 i_global = row_indices(iib)
989 ! check the process in the group holding this element
990 send_prow = int(real(i_global - 1, kind=dp)/part_row)
991 send_prow = max(0, send_prow)
992 send_prow = min(send_prow, ngroup_row - 1)
993 DO
994 itmp = get_limit(num_rows, ngroup_row, send_prow)
995 IF (i_global >= itmp(1) .AND. i_global <= itmp(2)) EXIT
996 IF (i_global < itmp(1)) send_prow = send_prow - 1
997 IF (i_global > itmp(2)) send_prow = send_prow + 1
998 END DO
999
1000 proc_send = group_grid_2_mepos(send_prow, send_pcol)
1001
1002 map_send_size(proc_send) = map_send_size(proc_send) + 1
1003
1004 END DO ! i_global
1005 END DO ! j_global
1006
1007 ! 2) loop over my local data of the array and define a map for the proc from which rec data
1008 ALLOCATE (map_rec_size(0:para_env%num_pe - 1))
1009 map_rec_size = 0
1010 dummy_proc = 0
1011 DO jjb = my_start_col, my_end_col
1012 rec_pcol = cp_fm_indxg2p(jjb, ncol_block, dummy_proc, &
1013 fm_mat%matrix_struct%first_p_pos(2), npcol)
1014 DO iib = my_start_row, my_end_row
1015 rec_prow = cp_fm_indxg2p(iib, nrow_block, dummy_proc, &
1016 fm_mat%matrix_struct%first_p_pos(1), nprow)
1017 proc_receive = grid_2_mepos(rec_prow, rec_pcol)
1018 map_rec_size(proc_receive) = map_rec_size(proc_receive) + 1
1019 END DO
1020 END DO
1021
1022 ! 3) check if the local data in the fm_mat has to be stored in the new array
1023 IF (map_rec_size(para_env%mepos) > 0) THEN
1024 DO jjb = 1, ncol_local
1025 j_global = col_indices(jjb)
1026 IF (j_global >= my_start_col .AND. j_global <= my_end_col) THEN
1027 DO iib = 1, nrow_local
1028 i_global = row_indices(iib)
1029 IF (i_global >= my_start_row .AND. i_global <= my_end_row) THEN
1030 mat2d(i_global - my_start_row + 1, j_global - my_start_col + 1) = fm_mat%local_data(iib, jjb)
1031 END IF
1032 END DO
1033 END IF
1034 END DO
1035 END IF
1036 CALL timestop(handle2)
1037
1038 ! 4) reorder data in the send_buffer
1039 CALL timeset(routinen//"_buffer_send", handle2)
1040 ! count the number of messages to send
1041 number_of_send = 0
1042 DO proc_shift = 1, para_env%num_pe - 1
1043 proc_send = modulo(para_env%mepos + proc_shift, para_env%num_pe)
1044 IF (map_send_size(proc_send) > 0) THEN
1045 number_of_send = number_of_send + 1
1046 END IF
1047 END DO
1048 ! allocate the structure that will hold the messages to be sent
1049 ALLOCATE (buffer_send(number_of_send))
1050 ! grid_ref_2_send_pos is an array (map) that given a pair
1051 ! (ref_send_prow,ref_send_pcol) returns
1052 ! the position in the buffer_send associated to that process
1053
1054 ALLOCATE (grid_ref_2_send_pos(0:ngroup_row - 1, 0:ngroup_col - 1))
1055 grid_ref_2_send_pos = 0
1056
1057 ! finalize the allocation of buffer_send with the actual size
1058 ! of each message (actual size is size_send_buffer)
1059 send_counter = 0
1060 DO proc_shift = 1, para_env%num_pe - 1
1061 proc_send = modulo(para_env%mepos + proc_shift, para_env%num_pe)
1062 size_send_buffer = map_send_size(proc_send)
1063 IF (map_send_size(proc_send) > 0) THEN
1064 send_counter = send_counter + 1
1065 ! allocate the sending buffer (msg)
1066 ALLOCATE (buffer_send(send_counter)%msg(size_send_buffer))
1067 buffer_send(send_counter)%msg = 0.0_dp
1068 buffer_send(send_counter)%proc = proc_send
1069 ! get the pointer to prow, pcol of the process that has
1070 ! to receive this message
1071 ref_send_prow = mepos_2_grid_group(1, proc_send)
1072 ref_send_pcol = mepos_2_grid_group(2, proc_send)
1073 ! save the rank of the process that has to receive this message
1074 grid_ref_2_send_pos(ref_send_prow, ref_send_pcol) = send_counter
1075 END IF
1076 END DO
1077
1078 ! loop over the locally held data and fill the buffer_send
1079 ! for doing that we need an array that keep track if the
1080 ! sequential increase of the index for each message
1081 ALLOCATE (iii_vet(number_of_send))
1082 iii_vet = 0
1083 DO jjb = 1, ncol_local
1084 j_global = col_indices(jjb)
1085 ! check the group holding this element
1086 ! dirty way, if someone has a better idea ...
1087 send_pcol = int(real(j_global - 1, kind=dp)/part_col)
1088 send_pcol = max(0, send_pcol)
1089 send_pcol = min(send_pcol, ngroup_col - 1)
1090 DO
1091 itmp = get_limit(num_cols, ngroup_col, send_pcol)
1092 IF (j_global >= itmp(1) .AND. j_global <= itmp(2)) EXIT
1093 IF (j_global < itmp(1)) send_pcol = send_pcol - 1
1094 IF (j_global > itmp(2)) send_pcol = send_pcol + 1
1095 END DO
1096
1097 DO iib = 1, nrow_local
1098 i_global = row_indices(iib)
1099 ! check the process in the group holding this element
1100 send_prow = int(real(i_global - 1, kind=dp)/part_row)
1101 send_prow = max(0, send_prow)
1102 send_prow = min(send_prow, ngroup_row - 1)
1103 DO
1104 itmp = get_limit(num_rows, ngroup_row, send_prow)
1105 IF (i_global >= itmp(1) .AND. i_global <= itmp(2)) EXIT
1106 IF (i_global < itmp(1)) send_prow = send_prow - 1
1107 IF (i_global > itmp(2)) send_prow = send_prow + 1
1108 END DO
1109 ! we don't need to send to ourselves
1110 IF (group_grid_2_mepos(send_prow, send_pcol) == para_env%mepos) cycle
1111
1112 send_counter = grid_ref_2_send_pos(send_prow, send_pcol)
1113 iii_vet(send_counter) = iii_vet(send_counter) + 1
1114 iii = iii_vet(send_counter)
1115 buffer_send(send_counter)%msg(iii) = fm_mat%local_data(iib, jjb)
1116 END DO
1117 END DO
1118
1119 DEALLOCATE (iii_vet)
1120 DEALLOCATE (grid_ref_2_send_pos)
1121 CALL timestop(handle2)
1122
1123 ! 5) similarly to what done for the buffer_send
1124 ! create the buffer for receive, post the message with irecv
1125 ! and send the messages with non-blocking
1126 CALL timeset(routinen//"_isendrecv", handle2)
1127 ! count the number of messages to be received
1128 number_of_rec = 0
1129 DO proc_shift = 1, para_env%num_pe - 1
1130 proc_receive = modulo(para_env%mepos - proc_shift, para_env%num_pe)
1131 IF (map_rec_size(proc_receive) > 0) THEN
1132 number_of_rec = number_of_rec + 1
1133 END IF
1134 END DO
1135
1136 ALLOCATE (buffer_rec(number_of_rec))
1137
1138 rec_counter = 0
1139 DO proc_shift = 1, para_env%num_pe - 1
1140 proc_receive = modulo(para_env%mepos - proc_shift, para_env%num_pe)
1141 size_rec_buffer = map_rec_size(proc_receive)
1142 IF (map_rec_size(proc_receive) > 0) THEN
1143 rec_counter = rec_counter + 1
1144 ! prepare the buffer for receive
1145 ALLOCATE (buffer_rec(rec_counter)%msg(size_rec_buffer))
1146 buffer_rec(rec_counter)%msg = 0.0_dp
1147 buffer_rec(rec_counter)%proc = proc_receive
1148 ! post the message to be received
1149 CALL para_env%irecv(buffer_rec(rec_counter)%msg, proc_receive, &
1150 buffer_rec(rec_counter)%msg_req)
1151 END IF
1152 END DO
1153
1154 ! send messages
1155 ALLOCATE (req_send(number_of_send))
1156 send_counter = 0
1157 DO proc_shift = 1, para_env%num_pe - 1
1158 proc_send = modulo(para_env%mepos + proc_shift, para_env%num_pe)
1159 IF (map_send_size(proc_send) > 0) THEN
1160 send_counter = send_counter + 1
1161 CALL para_env%isend(buffer_send(send_counter)%msg, proc_send, &
1162 buffer_send(send_counter)%msg_req)
1163 req_send(send_counter) = buffer_send(send_counter)%msg_req
1164 END IF
1165 END DO
1166 CALL timestop(handle2)
1167
1168 ! 6) fill the fm_mat matrix with the messages received from the
1169 ! other processes
1170 CALL timeset(routinen//"_fill", handle2)
1171 CALL cp_fm_get_info(matrix=fm_mat, &
1172 nrow_local=nrow_local, &
1173 ncol_local=ncol_local)
1174 ALLOCATE (sizes(2, 0:para_env%num_pe - 1))
1175 CALL para_env%allgather([nrow_local, ncol_local], sizes)
1176 iib = maxval(sizes(1, :))
1177 ALLOCATE (index_row_rec(iib))
1178 index_row_rec = 0
1179 rec_counter = 0
1180 DO proc_shift = 1, para_env%num_pe - 1
1181 proc_receive = modulo(para_env%mepos - proc_shift, para_env%num_pe)
1182 size_rec_buffer = map_rec_size(proc_receive)
1183
1184 IF (map_rec_size(proc_receive) > 0) THEN
1185 rec_counter = rec_counter + 1
1186
1187 rec_col_size = sizes(2, proc_receive)
1188 rec_row_size = sizes(1, proc_receive)
1189
1190 ! precompute the indeces of the rows
1191 num_rec_rows = 0
1192 DO iib = 1, rec_row_size
1193 i_global = cp_fm_indxl2g(iib, nrow_block, mepos_2_grid(1, proc_receive), &
1194 fm_mat%matrix_struct%first_p_pos(1), nprow)
1195 IF (i_global >= my_start_row .AND. i_global <= my_end_row) THEN
1196 num_rec_rows = num_rec_rows + 1
1197 index_row_rec(num_rec_rows) = i_global
1198 END IF
1199 END DO
1200
1201 CALL buffer_rec(rec_counter)%msg_req%wait()
1202
1203 iii = 0
1204 DO jjb = 1, rec_col_size
1205 j_global = cp_fm_indxl2g(jjb, ncol_block, mepos_2_grid(2, proc_receive), &
1206 fm_mat%matrix_struct%first_p_pos(2), npcol)
1207 IF (j_global >= my_start_col .AND. j_global <= my_end_col) THEN
1208 DO iib = 1, num_rec_rows
1209 i_global = index_row_rec(iib)
1210 iii = iii + 1
1211 mat2d(i_global - my_start_row + 1, j_global - my_start_col + 1) = buffer_rec(rec_counter)%msg(iii)
1212 END DO
1213 END IF
1214 END DO
1215
1216 DEALLOCATE (buffer_rec(rec_counter)%msg)
1217 END IF
1218 END DO
1219 DEALLOCATE (buffer_rec)
1220 DEALLOCATE (index_row_rec)
1221 CALL cp_fm_release(fm_mat)
1222 CALL timestop(handle2)
1223
1224 ! 7) Finally wait for all messeges to be sent
1225 CALL timeset(routinen//"_waitall", handle2)
1226 CALL mp_waitall(req_send(:))
1227 DO send_counter = 1, number_of_send
1228 DEALLOCATE (buffer_send(send_counter)%msg)
1229 END DO
1230 DEALLOCATE (buffer_send)
1231 CALL timestop(handle2)
1232
1233 CALL timestop(handle)
1234
1235 END SUBROUTINE fm2array
1236
1237! **************************************************************************************************
1238!> \brief redistribute 2D representation of 3d tensor to a set of dbcsr matrices
1239!> \param Gamma_2D ...
1240!> \param homo ...
1241!> \param virtual ...
1242!> \param dimen_ia ...
1243!> \param para_env_sub ...
1244!> \param my_ia_start ...
1245!> \param my_ia_end ...
1246!> \param my_group_L_size ...
1247!> \param gd_ia ...
1248!> \param dbcsr_Gamma ...
1249!> \param mo_coeff_o ...
1250! **************************************************************************************************
1251 SUBROUTINE create_dbcsr_gamma(Gamma_2D, homo, virtual, dimen_ia, para_env_sub, &
1252 my_ia_start, my_ia_end, my_group_L_size, &
1253 gd_ia, dbcsr_Gamma, mo_coeff_o)
1254 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gamma_2d
1255 INTEGER, INTENT(IN) :: homo, virtual, dimen_ia
1256 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_sub
1257 INTEGER, INTENT(IN) :: my_ia_start, my_ia_end, my_group_l_size
1258 TYPE(group_dist_d1_type), INTENT(IN) :: gd_ia
1259 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: dbcsr_gamma
1260 TYPE(dbcsr_type), INTENT(INOUT) :: mo_coeff_o
1261
1262 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_dbcsr_gamma'
1263
1264 INTEGER :: dummy_proc, handle, i_global, i_local, iaia, iib, iii, itmp(2), j_global, &
1265 j_local, jjb, jjj, kkb, mypcol, myprow, ncol_block, ncol_local, npcol, nprow, nrow_block, &
1266 nrow_local, number_of_rec, number_of_send, proc_receive, proc_send, proc_shift, &
1267 rec_counter, rec_iaia_end, rec_iaia_size, rec_iaia_start, rec_pcol, rec_prow, &
1268 ref_send_pcol, ref_send_prow, send_counter, send_pcol, send_prow, size_rec_buffer, &
1269 size_send_buffer
1270 INTEGER, ALLOCATABLE, DIMENSION(:) :: iii_vet, map_rec_size, map_send_size
1271 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: grid_2_mepos, grid_ref_2_send_pos, &
1272 indeces_map_my, mepos_2_grid
1273 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1274 REAL(kind=dp) :: part_ia
1275 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1276 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1277 TYPE(cp_fm_type) :: fm_ia
1278 TYPE(index_map), ALLOCATABLE, DIMENSION(:) :: indeces_rec
1279 TYPE(integ_mat_buffer_type), ALLOCATABLE, &
1280 DIMENSION(:) :: buffer_rec, buffer_send
1281 TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: req_send
1282
1283 CALL timeset(routinen, handle)
1284
1285 ! create sub blacs env
1286 NULLIFY (blacs_env)
1287 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env_sub)
1288
1289 ! create the fm_ia buffer matrix
1290 NULLIFY (fm_struct)
1291 CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=homo, &
1292 ncol_global=virtual, para_env=para_env_sub)
1293 CALL cp_fm_create(fm_ia, fm_struct, name="fm_ia")
1294 ! release structure
1295 CALL cp_fm_struct_release(fm_struct)
1296 ! release blacs_env
1297 CALL cp_blacs_env_release(blacs_env)
1298
1299 ! get array information
1300 CALL cp_fm_get_info(matrix=fm_ia, &
1301 nrow_local=nrow_local, &
1302 ncol_local=ncol_local, &
1303 row_indices=row_indices, &
1304 col_indices=col_indices, &
1305 nrow_block=nrow_block, &
1306 ncol_block=ncol_block)
1307 myprow = fm_ia%matrix_struct%context%mepos(1)
1308 mypcol = fm_ia%matrix_struct%context%mepos(2)
1309 nprow = fm_ia%matrix_struct%context%num_pe(1)
1310 npcol = fm_ia%matrix_struct%context%num_pe(2)
1311
1312 ! 0) create array containing the processes position and supporting infos
1313 ALLOCATE (grid_2_mepos(0:nprow - 1, 0:npcol - 1))
1314 grid_2_mepos = 0
1315 ALLOCATE (mepos_2_grid(2, 0:para_env_sub%num_pe - 1))
1316 ! fill the info array
1317 grid_2_mepos(myprow, mypcol) = para_env_sub%mepos
1318 ! sum infos
1319 CALL para_env_sub%sum(grid_2_mepos)
1320 CALL para_env_sub%allgather([myprow, mypcol], mepos_2_grid)
1321
1322 ! loop over local index range and define the sending map
1323 ALLOCATE (map_send_size(0:para_env_sub%num_pe - 1))
1324 map_send_size = 0
1325 dummy_proc = 0
1326 DO iaia = my_ia_start, my_ia_end
1327 i_global = (iaia - 1)/virtual + 1
1328 j_global = mod(iaia - 1, virtual) + 1
1329 send_prow = cp_fm_indxg2p(i_global, nrow_block, dummy_proc, &
1330 fm_ia%matrix_struct%first_p_pos(1), nprow)
1331 send_pcol = cp_fm_indxg2p(j_global, ncol_block, dummy_proc, &
1332 fm_ia%matrix_struct%first_p_pos(2), npcol)
1333 proc_send = grid_2_mepos(send_prow, send_pcol)
1334 map_send_size(proc_send) = map_send_size(proc_send) + 1
1335 END DO
1336
1337 ! loop over local data of fm_ia and define the receiving map
1338 ALLOCATE (map_rec_size(0:para_env_sub%num_pe - 1))
1339 map_rec_size = 0
1340 part_ia = real(dimen_ia, kind=dp)/real(para_env_sub%num_pe, kind=dp)
1341
1342 DO iib = 1, nrow_local
1343 i_global = row_indices(iib)
1344 DO jjb = 1, ncol_local
1345 j_global = col_indices(jjb)
1346 iaia = (i_global - 1)*virtual + j_global
1347 proc_receive = int(real(iaia - 1, kind=dp)/part_ia)
1348 proc_receive = max(0, proc_receive)
1349 proc_receive = min(proc_receive, para_env_sub%num_pe - 1)
1350 DO
1351 itmp = get_limit(dimen_ia, para_env_sub%num_pe, proc_receive)
1352 IF (iaia >= itmp(1) .AND. iaia <= itmp(2)) EXIT
1353 IF (iaia < itmp(1)) proc_receive = proc_receive - 1
1354 IF (iaia > itmp(2)) proc_receive = proc_receive + 1
1355 END DO
1356 map_rec_size(proc_receive) = map_rec_size(proc_receive) + 1
1357 END DO
1358 END DO
1359
1360 ! allocate the buffer for sending data
1361 number_of_send = 0
1362 DO proc_shift = 1, para_env_sub%num_pe - 1
1363 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
1364 IF (map_send_size(proc_send) > 0) THEN
1365 number_of_send = number_of_send + 1
1366 END IF
1367 END DO
1368 ! allocate the structure that will hold the messages to be sent
1369 ALLOCATE (buffer_send(number_of_send))
1370 ! and the map from the grid of processess to the message position
1371 ALLOCATE (grid_ref_2_send_pos(0:nprow - 1, 0:npcol - 1))
1372 grid_ref_2_send_pos = 0
1373 ! finally allocate each message
1374 send_counter = 0
1375 DO proc_shift = 1, para_env_sub%num_pe - 1
1376 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
1377 size_send_buffer = map_send_size(proc_send)
1378 IF (map_send_size(proc_send) > 0) THEN
1379 send_counter = send_counter + 1
1380 ! allocate the sending buffer (msg)
1381 ALLOCATE (buffer_send(send_counter)%msg(size_send_buffer))
1382 buffer_send(send_counter)%proc = proc_send
1383 ! get the pointer to prow, pcol of the process that has
1384 ! to receive this message
1385 ref_send_prow = mepos_2_grid(1, proc_send)
1386 ref_send_pcol = mepos_2_grid(2, proc_send)
1387 ! save the rank of the process that has to receive this message
1388 grid_ref_2_send_pos(ref_send_prow, ref_send_pcol) = send_counter
1389 END IF
1390 END DO
1391
1392 ! allocate the buffer for receiving data
1393 number_of_rec = 0
1394 DO proc_shift = 1, para_env_sub%num_pe - 1
1395 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1396 IF (map_rec_size(proc_receive) > 0) THEN
1397 number_of_rec = number_of_rec + 1
1398 END IF
1399 END DO
1400 ! allocate the structure that will hold the messages to be received
1401 ! and relative indeces
1402 ALLOCATE (buffer_rec(number_of_rec))
1403 ALLOCATE (indeces_rec(number_of_rec))
1404 ! finally allocate each message and fill the array of indeces
1405 rec_counter = 0
1406 DO proc_shift = 1, para_env_sub%num_pe - 1
1407 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1408 size_rec_buffer = map_rec_size(proc_receive)
1409 IF (map_rec_size(proc_receive) > 0) THEN
1410 rec_counter = rec_counter + 1
1411 ! prepare the buffer for receive
1412 ALLOCATE (buffer_rec(rec_counter)%msg(size_rec_buffer))
1413 buffer_rec(rec_counter)%proc = proc_receive
1414 ! create the indeces array
1415 ALLOCATE (indeces_rec(rec_counter)%map(2, size_rec_buffer))
1416 indeces_rec(rec_counter)%map = 0
1417 CALL get_group_dist(gd_ia, proc_receive, rec_iaia_start, rec_iaia_end, rec_iaia_size)
1418 iii = 0
1419 DO iaia = rec_iaia_start, rec_iaia_end
1420 i_global = (iaia - 1)/virtual + 1
1421 j_global = mod(iaia - 1, virtual) + 1
1422 rec_prow = cp_fm_indxg2p(i_global, nrow_block, dummy_proc, &
1423 fm_ia%matrix_struct%first_p_pos(1), nprow)
1424 rec_pcol = cp_fm_indxg2p(j_global, ncol_block, dummy_proc, &
1425 fm_ia%matrix_struct%first_p_pos(2), npcol)
1426 IF (grid_2_mepos(rec_prow, rec_pcol) /= para_env_sub%mepos) cycle
1427 iii = iii + 1
1428 i_local = cp_fm_indxg2l(i_global, nrow_block, dummy_proc, &
1429 fm_ia%matrix_struct%first_p_pos(1), nprow)
1430 j_local = cp_fm_indxg2l(j_global, ncol_block, dummy_proc, &
1431 fm_ia%matrix_struct%first_p_pos(2), npcol)
1432 indeces_rec(rec_counter)%map(1, iii) = i_local
1433 indeces_rec(rec_counter)%map(2, iii) = j_local
1434 END DO
1435 END IF
1436 END DO
1437 ! and create the index map for my local data
1438 IF (map_rec_size(para_env_sub%mepos) > 0) THEN
1439 size_rec_buffer = map_rec_size(para_env_sub%mepos)
1440 ALLOCATE (indeces_map_my(2, size_rec_buffer))
1441 indeces_map_my = 0
1442 iii = 0
1443 DO iaia = my_ia_start, my_ia_end
1444 i_global = (iaia - 1)/virtual + 1
1445 j_global = mod(iaia - 1, virtual) + 1
1446 rec_prow = cp_fm_indxg2p(i_global, nrow_block, dummy_proc, &
1447 fm_ia%matrix_struct%first_p_pos(1), nprow)
1448 rec_pcol = cp_fm_indxg2p(j_global, ncol_block, dummy_proc, &
1449 fm_ia%matrix_struct%first_p_pos(2), npcol)
1450 IF (grid_2_mepos(rec_prow, rec_pcol) /= para_env_sub%mepos) cycle
1451 iii = iii + 1
1452 i_local = cp_fm_indxg2l(i_global, nrow_block, dummy_proc, &
1453 fm_ia%matrix_struct%first_p_pos(1), nprow)
1454 j_local = cp_fm_indxg2l(j_global, ncol_block, dummy_proc, &
1455 fm_ia%matrix_struct%first_p_pos(2), npcol)
1456 indeces_map_my(1, iii) = i_local
1457 indeces_map_my(2, iii) = j_local
1458 END DO
1459 END IF
1460
1461 ! auxiliary vector of indeces for the send buffer
1462 ALLOCATE (iii_vet(number_of_send))
1463 ! vector for the send requests
1464 ALLOCATE (req_send(number_of_send))
1465 ! loop over auxiliary basis function and redistribute into a fm
1466 ! and then compy the fm into a dbcsr matrix
1467 DO kkb = 1, my_group_l_size
1468 ! zero the matries of the buffers and post the messages to be received
1469 CALL cp_fm_set_all(matrix=fm_ia, alpha=0.0_dp)
1470 rec_counter = 0
1471 DO proc_shift = 1, para_env_sub%num_pe - 1
1472 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1473 IF (map_rec_size(proc_receive) > 0) THEN
1474 rec_counter = rec_counter + 1
1475 buffer_rec(rec_counter)%msg = 0.0_dp
1476 CALL para_env_sub%irecv(buffer_rec(rec_counter)%msg, proc_receive, &
1477 buffer_rec(rec_counter)%msg_req)
1478 END IF
1479 END DO
1480 ! fill the sending buffer and send the messages
1481 DO send_counter = 1, number_of_send
1482 buffer_send(send_counter)%msg = 0.0_dp
1483 END DO
1484 iii_vet = 0
1485 jjj = 0
1486 DO iaia = my_ia_start, my_ia_end
1487 i_global = (iaia - 1)/virtual + 1
1488 j_global = mod(iaia - 1, virtual) + 1
1489 send_prow = cp_fm_indxg2p(i_global, nrow_block, dummy_proc, &
1490 fm_ia%matrix_struct%first_p_pos(1), nprow)
1491 send_pcol = cp_fm_indxg2p(j_global, ncol_block, dummy_proc, &
1492 fm_ia%matrix_struct%first_p_pos(2), npcol)
1493 proc_send = grid_2_mepos(send_prow, send_pcol)
1494 ! we don't need to send to ourselves
1495 IF (grid_2_mepos(send_prow, send_pcol) == para_env_sub%mepos) THEN
1496 ! filling fm_ia with local data
1497 jjj = jjj + 1
1498 i_local = indeces_map_my(1, jjj)
1499 j_local = indeces_map_my(2, jjj)
1500 fm_ia%local_data(i_local, j_local) = gamma_2d(iaia - my_ia_start + 1, kkb)
1501 ELSE
1502 send_counter = grid_ref_2_send_pos(send_prow, send_pcol)
1503 iii_vet(send_counter) = iii_vet(send_counter) + 1
1504 iii = iii_vet(send_counter)
1505 buffer_send(send_counter)%msg(iii) = gamma_2d(iaia - my_ia_start + 1, kkb)
1506 END IF
1507 END DO
1508 req_send = mp_request_null
1509 send_counter = 0
1510 DO proc_shift = 1, para_env_sub%num_pe - 1
1511 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
1512 IF (map_send_size(proc_send) > 0) THEN
1513 send_counter = send_counter + 1
1514 CALL para_env_sub%isend(buffer_send(send_counter)%msg, proc_send, &
1515 buffer_send(send_counter)%msg_req)
1516 req_send(send_counter) = buffer_send(send_counter)%msg_req
1517 END IF
1518 END DO
1519
1520 ! receive the massages and fill the fm_ia
1521 rec_counter = 0
1522 DO proc_shift = 1, para_env_sub%num_pe - 1
1523 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1524 size_rec_buffer = map_rec_size(proc_receive)
1525 IF (map_rec_size(proc_receive) > 0) THEN
1526 rec_counter = rec_counter + 1
1527 ! wait for the message
1528 CALL buffer_rec(rec_counter)%msg_req%wait()
1529 DO iii = 1, size_rec_buffer
1530 i_local = indeces_rec(rec_counter)%map(1, iii)
1531 j_local = indeces_rec(rec_counter)%map(2, iii)
1532 fm_ia%local_data(i_local, j_local) = buffer_rec(rec_counter)%msg(iii)
1533 END DO
1534 END IF
1535 END DO
1536
1537 ! wait all
1538 CALL mp_waitall(req_send(:))
1539
1540 ! now create the DBCSR matrix and copy fm_ia into it
1541 ALLOCATE (dbcsr_gamma(kkb)%matrix)
1542 CALL cp_dbcsr_m_by_n_from_template(dbcsr_gamma(kkb)%matrix, &
1543 template=mo_coeff_o, &
1544 m=homo, n=virtual, sym=dbcsr_type_no_symmetry)
1545 CALL copy_fm_to_dbcsr(fm_ia, dbcsr_gamma(kkb)%matrix, keep_sparsity=.false.)
1546
1547 END DO
1548
1549 ! deallocate stuff
1550 DEALLOCATE (gamma_2d)
1551 DEALLOCATE (iii_vet)
1552 DEALLOCATE (req_send)
1553 IF (map_rec_size(para_env_sub%mepos) > 0) THEN
1554 DEALLOCATE (indeces_map_my)
1555 END IF
1556 DO rec_counter = 1, number_of_rec
1557 DEALLOCATE (indeces_rec(rec_counter)%map)
1558 DEALLOCATE (buffer_rec(rec_counter)%msg)
1559 END DO
1560 DEALLOCATE (indeces_rec)
1561 DEALLOCATE (buffer_rec)
1562 DO send_counter = 1, number_of_send
1563 DEALLOCATE (buffer_send(send_counter)%msg)
1564 END DO
1565 DEALLOCATE (buffer_send)
1566 DEALLOCATE (map_send_size)
1567 DEALLOCATE (map_rec_size)
1568 DEALLOCATE (grid_2_mepos)
1569 DEALLOCATE (mepos_2_grid)
1570
1571 ! release buffer matrix
1572 CALL cp_fm_release(fm_ia)
1573
1574 CALL timestop(handle)
1575
1576 END SUBROUTINE create_dbcsr_gamma
1577
1578! **************************************************************************************************
1579!> \brief prepare array for redistribution
1580!> \param para_env ...
1581!> \param para_env_sub ...
1582!> \param ngroup ...
1583!> \param group_grid_2_mepos ...
1584!> \param mepos_2_grid_group ...
1585!> \param pos_info ...
1586! **************************************************************************************************
1587 SUBROUTINE prepare_redistribution(para_env, para_env_sub, ngroup, &
1588 group_grid_2_mepos, mepos_2_grid_group, &
1589 pos_info)
1590 TYPE(mp_para_env_type), INTENT(IN) :: para_env, para_env_sub
1591 INTEGER, INTENT(IN) :: ngroup
1592 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: group_grid_2_mepos, mepos_2_grid_group
1593 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT), &
1594 OPTIONAL :: pos_info
1595
1596 INTEGER :: i, pos_group, pos_sub
1597 INTEGER, ALLOCATABLE, DIMENSION(:) :: my_pos_info
1598
1599 ALLOCATE (my_pos_info(0:para_env%num_pe - 1))
1600 CALL para_env%allgather(para_env_sub%mepos, my_pos_info)
1601
1602 ALLOCATE (group_grid_2_mepos(0:para_env_sub%num_pe - 1, 0:ngroup - 1))
1603 group_grid_2_mepos = 0
1604 ALLOCATE (mepos_2_grid_group(2, 0:para_env%num_pe - 1))
1605 mepos_2_grid_group = 0
1606
1607 DO i = 0, para_env%num_pe - 1
1608 ! calculate position of the group
1609 pos_group = i/para_env_sub%num_pe
1610 ! calculate position in the subgroup
1611 pos_sub = my_pos_info(i)
1612 ! fill the map from the grid of groups to process
1613 group_grid_2_mepos(pos_sub, pos_group) = i
1614 ! and the opposite, from the global pos to the grid pos
1615 mepos_2_grid_group(1, i) = pos_sub
1616 mepos_2_grid_group(2, i) = pos_group
1617 END DO
1618
1619 IF (PRESENT(pos_info)) CALL move_alloc(my_pos_info, pos_info)
1620
1621 END SUBROUTINE prepare_redistribution
1622
1623END MODULE mp2_ri_grad_util
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:117
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
group_dist_types::create_group_dist
Definition group_dist_types.F:39
group_dist_types::get_group_dist
Definition group_dist_types.F:47
group_dist_types::release_group_dist
Definition group_dist_types.F:53
message_passing::mp_waitall
Definition message_passing.F:871
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_blacs_env::cp_blacs_env_release
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
Definition cp_blacs_env.F:282
cp_blacs_env::cp_blacs_env_create
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Definition cp_blacs_env.F:123
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_m_by_n_from_template
subroutine, public cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym, data_type)
Utility function to create an arbitrary shaped dbcsr matrix with the same processor grid as the templ...
Definition cp_dbcsr_operations.F:1105
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:118
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_geadd
subroutine, public cp_fm_geadd(alpha, trans, matrix_a, beta, matrix_b)
interface to BLACS geadd: matrix_b = beta*matrix_b + alpha*opt(matrix_a) where opt(matrix_a) can be e...
Definition cp_fm_basic_linalg.F:274
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_get
subroutine, public cp_fm_struct_get(fmstruct, para_env, context, descriptor, ncol_block, nrow_block, nrow_global, ncol_global, first_p_pos, row_indices, col_indices, nrow_local, ncol_local, nrow_locals, ncol_locals, local_leading_dimension)
returns the values of various attributes of the matrix structure
Definition cp_fm_struct.F:409
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_indxl2g
integer function, public cp_fm_indxl2g(indxloc, nb, iproc, isrcproc, nprocs)
wrapper to scalapack function INDXL2G that computes the global index of a distributed matrix entry po...
Definition cp_fm_types.F:2585
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1016
cp_fm_types::cp_fm_indxg2l
integer function, public cp_fm_indxg2l(indxglob, nb, iproc, isrcproc, nprocs)
wrapper to scalapack function INDXG2L that computes the local index of a distributed matrix entry poi...
Definition cp_fm_types.F:2525
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:535
cp_fm_types::cp_fm_indxg2p
integer function, public cp_fm_indxg2p(indxglob, nb, iproc, isrcproc, nprocs)
wrapper to scalapack function INDXG2P that computes the process coordinate which possesses the entry ...
Definition cp_fm_types.F:2466
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
group_dist_types
Types to describe group distributions.
Definition group_dist_types.F:14
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
libint_2c_3c
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition libint_2c_3c.F:14
libint_2c_3c::compare_potential_types
pure logical function, public compare_potential_types(potential1, potential2)
Helper function to compare libint_potential_types.
Definition libint_2c_3c.F:963
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
message_passing::mp_request_null
type(mp_request_type), parameter, public mp_request_null
Definition message_passing.F:814
mp2_ri_grad_util
Routines for calculating RI-MP2 gradients.
Definition mp2_ri_grad_util.F:13
mp2_ri_grad_util::complete_gamma
subroutine, public complete_gamma(mp2_env, b_ia_q, dimen_ri, homo, virtual, para_env, para_env_sub, ngroup, my_group_l_size, my_group_l_start, my_group_l_end, my_b_size, my_b_virtual_start, gd_array, gd_b_virtual, kspin)
complete the calculation of the Gamma matrices
Definition mp2_ri_grad_util.F:84
mp2_ri_grad_util::prepare_redistribution
subroutine, public prepare_redistribution(para_env, para_env_sub, ngroup, group_grid_2_mepos, mepos_2_grid_group, pos_info)
prepare array for redistribution
Definition mp2_ri_grad_util.F:1590
mp2_ri_grad_util::array2fm
subroutine, public array2fm(mat2d, fm_struct, my_start_row, my_end_row, my_start_col, my_end_col, gd_row, gd_col, group_grid_2_mepos, ngroup_row, ngroup_col, fm_mat, integ_group_size, color_group, do_release_mat)
redistribute local part of array to fm
Definition mp2_ri_grad_util.F:495
mp2_ri_grad_util::create_dbcsr_gamma
subroutine, public create_dbcsr_gamma(gamma_2d, homo, virtual, dimen_ia, para_env_sub, my_ia_start, my_ia_end, my_group_l_size, gd_ia, dbcsr_gamma, mo_coeff_o)
redistribute 2D representation of 3d tensor to a set of dbcsr matrices
Definition mp2_ri_grad_util.F:1254
mp2_ri_grad_util::fm2array
subroutine, public fm2array(mat2d, my_rows, my_start_row, my_end_row, my_cols, my_start_col, my_end_col, group_grid_2_mepos, mepos_2_grid_group, ngroup_row, ngroup_col, fm_mat)
redistribute fm to local part of array
Definition mp2_ri_grad_util.F:902
mp2_types
Types needed for MP2 calculations.
Definition mp2_types.F:14
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
util
All kind of helpful little routines.
Definition util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition util.F:333
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
group_dist_types::group_dist_d1_type
Definition group_dist_types.F:35
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
message_passing::mp_request_type
Definition message_passing.F:623
mp2_types::integ_mat_buffer_type
Definition mp2_types.F:377
mp2_types::mp2_type
Definition mp2_types.F:317