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 39308 : 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 39308 : CALL timeset(routineN, handle)
216 :
217 39308 : IF (PRESENT(border_points)) THEN
218 128 : border_size = border_points
219 : ELSE
220 : border_size = 0
221 : END IF
222 :
223 1926092 : ALLOCATE (desc)
224 :
225 39308 : CALL pw_grid%para%group%sync()
226 :
227 39308 : desc%pw => pw_grid
228 39308 : CALL pw_grid_retain(desc%pw)
229 :
230 511004 : desc%dh = pw_grid%dh
231 511004 : desc%dh_inv = pw_grid%dh_inv
232 39308 : desc%orthorhombic = pw_grid%orthorhombic
233 39308 : desc%ref_count = 1
234 :
235 39308 : 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 35080 : desc%group_size = pw_grid%para%group%num_pe
259 140320 : desc%npts = pw_grid%npts
260 140320 : desc%ngpts = PRODUCT(INT(desc%npts, KIND=int_8))
261 140320 : desc%lb = pw_grid%bounds(1, :)
262 140320 : desc%ub = pw_grid%bounds(2, :)
263 :
264 : ! this is the eventual border size
265 35080 : IF (border_size == 0) THEN
266 34952 : nmin = (input_settings%nsmax + 1)/2
267 34952 : 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 35080 : IF (input_settings%distribution_type == rsgrid_replicated) THEN
274 :
275 54280 : n_slices = 1
276 13570 : 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 86040 : n_slices = 1
286 21510 : ratio_best = -HUGE(ratio_best)
287 :
288 : ! don't allow distributions with more processors than real grid points
289 64530 : DO k = 1, MIN(desc%npts(3), desc%group_size)
290 150570 : DO j = 1, MIN(desc%npts(2), desc%group_size)
291 86040 : i = MIN(desc%npts(1), desc%group_size/(j*k))
292 344160 : n_slices_tmp = [i, j, k]
293 :
294 : ! we don't match the actual number of CPUs
295 344160 : 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 516498 : IF (.NOT. ALL(PACK(n_slices_tmp == input_settings%distribution_layout, &
300 : [-1, -1, -1] /= input_settings%distribution_layout) &
301 64530 : )) 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 64304 : overlap = .FALSE.
306 257216 : DO dir = 1, 3
307 257216 : IF (n_slices_tmp(dir) > 1) THEN
308 192912 : DO l = 0, n_slices_tmp(dir) - 1
309 128608 : lb = get_limit(desc%npts(dir), n_slices_tmp(dir), l)
310 192912 : IF (lb(2) - lb(1) + 1 + 2*nmin > desc%npts(dir)) overlap = .TRUE.
311 : END DO
312 : END IF
313 : END DO
314 64304 : 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 67536 : MERGE([0.0, 0.0, 0.0], 2*[1.06*nmin, 1.05*nmin, 1.03*nmin], n_slices_tmp == [1, 1, 1]))
323 52668 : IF (ratio > ratio_best) THEN
324 9622 : ratio_best = ratio
325 9622 : 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 21510 : IF (input_settings%distribution_type == rsgrid_automatic) THEN
334 83840 : volume = PRODUCT(REAL(desc%npts, KIND=dp))
335 : volume_dist = PRODUCT(REAL(desc%npts, KIND=dp)/n_slices + &
336 83840 : MERGE([0, 0, 0], 2*[nmin, nmin, nmin], n_slices == [1, 1, 1]))
337 20960 : IF (volume < volume_dist*input_settings%memory_factor) THEN
338 83840 : n_slices = 1
339 : END IF
340 : END IF
341 :
342 : END IF
343 :
344 140320 : desc%group_dim(:) = n_slices(:)
345 35080 : CALL desc%group%from_dup(pw_grid%para%group)
346 35080 : desc%group_size = desc%group%num_pe
347 35080 : desc%my_pos = desc%group%mepos
348 :
349 140142 : 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 34922 : desc%border = border_size
353 34922 : IF (border_size == 0) THEN
354 139688 : desc%perd = 1
355 : ELSE
356 0 : desc%perd = 0
357 : END IF
358 34922 : desc%distributed = .FALSE.
359 34922 : desc%parallel = .TRUE.
360 139688 : desc%group_coor(:) = 0
361 34922 : desc%my_virtual_pos = 0
362 :
363 104766 : ALLOCATE (desc%virtual2real(0:desc%group_size - 1))
364 104766 : ALLOCATE (desc%real2virtual(0:desc%group_size - 1))
365 : ! Start with no reordering
366 104766 : DO i = 0, desc%group_size - 1
367 69844 : desc%virtual2real(i) = i
368 104766 : 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 39308 : CALL timestop(handle)
469 :
470 39308 : END SUBROUTINE rs_grid_create_descriptor
471 :
472 : ! **************************************************************************************************
473 : !> \brief ...
474 : !> \param rs ...
475 : !> \param desc ...
476 : ! **************************************************************************************************
477 5900286 : 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 280966 : CALL timeset(routineN, handle)
487 :
488 280966 : rs%desc => desc
489 280966 : CALL rs_grid_retain_descriptor(rs%desc)
490 :
491 280966 : 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 1121468 : 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 1116760 : rs%lb_real = desc%lb
506 1116760 : rs%ub_real = desc%ub
507 1116760 : rs%lb_local = rs%lb_real - desc%border*(1 - desc%perd)
508 1116760 : rs%ub_local = rs%ub_real + desc%border*(1 - desc%perd)
509 1116760 : rs%npts_local = rs%ub_local - rs%lb_local + 1
510 1116760 : 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 280966 : 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 280966 : rs%lb_local(3):rs%ub_local(3)) => rs%buffer%host_buffer
526 :
527 842898 : ALLOCATE (rs%px(desc%npts(1)))
528 842898 : ALLOCATE (rs%py(desc%npts(2)))
529 842898 : ALLOCATE (rs%pz(desc%npts(3)))
530 :
531 280966 : CALL timestop(handle)
532 :
533 280966 : 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 16940 : 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 16940 : IF (rs%desc%parallel) THEN
586 16648 : IF (iounit > 0) THEN
587 : WRITE (iounit, '(/,A,T71,I10)') &
588 7539 : " RS_GRID| Information for grid number ", rs%desc%pw%id_nr
589 30156 : DO i = 1, 3
590 22617 : WRITE (iounit, '(A,I3,T30,2I8,T62,A,T71,I10)') " RS_GRID| Bounds ", &
591 52773 : i, rs%desc%lb(i), rs%desc%ub(i), "Points:", rs%desc%npts(i)
592 : END DO
593 7539 : IF (.NOT. rs%desc%distributed) THEN
594 7524 : WRITE (iounit, '(A)') " RS_GRID| Real space fully replicated"
595 : WRITE (iounit, '(A,T71,I10)') &
596 7524 : " 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 16648 : 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 16940 : END SUBROUTINE rs_grid_print
645 :
646 : ! **************************************************************************************************
647 : !> \brief ...
648 : !> \param rs ...
649 : !> \param pw ...
650 : ! **************************************************************************************************
651 1277334 : 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 1277334 : CALL timeset(routineN, handle2)
660 1277334 : CALL timeset(routineN//"_"//TRIM(ADJUSTL(cp_to_string(CEILING(pw%pw_grid%cutoff/10)*10))), handle)
661 :
662 1277334 : IF (.NOT. ASSOCIATED(rs%desc%pw, pw%pw_grid)) THEN
663 0 : CPABORT("Different rs and pw indentifiers")
664 : END IF
665 :
666 1277334 : IF (rs%desc%distributed) THEN
667 1852 : CALL transfer_rs2pw_distributed(rs, pw)
668 1275482 : ELSE IF (rs%desc%parallel) THEN
669 1059620 : CALL transfer_rs2pw_replicated(rs, pw)
670 : ELSE ! treat simple serial case locally
671 215862 : IF (rs%desc%border == 0) THEN
672 863448 : 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 1277334 : CALL timestop(handle)
685 1277334 : CALL timestop(handle2)
686 :
687 1277334 : END SUBROUTINE transfer_rs2pw
688 :
689 : ! **************************************************************************************************
690 : !> \brief ...
691 : !> \param rs ...
692 : !> \param pw ...
693 : ! **************************************************************************************************
694 1210005 : 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 1210005 : CALL timeset(routineN, handle2)
704 1210005 : CALL timeset(routineN//"_"//TRIM(ADJUSTL(cp_to_string(CEILING(pw%pw_grid%cutoff/10)*10))), handle)
705 :
706 1210005 : IF (.NOT. ASSOCIATED(rs%desc%pw, pw%pw_grid)) THEN
707 0 : CPABORT("Different rs and pw indentifiers")
708 : END IF
709 :
710 1210005 : IF (rs%desc%distributed) THEN
711 872 : CALL transfer_pw2rs_distributed(rs, pw)
712 1209133 : ELSE IF (rs%desc%parallel) THEN
713 956358 : 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 1210005 : CALL timestop(handle)
754 1210005 : CALL timestop(handle2)
755 :
756 1210005 : 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 1059620 : 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 1059620 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount
770 : INTEGER, DIMENSION(3) :: lb, ub
771 1059620 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: recvbuf, sendbuf, swaparray
772 :
773 : ASSOCIATE (np => pw%pw_grid%para%group%num_pe, bo => pw%pw_grid%para%bo(1:2, 1:3, 0:pw%pw_grid%para%group%num_pe - 1, 1), &
774 : pbo => pw%pw_grid%bounds, group => pw%pw_grid%para%group, mepos => pw%pw_grid%para%group%mepos, &
775 : grid => rs%r)
776 3178860 : ALLOCATE (rcount(0:np - 1))
777 3178860 : DO ip = 1, np
778 9536580 : rcount(ip - 1) = PRODUCT(bo(2, :, ip) - bo(1, :, ip) + 1)
779 : END DO
780 3178860 : nma = MAXVAL(rcount(0:np - 1))
781 4238480 : ALLOCATE (sendbuf(nma), recvbuf(nma))
782 38721679900 : sendbuf = 1.0E99_dp; recvbuf = 1.0E99_dp ! init mpi'ed buffers to silence warnings under valgrind
783 :
784 : !sample peak memory
785 1059620 : CALL m_memory()
786 :
787 1059620 : dest = MODULO(mepos + 1, np)
788 1059620 : source = MODULO(mepos - 1, np)
789 1059620 : sendbuf = 0.0_dp
790 :
791 2119240 : DO ip = 1, np
792 :
793 8476960 : lb = pbo(1, :) + bo(1, :, MODULO(mepos - ip, np) + 1) - 1
794 8476960 : ub = pbo(1, :) + bo(2, :, MODULO(mepos - ip, np) + 1) - 1
795 : ! this loop takes about the same time as the message passing call
796 : ! notice that the range of ix is only a small fraction of the first index of grid
797 : ! therefore it seems faster to have the second index as the innermost loop
798 : ! if this runs on many cpus
799 : ! tested on itanium, pentium4, opteron, ultrasparc...
800 8476960 : s = ub - lb + 1
801 52974652 : DO iz = lb(3), ub(3)
802 924041862 : DO ix = lb(1), ub(1)
803 871067210 : ii = (iz - lb(3))*s(1)*s(2) + (ix - lb(1)) + 1
804 39371298898 : DO iy = lb(2), ub(2)
805 38449376276 : sendbuf(ii) = sendbuf(ii) + grid(ix, iy, iz)
806 39320443486 : ii = ii + s(1)
807 : END DO
808 : END DO
809 : END DO
810 2119240 : IF (ip == np) EXIT
811 1059620 : CALL group%sendrecv(sendbuf, dest, recvbuf, source, 13)
812 1059620 : CALL MOVE_ALLOC(sendbuf, swaparray)
813 1059620 : CALL MOVE_ALLOC(recvbuf, sendbuf)
814 2119240 : CALL MOVE_ALLOC(swaparray, recvbuf)
815 : END DO
816 1059620 : nn = rcount(mepos)
817 : END ASSOCIATE
818 :
819 1059620 : CALL dcopy(nn, sendbuf, 1, pw%array, 1)
820 :
821 1059620 : DEALLOCATE (rcount)
822 1059620 : DEALLOCATE (sendbuf)
823 1059620 : DEALLOCATE (recvbuf)
824 :
825 1059620 : END SUBROUTINE transfer_rs2pw_replicated
826 :
827 : ! **************************************************************************************************
828 : !> \brief transfer from a planewave grid to a realspace grid
829 : !> \param rs ...
830 : !> \param pw ...
831 : ! **************************************************************************************************
832 956358 : SUBROUTINE transfer_pw2rs_replicated(rs, pw)
833 : TYPE(realspace_grid_type), INTENT(IN) :: rs
834 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
835 :
836 : INTEGER :: dest, i, ii, im, ip, ix, iy, iz, j, jm, &
837 : k, km, nma, nn, source
838 956358 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount
839 : INTEGER, DIMENSION(3) :: lb, ub
840 956358 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: recvbuf, sendbuf, swaparray
841 2869074 : TYPE(mp_request_type), DIMENSION(2) :: req
842 :
843 : ASSOCIATE (np => pw%pw_grid%para%group%num_pe, bo => pw%pw_grid%para%bo(1:2, 1:3, 0:pw%pw_grid%para%group%num_pe - 1, 1), &
844 : pbo => pw%pw_grid%bounds, group => pw%pw_grid%para%group, mepos => pw%pw_grid%para%group%mepos, &
845 : grid => rs%r)
846 2869074 : ALLOCATE (rcount(0:np - 1))
847 2869074 : DO ip = 1, np
848 8607222 : rcount(ip - 1) = PRODUCT(bo(2, :, ip) - bo(1, :, ip) + 1)
849 : END DO
850 2869074 : nma = MAXVAL(rcount(0:np - 1))
851 3825432 : ALLOCATE (sendbuf(nma), recvbuf(nma))
852 35812513400 : sendbuf = 1.0E99_dp; recvbuf = 1.0E99_dp ! init mpi'ed buffers to silence warnings under valgrind
853 :
854 : !sample peak memory
855 956358 : CALL m_memory()
856 :
857 956358 : nn = rcount(mepos)
858 956358 : CALL dcopy(nn, pw%array, 1, sendbuf, 1)
859 :
860 956358 : dest = MODULO(mepos + 1, np)
861 956358 : source = MODULO(mepos - 1, np)
862 :
863 2869074 : DO ip = 0, np - 1
864 : ! we must shift the buffer only np-1 times around
865 1912716 : IF (ip /= np - 1) THEN
866 : CALL group%isendrecv(sendbuf, dest, recvbuf, source, &
867 956358 : req(1), req(2), 13)
868 : END IF
869 7650864 : lb = pbo(1, :) + bo(1, :, MODULO(mepos - ip, np) + 1) - 1
870 7650864 : ub = pbo(1, :) + bo(2, :, MODULO(mepos - ip, np) + 1) - 1
871 1912716 : ii = 0
872 : ! this loop takes about the same time as the message passing call
873 : ! If I read the code correctly then:
874 48334500 : DO iz = lb(3), ub(3)
875 1641531764 : DO iy = lb(2), ub(2)
876 37204245900 : DO ix = lb(1), ub(1)
877 35564626852 : ii = ii + 1
878 37157824116 : grid(ix, iy, iz) = sendbuf(ii)
879 : END DO
880 : END DO
881 : END DO
882 1912716 : IF (ip /= np - 1) THEN
883 956358 : CALL mp_waitall(req)
884 : END IF
885 1912716 : CALL MOVE_ALLOC(sendbuf, swaparray)
886 1912716 : CALL MOVE_ALLOC(recvbuf, sendbuf)
887 2869074 : CALL MOVE_ALLOC(swaparray, recvbuf)
888 : END DO
889 1912716 : IF (rs%desc%border > 0) THEN
890 : !$OMP PARALLEL DO DEFAULT(NONE) &
891 : !$OMP PRIVATE(i,im,j,jm,k,km) &
892 0 : !$OMP SHARED(rs)
893 : DO k = rs%lb_local(3), rs%ub_local(3)
894 : IF (k < rs%lb_real(3)) THEN
895 : km = k + rs%desc%npts(3)
896 : ELSE IF (k > rs%ub_real(3)) THEN
897 : km = k - rs%desc%npts(3)
898 : ELSE
899 : km = k
900 : END IF
901 : DO j = rs%lb_local(2), rs%ub_local(2)
902 : IF (j < rs%lb_real(2)) THEN
903 : jm = j + rs%desc%npts(2)
904 : ELSE IF (j > rs%ub_real(2)) THEN
905 : jm = j - rs%desc%npts(2)
906 : ELSE
907 : jm = j
908 : END IF
909 : DO i = rs%lb_local(1), rs%ub_local(1)
910 : IF (i < rs%lb_real(1)) THEN
911 : im = i + rs%desc%npts(1)
912 : ELSE IF (i > rs%ub_real(1)) THEN
913 : im = i - rs%desc%npts(1)
914 : ELSE
915 : im = i
916 : END IF
917 : rs%r(i, j, k) = rs%r(im, jm, km)
918 : END DO
919 : END DO
920 : END DO
921 : !$OMP END PARALLEL DO
922 : END IF
923 : END ASSOCIATE
924 :
925 956358 : DEALLOCATE (rcount)
926 956358 : DEALLOCATE (sendbuf)
927 956358 : DEALLOCATE (recvbuf)
928 :
929 956358 : END SUBROUTINE transfer_pw2rs_replicated
930 :
931 : ! **************************************************************************************************
932 : !> \brief does the rs2pw transfer in the case where the rs grid is
933 : !> distributed (3D domain decomposition)
934 : !> \param rs ...
935 : !> \param pw ...
936 : !> \par History
937 : !> 12.2007 created [Matt Watkins]
938 : !> 9.2008 reduced amount of halo data sent [Iain Bethune]
939 : !> 10.2008 added non-blocking communication [Iain Bethune]
940 : !> 4.2009 added support for rank-reordering on the grid [Iain Bethune]
941 : !> 12.2009 added OMP and sparse alltoall [Iain Bethune]
942 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2008-2009 on behalf of the HECToR project
943 : !> \note
944 : !> the transfer is a two step procedure. For example, for the rs2pw transfer:
945 : !>
946 : !> 1) Halo-exchange in 3D so that the local part of the rs_grid contains the full data
947 : !> 2) an alltoall communication to redistribute the local rs_grid to the local pw_grid
948 : !>
949 : !> the halo exchange is most expensive on a large number of CPUs. Particular in this halo
950 : !> exchange is that the border region is rather large (e.g. 20 points) and that it might overlap
951 : !> with the central domain of several CPUs (i.e. next nearest neighbors)
952 : ! **************************************************************************************************
953 1852 : SUBROUTINE transfer_rs2pw_distributed(rs, pw)
954 : TYPE(realspace_grid_type), INTENT(IN) :: rs
955 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
956 :
957 : CHARACTER(LEN=200) :: error_string
958 : INTEGER :: completed, dest_down, dest_up, i, idir, j, k, lb, my_id, my_pw_rank, my_rs_rank, &
959 : n_shifts, nn, num_threads, position, source_down, source_up, ub, x, y, z
960 1852 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dshifts, recv_disps, recv_sizes, &
961 1852 : send_disps, send_sizes, ushifts
962 3704 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: bounds, recv_tasks, send_tasks
963 : INTEGER, DIMENSION(2) :: neighbours, pos
964 : INTEGER, DIMENSION(3) :: coords, lb_recv, lb_recv_down, lb_recv_up, lb_send, lb_send_down, &
965 : lb_send_up, ub_recv, ub_recv_down, ub_recv_up, ub_send, ub_send_down, ub_send_up
966 : LOGICAL, DIMENSION(3) :: halo_swapped
967 : REAL(KIND=dp) :: pw_sum, rs_sum
968 1852 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: recv_buf_3d_down, recv_buf_3d_up, &
969 1852 : send_buf_3d_down, send_buf_3d_up
970 3704 : TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: recv_bufs, send_bufs
971 1852 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_reqs, send_reqs
972 9260 : TYPE(mp_request_type), DIMENSION(4) :: req
973 :
974 1852 : num_threads = 1
975 1852 : my_id = 0
976 :
977 : ! safety check, to be removed once we're absolute sure the routine is correct
978 : IF (debug_this_module) THEN
979 : rs_sum = accurate_sum(rs%r)*ABS(det_3x3(rs%desc%dh))
980 : CALL rs%desc%group%sum(rs_sum)
981 : END IF
982 :
983 1852 : halo_swapped = .FALSE.
984 : ! We don't need to send the 'edges' of the halos that have already been sent
985 : ! Halos are contiguous in memory in z-direction only, so swap these first,
986 : ! and send less data in the y and x directions which are more expensive
987 :
988 7408 : DO idir = 3, 1, -1
989 :
990 5556 : IF (rs%desc%perd(idir) /= 1) THEN
991 :
992 14412 : ALLOCATE (dshifts(0:rs%desc%neighbours(idir)))
993 9608 : ALLOCATE (ushifts(0:rs%desc%neighbours(idir)))
994 :
995 4804 : ushifts = 0
996 4804 : dshifts = 0
997 :
998 : ! check that we don't try to send data to ourself
999 6656 : DO n_shifts = 1, MIN(rs%desc%neighbours(idir), rs%desc%group_dim(idir) - 1)
1000 :
1001 : ! need to take into account the possible varying widths of neighbouring cells
1002 : ! offset_up and offset_down hold the real size of the neighbouring cells
1003 1852 : position = MODULO(rs%desc%virtual_group_coor(idir) - n_shifts, rs%desc%group_dim(idir))
1004 1852 : neighbours = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), position)
1005 1852 : dshifts(n_shifts) = dshifts(n_shifts - 1) + (neighbours(2) - neighbours(1) + 1)
1006 :
1007 1852 : position = MODULO(rs%desc%virtual_group_coor(idir) + n_shifts, rs%desc%group_dim(idir))
1008 1852 : neighbours = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), position)
1009 1852 : ushifts(n_shifts) = ushifts(n_shifts - 1) + (neighbours(2) - neighbours(1) + 1)
1010 :
1011 : ! The border data has to be send/received from the neighbours
1012 : ! First we calculate the source and destination processes for the shift
1013 : ! We do both shifts at once to allow for more overlap of communication and buffer packing/unpacking
1014 :
1015 1852 : CALL cart_shift(rs, idir, -1*n_shifts, source_down, dest_down)
1016 :
1017 7408 : lb_send_down(:) = rs%lb_local(:)
1018 7408 : lb_recv_down(:) = rs%lb_local(:)
1019 7408 : ub_recv_down(:) = rs%ub_local(:)
1020 7408 : ub_send_down(:) = rs%ub_local(:)
1021 :
1022 1852 : IF (dshifts(n_shifts - 1) <= rs%desc%border) THEN
1023 1852 : ub_send_down(idir) = lb_send_down(idir) + rs%desc%border - 1 - dshifts(n_shifts - 1)
1024 : lb_send_down(idir) = MAX(lb_send_down(idir), &
1025 1852 : lb_send_down(idir) + rs%desc%border - dshifts(n_shifts))
1026 :
1027 1852 : ub_recv_down(idir) = ub_recv_down(idir) - rs%desc%border
1028 : lb_recv_down(idir) = MAX(lb_recv_down(idir) + rs%desc%border, &
1029 1852 : ub_recv_down(idir) - rs%desc%border + 1 + ushifts(n_shifts - 1))
1030 : ELSE
1031 0 : lb_send_down(idir) = 0
1032 0 : ub_send_down(idir) = -1
1033 0 : lb_recv_down(idir) = 0
1034 0 : ub_recv_down(idir) = -1
1035 : END IF
1036 :
1037 7408 : DO i = 1, 3
1038 7408 : IF (halo_swapped(i)) THEN
1039 554 : lb_send_down(i) = rs%lb_real(i)
1040 554 : ub_send_down(i) = rs%ub_real(i)
1041 554 : lb_recv_down(i) = rs%lb_real(i)
1042 554 : ub_recv_down(i) = rs%ub_real(i)
1043 : END IF
1044 : END DO
1045 :
1046 : ! post the receive
1047 0 : ALLOCATE (recv_buf_3d_down(lb_recv_down(1):ub_recv_down(1), &
1048 9260 : lb_recv_down(2):ub_recv_down(2), lb_recv_down(3):ub_recv_down(3)))
1049 1852 : CALL rs%desc%group%irecv(recv_buf_3d_down, source_down, req(1))
1050 :
1051 : ! now allocate, pack and send the send buffer
1052 7408 : nn = PRODUCT(ub_send_down - lb_send_down + 1)
1053 0 : ALLOCATE (send_buf_3d_down(lb_send_down(1):ub_send_down(1), &
1054 9260 : lb_send_down(2):ub_send_down(2), lb_send_down(3):ub_send_down(3)))
1055 :
1056 : !$OMP PARALLEL DEFAULT(NONE), &
1057 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1058 1852 : !$OMP SHARED(send_buf_3d_down,rs,lb_send_down,ub_send_down)
1059 : !$ num_threads = MIN(omp_get_max_threads(), ub_send_down(3) - lb_send_down(3) + 1)
1060 : !$ my_id = omp_get_thread_num()
1061 : IF (my_id < num_threads) THEN
1062 : lb = lb_send_down(3) + ((ub_send_down(3) - lb_send_down(3) + 1)*my_id)/num_threads
1063 : ub = lb_send_down(3) + ((ub_send_down(3) - lb_send_down(3) + 1)*(my_id + 1))/num_threads - 1
1064 :
1065 : send_buf_3d_down(lb_send_down(1):ub_send_down(1), lb_send_down(2):ub_send_down(2), &
1066 : lb:ub) = rs%r(lb_send_down(1):ub_send_down(1), &
1067 : lb_send_down(2):ub_send_down(2), lb:ub)
1068 : END IF
1069 : !$OMP END PARALLEL
1070 :
1071 1852 : CALL rs%desc%group%isend(send_buf_3d_down, dest_down, req(3))
1072 :
1073 : ! Now for the other direction
1074 1852 : CALL cart_shift(rs, idir, n_shifts, source_up, dest_up)
1075 :
1076 7408 : lb_send_up(:) = rs%lb_local(:)
1077 7408 : lb_recv_up(:) = rs%lb_local(:)
1078 7408 : ub_recv_up(:) = rs%ub_local(:)
1079 7408 : ub_send_up(:) = rs%ub_local(:)
1080 :
1081 1852 : IF (ushifts(n_shifts - 1) <= rs%desc%border) THEN
1082 :
1083 1852 : lb_send_up(idir) = ub_send_up(idir) - rs%desc%border + 1 + ushifts(n_shifts - 1)
1084 : ub_send_up(idir) = MIN(ub_send_up(idir), &
1085 1852 : ub_send_up(idir) - rs%desc%border + ushifts(n_shifts))
1086 :
1087 1852 : lb_recv_up(idir) = lb_recv_up(idir) + rs%desc%border
1088 : ub_recv_up(idir) = MIN(ub_recv_up(idir) - rs%desc%border, &
1089 1852 : lb_recv_up(idir) + rs%desc%border - 1 - dshifts(n_shifts - 1))
1090 : ELSE
1091 0 : lb_send_up(idir) = 0
1092 0 : ub_send_up(idir) = -1
1093 0 : lb_recv_up(idir) = 0
1094 0 : ub_recv_up(idir) = -1
1095 : END IF
1096 :
1097 7408 : DO i = 1, 3
1098 7408 : IF (halo_swapped(i)) THEN
1099 554 : lb_send_up(i) = rs%lb_real(i)
1100 554 : ub_send_up(i) = rs%ub_real(i)
1101 554 : lb_recv_up(i) = rs%lb_real(i)
1102 554 : ub_recv_up(i) = rs%ub_real(i)
1103 : END IF
1104 : END DO
1105 :
1106 : ! post the receive
1107 0 : ALLOCATE (recv_buf_3d_up(lb_recv_up(1):ub_recv_up(1), &
1108 9260 : lb_recv_up(2):ub_recv_up(2), lb_recv_up(3):ub_recv_up(3)))
1109 1852 : CALL rs%desc%group%irecv(recv_buf_3d_up, source_up, req(2))
1110 :
1111 : ! now allocate,pack and send the send buffer
1112 7408 : nn = PRODUCT(ub_send_up - lb_send_up + 1)
1113 0 : ALLOCATE (send_buf_3d_up(lb_send_up(1):ub_send_up(1), &
1114 9260 : lb_send_up(2):ub_send_up(2), lb_send_up(3):ub_send_up(3)))
1115 :
1116 : !$OMP PARALLEL DEFAULT(NONE), &
1117 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1118 1852 : !$OMP SHARED(send_buf_3d_up,rs,lb_send_up,ub_send_up)
1119 : !$ num_threads = MIN(omp_get_max_threads(), ub_send_up(3) - lb_send_up(3) + 1)
1120 : !$ my_id = omp_get_thread_num()
1121 : IF (my_id < num_threads) THEN
1122 : lb = lb_send_up(3) + ((ub_send_up(3) - lb_send_up(3) + 1)*my_id)/num_threads
1123 : ub = lb_send_up(3) + ((ub_send_up(3) - lb_send_up(3) + 1)*(my_id + 1))/num_threads - 1
1124 :
1125 : send_buf_3d_up(lb_send_up(1):ub_send_up(1), lb_send_up(2):ub_send_up(2), &
1126 : lb:ub) = rs%r(lb_send_up(1):ub_send_up(1), &
1127 : lb_send_up(2):ub_send_up(2), lb:ub)
1128 : END IF
1129 : !$OMP END PARALLEL
1130 :
1131 1852 : CALL rs%desc%group%isend(send_buf_3d_up, dest_up, req(4))
1132 :
1133 : ! wait for a recv to complete, then we can unpack
1134 :
1135 5556 : DO i = 1, 2
1136 :
1137 3704 : CALL mp_waitany(req(1:2), completed)
1138 :
1139 5556 : IF (completed == 1) THEN
1140 :
1141 : ! only some procs may need later shifts
1142 1852 : IF (ub_recv_down(idir) >= lb_recv_down(idir)) THEN
1143 : ! Sum the data in the RS Grid
1144 : !$OMP PARALLEL DEFAULT(NONE), &
1145 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1146 1852 : !$OMP SHARED(recv_buf_3d_down,rs,lb_recv_down,ub_recv_down)
1147 : !$ num_threads = MIN(omp_get_max_threads(), ub_recv_down(3) - lb_recv_down(3) + 1)
1148 : !$ my_id = omp_get_thread_num()
1149 : IF (my_id < num_threads) THEN
1150 : lb = lb_recv_down(3) + ((ub_recv_down(3) - lb_recv_down(3) + 1)*my_id)/num_threads
1151 : ub = lb_recv_down(3) + ((ub_recv_down(3) - lb_recv_down(3) + 1)*(my_id + 1))/num_threads - 1
1152 :
1153 : rs%r(lb_recv_down(1):ub_recv_down(1), &
1154 : lb_recv_down(2):ub_recv_down(2), lb:ub) = &
1155 : rs%r(lb_recv_down(1):ub_recv_down(1), &
1156 : lb_recv_down(2):ub_recv_down(2), lb:ub) + &
1157 : recv_buf_3d_down(:, :, lb:ub)
1158 : END IF
1159 : !$OMP END PARALLEL
1160 : END IF
1161 1852 : DEALLOCATE (recv_buf_3d_down)
1162 : ELSE
1163 :
1164 : ! only some procs may need later shifts
1165 1852 : IF (ub_recv_up(idir) >= lb_recv_up(idir)) THEN
1166 : ! Sum the data in the RS Grid
1167 : !$OMP PARALLEL DEFAULT(NONE), &
1168 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1169 1852 : !$OMP SHARED(recv_buf_3d_up,rs,lb_recv_up,ub_recv_up)
1170 : !$ num_threads = MIN(omp_get_max_threads(), ub_recv_up(3) - lb_recv_up(3) + 1)
1171 : !$ my_id = omp_get_thread_num()
1172 : IF (my_id < num_threads) THEN
1173 : lb = lb_recv_up(3) + ((ub_recv_up(3) - lb_recv_up(3) + 1)*my_id)/num_threads
1174 : ub = lb_recv_up(3) + ((ub_recv_up(3) - lb_recv_up(3) + 1)*(my_id + 1))/num_threads - 1
1175 :
1176 : rs%r(lb_recv_up(1):ub_recv_up(1), &
1177 : lb_recv_up(2):ub_recv_up(2), lb:ub) = &
1178 : rs%r(lb_recv_up(1):ub_recv_up(1), &
1179 : lb_recv_up(2):ub_recv_up(2), lb:ub) + &
1180 : recv_buf_3d_up(:, :, lb:ub)
1181 : END IF
1182 : !$OMP END PARALLEL
1183 : END IF
1184 1852 : DEALLOCATE (recv_buf_3d_up)
1185 : END IF
1186 :
1187 : END DO
1188 :
1189 : ! make sure the sends have completed before we deallocate
1190 :
1191 1852 : CALL mp_waitall(req(3:4))
1192 :
1193 1852 : DEALLOCATE (send_buf_3d_down)
1194 8508 : DEALLOCATE (send_buf_3d_up)
1195 : END DO
1196 :
1197 4804 : DEALLOCATE (dshifts)
1198 4804 : DEALLOCATE (ushifts)
1199 :
1200 : END IF
1201 :
1202 7408 : halo_swapped(idir) = .TRUE.
1203 :
1204 : END DO
1205 :
1206 : ! This is the real redistribution
1207 7408 : ALLOCATE (bounds(0:pw%pw_grid%para%group%num_pe - 1, 1:4))
1208 :
1209 : ! work out the pw grid points each proc holds
1210 5556 : DO i = 0, pw%pw_grid%para%group%num_pe - 1
1211 11112 : bounds(i, 1:2) = pw%pw_grid%para%bo(1:2, 1, i, 1)
1212 11112 : bounds(i, 3:4) = pw%pw_grid%para%bo(1:2, 2, i, 1)
1213 11112 : bounds(i, 1:2) = bounds(i, 1:2) - pw%pw_grid%npts(1)/2 - 1
1214 12964 : bounds(i, 3:4) = bounds(i, 3:4) - pw%pw_grid%npts(2)/2 - 1
1215 : END DO
1216 :
1217 7408 : ALLOCATE (send_tasks(0:pw%pw_grid%para%group%num_pe - 1, 1:6))
1218 5556 : ALLOCATE (send_sizes(0:pw%pw_grid%para%group%num_pe - 1))
1219 3704 : ALLOCATE (send_disps(0:pw%pw_grid%para%group%num_pe - 1))
1220 3704 : ALLOCATE (recv_tasks(0:pw%pw_grid%para%group%num_pe - 1, 1:6))
1221 3704 : ALLOCATE (recv_sizes(0:pw%pw_grid%para%group%num_pe - 1))
1222 3704 : ALLOCATE (recv_disps(0:pw%pw_grid%para%group%num_pe - 1))
1223 5556 : send_tasks(:, 1) = 1
1224 5556 : send_tasks(:, 2) = 0
1225 5556 : send_tasks(:, 3) = 1
1226 5556 : send_tasks(:, 4) = 0
1227 5556 : send_tasks(:, 5) = 1
1228 5556 : send_tasks(:, 6) = 0
1229 1852 : send_sizes = 0
1230 1852 : recv_sizes = 0
1231 :
1232 1852 : my_rs_rank = rs%desc%my_pos
1233 1852 : my_pw_rank = pw%pw_grid%para%group%mepos
1234 :
1235 : ! find the processors that should hold our data
1236 : ! should be part of the rs grid type
1237 : ! this is a loop over real ranks (i.e. the in-order cartesian ranks)
1238 : ! do the recv and send tasks in two separate loops which will
1239 : ! load balance better for OpenMP with large numbers of MPI tasks
1240 :
1241 : !$OMP PARALLEL DO DEFAULT(NONE), &
1242 : !$OMP PRIVATE(coords,idir,pos,lb_send,ub_send), &
1243 1852 : !$OMP SHARED(rs,bounds,my_rs_rank,recv_tasks,recv_sizes)
1244 : DO i = 0, rs%desc%group_size - 1
1245 :
1246 : coords(:) = rs%desc%rank2coord(:, rs%desc%real2virtual(i))
1247 : !calculate the rs grid points on each processor
1248 : !coords is the part of the grid that rank i actually holds
1249 : DO idir = 1, 3
1250 : pos(:) = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), coords(idir))
1251 : pos(:) = pos(:) - rs%desc%npts(idir)/2 - 1
1252 : lb_send(idir) = pos(1)
1253 : ub_send(idir) = pos(2)
1254 : END DO
1255 :
1256 : IF (lb_send(1) > bounds(my_rs_rank, 2)) CYCLE
1257 : IF (ub_send(1) < bounds(my_rs_rank, 1)) CYCLE
1258 : IF (lb_send(2) > bounds(my_rs_rank, 4)) CYCLE
1259 : IF (ub_send(2) < bounds(my_rs_rank, 3)) CYCLE
1260 :
1261 : recv_tasks(i, 1) = MAX(lb_send(1), bounds(my_rs_rank, 1))
1262 : recv_tasks(i, 2) = MIN(ub_send(1), bounds(my_rs_rank, 2))
1263 : recv_tasks(i, 3) = MAX(lb_send(2), bounds(my_rs_rank, 3))
1264 : recv_tasks(i, 4) = MIN(ub_send(2), bounds(my_rs_rank, 4))
1265 : recv_tasks(i, 5) = lb_send(3)
1266 : recv_tasks(i, 6) = ub_send(3)
1267 : recv_sizes(i) = (recv_tasks(i, 2) - recv_tasks(i, 1) + 1)* &
1268 : (recv_tasks(i, 4) - recv_tasks(i, 3) + 1)*(recv_tasks(i, 6) - recv_tasks(i, 5) + 1)
1269 :
1270 : END DO
1271 : !$OMP END PARALLEL DO
1272 :
1273 7408 : coords(:) = rs%desc%rank2coord(:, rs%desc%real2virtual(my_rs_rank))
1274 7408 : DO idir = 1, 3
1275 5556 : pos(:) = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), coords(idir))
1276 16668 : pos(:) = pos(:) - rs%desc%npts(idir)/2 - 1
1277 5556 : lb_send(idir) = pos(1)
1278 7408 : ub_send(idir) = pos(2)
1279 : END DO
1280 :
1281 1852 : lb_recv(:) = lb_send(:)
1282 1852 : ub_recv(:) = ub_send(:)
1283 : !$OMP PARALLEL DO DEFAULT(NONE), &
1284 1852 : !$OMP SHARED(pw,lb_send,ub_send,bounds,send_tasks,send_sizes)
1285 : DO j = 0, pw%pw_grid%para%group%num_pe - 1
1286 :
1287 : IF (lb_send(1) > bounds(j, 2)) CYCLE
1288 : IF (ub_send(1) < bounds(j, 1)) CYCLE
1289 : IF (lb_send(2) > bounds(j, 4)) CYCLE
1290 : IF (ub_send(2) < bounds(j, 3)) CYCLE
1291 :
1292 : send_tasks(j, 1) = MAX(lb_send(1), bounds(j, 1))
1293 : send_tasks(j, 2) = MIN(ub_send(1), bounds(j, 2))
1294 : send_tasks(j, 3) = MAX(lb_send(2), bounds(j, 3))
1295 : send_tasks(j, 4) = MIN(ub_send(2), bounds(j, 4))
1296 : send_tasks(j, 5) = lb_send(3)
1297 : send_tasks(j, 6) = ub_send(3)
1298 : send_sizes(j) = (send_tasks(j, 2) - send_tasks(j, 1) + 1)* &
1299 : (send_tasks(j, 4) - send_tasks(j, 3) + 1)*(send_tasks(j, 6) - send_tasks(j, 5) + 1)
1300 :
1301 : END DO
1302 : !$OMP END PARALLEL DO
1303 :
1304 1852 : send_disps(0) = 0
1305 1852 : recv_disps(0) = 0
1306 3704 : DO i = 1, pw%pw_grid%para%group%num_pe - 1
1307 1852 : send_disps(i) = send_disps(i - 1) + send_sizes(i - 1)
1308 3704 : recv_disps(i) = recv_disps(i - 1) + recv_sizes(i - 1)
1309 : END DO
1310 :
1311 12964 : CPASSERT(SUM(send_sizes) == PRODUCT(ub_recv - lb_recv + 1))
1312 :
1313 9260 : ALLOCATE (send_bufs(0:rs%desc%group_size - 1))
1314 11112 : ALLOCATE (recv_bufs(0:rs%desc%group_size - 1))
1315 :
1316 5556 : DO i = 0, rs%desc%group_size - 1
1317 3704 : IF (send_sizes(i) /= 0) THEN
1318 10374 : ALLOCATE (send_bufs(i)%array(send_sizes(i)))
1319 : ELSE
1320 246 : NULLIFY (send_bufs(i)%array)
1321 : END IF
1322 5556 : IF (recv_sizes(i) /= 0) THEN
1323 10374 : ALLOCATE (recv_bufs(i)%array(recv_sizes(i)))
1324 : ELSE
1325 246 : NULLIFY (recv_bufs(i)%array)
1326 : END IF
1327 : END DO
1328 :
1329 9260 : ALLOCATE (recv_reqs(0:rs%desc%group_size - 1))
1330 5556 : recv_reqs = mp_request_null
1331 :
1332 5556 : DO i = 0, rs%desc%group_size - 1
1333 5556 : IF (recv_sizes(i) /= 0) THEN
1334 3458 : CALL rs%desc%group%irecv(recv_bufs(i)%array, i, recv_reqs(i))
1335 : END IF
1336 : END DO
1337 :
1338 : ! do packing
1339 : !$OMP PARALLEL DO DEFAULT(NONE), &
1340 : !$OMP PRIVATE(k,z,y,x), &
1341 1852 : !$OMP SHARED(rs,send_tasks,send_bufs,send_disps)
1342 : DO i = 0, rs%desc%group_size - 1
1343 : k = 0
1344 : DO z = send_tasks(i, 5), send_tasks(i, 6)
1345 : DO y = send_tasks(i, 3), send_tasks(i, 4)
1346 : DO x = send_tasks(i, 1), send_tasks(i, 2)
1347 : k = k + 1
1348 : send_bufs(i)%array(k) = rs%r(x, y, z)
1349 : END DO
1350 : END DO
1351 : END DO
1352 : END DO
1353 : !$OMP END PARALLEL DO
1354 :
1355 9260 : ALLOCATE (send_reqs(0:rs%desc%group_size - 1))
1356 5556 : send_reqs = mp_request_null
1357 :
1358 5556 : DO i = 0, rs%desc%group_size - 1
1359 5556 : IF (send_sizes(i) /= 0) THEN
1360 3458 : CALL rs%desc%group%isend(send_bufs(i)%array, i, send_reqs(i))
1361 : END IF
1362 : END DO
1363 :
1364 : ! do unpacking
1365 : ! no OMP here so we can unpack each message as it arrives
1366 5556 : DO i = 0, rs%desc%group_size - 1
1367 3704 : IF (recv_sizes(i) == 0) CYCLE
1368 :
1369 3458 : CALL mp_waitany(recv_reqs, completed)
1370 3458 : k = 0
1371 143780 : DO z = recv_tasks(completed - 1, 5), recv_tasks(completed - 1, 6)
1372 10754290 : DO y = recv_tasks(completed - 1, 3), recv_tasks(completed - 1, 4)
1373 452245022 : DO x = recv_tasks(completed - 1, 1), recv_tasks(completed - 1, 2)
1374 441494436 : k = k + 1
1375 452106552 : pw%array(x, y, z) = recv_bufs(completed - 1)%array(k)
1376 : END DO
1377 : END DO
1378 : END DO
1379 : END DO
1380 :
1381 1852 : CALL mp_waitall(send_reqs)
1382 :
1383 1852 : DEALLOCATE (recv_reqs)
1384 1852 : DEALLOCATE (send_reqs)
1385 :
1386 5556 : DO i = 0, rs%desc%group_size - 1
1387 3704 : IF (ASSOCIATED(send_bufs(i)%array)) THEN
1388 3458 : DEALLOCATE (send_bufs(i)%array)
1389 : END IF
1390 5556 : IF (ASSOCIATED(recv_bufs(i)%array)) THEN
1391 3458 : DEALLOCATE (recv_bufs(i)%array)
1392 : END IF
1393 : END DO
1394 :
1395 1852 : DEALLOCATE (send_bufs)
1396 1852 : DEALLOCATE (recv_bufs)
1397 1852 : DEALLOCATE (send_tasks)
1398 1852 : DEALLOCATE (send_sizes)
1399 1852 : DEALLOCATE (send_disps)
1400 1852 : DEALLOCATE (recv_tasks)
1401 1852 : DEALLOCATE (recv_sizes)
1402 1852 : DEALLOCATE (recv_disps)
1403 :
1404 : IF (debug_this_module) THEN
1405 : ! safety check, to be removed once we're absolute sure the routine is correct
1406 : pw_sum = pw_integrate_function(pw)
1407 : IF (ABS(pw_sum - rs_sum)/MAX(1.0_dp, ABS(pw_sum), ABS(rs_sum)) > EPSILON(rs_sum)*1000) THEN
1408 : WRITE (error_string, '(A,6(1X,I4.4),3F25.16)') "rs_pw_transfer_distributed", &
1409 : rs%desc%npts, rs%desc%group_dim, pw_sum, rs_sum, ABS(pw_sum - rs_sum)
1410 : CALL cp_abort(__LOCATION__, &
1411 : error_string//" Please report this bug ... quick workaround: use "// &
1412 : "DISTRIBUTION_TYPE REPLICATED")
1413 : END IF
1414 : END IF
1415 :
1416 1852 : END SUBROUTINE transfer_rs2pw_distributed
1417 :
1418 : ! **************************************************************************************************
1419 : !> \brief does the pw2rs transfer in the case where the rs grid is
1420 : !> distributed (3D domain decomposition)
1421 : !> \param rs ...
1422 : !> \param pw ...
1423 : !> \par History
1424 : !> 12.2007 created [Matt Watkins]
1425 : !> 9.2008 reduced amount of halo data sent [Iain Bethune]
1426 : !> 10.2008 added non-blocking communication [Iain Bethune]
1427 : !> 4.2009 added support for rank-reordering on the grid [Iain Bethune]
1428 : !> 12.2009 added OMP and sparse alltoall [Iain Bethune]
1429 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2008-2009 on behalf of the HECToR project
1430 : !> \note
1431 : !> the transfer is a two step procedure. For example, for the rs2pw transfer:
1432 : !>
1433 : !> 1) Halo-exchange in 3D so that the local part of the rs_grid contains the full data
1434 : !> 2) an alltoall communication to redistribute the local rs_grid to the local pw_grid
1435 : !>
1436 : !> the halo exchange is most expensive on a large number of CPUs. Particular in this halo
1437 : !> exchange is that the border region is rather large (e.g. 20 points) and that it might overlap
1438 : !> with the central domain of several CPUs (i.e. next nearest neighbors)
1439 : ! **************************************************************************************************
1440 872 : SUBROUTINE transfer_pw2rs_distributed(rs, pw)
1441 : TYPE(realspace_grid_type), INTENT(IN) :: rs
1442 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
1443 :
1444 : INTEGER :: completed, dest_down, dest_up, i, idir, j, k, lb, my_id, my_pw_rank, my_rs_rank, &
1445 : n_shifts, nn, num_threads, position, source_down, source_up, ub, x, y, z
1446 872 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dshifts, recv_disps, recv_sizes, &
1447 872 : send_disps, send_sizes, ushifts
1448 1744 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: bounds, recv_tasks, send_tasks
1449 : INTEGER, DIMENSION(2) :: neighbours, pos
1450 : INTEGER, DIMENSION(3) :: coords, lb_recv, lb_recv_down, lb_recv_up, lb_send, lb_send_down, &
1451 : lb_send_up, ub_recv, ub_recv_down, ub_recv_up, ub_send, ub_send_down, ub_send_up
1452 : LOGICAL, DIMENSION(3) :: halo_swapped
1453 872 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: recv_buf_3d_down, recv_buf_3d_up, &
1454 872 : send_buf_3d_down, send_buf_3d_up
1455 1744 : TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: recv_bufs, send_bufs
1456 872 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_reqs, send_reqs
1457 4360 : TYPE(mp_request_type), DIMENSION(4) :: req
1458 :
1459 872 : num_threads = 1
1460 872 : my_id = 0
1461 :
1462 872 : CALL rs_grid_zero(rs)
1463 :
1464 : ! This is the real redistribution
1465 :
1466 3488 : ALLOCATE (bounds(0:pw%pw_grid%para%group%num_pe - 1, 1:4))
1467 :
1468 2616 : DO i = 0, pw%pw_grid%para%group%num_pe - 1
1469 5232 : bounds(i, 1:2) = pw%pw_grid%para%bo(1:2, 1, i, 1)
1470 5232 : bounds(i, 3:4) = pw%pw_grid%para%bo(1:2, 2, i, 1)
1471 5232 : bounds(i, 1:2) = bounds(i, 1:2) - pw%pw_grid%npts(1)/2 - 1
1472 6104 : bounds(i, 3:4) = bounds(i, 3:4) - pw%pw_grid%npts(2)/2 - 1
1473 : END DO
1474 :
1475 3488 : ALLOCATE (send_tasks(0:pw%pw_grid%para%group%num_pe - 1, 1:6))
1476 2616 : ALLOCATE (send_sizes(0:pw%pw_grid%para%group%num_pe - 1))
1477 1744 : ALLOCATE (send_disps(0:pw%pw_grid%para%group%num_pe - 1))
1478 1744 : ALLOCATE (recv_tasks(0:pw%pw_grid%para%group%num_pe - 1, 1:6))
1479 1744 : ALLOCATE (recv_sizes(0:pw%pw_grid%para%group%num_pe - 1))
1480 1744 : ALLOCATE (recv_disps(0:pw%pw_grid%para%group%num_pe - 1))
1481 :
1482 872 : send_tasks = 0
1483 2616 : send_tasks(:, 1) = 1
1484 2616 : send_tasks(:, 2) = 0
1485 2616 : send_tasks(:, 3) = 1
1486 2616 : send_tasks(:, 4) = 0
1487 2616 : send_tasks(:, 5) = 1
1488 2616 : send_tasks(:, 6) = 0
1489 872 : send_sizes = 0
1490 :
1491 872 : recv_tasks = 0
1492 2616 : recv_tasks(:, 1) = 1
1493 2616 : recv_tasks(:, 2) = 0
1494 2616 : send_tasks(:, 3) = 1
1495 2616 : send_tasks(:, 4) = 0
1496 2616 : send_tasks(:, 5) = 1
1497 2616 : send_tasks(:, 6) = 0
1498 872 : recv_sizes = 0
1499 :
1500 872 : my_rs_rank = rs%desc%my_pos
1501 872 : my_pw_rank = pw%pw_grid%para%group%mepos
1502 :
1503 : ! find the processors that should hold our data
1504 : ! should be part of the rs grid type
1505 : ! this is a loop over real ranks (i.e. the in-order cartesian ranks)
1506 : ! do the recv and send tasks in two separate loops which will
1507 : ! load balance better for OpenMP with large numbers of MPI tasks
1508 :
1509 : ! this is the reverse of rs2pw: what were the sends are now the recvs
1510 :
1511 : !$OMP PARALLEL DO DEFAULT(NONE), &
1512 : !$OMP PRIVATE(coords,idir,pos,lb_send,ub_send), &
1513 872 : !$OMP SHARED(rs,bounds,my_rs_rank,send_tasks,send_sizes,pw)
1514 : DO i = 0, pw%pw_grid%para%group%num_pe - 1
1515 :
1516 : coords(:) = rs%desc%rank2coord(:, rs%desc%real2virtual(i))
1517 : !calculate the real rs grid points on each processor
1518 : !coords is the part of the grid that rank i actually holds
1519 : DO idir = 1, 3
1520 : pos(:) = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), coords(idir))
1521 : pos(:) = pos(:) - rs%desc%npts(idir)/2 - 1
1522 : lb_send(idir) = pos(1)
1523 : ub_send(idir) = pos(2)
1524 : END DO
1525 :
1526 : IF (ub_send(1) < bounds(my_rs_rank, 1)) CYCLE
1527 : IF (lb_send(1) > bounds(my_rs_rank, 2)) CYCLE
1528 : IF (ub_send(2) < bounds(my_rs_rank, 3)) CYCLE
1529 : IF (lb_send(2) > bounds(my_rs_rank, 4)) CYCLE
1530 :
1531 : send_tasks(i, 1) = MAX(lb_send(1), bounds(my_rs_rank, 1))
1532 : send_tasks(i, 2) = MIN(ub_send(1), bounds(my_rs_rank, 2))
1533 : send_tasks(i, 3) = MAX(lb_send(2), bounds(my_rs_rank, 3))
1534 : send_tasks(i, 4) = MIN(ub_send(2), bounds(my_rs_rank, 4))
1535 : send_tasks(i, 5) = lb_send(3)
1536 : send_tasks(i, 6) = ub_send(3)
1537 : send_sizes(i) = (send_tasks(i, 2) - send_tasks(i, 1) + 1)* &
1538 : (send_tasks(i, 4) - send_tasks(i, 3) + 1)*(send_tasks(i, 6) - send_tasks(i, 5) + 1)
1539 :
1540 : END DO
1541 : !$OMP END PARALLEL DO
1542 :
1543 3488 : coords(:) = rs%desc%rank2coord(:, rs%desc%real2virtual(my_rs_rank))
1544 3488 : DO idir = 1, 3
1545 2616 : pos(:) = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), coords(idir))
1546 7848 : pos(:) = pos(:) - rs%desc%npts(idir)/2 - 1
1547 2616 : lb_send(idir) = pos(1)
1548 3488 : ub_send(idir) = pos(2)
1549 : END DO
1550 :
1551 872 : lb_recv(:) = lb_send(:)
1552 872 : ub_recv(:) = ub_send(:)
1553 :
1554 : !$OMP PARALLEL DO DEFAULT(NONE), &
1555 872 : !$OMP SHARED(pw,lb_send,ub_send,bounds,recv_tasks,recv_sizes)
1556 : DO j = 0, pw%pw_grid%para%group%num_pe - 1
1557 :
1558 : IF (ub_send(1) < bounds(j, 1)) CYCLE
1559 : IF (lb_send(1) > bounds(j, 2)) CYCLE
1560 : IF (ub_send(2) < bounds(j, 3)) CYCLE
1561 : IF (lb_send(2) > bounds(j, 4)) CYCLE
1562 :
1563 : recv_tasks(j, 1) = MAX(lb_send(1), bounds(j, 1))
1564 : recv_tasks(j, 2) = MIN(ub_send(1), bounds(j, 2))
1565 : recv_tasks(j, 3) = MAX(lb_send(2), bounds(j, 3))
1566 : recv_tasks(j, 4) = MIN(ub_send(2), bounds(j, 4))
1567 : recv_tasks(j, 5) = lb_send(3)
1568 : recv_tasks(j, 6) = ub_send(3)
1569 : recv_sizes(j) = (recv_tasks(j, 2) - recv_tasks(j, 1) + 1)* &
1570 : (recv_tasks(j, 4) - recv_tasks(j, 3) + 1)*(recv_tasks(j, 6) - recv_tasks(j, 5) + 1)
1571 :
1572 : END DO
1573 : !$OMP END PARALLEL DO
1574 :
1575 872 : send_disps(0) = 0
1576 872 : recv_disps(0) = 0
1577 1744 : DO i = 1, pw%pw_grid%para%group%num_pe - 1
1578 872 : send_disps(i) = send_disps(i - 1) + send_sizes(i - 1)
1579 1744 : recv_disps(i) = recv_disps(i - 1) + recv_sizes(i - 1)
1580 : END DO
1581 :
1582 6104 : CPASSERT(SUM(recv_sizes) == PRODUCT(ub_recv - lb_recv + 1))
1583 :
1584 4360 : ALLOCATE (send_bufs(0:rs%desc%group_size - 1))
1585 5232 : ALLOCATE (recv_bufs(0:rs%desc%group_size - 1))
1586 :
1587 2616 : DO i = 0, rs%desc%group_size - 1
1588 1744 : IF (send_sizes(i) /= 0) THEN
1589 4986 : ALLOCATE (send_bufs(i)%array(send_sizes(i)))
1590 : ELSE
1591 82 : NULLIFY (send_bufs(i)%array)
1592 : END IF
1593 2616 : IF (recv_sizes(i) /= 0) THEN
1594 4986 : ALLOCATE (recv_bufs(i)%array(recv_sizes(i)))
1595 : ELSE
1596 82 : NULLIFY (recv_bufs(i)%array)
1597 : END IF
1598 : END DO
1599 :
1600 4360 : ALLOCATE (recv_reqs(0:rs%desc%group_size - 1))
1601 2616 : recv_reqs = mp_request_null
1602 :
1603 2616 : DO i = 0, rs%desc%group_size - 1
1604 2616 : IF (recv_sizes(i) /= 0) THEN
1605 1662 : CALL rs%desc%group%irecv(recv_bufs(i)%array, i, recv_reqs(i))
1606 : END IF
1607 : END DO
1608 :
1609 : ! do packing
1610 : !$OMP PARALLEL DO DEFAULT(NONE), &
1611 : !$OMP PRIVATE(k,z,y,x), &
1612 872 : !$OMP SHARED(pw,rs,send_tasks,send_bufs,send_disps)
1613 : DO i = 0, rs%desc%group_size - 1
1614 : k = 0
1615 : DO z = send_tasks(i, 5), send_tasks(i, 6)
1616 : DO y = send_tasks(i, 3), send_tasks(i, 4)
1617 : DO x = send_tasks(i, 1), send_tasks(i, 2)
1618 : k = k + 1
1619 : send_bufs(i)%array(k) = pw%array(x, y, z)
1620 : END DO
1621 : END DO
1622 : END DO
1623 : END DO
1624 : !$OMP END PARALLEL DO
1625 :
1626 4360 : ALLOCATE (send_reqs(0:rs%desc%group_size - 1))
1627 2616 : send_reqs = mp_request_null
1628 :
1629 2616 : DO i = 0, rs%desc%group_size - 1
1630 2616 : IF (send_sizes(i) /= 0) THEN
1631 1662 : CALL rs%desc%group%isend(send_bufs(i)%array, i, send_reqs(i))
1632 : END IF
1633 : END DO
1634 :
1635 : ! do unpacking
1636 : ! no OMP here so we can unpack each message as it arrives
1637 :
1638 2616 : DO i = 0, rs%desc%group_size - 1
1639 1744 : IF (recv_sizes(i) == 0) CYCLE
1640 :
1641 1662 : CALL mp_waitany(recv_reqs, completed)
1642 1662 : k = 0
1643 66556 : DO z = recv_tasks(completed - 1, 5), recv_tasks(completed - 1, 6)
1644 4763300 : DO y = recv_tasks(completed - 1, 3), recv_tasks(completed - 1, 4)
1645 194725231 : DO x = recv_tasks(completed - 1, 1), recv_tasks(completed - 1, 2)
1646 189963675 : k = k + 1
1647 194661209 : rs%r(x, y, z) = recv_bufs(completed - 1)%array(k)
1648 : END DO
1649 : END DO
1650 : END DO
1651 : END DO
1652 :
1653 872 : CALL mp_waitall(send_reqs)
1654 :
1655 872 : DEALLOCATE (recv_reqs)
1656 872 : DEALLOCATE (send_reqs)
1657 :
1658 2616 : DO i = 0, rs%desc%group_size - 1
1659 1744 : IF (ASSOCIATED(send_bufs(i)%array)) THEN
1660 1662 : DEALLOCATE (send_bufs(i)%array)
1661 : END IF
1662 2616 : IF (ASSOCIATED(recv_bufs(i)%array)) THEN
1663 1662 : DEALLOCATE (recv_bufs(i)%array)
1664 : END IF
1665 : END DO
1666 :
1667 872 : DEALLOCATE (send_bufs)
1668 872 : DEALLOCATE (recv_bufs)
1669 872 : DEALLOCATE (send_tasks)
1670 872 : DEALLOCATE (send_sizes)
1671 872 : DEALLOCATE (send_disps)
1672 872 : DEALLOCATE (recv_tasks)
1673 872 : DEALLOCATE (recv_sizes)
1674 872 : DEALLOCATE (recv_disps)
1675 :
1676 : ! now pass wings around
1677 872 : halo_swapped = .FALSE.
1678 :
1679 3488 : DO idir = 1, 3
1680 :
1681 2616 : IF (rs%desc%perd(idir) /= 1) THEN
1682 :
1683 5568 : ALLOCATE (dshifts(0:rs%desc%neighbours(idir)))
1684 3712 : ALLOCATE (ushifts(0:rs%desc%neighbours(idir)))
1685 1856 : ushifts = 0
1686 1856 : dshifts = 0
1687 :
1688 3712 : DO n_shifts = 1, rs%desc%neighbours(idir)
1689 :
1690 : ! need to take into account the possible varying widths of neighbouring cells
1691 : ! ushifts and dshifts hold the real size of the neighbouring cells
1692 :
1693 1856 : position = MODULO(rs%desc%virtual_group_coor(idir) - n_shifts, rs%desc%group_dim(idir))
1694 1856 : neighbours = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), position)
1695 1856 : dshifts(n_shifts) = dshifts(n_shifts - 1) + (neighbours(2) - neighbours(1) + 1)
1696 :
1697 1856 : position = MODULO(rs%desc%virtual_group_coor(idir) + n_shifts, rs%desc%group_dim(idir))
1698 1856 : neighbours = get_limit(rs%desc%npts(idir), rs%desc%group_dim(idir), position)
1699 1856 : ushifts(n_shifts) = ushifts(n_shifts - 1) + (neighbours(2) - neighbours(1) + 1)
1700 :
1701 : ! The border data has to be send/received from the neighbors
1702 : ! First we calculate the source and destination processes for the shift
1703 : ! The first shift is "downwards"
1704 :
1705 1856 : CALL cart_shift(rs, idir, -1*n_shifts, source_down, dest_down)
1706 :
1707 7424 : lb_send_down(:) = rs%lb_local(:)
1708 7424 : ub_send_down(:) = rs%ub_local(:)
1709 7424 : lb_recv_down(:) = rs%lb_local(:)
1710 7424 : ub_recv_down(:) = rs%ub_local(:)
1711 :
1712 1856 : IF (dshifts(n_shifts - 1) <= rs%desc%border) THEN
1713 1856 : lb_send_down(idir) = lb_send_down(idir) + rs%desc%border
1714 : ub_send_down(idir) = MIN(ub_send_down(idir) - rs%desc%border, &
1715 1856 : lb_send_down(idir) + rs%desc%border - 1 - dshifts(n_shifts - 1))
1716 :
1717 1856 : lb_recv_down(idir) = ub_recv_down(idir) - rs%desc%border + 1 + ushifts(n_shifts - 1)
1718 : ub_recv_down(idir) = MIN(ub_recv_down(idir), &
1719 1856 : ub_recv_down(idir) - rs%desc%border + ushifts(n_shifts))
1720 : ELSE
1721 0 : lb_send_down(idir) = 0
1722 0 : ub_send_down(idir) = -1
1723 0 : lb_recv_down(idir) = 0
1724 0 : ub_recv_down(idir) = -1
1725 : END IF
1726 :
1727 7424 : DO i = 1, 3
1728 7424 : IF (.NOT. (halo_swapped(i) .OR. i == idir)) THEN
1729 1560 : lb_send_down(i) = rs%lb_real(i)
1730 1560 : ub_send_down(i) = rs%ub_real(i)
1731 1560 : lb_recv_down(i) = rs%lb_real(i)
1732 1560 : ub_recv_down(i) = rs%ub_real(i)
1733 : END IF
1734 : END DO
1735 :
1736 : ! allocate the recv buffer
1737 7424 : nn = PRODUCT(ub_recv_down - lb_recv_down + 1)
1738 0 : ALLOCATE (recv_buf_3d_down(lb_recv_down(1):ub_recv_down(1), &
1739 9280 : lb_recv_down(2):ub_recv_down(2), lb_recv_down(3):ub_recv_down(3)))
1740 :
1741 : ! recv buffer is now ready, so post the receive
1742 1856 : CALL rs%desc%group%irecv(recv_buf_3d_down, source_down, req(1))
1743 :
1744 : ! now allocate,pack and send the send buffer
1745 7424 : nn = PRODUCT(ub_send_down - lb_send_down + 1)
1746 0 : ALLOCATE (send_buf_3d_down(lb_send_down(1):ub_send_down(1), &
1747 9280 : lb_send_down(2):ub_send_down(2), lb_send_down(3):ub_send_down(3)))
1748 :
1749 : !$OMP PARALLEL DEFAULT(NONE), &
1750 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1751 1856 : !$OMP SHARED(send_buf_3d_down,rs,lb_send_down,ub_send_down)
1752 : !$ num_threads = MIN(omp_get_max_threads(), ub_send_down(3) - lb_send_down(3) + 1)
1753 : !$ my_id = omp_get_thread_num()
1754 : IF (my_id < num_threads) THEN
1755 : lb = lb_send_down(3) + ((ub_send_down(3) - lb_send_down(3) + 1)*my_id)/num_threads
1756 : ub = lb_send_down(3) + ((ub_send_down(3) - lb_send_down(3) + 1)*(my_id + 1))/num_threads - 1
1757 :
1758 : send_buf_3d_down(lb_send_down(1):ub_send_down(1), lb_send_down(2):ub_send_down(2), &
1759 : lb:ub) = rs%r(lb_send_down(1):ub_send_down(1), &
1760 : lb_send_down(2):ub_send_down(2), lb:ub)
1761 : END IF
1762 : !$OMP END PARALLEL
1763 :
1764 1856 : CALL rs%desc%group%isend(send_buf_3d_down, dest_down, req(3))
1765 :
1766 : ! Now for the other direction
1767 :
1768 1856 : CALL cart_shift(rs, idir, n_shifts, source_up, dest_up)
1769 :
1770 7424 : lb_send_up(:) = rs%lb_local(:)
1771 7424 : ub_send_up(:) = rs%ub_local(:)
1772 7424 : lb_recv_up(:) = rs%lb_local(:)
1773 7424 : ub_recv_up(:) = rs%ub_local(:)
1774 :
1775 1856 : IF (ushifts(n_shifts - 1) <= rs%desc%border) THEN
1776 1856 : ub_send_up(idir) = ub_send_up(idir) - rs%desc%border
1777 : lb_send_up(idir) = MAX(lb_send_up(idir) + rs%desc%border, &
1778 1856 : ub_send_up(idir) - rs%desc%border + 1 + ushifts(n_shifts - 1))
1779 :
1780 1856 : ub_recv_up(idir) = lb_recv_up(idir) + rs%desc%border - 1 - dshifts(n_shifts - 1)
1781 : lb_recv_up(idir) = MAX(lb_recv_up(idir), &
1782 1856 : lb_recv_up(idir) + rs%desc%border - dshifts(n_shifts))
1783 : ELSE
1784 0 : lb_send_up(idir) = 0
1785 0 : ub_send_up(idir) = -1
1786 0 : lb_recv_up(idir) = 0
1787 0 : ub_recv_up(idir) = -1
1788 : END IF
1789 :
1790 7424 : DO i = 1, 3
1791 7424 : IF (.NOT. (halo_swapped(i) .OR. i == idir)) THEN
1792 1560 : lb_send_up(i) = rs%lb_real(i)
1793 1560 : ub_send_up(i) = rs%ub_real(i)
1794 1560 : lb_recv_up(i) = rs%lb_real(i)
1795 1560 : ub_recv_up(i) = rs%ub_real(i)
1796 : END IF
1797 : END DO
1798 :
1799 : ! allocate the recv buffer
1800 7424 : nn = PRODUCT(ub_recv_up - lb_recv_up + 1)
1801 0 : ALLOCATE (recv_buf_3d_up(lb_recv_up(1):ub_recv_up(1), &
1802 9280 : lb_recv_up(2):ub_recv_up(2), lb_recv_up(3):ub_recv_up(3)))
1803 :
1804 : ! recv buffer is now ready, so post the receive
1805 :
1806 1856 : CALL rs%desc%group%irecv(recv_buf_3d_up, source_up, req(2))
1807 :
1808 : ! now allocate,pack and send the send buffer
1809 7424 : nn = PRODUCT(ub_send_up - lb_send_up + 1)
1810 0 : ALLOCATE (send_buf_3d_up(lb_send_up(1):ub_send_up(1), &
1811 9280 : lb_send_up(2):ub_send_up(2), lb_send_up(3):ub_send_up(3)))
1812 :
1813 : !$OMP PARALLEL DEFAULT(NONE), &
1814 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1815 1856 : !$OMP SHARED(send_buf_3d_up,rs,lb_send_up,ub_send_up)
1816 : !$ num_threads = MIN(omp_get_max_threads(), ub_send_up(3) - lb_send_up(3) + 1)
1817 : !$ my_id = omp_get_thread_num()
1818 : IF (my_id < num_threads) THEN
1819 : lb = lb_send_up(3) + ((ub_send_up(3) - lb_send_up(3) + 1)*my_id)/num_threads
1820 : ub = lb_send_up(3) + ((ub_send_up(3) - lb_send_up(3) + 1)*(my_id + 1))/num_threads - 1
1821 :
1822 : send_buf_3d_up(lb_send_up(1):ub_send_up(1), lb_send_up(2):ub_send_up(2), &
1823 : lb:ub) = rs%r(lb_send_up(1):ub_send_up(1), &
1824 : lb_send_up(2):ub_send_up(2), lb:ub)
1825 : END IF
1826 : !$OMP END PARALLEL
1827 :
1828 1856 : CALL rs%desc%group%isend(send_buf_3d_up, dest_up, req(4))
1829 :
1830 : ! wait for a recv to complete, then we can unpack
1831 :
1832 5568 : DO i = 1, 2
1833 :
1834 3712 : CALL mp_waitany(req(1:2), completed)
1835 :
1836 5568 : IF (completed == 1) THEN
1837 :
1838 : ! only some procs may need later shifts
1839 1856 : IF (ub_recv_down(idir) >= lb_recv_down(idir)) THEN
1840 :
1841 : ! Add the data to the RS Grid
1842 : !$OMP PARALLEL DEFAULT(NONE), &
1843 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1844 1856 : !$OMP SHARED(recv_buf_3d_down,rs,lb_recv_down,ub_recv_down)
1845 : !$ num_threads = MIN(omp_get_max_threads(), ub_recv_down(3) - lb_recv_down(3) + 1)
1846 : !$ my_id = omp_get_thread_num()
1847 : IF (my_id < num_threads) THEN
1848 : lb = lb_recv_down(3) + ((ub_recv_down(3) - lb_recv_down(3) + 1)*my_id)/num_threads
1849 : ub = lb_recv_down(3) + ((ub_recv_down(3) - lb_recv_down(3) + 1)*(my_id + 1))/num_threads - 1
1850 :
1851 : rs%r(lb_recv_down(1):ub_recv_down(1), lb_recv_down(2):ub_recv_down(2), &
1852 : lb:ub) = recv_buf_3d_down(:, :, lb:ub)
1853 : END IF
1854 : !$OMP END PARALLEL
1855 : END IF
1856 :
1857 1856 : DEALLOCATE (recv_buf_3d_down)
1858 : ELSE
1859 :
1860 : ! only some procs may need later shifts
1861 1856 : IF (ub_recv_up(idir) >= lb_recv_up(idir)) THEN
1862 :
1863 : ! Add the data to the RS Grid
1864 : !$OMP PARALLEL DEFAULT(NONE), &
1865 : !$OMP PRIVATE(lb,ub,my_id,NUM_THREADS), &
1866 1856 : !$OMP SHARED(recv_buf_3d_up,rs,lb_recv_up,ub_recv_up)
1867 : !$ num_threads = MIN(omp_get_max_threads(), ub_recv_up(3) - lb_recv_up(3) + 1)
1868 : !$ my_id = omp_get_thread_num()
1869 : IF (my_id < num_threads) THEN
1870 : lb = lb_recv_up(3) + ((ub_recv_up(3) - lb_recv_up(3) + 1)*my_id)/num_threads
1871 : ub = lb_recv_up(3) + ((ub_recv_up(3) - lb_recv_up(3) + 1)*(my_id + 1))/num_threads - 1
1872 :
1873 : rs%r(lb_recv_up(1):ub_recv_up(1), lb_recv_up(2):ub_recv_up(2), &
1874 : lb:ub) = recv_buf_3d_up(:, :, lb:ub)
1875 : END IF
1876 : !$OMP END PARALLEL
1877 : END IF
1878 :
1879 1856 : DEALLOCATE (recv_buf_3d_up)
1880 : END IF
1881 : END DO
1882 :
1883 1856 : CALL mp_waitall(req(3:4))
1884 :
1885 1856 : DEALLOCATE (send_buf_3d_down)
1886 5568 : DEALLOCATE (send_buf_3d_up)
1887 : END DO
1888 :
1889 1856 : DEALLOCATE (ushifts)
1890 1856 : DEALLOCATE (dshifts)
1891 : END IF
1892 :
1893 3488 : halo_swapped(idir) = .TRUE.
1894 :
1895 : END DO
1896 :
1897 872 : END SUBROUTINE transfer_pw2rs_distributed
1898 :
1899 : ! **************************************************************************************************
1900 : !> \brief Initialize grid to zero
1901 : !> \param rs ...
1902 : !> \par History
1903 : !> none
1904 : !> \author JGH (23-Mar-2002)
1905 : ! **************************************************************************************************
1906 405692 : SUBROUTINE rs_grid_zero(rs)
1907 :
1908 : TYPE(realspace_grid_type), INTENT(IN) :: rs
1909 :
1910 : CHARACTER(len=*), PARAMETER :: routineN = 'rs_grid_zero'
1911 :
1912 : INTEGER :: handle, i, j, k, l(3), u(3)
1913 :
1914 405692 : CALL timeset(routineN, handle)
1915 1217076 : l(1) = LBOUND(rs%r, 1); l(2) = LBOUND(rs%r, 2); l(3) = LBOUND(rs%r, 3)
1916 1217076 : u(1) = UBOUND(rs%r, 1); u(2) = UBOUND(rs%r, 2); u(3) = UBOUND(rs%r, 3)
1917 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(3) &
1918 : !$OMP PRIVATE(i,j,k) &
1919 405692 : !$OMP SHARED(rs,l,u)
1920 : DO k = l(3), u(3)
1921 : DO j = l(2), u(2)
1922 : DO i = l(1), u(1)
1923 : rs%r(i, j, k) = 0.0_dp
1924 : END DO
1925 : END DO
1926 : END DO
1927 : !$OMP END PARALLEL DO
1928 405692 : CALL timestop(handle)
1929 :
1930 405692 : END SUBROUTINE rs_grid_zero
1931 :
1932 : ! **************************************************************************************************
1933 : !> \brief rs1(i) = rs1(i) + rs2(i)*rs3(i)
1934 : !> \param rs1 ...
1935 : !> \param rs2 ...
1936 : !> \param rs3 ...
1937 : !> \param scalar ...
1938 : !> \par History
1939 : !> none
1940 : !> \author
1941 : ! **************************************************************************************************
1942 1404 : SUBROUTINE rs_grid_mult_and_add(rs1, rs2, rs3, scalar)
1943 :
1944 : TYPE(realspace_grid_type), INTENT(IN) :: rs1, rs2, rs3
1945 : REAL(dp), INTENT(IN) :: scalar
1946 :
1947 : CHARACTER(len=*), PARAMETER :: routineN = 'rs_grid_mult_and_add'
1948 :
1949 : INTEGER :: handle, i, j, k, l(3), u(3)
1950 :
1951 : !-----------------------------------------------------------------------------!
1952 :
1953 1404 : CALL timeset(routineN, handle)
1954 1404 : IF (scalar /= 0.0_dp) THEN
1955 4212 : l(1) = LBOUND(rs1%r, 1); l(2) = LBOUND(rs1%r, 2); l(3) = LBOUND(rs1%r, 3)
1956 4212 : u(1) = UBOUND(rs1%r, 1); u(2) = UBOUND(rs1%r, 2); u(3) = UBOUND(rs1%r, 3)
1957 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(3) &
1958 : !$OMP PRIVATE(i,j,k) &
1959 1404 : !$OMP SHARED(rs1,rs2,rs3,scalar,l,u)
1960 : DO k = l(3), u(3)
1961 : DO j = l(2), u(2)
1962 : DO i = l(1), u(1)
1963 : rs1%r(i, j, k) = rs1%r(i, j, k) + scalar*rs2%r(i, j, k)*rs3%r(i, j, k)
1964 : END DO
1965 : END DO
1966 : END DO
1967 : !$OMP END PARALLEL DO
1968 : END IF
1969 1404 : CALL timestop(handle)
1970 1404 : END SUBROUTINE rs_grid_mult_and_add
1971 :
1972 : ! **************************************************************************************************
1973 : !> \brief Set box matrix info for real space grid
1974 : !> This is needed for variable cell simulations
1975 : !> \param pw_grid ...
1976 : !> \param rs ...
1977 : !> \par History
1978 : !> none
1979 : !> \author JGH (15-May-2007)
1980 : ! **************************************************************************************************
1981 204206 : SUBROUTINE rs_grid_set_box(pw_grid, rs)
1982 :
1983 : TYPE(pw_grid_type), INTENT(IN), TARGET :: pw_grid
1984 : TYPE(realspace_grid_type), INTENT(IN) :: rs
1985 :
1986 204206 : CPASSERT(ASSOCIATED(rs%desc%pw, pw_grid))
1987 2654678 : rs%desc%dh = pw_grid%dh
1988 2654678 : rs%desc%dh_inv = pw_grid%dh_inv
1989 :
1990 204206 : END SUBROUTINE rs_grid_set_box
1991 :
1992 : ! **************************************************************************************************
1993 : !> \brief retains the given rs grid descriptor (see doc/ReferenceCounting.html)
1994 : !> \param rs_desc the grid descriptor to retain
1995 : !> \par History
1996 : !> 04.2009 created [Iain Bethune]
1997 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
1998 : ! **************************************************************************************************
1999 283454 : SUBROUTINE rs_grid_retain_descriptor(rs_desc)
2000 : TYPE(realspace_grid_desc_type), INTENT(INOUT) :: rs_desc
2001 :
2002 283454 : CPASSERT(rs_desc%ref_count > 0)
2003 283454 : rs_desc%ref_count = rs_desc%ref_count + 1
2004 283454 : END SUBROUTINE rs_grid_retain_descriptor
2005 :
2006 : ! **************************************************************************************************
2007 : !> \brief releases the given rs grid (see doc/ReferenceCounting.html)
2008 : !> \param rs_grid the rs grid to release
2009 : !> \par History
2010 : !> 03.2003 created [fawzi]
2011 : !> \author fawzi
2012 : ! **************************************************************************************************
2013 282718 : SUBROUTINE rs_grid_release(rs_grid)
2014 : TYPE(realspace_grid_type), INTENT(INOUT) :: rs_grid
2015 :
2016 282718 : CALL rs_grid_release_descriptor(rs_grid%desc)
2017 :
2018 282718 : CALL offload_free_buffer(rs_grid%buffer)
2019 282718 : NULLIFY (rs_grid%r)
2020 :
2021 282718 : IF (ALLOCATED(rs_grid%px)) DEALLOCATE (rs_grid%px)
2022 282718 : IF (ALLOCATED(rs_grid%py)) DEALLOCATE (rs_grid%py)
2023 282718 : IF (ALLOCATED(rs_grid%pz)) DEALLOCATE (rs_grid%pz)
2024 282718 : END SUBROUTINE rs_grid_release
2025 :
2026 : ! **************************************************************************************************
2027 : !> \brief releases the given rs grid descriptor (see doc/ReferenceCounting.html)
2028 : !> \param rs_desc the rs grid descriptor to release
2029 : !> \par History
2030 : !> 04.2009 created [Iain Bethune]
2031 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
2032 : ! **************************************************************************************************
2033 326389 : SUBROUTINE rs_grid_release_descriptor(rs_desc)
2034 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
2035 :
2036 326389 : IF (ASSOCIATED(rs_desc)) THEN
2037 322762 : CPASSERT(rs_desc%ref_count > 0)
2038 322762 : rs_desc%ref_count = rs_desc%ref_count - 1
2039 322762 : IF (rs_desc%ref_count == 0) THEN
2040 :
2041 39308 : CALL pw_grid_release(rs_desc%pw)
2042 :
2043 39308 : IF (rs_desc%parallel) THEN
2044 : ! release the group communicator
2045 35080 : CALL rs_desc%group%free()
2046 :
2047 35080 : DEALLOCATE (rs_desc%virtual2real)
2048 35080 : DEALLOCATE (rs_desc%real2virtual)
2049 : END IF
2050 :
2051 39308 : IF (rs_desc%distributed) THEN
2052 158 : DEALLOCATE (rs_desc%rank2coord)
2053 158 : DEALLOCATE (rs_desc%coord2rank)
2054 158 : DEALLOCATE (rs_desc%lb_global)
2055 158 : DEALLOCATE (rs_desc%ub_global)
2056 158 : DEALLOCATE (rs_desc%x2coord)
2057 158 : DEALLOCATE (rs_desc%y2coord)
2058 158 : DEALLOCATE (rs_desc%z2coord)
2059 : END IF
2060 :
2061 39308 : DEALLOCATE (rs_desc)
2062 : END IF
2063 : END IF
2064 326389 : NULLIFY (rs_desc)
2065 326389 : END SUBROUTINE rs_grid_release_descriptor
2066 :
2067 : ! **************************************************************************************************
2068 : !> \brief emulates the function of an MPI_cart_shift operation, but the shift is
2069 : !> done in virtual coordinates, and the corresponding real ranks are returned
2070 : !> \param rs_grid ...
2071 : !> \param dir ...
2072 : !> \param disp ...
2073 : !> \param source ...
2074 : !> \param dest ...
2075 : !> \par History
2076 : !> 04.2009 created [Iain Bethune]
2077 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
2078 : ! **************************************************************************************************
2079 7416 : PURE SUBROUTINE cart_shift(rs_grid, dir, disp, source, dest)
2080 :
2081 : TYPE(realspace_grid_type), INTENT(IN) :: rs_grid
2082 : INTEGER, INTENT(IN) :: dir, disp
2083 : INTEGER, INTENT(OUT) :: source, dest
2084 :
2085 : INTEGER, DIMENSION(3) :: shift_coords
2086 :
2087 29664 : shift_coords = rs_grid%desc%virtual_group_coor
2088 7416 : shift_coords(dir) = MODULO(shift_coords(dir) + disp, rs_grid%desc%group_dim(dir))
2089 7416 : dest = rs_grid%desc%virtual2real(rs_grid%desc%coord2rank(shift_coords(1), shift_coords(2), shift_coords(3)))
2090 29664 : shift_coords = rs_grid%desc%virtual_group_coor
2091 7416 : shift_coords(dir) = MODULO(shift_coords(dir) - disp, rs_grid%desc%group_dim(dir))
2092 7416 : source = rs_grid%desc%virtual2real(rs_grid%desc%coord2rank(shift_coords(1), shift_coords(2), shift_coords(3)))
2093 :
2094 7416 : END SUBROUTINE cart_shift
2095 :
2096 : ! **************************************************************************************************
2097 : !> \brief returns the maximum number of points in the local grid of any process
2098 : !> to account for the case where the grid may later be reordered
2099 : !> \param desc ...
2100 : !> \return ...
2101 : !> \par History
2102 : !> 10.2011 created [Iain Bethune]
2103 : ! **************************************************************************************************
2104 0 : FUNCTION rs_grid_max_ngpts(desc) RESULT(max_ngpts)
2105 : TYPE(realspace_grid_desc_type), INTENT(IN) :: desc
2106 : INTEGER :: max_ngpts
2107 :
2108 : CHARACTER(len=*), PARAMETER :: routineN = 'rs_grid_max_ngpts'
2109 :
2110 : INTEGER :: handle, i
2111 : INTEGER, DIMENSION(3) :: lb, ub
2112 :
2113 0 : CALL timeset(routineN, handle)
2114 :
2115 0 : max_ngpts = 0
2116 0 : IF ((desc%pw%para%mode == PW_MODE_LOCAL) .OR. &
2117 : (ALL(desc%group_dim == 1))) THEN
2118 0 : CPASSERT(PRODUCT(INT(desc%npts, KIND=int_8)) < HUGE(1))
2119 0 : max_ngpts = PRODUCT(desc%npts)
2120 : ELSE
2121 0 : DO i = 0, desc%group_size - 1
2122 0 : lb = desc%lb_global(:, i)
2123 0 : ub = desc%ub_global(:, i)
2124 0 : lb = lb - desc%border*(1 - desc%perd)
2125 0 : ub = ub + desc%border*(1 - desc%perd)
2126 0 : CPASSERT(PRODUCT(INT(ub - lb + 1, KIND=int_8)) < HUGE(1))
2127 0 : max_ngpts = MAX(max_ngpts, PRODUCT(ub - lb + 1))
2128 : END DO
2129 : END IF
2130 :
2131 0 : CALL timestop(handle)
2132 :
2133 0 : END FUNCTION rs_grid_max_ngpts
2134 :
2135 : ! **************************************************************************************************
2136 : !> \brief ...
2137 : !> \param rs_grid ...
2138 : !> \param h_inv ...
2139 : !> \param ra ...
2140 : !> \param offset ...
2141 : !> \param group_size ...
2142 : !> \param my_pos ...
2143 : !> \return ...
2144 : ! **************************************************************************************************
2145 1625923 : PURE LOGICAL FUNCTION map_gaussian_here(rs_grid, h_inv, ra, offset, group_size, my_pos) RESULT(res)
2146 : TYPE(realspace_grid_type), INTENT(IN) :: rs_grid
2147 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: h_inv
2148 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra
2149 : INTEGER, INTENT(IN), OPTIONAL :: offset, group_size, my_pos
2150 :
2151 : INTEGER :: dir, lb(3), location(3), tp(3), ub(3)
2152 :
2153 1625923 : res = .FALSE.
2154 :
2155 6503684 : IF (.NOT. ALL(rs_grid%desc%perd == 1)) THEN
2156 32 : DO dir = 1, 3
2157 : ! bounds of local grid (i.e. removing the 'wings'), if periodic
2158 96 : tp(dir) = FLOOR(DOT_PRODUCT(h_inv(dir, :), ra)*rs_grid%desc%npts(dir))
2159 24 : tp(dir) = MODULO(tp(dir), rs_grid%desc%npts(dir))
2160 24 : IF (rs_grid%desc%perd(dir) /= 1) THEN
2161 8 : lb(dir) = rs_grid%lb_local(dir) + rs_grid%desc%border
2162 8 : ub(dir) = rs_grid%ub_local(dir) - rs_grid%desc%border
2163 : ELSE
2164 16 : lb(dir) = rs_grid%lb_local(dir)
2165 16 : ub(dir) = rs_grid%ub_local(dir)
2166 : END IF
2167 : ! distributed grid, only map if it is local to the grid
2168 32 : location(dir) = tp(dir) + rs_grid%desc%lb(dir)
2169 : END DO
2170 60 : IF (ALL(lb(:) <= location(:)) .AND. ALL(location(:) <= ub(:))) THEN
2171 4 : res = .TRUE.
2172 : END IF
2173 : ELSE
2174 1625915 : IF (PRESENT(offset) .AND. PRESENT(group_size) .AND. PRESENT(my_pos)) THEN
2175 : ! not distributed, just a round-robin distribution over the full set of CPUs
2176 1625915 : IF (MODULO(offset, group_size) == my_pos) res = .TRUE.
2177 : END IF
2178 : END IF
2179 :
2180 1625923 : END FUNCTION map_gaussian_here
2181 :
2182 0 : END MODULE realspace_grid_types
|
