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 : !> \note
10 : !> Basic type for real space grid methods
11 : !> \par History
12 : !> JGH (22-May-2002) : New routine rs_grid_zero
13 : !> JGH (12-Jun-2002) : Bug fix for mpi groups
14 : !> JGH (19-Jun-2003) : Added routine for task distribution
15 : !> JGH (23-Nov-2003) : Added routine for task loop separation
16 : !> \author JGH (18-Mar-2001)
17 : ! **************************************************************************************************
18 : MODULE realspace_grid_types
19 : USE cp_array_utils, ONLY: cp_1d_r_p_type
20 : USE cp_log_handling, ONLY: cp_to_string
21 : USE kahan_sum, ONLY: accurate_sum
22 : USE kinds, ONLY: dp,&
23 : int_8
24 : USE machine, ONLY: m_memory
25 : USE mathlib, ONLY: det_3x3
26 : USE message_passing, ONLY: mp_comm_null,&
27 : mp_comm_type,&
28 : mp_request_null,&
29 : mp_request_type,&
30 : mp_waitall,&
31 : mp_waitany
32 : USE offload_api, ONLY: offload_buffer_type,&
33 : offload_create_buffer,&
34 : offload_free_buffer
35 : USE pw_grid_types, ONLY: PW_MODE_LOCAL,&
36 : pw_grid_type
37 : USE pw_grids, ONLY: pw_grid_release,&
38 : pw_grid_retain
39 : USE pw_methods, ONLY: pw_integrate_function
40 : USE pw_types, ONLY: pw_r3d_rs_type
41 : USE util, ONLY: get_limit
42 :
43 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
44 :
45 : #include "../base/base_uses.f90"
46 :
47 : IMPLICIT NONE
48 :
49 : PRIVATE
50 : PUBLIC :: realspace_grid_type, &
51 : realspace_grid_desc_type, &
52 : realspace_grid_p_type, &
53 : realspace_grid_desc_p_type, &
54 : realspace_grid_input_type
55 :
56 : PUBLIC :: transfer_rs2pw, &
57 : transfer_pw2rs, &
58 : rs_grid_zero, &
59 : rs_grid_set_box, &
60 : rs_grid_create, &
61 : rs_grid_create_descriptor, &
62 : rs_grid_retain_descriptor, &
63 : rs_grid_release, &
64 : rs_grid_release_descriptor, &
65 : rs_grid_reorder_ranks, &
66 : rs_grid_print, &
67 : rs_grid_locate_rank, &
68 : rs_grid_max_ngpts, &
69 : rs_grid_mult_and_add, &
70 : map_gaussian_here
71 :
72 : INTEGER, PARAMETER, PUBLIC :: rsgrid_distributed = 0, &
73 : rsgrid_replicated = 1, &
74 : rsgrid_automatic = 2
75 :
76 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
77 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'realspace_grid_types'
78 :
79 : ! **************************************************************************************************
80 : TYPE realspace_grid_input_type
81 : INTEGER :: distribution_type = rsgrid_replicated
82 : INTEGER :: distribution_layout(3) = -1
83 : REAL(KIND=dp) :: memory_factor = 0.0_dp
84 : LOGICAL :: lock_distribution = .FALSE.
85 : INTEGER :: nsmax = -1
86 : REAL(KIND=dp) :: halo_reduction_factor = 1.0_dp
87 : END TYPE realspace_grid_input_type
88 :
89 : ! **************************************************************************************************
90 : TYPE realspace_grid_desc_type
91 : TYPE(pw_grid_type), POINTER :: pw => NULL() ! the pw grid
92 :
93 : INTEGER :: ref_count = 0 ! reference count
94 :
95 : INTEGER(int_8) :: ngpts = 0_int_8 ! # grid points
96 : INTEGER, DIMENSION(3) :: npts = 0 ! # grid points per dimension
97 : INTEGER, DIMENSION(3) :: lb = 0 ! lower bounds
98 : INTEGER, DIMENSION(3) :: ub = 0 ! upper bounds
99 :
100 : INTEGER :: border = 0 ! border points
101 :
102 : INTEGER, DIMENSION(3) :: perd = -1 ! periodicity enforced
103 : REAL(KIND=dp), DIMENSION(3, 3) :: dh = 0.0_dp ! incremental grid matrix
104 : REAL(KIND=dp), DIMENSION(3, 3) :: dh_inv = 0.0_dp ! inverse incremental grid matrix
105 : LOGICAL :: orthorhombic = .TRUE. ! grid symmetry
106 :
107 : LOGICAL :: parallel = .TRUE. ! whether the corresponding pw grid is distributed
108 : LOGICAL :: distributed = .TRUE. ! whether the rs grid is distributed
109 : ! these MPI related quantities are only meaningful depending on how the grid has been laid out
110 : ! they are most useful for fully distributed grids, where they reflect the topology of the grid
111 : TYPE(mp_comm_type) :: group = mp_comm_null
112 : INTEGER :: my_pos = -1
113 : INTEGER :: group_size = 0
114 : INTEGER, DIMENSION(3) :: group_dim = -1
115 : INTEGER, DIMENSION(3) :: group_coor = -1
116 : INTEGER, DIMENSION(3) :: neighbours = -1
117 : ! only meaningful on distributed grids
118 : ! a list of bounds for each CPU
119 : INTEGER, DIMENSION(:, :), ALLOCATABLE :: lb_global
120 : INTEGER, DIMENSION(:, :), ALLOCATABLE :: ub_global
121 : ! a mapping from linear rank to 3d coord
122 : INTEGER, DIMENSION(:, :), ALLOCATABLE :: rank2coord
123 : INTEGER, DIMENSION(:, :, :), ALLOCATABLE :: coord2rank
124 : ! a mapping from index to rank (which allows to figure out easily on which rank a given point of the grid is)
125 : INTEGER, DIMENSION(:), ALLOCATABLE :: x2coord
126 : INTEGER, DIMENSION(:), ALLOCATABLE :: y2coord
127 : INTEGER, DIMENSION(:), ALLOCATABLE :: z2coord
128 :
129 : INTEGER :: my_virtual_pos = -1
130 : INTEGER, DIMENSION(3) :: virtual_group_coor = -1
131 :
132 : INTEGER, DIMENSION(:), ALLOCATABLE :: virtual2real, real2virtual
133 :
134 : END TYPE realspace_grid_desc_type
135 :
136 : TYPE realspace_grid_type
137 :
138 : TYPE(realspace_grid_desc_type), POINTER :: desc => NULL()
139 :
140 : INTEGER :: ngpts_local = -1 ! local dimensions
141 : INTEGER, DIMENSION(3) :: npts_local = -1
142 : INTEGER, DIMENSION(3) :: lb_local = -1
143 : INTEGER, DIMENSION(3) :: ub_local = -1
144 : INTEGER, DIMENSION(3) :: lb_real = -1 ! lower bounds of the real local data
145 : INTEGER, DIMENSION(3) :: ub_real = -1 ! upper bounds of the real local data
146 :
147 : INTEGER, DIMENSION(:), ALLOCATABLE :: px, py, pz ! index translators
148 : TYPE(offload_buffer_type) :: buffer = offload_buffer_type() ! owner of the grid's memory
149 : REAL(KIND=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: r => NULL() ! the grid (pointer to buffer%host_buffer)
150 :
151 : END TYPE realspace_grid_type
152 :
153 : ! **************************************************************************************************
154 : TYPE realspace_grid_p_type
155 : TYPE(realspace_grid_type), POINTER :: rs_grid => NULL()
156 : END TYPE realspace_grid_p_type
157 :
158 : TYPE realspace_grid_desc_p_type
159 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc => NULL()
160 : END TYPE realspace_grid_desc_p_type
161 :
162 : CONTAINS
163 :
164 : ! **************************************************************************************************
165 : !> \brief returns the 1D rank of the task which is a cartesian shift away from 1D rank rank_in
166 : !> only possible if rs_grid is a distributed grid
167 : !> \param rs_desc ...
168 : !> \param rank_in ...
169 : !> \param shift ...
170 : !> \return ...
171 : ! **************************************************************************************************
172 2346 : PURE FUNCTION rs_grid_locate_rank(rs_desc, rank_in, shift) RESULT(rank_out)
173 : TYPE(realspace_grid_desc_type), INTENT(IN) :: rs_desc
174 : INTEGER, INTENT(IN) :: rank_in
175 : INTEGER, DIMENSION(3), INTENT(IN) :: shift
176 : INTEGER :: rank_out
177 :
178 : INTEGER :: coord(3)
179 :
180 9384 : coord = MODULO(rs_desc%rank2coord(:, rank_in) + shift, rs_desc%group_dim)
181 2346 : rank_out = rs_desc%coord2rank(coord(1), coord(2), coord(3))
182 2346 : END FUNCTION rs_grid_locate_rank
183 :
184 : ! **************************************************************************************************
185 : !> \brief Determine the setup of real space grids - this is divided up into the
186 : !> creation of a descriptor and the actual grid itself (see rs_grid_create)
187 : !> \param desc ...
188 : !> \param pw_grid ...
189 : !> \param input_settings ...
190 : !> \param border_points ...
191 : !> \par History
192 : !> JGH (08-Jun-2003) : nsmax <= 0 indicates fully replicated grid
193 : !> Iain Bethune (05-Sep-2008) : modified cut heuristic
194 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2008 on behalf of the HECToR project
195 : !> - Create a descriptor for realspace grids with a number of border
196 : !> points as exactly given by the optional argument border_points.
197 : !> These grids are always distributed.
198 : !> (27.11.2013, Matthias Krack)
199 : !> \author JGH (18-Mar-2001)
200 : ! **************************************************************************************************
201 41330 : SUBROUTINE rs_grid_create_descriptor(desc, pw_grid, input_settings, border_points)
202 : TYPE(realspace_grid_desc_type), POINTER :: desc
203 : TYPE(pw_grid_type), INTENT(INOUT), TARGET :: pw_grid
204 : TYPE(realspace_grid_input_type), INTENT(IN) :: input_settings
205 : INTEGER, INTENT(IN), OPTIONAL :: border_points
206 :
207 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rs_grid_create_descriptor'
208 :
209 : INTEGER :: border_size, dir, handle, i, j, k, l, &
210 : lb(2), min_npts_real, n_slices(3), &
211 : n_slices_tmp(3), nmin
212 : LOGICAL :: overlap
213 : REAL(KIND=dp) :: ratio, ratio_best, volume, volume_dist
214 :
215 41330 : CALL timeset(routineN, handle)
216 :
217 41330 : IF (PRESENT(border_points)) THEN
218 128 : border_size = border_points
219 : ELSE
220 : border_size = 0
221 : END IF
222 :
223 2025170 : ALLOCATE (desc)
224 :
225 41330 : CALL pw_grid%para%group%sync()
226 :
227 41330 : desc%pw => pw_grid
228 41330 : CALL pw_grid_retain(desc%pw)
229 :
230 537290 : desc%dh = pw_grid%dh
231 537290 : desc%dh_inv = pw_grid%dh_inv
232 41330 : desc%orthorhombic = pw_grid%orthorhombic
233 41330 : desc%ref_count = 1
234 :
235 41330 : IF (pw_grid%para%mode == PW_MODE_LOCAL) THEN
236 : ! The corresponding group has dimension 1
237 : ! All operations will be done locally
238 16912 : desc%npts = pw_grid%npts
239 16912 : desc%ngpts = PRODUCT(INT(desc%npts, KIND=int_8))
240 16912 : desc%lb = pw_grid%bounds(1, :)
241 16912 : desc%ub = pw_grid%bounds(2, :)
242 4228 : desc%border = border_size
243 4228 : IF (border_size == 0) THEN
244 16912 : desc%perd = 1
245 : ELSE
246 0 : desc%perd = 0
247 : END IF
248 4228 : desc%parallel = .FALSE.
249 4228 : desc%distributed = .FALSE.
250 4228 : desc%group = mp_comm_null
251 4228 : desc%group_size = 1
252 16912 : desc%group_dim = 1
253 16912 : desc%group_coor = 0
254 4228 : desc%my_pos = 0
255 : ELSE
256 : ! group size of desc grid
257 : ! global grid dimensions are still the same
258 37102 : desc%group_size = pw_grid%para%group%num_pe
259 148408 : desc%npts = pw_grid%npts
260 148408 : desc%ngpts = PRODUCT(INT(desc%npts, KIND=int_8))
261 148408 : desc%lb = pw_grid%bounds(1, :)
262 148408 : desc%ub = pw_grid%bounds(2, :)
263 :
264 : ! this is the eventual border size
265 37102 : IF (border_size == 0) THEN
266 36974 : nmin = (input_settings%nsmax + 1)/2
267 36974 : nmin = MAX(0, NINT(nmin*input_settings%halo_reduction_factor))
268 : ELSE
269 : ! Set explicitly the requested border size
270 : nmin = border_size
271 : END IF
272 :
273 37102 : IF (input_settings%distribution_type == rsgrid_replicated) THEN
274 :
275 57992 : n_slices = 1
276 14498 : IF (border_size > 0) THEN
277 : CALL cp_abort(__LOCATION__, &
278 : "An explicit border size > 0 is not yet working for "// &
279 : "replicated realspace grids. Request DISTRIBUTION_TYPE "// &
280 0 : "distributed for RS_GRID explicitly.")
281 : END IF
282 :
283 : ELSE
284 :
285 90416 : n_slices = 1
286 22604 : ratio_best = -HUGE(ratio_best)
287 :
288 : ! don't allow distributions with more processors than real grid points
289 67812 : DO k = 1, MIN(desc%npts(3), desc%group_size)
290 158228 : DO j = 1, MIN(desc%npts(2), desc%group_size)
291 90416 : i = MIN(desc%npts(1), desc%group_size/(j*k))
292 361664 : n_slices_tmp = [i, j, k]
293 :
294 : ! we don't match the actual number of CPUs
295 361664 : IF (PRODUCT(n_slices_tmp) /= desc%group_size) CYCLE
296 :
297 : ! we see if there has been a input constraint
298 : ! i.e. if the layout is not -1 we need to fullfil it
299 542754 : IF (.NOT. ALL(PACK(n_slices_tmp == input_settings%distribution_layout, &
300 : [-1, -1, -1] /= input_settings%distribution_layout) &
301 67812 : )) CYCLE
302 :
303 : ! We can not work with a grid that has more local than global grid points.
304 : ! This can happen when a halo region wraps around and overlaps with the other halo.
305 67586 : overlap = .FALSE.
306 270344 : DO dir = 1, 3
307 270344 : IF (n_slices_tmp(dir) > 1) THEN
308 202758 : DO l = 0, n_slices_tmp(dir) - 1
309 135172 : lb = get_limit(desc%npts(dir), n_slices_tmp(dir), l)
310 202758 : IF (lb(2) - lb(1) + 1 + 2*nmin > desc%npts(dir)) overlap = .TRUE.
311 : END DO
312 : END IF
313 : END DO
314 67586 : IF (overlap) CYCLE
315 :
316 : ! a heuristic optimisation to reduce the memory usage
317 : ! we go for the smallest local to real volume
318 : ! volume of the box without the wings / volume of the box with the wings
319 : ! with prefactodesc to promote less cuts in Z dimension
320 : ratio = PRODUCT(REAL(desc%npts, KIND=dp)/n_slices_tmp)/ &
321 : PRODUCT(REAL(desc%npts, KIND=dp)/n_slices_tmp + &
322 71414 : MERGE([0.0_dp, 0.0_dp, 0.0_dp], 2*[1.06_dp*nmin, 1.05_dp*nmin, 1.03_dp*nmin], n_slices_tmp == [1, 1, 1]))
323 55410 : IF (ratio > ratio_best) THEN
324 10172 : ratio_best = ratio
325 10172 : n_slices = n_slices_tmp
326 : END IF
327 :
328 : END DO
329 : END DO
330 :
331 : ! if automatic we can still decide this is a replicated grid
332 : ! if the memory gain (or the gain is messages) is too small.
333 22604 : IF (input_settings%distribution_type == rsgrid_automatic) THEN
334 88216 : volume = PRODUCT(REAL(desc%npts, KIND=dp))
335 : volume_dist = PRODUCT(REAL(desc%npts, KIND=dp)/n_slices + &
336 88216 : MERGE([0, 0, 0], 2*[nmin, nmin, nmin], n_slices == [1, 1, 1]))
337 22054 : IF (volume < volume_dist*input_settings%memory_factor) THEN
338 88216 : n_slices = 1
339 : END IF
340 : END IF
341 :
342 : END IF
343 :
344 148408 : desc%group_dim(:) = n_slices(:)
345 37102 : CALL desc%group%from_dup(pw_grid%para%group)
346 37102 : desc%group_size = desc%group%num_pe
347 37102 : desc%my_pos = desc%group%mepos
348 :
349 148230 : IF (ALL(n_slices == 1)) THEN
350 : ! CASE 1 : only one slice: we do not need overlapping regions and special
351 : ! recombination of the total density
352 36944 : desc%border = border_size
353 36944 : IF (border_size == 0) THEN
354 147776 : desc%perd = 1
355 : ELSE
356 0 : desc%perd = 0
357 : END IF
358 36944 : desc%distributed = .FALSE.
359 36944 : desc%parallel = .TRUE.
360 147776 : desc%group_coor(:) = 0
361 36944 : desc%my_virtual_pos = 0
362 :
363 110832 : ALLOCATE (desc%virtual2real(0:desc%group_size - 1))
364 110832 : ALLOCATE (desc%real2virtual(0:desc%group_size - 1))
365 : ! Start with no reordering
366 110832 : DO i = 0, desc%group_size - 1
367 73888 : desc%virtual2real(i) = i
368 110832 : desc%real2virtual(i) = i
369 : END DO
370 : ELSE
371 : ! CASE 2 : general case
372 : ! periodicity is no longer enforced arbritary directions
373 158 : IF (border_size == 0) THEN
374 120 : desc%perd = 1
375 120 : DO dir = 1, 3
376 120 : IF (n_slices(dir) > 1) desc%perd(dir) = 0
377 : END DO
378 : ELSE
379 512 : desc%perd(:) = 0
380 : END IF
381 : ! we keep a border of nmin points
382 158 : desc%border = nmin
383 : ! we are going parallel on the real space grid
384 158 : desc%parallel = .TRUE.
385 158 : desc%distributed = .TRUE.
386 :
387 : ! set up global info about the distribution
388 474 : ALLOCATE (desc%rank2coord(3, 0:desc%group_size - 1))
389 790 : ALLOCATE (desc%coord2rank(0:desc%group_dim(1) - 1, 0:desc%group_dim(2) - 1, 0:desc%group_dim(3) - 1))
390 474 : ALLOCATE (desc%lb_global(3, 0:desc%group_size - 1))
391 474 : ALLOCATE (desc%ub_global(3, 0:desc%group_size - 1))
392 474 : ALLOCATE (desc%x2coord(desc%lb(1):desc%ub(1)))
393 474 : ALLOCATE (desc%y2coord(desc%lb(2):desc%ub(2)))
394 474 : ALLOCATE (desc%z2coord(desc%lb(3):desc%ub(3)))
395 :
396 474 : DO i = 0, desc%group_size - 1
397 : ! Calculate coordinates in a row-major order (to be SMP-friendly)
398 316 : desc%rank2coord(1, i) = i/(desc%group_dim(2)*desc%group_dim(3))
399 : desc%rank2coord(2, i) = MODULO(i, desc%group_dim(2)*desc%group_dim(3)) &
400 316 : /desc%group_dim(3)
401 316 : desc%rank2coord(3, i) = MODULO(i, desc%group_dim(3))
402 :
403 316 : IF (i == desc%my_pos) THEN
404 632 : desc%group_coor = desc%rank2coord(:, i)
405 : END IF
406 :
407 316 : desc%coord2rank(desc%rank2coord(1, i), desc%rank2coord(2, i), desc%rank2coord(3, i)) = i
408 : ! the lb_global and ub_global correspond to lb_real and ub_real of each task
409 1264 : desc%lb_global(:, i) = desc%lb
410 1264 : desc%ub_global(:, i) = desc%ub
411 1422 : DO dir = 1, 3
412 1264 : IF (desc%group_dim(dir) > 1) THEN
413 316 : lb = get_limit(desc%npts(dir), desc%group_dim(dir), desc%rank2coord(dir, i))
414 316 : desc%lb_global(dir, i) = lb(1) + desc%lb(dir) - 1
415 316 : desc%ub_global(dir, i) = lb(2) + desc%lb(dir) - 1
416 : END IF
417 : END DO
418 : END DO
419 :
420 : ! map a grid point to a CPU coord
421 632 : DO dir = 1, 3
422 1264 : DO l = 0, desc%group_dim(dir) - 1
423 632 : IF (desc%group_dim(dir) > 1) THEN
424 316 : lb = get_limit(desc%npts(dir), desc%group_dim(dir), l)
425 948 : lb = lb + desc%lb(dir) - 1
426 : ELSE
427 316 : lb(1) = desc%lb(dir)
428 316 : lb(2) = desc%ub(dir)
429 : END IF
430 474 : SELECT CASE (dir)
431 : CASE (1)
432 11696 : desc%x2coord(lb(1):lb(2)) = l
433 : CASE (2)
434 12104 : desc%y2coord(lb(1):lb(2)) = l
435 : CASE (3)
436 12200 : desc%z2coord(lb(1):lb(2)) = l
437 : END SELECT
438 : END DO
439 : END DO
440 :
441 : ! an upper bound for the number of neighbours the border is overlapping with
442 632 : DO dir = 1, 3
443 474 : desc%neighbours(dir) = 0
444 632 : IF ((n_slices(dir) > 1) .OR. (border_size > 0)) THEN
445 414 : min_npts_real = HUGE(0)
446 986 : DO l = 0, n_slices(dir) - 1
447 572 : lb = get_limit(desc%npts(dir), n_slices(dir), l)
448 986 : min_npts_real = MIN(lb(2) - lb(1) + 1, min_npts_real)
449 : END DO
450 414 : desc%neighbours(dir) = (desc%border + min_npts_real - 1)/min_npts_real
451 : END IF
452 : END DO
453 :
454 474 : ALLOCATE (desc%virtual2real(0:desc%group_size - 1))
455 474 : ALLOCATE (desc%real2virtual(0:desc%group_size - 1))
456 : ! Start with no reordering
457 474 : DO i = 0, desc%group_size - 1
458 316 : desc%virtual2real(i) = i
459 474 : desc%real2virtual(i) = i
460 : END DO
461 :
462 158 : desc%my_virtual_pos = desc%real2virtual(desc%my_pos)
463 632 : desc%virtual_group_coor(:) = desc%rank2coord(:, desc%my_virtual_pos)
464 :
465 : END IF
466 : END IF
467 :
468 41330 : CALL timestop(handle)
469 :
470 41330 : END SUBROUTINE rs_grid_create_descriptor
471 :
472 : ! **************************************************************************************************
473 : !> \brief ...
474 : !> \param rs ...
475 : !> \param desc ...
476 : ! **************************************************************************************************
477 6066438 : SUBROUTINE rs_grid_create(rs, desc)
478 : TYPE(realspace_grid_type), INTENT(OUT) :: rs
479 : TYPE(realspace_grid_desc_type), INTENT(INOUT), &
480 : TARGET :: desc
481 :
482 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rs_grid_create'
483 :
484 : INTEGER :: handle
485 :
486 288878 : CALL timeset(routineN, handle)
487 :
488 288878 : rs%desc => desc
489 288878 : CALL rs_grid_retain_descriptor(rs%desc)
490 :
491 288878 : IF (desc%pw%para%mode == PW_MODE_LOCAL) THEN
492 : ! The corresponding group has dimension 1
493 : ! All operations will be done locally
494 70976 : rs%lb_real = desc%lb
495 70976 : rs%ub_real = desc%ub
496 70976 : rs%lb_local = rs%lb_real - desc%border*(1 - desc%perd)
497 70976 : rs%ub_local = rs%ub_real + desc%border*(1 - desc%perd)
498 70976 : rs%npts_local = rs%ub_local - rs%lb_local + 1
499 70976 : rs%ngpts_local = PRODUCT(rs%npts_local)
500 : END IF
501 :
502 1153116 : IF (ALL(rs%desc%group_dim == 1)) THEN
503 : ! CASE 1 : only one slice: we do not need overlapping regions and special
504 : ! recombination of the total density
505 1148408 : rs%lb_real = desc%lb
506 1148408 : rs%ub_real = desc%ub
507 1148408 : rs%lb_local = rs%lb_real - desc%border*(1 - desc%perd)
508 1148408 : rs%ub_local = rs%ub_real + desc%border*(1 - desc%perd)
509 1148408 : rs%npts_local = rs%ub_local - rs%lb_local + 1
510 1148408 : rs%ngpts_local = PRODUCT(rs%npts_local)
511 : ELSE
512 : ! CASE 2 : general case
513 : ! extract some more derived quantities about the local grid
514 7104 : rs%lb_real = desc%lb_global(:, desc%my_virtual_pos)
515 7104 : rs%ub_real = desc%ub_global(:, desc%my_virtual_pos)
516 7104 : rs%lb_local = rs%lb_real - desc%border*(1 - desc%perd)
517 7104 : rs%ub_local = rs%ub_real + desc%border*(1 - desc%perd)
518 7104 : rs%npts_local = rs%ub_local - rs%lb_local + 1
519 7104 : rs%ngpts_local = PRODUCT(rs%npts_local)
520 : END IF
521 :
522 288878 : CALL offload_create_buffer(rs%ngpts_local, rs%buffer)
523 : rs%r(rs%lb_local(1):rs%ub_local(1), &
524 : rs%lb_local(2):rs%ub_local(2), &
525 288878 : rs%lb_local(3):rs%ub_local(3)) => rs%buffer%host_buffer
526 :
527 866634 : ALLOCATE (rs%px(desc%npts(1)))
528 866634 : ALLOCATE (rs%py(desc%npts(2)))
529 866634 : ALLOCATE (rs%pz(desc%npts(3)))
530 :
531 288878 : CALL timestop(handle)
532 :
533 288878 : END SUBROUTINE rs_grid_create
534 :
535 : ! **************************************************************************************************
536 : !> \brief Defines a new ordering of ranks on this realspace grid, recalculating
537 : !> the data bounds and reallocating the grid. As a result, each MPI process
538 : !> now has a real rank (i.e., its rank in the MPI communicator from the pw grid)
539 : !> and a virtual rank (the rank of the process where the data now owned by this
540 : !> process would reside in an ordinary cartesian distribution).
541 : !> NB. Since the grid size required may change, the caller should be sure to release
542 : !> and recreate the corresponding rs_grids
543 : !> The desc%real2virtual and desc%virtual2real arrays can be used to map
544 : !> a physical rank to the 'rank' of data owned by that process and vice versa
545 : !> \param desc ...
546 : !> \param real2virtual ...
547 : !> \par History
548 : !> 04-2009 created [Iain Bethune]
549 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
550 : ! **************************************************************************************************
551 6 : PURE SUBROUTINE rs_grid_reorder_ranks(desc, real2virtual)
552 :
553 : TYPE(realspace_grid_desc_type), INTENT(INOUT) :: desc
554 : INTEGER, DIMENSION(:), INTENT(IN) :: real2virtual
555 :
556 : INTEGER :: i
557 :
558 18 : desc%real2virtual(:) = real2virtual
559 :
560 18 : DO i = 0, desc%group_size - 1
561 18 : desc%virtual2real(desc%real2virtual(i)) = i
562 : END DO
563 :
564 6 : desc%my_virtual_pos = desc%real2virtual(desc%my_pos)
565 :
566 12 : IF (.NOT. ALL(desc%group_dim == 1)) THEN
567 24 : desc%virtual_group_coor(:) = desc%rank2coord(:, desc%my_virtual_pos)
568 : END IF
569 :
570 6 : END SUBROUTINE rs_grid_reorder_ranks
571 :
572 : ! **************************************************************************************************
573 : !> \brief Print information on grids to output
574 : !> \param rs ...
575 : !> \param iounit ...
576 : !> \author JGH (17-May-2007)
577 : ! **************************************************************************************************
578 17512 : SUBROUTINE rs_grid_print(rs, iounit)
579 : TYPE(realspace_grid_type), INTENT(IN) :: rs
580 : INTEGER, INTENT(in) :: iounit
581 :
582 : INTEGER :: dir, i, nn
583 : REAL(KIND=dp) :: pp(3)
584 :
585 17512 : IF (rs%desc%parallel) THEN
586 17220 : IF (iounit > 0) THEN
587 : WRITE (iounit, '(/,A,T71,I10)') &
588 7767 : " RS_GRID| Information for grid number ", rs%desc%pw%id_nr
589 31068 : DO i = 1, 3
590 23301 : WRITE (iounit, '(A,I3,T30,2I8,T62,A,T71,I10)') " RS_GRID| Bounds ", &
591 54369 : i, rs%desc%lb(i), rs%desc%ub(i), "Points:", rs%desc%npts(i)
592 : END DO
593 7767 : IF (.NOT. rs%desc%distributed) THEN
594 7752 : WRITE (iounit, '(A)') " RS_GRID| Real space fully replicated"
595 : WRITE (iounit, '(A,T71,I10)') &
596 7752 : " RS_GRID| Group size ", rs%desc%group_dim(2)
597 : ELSE
598 60 : DO dir = 1, 3
599 60 : IF (rs%desc%perd(dir) /= 1) THEN
600 : WRITE (iounit, '(A,T71,I3,A)') &
601 15 : " RS_GRID| Real space distribution over ", rs%desc%group_dim(dir), " groups"
602 : WRITE (iounit, '(A,T71,I10)') &
603 15 : " RS_GRID| Real space distribution along direction ", dir
604 : WRITE (iounit, '(A,T71,I10)') &
605 15 : " RS_GRID| Border size ", rs%desc%border
606 : END IF
607 : END DO
608 : END IF
609 : END IF
610 17220 : IF (rs%desc%distributed) THEN
611 120 : DO dir = 1, 3
612 120 : IF (rs%desc%perd(dir) /= 1) THEN
613 30 : nn = rs%npts_local(dir)
614 30 : CALL rs%desc%group%sum(nn)
615 120 : pp(1) = REAL(nn, KIND=dp)/REAL(PRODUCT(rs%desc%group_dim), KIND=dp)
616 30 : nn = rs%npts_local(dir)
617 30 : CALL rs%desc%group%max(nn)
618 30 : pp(2) = REAL(nn, KIND=dp)
619 30 : nn = rs%npts_local(dir)
620 30 : CALL rs%desc%group%min(nn)
621 30 : pp(3) = REAL(nn, KIND=dp)
622 30 : IF (iounit > 0) THEN
623 15 : WRITE (iounit, '(A,T48,A)') " RS_GRID| Distribution", &
624 30 : " Average Max Min"
625 15 : WRITE (iounit, '(A,T45,F12.1,2I12)') " RS_GRID| Planes ", &
626 30 : pp(1), NINT(pp(2)), NINT(pp(3))
627 : END IF
628 : END IF
629 : END DO
630 : ! WRITE ( iounit, '(/)' )
631 : END IF
632 : ELSE
633 292 : IF (iounit > 0) THEN
634 : WRITE (iounit, '(/,A,T71,I10)') &
635 180 : " RS_GRID| Information for grid number ", rs%desc%pw%id_nr
636 720 : DO i = 1, 3
637 540 : WRITE (iounit, '(A,I3,T30,2I8,T62,A,T71,I10)') " RS_GRID| Bounds ", &
638 1260 : i, rs%desc%lb(i), rs%desc%ub(i), "Points:", rs%desc%npts(i)
639 : END DO
640 : ! WRITE ( iounit, '(/)' )
641 : END IF
642 : END IF
643 :
644 17512 : END SUBROUTINE rs_grid_print
645 :
646 : ! **************************************************************************************************
647 : !> \brief ...
648 : !> \param rs ...
649 : !> \param pw ...
650 : ! **************************************************************************************************
651 1368160 : SUBROUTINE transfer_rs2pw(rs, pw)
652 : TYPE(realspace_grid_type), INTENT(IN) :: rs
653 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw
654 :
655 : CHARACTER(len=*), PARAMETER :: routineN = 'transfer_rs2pw'
656 :
657 : INTEGER :: handle, handle2, i
658 :
659 1368160 : CALL timeset(routineN, handle2)
660 1368160 : CALL timeset(routineN//"_"//TRIM(ADJUSTL(cp_to_string(CEILING(pw%pw_grid%cutoff/10)*10))), handle)
661 :
662 1368160 : IF (.NOT. ASSOCIATED(rs%desc%pw, pw%pw_grid)) THEN
663 0 : CPABORT("Different rs and pw indentifiers")
664 : END IF
665 :
666 1368160 : IF (rs%desc%distributed) THEN
667 1852 : CALL transfer_rs2pw_distributed(rs, pw)
668 1366308 : ELSE IF (rs%desc%parallel) THEN
669 1150150 : CALL transfer_rs2pw_replicated(rs, pw)
670 : ELSE ! treat simple serial case locally
671 216158 : IF (rs%desc%border == 0) THEN
672 864632 : CALL dcopy(SIZE(rs%r), rs%r, 1, pw%array, 1)
673 : ELSE
674 0 : CPASSERT(LBOUND(pw%array, 3) == rs%lb_real(3))
675 0 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(pw,rs)
676 : DO i = rs%lb_real(3), rs%ub_real(3)
677 : pw%array(:, :, i) = rs%r(rs%lb_real(1):rs%ub_real(1), &
678 : rs%lb_real(2):rs%ub_real(2), i)
679 : END DO
680 : !$OMP END PARALLEL DO
681 : END IF
682 : END IF
683 :
684 1368160 : CALL timestop(handle)
685 1368160 : CALL timestop(handle2)
686 :
687 1368160 : END SUBROUTINE transfer_rs2pw
688 :
689 : ! **************************************************************************************************
690 : !> \brief ...
691 : !> \param rs ...
692 : !> \param pw ...
693 : ! **************************************************************************************************
694 1293163 : SUBROUTINE transfer_pw2rs(rs, pw)
695 :
696 : TYPE(realspace_grid_type), INTENT(IN) :: rs
697 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
698 :
699 : CHARACTER(len=*), PARAMETER :: routineN = 'transfer_pw2rs'
700 :
701 : INTEGER :: handle, handle2, i, im, j, jm, k, km
702 :
703 1293163 : CALL timeset(routineN, handle2)
704 1293163 : CALL timeset(routineN//"_"//TRIM(ADJUSTL(cp_to_string(CEILING(pw%pw_grid%cutoff/10)*10))), handle)
705 :
706 1293163 : IF (.NOT. ASSOCIATED(rs%desc%pw, pw%pw_grid)) THEN
707 0 : CPABORT("Different rs and pw indentifiers")
708 : END IF
709 :
710 1293163 : IF (rs%desc%distributed) THEN
711 872 : CALL transfer_pw2rs_distributed(rs, pw)
712 1292291 : ELSE IF (rs%desc%parallel) THEN
713 1039516 : CALL transfer_pw2rs_replicated(rs, pw)
714 : ELSE ! treat simple serial case locally
715 252775 : IF (rs%desc%border == 0) THEN
716 1011100 : CALL dcopy(SIZE(rs%r), pw%array, 1, rs%r, 1)
717 : ELSE
718 : !$OMP PARALLEL DO DEFAULT(NONE) &
719 : !$OMP PRIVATE(i,im,j,jm,k,km) &
720 0 : !$OMP SHARED(pw,rs)
721 : DO k = rs%lb_local(3), rs%ub_local(3)
722 : IF (k < rs%lb_real(3)) THEN
723 : km = k + rs%desc%npts(3)
724 : ELSE IF (k > rs%ub_real(3)) THEN
725 : km = k - rs%desc%npts(3)
726 : ELSE
727 : km = k
728 : END IF
729 : DO j = rs%lb_local(2), rs%ub_local(2)
730 : IF (j < rs%lb_real(2)) THEN
731 : jm = j + rs%desc%npts(2)
732 : ELSE IF (j > rs%ub_real(2)) THEN
733 : jm = j - rs%desc%npts(2)
734 : ELSE
735 : jm = j
736 : END IF
737 : DO i = rs%lb_local(1), rs%ub_local(1)
738 : IF (i < rs%lb_real(1)) THEN
739 : im = i + rs%desc%npts(1)
740 : ELSE IF (i > rs%ub_real(1)) THEN
741 : im = i - rs%desc%npts(1)
742 : ELSE
743 : im = i
744 : END IF
745 : rs%r(i, j, k) = pw%array(im, jm, km)
746 : END DO
747 : END DO
748 : END DO
749 : !$OMP END PARALLEL DO
750 : END IF
751 : END IF
752 :
753 1293163 : CALL timestop(handle)
754 1293163 : CALL timestop(handle2)
755 :
756 1293163 : END SUBROUTINE transfer_pw2rs
757 :
758 : ! **************************************************************************************************
759 : !> \brief transfer from a realspace grid to a planewave grid
760 : !> \param rs ...
761 : !> \param pw ...
762 : ! **************************************************************************************************
763 1150150 : SUBROUTINE transfer_rs2pw_replicated(rs, pw)
764 : TYPE(realspace_grid_type), INTENT(IN) :: rs
765 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw
766 :
767 : INTEGER :: dest, ii, ip, ix, iy, iz, nma, nn, s(3), &
768 : source
769 1150150 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount
770 : INTEGER, DIMENSION(3) :: lb, ub
771 1150150 : INTEGER, DIMENSION(:, :, :, :), POINTER :: bounds_by_rank
772 1150150 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: recvbuf, sendbuf, swaparray
773 :
774 : ! Avoid nested component array sections: GCC 16 can generate invalid bounds checks.
775 1150150 : bounds_by_rank => pw%pw_grid%para%bo
776 : ASSOCIATE (np => pw%pw_grid%para%group%num_pe, bo => bounds_by_rank(:, :, :, 1), &
777 : pbo => pw%pw_grid%bounds, group => pw%pw_grid%para%group, mepos => pw%pw_grid%para%group%mepos, &
778 : grid => rs%r)
779 3450450 : ALLOCATE (rcount(0:np - 1))
780 3450450 : DO ip = 1, np
781 10351350 : rcount(ip - 1) = PRODUCT(bo(2, :, ip) - bo(1, :, ip) + 1)
782 : END DO
783 3450450 : nma = MAXVAL(rcount(0:np - 1))
784 4600600 : ALLOCATE (sendbuf(nma), recvbuf(nma))
785 43424367600 : sendbuf = 1.0E99_dp; recvbuf = 1.0E99_dp ! init mpi'ed buffers to silence warnings under valgrind
786 :
787 : !sample peak memory
788 1150150 : CALL m_memory()
789 :
790 1150150 : dest = MODULO(mepos + 1, np)
791 1150150 : source = MODULO(mepos - 1, np)
792 1150150 : sendbuf = 0.0_dp
793 :
794 2300300 : DO ip = 1, np
795 :
796 9201200 : lb = pbo(1, :) + bo(1, :, MODULO(mepos - ip, np) + 1) - 1
797 9201200 : ub = pbo(1, :) + bo(2, :, MODULO(mepos - ip, np) + 1) - 1
798 : ! this loop takes about the same time as the message passing call
799 : ! notice that the range of ix is only a small fraction of the first index of grid
800 : ! therefore it seems faster to have the second index as the innermost loop
801 : ! if this runs on many cpus
802 : ! tested on itanium, pentium4, opteron, ultrasparc...
803 9201200 : s = ub - lb + 1
804 57770480 : DO iz = lb(3), ub(3)
805 1021915492 : DO ix = lb(1), ub(1)
806 964145012 : ii = (iz - lb(3))*s(1)*s(2) + (ix - lb(1)) + 1
807 44139472380 : DO iy = lb(2), ub(2)
808 43119857188 : sendbuf(ii) = sendbuf(ii) + grid(ix, iy, iz)
809 44084002200 : ii = ii + s(1)
810 : END DO
811 : END DO
812 : END DO
813 2300300 : IF (ip == np) EXIT
814 1150150 : CALL group%sendrecv(sendbuf, dest, recvbuf, source, 13)
815 1150150 : CALL MOVE_ALLOC(sendbuf, swaparray)
816 1150150 : CALL MOVE_ALLOC(recvbuf, sendbuf)
817 2300300 : CALL MOVE_ALLOC(swaparray, recvbuf)
818 : END DO
819 1150150 : nn = rcount(mepos)
820 : END ASSOCIATE
821 :
822 1150150 : CALL dcopy(nn, sendbuf, 1, pw%array, 1)
823 :
824 1150150 : DEALLOCATE (rcount)
825 1150150 : DEALLOCATE (sendbuf)
826 1150150 : DEALLOCATE (recvbuf)
827 :
828 1150150 : END SUBROUTINE transfer_rs2pw_replicated
829 :
830 : ! **************************************************************************************************
831 : !> \brief transfer from a planewave grid to a realspace grid
832 : !> \param rs ...
833 : !> \param pw ...
834 : ! **************************************************************************************************
835 1039516 : SUBROUTINE transfer_pw2rs_replicated(rs, pw)
836 : TYPE(realspace_grid_type), INTENT(IN) :: rs
837 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
838 :
839 : INTEGER :: dest, i, ii, im, ip, ix, iy, iz, j, jm, &
840 : k, km, nma, nn, source
841 1039516 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount
842 : INTEGER, DIMENSION(3) :: lb, ub
843 1039516 : INTEGER, DIMENSION(:, :, :, :), POINTER :: bounds_by_rank
844 1039516 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: recvbuf, sendbuf, swaparray
845 3118548 : TYPE(mp_request_type), DIMENSION(2) :: req
846 :
847 : ! Keep the reverse transfer on the same compiler-safe array-section path.
848 1039516 : bounds_by_rank => pw%pw_grid%para%bo
849 : ASSOCIATE (np => pw%pw_grid%para%group%num_pe, bo => bounds_by_rank(:, :, :, 1), &
850 : pbo => pw%pw_grid%bounds, group => pw%pw_grid%para%group, mepos => pw%pw_grid%para%group%mepos, &
851 : grid => rs%r)
852 3118548 : ALLOCATE (rcount(0:np - 1))
853 3118548 : DO ip = 1, np
854 9355644 : rcount(ip - 1) = PRODUCT(bo(2, :, ip) - bo(1, :, ip) + 1)
855 : END DO
856 3118548 : nma = MAXVAL(rcount(0:np - 1))
857 4158064 : ALLOCATE (sendbuf(nma), recvbuf(nma))
858 39984648968 : sendbuf = 1.0E99_dp; recvbuf = 1.0E99_dp ! init mpi'ed buffers to silence warnings under valgrind
859 :
860 : !sample peak memory
861 1039516 : CALL m_memory()
862 :
863 1039516 : nn = rcount(mepos)
864 1039516 : CALL dcopy(nn, pw%array, 1, sendbuf, 1)
865 :
866 1039516 : dest = MODULO(mepos + 1, np)
867 1039516 : source = MODULO(mepos - 1, np)
868 :
869 3118548 : DO ip = 0, np - 1
870 : ! we must shift the buffer only np-1 times around
871 2079032 : IF (ip /= np - 1) THEN
872 : CALL group%isendrecv(sendbuf, dest, recvbuf, source, &
873 1039516 : req(1), req(2), 13)
874 : END IF
875 8316128 : lb = pbo(1, :) + bo(1, :, MODULO(mepos - ip, np) + 1) - 1
876 8316128 : ub = pbo(1, :) + bo(2, :, MODULO(mepos - ip, np) + 1) - 1
877 2079032 : ii = 0
878 : ! this loop takes about the same time as the message passing call
879 : ! If I read the code correctly then:
880 52822572 : DO iz = lb(3), ub(3)
881 1813325704 : DO iy = lb(2), ub(2)
882 41519440252 : DO ix = lb(1), ub(1)
883 39708193580 : ii = ii + 1
884 41468696712 : grid(ix, iy, iz) = sendbuf(ii)
885 : END DO
886 : END DO
887 : END DO
888 2079032 : IF (ip /= np - 1) THEN
889 1039516 : CALL mp_waitall(req)
890 : END IF
891 2079032 : CALL MOVE_ALLOC(sendbuf, swaparray)
892 2079032 : CALL MOVE_ALLOC(recvbuf, sendbuf)
893 3118548 : CALL MOVE_ALLOC(swaparray, recvbuf)
894 : END DO
895 2079032 : IF (rs%desc%border > 0) THEN
896 : !$OMP PARALLEL DO DEFAULT(NONE) &
897 : !$OMP PRIVATE(i,im,j,jm,k,km) &
898 0 : !$OMP SHARED(rs)
899 : DO k = rs%lb_local(3), rs%ub_local(3)
900 : IF (k < rs%lb_real(3)) THEN
901 : km = k + rs%desc%npts(3)
902 : ELSE IF (k > rs%ub_real(3)) THEN
903 : km = k - rs%desc%npts(3)
904 : ELSE
905 : km = k
906 : END IF
907 : DO j = rs%lb_local(2), rs%ub_local(2)
908 : IF (j < rs%lb_real(2)) THEN
909 : jm = j + rs%desc%npts(2)
910 : ELSE IF (j > rs%ub_real(2)) THEN
911 : jm = j - rs%desc%npts(2)
912 : ELSE
913 : jm = j
914 : END IF
915 : DO i = rs%lb_local(1), rs%ub_local(1)
916 : IF (i < rs%lb_real(1)) THEN
917 : im = i + rs%desc%npts(1)
918 : ELSE IF (i > rs%ub_real(1)) THEN
919 : im = i - rs%desc%npts(1)
920 : ELSE
921 : im = i
922 : END IF
923 : rs%r(i, j, k) = rs%r(im, jm, km)
924 : END DO
925 : END DO
926 : END DO
927 : !$OMP END PARALLEL DO
928 : END IF
929 : END ASSOCIATE
930 :
931 1039516 : DEALLOCATE (rcount)
932 1039516 : DEALLOCATE (sendbuf)
933 1039516 : DEALLOCATE (recvbuf)
934 :
935 1039516 : END SUBROUTINE transfer_pw2rs_replicated
936 :
937 : ! **************************************************************************************************
938 : !> \brief does the rs2pw transfer in the case where the rs grid is
939 : !> distributed (3D domain decomposition)
940 : !> \param rs ...
941 : !> \param pw ...
942 : !> \par History
943 : !> 12.2007 created [Matt Watkins]
944 : !> 9.2008 reduced amount of halo data sent [Iain Bethune]
945 : !> 10.2008 added non-blocking communication [Iain Bethune]
946 : !> 4.2009 added support for rank-reordering on the grid [Iain Bethune]
947 : !> 12.2009 added OMP and sparse alltoall [Iain Bethune]
948 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2008-2009 on behalf of the HECToR project
949 : !> \note
950 : !> the transfer is a two step procedure. For example, for the rs2pw transfer:
951 : !>
952 : !> 1) Halo-exchange in 3D so that the local part of the rs_grid contains the full data
953 : !> 2) an alltoall communication to redistribute the local rs_grid to the local pw_grid
954 : !>
955 : !> the halo exchange is most expensive on a large number of CPUs. Particular in this halo
956 : !> exchange is that the border region is rather large (e.g. 20 points) and that it might overlap
957 : !> with the central domain of several CPUs (i.e. next nearest neighbors)
958 : ! **************************************************************************************************
959 1852 : SUBROUTINE transfer_rs2pw_distributed(rs, pw)
960 : TYPE(realspace_grid_type), INTENT(IN) :: rs
961 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
962 :
963 : CHARACTER(LEN=200) :: error_string
964 : INTEGER :: completed, dest_down, dest_up, i, idir, j, k, lb, my_id, my_pw_rank, my_rs_rank, &
965 : n_shifts, nn, num_threads, position, source_down, source_up, ub, x, y, z
966 1852 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dshifts, recv_disps, recv_sizes, &
967 1852 : send_disps, send_sizes, ushifts
968 3704 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: bounds, recv_tasks, send_tasks
969 : INTEGER, DIMENSION(2) :: neighbours, pos
970 : INTEGER, DIMENSION(3) :: coords, lb_recv, lb_recv_down, lb_recv_up, lb_send, lb_send_down, &
971 : lb_send_up, ub_recv, ub_recv_down, ub_recv_up, ub_send, ub_send_down, ub_send_up
972 : LOGICAL, DIMENSION(3) :: halo_swapped
973 : REAL(KIND=dp) :: pw_sum, rs_sum
974 1852 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: recv_buf_3d_down, recv_buf_3d_up, &
975 1852 : send_buf_3d_down, send_buf_3d_up
976 3704 : TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: recv_bufs, send_bufs
977 1852 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_reqs, send_reqs
978 9260 : TYPE(mp_request_type), DIMENSION(4) :: req
979 :
980 1852 : num_threads = 1
981 1852 : my_id = 0
982 :
983 : ! safety check, to be removed once we're absolute sure the routine is correct
984 : IF (debug_this_module) THEN
985 : rs_sum = accurate_sum(rs%r)*ABS(det_3x3(rs%desc%dh))
986 : CALL rs%desc%group%sum(rs_sum)
987 : END IF
988 :
989 1852 : halo_swapped = .FALSE.
990 : ! We don't need to send the 'edges' of the halos that have already been sent
991 : ! Halos are contiguous in memory in z-direction only, so swap these first,
992 : ! and send less data in the y and x directions which are more expensive
993 :
994 7408 : DO idir = 3, 1, -1
995 :
996 5556 : IF (rs%desc%perd(idir) /= 1) THEN
997 :
998 14412 : ALLOCATE (dshifts(0:rs%desc%neighbours(idir)))
999 9608 : ALLOCATE (ushifts(0:rs%desc%neighbours(idir)))
1000 :
1001 4804 : ushifts = 0
1002 4804 : dshifts = 0
1003 :
1004 : ! check that we don't try to send data to ourself
1005 6656 : DO n_shifts = 1, MIN(rs%desc%neighbours(idir), rs%desc%group_dim(idir) - 1)
1006 :
1007 : ! need to take into account the possible varying widths of neighbouring cells
1008 : ! offset_up and offset_down hold the real size of the neighbouring cells
1009 1852 : position = MODULO(rs%desc%virtual_group_coor(idir) - n_shifts, rs%desc%group_dim(idir))
1010 1852 : neighbours = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), position)
1011 1852 : dshifts(n_shifts) = dshifts(n_shifts - 1) + (neighbours(2) - neighbours(1) + 1)
1012 :
1013 1852 : position = MODULO(rs%desc%virtual_group_coor(idir) + n_shifts, rs%desc%group_dim(idir))
1014 1852 : neighbours = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), position)
1015 1852 : ushifts(n_shifts) = ushifts(n_shifts - 1) + (neighbours(2) - neighbours(1) + 1)
1016 :
1017 : ! The border data has to be send/received from the neighbours
1018 : ! First we calculate the source and destination processes for the shift
1019 : ! We do both shifts at once to allow for more overlap of communication and buffer packing/unpacking
1020 :
1021 1852 : CALL cart_shift(rs, idir, -1*n_shifts, source_down, dest_down)
1022 :
1023 7408 : lb_send_down(:) = rs%lb_local(:)
1024 7408 : lb_recv_down(:) = rs%lb_local(:)
1025 7408 : ub_recv_down(:) = rs%ub_local(:)
1026 7408 : ub_send_down(:) = rs%ub_local(:)
1027 :
1028 1852 : IF (dshifts(n_shifts - 1) <= rs%desc%border) THEN
1029 1852 : ub_send_down(idir) = lb_send_down(idir) + rs%desc%border - 1 - dshifts(n_shifts - 1)
1030 : lb_send_down(idir) = MAX(lb_send_down(idir), &
1031 1852 : lb_send_down(idir) + rs%desc%border - dshifts(n_shifts))
1032 :
1033 1852 : ub_recv_down(idir) = ub_recv_down(idir) - rs%desc%border
1034 : lb_recv_down(idir) = MAX(lb_recv_down(idir) + rs%desc%border, &
1035 1852 : ub_recv_down(idir) - rs%desc%border + 1 + ushifts(n_shifts - 1))
1036 : ELSE
1037 0 : lb_send_down(idir) = 0
1038 0 : ub_send_down(idir) = -1
1039 0 : lb_recv_down(idir) = 0
1040 0 : ub_recv_down(idir) = -1
1041 : END IF
1042 :
1043 7408 : DO i = 1, 3
1044 7408 : IF (halo_swapped(i)) THEN
1045 554 : lb_send_down(i) = rs%lb_real(i)
1046 554 : ub_send_down(i) = rs%ub_real(i)
1047 554 : lb_recv_down(i) = rs%lb_real(i)
1048 554 : ub_recv_down(i) = rs%ub_real(i)
1049 : END IF
1050 : END DO
1051 :
1052 : ! post the receive
1053 0 : ALLOCATE (recv_buf_3d_down(lb_recv_down(1):ub_recv_down(1), &
1054 9260 : lb_recv_down(2):ub_recv_down(2), lb_recv_down(3):ub_recv_down(3)))
1055 1852 : CALL rs%desc%group%irecv(recv_buf_3d_down, source_down, req(1))
1056 :
1057 : ! now allocate, pack and send the send buffer
1058 7408 : nn = PRODUCT(ub_send_down - lb_send_down + 1)
1059 0 : ALLOCATE (send_buf_3d_down(lb_send_down(1):ub_send_down(1), &
1060 9260 : lb_send_down(2):ub_send_down(2), lb_send_down(3):ub_send_down(3)))
1061 :
1062 : !$OMP PARALLEL DEFAULT(NONE), &
1063 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1064 1852 : !$OMP SHARED(send_buf_3d_down,rs,lb_send_down,ub_send_down)
1065 : !$ num_threads = MIN(omp_get_max_threads(), ub_send_down(3) - lb_send_down(3) + 1)
1066 : !$ my_id = omp_get_thread_num()
1067 : IF (my_id < num_threads) THEN
1068 : lb = lb_send_down(3) + ((ub_send_down(3) - lb_send_down(3) + 1)*my_id)/num_threads
1069 : ub = lb_send_down(3) + ((ub_send_down(3) - lb_send_down(3) + 1)*(my_id + 1))/num_threads - 1
1070 :
1071 : send_buf_3d_down(lb_send_down(1):ub_send_down(1), lb_send_down(2):ub_send_down(2), &
1072 : lb:ub) = rs%r(lb_send_down(1):ub_send_down(1), &
1073 : lb_send_down(2):ub_send_down(2), lb:ub)
1074 : END IF
1075 : !$OMP END PARALLEL
1076 :
1077 1852 : CALL rs%desc%group%isend(send_buf_3d_down, dest_down, req(3))
1078 :
1079 : ! Now for the other direction
1080 1852 : CALL cart_shift(rs, idir, n_shifts, source_up, dest_up)
1081 :
1082 7408 : lb_send_up(:) = rs%lb_local(:)
1083 7408 : lb_recv_up(:) = rs%lb_local(:)
1084 7408 : ub_recv_up(:) = rs%ub_local(:)
1085 7408 : ub_send_up(:) = rs%ub_local(:)
1086 :
1087 1852 : IF (ushifts(n_shifts - 1) <= rs%desc%border) THEN
1088 :
1089 1852 : lb_send_up(idir) = ub_send_up(idir) - rs%desc%border + 1 + ushifts(n_shifts - 1)
1090 : ub_send_up(idir) = MIN(ub_send_up(idir), &
1091 1852 : ub_send_up(idir) - rs%desc%border + ushifts(n_shifts))
1092 :
1093 1852 : lb_recv_up(idir) = lb_recv_up(idir) + rs%desc%border
1094 : ub_recv_up(idir) = MIN(ub_recv_up(idir) - rs%desc%border, &
1095 1852 : lb_recv_up(idir) + rs%desc%border - 1 - dshifts(n_shifts - 1))
1096 : ELSE
1097 0 : lb_send_up(idir) = 0
1098 0 : ub_send_up(idir) = -1
1099 0 : lb_recv_up(idir) = 0
1100 0 : ub_recv_up(idir) = -1
1101 : END IF
1102 :
1103 7408 : DO i = 1, 3
1104 7408 : IF (halo_swapped(i)) THEN
1105 554 : lb_send_up(i) = rs%lb_real(i)
1106 554 : ub_send_up(i) = rs%ub_real(i)
1107 554 : lb_recv_up(i) = rs%lb_real(i)
1108 554 : ub_recv_up(i) = rs%ub_real(i)
1109 : END IF
1110 : END DO
1111 :
1112 : ! post the receive
1113 0 : ALLOCATE (recv_buf_3d_up(lb_recv_up(1):ub_recv_up(1), &
1114 9260 : lb_recv_up(2):ub_recv_up(2), lb_recv_up(3):ub_recv_up(3)))
1115 1852 : CALL rs%desc%group%irecv(recv_buf_3d_up, source_up, req(2))
1116 :
1117 : ! now allocate,pack and send the send buffer
1118 7408 : nn = PRODUCT(ub_send_up - lb_send_up + 1)
1119 0 : ALLOCATE (send_buf_3d_up(lb_send_up(1):ub_send_up(1), &
1120 9260 : lb_send_up(2):ub_send_up(2), lb_send_up(3):ub_send_up(3)))
1121 :
1122 : !$OMP PARALLEL DEFAULT(NONE), &
1123 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1124 1852 : !$OMP SHARED(send_buf_3d_up,rs,lb_send_up,ub_send_up)
1125 : !$ num_threads = MIN(omp_get_max_threads(), ub_send_up(3) - lb_send_up(3) + 1)
1126 : !$ my_id = omp_get_thread_num()
1127 : IF (my_id < num_threads) THEN
1128 : lb = lb_send_up(3) + ((ub_send_up(3) - lb_send_up(3) + 1)*my_id)/num_threads
1129 : ub = lb_send_up(3) + ((ub_send_up(3) - lb_send_up(3) + 1)*(my_id + 1))/num_threads - 1
1130 :
1131 : send_buf_3d_up(lb_send_up(1):ub_send_up(1), lb_send_up(2):ub_send_up(2), &
1132 : lb:ub) = rs%r(lb_send_up(1):ub_send_up(1), &
1133 : lb_send_up(2):ub_send_up(2), lb:ub)
1134 : END IF
1135 : !$OMP END PARALLEL
1136 :
1137 1852 : CALL rs%desc%group%isend(send_buf_3d_up, dest_up, req(4))
1138 :
1139 : ! wait for a recv to complete, then we can unpack
1140 :
1141 5556 : DO i = 1, 2
1142 :
1143 3704 : CALL mp_waitany(req(1:2), completed)
1144 :
1145 5556 : IF (completed == 1) THEN
1146 :
1147 : ! only some procs may need later shifts
1148 1852 : IF (ub_recv_down(idir) >= lb_recv_down(idir)) THEN
1149 : ! Sum the data in the RS Grid
1150 : !$OMP PARALLEL DEFAULT(NONE), &
1151 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1152 1852 : !$OMP SHARED(recv_buf_3d_down,rs,lb_recv_down,ub_recv_down)
1153 : !$ num_threads = MIN(omp_get_max_threads(), ub_recv_down(3) - lb_recv_down(3) + 1)
1154 : !$ my_id = omp_get_thread_num()
1155 : IF (my_id < num_threads) THEN
1156 : lb = lb_recv_down(3) + ((ub_recv_down(3) - lb_recv_down(3) + 1)*my_id)/num_threads
1157 : ub = lb_recv_down(3) + ((ub_recv_down(3) - lb_recv_down(3) + 1)*(my_id + 1))/num_threads - 1
1158 :
1159 : rs%r(lb_recv_down(1):ub_recv_down(1), &
1160 : lb_recv_down(2):ub_recv_down(2), lb:ub) = &
1161 : rs%r(lb_recv_down(1):ub_recv_down(1), &
1162 : lb_recv_down(2):ub_recv_down(2), lb:ub) + &
1163 : recv_buf_3d_down(:, :, lb:ub)
1164 : END IF
1165 : !$OMP END PARALLEL
1166 : END IF
1167 1852 : DEALLOCATE (recv_buf_3d_down)
1168 : ELSE
1169 :
1170 : ! only some procs may need later shifts
1171 1852 : IF (ub_recv_up(idir) >= lb_recv_up(idir)) THEN
1172 : ! Sum the data in the RS Grid
1173 : !$OMP PARALLEL DEFAULT(NONE), &
1174 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1175 1852 : !$OMP SHARED(recv_buf_3d_up,rs,lb_recv_up,ub_recv_up)
1176 : !$ num_threads = MIN(omp_get_max_threads(), ub_recv_up(3) - lb_recv_up(3) + 1)
1177 : !$ my_id = omp_get_thread_num()
1178 : IF (my_id < num_threads) THEN
1179 : lb = lb_recv_up(3) + ((ub_recv_up(3) - lb_recv_up(3) + 1)*my_id)/num_threads
1180 : ub = lb_recv_up(3) + ((ub_recv_up(3) - lb_recv_up(3) + 1)*(my_id + 1))/num_threads - 1
1181 :
1182 : rs%r(lb_recv_up(1):ub_recv_up(1), &
1183 : lb_recv_up(2):ub_recv_up(2), lb:ub) = &
1184 : rs%r(lb_recv_up(1):ub_recv_up(1), &
1185 : lb_recv_up(2):ub_recv_up(2), lb:ub) + &
1186 : recv_buf_3d_up(:, :, lb:ub)
1187 : END IF
1188 : !$OMP END PARALLEL
1189 : END IF
1190 1852 : DEALLOCATE (recv_buf_3d_up)
1191 : END IF
1192 :
1193 : END DO
1194 :
1195 : ! make sure the sends have completed before we deallocate
1196 :
1197 1852 : CALL mp_waitall(req(3:4))
1198 :
1199 1852 : DEALLOCATE (send_buf_3d_down)
1200 8508 : DEALLOCATE (send_buf_3d_up)
1201 : END DO
1202 :
1203 4804 : DEALLOCATE (dshifts)
1204 4804 : DEALLOCATE (ushifts)
1205 :
1206 : END IF
1207 :
1208 7408 : halo_swapped(idir) = .TRUE.
1209 :
1210 : END DO
1211 :
1212 : ! This is the real redistribution
1213 7408 : ALLOCATE (bounds(0:pw%pw_grid%para%group%num_pe - 1, 1:4))
1214 :
1215 : ! work out the pw grid points each proc holds
1216 5556 : DO i = 0, pw%pw_grid%para%group%num_pe - 1
1217 11112 : bounds(i, 1:2) = pw%pw_grid%para%bo(1:2, 1, i, 1)
1218 11112 : bounds(i, 3:4) = pw%pw_grid%para%bo(1:2, 2, i, 1)
1219 11112 : bounds(i, 1:2) = bounds(i, 1:2) - pw%pw_grid%npts(1)/2 - 1
1220 12964 : bounds(i, 3:4) = bounds(i, 3:4) - pw%pw_grid%npts(2)/2 - 1
1221 : END DO
1222 :
1223 7408 : ALLOCATE (send_tasks(0:pw%pw_grid%para%group%num_pe - 1, 1:6))
1224 5556 : ALLOCATE (send_sizes(0:pw%pw_grid%para%group%num_pe - 1))
1225 3704 : ALLOCATE (send_disps(0:pw%pw_grid%para%group%num_pe - 1))
1226 3704 : ALLOCATE (recv_tasks(0:pw%pw_grid%para%group%num_pe - 1, 1:6))
1227 3704 : ALLOCATE (recv_sizes(0:pw%pw_grid%para%group%num_pe - 1))
1228 3704 : ALLOCATE (recv_disps(0:pw%pw_grid%para%group%num_pe - 1))
1229 5556 : send_tasks(:, 1) = 1
1230 5556 : send_tasks(:, 2) = 0
1231 5556 : send_tasks(:, 3) = 1
1232 5556 : send_tasks(:, 4) = 0
1233 5556 : send_tasks(:, 5) = 1
1234 5556 : send_tasks(:, 6) = 0
1235 1852 : send_sizes = 0
1236 1852 : recv_sizes = 0
1237 :
1238 1852 : my_rs_rank = rs%desc%my_pos
1239 1852 : my_pw_rank = pw%pw_grid%para%group%mepos
1240 :
1241 : ! find the processors that should hold our data
1242 : ! should be part of the rs grid type
1243 : ! this is a loop over real ranks (i.e. the in-order cartesian ranks)
1244 : ! do the recv and send tasks in two separate loops which will
1245 : ! load balance better for OpenMP with large numbers of MPI tasks
1246 :
1247 : !$OMP PARALLEL DO DEFAULT(NONE), &
1248 : !$OMP PRIVATE(coords,idir,pos,lb_send,ub_send), &
1249 1852 : !$OMP SHARED(rs,bounds,my_rs_rank,recv_tasks,recv_sizes)
1250 : DO i = 0, rs%desc%group_size - 1
1251 :
1252 : coords(:) = rs%desc%rank2coord(:, rs%desc%real2virtual(i))
1253 : !calculate the rs grid points on each processor
1254 : !coords is the part of the grid that rank i actually holds
1255 : DO idir = 1, 3
1256 : pos(:) = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), coords(idir))
1257 : pos(:) = pos(:) - rs%desc%npts(idir)/2 - 1
1258 : lb_send(idir) = pos(1)
1259 : ub_send(idir) = pos(2)
1260 : END DO
1261 :
1262 : IF (lb_send(1) > bounds(my_rs_rank, 2)) CYCLE
1263 : IF (ub_send(1) < bounds(my_rs_rank, 1)) CYCLE
1264 : IF (lb_send(2) > bounds(my_rs_rank, 4)) CYCLE
1265 : IF (ub_send(2) < bounds(my_rs_rank, 3)) CYCLE
1266 :
1267 : recv_tasks(i, 1) = MAX(lb_send(1), bounds(my_rs_rank, 1))
1268 : recv_tasks(i, 2) = MIN(ub_send(1), bounds(my_rs_rank, 2))
1269 : recv_tasks(i, 3) = MAX(lb_send(2), bounds(my_rs_rank, 3))
1270 : recv_tasks(i, 4) = MIN(ub_send(2), bounds(my_rs_rank, 4))
1271 : recv_tasks(i, 5) = lb_send(3)
1272 : recv_tasks(i, 6) = ub_send(3)
1273 : recv_sizes(i) = (recv_tasks(i, 2) - recv_tasks(i, 1) + 1)* &
1274 : (recv_tasks(i, 4) - recv_tasks(i, 3) + 1)*(recv_tasks(i, 6) - recv_tasks(i, 5) + 1)
1275 :
1276 : END DO
1277 : !$OMP END PARALLEL DO
1278 :
1279 7408 : coords(:) = rs%desc%rank2coord(:, rs%desc%real2virtual(my_rs_rank))
1280 7408 : DO idir = 1, 3
1281 5556 : pos(:) = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), coords(idir))
1282 16668 : pos(:) = pos(:) - rs%desc%npts(idir)/2 - 1
1283 5556 : lb_send(idir) = pos(1)
1284 7408 : ub_send(idir) = pos(2)
1285 : END DO
1286 :
1287 1852 : lb_recv(:) = lb_send(:)
1288 1852 : ub_recv(:) = ub_send(:)
1289 : !$OMP PARALLEL DO DEFAULT(NONE), &
1290 1852 : !$OMP SHARED(pw,lb_send,ub_send,bounds,send_tasks,send_sizes)
1291 : DO j = 0, pw%pw_grid%para%group%num_pe - 1
1292 :
1293 : IF (lb_send(1) > bounds(j, 2)) CYCLE
1294 : IF (ub_send(1) < bounds(j, 1)) CYCLE
1295 : IF (lb_send(2) > bounds(j, 4)) CYCLE
1296 : IF (ub_send(2) < bounds(j, 3)) CYCLE
1297 :
1298 : send_tasks(j, 1) = MAX(lb_send(1), bounds(j, 1))
1299 : send_tasks(j, 2) = MIN(ub_send(1), bounds(j, 2))
1300 : send_tasks(j, 3) = MAX(lb_send(2), bounds(j, 3))
1301 : send_tasks(j, 4) = MIN(ub_send(2), bounds(j, 4))
1302 : send_tasks(j, 5) = lb_send(3)
1303 : send_tasks(j, 6) = ub_send(3)
1304 : send_sizes(j) = (send_tasks(j, 2) - send_tasks(j, 1) + 1)* &
1305 : (send_tasks(j, 4) - send_tasks(j, 3) + 1)*(send_tasks(j, 6) - send_tasks(j, 5) + 1)
1306 :
1307 : END DO
1308 : !$OMP END PARALLEL DO
1309 :
1310 1852 : send_disps(0) = 0
1311 1852 : recv_disps(0) = 0
1312 3704 : DO i = 1, pw%pw_grid%para%group%num_pe - 1
1313 1852 : send_disps(i) = send_disps(i - 1) + send_sizes(i - 1)
1314 3704 : recv_disps(i) = recv_disps(i - 1) + recv_sizes(i - 1)
1315 : END DO
1316 :
1317 12964 : CPASSERT(SUM(send_sizes) == PRODUCT(ub_recv - lb_recv + 1))
1318 :
1319 9260 : ALLOCATE (send_bufs(0:rs%desc%group_size - 1))
1320 11112 : ALLOCATE (recv_bufs(0:rs%desc%group_size - 1))
1321 :
1322 5556 : DO i = 0, rs%desc%group_size - 1
1323 3704 : IF (send_sizes(i) /= 0) THEN
1324 10374 : ALLOCATE (send_bufs(i)%array(send_sizes(i)))
1325 : ELSE
1326 246 : NULLIFY (send_bufs(i)%array)
1327 : END IF
1328 5556 : IF (recv_sizes(i) /= 0) THEN
1329 10374 : ALLOCATE (recv_bufs(i)%array(recv_sizes(i)))
1330 : ELSE
1331 246 : NULLIFY (recv_bufs(i)%array)
1332 : END IF
1333 : END DO
1334 :
1335 9260 : ALLOCATE (recv_reqs(0:rs%desc%group_size - 1))
1336 5556 : recv_reqs = mp_request_null
1337 :
1338 5556 : DO i = 0, rs%desc%group_size - 1
1339 5556 : IF (recv_sizes(i) /= 0) THEN
1340 3458 : CALL rs%desc%group%irecv(recv_bufs(i)%array, i, recv_reqs(i))
1341 : END IF
1342 : END DO
1343 :
1344 : ! do packing
1345 : !$OMP PARALLEL DO DEFAULT(NONE), &
1346 : !$OMP PRIVATE(k,z,y,x), &
1347 1852 : !$OMP SHARED(rs,send_tasks,send_bufs,send_disps)
1348 : DO i = 0, rs%desc%group_size - 1
1349 : k = 0
1350 : DO z = send_tasks(i, 5), send_tasks(i, 6)
1351 : DO y = send_tasks(i, 3), send_tasks(i, 4)
1352 : DO x = send_tasks(i, 1), send_tasks(i, 2)
1353 : k = k + 1
1354 : send_bufs(i)%array(k) = rs%r(x, y, z)
1355 : END DO
1356 : END DO
1357 : END DO
1358 : END DO
1359 : !$OMP END PARALLEL DO
1360 :
1361 9260 : ALLOCATE (send_reqs(0:rs%desc%group_size - 1))
1362 5556 : send_reqs = mp_request_null
1363 :
1364 5556 : DO i = 0, rs%desc%group_size - 1
1365 5556 : IF (send_sizes(i) /= 0) THEN
1366 3458 : CALL rs%desc%group%isend(send_bufs(i)%array, i, send_reqs(i))
1367 : END IF
1368 : END DO
1369 :
1370 : ! do unpacking
1371 : ! no OMP here so we can unpack each message as it arrives
1372 5556 : DO i = 0, rs%desc%group_size - 1
1373 3704 : IF (recv_sizes(i) == 0) CYCLE
1374 :
1375 3458 : CALL mp_waitany(recv_reqs, completed)
1376 3458 : k = 0
1377 143780 : DO z = recv_tasks(completed - 1, 5), recv_tasks(completed - 1, 6)
1378 10754290 : DO y = recv_tasks(completed - 1, 3), recv_tasks(completed - 1, 4)
1379 452245022 : DO x = recv_tasks(completed - 1, 1), recv_tasks(completed - 1, 2)
1380 441494436 : k = k + 1
1381 452106552 : pw%array(x, y, z) = recv_bufs(completed - 1)%array(k)
1382 : END DO
1383 : END DO
1384 : END DO
1385 : END DO
1386 :
1387 1852 : CALL mp_waitall(send_reqs)
1388 :
1389 1852 : DEALLOCATE (recv_reqs)
1390 1852 : DEALLOCATE (send_reqs)
1391 :
1392 5556 : DO i = 0, rs%desc%group_size - 1
1393 3704 : IF (ASSOCIATED(send_bufs(i)%array)) THEN
1394 3458 : DEALLOCATE (send_bufs(i)%array)
1395 : END IF
1396 5556 : IF (ASSOCIATED(recv_bufs(i)%array)) THEN
1397 3458 : DEALLOCATE (recv_bufs(i)%array)
1398 : END IF
1399 : END DO
1400 :
1401 1852 : DEALLOCATE (send_bufs)
1402 1852 : DEALLOCATE (recv_bufs)
1403 1852 : DEALLOCATE (send_tasks)
1404 1852 : DEALLOCATE (send_sizes)
1405 1852 : DEALLOCATE (send_disps)
1406 1852 : DEALLOCATE (recv_tasks)
1407 1852 : DEALLOCATE (recv_sizes)
1408 1852 : DEALLOCATE (recv_disps)
1409 :
1410 : IF (debug_this_module) THEN
1411 : ! safety check, to be removed once we're absolute sure the routine is correct
1412 : pw_sum = pw_integrate_function(pw)
1413 : IF (ABS(pw_sum - rs_sum)/MAX(1.0_dp, ABS(pw_sum), ABS(rs_sum)) > EPSILON(rs_sum)*1000) THEN
1414 : WRITE (error_string, '(A,6(1X,I4.4),3F25.16)') "rs_pw_transfer_distributed", &
1415 : rs%desc%npts, rs%desc%group_dim, pw_sum, rs_sum, ABS(pw_sum - rs_sum)
1416 : CALL cp_abort(__LOCATION__, &
1417 : error_string//" Please report this bug ... quick workaround: use "// &
1418 : "DISTRIBUTION_TYPE REPLICATED")
1419 : END IF
1420 : END IF
1421 :
1422 1852 : END SUBROUTINE transfer_rs2pw_distributed
1423 :
1424 : ! **************************************************************************************************
1425 : !> \brief does the pw2rs transfer in the case where the rs grid is
1426 : !> distributed (3D domain decomposition)
1427 : !> \param rs ...
1428 : !> \param pw ...
1429 : !> \par History
1430 : !> 12.2007 created [Matt Watkins]
1431 : !> 9.2008 reduced amount of halo data sent [Iain Bethune]
1432 : !> 10.2008 added non-blocking communication [Iain Bethune]
1433 : !> 4.2009 added support for rank-reordering on the grid [Iain Bethune]
1434 : !> 12.2009 added OMP and sparse alltoall [Iain Bethune]
1435 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2008-2009 on behalf of the HECToR project
1436 : !> \note
1437 : !> the transfer is a two step procedure. For example, for the rs2pw transfer:
1438 : !>
1439 : !> 1) Halo-exchange in 3D so that the local part of the rs_grid contains the full data
1440 : !> 2) an alltoall communication to redistribute the local rs_grid to the local pw_grid
1441 : !>
1442 : !> the halo exchange is most expensive on a large number of CPUs. Particular in this halo
1443 : !> exchange is that the border region is rather large (e.g. 20 points) and that it might overlap
1444 : !> with the central domain of several CPUs (i.e. next nearest neighbors)
1445 : ! **************************************************************************************************
1446 872 : SUBROUTINE transfer_pw2rs_distributed(rs, pw)
1447 : TYPE(realspace_grid_type), INTENT(IN) :: rs
1448 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
1449 :
1450 : INTEGER :: completed, dest_down, dest_up, i, idir, j, k, lb, my_id, my_pw_rank, my_rs_rank, &
1451 : n_shifts, nn, num_threads, position, source_down, source_up, ub, x, y, z
1452 872 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dshifts, recv_disps, recv_sizes, &
1453 872 : send_disps, send_sizes, ushifts
1454 1744 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: bounds, recv_tasks, send_tasks
1455 : INTEGER, DIMENSION(2) :: neighbours, pos
1456 : INTEGER, DIMENSION(3) :: coords, lb_recv, lb_recv_down, lb_recv_up, lb_send, lb_send_down, &
1457 : lb_send_up, ub_recv, ub_recv_down, ub_recv_up, ub_send, ub_send_down, ub_send_up
1458 : LOGICAL, DIMENSION(3) :: halo_swapped
1459 872 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: recv_buf_3d_down, recv_buf_3d_up, &
1460 872 : send_buf_3d_down, send_buf_3d_up
1461 1744 : TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: recv_bufs, send_bufs
1462 872 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_reqs, send_reqs
1463 4360 : TYPE(mp_request_type), DIMENSION(4) :: req
1464 :
1465 872 : num_threads = 1
1466 872 : my_id = 0
1467 :
1468 872 : CALL rs_grid_zero(rs)
1469 :
1470 : ! This is the real redistribution
1471 :
1472 3488 : ALLOCATE (bounds(0:pw%pw_grid%para%group%num_pe - 1, 1:4))
1473 :
1474 2616 : DO i = 0, pw%pw_grid%para%group%num_pe - 1
1475 5232 : bounds(i, 1:2) = pw%pw_grid%para%bo(1:2, 1, i, 1)
1476 5232 : bounds(i, 3:4) = pw%pw_grid%para%bo(1:2, 2, i, 1)
1477 5232 : bounds(i, 1:2) = bounds(i, 1:2) - pw%pw_grid%npts(1)/2 - 1
1478 6104 : bounds(i, 3:4) = bounds(i, 3:4) - pw%pw_grid%npts(2)/2 - 1
1479 : END DO
1480 :
1481 3488 : ALLOCATE (send_tasks(0:pw%pw_grid%para%group%num_pe - 1, 1:6))
1482 2616 : ALLOCATE (send_sizes(0:pw%pw_grid%para%group%num_pe - 1))
1483 1744 : ALLOCATE (send_disps(0:pw%pw_grid%para%group%num_pe - 1))
1484 1744 : ALLOCATE (recv_tasks(0:pw%pw_grid%para%group%num_pe - 1, 1:6))
1485 1744 : ALLOCATE (recv_sizes(0:pw%pw_grid%para%group%num_pe - 1))
1486 1744 : ALLOCATE (recv_disps(0:pw%pw_grid%para%group%num_pe - 1))
1487 :
1488 872 : send_tasks = 0
1489 2616 : send_tasks(:, 1) = 1
1490 2616 : send_tasks(:, 2) = 0
1491 2616 : send_tasks(:, 3) = 1
1492 2616 : send_tasks(:, 4) = 0
1493 2616 : send_tasks(:, 5) = 1
1494 2616 : send_tasks(:, 6) = 0
1495 872 : send_sizes = 0
1496 :
1497 872 : recv_tasks = 0
1498 2616 : recv_tasks(:, 1) = 1
1499 2616 : recv_tasks(:, 2) = 0
1500 2616 : send_tasks(:, 3) = 1
1501 2616 : send_tasks(:, 4) = 0
1502 2616 : send_tasks(:, 5) = 1
1503 2616 : send_tasks(:, 6) = 0
1504 872 : recv_sizes = 0
1505 :
1506 872 : my_rs_rank = rs%desc%my_pos
1507 872 : my_pw_rank = pw%pw_grid%para%group%mepos
1508 :
1509 : ! find the processors that should hold our data
1510 : ! should be part of the rs grid type
1511 : ! this is a loop over real ranks (i.e. the in-order cartesian ranks)
1512 : ! do the recv and send tasks in two separate loops which will
1513 : ! load balance better for OpenMP with large numbers of MPI tasks
1514 :
1515 : ! this is the reverse of rs2pw: what were the sends are now the recvs
1516 :
1517 : !$OMP PARALLEL DO DEFAULT(NONE), &
1518 : !$OMP PRIVATE(coords,idir,pos,lb_send,ub_send), &
1519 872 : !$OMP SHARED(rs,bounds,my_rs_rank,send_tasks,send_sizes,pw)
1520 : DO i = 0, pw%pw_grid%para%group%num_pe - 1
1521 :
1522 : coords(:) = rs%desc%rank2coord(:, rs%desc%real2virtual(i))
1523 : !calculate the real rs grid points on each processor
1524 : !coords is the part of the grid that rank i actually holds
1525 : DO idir = 1, 3
1526 : pos(:) = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), coords(idir))
1527 : pos(:) = pos(:) - rs%desc%npts(idir)/2 - 1
1528 : lb_send(idir) = pos(1)
1529 : ub_send(idir) = pos(2)
1530 : END DO
1531 :
1532 : IF (ub_send(1) < bounds(my_rs_rank, 1)) CYCLE
1533 : IF (lb_send(1) > bounds(my_rs_rank, 2)) CYCLE
1534 : IF (ub_send(2) < bounds(my_rs_rank, 3)) CYCLE
1535 : IF (lb_send(2) > bounds(my_rs_rank, 4)) CYCLE
1536 :
1537 : send_tasks(i, 1) = MAX(lb_send(1), bounds(my_rs_rank, 1))
1538 : send_tasks(i, 2) = MIN(ub_send(1), bounds(my_rs_rank, 2))
1539 : send_tasks(i, 3) = MAX(lb_send(2), bounds(my_rs_rank, 3))
1540 : send_tasks(i, 4) = MIN(ub_send(2), bounds(my_rs_rank, 4))
1541 : send_tasks(i, 5) = lb_send(3)
1542 : send_tasks(i, 6) = ub_send(3)
1543 : send_sizes(i) = (send_tasks(i, 2) - send_tasks(i, 1) + 1)* &
1544 : (send_tasks(i, 4) - send_tasks(i, 3) + 1)*(send_tasks(i, 6) - send_tasks(i, 5) + 1)
1545 :
1546 : END DO
1547 : !$OMP END PARALLEL DO
1548 :
1549 3488 : coords(:) = rs%desc%rank2coord(:, rs%desc%real2virtual(my_rs_rank))
1550 3488 : DO idir = 1, 3
1551 2616 : pos(:) = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), coords(idir))
1552 7848 : pos(:) = pos(:) - rs%desc%npts(idir)/2 - 1
1553 2616 : lb_send(idir) = pos(1)
1554 3488 : ub_send(idir) = pos(2)
1555 : END DO
1556 :
1557 872 : lb_recv(:) = lb_send(:)
1558 872 : ub_recv(:) = ub_send(:)
1559 :
1560 : !$OMP PARALLEL DO DEFAULT(NONE), &
1561 872 : !$OMP SHARED(pw,lb_send,ub_send,bounds,recv_tasks,recv_sizes)
1562 : DO j = 0, pw%pw_grid%para%group%num_pe - 1
1563 :
1564 : IF (ub_send(1) < bounds(j, 1)) CYCLE
1565 : IF (lb_send(1) > bounds(j, 2)) CYCLE
1566 : IF (ub_send(2) < bounds(j, 3)) CYCLE
1567 : IF (lb_send(2) > bounds(j, 4)) CYCLE
1568 :
1569 : recv_tasks(j, 1) = MAX(lb_send(1), bounds(j, 1))
1570 : recv_tasks(j, 2) = MIN(ub_send(1), bounds(j, 2))
1571 : recv_tasks(j, 3) = MAX(lb_send(2), bounds(j, 3))
1572 : recv_tasks(j, 4) = MIN(ub_send(2), bounds(j, 4))
1573 : recv_tasks(j, 5) = lb_send(3)
1574 : recv_tasks(j, 6) = ub_send(3)
1575 : recv_sizes(j) = (recv_tasks(j, 2) - recv_tasks(j, 1) + 1)* &
1576 : (recv_tasks(j, 4) - recv_tasks(j, 3) + 1)*(recv_tasks(j, 6) - recv_tasks(j, 5) + 1)
1577 :
1578 : END DO
1579 : !$OMP END PARALLEL DO
1580 :
1581 872 : send_disps(0) = 0
1582 872 : recv_disps(0) = 0
1583 1744 : DO i = 1, pw%pw_grid%para%group%num_pe - 1
1584 872 : send_disps(i) = send_disps(i - 1) + send_sizes(i - 1)
1585 1744 : recv_disps(i) = recv_disps(i - 1) + recv_sizes(i - 1)
1586 : END DO
1587 :
1588 6104 : CPASSERT(SUM(recv_sizes) == PRODUCT(ub_recv - lb_recv + 1))
1589 :
1590 4360 : ALLOCATE (send_bufs(0:rs%desc%group_size - 1))
1591 5232 : ALLOCATE (recv_bufs(0:rs%desc%group_size - 1))
1592 :
1593 2616 : DO i = 0, rs%desc%group_size - 1
1594 1744 : IF (send_sizes(i) /= 0) THEN
1595 4986 : ALLOCATE (send_bufs(i)%array(send_sizes(i)))
1596 : ELSE
1597 82 : NULLIFY (send_bufs(i)%array)
1598 : END IF
1599 2616 : IF (recv_sizes(i) /= 0) THEN
1600 4986 : ALLOCATE (recv_bufs(i)%array(recv_sizes(i)))
1601 : ELSE
1602 82 : NULLIFY (recv_bufs(i)%array)
1603 : END IF
1604 : END DO
1605 :
1606 4360 : ALLOCATE (recv_reqs(0:rs%desc%group_size - 1))
1607 2616 : recv_reqs = mp_request_null
1608 :
1609 2616 : DO i = 0, rs%desc%group_size - 1
1610 2616 : IF (recv_sizes(i) /= 0) THEN
1611 1662 : CALL rs%desc%group%irecv(recv_bufs(i)%array, i, recv_reqs(i))
1612 : END IF
1613 : END DO
1614 :
1615 : ! do packing
1616 : !$OMP PARALLEL DO DEFAULT(NONE), &
1617 : !$OMP PRIVATE(k,z,y,x), &
1618 872 : !$OMP SHARED(pw,rs,send_tasks,send_bufs,send_disps)
1619 : DO i = 0, rs%desc%group_size - 1
1620 : k = 0
1621 : DO z = send_tasks(i, 5), send_tasks(i, 6)
1622 : DO y = send_tasks(i, 3), send_tasks(i, 4)
1623 : DO x = send_tasks(i, 1), send_tasks(i, 2)
1624 : k = k + 1
1625 : send_bufs(i)%array(k) = pw%array(x, y, z)
1626 : END DO
1627 : END DO
1628 : END DO
1629 : END DO
1630 : !$OMP END PARALLEL DO
1631 :
1632 4360 : ALLOCATE (send_reqs(0:rs%desc%group_size - 1))
1633 2616 : send_reqs = mp_request_null
1634 :
1635 2616 : DO i = 0, rs%desc%group_size - 1
1636 2616 : IF (send_sizes(i) /= 0) THEN
1637 1662 : CALL rs%desc%group%isend(send_bufs(i)%array, i, send_reqs(i))
1638 : END IF
1639 : END DO
1640 :
1641 : ! do unpacking
1642 : ! no OMP here so we can unpack each message as it arrives
1643 :
1644 2616 : DO i = 0, rs%desc%group_size - 1
1645 1744 : IF (recv_sizes(i) == 0) CYCLE
1646 :
1647 1662 : CALL mp_waitany(recv_reqs, completed)
1648 1662 : k = 0
1649 66556 : DO z = recv_tasks(completed - 1, 5), recv_tasks(completed - 1, 6)
1650 4763300 : DO y = recv_tasks(completed - 1, 3), recv_tasks(completed - 1, 4)
1651 194725231 : DO x = recv_tasks(completed - 1, 1), recv_tasks(completed - 1, 2)
1652 189963675 : k = k + 1
1653 194661209 : rs%r(x, y, z) = recv_bufs(completed - 1)%array(k)
1654 : END DO
1655 : END DO
1656 : END DO
1657 : END DO
1658 :
1659 872 : CALL mp_waitall(send_reqs)
1660 :
1661 872 : DEALLOCATE (recv_reqs)
1662 872 : DEALLOCATE (send_reqs)
1663 :
1664 2616 : DO i = 0, rs%desc%group_size - 1
1665 1744 : IF (ASSOCIATED(send_bufs(i)%array)) THEN
1666 1662 : DEALLOCATE (send_bufs(i)%array)
1667 : END IF
1668 2616 : IF (ASSOCIATED(recv_bufs(i)%array)) THEN
1669 1662 : DEALLOCATE (recv_bufs(i)%array)
1670 : END IF
1671 : END DO
1672 :
1673 872 : DEALLOCATE (send_bufs)
1674 872 : DEALLOCATE (recv_bufs)
1675 872 : DEALLOCATE (send_tasks)
1676 872 : DEALLOCATE (send_sizes)
1677 872 : DEALLOCATE (send_disps)
1678 872 : DEALLOCATE (recv_tasks)
1679 872 : DEALLOCATE (recv_sizes)
1680 872 : DEALLOCATE (recv_disps)
1681 :
1682 : ! now pass wings around
1683 872 : halo_swapped = .FALSE.
1684 :
1685 3488 : DO idir = 1, 3
1686 :
1687 2616 : IF (rs%desc%perd(idir) /= 1) THEN
1688 :
1689 5568 : ALLOCATE (dshifts(0:rs%desc%neighbours(idir)))
1690 3712 : ALLOCATE (ushifts(0:rs%desc%neighbours(idir)))
1691 1856 : ushifts = 0
1692 1856 : dshifts = 0
1693 :
1694 3712 : DO n_shifts = 1, rs%desc%neighbours(idir)
1695 :
1696 : ! need to take into account the possible varying widths of neighbouring cells
1697 : ! ushifts and dshifts hold the real size of the neighbouring cells
1698 :
1699 1856 : position = MODULO(rs%desc%virtual_group_coor(idir) - n_shifts, rs%desc%group_dim(idir))
1700 1856 : neighbours = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), position)
1701 1856 : dshifts(n_shifts) = dshifts(n_shifts - 1) + (neighbours(2) - neighbours(1) + 1)
1702 :
1703 1856 : position = MODULO(rs%desc%virtual_group_coor(idir) + n_shifts, rs%desc%group_dim(idir))
1704 1856 : neighbours = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), position)
1705 1856 : ushifts(n_shifts) = ushifts(n_shifts - 1) + (neighbours(2) - neighbours(1) + 1)
1706 :
1707 : ! The border data has to be send/received from the neighbors
1708 : ! First we calculate the source and destination processes for the shift
1709 : ! The first shift is "downwards"
1710 :
1711 1856 : CALL cart_shift(rs, idir, -1*n_shifts, source_down, dest_down)
1712 :
1713 7424 : lb_send_down(:) = rs%lb_local(:)
1714 7424 : ub_send_down(:) = rs%ub_local(:)
1715 7424 : lb_recv_down(:) = rs%lb_local(:)
1716 7424 : ub_recv_down(:) = rs%ub_local(:)
1717 :
1718 1856 : IF (dshifts(n_shifts - 1) <= rs%desc%border) THEN
1719 1856 : lb_send_down(idir) = lb_send_down(idir) + rs%desc%border
1720 : ub_send_down(idir) = MIN(ub_send_down(idir) - rs%desc%border, &
1721 1856 : lb_send_down(idir) + rs%desc%border - 1 - dshifts(n_shifts - 1))
1722 :
1723 1856 : lb_recv_down(idir) = ub_recv_down(idir) - rs%desc%border + 1 + ushifts(n_shifts - 1)
1724 : ub_recv_down(idir) = MIN(ub_recv_down(idir), &
1725 1856 : ub_recv_down(idir) - rs%desc%border + ushifts(n_shifts))
1726 : ELSE
1727 0 : lb_send_down(idir) = 0
1728 0 : ub_send_down(idir) = -1
1729 0 : lb_recv_down(idir) = 0
1730 0 : ub_recv_down(idir) = -1
1731 : END IF
1732 :
1733 7424 : DO i = 1, 3
1734 7424 : IF (.NOT. (halo_swapped(i) .OR. i == idir)) THEN
1735 1560 : lb_send_down(i) = rs%lb_real(i)
1736 1560 : ub_send_down(i) = rs%ub_real(i)
1737 1560 : lb_recv_down(i) = rs%lb_real(i)
1738 1560 : ub_recv_down(i) = rs%ub_real(i)
1739 : END IF
1740 : END DO
1741 :
1742 : ! allocate the recv buffer
1743 7424 : nn = PRODUCT(ub_recv_down - lb_recv_down + 1)
1744 0 : ALLOCATE (recv_buf_3d_down(lb_recv_down(1):ub_recv_down(1), &
1745 9280 : lb_recv_down(2):ub_recv_down(2), lb_recv_down(3):ub_recv_down(3)))
1746 :
1747 : ! recv buffer is now ready, so post the receive
1748 1856 : CALL rs%desc%group%irecv(recv_buf_3d_down, source_down, req(1))
1749 :
1750 : ! now allocate,pack and send the send buffer
1751 7424 : nn = PRODUCT(ub_send_down - lb_send_down + 1)
1752 0 : ALLOCATE (send_buf_3d_down(lb_send_down(1):ub_send_down(1), &
1753 9280 : lb_send_down(2):ub_send_down(2), lb_send_down(3):ub_send_down(3)))
1754 :
1755 : !$OMP PARALLEL DEFAULT(NONE), &
1756 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1757 1856 : !$OMP SHARED(send_buf_3d_down,rs,lb_send_down,ub_send_down)
1758 : !$ num_threads = MIN(omp_get_max_threads(), ub_send_down(3) - lb_send_down(3) + 1)
1759 : !$ my_id = omp_get_thread_num()
1760 : IF (my_id < num_threads) THEN
1761 : lb = lb_send_down(3) + ((ub_send_down(3) - lb_send_down(3) + 1)*my_id)/num_threads
1762 : ub = lb_send_down(3) + ((ub_send_down(3) - lb_send_down(3) + 1)*(my_id + 1))/num_threads - 1
1763 :
1764 : send_buf_3d_down(lb_send_down(1):ub_send_down(1), lb_send_down(2):ub_send_down(2), &
1765 : lb:ub) = rs%r(lb_send_down(1):ub_send_down(1), &
1766 : lb_send_down(2):ub_send_down(2), lb:ub)
1767 : END IF
1768 : !$OMP END PARALLEL
1769 :
1770 1856 : CALL rs%desc%group%isend(send_buf_3d_down, dest_down, req(3))
1771 :
1772 : ! Now for the other direction
1773 :
1774 1856 : CALL cart_shift(rs, idir, n_shifts, source_up, dest_up)
1775 :
1776 7424 : lb_send_up(:) = rs%lb_local(:)
1777 7424 : ub_send_up(:) = rs%ub_local(:)
1778 7424 : lb_recv_up(:) = rs%lb_local(:)
1779 7424 : ub_recv_up(:) = rs%ub_local(:)
1780 :
1781 1856 : IF (ushifts(n_shifts - 1) <= rs%desc%border) THEN
1782 1856 : ub_send_up(idir) = ub_send_up(idir) - rs%desc%border
1783 : lb_send_up(idir) = MAX(lb_send_up(idir) + rs%desc%border, &
1784 1856 : ub_send_up(idir) - rs%desc%border + 1 + ushifts(n_shifts - 1))
1785 :
1786 1856 : ub_recv_up(idir) = lb_recv_up(idir) + rs%desc%border - 1 - dshifts(n_shifts - 1)
1787 : lb_recv_up(idir) = MAX(lb_recv_up(idir), &
1788 1856 : lb_recv_up(idir) + rs%desc%border - dshifts(n_shifts))
1789 : ELSE
1790 0 : lb_send_up(idir) = 0
1791 0 : ub_send_up(idir) = -1
1792 0 : lb_recv_up(idir) = 0
1793 0 : ub_recv_up(idir) = -1
1794 : END IF
1795 :
1796 7424 : DO i = 1, 3
1797 7424 : IF (.NOT. (halo_swapped(i) .OR. i == idir)) THEN
1798 1560 : lb_send_up(i) = rs%lb_real(i)
1799 1560 : ub_send_up(i) = rs%ub_real(i)
1800 1560 : lb_recv_up(i) = rs%lb_real(i)
1801 1560 : ub_recv_up(i) = rs%ub_real(i)
1802 : END IF
1803 : END DO
1804 :
1805 : ! allocate the recv buffer
1806 7424 : nn = PRODUCT(ub_recv_up - lb_recv_up + 1)
1807 0 : ALLOCATE (recv_buf_3d_up(lb_recv_up(1):ub_recv_up(1), &
1808 9280 : lb_recv_up(2):ub_recv_up(2), lb_recv_up(3):ub_recv_up(3)))
1809 :
1810 : ! recv buffer is now ready, so post the receive
1811 :
1812 1856 : CALL rs%desc%group%irecv(recv_buf_3d_up, source_up, req(2))
1813 :
1814 : ! now allocate,pack and send the send buffer
1815 7424 : nn = PRODUCT(ub_send_up - lb_send_up + 1)
1816 0 : ALLOCATE (send_buf_3d_up(lb_send_up(1):ub_send_up(1), &
1817 9280 : lb_send_up(2):ub_send_up(2), lb_send_up(3):ub_send_up(3)))
1818 :
1819 : !$OMP PARALLEL DEFAULT(NONE), &
1820 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1821 1856 : !$OMP SHARED(send_buf_3d_up,rs,lb_send_up,ub_send_up)
1822 : !$ num_threads = MIN(omp_get_max_threads(), ub_send_up(3) - lb_send_up(3) + 1)
1823 : !$ my_id = omp_get_thread_num()
1824 : IF (my_id < num_threads) THEN
1825 : lb = lb_send_up(3) + ((ub_send_up(3) - lb_send_up(3) + 1)*my_id)/num_threads
1826 : ub = lb_send_up(3) + ((ub_send_up(3) - lb_send_up(3) + 1)*(my_id + 1))/num_threads - 1
1827 :
1828 : send_buf_3d_up(lb_send_up(1):ub_send_up(1), lb_send_up(2):ub_send_up(2), &
1829 : lb:ub) = rs%r(lb_send_up(1):ub_send_up(1), &
1830 : lb_send_up(2):ub_send_up(2), lb:ub)
1831 : END IF
1832 : !$OMP END PARALLEL
1833 :
1834 1856 : CALL rs%desc%group%isend(send_buf_3d_up, dest_up, req(4))
1835 :
1836 : ! wait for a recv to complete, then we can unpack
1837 :
1838 5568 : DO i = 1, 2
1839 :
1840 3712 : CALL mp_waitany(req(1:2), completed)
1841 :
1842 5568 : IF (completed == 1) THEN
1843 :
1844 : ! only some procs may need later shifts
1845 1856 : IF (ub_recv_down(idir) >= lb_recv_down(idir)) THEN
1846 :
1847 : ! Add the data to the RS Grid
1848 : !$OMP PARALLEL DEFAULT(NONE), &
1849 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1850 1856 : !$OMP SHARED(recv_buf_3d_down,rs,lb_recv_down,ub_recv_down)
1851 : !$ num_threads = MIN(omp_get_max_threads(), ub_recv_down(3) - lb_recv_down(3) + 1)
1852 : !$ my_id = omp_get_thread_num()
1853 : IF (my_id < num_threads) THEN
1854 : lb = lb_recv_down(3) + ((ub_recv_down(3) - lb_recv_down(3) + 1)*my_id)/num_threads
1855 : ub = lb_recv_down(3) + ((ub_recv_down(3) - lb_recv_down(3) + 1)*(my_id + 1))/num_threads - 1
1856 :
1857 : rs%r(lb_recv_down(1):ub_recv_down(1), lb_recv_down(2):ub_recv_down(2), &
1858 : lb:ub) = recv_buf_3d_down(:, :, lb:ub)
1859 : END IF
1860 : !$OMP END PARALLEL
1861 : END IF
1862 :
1863 1856 : DEALLOCATE (recv_buf_3d_down)
1864 : ELSE
1865 :
1866 : ! only some procs may need later shifts
1867 1856 : IF (ub_recv_up(idir) >= lb_recv_up(idir)) THEN
1868 :
1869 : ! Add the data to the RS Grid
1870 : !$OMP PARALLEL DEFAULT(NONE), &
1871 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1872 1856 : !$OMP SHARED(recv_buf_3d_up,rs,lb_recv_up,ub_recv_up)
1873 : !$ num_threads = MIN(omp_get_max_threads(), ub_recv_up(3) - lb_recv_up(3) + 1)
1874 : !$ my_id = omp_get_thread_num()
1875 : IF (my_id < num_threads) THEN
1876 : lb = lb_recv_up(3) + ((ub_recv_up(3) - lb_recv_up(3) + 1)*my_id)/num_threads
1877 : ub = lb_recv_up(3) + ((ub_recv_up(3) - lb_recv_up(3) + 1)*(my_id + 1))/num_threads - 1
1878 :
1879 : rs%r(lb_recv_up(1):ub_recv_up(1), lb_recv_up(2):ub_recv_up(2), &
1880 : lb:ub) = recv_buf_3d_up(:, :, lb:ub)
1881 : END IF
1882 : !$OMP END PARALLEL
1883 : END IF
1884 :
1885 1856 : DEALLOCATE (recv_buf_3d_up)
1886 : END IF
1887 : END DO
1888 :
1889 1856 : CALL mp_waitall(req(3:4))
1890 :
1891 1856 : DEALLOCATE (send_buf_3d_down)
1892 5568 : DEALLOCATE (send_buf_3d_up)
1893 : END DO
1894 :
1895 1856 : DEALLOCATE (ushifts)
1896 1856 : DEALLOCATE (dshifts)
1897 : END IF
1898 :
1899 3488 : halo_swapped(idir) = .TRUE.
1900 :
1901 : END DO
1902 :
1903 872 : END SUBROUTINE transfer_pw2rs_distributed
1904 :
1905 : ! **************************************************************************************************
1906 : !> \brief Initialize grid to zero
1907 : !> \param rs ...
1908 : !> \par History
1909 : !> none
1910 : !> \author JGH (23-Mar-2002)
1911 : ! **************************************************************************************************
1912 412122 : SUBROUTINE rs_grid_zero(rs)
1913 :
1914 : TYPE(realspace_grid_type), INTENT(IN) :: rs
1915 :
1916 : CHARACTER(len=*), PARAMETER :: routineN = 'rs_grid_zero'
1917 :
1918 : INTEGER :: handle, i, j, k, l(3), u(3)
1919 :
1920 412122 : CALL timeset(routineN, handle)
1921 1236366 : l(1) = LBOUND(rs%r, 1); l(2) = LBOUND(rs%r, 2); l(3) = LBOUND(rs%r, 3)
1922 1236366 : u(1) = UBOUND(rs%r, 1); u(2) = UBOUND(rs%r, 2); u(3) = UBOUND(rs%r, 3)
1923 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(3) &
1924 : !$OMP PRIVATE(i,j,k) &
1925 412122 : !$OMP SHARED(rs,l,u)
1926 : DO k = l(3), u(3)
1927 : DO j = l(2), u(2)
1928 : DO i = l(1), u(1)
1929 : rs%r(i, j, k) = 0.0_dp
1930 : END DO
1931 : END DO
1932 : END DO
1933 : !$OMP END PARALLEL DO
1934 412122 : CALL timestop(handle)
1935 :
1936 412122 : END SUBROUTINE rs_grid_zero
1937 :
1938 : ! **************************************************************************************************
1939 : !> \brief rs1(i) = rs1(i) + rs2(i)*rs3(i)
1940 : !> \param rs1 ...
1941 : !> \param rs2 ...
1942 : !> \param rs3 ...
1943 : !> \param scalar ...
1944 : !> \par History
1945 : !> none
1946 : !> \author
1947 : ! **************************************************************************************************
1948 1404 : SUBROUTINE rs_grid_mult_and_add(rs1, rs2, rs3, scalar)
1949 :
1950 : TYPE(realspace_grid_type), INTENT(IN) :: rs1, rs2, rs3
1951 : REAL(dp), INTENT(IN) :: scalar
1952 :
1953 : CHARACTER(len=*), PARAMETER :: routineN = 'rs_grid_mult_and_add'
1954 :
1955 : INTEGER :: handle, i, j, k, l(3), u(3)
1956 :
1957 : !-----------------------------------------------------------------------------!
1958 :
1959 1404 : CALL timeset(routineN, handle)
1960 1404 : IF (scalar /= 0.0_dp) THEN
1961 4212 : l(1) = LBOUND(rs1%r, 1); l(2) = LBOUND(rs1%r, 2); l(3) = LBOUND(rs1%r, 3)
1962 4212 : u(1) = UBOUND(rs1%r, 1); u(2) = UBOUND(rs1%r, 2); u(3) = UBOUND(rs1%r, 3)
1963 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(3) &
1964 : !$OMP PRIVATE(i,j,k) &
1965 1404 : !$OMP SHARED(rs1,rs2,rs3,scalar,l,u)
1966 : DO k = l(3), u(3)
1967 : DO j = l(2), u(2)
1968 : DO i = l(1), u(1)
1969 : rs1%r(i, j, k) = rs1%r(i, j, k) + scalar*rs2%r(i, j, k)*rs3%r(i, j, k)
1970 : END DO
1971 : END DO
1972 : END DO
1973 : !$OMP END PARALLEL DO
1974 : END IF
1975 1404 : CALL timestop(handle)
1976 1404 : END SUBROUTINE rs_grid_mult_and_add
1977 :
1978 : ! **************************************************************************************************
1979 : !> \brief Set box matrix info for real space grid
1980 : !> This is needed for variable cell simulations
1981 : !> \param pw_grid ...
1982 : !> \param rs ...
1983 : !> \par History
1984 : !> none
1985 : !> \author JGH (15-May-2007)
1986 : ! **************************************************************************************************
1987 207980 : SUBROUTINE rs_grid_set_box(pw_grid, rs)
1988 :
1989 : TYPE(pw_grid_type), INTENT(IN), TARGET :: pw_grid
1990 : TYPE(realspace_grid_type), INTENT(IN) :: rs
1991 :
1992 207980 : CPASSERT(ASSOCIATED(rs%desc%pw, pw_grid))
1993 2703740 : rs%desc%dh = pw_grid%dh
1994 2703740 : rs%desc%dh_inv = pw_grid%dh_inv
1995 :
1996 207980 : END SUBROUTINE rs_grid_set_box
1997 :
1998 : ! **************************************************************************************************
1999 : !> \brief retains the given rs grid descriptor (see doc/ReferenceCounting.html)
2000 : !> \param rs_desc the grid descriptor to retain
2001 : !> \par History
2002 : !> 04.2009 created [Iain Bethune]
2003 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
2004 : ! **************************************************************************************************
2005 291514 : SUBROUTINE rs_grid_retain_descriptor(rs_desc)
2006 : TYPE(realspace_grid_desc_type), INTENT(INOUT) :: rs_desc
2007 :
2008 291514 : CPASSERT(rs_desc%ref_count > 0)
2009 291514 : rs_desc%ref_count = rs_desc%ref_count + 1
2010 291514 : END SUBROUTINE rs_grid_retain_descriptor
2011 :
2012 : ! **************************************************************************************************
2013 : !> \brief releases the given rs grid (see doc/ReferenceCounting.html)
2014 : !> \param rs_grid the rs grid to release
2015 : !> \par History
2016 : !> 03.2003 created [fawzi]
2017 : !> \author fawzi
2018 : ! **************************************************************************************************
2019 290630 : SUBROUTINE rs_grid_release(rs_grid)
2020 : TYPE(realspace_grid_type), INTENT(INOUT) :: rs_grid
2021 :
2022 290630 : CALL rs_grid_release_descriptor(rs_grid%desc)
2023 :
2024 290630 : CALL offload_free_buffer(rs_grid%buffer)
2025 290630 : NULLIFY (rs_grid%r)
2026 :
2027 290630 : IF (ALLOCATED(rs_grid%px)) DEALLOCATE (rs_grid%px)
2028 290630 : IF (ALLOCATED(rs_grid%py)) DEALLOCATE (rs_grid%py)
2029 290630 : IF (ALLOCATED(rs_grid%pz)) DEALLOCATE (rs_grid%pz)
2030 290630 : END SUBROUTINE rs_grid_release
2031 :
2032 : ! **************************************************************************************************
2033 : !> \brief releases the given rs grid descriptor (see doc/ReferenceCounting.html)
2034 : !> \param rs_desc the rs grid descriptor to release
2035 : !> \par History
2036 : !> 04.2009 created [Iain Bethune]
2037 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
2038 : ! **************************************************************************************************
2039 336597 : SUBROUTINE rs_grid_release_descriptor(rs_desc)
2040 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
2041 :
2042 336597 : IF (ASSOCIATED(rs_desc)) THEN
2043 332844 : CPASSERT(rs_desc%ref_count > 0)
2044 332844 : rs_desc%ref_count = rs_desc%ref_count - 1
2045 332844 : IF (rs_desc%ref_count == 0) THEN
2046 :
2047 41330 : CALL pw_grid_release(rs_desc%pw)
2048 :
2049 41330 : IF (rs_desc%parallel) THEN
2050 : ! release the group communicator
2051 37102 : CALL rs_desc%group%free()
2052 :
2053 37102 : DEALLOCATE (rs_desc%virtual2real)
2054 37102 : DEALLOCATE (rs_desc%real2virtual)
2055 : END IF
2056 :
2057 41330 : IF (rs_desc%distributed) THEN
2058 158 : DEALLOCATE (rs_desc%rank2coord)
2059 158 : DEALLOCATE (rs_desc%coord2rank)
2060 158 : DEALLOCATE (rs_desc%lb_global)
2061 158 : DEALLOCATE (rs_desc%ub_global)
2062 158 : DEALLOCATE (rs_desc%x2coord)
2063 158 : DEALLOCATE (rs_desc%y2coord)
2064 158 : DEALLOCATE (rs_desc%z2coord)
2065 : END IF
2066 :
2067 41330 : DEALLOCATE (rs_desc)
2068 : END IF
2069 : END IF
2070 336597 : NULLIFY (rs_desc)
2071 336597 : END SUBROUTINE rs_grid_release_descriptor
2072 :
2073 : ! **************************************************************************************************
2074 : !> \brief emulates the function of an MPI_cart_shift operation, but the shift is
2075 : !> done in virtual coordinates, and the corresponding real ranks are returned
2076 : !> \param rs_grid ...
2077 : !> \param dir ...
2078 : !> \param disp ...
2079 : !> \param source ...
2080 : !> \param dest ...
2081 : !> \par History
2082 : !> 04.2009 created [Iain Bethune]
2083 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
2084 : ! **************************************************************************************************
2085 7416 : PURE SUBROUTINE cart_shift(rs_grid, dir, disp, source, dest)
2086 :
2087 : TYPE(realspace_grid_type), INTENT(IN) :: rs_grid
2088 : INTEGER, INTENT(IN) :: dir, disp
2089 : INTEGER, INTENT(OUT) :: source, dest
2090 :
2091 : INTEGER, DIMENSION(3) :: shift_coords
2092 :
2093 29664 : shift_coords = rs_grid%desc%virtual_group_coor
2094 7416 : shift_coords(dir) = MODULO(shift_coords(dir) + disp, rs_grid%desc%group_dim(dir))
2095 7416 : dest = rs_grid%desc%virtual2real(rs_grid%desc%coord2rank(shift_coords(1), shift_coords(2), shift_coords(3)))
2096 29664 : shift_coords = rs_grid%desc%virtual_group_coor
2097 7416 : shift_coords(dir) = MODULO(shift_coords(dir) - disp, rs_grid%desc%group_dim(dir))
2098 7416 : source = rs_grid%desc%virtual2real(rs_grid%desc%coord2rank(shift_coords(1), shift_coords(2), shift_coords(3)))
2099 :
2100 7416 : END SUBROUTINE cart_shift
2101 :
2102 : ! **************************************************************************************************
2103 : !> \brief returns the maximum number of points in the local grid of any process
2104 : !> to account for the case where the grid may later be reordered
2105 : !> \param desc ...
2106 : !> \return ...
2107 : !> \par History
2108 : !> 10.2011 created [Iain Bethune]
2109 : ! **************************************************************************************************
2110 0 : FUNCTION rs_grid_max_ngpts(desc) RESULT(max_ngpts)
2111 : TYPE(realspace_grid_desc_type), INTENT(IN) :: desc
2112 : INTEGER :: max_ngpts
2113 :
2114 : CHARACTER(len=*), PARAMETER :: routineN = 'rs_grid_max_ngpts'
2115 :
2116 : INTEGER :: handle, i
2117 : INTEGER, DIMENSION(3) :: lb, ub
2118 :
2119 0 : CALL timeset(routineN, handle)
2120 :
2121 0 : max_ngpts = 0
2122 0 : IF ((desc%pw%para%mode == PW_MODE_LOCAL) .OR. &
2123 : (ALL(desc%group_dim == 1))) THEN
2124 0 : CPASSERT(PRODUCT(INT(desc%npts, KIND=int_8)) < HUGE(1))
2125 0 : max_ngpts = PRODUCT(desc%npts)
2126 : ELSE
2127 0 : DO i = 0, desc%group_size - 1
2128 0 : lb = desc%lb_global(:, i)
2129 0 : ub = desc%ub_global(:, i)
2130 0 : lb = lb - desc%border*(1 - desc%perd)
2131 0 : ub = ub + desc%border*(1 - desc%perd)
2132 0 : CPASSERT(PRODUCT(INT(ub - lb + 1, KIND=int_8)) < HUGE(1))
2133 0 : max_ngpts = MAX(max_ngpts, PRODUCT(ub - lb + 1))
2134 : END DO
2135 : END IF
2136 :
2137 0 : CALL timestop(handle)
2138 :
2139 0 : END FUNCTION rs_grid_max_ngpts
2140 :
2141 : ! **************************************************************************************************
2142 : !> \brief ...
2143 : !> \param rs_grid ...
2144 : !> \param h_inv ...
2145 : !> \param ra ...
2146 : !> \param offset ...
2147 : !> \param group_size ...
2148 : !> \param my_pos ...
2149 : !> \return ...
2150 : ! **************************************************************************************************
2151 1611955 : PURE LOGICAL FUNCTION map_gaussian_here(rs_grid, h_inv, ra, offset, group_size, my_pos) RESULT(res)
2152 : TYPE(realspace_grid_type), INTENT(IN) :: rs_grid
2153 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: h_inv
2154 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra
2155 : INTEGER, INTENT(IN), OPTIONAL :: offset, group_size, my_pos
2156 :
2157 : INTEGER :: dir, lb(3), location(3), tp(3), ub(3)
2158 :
2159 1611955 : res = .FALSE.
2160 :
2161 6447812 : IF (.NOT. ALL(rs_grid%desc%perd == 1)) THEN
2162 32 : DO dir = 1, 3
2163 : ! bounds of local grid (i.e. removing the 'wings'), if periodic
2164 96 : tp(dir) = FLOOR(DOT_PRODUCT(h_inv(dir, :), ra)*rs_grid%desc%npts(dir))
2165 24 : tp(dir) = MODULO(tp(dir), rs_grid%desc%npts(dir))
2166 24 : IF (rs_grid%desc%perd(dir) /= 1) THEN
2167 8 : lb(dir) = rs_grid%lb_local(dir) + rs_grid%desc%border
2168 8 : ub(dir) = rs_grid%ub_local(dir) - rs_grid%desc%border
2169 : ELSE
2170 16 : lb(dir) = rs_grid%lb_local(dir)
2171 16 : ub(dir) = rs_grid%ub_local(dir)
2172 : END IF
2173 : ! distributed grid, only map if it is local to the grid
2174 32 : location(dir) = tp(dir) + rs_grid%desc%lb(dir)
2175 : END DO
2176 60 : IF (ALL(lb(:) <= location(:)) .AND. ALL(location(:) <= ub(:))) THEN
2177 4 : res = .TRUE.
2178 : END IF
2179 : ELSE
2180 1611947 : IF (PRESENT(offset) .AND. PRESENT(group_size) .AND. PRESENT(my_pos)) THEN
2181 : ! not distributed, just a round-robin distribution over the full set of CPUs
2182 1611947 : IF (MODULO(offset, group_size) == my_pos) res = .TRUE.
2183 : END IF
2184 : END IF
2185 :
2186 1611955 : END FUNCTION map_gaussian_here
2187 :
2188 0 : END MODULE realspace_grid_types
|
