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 : !> \par History
10 : !> JGH (30-Nov-2000): ESSL FFT Library added
11 : !> JGH (05-Jan-2001): Added SGI library FFT
12 : !> JGH (14-Jan-2001): Added parallel 3d FFT
13 : !> JGH (10-Feb-2006): New interface type
14 : !> JGH (31-Mar-2008): Remove local allocates and reshapes (performance)
15 : !> Possible problems can be related with setting arrays
16 : !> not to zero
17 : !> Some interfaces could be further simplified by avoiding
18 : !> an initial copy. However, this assumes contiguous arrays
19 : !> IAB (15-Oct-2008): Moved mp_cart_sub calls out of cube_tranpose_* and into
20 : !> fft_scratch type, reducing number of calls dramatically
21 : !> IAB (05-Dec-2008): Moved all other non-essential MPI calls into scratch type
22 : !> IAB (09-Jan-2009): Added fft_plan_type to store FFT data, including cached FFTW plans
23 : !> IAB (13-Feb-2009): Extended plan caching to serial 3D FFT (fft3d_s)
24 : !> IAB (09-Oct-2009): Added OpenMP directives to parallel 3D FFT
25 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2008-2009 on behalf of the HECToR project
26 : !> HFP (17-Oct-2024): Thread-safety insurance (CPASSERT), and OMP ATOMIC (tick_fft_pool)
27 : !> \author JGH
28 : ! **************************************************************************************************
29 : MODULE fft_tools
30 : USE ISO_C_BINDING, ONLY: C_F_POINTER,&
31 : C_LOC,&
32 : C_PTR,&
33 : C_SIZE_T
34 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
35 : USE fft_lib, ONLY: &
36 : fft_1d, fft_3d, fft_alloc, fft_create_plan_1d, fft_create_plan_3d, fft_dealloc, &
37 : fft_destroy_plan, fft_do_cleanup, fft_do_init, fft_get_lengths
38 : USE fft_plan, ONLY: fft_plan_type
39 : USE kinds, ONLY: dp,&
40 : dp_size,&
41 : sp
42 : USE mathconstants, ONLY: z_zero
43 : USE message_passing, ONLY: mp_cart_type,&
44 : mp_comm_null,&
45 : mp_comm_type,&
46 : mp_request_type,&
47 : mp_waitall
48 : USE offload_api, ONLY: offload_free_pinned_mem,&
49 : offload_malloc_pinned_mem
50 :
51 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_in_parallel
52 :
53 : #include "../base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fft_tools'
58 :
59 : ! Types for the pool of scratch data needed in FFT routines
60 : ! keep the subroutine "is_equal" up-to-date
61 : ! needs a default initialization
62 : TYPE fft_scratch_sizes
63 : INTEGER :: nx = 0, ny = 0, nz = 0
64 : INTEGER :: lmax = 0, mmax = 0, nmax = 0
65 : INTEGER :: mx1 = 0, mx2 = 0, mx3 = 0
66 : INTEGER :: my1 = 0, my2 = 0, my3 = 0
67 : INTEGER :: mz1 = 0, mz2 = 0, mz3 = 0
68 : INTEGER :: mcz1 = 0, mcz2 = 0, mcy3 = 0, mcx2 = 0
69 : INTEGER :: lg = 0, mg = 0
70 : INTEGER :: nbx = 0, nbz = 0
71 : INTEGER :: nmray = 0, nyzray = 0
72 : TYPE(mp_cart_type) :: rs_group = mp_cart_type()
73 : INTEGER, DIMENSION(2) :: g_pos = 0, r_pos = 0, r_dim = 0
74 : INTEGER :: numtask = 0
75 : END TYPE fft_scratch_sizes
76 :
77 : TYPE fft_scratch_type
78 : INTEGER :: fft_scratch_id = 0
79 : INTEGER :: tf_type = -1
80 : LOGICAL :: in_use = .TRUE.
81 : TYPE(mp_comm_type) :: group = mp_comm_type()
82 : INTEGER, DIMENSION(3) :: nfft = -1
83 : ! to be used in cube_transpose_* routines
84 : TYPE(mp_cart_type), DIMENSION(2) :: cart_sub_comm = mp_cart_type()
85 : INTEGER, DIMENSION(2) :: dim = -1, pos = -1
86 : ! to be used in fft3d_s
87 : COMPLEX(KIND=dp), DIMENSION(:, :, :), POINTER, CONTIGUOUS &
88 : :: ziptr => NULL(), zoptr => NULL()
89 : ! to be used in fft3d_ps : block distribution
90 : COMPLEX(KIND=dp), DIMENSION(:, :), CONTIGUOUS, POINTER &
91 : :: p1buf => NULL(), p2buf => NULL(), p3buf => NULL(), p4buf => NULL(), &
92 : p5buf => NULL(), p6buf => NULL(), p7buf => NULL()
93 : ! to be used in fft3d_ps : plane distribution
94 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER, CONTIGUOUS &
95 : :: r1buf => NULL(), r2buf => NULL()
96 : COMPLEX(KIND=dp), DIMENSION(:, :, :), POINTER, CONTIGUOUS &
97 : :: tbuf => NULL()
98 : ! to be used in fft3d_pb
99 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER, CONTIGUOUS &
100 : :: a1buf => NULL(), a2buf => NULL(), a3buf => NULL(), &
101 : a4buf => NULL(), a5buf => NULL(), a6buf => NULL()
102 : ! to be used in communication routines
103 : INTEGER, DIMENSION(:), CONTIGUOUS, POINTER :: scount => NULL(), rcount => NULL(), sdispl => NULL(), rdispl => NULL()
104 : INTEGER, DIMENSION(:, :), CONTIGUOUS, POINTER :: pgcube => NULL()
105 : INTEGER, DIMENSION(:), CONTIGUOUS, POINTER :: xzcount => NULL(), yzcount => NULL(), xzdispl => NULL(), yzdispl => NULL()
106 : INTEGER :: in = 0, mip = -1
107 : REAL(KIND=dp) :: rsratio = 1.0_dp
108 : COMPLEX(KIND=dp), DIMENSION(:), POINTER, CONTIGUOUS &
109 : :: xzbuf => NULL(), yzbuf => NULL()
110 : COMPLEX(KIND=sp), DIMENSION(:), POINTER, CONTIGUOUS &
111 : :: xzbuf_sgl => NULL(), yzbuf_sgl => NULL()
112 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER, CONTIGUOUS &
113 : :: rbuf1 => NULL(), rbuf2 => NULL(), rbuf3 => NULL(), rbuf4 => NULL(), &
114 : rbuf5 => NULL(), rbuf6 => NULL(), rr => NULL()
115 : COMPLEX(KIND=sp), DIMENSION(:, :), POINTER, CONTIGUOUS &
116 : :: ss => NULL(), tt => NULL()
117 : INTEGER, DIMENSION(:, :), POINTER, CONTIGUOUS :: pgrid => NULL()
118 : INTEGER, DIMENSION(:), POINTER, CONTIGUOUS :: xcor => NULL(), zcor => NULL(), pzcoord => NULL()
119 : TYPE(fft_scratch_sizes) :: sizes = fft_scratch_sizes()
120 : TYPE(fft_plan_type), DIMENSION(6) :: fft_plan = fft_plan_type()
121 : INTEGER :: last_tick = -1
122 : END TYPE fft_scratch_type
123 :
124 : TYPE fft_scratch_pool_type
125 : TYPE(fft_scratch_type), POINTER :: fft_scratch => NULL()
126 : TYPE(fft_scratch_pool_type), POINTER :: fft_scratch_next => NULL()
127 : END TYPE fft_scratch_pool_type
128 :
129 : INTEGER, SAVE :: init_fft_pool = 0
130 : ! the clock for fft pool. Allows to identify the least recently used scratch
131 : INTEGER, SAVE :: tick_fft_pool = 0
132 : ! limit the number of scratch pools to fft_pool_scratch_limit.
133 : INTEGER, SAVE :: fft_pool_scratch_limit = 15
134 : TYPE(fft_scratch_pool_type), POINTER, SAVE :: fft_scratch_first
135 : ! END of types for the pool of scratch data needed in FFT routines
136 :
137 : PRIVATE
138 : PUBLIC :: init_fft, fft3d, finalize_fft
139 : PUBLIC :: fft_alloc, fft_dealloc
140 : PUBLIC :: init_fft_scratch_pool
141 : PUBLIC :: fft_radix_operations, fft_1d_many
142 : PUBLIC :: FWFFT, BWFFT
143 : PUBLIC :: FFT_RADIX_CLOSEST, FFT_RADIX_NEXT
144 : PUBLIC :: FFT_RADIX_NEXT_ODD
145 :
146 : INTEGER, PARAMETER :: FWFFT = +1, BWFFT = -1
147 : INTEGER, PARAMETER :: FFT_RADIX_CLOSEST = 493, FFT_RADIX_NEXT = 494
148 : INTEGER, PARAMETER :: FFT_RADIX_ALLOWED = 495, FFT_RADIX_DISALLOWED = 496
149 : INTEGER, PARAMETER :: FFT_RADIX_NEXT_ODD = 497
150 :
151 : REAL(KIND=dp), PARAMETER :: ratio_sparse_alltoall = 0.5_dp
152 :
153 : ! these saved variables are FFT globals
154 : LOGICAL, SAVE :: alltoall_sgl = .FALSE.
155 :
156 : ! these are only needed for pw_gpu (-D__OFFLOAD)
157 : PUBLIC :: get_fft_scratch, release_fft_scratch
158 : PUBLIC :: cube_transpose_1, cube_transpose_2
159 : PUBLIC :: yz_to_x, x_to_yz, xz_to_yz, yz_to_xz
160 : PUBLIC :: fft_scratch_sizes, fft_scratch_type
161 :
162 : INTERFACE fft3d
163 : MODULE PROCEDURE fft3d_s, fft3d_ps, fft3d_pb
164 : END INTERFACE
165 :
166 : ! **************************************************************************************************
167 :
168 : CONTAINS
169 :
170 : ! **************************************************************************************************
171 : !> \brief ...
172 : !> \param fftlib ...
173 : !> \param alltoall ...
174 : !> \param pool_limit ...
175 : !> \param wisdom_file ...
176 : !> \param plan_style ...
177 : !> \author JGH
178 : ! **************************************************************************************************
179 11603 : SUBROUTINE init_fft(fftlib, alltoall, pool_limit, wisdom_file, &
180 : plan_style)
181 :
182 : CHARACTER(LEN=*), INTENT(IN) :: fftlib
183 : LOGICAL, INTENT(IN) :: alltoall
184 : INTEGER, INTENT(IN) :: pool_limit
185 : CHARACTER(LEN=*), INTENT(IN) :: wisdom_file
186 : INTEGER, INTENT(IN) :: plan_style
187 :
188 11603 : alltoall_sgl = alltoall
189 11603 : fft_pool_scratch_limit = pool_limit
190 :
191 11603 : CALL release_fft_scratch_pool()
192 :
193 11603 : CALL fft_do_init(fftlib, plan_style, wisdom_file)
194 :
195 : ! setup the FFT scratch pool, if one is associated, clear first
196 11603 : CALL init_fft_scratch_pool()
197 :
198 11603 : END SUBROUTINE init_fft
199 :
200 : ! **************************************************************************************************
201 : !> \brief does whatever is needed to finalize the current fft setup
202 : !> \param para_env ...
203 : !> \param wisdom_file ...
204 : !> \par History
205 : !> 10.2007 created [Joost VandeVondele]
206 : ! **************************************************************************************************
207 11393 : SUBROUTINE finalize_fft(para_env, wisdom_file)
208 : CLASS(mp_comm_type) :: para_env
209 : CHARACTER(LEN=*), INTENT(IN) :: wisdom_file
210 :
211 : ! release the FFT scratch pool
212 :
213 11393 : CALL release_fft_scratch_pool()
214 :
215 : ! finalize fft libs
216 :
217 11393 : CALL fft_do_cleanup(wisdom_file, para_env%is_source())
218 :
219 11393 : END SUBROUTINE finalize_fft
220 :
221 : ! **************************************************************************************************
222 : !> \brief Determine the allowed lengths of FFT's '''
223 : !> \param radix_in ...
224 : !> \param radix_out ...
225 : !> \param operation ...
226 : !> \par History
227 : !> new library structure (JGH)
228 : !> \author Ari Seitsonen
229 : ! **************************************************************************************************
230 269191 : SUBROUTINE fft_radix_operations(radix_in, radix_out, operation)
231 :
232 : INTEGER, INTENT(IN) :: radix_in
233 : INTEGER, INTENT(OUT) :: radix_out
234 : INTEGER, INTENT(IN) :: operation
235 :
236 : INTEGER, PARAMETER :: number_of_small_radices = 81
237 : INTEGER, DIMENSION(number_of_small_radices), PARAMETER :: small_radices = [2, 4, 6, 8, 9, 12,&
238 : 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 54, 60, 64, 72, 75, 80, 81, 90, 96, &
239 : 100, 108, 120, 125, 128, 135, 144, 150, 160, 162, 180, 192, 200, 216, 225, 240, 243, 256, &
240 : 270, 288, 300, 320, 324, 360, 375, 384, 400, 405, 432, 450, 480, 486, 500, 512, 540, 576, &
241 : 600, 625, 640, 648, 675, 720, 729, 750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024]
242 :
243 : INTEGER :: i, iloc, ldata, trial_number
244 269191 : INTEGER, ALLOCATABLE, DIMENSION(:) :: DATA
245 :
246 269191 : ldata = 1024
247 269191 : ALLOCATE (DATA(ldata))
248 275920775 : DATA = -1
249 :
250 22073662 : DATA(1:number_of_small_radices) = small_radices(1:number_of_small_radices)
251 :
252 : ! First run: try smaller radices
253 : ! Second run: Fetch extended list of radices and try these
254 269191 : DO trial_number = 1, 2
255 269191 : IF (trial_number == 2) CALL fft_get_lengths(DATA, ldata)
256 269191 : iloc = 0
257 3465865 : DO i = 1, ldata
258 3465865 : IF (DATA(i) == radix_in) THEN
259 : iloc = i
260 : EXIT
261 : ELSE
262 3320319 : IF (OPERATION == FFT_RADIX_ALLOWED) THEN
263 : CYCLE
264 3320319 : ELSE IF (DATA(i) > radix_in) THEN
265 : iloc = i
266 : EXIT
267 : END IF
268 : END IF
269 : END DO
270 :
271 269191 : IF (iloc == 0) THEN
272 0 : IF (trial_number == 1) THEN
273 : ! We haven't found any suitable size in the list of smaller radices
274 : ! because the requested radix is too large (large cell or large cutoff)
275 : ! So, try again with the larger radices
276 : CYCLE
277 : ELSE
278 0 : IF (OPERATION == FFT_RADIX_CLOSEST) THEN
279 0 : radix_out = MAXVAL(DATA)
280 : ELSE
281 : CALL cp_abort(__LOCATION__, &
282 0 : "No suitable FFT lengths found! Either enable FFTW or reduce the grid size!")
283 : END IF
284 0 : RETURN
285 : END IF
286 :
287 : END IF
288 :
289 269191 : IF (OPERATION == FFT_RADIX_ALLOWED) THEN
290 0 : IF (DATA(iloc) == radix_in) THEN
291 0 : radix_out = FFT_RADIX_ALLOWED
292 : ELSE
293 0 : radix_out = FFT_RADIX_DISALLOWED
294 : END IF
295 :
296 269191 : ELSE IF (OPERATION == FFT_RADIX_CLOSEST) THEN
297 288 : IF (DATA(iloc) == radix_in) THEN
298 120 : radix_out = DATA(iloc)
299 : ELSE
300 168 : IF (ABS(DATA(iloc - 1) - radix_in) <= &
301 : ABS(DATA(iloc) - radix_in)) THEN
302 162 : radix_out = DATA(iloc - 1)
303 : ELSE
304 6 : radix_out = DATA(iloc)
305 : END IF
306 : END IF
307 :
308 268903 : ELSE IF (OPERATION == FFT_RADIX_NEXT) THEN
309 266881 : radix_out = DATA(iloc)
310 :
311 2022 : ELSE IF (OPERATION == FFT_RADIX_NEXT_ODD) THEN
312 3476 : DO i = iloc, ldata
313 3476 : IF (MOD(DATA(i), 2) == 1) THEN
314 2022 : radix_out = DATA(i)
315 2022 : EXIT
316 : END IF
317 : END DO
318 2022 : IF (MOD(radix_out, 2) == 0) THEN
319 : ! No larger odd radix found
320 0 : IF (trial_number == 1) THEN
321 : CYCLE
322 : ELSE
323 0 : radix_out = radix_in + 1 - MOD(radix_in, 2)
324 : END IF
325 : END IF
326 :
327 : END IF
328 0 : EXIT
329 : END DO
330 :
331 269191 : DEALLOCATE (DATA)
332 :
333 0 : END SUBROUTINE fft_radix_operations
334 :
335 : ! **************************************************************************************************
336 : !> \brief Performs m 1-D forward FFT-s of size n.
337 : !> \param sign ...
338 : !> \param n size of the FFT
339 : !> \param m number of FFT-s
340 : !> \param trans_in transposition state of the input array
341 : !> \param trans_out transposition state of the output array
342 : !> \param ldx_in ...
343 : !> \param ldx_out ...
344 : !> \param zin input array
345 : !> \param zout output array
346 : !> \param scale scaling factor
347 : !> \param stat status of the operation, non-zero code indicates an error
348 : ! **************************************************************************************************
349 1884 : SUBROUTINE fft_1d_many(sign, n, m, trans_in, trans_out, ldx_in, ldx_out, zin, zout, scale, stat)
350 : INTEGER, INTENT(in) :: sign, n, m
351 : LOGICAL, INTENT(in) :: trans_in, trans_out
352 : INTEGER, INTENT(in) :: ldx_in, ldx_out
353 : COMPLEX(kind=dp), DIMENSION(*), INTENT(inout) :: zin, zout
354 : REAL(kind=dp), INTENT(in) :: scale
355 : INTEGER, INTENT(out) :: stat
356 :
357 : CHARACTER(len=*), PARAMETER :: routineN = 'fft_1d_many'
358 :
359 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:), &
360 1884 : POINTER :: ztmp
361 : INTEGER :: handle
362 : TYPE(fft_plan_type) :: fft_plan
363 :
364 1884 : CALL timeset(routineN, handle)
365 :
366 1884 : IF (sign /= FWFFT .AND. sign /= BWFFT) THEN
367 0 : CPABORT("Unknown FFT-direction!")
368 : END IF
369 :
370 1884 : NULLIFY (ztmp)
371 1884 : IF (trans_in) THEN
372 0 : CALL fft_alloc(ztmp, [ldx_in*n])
373 : ELSE
374 3768 : CALL fft_alloc(ztmp, [ldx_in*m])
375 : END IF
376 1884 : CALL fft_create_plan_1d(fft_plan, sign, trans_in, trans_out, ldx_in, ldx_out, n, m, ztmp, zout)
377 1884 : CALL fft_dealloc(ztmp)
378 1884 : CALL fft_1d(fft_plan, zin, zout, scale, stat)
379 1884 : CALL fft_destroy_plan(fft_plan)
380 :
381 1884 : CALL timestop(handle)
382 7536 : END SUBROUTINE fft_1d_many
383 :
384 : ! **************************************************************************************************
385 : !> \brief Calls the 3D-FFT function from the initialized library
386 : !> \param fsign ...
387 : !> \param n ...
388 : !> \param zin ...
389 : !> \param zout ...
390 : !> \param status ...
391 : !> \param debug ...
392 : !> \par History
393 : !> none
394 : !> \author JGH
395 : ! **************************************************************************************************
396 678145 : SUBROUTINE fft3d_s(fsign, n, zin, zout, status, debug)
397 :
398 : INTEGER, INTENT(IN) :: fsign
399 : INTEGER, DIMENSION(:), INTENT(INOUT) :: n
400 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
401 : INTENT(INOUT) :: zin
402 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
403 : INTENT(INOUT), OPTIONAL, TARGET :: zout
404 : INTEGER, INTENT(OUT), OPTIONAL :: status
405 : LOGICAL, INTENT(IN), OPTIONAL :: debug
406 :
407 : CHARACTER(len=*), PARAMETER :: routineN = 'fft3d_s'
408 :
409 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
410 678145 : POINTER :: zoptr
411 : COMPLEX(KIND=dp), DIMENSION(1, 1, 1), TARGET :: zdum
412 : INTEGER :: handle, ld(3), lo(3), output_unit, sign, &
413 : stat
414 : LOGICAL :: fft_in_place, test
415 : REAL(KIND=dp) :: in_sum, norm, out_sum
416 : TYPE(fft_scratch_type), POINTER :: fft_scratch
417 :
418 678145 : CALL timeset(routineN, handle)
419 678145 : output_unit = cp_logger_get_default_io_unit()
420 :
421 678145 : IF (fsign == FWFFT) THEN
422 1274016 : norm = 1.0_dp/REAL(PRODUCT(n), KIND=dp)
423 359641 : ELSE IF (fsign == BWFFT) THEN
424 359641 : norm = 1.0_dp
425 : ELSE
426 0 : CPABORT("Unknown FFT direction!")
427 : END IF
428 :
429 678145 : IF (PRESENT(debug)) THEN
430 618424 : test = debug
431 : ELSE
432 : test = .FALSE.
433 : END IF
434 :
435 678145 : IF (PRESENT(zout)) THEN
436 : fft_in_place = .FALSE.
437 : ELSE
438 670443 : fft_in_place = .TRUE.
439 : END IF
440 :
441 678145 : IF (test) THEN
442 0 : in_sum = SUM(ABS(zin))
443 : END IF
444 :
445 678145 : ld(1) = SIZE(zin, 1)
446 678145 : ld(2) = SIZE(zin, 2)
447 678145 : ld(3) = SIZE(zin, 3)
448 :
449 678145 : IF (n(1) /= ld(1) .OR. n(2) /= ld(2) .OR. n(3) /= ld(3)) THEN
450 0 : CPABORT("Size and dimension (zin) have to be the same.")
451 : END IF
452 :
453 678145 : sign = fsign
454 678145 : CALL get_fft_scratch(fft_scratch, tf_type=400, n=n)
455 :
456 678145 : IF (fft_in_place) THEN
457 670443 : zoptr => zdum
458 670443 : IF (fsign == FWFFT) THEN
459 313268 : CALL fft_3d(fft_scratch%fft_plan(1), norm, zin, zoptr, stat)
460 : ELSE
461 357175 : CALL fft_3d(fft_scratch%fft_plan(2), norm, zin, zoptr, stat)
462 : END IF
463 : ELSE
464 7702 : IF (fsign == FWFFT) THEN
465 5236 : CALL fft_3d(fft_scratch%fft_plan(3), norm, zin, zout, stat)
466 : ELSE
467 2466 : CALL fft_3d(fft_scratch%fft_plan(4), norm, zin, zout, stat)
468 : END IF
469 : END IF
470 :
471 678145 : CALL release_fft_scratch(fft_scratch)
472 :
473 678145 : IF (PRESENT(zout)) THEN
474 7702 : lo(1) = SIZE(zout, 1)
475 7702 : lo(2) = SIZE(zout, 2)
476 7702 : lo(3) = SIZE(zout, 3)
477 7702 : IF (n(1) /= lo(1) .OR. n(2) /= lo(2) .OR. n(3) /= lo(3)) THEN
478 0 : CPABORT("Size and dimension (zout) have to be the same.")
479 : END IF
480 : END IF
481 :
482 678145 : IF (PRESENT(status)) THEN
483 19283 : status = stat
484 : END IF
485 :
486 678145 : IF (test .AND. output_unit > 0) THEN
487 0 : IF (PRESENT(zout)) THEN
488 0 : out_sum = SUM(ABS(zout))
489 0 : WRITE (output_unit, '(A)') " Out of place 3D FFT (local) : fft3d_s"
490 0 : WRITE (output_unit, '(A,T60,3I7)') " Transform lengths ", n
491 0 : WRITE (output_unit, '(A,T60,3I7)') " Input array dimensions ", ld
492 0 : WRITE (output_unit, '(A,T60,3I7)') " Output array dimensions ", lo
493 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of input data ", in_sum
494 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of output data ", out_sum
495 : ELSE
496 0 : out_sum = SUM(ABS(zin))
497 0 : WRITE (output_unit, '(A)') " In place 3D FFT (local) : fft3d_s"
498 0 : WRITE (output_unit, '(A,T60,3I7)') " Transform lengths ", n
499 0 : WRITE (output_unit, '(A,T60,3I7)') " Input/output array dimensions ", ld
500 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of input data ", in_sum
501 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of output data ", out_sum
502 : END IF
503 : END IF
504 :
505 678145 : CALL timestop(handle)
506 :
507 678145 : END SUBROUTINE fft3d_s
508 :
509 : ! **************************************************************************************************
510 : !> \brief ...
511 : !> \param fsign ...
512 : !> \param n ...
513 : !> \param cin ...
514 : !> \param gin ...
515 : !> \param rs_group ...
516 : !> \param yzp ...
517 : !> \param nyzray ...
518 : !> \param bo ...
519 : !> \param status ...
520 : !> \param debug ...
521 : ! **************************************************************************************************
522 3529842 : SUBROUTINE fft3d_ps(fsign, n, cin, gin, rs_group, yzp, nyzray, &
523 3529842 : bo, status, debug)
524 :
525 : INTEGER, INTENT(IN) :: fsign
526 : INTEGER, DIMENSION(:), INTENT(IN) :: n
527 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
528 : INTENT(INOUT) :: cin
529 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
530 : INTENT(INOUT) :: gin
531 : TYPE(mp_cart_type), INTENT(IN) :: rs_group
532 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), &
533 : INTENT(IN) :: yzp
534 : INTEGER, CONTIGUOUS, DIMENSION(0:), INTENT(IN) :: nyzray
535 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:, :), &
536 : INTENT(IN) :: bo
537 : INTEGER, INTENT(OUT), OPTIONAL :: status
538 : LOGICAL, INTENT(IN), OPTIONAL :: debug
539 :
540 : CHARACTER(len=*), PARAMETER :: routineN = 'fft3d_ps'
541 :
542 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
543 3529842 : POINTER :: pbuf, qbuf, rbuf, sbuf
544 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
545 3529842 : POINTER :: tbuf
546 : INTEGER :: g_pos, handle, lg, lmax, mcx2, mcz1, mcz2, mg, mmax, mx1, mx2, my1, mz2, n1, n2, &
547 : nmax, numtask, nx, ny, nz, output_unit, r_dim(2), r_pos(2), rp, sign, stat
548 3529842 : INTEGER, ALLOCATABLE, DIMENSION(:) :: p2p
549 : LOGICAL :: test
550 : REAL(KIND=dp) :: norm, sum_data
551 24708894 : TYPE(fft_scratch_sizes) :: fft_scratch_size
552 : TYPE(fft_scratch_type), POINTER :: fft_scratch
553 :
554 3529842 : CALL timeset(routineN, handle)
555 3529842 : output_unit = cp_logger_get_default_io_unit()
556 :
557 3529842 : IF (PRESENT(debug)) THEN
558 3529842 : test = debug
559 : ELSE
560 : test = .FALSE.
561 : END IF
562 :
563 3529842 : g_pos = rs_group%mepos
564 3529842 : numtask = rs_group%num_pe
565 10589526 : r_dim = rs_group%num_pe_cart
566 10589526 : r_pos = rs_group%mepos_cart
567 :
568 3529842 : IF (fsign == FWFFT) THEN
569 7086120 : norm = 1.0_dp/REAL(PRODUCT(n), KIND=dp)
570 1758312 : ELSE IF (fsign == BWFFT) THEN
571 1758312 : norm = 1.0_dp
572 : ELSE
573 0 : CPABORT("Unknown FFT direction!")
574 : END IF
575 :
576 3529842 : sign = fsign
577 :
578 3529842 : lg = SIZE(gin, 1)
579 3529842 : mg = SIZE(gin, 2)
580 :
581 3529842 : nx = SIZE(cin, 1)
582 3529842 : ny = SIZE(cin, 2)
583 3529842 : nz = SIZE(cin, 3)
584 :
585 3529842 : IF (mg == 0) THEN
586 : mmax = 1
587 : ELSE
588 3529842 : mmax = mg
589 : END IF
590 3529842 : lmax = MAX(lg, (nx*ny*nz)/mmax + 1)
591 :
592 10589526 : ALLOCATE (p2p(0:numtask - 1))
593 :
594 3529842 : CALL rs_group%rank_compare(rs_group, p2p)
595 :
596 3529842 : rp = p2p(g_pos)
597 3529842 : mx1 = bo(2, 1, rp, 1) - bo(1, 1, rp, 1) + 1
598 3529842 : my1 = bo(2, 2, rp, 1) - bo(1, 2, rp, 1) + 1
599 3529842 : mx2 = bo(2, 1, rp, 2) - bo(1, 1, rp, 2) + 1
600 3529842 : mz2 = bo(2, 3, rp, 2) - bo(1, 3, rp, 2) + 1
601 :
602 10589526 : n1 = MAXVAL(bo(2, 1, :, 1) - bo(1, 1, :, 1) + 1)
603 10589526 : n2 = MAXVAL(bo(2, 2, :, 1) - bo(1, 2, :, 1) + 1)
604 3529842 : nmax = MAX((2*n2)/numtask, 2)*mx2*mz2
605 10589526 : nmax = MAX(nmax, n1*MAXVAL(nyzray))
606 10589526 : n1 = MAXVAL(bo(2, 1, :, 2))
607 10589526 : n2 = MAXVAL(bo(2, 3, :, 2))
608 :
609 3529842 : fft_scratch_size%nx = nx
610 3529842 : fft_scratch_size%ny = ny
611 3529842 : fft_scratch_size%nz = nz
612 3529842 : fft_scratch_size%lmax = lmax
613 3529842 : fft_scratch_size%mmax = mmax
614 3529842 : fft_scratch_size%mx1 = mx1
615 3529842 : fft_scratch_size%mx2 = mx2
616 3529842 : fft_scratch_size%my1 = my1
617 3529842 : fft_scratch_size%mz2 = mz2
618 3529842 : fft_scratch_size%lg = lg
619 3529842 : fft_scratch_size%mg = mg
620 3529842 : fft_scratch_size%nbx = n1
621 3529842 : fft_scratch_size%nbz = n2
622 10589526 : mcz1 = MAXVAL(bo(2, 3, :, 1) - bo(1, 3, :, 1) + 1)
623 10589526 : mcx2 = MAXVAL(bo(2, 1, :, 2) - bo(1, 1, :, 2) + 1)
624 10589526 : mcz2 = MAXVAL(bo(2, 3, :, 2) - bo(1, 3, :, 2) + 1)
625 3529842 : fft_scratch_size%mcz1 = mcz1
626 3529842 : fft_scratch_size%mcx2 = mcx2
627 3529842 : fft_scratch_size%mcz2 = mcz2
628 3529842 : fft_scratch_size%nmax = nmax
629 10589526 : fft_scratch_size%nmray = MAXVAL(nyzray)
630 3529842 : fft_scratch_size%nyzray = nyzray(g_pos)
631 3529842 : fft_scratch_size%rs_group = rs_group
632 10589526 : fft_scratch_size%g_pos = g_pos
633 10589526 : fft_scratch_size%r_pos = r_pos
634 10589526 : fft_scratch_size%r_dim = r_dim
635 3529842 : fft_scratch_size%numtask = numtask
636 :
637 3529842 : IF (test) THEN
638 8 : IF (g_pos == 0 .AND. output_unit > 0) THEN
639 4 : WRITE (output_unit, '(A)') " Parallel 3D FFT : fft3d_ps"
640 4 : WRITE (output_unit, '(A,T60,3I7)') " Transform lengths ", n
641 4 : WRITE (output_unit, '(A,T67,2I7)') " Array dimensions (gin) ", lg, mg
642 4 : WRITE (output_unit, '(A,T60,3I7)') " Array dimensions (cin) ", nx, ny, nz
643 : END IF
644 : END IF
645 :
646 3529842 : IF (r_dim(2) > 1) THEN
647 :
648 : !
649 : ! real space is distributed over x and y coordinate
650 : ! we have two stages of communication
651 : !
652 :
653 0 : IF (r_dim(1) == 1) THEN
654 0 : CPABORT("This processor distribution is not supported.")
655 : END IF
656 0 : CALL get_fft_scratch(fft_scratch, tf_type=300, n=n, fft_sizes=fft_scratch_size)
657 :
658 0 : IF (sign == FWFFT) THEN
659 : ! cin -> gin
660 :
661 0 : IF (test) THEN
662 0 : sum_data = ABS(SUM(cin))
663 0 : CALL rs_group%sum(sum_data)
664 0 : IF (g_pos == 0 .AND. output_unit > 0) THEN
665 0 : WRITE (output_unit, '(A)') " Two step communication algorithm "
666 0 : WRITE (output_unit, '(A,T60,3I7)') " Transform Z ", n(3), mx1*my1
667 0 : WRITE (output_unit, '(A,T60,3I7)') " Transform Y ", n(2), mx2*mz2
668 0 : WRITE (output_unit, '(A,T67,2I7)') " Transform X ", n(1), nyzray(g_pos)
669 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(1) ", sum_data
670 : END IF
671 : END IF
672 :
673 0 : pbuf => fft_scratch%p1buf
674 0 : qbuf => fft_scratch%p2buf
675 :
676 : ! FFT along z
677 0 : CALL fft_1d(fft_scratch%fft_plan(1), cin, qbuf, norm, stat)
678 :
679 0 : rbuf => fft_scratch%p3buf
680 :
681 0 : IF (test) THEN
682 0 : sum_data = ABS(SUM(qbuf))
683 0 : CALL rs_group%sum(sum_data)
684 0 : IF (g_pos == 0 .AND. output_unit > 0) THEN
685 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(2) T", sum_data
686 : END IF
687 : END IF
688 :
689 : ! Exchange data ( transpose of matrix )
690 0 : CALL cube_transpose_2(qbuf, bo(:, :, :, 1), bo(:, :, :, 2), rbuf, fft_scratch)
691 :
692 0 : IF (test) THEN
693 0 : sum_data = ABS(SUM(rbuf))
694 0 : CALL rs_group%sum(sum_data)
695 0 : IF (g_pos == 0 .AND. output_unit > 0) THEN
696 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(3) T", sum_data
697 : END IF
698 : END IF
699 :
700 0 : pbuf => fft_scratch%p4buf
701 :
702 : ! FFT along y
703 0 : CALL fft_1d(fft_scratch%fft_plan(2), rbuf, pbuf, 1.0_dp, stat)
704 :
705 0 : qbuf => fft_scratch%p5buf
706 :
707 0 : IF (test) THEN
708 0 : sum_data = ABS(SUM(pbuf))
709 0 : CALL rs_group%sum(sum_data)
710 0 : IF (g_pos == 0 .AND. output_unit > 0) THEN
711 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(4) TS", sum_data
712 : END IF
713 : END IF
714 :
715 : ! Exchange data ( transpose of matrix ) and sort
716 : CALL xz_to_yz(pbuf, rs_group, r_dim, g_pos, p2p, yzp, nyzray, &
717 0 : bo(:, :, :, 2), qbuf, fft_scratch)
718 :
719 0 : IF (test) THEN
720 0 : sum_data = ABS(SUM(qbuf))
721 0 : CALL rs_group%sum(sum_data)
722 0 : IF (g_pos == 0 .AND. output_unit > 0) THEN
723 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(5) TS", sum_data
724 : END IF
725 : END IF
726 :
727 : ! FFT along x
728 0 : CALL fft_1d(fft_scratch%fft_plan(3), qbuf, gin, 1.0_dp, stat)
729 :
730 0 : IF (test) THEN
731 0 : sum_data = ABS(SUM(gin))
732 0 : CALL rs_group%sum(sum_data)
733 0 : IF (g_pos == 0 .AND. output_unit > 0) THEN
734 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(6) ", sum_data
735 : END IF
736 : END IF
737 :
738 0 : ELSE IF (sign == BWFFT) THEN
739 : ! gin -> cin
740 :
741 0 : IF (test) THEN
742 0 : sum_data = ABS(SUM(gin))
743 0 : CALL rs_group%sum(sum_data)
744 0 : IF (g_pos == 0 .AND. output_unit > 0) THEN
745 0 : WRITE (output_unit, '(A)') " Two step communication algorithm "
746 0 : WRITE (output_unit, '(A,T67,2I7)') " Transform X ", n(1), nyzray(g_pos)
747 0 : WRITE (output_unit, '(A,T60,3I7)') " Transform Y ", n(2), mx2*mz2
748 0 : WRITE (output_unit, '(A,T60,3I7)') " Transform Z ", n(3), mx1*my1
749 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(1) ", sum_data
750 : END IF
751 : END IF
752 :
753 0 : pbuf => fft_scratch%p7buf
754 :
755 : ! FFT along x
756 0 : CALL fft_1d(fft_scratch%fft_plan(4), gin, pbuf, norm, stat)
757 :
758 0 : qbuf => fft_scratch%p4buf
759 :
760 0 : IF (test) THEN
761 0 : sum_data = ABS(SUM(pbuf))
762 0 : CALL rs_group%sum(sum_data)
763 0 : IF (g_pos == 0 .AND. output_unit > 0) THEN
764 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(2) TS", sum_data
765 : END IF
766 : END IF
767 :
768 : ! Exchange data ( transpose of matrix ) and sort
769 : CALL yz_to_xz(pbuf, rs_group, r_dim, g_pos, p2p, yzp, nyzray, &
770 0 : bo(:, :, :, 2), qbuf, fft_scratch)
771 :
772 0 : IF (test) THEN
773 0 : sum_data = ABS(SUM(qbuf))
774 0 : CALL rs_group%sum(sum_data)
775 0 : IF (g_pos == 0 .AND. output_unit > 0) THEN
776 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(3) TS", sum_data
777 : END IF
778 : END IF
779 :
780 0 : rbuf => fft_scratch%p3buf
781 :
782 : ! FFT along y
783 0 : CALL fft_1d(fft_scratch%fft_plan(5), qbuf, rbuf, 1.0_dp, stat)
784 :
785 0 : pbuf => fft_scratch%p2buf
786 :
787 0 : IF (test) THEN
788 0 : sum_data = ABS(SUM(rbuf))
789 0 : CALL rs_group%sum(sum_data)
790 0 : IF (g_pos == 0 .AND. output_unit > 0) THEN
791 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(4) T", sum_data
792 : END IF
793 : END IF
794 :
795 : ! Exchange data ( transpose of matrix )
796 0 : CALL cube_transpose_1(rbuf, bo(:, :, :, 2), bo(:, :, :, 1), pbuf, fft_scratch)
797 :
798 0 : IF (test) THEN
799 0 : sum_data = ABS(SUM(pbuf))
800 0 : CALL rs_group%sum(sum_data)
801 0 : IF (g_pos == 0 .AND. output_unit > 0) THEN
802 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(5) T", sum_data
803 : END IF
804 : END IF
805 :
806 0 : qbuf => fft_scratch%p1buf
807 :
808 : ! FFT along z
809 0 : CALL fft_1d(fft_scratch%fft_plan(6), pbuf, cin, 1.0_dp, stat)
810 :
811 0 : IF (test) THEN
812 0 : sum_data = ABS(SUM(cin))
813 0 : CALL rs_group%sum(sum_data)
814 0 : IF (g_pos == 0 .AND. output_unit > 0) THEN
815 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(6) ", sum_data
816 : END IF
817 : END IF
818 :
819 : ELSE
820 :
821 0 : CPABORT("Illegal fsign parameter.")
822 :
823 : END IF
824 :
825 0 : CALL release_fft_scratch(fft_scratch)
826 :
827 : ELSE
828 :
829 : !
830 : ! real space is only distributed over x coordinate
831 : ! we have one stage of communication, after the transform of
832 : ! direction x
833 : !
834 :
835 3529842 : CALL get_fft_scratch(fft_scratch, tf_type=200, n=n, fft_sizes=fft_scratch_size)
836 :
837 3529842 : sbuf => fft_scratch%r1buf
838 3529842 : tbuf => fft_scratch%tbuf
839 :
840 3529842 : sbuf = z_zero
841 3529842 : tbuf = z_zero
842 :
843 3529842 : IF (sign == FWFFT) THEN
844 : ! cin -> gin
845 :
846 1771530 : IF (test) THEN
847 9284 : sum_data = ABS(SUM(cin))
848 4 : CALL rs_group%sum(sum_data)
849 4 : IF (g_pos == 0 .AND. output_unit > 0) THEN
850 2 : WRITE (output_unit, '(A)') " One step communication algorithm "
851 2 : WRITE (output_unit, '(A,T60,3I7)') " Transform YZ ", n(2), n(3), nx
852 2 : WRITE (output_unit, '(A,T67,2I7)') " Transform X ", n(1), nyzray(g_pos)
853 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(1) ", sum_data
854 : END IF
855 : END IF
856 :
857 : ! FFT along y and z
858 1771530 : CALL fft_1d(fft_scratch%fft_plan(1), cin, sbuf, 1._dp, stat)
859 1771530 : CALL fft_1d(fft_scratch%fft_plan(2), sbuf, tbuf, 1._dp, stat)
860 :
861 1771530 : IF (test) THEN
862 8740 : sum_data = ABS(SUM(tbuf))
863 4 : CALL rs_group%sum(sum_data)
864 4 : IF (g_pos == 0 .AND. output_unit > 0) THEN
865 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(2) TS", sum_data
866 : END IF
867 : END IF
868 :
869 : ! Exchange data ( transpose of matrix ) and sort
870 : CALL yz_to_x(tbuf, rs_group, g_pos, p2p, yzp, nyzray, &
871 1771530 : bo(:, :, :, 2), sbuf, fft_scratch)
872 :
873 1771530 : IF (test) THEN
874 8776 : sum_data = ABS(SUM(sbuf))
875 4 : CALL rs_group%sum(sum_data)
876 4 : IF (g_pos == 0 .AND. output_unit > 0) THEN
877 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(3) TS", sum_data
878 : END IF
879 : END IF
880 : ! FFT along x
881 1771530 : CALL fft_1d(fft_scratch%fft_plan(3), sbuf, gin, norm, stat)
882 :
883 1771530 : IF (test) THEN
884 8708 : sum_data = ABS(SUM(gin))
885 4 : CALL rs_group%sum(sum_data)
886 4 : IF (g_pos == 0 .AND. output_unit > 0) THEN
887 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(4) ", sum_data
888 : END IF
889 : END IF
890 :
891 1758312 : ELSE IF (sign == BWFFT) THEN
892 : ! gin -> cin
893 :
894 1758312 : IF (test) THEN
895 8708 : sum_data = ABS(SUM(gin))
896 4 : CALL rs_group%sum(sum_data)
897 4 : IF (g_pos == 0 .AND. output_unit > 0) THEN
898 2 : WRITE (output_unit, '(A)') " One step communication algorithm "
899 2 : WRITE (output_unit, '(A,T67,2I7)') " Transform X ", n(1), nyzray(g_pos)
900 2 : WRITE (output_unit, '(A,T60,3I7)') " Transform YZ ", n(2), n(3), nx
901 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(1) ", sum_data
902 : END IF
903 : END IF
904 :
905 : ! FFT along x
906 1758312 : CALL fft_1d(fft_scratch%fft_plan(4), gin, sbuf, norm, stat)
907 :
908 1758312 : IF (test) THEN
909 8776 : sum_data = ABS(SUM(sbuf))
910 4 : CALL rs_group%sum(sum_data)
911 4 : IF (g_pos == 0 .AND. output_unit > 0) THEN
912 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(2) TS", sum_data
913 : END IF
914 : END IF
915 :
916 : ! Exchange data ( transpose of matrix ) and sort
917 : CALL x_to_yz(sbuf, rs_group, g_pos, p2p, yzp, nyzray, &
918 1758312 : bo(:, :, :, 2), tbuf, fft_scratch)
919 :
920 1758312 : IF (test) THEN
921 8740 : sum_data = ABS(SUM(tbuf))
922 4 : CALL rs_group%sum(sum_data)
923 4 : IF (g_pos == 0 .AND. output_unit > 0) THEN
924 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(3) TS", sum_data
925 : END IF
926 : END IF
927 :
928 : ! FFT along y and z
929 1758312 : CALL fft_1d(fft_scratch%fft_plan(5), tbuf, sbuf, 1._dp, stat)
930 1758312 : CALL fft_1d(fft_scratch%fft_plan(6), sbuf, cin, 1._dp, stat)
931 :
932 1758312 : IF (test) THEN
933 9284 : sum_data = ABS(SUM(cin))
934 4 : CALL rs_group%sum(sum_data)
935 4 : IF (g_pos == 0 .AND. output_unit > 0) THEN
936 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(4) ", sum_data
937 : END IF
938 : END IF
939 : ELSE
940 0 : CPABORT("Illegal fsign parameter.")
941 : END IF
942 :
943 3529842 : CALL release_fft_scratch(fft_scratch)
944 :
945 : END IF
946 :
947 3529842 : DEALLOCATE (p2p)
948 :
949 3529842 : IF (PRESENT(status)) THEN
950 0 : status = stat
951 : END IF
952 3529842 : CALL timestop(handle)
953 :
954 7059684 : END SUBROUTINE fft3d_ps
955 :
956 : ! **************************************************************************************************
957 : !> \brief ...
958 : !> \param fsign ...
959 : !> \param n ...
960 : !> \param zin ...
961 : !> \param gin ...
962 : !> \param group ...
963 : !> \param bo ...
964 : !> \param status ...
965 : !> \param debug ...
966 : ! **************************************************************************************************
967 208 : SUBROUTINE fft3d_pb(fsign, n, zin, gin, group, bo, status, debug)
968 :
969 : INTEGER, INTENT(IN) :: fsign
970 : INTEGER, DIMENSION(3), INTENT(IN) :: n
971 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
972 : INTENT(INOUT) :: zin
973 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
974 : INTENT(INOUT) :: gin
975 : TYPE(mp_cart_type), INTENT(IN) :: group
976 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:, :), &
977 : INTENT(IN) :: bo
978 : INTEGER, INTENT(OUT), OPTIONAL :: status
979 : LOGICAL, INTENT(IN), OPTIONAL :: debug
980 :
981 : CHARACTER(len=*), PARAMETER :: routineN = 'fft3d_pb'
982 :
983 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
984 208 : POINTER :: abuf, bbuf
985 : INTEGER :: handle, lg(2), lz(3), mcx2, mcy3, mcz1, &
986 : mcz2, mx1, mx2, mx3, my1, my2, my3, &
987 : my_pos, mz1, mz2, mz3, output_unit, &
988 : sign, stat
989 : INTEGER, DIMENSION(2) :: dim
990 : LOGICAL :: test
991 : REAL(KIND=dp) :: norm, sum_data
992 1456 : TYPE(fft_scratch_sizes) :: fft_scratch_size
993 : TYPE(fft_scratch_type), POINTER :: fft_scratch
994 :
995 : !------------------------------------------------------------------------------
996 : ! "Real Space" 1) xyZ or 1) xYZ
997 : ! 2) xYz or not used
998 : ! "G Space" 3) Xyz or 3) XYz
999 : !
1000 : ! There is one communicator (2-dimensional) for all distributions
1001 : ! np = n1 * n2, where np is the total number of processors
1002 : ! If n2 = 1, we have the second case and only one transpose step is needed
1003 : !
1004 : ! Assignment of dimensions to axis for different steps
1005 : ! First case: 1) n1=x; n2=y
1006 : ! 2) n1=x; n2=z
1007 : ! 3) n1=y; n2=z
1008 : ! Second case 1) n1=x
1009 : ! 3) n1=z
1010 : !
1011 : ! The more general case with two communicators for the initial and final
1012 : ! distribution is not covered.
1013 : !------------------------------------------------------------------------------
1014 :
1015 208 : CALL timeset(routineN, handle)
1016 208 : output_unit = cp_logger_get_default_io_unit()
1017 :
1018 624 : dim = group%num_pe_cart
1019 208 : my_pos = group%mepos
1020 :
1021 208 : IF (PRESENT(debug)) THEN
1022 208 : test = debug
1023 : ELSE
1024 : test = .FALSE.
1025 : END IF
1026 :
1027 208 : IF (fsign == FWFFT) THEN
1028 416 : norm = 1.0_dp/REAL(PRODUCT(n), KIND=dp)
1029 104 : ELSE IF (fsign == BWFFT) THEN
1030 104 : norm = 1.0_dp
1031 : ELSE
1032 0 : CPABORT("Unknown FFT direction!")
1033 : END IF
1034 :
1035 208 : sign = fsign
1036 :
1037 208 : IF (test) THEN
1038 8 : lg(1) = SIZE(gin, 1)
1039 8 : lg(2) = SIZE(gin, 2)
1040 8 : lz(1) = SIZE(zin, 1)
1041 8 : lz(2) = SIZE(zin, 2)
1042 8 : lz(3) = SIZE(zin, 3)
1043 8 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1044 4 : WRITE (output_unit, '(A)') " Parallel 3D FFT : fft3d_pb"
1045 4 : WRITE (output_unit, '(A,T60,3I7)') " Transform lengths ", n
1046 4 : WRITE (output_unit, '(A,T67,2I7)') " Array dimensions (gin) ", lg
1047 4 : WRITE (output_unit, '(A,T60,3I7)') " Array dimensions (cin) ", lz
1048 : END IF
1049 : END IF
1050 :
1051 208 : mx1 = bo(2, 1, my_pos, 1) - bo(1, 1, my_pos, 1) + 1
1052 208 : my1 = bo(2, 2, my_pos, 1) - bo(1, 2, my_pos, 1) + 1
1053 208 : mz1 = bo(2, 3, my_pos, 1) - bo(1, 3, my_pos, 1) + 1
1054 208 : mx2 = bo(2, 1, my_pos, 2) - bo(1, 1, my_pos, 2) + 1
1055 208 : my2 = bo(2, 2, my_pos, 2) - bo(1, 2, my_pos, 2) + 1
1056 208 : mz2 = bo(2, 3, my_pos, 2) - bo(1, 3, my_pos, 2) + 1
1057 208 : mx3 = bo(2, 1, my_pos, 3) - bo(1, 1, my_pos, 3) + 1
1058 208 : my3 = bo(2, 2, my_pos, 3) - bo(1, 2, my_pos, 3) + 1
1059 208 : mz3 = bo(2, 3, my_pos, 3) - bo(1, 3, my_pos, 3) + 1
1060 208 : fft_scratch_size%mx1 = mx1
1061 208 : fft_scratch_size%mx2 = mx2
1062 208 : fft_scratch_size%mx3 = mx3
1063 208 : fft_scratch_size%my1 = my1
1064 208 : fft_scratch_size%my2 = my2
1065 208 : fft_scratch_size%my3 = my3
1066 208 : fft_scratch_size%mz1 = mz1
1067 208 : fft_scratch_size%mz2 = mz2
1068 208 : fft_scratch_size%mz3 = mz3
1069 624 : mcz1 = MAXVAL(bo(2, 3, :, 1) - bo(1, 3, :, 1) + 1)
1070 624 : mcx2 = MAXVAL(bo(2, 1, :, 2) - bo(1, 1, :, 2) + 1)
1071 624 : mcz2 = MAXVAL(bo(2, 3, :, 2) - bo(1, 3, :, 2) + 1)
1072 624 : mcy3 = MAXVAL(bo(2, 2, :, 3) - bo(1, 2, :, 3) + 1)
1073 208 : fft_scratch_size%mcz1 = mcz1
1074 208 : fft_scratch_size%mcx2 = mcx2
1075 208 : fft_scratch_size%mcz2 = mcz2
1076 208 : fft_scratch_size%mcy3 = mcy3
1077 208 : fft_scratch_size%rs_group = group
1078 624 : fft_scratch_size%g_pos = my_pos
1079 208 : fft_scratch_size%numtask = DIM(1)*DIM(2)
1080 :
1081 208 : IF (DIM(1) > 1 .AND. DIM(2) > 1) THEN
1082 :
1083 : !
1084 : ! First case; two stages of communication
1085 : !
1086 :
1087 0 : CALL get_fft_scratch(fft_scratch, tf_type=100, n=n, fft_sizes=fft_scratch_size)
1088 :
1089 0 : IF (sign == FWFFT) THEN
1090 : ! Stage 1 -> 3
1091 :
1092 0 : bbuf => fft_scratch%a2buf
1093 :
1094 0 : IF (test) THEN
1095 0 : sum_data = ABS(SUM(zin))
1096 0 : CALL group%sum(sum_data)
1097 0 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1098 0 : WRITE (output_unit, '(A)') " Two step communication algorithm "
1099 0 : WRITE (output_unit, '(A,T67,2I7)') " Transform Z ", n(3), mx1*my1
1100 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(1) ", sum_data
1101 : END IF
1102 : END IF
1103 :
1104 : ! FFT along z
1105 0 : CALL fft_1d(fft_scratch%fft_plan(1), zin, bbuf, norm, stat)
1106 :
1107 0 : abuf => fft_scratch%a3buf
1108 :
1109 0 : IF (test) THEN
1110 0 : sum_data = ABS(SUM(bbuf))
1111 0 : CALL group%sum(sum_data)
1112 0 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1113 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(2) T", sum_data
1114 : END IF
1115 : END IF
1116 :
1117 0 : CALL cube_transpose_2(bbuf, bo(:, :, :, 1), bo(:, :, :, 2), abuf, fft_scratch)
1118 :
1119 0 : bbuf => fft_scratch%a4buf
1120 :
1121 0 : IF (test) THEN
1122 0 : sum_data = ABS(SUM(abuf))
1123 0 : CALL group%sum(sum_data)
1124 0 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1125 0 : WRITE (output_unit, '(A,T67,2I7)') " Transform Y ", n(2), mx2*mz2
1126 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(3) ", sum_data
1127 : END IF
1128 : END IF
1129 :
1130 : ! FFT along y
1131 0 : CALL fft_1d(fft_scratch%fft_plan(2), abuf, bbuf, 1.0_dp, stat)
1132 :
1133 0 : abuf => fft_scratch%a5buf
1134 :
1135 0 : IF (test) THEN
1136 0 : sum_data = ABS(SUM(bbuf))
1137 0 : CALL group%sum(sum_data)
1138 0 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1139 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(4) T", sum_data
1140 : END IF
1141 : END IF
1142 :
1143 0 : CALL cube_transpose_4(bbuf, bo(:, :, :, 2), bo(:, :, :, 3), abuf, fft_scratch)
1144 :
1145 0 : IF (test) THEN
1146 0 : sum_data = ABS(SUM(abuf))
1147 0 : CALL group%sum(sum_data)
1148 0 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1149 0 : WRITE (output_unit, '(A,T67,2I7)') " Transform X ", n(1), my3*mz3
1150 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(5) ", sum_data
1151 : END IF
1152 : END IF
1153 :
1154 : ! FFT along x
1155 0 : CALL fft_1d(fft_scratch%fft_plan(3), abuf, gin, 1.0_dp, stat)
1156 :
1157 0 : IF (test) THEN
1158 0 : sum_data = ABS(SUM(gin))
1159 0 : CALL group%sum(sum_data)
1160 0 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1161 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(6) ", sum_data
1162 : END IF
1163 : END IF
1164 :
1165 0 : ELSE IF (sign == BWFFT) THEN
1166 : ! Stage 3 -> 1
1167 :
1168 0 : bbuf => fft_scratch%a5buf
1169 :
1170 0 : IF (test) THEN
1171 0 : sum_data = ABS(SUM(gin))
1172 0 : CALL group%sum(sum_data)
1173 0 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1174 0 : WRITE (output_unit, '(A)') " Two step communication algorithm "
1175 0 : WRITE (output_unit, '(A,T67,2I7)') " Transform X ", n(1), my3*mz3
1176 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(1) ", sum_data
1177 : END IF
1178 : END IF
1179 :
1180 : ! FFT along x
1181 0 : CALL fft_1d(fft_scratch%fft_plan(4), gin, bbuf, 1.0_dp, stat)
1182 :
1183 0 : abuf => fft_scratch%a4buf
1184 :
1185 0 : IF (test) THEN
1186 0 : sum_data = ABS(SUM(bbuf))
1187 0 : CALL group%sum(sum_data)
1188 0 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1189 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(2) T", sum_data
1190 : END IF
1191 : END IF
1192 :
1193 0 : CALL cube_transpose_3(bbuf, bo(:, :, :, 3), bo(:, :, :, 2), abuf, fft_scratch)
1194 :
1195 0 : bbuf => fft_scratch%a3buf
1196 :
1197 0 : IF (test) THEN
1198 0 : sum_data = ABS(SUM(abuf))
1199 0 : CALL group%sum(sum_data)
1200 0 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1201 0 : WRITE (output_unit, '(A,T67,2I7)') " Transform Y ", n(2), mx2*mz2
1202 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(3) ", sum_data
1203 : END IF
1204 : END IF
1205 :
1206 : ! FFT along y
1207 0 : CALL fft_1d(fft_scratch%fft_plan(5), abuf, bbuf, 1.0_dp, stat)
1208 :
1209 0 : abuf => fft_scratch%a2buf
1210 :
1211 0 : IF (test) THEN
1212 0 : sum_data = ABS(SUM(bbuf))
1213 0 : CALL group%sum(sum_data)
1214 0 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1215 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(4) T", sum_data
1216 : END IF
1217 : END IF
1218 :
1219 0 : CALL cube_transpose_1(bbuf, bo(:, :, :, 2), bo(:, :, :, 1), abuf, fft_scratch)
1220 :
1221 0 : IF (test) THEN
1222 0 : sum_data = ABS(SUM(abuf))
1223 0 : CALL group%sum(sum_data)
1224 0 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1225 0 : WRITE (output_unit, '(A,T67,2I7)') " Transform Z ", n(3), mx1*my1
1226 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(5) ", sum_data
1227 : END IF
1228 : END IF
1229 :
1230 : ! FFT along z
1231 0 : CALL fft_1d(fft_scratch%fft_plan(6), abuf, zin, norm, stat)
1232 :
1233 0 : IF (test) THEN
1234 0 : sum_data = ABS(SUM(zin))
1235 0 : CALL group%sum(sum_data)
1236 0 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1237 0 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(6) ", sum_data
1238 : END IF
1239 : END IF
1240 :
1241 : ELSE
1242 0 : CPABORT("Illegal fsign parameter.")
1243 : END IF
1244 :
1245 0 : CALL release_fft_scratch(fft_scratch)
1246 :
1247 208 : ELSE IF (DIM(2) == 1) THEN
1248 :
1249 : !
1250 : ! Second case; one stage of communication
1251 : !
1252 :
1253 208 : CALL get_fft_scratch(fft_scratch, tf_type=101, n=n, fft_sizes=fft_scratch_size)
1254 :
1255 208 : IF (sign == FWFFT) THEN
1256 : ! Stage 1 -> 3
1257 :
1258 104 : IF (test) THEN
1259 9284 : sum_data = ABS(SUM(zin))
1260 4 : CALL group%sum(sum_data)
1261 4 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1262 2 : WRITE (output_unit, '(A)') " one step communication algorithm "
1263 2 : WRITE (output_unit, '(A,T67,2I7)') " Transform Z ", n(3), mx1*my1
1264 2 : WRITE (output_unit, '(A,T67,2I7)') " Transform Y ", n(2), mx1*mz1
1265 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(1) ", sum_data
1266 : END IF
1267 : END IF
1268 :
1269 104 : abuf => fft_scratch%a3buf
1270 104 : bbuf => fft_scratch%a4buf
1271 : ! FFT along z and y
1272 104 : CALL fft_1d(fft_scratch%fft_plan(1), zin, abuf, norm, stat)
1273 104 : CALL fft_1d(fft_scratch%fft_plan(2), abuf, bbuf, 1.0_dp, stat)
1274 :
1275 104 : abuf => fft_scratch%a5buf
1276 :
1277 104 : IF (test) THEN
1278 8708 : sum_data = ABS(SUM(bbuf))
1279 4 : CALL group%sum(sum_data)
1280 4 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1281 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(2) T", sum_data
1282 : END IF
1283 : END IF
1284 :
1285 104 : CALL cube_transpose_6(bbuf, group, bo(:, :, :, 1), bo(:, :, :, 3), abuf, fft_scratch)
1286 :
1287 104 : IF (test) THEN
1288 8260 : sum_data = ABS(SUM(abuf))
1289 4 : CALL group%sum(sum_data)
1290 4 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1291 2 : WRITE (output_unit, '(A,T67,2I7)') " Transform X ", n(1), my3*mz3
1292 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(3) ", sum_data
1293 : END IF
1294 : END IF
1295 :
1296 : ! FFT along x
1297 104 : CALL fft_1d(fft_scratch%fft_plan(3), abuf, gin, 1.0_dp, stat)
1298 :
1299 104 : IF (test) THEN
1300 8708 : sum_data = ABS(SUM(gin))
1301 4 : CALL group%sum(sum_data)
1302 4 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1303 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(4) ", sum_data
1304 : END IF
1305 : END IF
1306 :
1307 104 : ELSE IF (sign == BWFFT) THEN
1308 : ! Stage 3 -> 1
1309 :
1310 104 : IF (test) THEN
1311 8708 : sum_data = ABS(SUM(gin))
1312 4 : CALL group%sum(sum_data)
1313 4 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1314 2 : WRITE (output_unit, '(A)') " one step communication algorithm "
1315 2 : WRITE (output_unit, '(A,T67,2I7)') " Transform X ", n(1), my3*mz3
1316 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(1) ", sum_data
1317 : END IF
1318 : END IF
1319 :
1320 104 : bbuf => fft_scratch%a5buf
1321 :
1322 : ! FFT along x
1323 104 : CALL fft_1d(fft_scratch%fft_plan(4), gin, bbuf, 1.0_dp, stat)
1324 :
1325 104 : abuf => fft_scratch%a4buf
1326 :
1327 104 : IF (test) THEN
1328 8260 : sum_data = ABS(SUM(bbuf))
1329 4 : CALL group%sum(sum_data)
1330 4 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1331 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(2) T", sum_data
1332 : END IF
1333 : END IF
1334 :
1335 104 : CALL cube_transpose_5(bbuf, group, bo(:, :, :, 3), bo(:, :, :, 1), abuf, fft_scratch)
1336 :
1337 104 : bbuf => fft_scratch%a3buf
1338 :
1339 104 : IF (test) THEN
1340 8708 : sum_data = ABS(SUM(abuf))
1341 4 : CALL group%sum(sum_data)
1342 4 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1343 2 : WRITE (output_unit, '(A,T67,2I7)') " Transform Y ", n(2), mx1*mz1
1344 2 : WRITE (output_unit, '(A,T67,2I7)') " Transform Z ", n(3), mx1*my1
1345 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(3) ", sum_data
1346 : END IF
1347 : END IF
1348 :
1349 : ! FFT along y
1350 104 : CALL fft_1d(fft_scratch%fft_plan(5), abuf, bbuf, 1.0_dp, stat)
1351 :
1352 : ! FFT along z
1353 104 : CALL fft_1d(fft_scratch%fft_plan(6), bbuf, zin, norm, stat)
1354 :
1355 104 : IF (test) THEN
1356 9284 : sum_data = ABS(SUM(zin))
1357 4 : CALL group%sum(sum_data)
1358 4 : IF (my_pos == 0 .AND. output_unit > 0) THEN
1359 2 : WRITE (output_unit, '(A,T61,E20.14)') " Sum of data(4) ", sum_data
1360 : END IF
1361 : END IF
1362 :
1363 : ELSE
1364 0 : CPABORT("Illegal fsign parameter.")
1365 : END IF
1366 :
1367 208 : CALL release_fft_scratch(fft_scratch)
1368 :
1369 : ELSE
1370 :
1371 0 : CPABORT("Partition not implemented.")
1372 :
1373 : END IF
1374 :
1375 208 : IF (PRESENT(status)) THEN
1376 0 : status = stat
1377 : END IF
1378 :
1379 208 : CALL timestop(handle)
1380 :
1381 416 : END SUBROUTINE fft3d_pb
1382 :
1383 : ! **************************************************************************************************
1384 : !> \brief ...
1385 : !> \param sb ...
1386 : !> \param group ...
1387 : !> \param my_pos ...
1388 : !> \param p2p ...
1389 : !> \param yzp ...
1390 : !> \param nray ...
1391 : !> \param bo ...
1392 : !> \param tb ...
1393 : !> \param fft_scratch ...
1394 : !> \par History
1395 : !> 15. Feb. 2006 : single precision all_to_all
1396 : !> \author JGH (14-Jan-2001)
1397 : ! **************************************************************************************************
1398 1758312 : SUBROUTINE x_to_yz(sb, group, my_pos, p2p, yzp, nray, bo, tb, fft_scratch)
1399 :
1400 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1401 : INTENT(IN) :: sb
1402 :
1403 : CLASS(mp_comm_type), INTENT(IN) :: group
1404 : INTEGER, INTENT(IN) :: my_pos
1405 : INTEGER, CONTIGUOUS, DIMENSION(0:), INTENT(IN) :: p2p
1406 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), &
1407 : INTENT(IN) :: yzp
1408 : INTEGER, CONTIGUOUS, DIMENSION(0:), INTENT(IN) :: nray
1409 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), &
1410 : INTENT(IN) :: bo
1411 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
1412 : INTENT(INOUT) :: tb
1413 : TYPE(fft_scratch_type), INTENT(IN) :: fft_scratch
1414 :
1415 : CHARACTER(len=*), PARAMETER :: routineN = 'x_to_yz'
1416 :
1417 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1418 1758312 : POINTER :: rr
1419 : COMPLEX(KIND=sp), CONTIGUOUS, DIMENSION(:, :), &
1420 1758312 : POINTER :: ss, tt
1421 : INTEGER :: handle, ip, ir, ix, ixx, iy, iz, mpr, &
1422 : nm, np, nr, nx
1423 1758312 : INTEGER, CONTIGUOUS, DIMENSION(:), POINTER :: rcount, rdispl, scount, sdispl
1424 :
1425 1758312 : CALL timeset(routineN, handle)
1426 :
1427 1758312 : np = SIZE(p2p)
1428 1758312 : scount => fft_scratch%scount
1429 1758312 : rcount => fft_scratch%rcount
1430 1758312 : sdispl => fft_scratch%sdispl
1431 1758312 : rdispl => fft_scratch%rdispl
1432 :
1433 1758312 : IF (alltoall_sgl) THEN
1434 230 : ss => fft_scratch%ss
1435 230 : tt => fft_scratch%tt
1436 3441099 : ss(:, :) = CMPLX(sb(:, :), KIND=sp)
1437 230 : tt(:, :) = 0._sp
1438 : ELSE
1439 1758082 : rr => fft_scratch%rr
1440 : END IF
1441 :
1442 1758312 : mpr = p2p(my_pos)
1443 5274936 : nm = MAXVAL(nray(0:np - 1))
1444 1758312 : nr = nray(my_pos)
1445 : !$OMP PARALLEL DO DEFAULT(NONE), &
1446 : !$OMP PRIVATE(ix,nx), &
1447 1758312 : !$OMP SHARED(np,p2p,bo,nr,scount,sdispl)
1448 : DO ip = 0, np - 1
1449 : ix = p2p(ip)
1450 : nx = bo(2, 1, ix) - bo(1, 1, ix) + 1
1451 : scount(ip) = nr*nx
1452 : sdispl(ip) = nr*(bo(1, 1, ix) - 1)
1453 : END DO
1454 : !$OMP END PARALLEL DO
1455 1758312 : nx = bo(2, 1, mpr) - bo(1, 1, mpr) + 1
1456 : !$OMP PARALLEL DO DEFAULT(NONE), &
1457 : !$OMP PRIVATE(nr), &
1458 1758312 : !$OMP SHARED(np,nray,nx,rcount,rdispl,nm)
1459 : DO ip = 0, np - 1
1460 : nr = nray(ip)
1461 : rcount(ip) = nr*nx
1462 : rdispl(ip) = nm*nx*ip
1463 : END DO
1464 : !$OMP END PARALLEL DO
1465 1758312 : IF (alltoall_sgl) THEN
1466 230 : CALL group%alltoall(ss, scount, sdispl, tt, rcount, rdispl)
1467 : ELSE
1468 1758082 : CALL group%alltoall(sb, scount, sdispl, rr, rcount, rdispl)
1469 : END IF
1470 :
1471 1758312 : nx = bo(2, 1, mpr) - bo(1, 1, mpr) + 1
1472 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(2) &
1473 : !$OMP PRIVATE(ixx,ir,iy,iz,ix) &
1474 1758312 : !$OMP SHARED(np,nray,nx,alltoall_sgl,yzp,tt,rr,tb)
1475 : DO ip = 0, np - 1
1476 : DO ix = 1, nx
1477 : ixx = nray(ip)*(ix - 1)
1478 : IF (alltoall_sgl) THEN
1479 : DO ir = 1, nray(ip)
1480 : iy = yzp(1, ir, ip)
1481 : iz = yzp(2, ir, ip)
1482 : tb(iy, iz, ix) = tt(ir + ixx, ip)
1483 : END DO
1484 : ELSE
1485 : DO ir = 1, nray(ip)
1486 : iy = yzp(1, ir, ip)
1487 : iz = yzp(2, ir, ip)
1488 : tb(iy, iz, ix) = rr(ir + ixx, ip)
1489 : END DO
1490 : END IF
1491 : END DO
1492 : END DO
1493 : !$OMP END PARALLEL DO
1494 :
1495 1758312 : CALL timestop(handle)
1496 :
1497 1758312 : END SUBROUTINE x_to_yz
1498 :
1499 : ! **************************************************************************************************
1500 : !> \brief ...
1501 : !> \param tb ...
1502 : !> \param group ...
1503 : !> \param my_pos ...
1504 : !> \param p2p ...
1505 : !> \param yzp ...
1506 : !> \param nray ...
1507 : !> \param bo ...
1508 : !> \param sb ...
1509 : !> \param fft_scratch ...
1510 : !> \par History
1511 : !> 15. Feb. 2006 : single precision all_to_all
1512 : !> \author JGH (14-Jan-2001)
1513 : ! **************************************************************************************************
1514 1771530 : SUBROUTINE yz_to_x(tb, group, my_pos, p2p, yzp, nray, bo, sb, fft_scratch)
1515 :
1516 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
1517 : INTENT(IN) :: tb
1518 :
1519 : CLASS(mp_comm_type), INTENT(IN) :: group
1520 : INTEGER, INTENT(IN) :: my_pos
1521 : INTEGER, CONTIGUOUS, DIMENSION(0:), INTENT(IN) :: p2p
1522 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), &
1523 : INTENT(IN) :: yzp
1524 : INTEGER, CONTIGUOUS, DIMENSION(0:), INTENT(IN) :: nray
1525 : INTEGER, DIMENSION(:, :, 0:), INTENT(IN) :: bo
1526 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1527 : INTENT(INOUT) :: sb
1528 : TYPE(fft_scratch_type), INTENT(IN) :: fft_scratch
1529 :
1530 : CHARACTER(len=*), PARAMETER :: routineN = 'yz_to_x'
1531 :
1532 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1533 1771530 : POINTER :: rr
1534 : COMPLEX(KIND=sp), CONTIGUOUS, DIMENSION(:, :), &
1535 1771530 : POINTER :: ss, tt
1536 : INTEGER :: handle, ip, ir, ix, ixx, iy, iz, mpr, &
1537 : nm, np, nr, nx
1538 1771530 : INTEGER, CONTIGUOUS, DIMENSION(:), POINTER :: rcount, rdispl, scount, sdispl
1539 :
1540 1771530 : CALL timeset(routineN, handle)
1541 :
1542 1771530 : np = SIZE(p2p)
1543 1771530 : mpr = p2p(my_pos)
1544 1771530 : scount => fft_scratch%scount
1545 1771530 : rcount => fft_scratch%rcount
1546 1771530 : sdispl => fft_scratch%sdispl
1547 1771530 : rdispl => fft_scratch%rdispl
1548 :
1549 1771530 : IF (alltoall_sgl) THEN
1550 238 : ss => fft_scratch%ss
1551 238 : tt => fft_scratch%tt
1552 238 : ss = 0._sp
1553 238 : tt = 0._sp
1554 : ELSE
1555 1771292 : rr => fft_scratch%rr
1556 : END IF
1557 :
1558 1771530 : nx = bo(2, 1, mpr) - bo(1, 1, mpr) + 1
1559 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(2) &
1560 : !$OMP PRIVATE(ip, ixx, ir, iy, iz, ix) &
1561 1771530 : !$OMP SHARED(np,nray,nx,alltoall_sgl,yzp,tb,tt,rr)
1562 : DO ip = 0, np - 1
1563 : DO ix = 1, nx
1564 : ixx = nray(ip)*(ix - 1)
1565 : IF (alltoall_sgl) THEN
1566 : DO ir = 1, nray(ip)
1567 : iy = yzp(1, ir, ip)
1568 : iz = yzp(2, ir, ip)
1569 : tt(ir + ixx, ip) = CMPLX(tb(iy, iz, ix), KIND=sp)
1570 : END DO
1571 : ELSE
1572 : DO ir = 1, nray(ip)
1573 : iy = yzp(1, ir, ip)
1574 : iz = yzp(2, ir, ip)
1575 : rr(ir + ixx, ip) = tb(iy, iz, ix)
1576 : END DO
1577 : END IF
1578 : END DO
1579 : END DO
1580 : !$OMP END PARALLEL DO
1581 5314590 : nm = MAXVAL(nray(0:np - 1))
1582 1771530 : nr = nray(my_pos)
1583 : !$OMP PARALLEL DO DEFAULT(NONE), &
1584 : !$OMP PRIVATE(ix,nx), &
1585 1771530 : !$OMP SHARED(np,p2p,bo,rcount,rdispl,nr)
1586 : DO ip = 0, np - 1
1587 : ix = p2p(ip)
1588 : nx = bo(2, 1, ix) - bo(1, 1, ix) + 1
1589 : rcount(ip) = nr*nx
1590 : rdispl(ip) = nr*(bo(1, 1, ix) - 1)
1591 : END DO
1592 : !$OMP END PARALLEL DO
1593 1771530 : nx = bo(2, 1, mpr) - bo(1, 1, mpr) + 1
1594 : !$OMP PARALLEL DO DEFAULT(NONE), &
1595 : !$OMP PRIVATE(nr), &
1596 1771530 : !$OMP SHARED(np,nray,scount,sdispl,nx,nm)
1597 : DO ip = 0, np - 1
1598 : nr = nray(ip)
1599 : scount(ip) = nr*nx
1600 : sdispl(ip) = nm*nx*ip
1601 : END DO
1602 : !$OMP END PARALLEL DO
1603 :
1604 1771530 : IF (alltoall_sgl) THEN
1605 238 : CALL group%alltoall(tt, scount, sdispl, ss, rcount, rdispl)
1606 3703835 : sb = ss
1607 : ELSE
1608 1771292 : CALL group%alltoall(rr, scount, sdispl, sb, rcount, rdispl)
1609 : END IF
1610 :
1611 1771530 : CALL timestop(handle)
1612 :
1613 1771530 : END SUBROUTINE yz_to_x
1614 :
1615 : ! **************************************************************************************************
1616 : !> \brief ...
1617 : !> \param sb ...
1618 : !> \param group ...
1619 : !> \param dims ...
1620 : !> \param my_pos ...
1621 : !> \param p2p ...
1622 : !> \param yzp ...
1623 : !> \param nray ...
1624 : !> \param bo ...
1625 : !> \param tb ...
1626 : !> \param fft_scratch ...
1627 : !> \par History
1628 : !> 15. Feb. 2006 : single precision all_to_all
1629 : !> \author JGH (18-Jan-2001)
1630 : ! **************************************************************************************************
1631 0 : SUBROUTINE yz_to_xz(sb, group, dims, my_pos, p2p, yzp, nray, bo, tb, fft_scratch)
1632 :
1633 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1634 : INTENT(IN) :: sb
1635 :
1636 : CLASS(mp_comm_type), INTENT(IN) :: group
1637 : INTEGER, DIMENSION(2), INTENT(IN) :: dims
1638 : INTEGER, INTENT(IN) :: my_pos
1639 : INTEGER, CONTIGUOUS, DIMENSION(0:), INTENT(IN) :: p2p
1640 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), INTENT(IN) :: yzp
1641 : INTEGER, CONTIGUOUS, DIMENSION(0:), INTENT(IN) :: nray
1642 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), INTENT(IN) :: bo
1643 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT), CONTIGUOUS :: tb
1644 : TYPE(fft_scratch_type), INTENT(INOUT) :: fft_scratch
1645 :
1646 : CHARACTER(len=*), PARAMETER :: routineN = 'yz_to_xz'
1647 :
1648 0 : COMPLEX(KIND=dp), DIMENSION(:), POINTER, CONTIGUOUS :: xzbuf, yzbuf
1649 0 : COMPLEX(KIND=sp), DIMENSION(:), POINTER, CONTIGUOUS :: xzbuf_sgl, yzbuf_sgl
1650 : INTEGER :: handle, icrs, ip, ipl, ipr, ir, ix, iz, &
1651 : jj, jx, jy, jz, myx, myz, np, npx, &
1652 : npz, nx, nz, rs_pos
1653 0 : INTEGER, DIMENSION(:), POINTER, CONTIGUOUS :: pzcoord, rcount, rdispl, scount, sdispl, &
1654 0 : xcor, zcor
1655 0 : INTEGER, DIMENSION(:, :), CONTIGUOUS, POINTER :: pgrid
1656 :
1657 0 : CALL timeset(routineN, handle)
1658 :
1659 0 : np = SIZE(p2p)
1660 :
1661 0 : rs_pos = p2p(my_pos)
1662 :
1663 0 : IF (alltoall_sgl) THEN
1664 0 : yzbuf_sgl => fft_scratch%yzbuf_sgl
1665 0 : xzbuf_sgl => fft_scratch%xzbuf_sgl
1666 : ELSE
1667 0 : yzbuf => fft_scratch%yzbuf
1668 0 : xzbuf => fft_scratch%xzbuf
1669 : END IF
1670 0 : npx = dims(1)
1671 0 : npz = dims(2)
1672 0 : pgrid => fft_scratch%pgrid
1673 0 : xcor => fft_scratch%xcor
1674 0 : zcor => fft_scratch%zcor
1675 0 : pzcoord => fft_scratch%pzcoord
1676 0 : scount => fft_scratch%scount
1677 0 : rcount => fft_scratch%rcount
1678 0 : sdispl => fft_scratch%sdispl
1679 0 : rdispl => fft_scratch%rdispl
1680 :
1681 0 : nx = SIZE(sb, 2)
1682 :
1683 : ! If the send and recv counts are not already cached, then
1684 : ! calculate and store them
1685 0 : IF (fft_scratch%in == 0) THEN
1686 :
1687 0 : scount = 0
1688 :
1689 0 : DO ix = 0, npx - 1
1690 0 : ip = pgrid(ix, 0)
1691 0 : xcor(bo(1, 1, ip):bo(2, 1, ip)) = ix
1692 : END DO
1693 0 : DO iz = 0, npz - 1
1694 0 : ip = pgrid(0, iz)
1695 0 : zcor(bo(1, 3, ip):bo(2, 3, ip)) = iz
1696 : END DO
1697 0 : DO jx = 1, nx
1698 0 : IF (alltoall_sgl) THEN
1699 0 : DO ir = 1, nray(my_pos)
1700 0 : jy = yzp(1, ir, my_pos)
1701 0 : jz = yzp(2, ir, my_pos)
1702 0 : ip = pgrid(xcor(jx), zcor(jz))
1703 0 : scount(ip) = scount(ip) + 1
1704 : END DO
1705 : ELSE
1706 0 : DO ir = 1, nray(my_pos)
1707 0 : jy = yzp(1, ir, my_pos)
1708 0 : jz = yzp(2, ir, my_pos)
1709 0 : ip = pgrid(xcor(jx), zcor(jz))
1710 0 : scount(ip) = scount(ip) + 1
1711 : END DO
1712 : END IF
1713 : END DO
1714 :
1715 0 : CALL group%alltoall(scount, rcount, 1)
1716 0 : fft_scratch%yzcount = scount
1717 0 : fft_scratch%xzcount = rcount
1718 :
1719 : ! Work out the correct displacements in the buffers
1720 0 : sdispl(0) = 0
1721 0 : rdispl(0) = 0
1722 0 : DO ip = 1, np - 1
1723 0 : sdispl(ip) = sdispl(ip - 1) + scount(ip - 1)
1724 0 : rdispl(ip) = rdispl(ip - 1) + rcount(ip - 1)
1725 : END DO
1726 :
1727 0 : fft_scratch%yzdispl = sdispl
1728 0 : fft_scratch%xzdispl = rdispl
1729 :
1730 0 : icrs = 0
1731 0 : DO ip = 0, np - 1
1732 0 : IF (scount(ip) /= 0) icrs = icrs + 1
1733 0 : IF (rcount(ip) /= 0) icrs = icrs + 1
1734 : END DO
1735 0 : CALL group%sum(icrs)
1736 0 : fft_scratch%rsratio = REAL(icrs, KIND=dp)/(REAL(2*np, KIND=dp)*REAL(np, KIND=dp))
1737 :
1738 0 : fft_scratch%in = 1
1739 : ELSE
1740 0 : scount = fft_scratch%yzcount
1741 0 : rcount = fft_scratch%xzcount
1742 0 : sdispl = fft_scratch%yzdispl
1743 0 : rdispl = fft_scratch%xzdispl
1744 : END IF
1745 :
1746 : ! Do the actual packing
1747 : !$OMP PARALLEL DO DEFAULT(NONE), &
1748 : !$OMP PRIVATE(ipl,jj,nx,ir,jx,jy,jz),&
1749 : !$OMP SHARED(np,p2p,pzcoord,bo,nray,yzp,zcor),&
1750 : !$OMP SHARED(yzbuf,sb,scount,sdispl,my_pos),&
1751 0 : !$OMP SHARED(yzbuf_sgl,alltoall_sgl)
1752 : DO ip = 0, np - 1
1753 : IF (scount(ip) == 0) CYCLE
1754 : ipl = p2p(ip)
1755 : jj = 0
1756 : nx = bo(2, 1, ipl) - bo(1, 1, ipl) + 1
1757 : DO ir = 1, nray(my_pos)
1758 : jz = yzp(2, ir, my_pos)
1759 : IF (zcor(jz) == pzcoord(ipl)) THEN
1760 : jj = jj + 1
1761 : jy = yzp(1, ir, my_pos)
1762 : IF (alltoall_sgl) THEN
1763 : DO jx = 0, nx - 1
1764 : yzbuf_sgl(sdispl(ip) + jj + jx*scount(ip)/nx) = CMPLX(sb(ir, jx + bo(1, 1, ipl)), KIND=sp)
1765 : END DO
1766 : ELSE
1767 : DO jx = 0, nx - 1
1768 : yzbuf(sdispl(ip) + jj + jx*scount(ip)/nx) = sb(ir, jx + bo(1, 1, ipl))
1769 : END DO
1770 : END IF
1771 : END IF
1772 : END DO
1773 : END DO
1774 : !$OMP END PARALLEL DO
1775 :
1776 0 : IF (alltoall_sgl) THEN
1777 0 : CALL group%alltoall(yzbuf_sgl, scount, sdispl, xzbuf_sgl, rcount, rdispl)
1778 : ELSE
1779 0 : IF (fft_scratch%rsratio < ratio_sparse_alltoall) THEN
1780 0 : CALL sparse_alltoall(yzbuf, scount, sdispl, xzbuf, rcount, rdispl, group)
1781 : ELSE
1782 0 : CALL group%alltoall(yzbuf, scount, sdispl, xzbuf, rcount, rdispl)
1783 : END IF
1784 : END IF
1785 :
1786 0 : myx = fft_scratch%sizes%r_pos(1)
1787 0 : myz = fft_scratch%sizes%r_pos(2)
1788 0 : nz = bo(2, 3, rs_pos) - bo(1, 3, rs_pos) + 1
1789 :
1790 : !$OMP PARALLEL DO DEFAULT(NONE), &
1791 : !$OMP PRIVATE(ipr,jj,ir,jx,jy,jz),&
1792 : !$OMP SHARED(tb,np,p2p,bo,rs_pos,nray),&
1793 : !$OMP SHARED(yzp,alltoall_sgl,zcor,myz),&
1794 0 : !$OMP SHARED(xzbuf,xzbuf_sgl,nz,rdispl)
1795 : DO ip = 0, np - 1
1796 : ipr = p2p(ip)
1797 : jj = 0
1798 : DO jx = 0, bo(2, 1, rs_pos) - bo(1, 1, rs_pos)
1799 : DO ir = 1, nray(ip)
1800 : jz = yzp(2, ir, ip)
1801 : IF (alltoall_sgl) THEN
1802 : IF (zcor(jz) == myz) THEN
1803 : jj = jj + 1
1804 : jy = yzp(1, ir, ip)
1805 : jz = jz - bo(1, 3, rs_pos) + 1
1806 : tb(jy, jz + jx*nz) = xzbuf_sgl(jj + rdispl(ipr))
1807 : END IF
1808 : ELSE
1809 : IF (zcor(jz) == myz) THEN
1810 : jj = jj + 1
1811 : jy = yzp(1, ir, ip)
1812 : jz = jz - bo(1, 3, rs_pos) + 1
1813 : tb(jy, jz + jx*nz) = xzbuf(jj + rdispl(ipr))
1814 : END IF
1815 : END IF
1816 : END DO
1817 : END DO
1818 : END DO
1819 : !$OMP END PARALLEL DO
1820 :
1821 0 : CALL timestop(handle)
1822 :
1823 0 : END SUBROUTINE yz_to_xz
1824 :
1825 : ! **************************************************************************************************
1826 : !> \brief ...
1827 : !> \param sb ...
1828 : !> \param group ...
1829 : !> \param dims ...
1830 : !> \param my_pos ...
1831 : !> \param p2p ...
1832 : !> \param yzp ...
1833 : !> \param nray ...
1834 : !> \param bo ...
1835 : !> \param tb ...
1836 : !> \param fft_scratch ...
1837 : !> \par History
1838 : !> 15. Feb. 2006 : single precision all_to_all
1839 : !> \author JGH (19-Jan-2001)
1840 : ! **************************************************************************************************
1841 0 : SUBROUTINE xz_to_yz(sb, group, dims, my_pos, p2p, yzp, nray, bo, tb, fft_scratch)
1842 :
1843 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1844 : INTENT(IN) :: sb
1845 :
1846 : CLASS(mp_comm_type), INTENT(IN) :: group
1847 : INTEGER, DIMENSION(2), INTENT(IN) :: dims
1848 : INTEGER, INTENT(IN) :: my_pos
1849 : INTEGER, CONTIGUOUS, DIMENSION(0:), INTENT(IN) :: p2p
1850 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), INTENT(IN) :: yzp
1851 : INTEGER, CONTIGUOUS, DIMENSION(0:), INTENT(IN) :: nray
1852 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), INTENT(IN) :: bo
1853 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT), CONTIGUOUS :: tb
1854 : TYPE(fft_scratch_type), INTENT(INOUT) :: fft_scratch
1855 :
1856 : CHARACTER(len=*), PARAMETER :: routineN = 'xz_to_yz'
1857 :
1858 0 : COMPLEX(KIND=dp), DIMENSION(:), POINTER, CONTIGUOUS :: xzbuf, yzbuf
1859 0 : COMPLEX(KIND=sp), DIMENSION(:), POINTER, CONTIGUOUS :: xzbuf_sgl, yzbuf_sgl
1860 : INTEGER :: handle, icrs, ip, ipl, ir, ix, ixx, iz, &
1861 : jj, jx, jy, jz, mp, myx, myz, np, npx, &
1862 : npz, nx, nz
1863 0 : INTEGER, DIMENSION(:), POINTER, CONTIGUOUS :: pzcoord, rcount, rdispl, scount, sdispl, &
1864 0 : xcor, zcor
1865 0 : INTEGER, DIMENSION(:, :), CONTIGUOUS, POINTER :: pgrid
1866 :
1867 0 : CALL timeset(routineN, handle)
1868 :
1869 0 : np = SIZE(p2p)
1870 :
1871 0 : IF (alltoall_sgl) THEN
1872 0 : yzbuf_sgl => fft_scratch%yzbuf_sgl
1873 0 : xzbuf_sgl => fft_scratch%xzbuf_sgl
1874 : ELSE
1875 0 : yzbuf => fft_scratch%yzbuf
1876 0 : xzbuf => fft_scratch%xzbuf
1877 : END IF
1878 0 : npx = dims(1)
1879 0 : npz = dims(2)
1880 0 : pgrid => fft_scratch%pgrid
1881 0 : xcor => fft_scratch%xcor
1882 0 : zcor => fft_scratch%zcor
1883 0 : pzcoord => fft_scratch%pzcoord
1884 0 : scount => fft_scratch%scount
1885 0 : rcount => fft_scratch%rcount
1886 0 : sdispl => fft_scratch%sdispl
1887 0 : rdispl => fft_scratch%rdispl
1888 :
1889 : ! If the send and recv counts are not already cached, then
1890 : ! calculate and store them
1891 0 : IF (fft_scratch%in == 0) THEN
1892 :
1893 0 : rcount = 0
1894 0 : nx = MAXVAL(bo(2, 1, :))
1895 :
1896 0 : DO ix = 0, npx - 1
1897 0 : ip = pgrid(ix, 0)
1898 0 : xcor(bo(1, 1, ip):bo(2, 1, ip)) = ix
1899 : END DO
1900 0 : DO iz = 0, npz - 1
1901 0 : ip = pgrid(0, iz)
1902 0 : zcor(bo(1, 3, ip):bo(2, 3, ip)) = iz
1903 : END DO
1904 0 : DO jx = 1, nx
1905 0 : DO ir = 1, nray(my_pos)
1906 0 : jy = yzp(1, ir, my_pos)
1907 0 : jz = yzp(2, ir, my_pos)
1908 0 : ip = pgrid(xcor(jx), zcor(jz))
1909 0 : rcount(ip) = rcount(ip) + 1
1910 : END DO
1911 : END DO
1912 :
1913 0 : CALL group%alltoall(rcount, scount, 1)
1914 0 : fft_scratch%xzcount = scount
1915 0 : fft_scratch%yzcount = rcount
1916 :
1917 : ! Work out the correct displacements in the buffers
1918 0 : sdispl(0) = 0
1919 0 : rdispl(0) = 0
1920 0 : DO ip = 1, np - 1
1921 0 : sdispl(ip) = sdispl(ip - 1) + scount(ip - 1)
1922 0 : rdispl(ip) = rdispl(ip - 1) + rcount(ip - 1)
1923 : END DO
1924 :
1925 0 : fft_scratch%xzdispl = sdispl
1926 0 : fft_scratch%yzdispl = rdispl
1927 :
1928 0 : icrs = 0
1929 0 : DO ip = 0, np - 1
1930 0 : IF (scount(ip) /= 0) icrs = icrs + 1
1931 0 : IF (rcount(ip) /= 0) icrs = icrs + 1
1932 : END DO
1933 0 : CALL group%sum(icrs)
1934 0 : fft_scratch%rsratio = REAL(icrs, KIND=dp)/(REAL(2*np, KIND=dp)*REAL(np, KIND=dp))
1935 :
1936 0 : fft_scratch%in = 1
1937 : ELSE
1938 0 : scount = fft_scratch%xzcount
1939 0 : rcount = fft_scratch%yzcount
1940 0 : sdispl = fft_scratch%xzdispl
1941 0 : rdispl = fft_scratch%yzdispl
1942 : END IF
1943 :
1944 : ! Now do the actual packing
1945 0 : myx = fft_scratch%sizes%r_pos(1)
1946 0 : myz = fft_scratch%sizes%r_pos(2)
1947 0 : mp = p2p(my_pos)
1948 0 : nz = bo(2, 3, mp) - bo(1, 3, mp) + 1
1949 0 : nx = bo(2, 1, mp) - bo(1, 1, mp) + 1
1950 :
1951 : !$OMP PARALLEL DO DEFAULT(NONE), &
1952 : !$OMP PRIVATE(jj,ipl,ir,jx,jy,jz,ixx),&
1953 : !$OMP SHARED(np,p2p,nray,yzp,zcor,myz,bo,mp),&
1954 : !$OMP SHARED(alltoall_sgl,nx,scount,sdispl),&
1955 0 : !$OMP SHARED(xzbuf,xzbuf_sgl,sb,nz)
1956 : DO ip = 0, np - 1
1957 : jj = 0
1958 : ipl = p2p(ip)
1959 : DO ir = 1, nray(ip)
1960 : jz = yzp(2, ir, ip)
1961 : IF (zcor(jz) == myz) THEN
1962 : jj = jj + 1
1963 : jy = yzp(1, ir, ip)
1964 : jz = yzp(2, ir, ip) - bo(1, 3, mp) + 1
1965 : IF (alltoall_sgl) THEN
1966 : DO jx = 0, nx - 1
1967 : ixx = jj + jx*scount(ipl)/nx
1968 : xzbuf_sgl(ixx + sdispl(ipl)) = CMPLX(sb(jy, jz + jx*nz), KIND=sp)
1969 : END DO
1970 : ELSE
1971 : DO jx = 0, nx - 1
1972 : ixx = jj + jx*scount(ipl)/nx
1973 : xzbuf(ixx + sdispl(ipl)) = sb(jy, jz + jx*nz)
1974 : END DO
1975 : END IF
1976 : END IF
1977 : END DO
1978 : END DO
1979 : !$OMP END PARALLEL DO
1980 :
1981 0 : IF (alltoall_sgl) THEN
1982 0 : CALL group%alltoall(xzbuf_sgl, scount, sdispl, yzbuf_sgl, rcount, rdispl)
1983 : ELSE
1984 0 : IF (fft_scratch%rsratio < ratio_sparse_alltoall) THEN
1985 0 : CALL sparse_alltoall(xzbuf, scount, sdispl, yzbuf, rcount, rdispl, group)
1986 : ELSE
1987 0 : CALL group%alltoall(xzbuf, scount, sdispl, yzbuf, rcount, rdispl)
1988 : END IF
1989 : END IF
1990 :
1991 : !$OMP PARALLEL DO DEFAULT(NONE), &
1992 : !$OMP PRIVATE(ipl,jj,nx,ir,jx,jy,jz),&
1993 : !$OMP SHARED(p2p,pzcoord,bo,nray,my_pos,yzp),&
1994 : !$OMP SHARED(rcount,rdispl,tb,yzbuf,zcor),&
1995 0 : !$OMP SHARED(yzbuf_sgl,alltoall_sgl,np)
1996 : DO ip = 0, np - 1
1997 : IF (rcount(ip) == 0) CYCLE
1998 : ipl = p2p(ip)
1999 : jj = 0
2000 : nx = bo(2, 1, ipl) - bo(1, 1, ipl) + 1
2001 : DO ir = 1, nray(my_pos)
2002 : jz = yzp(2, ir, my_pos)
2003 : IF (zcor(jz) == pzcoord(ipl)) THEN
2004 : jj = jj + 1
2005 : jy = yzp(1, ir, my_pos)
2006 : IF (alltoall_sgl) THEN
2007 : DO jx = 0, nx - 1
2008 : tb(ir, jx + bo(1, 1, ipl)) = yzbuf_sgl(rdispl(ip) + jj + jx*rcount(ip)/nx)
2009 : END DO
2010 : ELSE
2011 : DO jx = 0, nx - 1
2012 : tb(ir, jx + bo(1, 1, ipl)) = yzbuf(rdispl(ip) + jj + jx*rcount(ip)/nx)
2013 : END DO
2014 : END IF
2015 : END IF
2016 : END DO
2017 : END DO
2018 : !$OMP END PARALLEL DO
2019 :
2020 0 : CALL timestop(handle)
2021 :
2022 0 : END SUBROUTINE xz_to_yz
2023 :
2024 : ! **************************************************************************************************
2025 : !> \brief ...
2026 : !> \param cin ...
2027 : !> \param boin ...
2028 : !> \param boout ...
2029 : !> \param sout ...
2030 : !> \param fft_scratch ...
2031 : !> \par History
2032 : !> none
2033 : !> \author JGH (20-Jan-2001)
2034 : ! **************************************************************************************************
2035 0 : SUBROUTINE cube_transpose_1(cin, boin, boout, sout, fft_scratch)
2036 :
2037 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2038 : INTENT(IN) :: cin
2039 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), &
2040 : INTENT(IN) :: boin, boout
2041 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2042 : INTENT(OUT) :: sout
2043 : TYPE(fft_scratch_type), INTENT(IN) :: fft_scratch
2044 :
2045 : CHARACTER(len=*), PARAMETER :: routineN = 'cube_transpose_1'
2046 :
2047 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2048 0 : POINTER :: rbuf
2049 : INTEGER :: handle, ip, ipl, ir, is, ixy, iz, mip, &
2050 : mz, np, nx, ny, nz
2051 0 : INTEGER, CONTIGUOUS, DIMENSION(:), POINTER :: rcount, rdispl, scount, sdispl
2052 0 : INTEGER, CONTIGUOUS, DIMENSION(:, :), POINTER :: pgrid
2053 : INTEGER, DIMENSION(2) :: dim, pos
2054 :
2055 0 : CALL timeset(routineN, handle)
2056 :
2057 0 : mip = fft_scratch%mip
2058 0 : dim = fft_scratch%dim
2059 0 : pos = fft_scratch%pos
2060 0 : scount => fft_scratch%scount
2061 0 : rcount => fft_scratch%rcount
2062 0 : sdispl => fft_scratch%sdispl
2063 0 : rdispl => fft_scratch%rdispl
2064 0 : pgrid => fft_scratch%pgcube
2065 0 : np = DIM(2)
2066 :
2067 0 : nx = boin(2, 1, mip) - boin(1, 1, mip) + 1
2068 0 : nz = boin(2, 3, mip) - boin(1, 3, mip) + 1
2069 :
2070 : !$OMP PARALLEL DO DEFAULT(NONE), &
2071 : !$OMP PRIVATE(ipl,ny), &
2072 0 : !$OMP SHARED(np,pgrid,boout,scount,sdispl,nx,nz)
2073 : DO ip = 0, np - 1
2074 : ipl = pgrid(ip, 2)
2075 : ny = boout(2, 2, ipl) - boout(1, 2, ipl) + 1
2076 : scount(ip) = nx*nz*ny
2077 : sdispl(ip) = nx*nz*(boout(1, 2, ipl) - 1)
2078 : END DO
2079 : !$OMP END PARALLEL DO
2080 0 : ny = boout(2, 2, mip) - boout(1, 2, mip) + 1
2081 0 : mz = MAXVAL(boin(2, 3, :) - boin(1, 3, :) + 1)
2082 : !$OMP PARALLEL DO DEFAULT(NONE), &
2083 : !$OMP PRIVATE(ipl,nz), &
2084 0 : !$OMP SHARED(np,pgrid,boin,nx,ny,rcount,rdispl,mz)
2085 : DO ip = 0, np - 1
2086 : ipl = pgrid(ip, 2)
2087 : nz = boin(2, 3, ipl) - boin(1, 3, ipl) + 1
2088 : rcount(ip) = nx*nz*ny
2089 : rdispl(ip) = nx*ny*mz*ip
2090 : END DO
2091 : !$OMP END PARALLEL DO
2092 :
2093 0 : rbuf => fft_scratch%rbuf1
2094 :
2095 0 : CALL fft_scratch%cart_sub_comm(2)%alltoall(cin, scount, sdispl, rbuf, rcount, rdispl)
2096 :
2097 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(2) &
2098 : !$OMP PRIVATE(ip,ipl,nz,iz,is,ir) &
2099 0 : !$OMP SHARED(nx,ny,np,pgrid,boin,sout,rbuf)
2100 : DO ixy = 1, nx*ny
2101 : DO ip = 0, np - 1
2102 : ipl = pgrid(ip, 2)
2103 : nz = boin(2, 3, ipl) - boin(1, 3, ipl) + 1
2104 : DO iz = 1, nz
2105 : is = boin(1, 3, ipl) + iz - 1
2106 : ir = iz + nz*(ixy - 1)
2107 : sout(is, ixy) = rbuf(ir, ip)
2108 : END DO
2109 : END DO
2110 : END DO
2111 : !$OMP END PARALLEL DO
2112 :
2113 0 : CALL timestop(handle)
2114 :
2115 0 : END SUBROUTINE cube_transpose_1
2116 :
2117 : ! **************************************************************************************************
2118 : !> \brief ...
2119 : !> \param cin ...
2120 : !> \param boin ...
2121 : !> \param boout ...
2122 : !> \param sout ...
2123 : !> \param fft_scratch ...
2124 : ! **************************************************************************************************
2125 0 : SUBROUTINE cube_transpose_2(cin, boin, boout, sout, fft_scratch)
2126 :
2127 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2128 : INTENT(IN) :: cin
2129 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), &
2130 : INTENT(IN) :: boin, boout
2131 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2132 : INTENT(OUT) :: sout
2133 : TYPE(fft_scratch_type), INTENT(IN) :: fft_scratch
2134 :
2135 : CHARACTER(len=*), PARAMETER :: routineN = 'cube_transpose_2'
2136 :
2137 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2138 0 : POINTER :: rbuf
2139 : INTEGER :: handle, ip, ipl, ir, ixy, iz, mip, mz, &
2140 : np, nx, ny, nz
2141 0 : INTEGER, CONTIGUOUS, DIMENSION(:), POINTER :: rcount, rdispl, scount, sdispl
2142 0 : INTEGER, CONTIGUOUS, DIMENSION(:, :), POINTER :: pgrid
2143 : INTEGER, DIMENSION(2) :: dim, pos
2144 : TYPE(mp_comm_type) :: sub_group
2145 :
2146 0 : CALL timeset(routineN, handle)
2147 :
2148 0 : sub_group = fft_scratch%cart_sub_comm(2)
2149 0 : mip = fft_scratch%mip
2150 0 : dim = fft_scratch%dim
2151 0 : pos = fft_scratch%pos
2152 0 : scount => fft_scratch%scount
2153 0 : rcount => fft_scratch%rcount
2154 0 : sdispl => fft_scratch%sdispl
2155 0 : rdispl => fft_scratch%rdispl
2156 0 : pgrid => fft_scratch%pgcube
2157 0 : np = DIM(2)
2158 :
2159 0 : nx = boin(2, 1, mip) - boin(1, 1, mip) + 1
2160 0 : ny = boin(2, 2, mip) - boin(1, 2, mip) + 1
2161 0 : mz = MAXVAL(boout(2, 3, :) - boout(1, 3, :) + 1)
2162 :
2163 0 : rbuf => fft_scratch%rbuf2
2164 :
2165 : !$OMP PARALLEL DEFAULT(NONE), &
2166 : !$OMP PRIVATE(ip,ipl,nz,iz,ir), &
2167 0 : !$OMP SHARED(nx,ny,np,pgrid,boout,rbuf,cin,scount,sdispl,mz)
2168 : !$OMP DO COLLAPSE(2)
2169 : DO ixy = 1, nx*ny
2170 : DO ip = 0, np - 1
2171 : ipl = pgrid(ip, 2)
2172 : nz = boout(2, 3, ipl) - boout(1, 3, ipl) + 1
2173 : DO iz = boout(1, 3, ipl), boout(2, 3, ipl)
2174 : ir = iz - boout(1, 3, ipl) + 1 + (ixy - 1)*nz
2175 : rbuf(ir, ip) = cin(iz, ixy)
2176 : END DO
2177 : END DO
2178 : END DO
2179 : !$OMP END DO
2180 : !$OMP DO
2181 : DO ip = 0, np - 1
2182 : ipl = pgrid(ip, 2)
2183 : nz = boout(2, 3, ipl) - boout(1, 3, ipl) + 1
2184 : scount(ip) = nx*ny*nz
2185 : sdispl(ip) = nx*ny*mz*ip
2186 : END DO
2187 : !$OMP END DO
2188 : !$OMP END PARALLEL
2189 0 : nz = boout(2, 3, mip) - boout(1, 3, mip) + 1
2190 : !$OMP PARALLEL DO DEFAULT(NONE), &
2191 : !$OMP PRIVATE(ipl,ny), &
2192 0 : !$OMP SHARED(np,pgrid,boin,nx,nz,rcount,rdispl)
2193 : DO ip = 0, np - 1
2194 : ipl = pgrid(ip, 2)
2195 : ny = boin(2, 2, ipl) - boin(1, 2, ipl) + 1
2196 : rcount(ip) = nx*ny*nz
2197 : rdispl(ip) = nx*nz*(boin(1, 2, ipl) - 1)
2198 : END DO
2199 : !$OMP END PARALLEL DO
2200 :
2201 0 : CALL sub_group%alltoall(rbuf, scount, sdispl, sout, rcount, rdispl)
2202 :
2203 0 : CALL timestop(handle)
2204 :
2205 0 : END SUBROUTINE cube_transpose_2
2206 :
2207 : ! **************************************************************************************************
2208 : !> \brief ...
2209 : !> \param cin ...
2210 : !> \param boin ...
2211 : !> \param boout ...
2212 : !> \param sout ...
2213 : !> \param fft_scratch ...
2214 : ! **************************************************************************************************
2215 0 : SUBROUTINE cube_transpose_3(cin, boin, boout, sout, fft_scratch)
2216 :
2217 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2218 : INTENT(IN) :: cin
2219 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), &
2220 : INTENT(IN) :: boin, boout
2221 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2222 : INTENT(OUT) :: sout
2223 : TYPE(fft_scratch_type), INTENT(IN) :: fft_scratch
2224 :
2225 : CHARACTER(len=*), PARAMETER :: routineN = 'cube_transpose_3'
2226 :
2227 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2228 0 : POINTER :: rbuf
2229 : INTEGER :: handle, ip, ipl, ir, is, ixz, iy, lb, &
2230 : mip, my, my_id, np, num_threads, nx, &
2231 : ny, nz, ub
2232 0 : INTEGER, CONTIGUOUS, DIMENSION(:), POINTER :: rcount, rdispl, scount, sdispl
2233 0 : INTEGER, CONTIGUOUS, DIMENSION(:, :), POINTER :: pgrid
2234 : INTEGER, DIMENSION(2) :: dim, pos
2235 : TYPE(mp_comm_type) :: sub_group
2236 :
2237 0 : CALL timeset(routineN, handle)
2238 :
2239 0 : sub_group = fft_scratch%cart_sub_comm(1)
2240 0 : mip = fft_scratch%mip
2241 0 : dim = fft_scratch%dim
2242 0 : pos = fft_scratch%pos
2243 0 : np = DIM(1)
2244 0 : scount => fft_scratch%scount
2245 0 : rcount => fft_scratch%rcount
2246 0 : sdispl => fft_scratch%sdispl
2247 0 : rdispl => fft_scratch%rdispl
2248 0 : pgrid => fft_scratch%pgcube
2249 :
2250 0 : ny = boin(2, 2, mip) - boin(1, 2, mip) + 1
2251 0 : nz = boin(2, 3, mip) - boin(1, 3, mip) + 1
2252 : !$OMP PARALLEL DO DEFAULT(NONE), &
2253 : !$OMP PRIVATE(ipl, nx), &
2254 0 : !$OMP SHARED(np,pgrid,boout,ny,nz,scount,sdispl)
2255 : DO ip = 0, np - 1
2256 : ipl = pgrid(ip, 1)
2257 : nx = boout(2, 1, ipl) - boout(1, 1, ipl) + 1
2258 : scount(ip) = nx*nz*ny
2259 : sdispl(ip) = ny*nz*(boout(1, 1, ipl) - 1)
2260 : END DO
2261 : !$OMP END PARALLEL DO
2262 0 : nx = boout(2, 1, mip) - boout(1, 1, mip) + 1
2263 0 : my = MAXVAL(boin(2, 2, :) - boin(1, 2, :) + 1)
2264 : !$OMP PARALLEL DO DEFAULT(NONE), &
2265 : !$OMP PRIVATE(ipl, ny), &
2266 0 : !$OMP SHARED(np,pgrid,boin,nx,nz,my,rcount,rdispl)
2267 : DO ip = 0, np - 1
2268 : ipl = pgrid(ip, 1)
2269 : ny = boin(2, 2, ipl) - boin(1, 2, ipl) + 1
2270 : rcount(ip) = nx*nz*ny
2271 : rdispl(ip) = nx*my*nz*ip
2272 : END DO
2273 : !$OMP END PARALLEL DO
2274 :
2275 0 : rbuf => fft_scratch%rbuf3
2276 0 : num_threads = 1
2277 0 : my_id = 0
2278 : !$OMP PARALLEL DEFAULT(NONE), &
2279 : !$OMP PRIVATE(NUM_THREADS, my_id, lb, ub) &
2280 0 : !$OMP SHARED(rbuf)
2281 : !$ num_threads = MIN(omp_get_max_threads(), SIZE(rbuf, 2))
2282 : !$ my_id = omp_get_thread_num()
2283 : IF (my_id < num_threads) THEN
2284 : lb = (SIZE(rbuf, 2)*my_id)/num_threads
2285 : ub = (SIZE(rbuf, 2)*(my_id + 1))/num_threads - 1
2286 : rbuf(:, lb:ub) = 0.0_dp
2287 : END IF
2288 : !$OMP END PARALLEL
2289 :
2290 0 : CALL sub_group%alltoall(cin, scount, sdispl, rbuf, rcount, rdispl)
2291 :
2292 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(2) &
2293 : !$OMP PRIVATE(ip,ipl,ny,iy,is,ir) &
2294 0 : !$OMP SHARED(nx,nz,np,pgrid,boin,rbuf,sout)
2295 : DO ixz = 1, nx*nz
2296 : DO ip = 0, np - 1
2297 : ipl = pgrid(ip, 1)
2298 : ny = boin(2, 2, ipl) - boin(1, 2, ipl) + 1
2299 : DO iy = 1, ny
2300 : is = boin(1, 2, ipl) + iy - 1
2301 : ir = iy + ny*(ixz - 1)
2302 : sout(is, ixz) = rbuf(ir, ip)
2303 : END DO
2304 : END DO
2305 : END DO
2306 : !$OMP END PARALLEL DO
2307 :
2308 0 : CALL timestop(handle)
2309 :
2310 0 : END SUBROUTINE cube_transpose_3
2311 :
2312 : ! **************************************************************************************************
2313 : !> \brief ...
2314 : !> \param cin ...
2315 : !> \param boin ...
2316 : !> \param boout ...
2317 : !> \param sout ...
2318 : !> \param fft_scratch ...
2319 : ! **************************************************************************************************
2320 0 : SUBROUTINE cube_transpose_4(cin, boin, boout, sout, fft_scratch)
2321 :
2322 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2323 : INTENT(IN) :: cin
2324 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), &
2325 : INTENT(IN) :: boin, boout
2326 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2327 : INTENT(OUT) :: sout
2328 : TYPE(fft_scratch_type), INTENT(IN) :: fft_scratch
2329 :
2330 : CHARACTER(len=*), PARAMETER :: routineN = 'cube_transpose_4'
2331 :
2332 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2333 0 : POINTER :: rbuf
2334 : INTEGER :: handle, ip, ipl, ir, iy, izx, lb, mip, &
2335 : my, my_id, np, num_threads, nx, ny, &
2336 : nz, ub
2337 0 : INTEGER, CONTIGUOUS, DIMENSION(:), POINTER :: rcount, rdispl, scount, sdispl
2338 0 : INTEGER, CONTIGUOUS, DIMENSION(:, :), POINTER :: pgrid
2339 : INTEGER, DIMENSION(2) :: dim, pos
2340 : TYPE(mp_comm_type) :: sub_group
2341 :
2342 0 : CALL timeset(routineN, handle)
2343 :
2344 0 : sub_group = fft_scratch%cart_sub_comm(1)
2345 0 : mip = fft_scratch%mip
2346 0 : dim = fft_scratch%dim
2347 0 : pos = fft_scratch%pos
2348 0 : np = DIM(1)
2349 0 : scount => fft_scratch%scount
2350 0 : rcount => fft_scratch%rcount
2351 0 : sdispl => fft_scratch%sdispl
2352 0 : rdispl => fft_scratch%rdispl
2353 0 : pgrid => fft_scratch%pgcube
2354 :
2355 0 : nx = boin(2, 1, mip) - boin(1, 1, mip) + 1
2356 0 : nz = boin(2, 3, mip) - boin(1, 3, mip) + 1
2357 0 : my = MAXVAL(boout(2, 2, :) - boout(1, 2, :) + 1)
2358 :
2359 0 : rbuf => fft_scratch%rbuf4
2360 0 : num_threads = 1
2361 0 : my_id = 0
2362 : !$OMP PARALLEL DEFAULT(NONE), &
2363 : !$OMP PRIVATE(NUM_THREADS,my_id,lb,ub,ip,ipl,ny,iy,ir), &
2364 0 : !$OMP SHARED(rbuf,nz,nx,np,pgrid,boout,cin,my,scount,sdispl)
2365 : !$ num_threads = MIN(omp_get_max_threads(), SIZE(rbuf, 2))
2366 : !$ my_id = omp_get_thread_num()
2367 : IF (my_id < num_threads) THEN
2368 : lb = (SIZE(rbuf, 2)*my_id)/num_threads
2369 : ub = (SIZE(rbuf, 2)*(my_id + 1))/num_threads - 1
2370 : rbuf(:, lb:ub) = 0.0_dp
2371 : END IF
2372 : !$OMP BARRIER
2373 :
2374 : !$OMP DO COLLAPSE(2)
2375 : DO izx = 1, nz*nx
2376 : DO ip = 0, np - 1
2377 : ipl = pgrid(ip, 1)
2378 : ny = boout(2, 2, ipl) - boout(1, 2, ipl) + 1
2379 : DO iy = boout(1, 2, ipl), boout(2, 2, ipl)
2380 : ir = iy - boout(1, 2, ipl) + 1 + (izx - 1)*ny
2381 : rbuf(ir, ip) = cin(iy, izx)
2382 : END DO
2383 : END DO
2384 : END DO
2385 : !$OMP END DO
2386 : !$OMP DO
2387 : DO ip = 0, np - 1
2388 : ipl = pgrid(ip, 1)
2389 : ny = boout(2, 2, ipl) - boout(1, 2, ipl) + 1
2390 : scount(ip) = nx*ny*nz
2391 : sdispl(ip) = nx*nz*my*ip
2392 : END DO
2393 : !$OMP END DO
2394 : !$OMP END PARALLEL
2395 0 : ny = boout(2, 2, mip) - boout(1, 2, mip) + 1
2396 : !$OMP PARALLEL DO DEFAULT(NONE), &
2397 : !$OMP PRIVATE(ipl,nx), &
2398 0 : !$OMP SHARED(np,pgrid,boin,rcount,rdispl,ny,nz)
2399 : DO ip = 0, np - 1
2400 : ipl = pgrid(ip, 1)
2401 : nx = boin(2, 1, ipl) - boin(1, 1, ipl) + 1
2402 : rcount(ip) = nx*ny*nz
2403 : rdispl(ip) = ny*nz*(boin(1, 1, ipl) - 1)
2404 : END DO
2405 : !$OMP END PARALLEL DO
2406 :
2407 0 : CALL sub_group%alltoall(rbuf, scount, sdispl, sout, rcount, rdispl)
2408 :
2409 0 : CALL timestop(handle)
2410 :
2411 0 : END SUBROUTINE cube_transpose_4
2412 :
2413 : ! **************************************************************************************************
2414 : !> \brief ...
2415 : !> \param cin ...
2416 : !> \param group ...
2417 : !> \param boin ...
2418 : !> \param boout ...
2419 : !> \param sout ...
2420 : !> \param fft_scratch ...
2421 : ! **************************************************************************************************
2422 104 : SUBROUTINE cube_transpose_5(cin, group, boin, boout, sout, fft_scratch)
2423 :
2424 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2425 : INTENT(IN) :: cin
2426 :
2427 : CLASS(mp_comm_type), INTENT(IN) :: group
2428 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), INTENT(IN) :: boin, boout
2429 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT), CONTIGUOUS :: sout
2430 : TYPE(fft_scratch_type), INTENT(IN) :: fft_scratch
2431 :
2432 : CHARACTER(len=*), PARAMETER :: routineN = 'cube_transpose_5'
2433 :
2434 104 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER, CONTIGUOUS :: rbuf
2435 : INTEGER :: handle, ip, ir, is, ixz, iy, lb, mip, &
2436 : my, my_id, np, num_threads, nx, ny, &
2437 : nz, ub
2438 104 : INTEGER, DIMENSION(:), POINTER, CONTIGUOUS :: rcount, rdispl, scount, sdispl
2439 :
2440 104 : CALL timeset(routineN, handle)
2441 :
2442 104 : np = fft_scratch%sizes%numtask
2443 104 : mip = fft_scratch%mip
2444 104 : scount => fft_scratch%scount
2445 104 : rcount => fft_scratch%rcount
2446 104 : sdispl => fft_scratch%sdispl
2447 104 : rdispl => fft_scratch%rdispl
2448 :
2449 104 : ny = boin(2, 2, mip) - boin(1, 2, mip) + 1
2450 104 : nz = boin(2, 3, mip) - boin(1, 3, mip) + 1
2451 : !$OMP PARALLEL DO DEFAULT(NONE), &
2452 : !$OMP PRIVATE(nx), &
2453 104 : !$OMP SHARED(np,boout,ny,nz,scount,sdispl)
2454 : DO ip = 0, np - 1
2455 : nx = boout(2, 1, ip) - boout(1, 1, ip) + 1
2456 : scount(ip) = nx*nz*ny
2457 : sdispl(ip) = ny*nz*(boout(1, 1, ip) - 1)
2458 : END DO
2459 : !$OMP END PARALLEL DO
2460 104 : nx = boout(2, 1, mip) - boout(1, 1, mip) + 1
2461 312 : my = MAXVAL(boin(2, 2, :) - boin(1, 2, :) + 1)
2462 : !$OMP PARALLEL DO DEFAULT(NONE), &
2463 : !$OMP PRIVATE(ny), &
2464 104 : !$OMP SHARED(np,boin,nx,nz,rcount,rdispl,my)
2465 : DO ip = 0, np - 1
2466 : ny = boin(2, 2, ip) - boin(1, 2, ip) + 1
2467 : rcount(ip) = nx*nz*ny
2468 : rdispl(ip) = nx*my*nz*ip
2469 : END DO
2470 : !$OMP END PARALLEL DO
2471 :
2472 104 : rbuf => fft_scratch%rbuf5
2473 104 : num_threads = 1
2474 104 : my_id = 0
2475 : !$OMP PARALLEL DEFAULT(NONE), &
2476 : !$OMP PRIVATE(NUM_THREADS, my_id, lb, ub), &
2477 104 : !$OMP SHARED(rbuf)
2478 : !$ num_threads = MIN(omp_get_max_threads(), SIZE(rbuf, 2))
2479 : !$ my_id = omp_get_thread_num()
2480 : IF (my_id < num_threads) THEN
2481 : lb = (SIZE(rbuf, 2)*my_id)/num_threads
2482 : ub = (SIZE(rbuf, 2)*(my_id + 1))/num_threads - 1
2483 : rbuf(:, lb:ub) = 0.0_dp
2484 : END IF
2485 : !$OMP END PARALLEL
2486 :
2487 104 : CALL group%alltoall(cin, scount, sdispl, rbuf, rcount, rdispl)
2488 :
2489 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(2) &
2490 : !$OMP PRIVATE(ip,ny,iy,is,ir) &
2491 104 : !$OMP SHARED(nx,nz,np,boin,sout,rbuf)
2492 : DO ixz = 1, nx*nz
2493 : DO ip = 0, np - 1
2494 : ny = boin(2, 2, ip) - boin(1, 2, ip) + 1
2495 : DO iy = 1, ny
2496 : is = boin(1, 2, ip) + iy - 1
2497 : ir = iy + ny*(ixz - 1)
2498 : sout(is, ixz) = rbuf(ir, ip)
2499 : END DO
2500 : END DO
2501 : END DO
2502 : !$OMP END PARALLEL DO
2503 :
2504 104 : CALL timestop(handle)
2505 :
2506 104 : END SUBROUTINE cube_transpose_5
2507 :
2508 : ! **************************************************************************************************
2509 : !> \brief ...
2510 : !> \param cin ...
2511 : !> \param group ...
2512 : !> \param boin ...
2513 : !> \param boout ...
2514 : !> \param sout ...
2515 : !> \param fft_scratch ...
2516 : ! **************************************************************************************************
2517 104 : SUBROUTINE cube_transpose_6(cin, group, boin, boout, sout, fft_scratch)
2518 :
2519 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
2520 : INTENT(IN) :: cin
2521 :
2522 : CLASS(mp_comm_type), INTENT(IN) :: group
2523 : INTEGER, CONTIGUOUS, DIMENSION(:, :, 0:), INTENT(IN) :: boin, boout
2524 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT), CONTIGUOUS :: sout
2525 : TYPE(fft_scratch_type), INTENT(IN) :: fft_scratch
2526 :
2527 : CHARACTER(len=*), PARAMETER :: routineN = 'cube_transpose_6'
2528 :
2529 104 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER, CONTIGUOUS :: rbuf
2530 : INTEGER :: handle, ip, ir, iy, izx, lb, mip, my, &
2531 : my_id, np, num_threads, nx, ny, nz, ub
2532 104 : INTEGER, DIMENSION(:), POINTER, CONTIGUOUS :: rcount, rdispl, scount, sdispl
2533 :
2534 104 : CALL timeset(routineN, handle)
2535 :
2536 104 : np = fft_scratch%sizes%numtask
2537 104 : mip = fft_scratch%mip
2538 104 : scount => fft_scratch%scount
2539 104 : rcount => fft_scratch%rcount
2540 104 : sdispl => fft_scratch%sdispl
2541 104 : rdispl => fft_scratch%rdispl
2542 :
2543 104 : nx = boin(2, 1, mip) - boin(1, 1, mip) + 1
2544 104 : nz = boin(2, 3, mip) - boin(1, 3, mip) + 1
2545 312 : my = MAXVAL(boout(2, 2, :) - boout(1, 2, :) + 1)
2546 :
2547 104 : rbuf => fft_scratch%rbuf5
2548 104 : num_threads = 1
2549 104 : my_id = 0
2550 : !$OMP PARALLEL DEFAULT(NONE), &
2551 : !$OMP PRIVATE(NUM_THREADS,my_id,lb,ub,ip,ny,iy,ir), &
2552 104 : !$OMP SHARED(rbuf,nx,nz,np,boout,cin,my,scount,sdispl)
2553 : !$ num_threads = MIN(omp_get_max_threads(), SIZE(rbuf, 2))
2554 : !$ my_id = omp_get_thread_num()
2555 : IF (my_id < num_threads) THEN
2556 : lb = (SIZE(rbuf, 2)*my_id)/num_threads
2557 : ub = (SIZE(rbuf, 2)*(my_id + 1))/num_threads - 1
2558 : rbuf(:, lb:ub) = 0.0_dp
2559 : END IF
2560 : !$OMP BARRIER
2561 :
2562 : !$OMP DO COLLAPSE(2)
2563 : DO izx = 1, nz*nx
2564 : DO ip = 0, np - 1
2565 : ny = boout(2, 2, ip) - boout(1, 2, ip) + 1
2566 : DO iy = boout(1, 2, ip), boout(2, 2, ip)
2567 : ir = iy - boout(1, 2, ip) + 1 + (izx - 1)*ny
2568 : rbuf(ir, ip) = cin(iy, izx)
2569 : END DO
2570 : END DO
2571 : END DO
2572 : !$OMP END DO
2573 : !$OMP DO
2574 : DO ip = 0, np - 1
2575 : ny = boout(2, 2, ip) - boout(1, 2, ip) + 1
2576 : scount(ip) = nx*ny*nz
2577 : sdispl(ip) = nx*nz*my*ip
2578 : END DO
2579 : !$OMP END DO
2580 : !$OMP END PARALLEL
2581 104 : ny = boout(2, 2, mip) - boout(1, 2, mip) + 1
2582 : !$OMP PARALLEL DO DEFAULT(NONE), &
2583 : !$OMP PRIVATE(nx), &
2584 104 : !$OMP SHARED(np,boin,rcount,rdispl,nz,ny)
2585 : DO ip = 0, np - 1
2586 : nx = boin(2, 1, ip) - boin(1, 1, ip) + 1
2587 : rcount(ip) = nx*ny*nz
2588 : rdispl(ip) = ny*nz*(boin(1, 1, ip) - 1)
2589 : END DO
2590 : !$OMP END PARALLEL DO
2591 :
2592 104 : CALL group%alltoall(rbuf, scount, sdispl, sout, rcount, rdispl)
2593 :
2594 104 : CALL timestop(handle)
2595 :
2596 104 : END SUBROUTINE cube_transpose_6
2597 :
2598 : ! **************************************************************************************************
2599 : !> \brief ...
2600 : ! **************************************************************************************************
2601 25089 : SUBROUTINE init_fft_scratch_pool()
2602 :
2603 25089 : CALL release_fft_scratch_pool()
2604 :
2605 : ! Allocate first scratch and mark it as used
2606 25089 : ALLOCATE (fft_scratch_first)
2607 727581 : ALLOCATE (fft_scratch_first%fft_scratch)
2608 : ! this is a very special scratch, it seems, we always keep it 'most - recent' so we will never delete it
2609 25089 : fft_scratch_first%fft_scratch%last_tick = HUGE(fft_scratch_first%fft_scratch%last_tick)
2610 :
2611 25089 : init_fft_pool = init_fft_pool + 1
2612 :
2613 25089 : END SUBROUTINE init_fft_scratch_pool
2614 :
2615 : ! **************************************************************************************************
2616 : !> \brief ...
2617 : !> \param fft_scratch ...
2618 : ! **************************************************************************************************
2619 75697 : SUBROUTINE deallocate_fft_scratch_type(fft_scratch)
2620 : TYPE(fft_scratch_type), INTENT(INOUT) :: fft_scratch
2621 :
2622 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2623 : INTEGER :: ierr
2624 : COMPLEX(KIND=dp), POINTER :: dummy_ptr_z
2625 : #endif
2626 :
2627 : ! deallocate structures
2628 75697 : IF (ASSOCIATED(fft_scratch%ziptr)) THEN
2629 15871 : CALL fft_dealloc(fft_scratch%ziptr)
2630 : END IF
2631 75697 : IF (ASSOCIATED(fft_scratch%zoptr)) THEN
2632 15871 : CALL fft_dealloc(fft_scratch%zoptr)
2633 : END IF
2634 75697 : IF (ASSOCIATED(fft_scratch%p1buf)) THEN
2635 0 : CALL fft_dealloc(fft_scratch%p1buf)
2636 : END IF
2637 75697 : IF (ASSOCIATED(fft_scratch%p2buf)) THEN
2638 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2639 : dummy_ptr_z => fft_scratch%p2buf(1, 1)
2640 : ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
2641 : #else
2642 0 : CALL fft_dealloc(fft_scratch%p2buf)
2643 : #endif
2644 : END IF
2645 75697 : IF (ASSOCIATED(fft_scratch%p3buf)) THEN
2646 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2647 : dummy_ptr_z => fft_scratch%p3buf(1, 1)
2648 : ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
2649 : #else
2650 0 : CALL fft_dealloc(fft_scratch%p3buf)
2651 : #endif
2652 : END IF
2653 75697 : IF (ASSOCIATED(fft_scratch%p4buf)) THEN
2654 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2655 : dummy_ptr_z => fft_scratch%p4buf(1, 1)
2656 : ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
2657 : #else
2658 0 : CALL fft_dealloc(fft_scratch%p4buf)
2659 : #endif
2660 : END IF
2661 75697 : IF (ASSOCIATED(fft_scratch%p5buf)) THEN
2662 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2663 : dummy_ptr_z => fft_scratch%p5buf(1, 1)
2664 : ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
2665 : #else
2666 0 : CALL fft_dealloc(fft_scratch%p5buf)
2667 : #endif
2668 : END IF
2669 75697 : IF (ASSOCIATED(fft_scratch%p6buf)) THEN
2670 0 : CALL fft_dealloc(fft_scratch%p6buf)
2671 : END IF
2672 75697 : IF (ASSOCIATED(fft_scratch%p7buf)) THEN
2673 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2674 : dummy_ptr_z => fft_scratch%p7buf(1, 1)
2675 : ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
2676 : #else
2677 0 : CALL fft_dealloc(fft_scratch%p7buf)
2678 : #endif
2679 : END IF
2680 75697 : IF (ASSOCIATED(fft_scratch%r1buf)) THEN
2681 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2682 : dummy_ptr_z => fft_scratch%r1buf(1, 1)
2683 : ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
2684 : #else
2685 34729 : CALL fft_dealloc(fft_scratch%r1buf)
2686 : #endif
2687 : END IF
2688 75697 : IF (ASSOCIATED(fft_scratch%r2buf)) THEN
2689 34729 : CALL fft_dealloc(fft_scratch%r2buf)
2690 : END IF
2691 75697 : IF (ASSOCIATED(fft_scratch%tbuf)) THEN
2692 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2693 : dummy_ptr_z => fft_scratch%tbuf(1, 1, 1)
2694 : ierr = offload_free_pinned_mem(c_loc(dummy_ptr_z))
2695 : #else
2696 34729 : CALL fft_dealloc(fft_scratch%tbuf)
2697 : #endif
2698 : END IF
2699 75697 : IF (ASSOCIATED(fft_scratch%a1buf)) THEN
2700 8 : CALL fft_dealloc(fft_scratch%a1buf)
2701 : END IF
2702 75697 : IF (ASSOCIATED(fft_scratch%a2buf)) THEN
2703 8 : CALL fft_dealloc(fft_scratch%a2buf)
2704 : END IF
2705 75697 : IF (ASSOCIATED(fft_scratch%a3buf)) THEN
2706 8 : CALL fft_dealloc(fft_scratch%a3buf)
2707 : END IF
2708 75697 : IF (ASSOCIATED(fft_scratch%a4buf)) THEN
2709 8 : CALL fft_dealloc(fft_scratch%a4buf)
2710 : END IF
2711 75697 : IF (ASSOCIATED(fft_scratch%a5buf)) THEN
2712 8 : CALL fft_dealloc(fft_scratch%a5buf)
2713 : END IF
2714 75697 : IF (ASSOCIATED(fft_scratch%a6buf)) THEN
2715 8 : CALL fft_dealloc(fft_scratch%a6buf)
2716 : END IF
2717 75697 : IF (ASSOCIATED(fft_scratch%scount)) THEN
2718 0 : DEALLOCATE (fft_scratch%scount, fft_scratch%rcount, &
2719 34737 : fft_scratch%sdispl, fft_scratch%rdispl)
2720 : END IF
2721 75697 : IF (ASSOCIATED(fft_scratch%rr)) THEN
2722 34721 : DEALLOCATE (fft_scratch%rr)
2723 : END IF
2724 75697 : IF (ASSOCIATED(fft_scratch%xzbuf)) THEN
2725 0 : DEALLOCATE (fft_scratch%xzbuf)
2726 : END IF
2727 75697 : IF (ASSOCIATED(fft_scratch%yzbuf)) THEN
2728 0 : DEALLOCATE (fft_scratch%yzbuf)
2729 : END IF
2730 75697 : IF (ASSOCIATED(fft_scratch%xzbuf_sgl)) THEN
2731 0 : DEALLOCATE (fft_scratch%xzbuf_sgl)
2732 : END IF
2733 75697 : IF (ASSOCIATED(fft_scratch%yzbuf_sgl)) THEN
2734 0 : DEALLOCATE (fft_scratch%yzbuf_sgl)
2735 : END IF
2736 75697 : IF (ASSOCIATED(fft_scratch%ss)) THEN
2737 8 : DEALLOCATE (fft_scratch%ss)
2738 : END IF
2739 75697 : IF (ASSOCIATED(fft_scratch%tt)) THEN
2740 8 : DEALLOCATE (fft_scratch%tt)
2741 : END IF
2742 75697 : IF (ASSOCIATED(fft_scratch%pgrid)) THEN
2743 0 : DEALLOCATE (fft_scratch%pgrid)
2744 : END IF
2745 75697 : IF (ASSOCIATED(fft_scratch%pgcube)) THEN
2746 34737 : DEALLOCATE (fft_scratch%pgcube)
2747 : END IF
2748 75697 : IF (ASSOCIATED(fft_scratch%xcor)) THEN
2749 0 : DEALLOCATE (fft_scratch%xcor, fft_scratch%zcor)
2750 : END IF
2751 75697 : IF (ASSOCIATED(fft_scratch%pzcoord)) THEN
2752 0 : DEALLOCATE (fft_scratch%pzcoord)
2753 : END IF
2754 75697 : IF (ASSOCIATED(fft_scratch%xzcount)) THEN
2755 0 : DEALLOCATE (fft_scratch%xzcount, fft_scratch%yzcount)
2756 0 : DEALLOCATE (fft_scratch%xzdispl, fft_scratch%yzdispl)
2757 0 : fft_scratch%in = 0
2758 0 : fft_scratch%rsratio = 1._dp
2759 : END IF
2760 75697 : IF (ASSOCIATED(fft_scratch%rbuf1)) THEN
2761 0 : DEALLOCATE (fft_scratch%rbuf1)
2762 : END IF
2763 75697 : IF (ASSOCIATED(fft_scratch%rbuf2)) THEN
2764 0 : DEALLOCATE (fft_scratch%rbuf2)
2765 : END IF
2766 75697 : IF (ASSOCIATED(fft_scratch%rbuf3)) THEN
2767 0 : DEALLOCATE (fft_scratch%rbuf3)
2768 : END IF
2769 75697 : IF (ASSOCIATED(fft_scratch%rbuf4)) THEN
2770 0 : DEALLOCATE (fft_scratch%rbuf4)
2771 : END IF
2772 75697 : IF (ASSOCIATED(fft_scratch%rbuf5)) THEN
2773 8 : DEALLOCATE (fft_scratch%rbuf5)
2774 : END IF
2775 75697 : IF (ASSOCIATED(fft_scratch%rbuf6)) THEN
2776 8 : DEALLOCATE (fft_scratch%rbuf6)
2777 : END IF
2778 :
2779 75697 : IF (fft_scratch%cart_sub_comm(1) /= mp_comm_null) THEN
2780 0 : CALL fft_scratch%cart_sub_comm(1)%free()
2781 : END IF
2782 75697 : IF (fft_scratch%cart_sub_comm(2) /= mp_comm_null) THEN
2783 0 : CALL fft_scratch%cart_sub_comm(2)%free()
2784 : END IF
2785 :
2786 75697 : CALL fft_destroy_plan(fft_scratch%fft_plan(1))
2787 75697 : CALL fft_destroy_plan(fft_scratch%fft_plan(2))
2788 75697 : CALL fft_destroy_plan(fft_scratch%fft_plan(3))
2789 75697 : CALL fft_destroy_plan(fft_scratch%fft_plan(4))
2790 75697 : CALL fft_destroy_plan(fft_scratch%fft_plan(5))
2791 75697 : CALL fft_destroy_plan(fft_scratch%fft_plan(6))
2792 :
2793 75697 : END SUBROUTINE deallocate_fft_scratch_type
2794 :
2795 : ! **************************************************************************************************
2796 : !> \brief ...
2797 : ! **************************************************************************************************
2798 48085 : SUBROUTINE release_fft_scratch_pool()
2799 : TYPE(fft_scratch_pool_type), POINTER :: fft_scratch, fft_scratch_current
2800 :
2801 48085 : !$ CPASSERT(.NOT. omp_in_parallel() .OR. 0 == omp_get_thread_num())
2802 :
2803 48085 : IF (init_fft_pool == 0) NULLIFY (fft_scratch_first)
2804 :
2805 48085 : fft_scratch => fft_scratch_first
2806 75564 : DO
2807 123649 : IF (ASSOCIATED(fft_scratch)) THEN
2808 75564 : fft_scratch_current => fft_scratch
2809 75564 : fft_scratch => fft_scratch_current%fft_scratch_next
2810 75564 : NULLIFY (fft_scratch_current%fft_scratch_next)
2811 :
2812 75564 : CALL deallocate_fft_scratch_type(fft_scratch_current%fft_scratch)
2813 :
2814 302256 : DEALLOCATE (fft_scratch_current%fft_scratch)
2815 75564 : DEALLOCATE (fft_scratch_current)
2816 : ELSE
2817 : EXIT
2818 : END IF
2819 : END DO
2820 :
2821 48085 : init_fft_pool = 0
2822 :
2823 48085 : END SUBROUTINE release_fft_scratch_pool
2824 :
2825 : ! **************************************************************************************************
2826 : !> \brief ...
2827 : ! **************************************************************************************************
2828 4208195 : SUBROUTINE resize_fft_scratch_pool()
2829 :
2830 : INTEGER :: last_tick, nscratch
2831 : TYPE(fft_scratch_pool_type), POINTER :: fft_scratch_current, fft_scratch_old
2832 :
2833 4208195 : nscratch = 0
2834 :
2835 4208195 : last_tick = HUGE(last_tick)
2836 4208195 : NULLIFY (fft_scratch_old)
2837 :
2838 : ! start at the global pool, count, and find a deletion candidate
2839 4208195 : fft_scratch_current => fft_scratch_first
2840 22557384 : DO
2841 26765579 : IF (ASSOCIATED(fft_scratch_current)) THEN
2842 22557384 : nscratch = nscratch + 1
2843 : ! is this a candidate for deletion (i.e. least recently used, and not in use)
2844 22557384 : IF (.NOT. fft_scratch_current%fft_scratch%in_use) THEN
2845 18349189 : IF (fft_scratch_current%fft_scratch%last_tick < last_tick) THEN
2846 8688439 : last_tick = fft_scratch_current%fft_scratch%last_tick
2847 8688439 : fft_scratch_old => fft_scratch_current
2848 : END IF
2849 : END IF
2850 22557384 : fft_scratch_current => fft_scratch_current%fft_scratch_next
2851 : ELSE
2852 : EXIT
2853 : END IF
2854 : END DO
2855 :
2856 : ! we should delete a scratch
2857 4208195 : IF (nscratch > fft_pool_scratch_limit) THEN
2858 : ! note that we never deallocate the first (special) element of the list
2859 133 : IF (ASSOCIATED(fft_scratch_old)) THEN
2860 : fft_scratch_current => fft_scratch_first
2861 : DO
2862 2261 : IF (ASSOCIATED(fft_scratch_current)) THEN
2863 : ! should we delete the next in the list?
2864 2128 : IF (ASSOCIATED(fft_scratch_current%fft_scratch_next, fft_scratch_old)) THEN
2865 : ! fix the linked list
2866 133 : fft_scratch_current%fft_scratch_next => fft_scratch_old%fft_scratch_next
2867 :
2868 : ! deallocate the element
2869 133 : CALL deallocate_fft_scratch_type(fft_scratch_old%fft_scratch)
2870 532 : DEALLOCATE (fft_scratch_old%fft_scratch)
2871 133 : DEALLOCATE (fft_scratch_old)
2872 :
2873 : ELSE
2874 : fft_scratch_current => fft_scratch_current%fft_scratch_next
2875 : END IF
2876 : ELSE
2877 : EXIT
2878 : END IF
2879 : END DO
2880 :
2881 : ELSE
2882 0 : CPWARN("The number of the scratches exceeded the limit, but none could be deallocated")
2883 : END IF
2884 : END IF
2885 :
2886 4208195 : END SUBROUTINE resize_fft_scratch_pool
2887 :
2888 : ! **************************************************************************************************
2889 : !> \brief ...
2890 : !> \param fft_scratch ...
2891 : !> \param tf_type ...
2892 : !> \param n ...
2893 : !> \param fft_sizes ...
2894 : ! **************************************************************************************************
2895 4208195 : SUBROUTINE get_fft_scratch(fft_scratch, tf_type, n, fft_sizes)
2896 : TYPE(fft_scratch_type), POINTER :: fft_scratch
2897 : INTEGER, INTENT(IN) :: tf_type
2898 : INTEGER, DIMENSION(:), INTENT(IN) :: n
2899 : TYPE(fft_scratch_sizes), INTENT(IN), &
2900 : OPTIONAL :: fft_sizes
2901 :
2902 : CHARACTER(len=*), PARAMETER :: routineN = 'get_fft_scratch'
2903 :
2904 : INTEGER :: coord(2), DIM(2), handle, i, ix, iz, lg, lmax, m1, m2, &
2905 : mcx2, mcy3, mcz1, mcz2, mg, mmax, mx1, mx2, my1, my3, mz1, mz2, mz3, &
2906 : nbx, nbz, nm, nmax, nmray, np, nx, ny, nyzray, nz, pos(2)
2907 : INTEGER, DIMENSION(3) :: pcoord
2908 : LOGICAL :: equal
2909 : LOGICAL, DIMENSION(2) :: dims
2910 : TYPE(fft_scratch_pool_type), POINTER :: fft_scratch_current, &
2911 : fft_scratch_last, &
2912 : fft_scratch_new
2913 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
2914 : INTEGER :: ierr
2915 : INTEGER(KIND=C_SIZE_T) :: length
2916 : TYPE(C_PTR) :: cptr_r1buf, cptr_tbuf, &
2917 : cptr_p2buf, cptr_p3buf, cptr_p4buf, cptr_p5buf, cptr_p7buf
2918 : #endif
2919 4208195 : CALL timeset(routineN, handle)
2920 :
2921 : ! this is the place to check that the scratch_pool does not grow without limits
2922 : ! before we add a new scratch check the size of the pool and release some of the list if needed
2923 4208195 : CALL resize_fft_scratch_pool()
2924 :
2925 : ! get the required scratch
2926 4208195 : !$OMP ATOMIC
2927 : tick_fft_pool = tick_fft_pool + 1
2928 4208195 : fft_scratch_current => fft_scratch_first
2929 : DO
2930 12988803 : IF (ASSOCIATED(fft_scratch_current)) THEN
2931 12938195 : IF (fft_scratch_current%fft_scratch%in_use) THEN
2932 4208195 : fft_scratch_last => fft_scratch_current
2933 4208195 : fft_scratch_current => fft_scratch_current%fft_scratch_next
2934 4208195 : CYCLE
2935 : END IF
2936 8730000 : IF (tf_type /= fft_scratch_current%fft_scratch%tf_type) THEN
2937 458894 : fft_scratch_last => fft_scratch_current
2938 458894 : fft_scratch_current => fft_scratch_current%fft_scratch_next
2939 458894 : CYCLE
2940 : END IF
2941 21015458 : IF (.NOT. ALL(n == fft_scratch_current%fft_scratch%nfft)) THEN
2942 4023714 : fft_scratch_last => fft_scratch_current
2943 4023714 : fft_scratch_current => fft_scratch_current%fft_scratch_next
2944 4023714 : CYCLE
2945 : END IF
2946 4247392 : IF (PRESENT(fft_sizes)) THEN
2947 3585118 : IF (fft_sizes%rs_group /= fft_scratch_current%fft_scratch%group) THEN
2948 89485 : fft_scratch_last => fft_scratch_current
2949 89485 : fft_scratch_current => fft_scratch_current%fft_scratch_next
2950 89485 : CYCLE
2951 : END IF
2952 3495633 : CALL is_equal(fft_sizes, fft_scratch_current%fft_scratch%sizes, equal)
2953 3495633 : IF (.NOT. equal) THEN
2954 320 : fft_scratch_last => fft_scratch_current
2955 320 : fft_scratch_current => fft_scratch_current%fft_scratch_next
2956 320 : CYCLE
2957 : END IF
2958 : END IF
2959 : ! Success
2960 4157587 : fft_scratch => fft_scratch_current%fft_scratch
2961 4157587 : fft_scratch_current%fft_scratch%in_use = .TRUE.
2962 4157587 : EXIT
2963 : ELSE
2964 : ! We cannot find the scratch type in this pool
2965 : ! Generate a new scratch set
2966 50608 : ALLOCATE (fft_scratch_new)
2967 1315808 : ALLOCATE (fft_scratch_new%fft_scratch)
2968 :
2969 50608 : IF (tf_type /= 400) THEN
2970 34737 : fft_scratch_new%fft_scratch%sizes = fft_sizes
2971 34737 : np = fft_sizes%numtask
2972 : ALLOCATE (fft_scratch_new%fft_scratch%scount(0:np - 1), fft_scratch_new%fft_scratch%rcount(0:np - 1), &
2973 : fft_scratch_new%fft_scratch%sdispl(0:np - 1), fft_scratch_new%fft_scratch%rdispl(0:np - 1), &
2974 312633 : fft_scratch_new%fft_scratch%pgcube(0:np - 1, 2))
2975 : END IF
2976 :
2977 0 : SELECT CASE (tf_type)
2978 : CASE DEFAULT
2979 0 : CPABORT("Invalid scratch type.")
2980 : CASE (100) ! fft3d_pb: full cube distribution
2981 0 : CPASSERT(PRESENT(fft_sizes))
2982 0 : mx1 = fft_sizes%mx1
2983 0 : my1 = fft_sizes%my1
2984 0 : mx2 = fft_sizes%mx2
2985 0 : mz2 = fft_sizes%mz2
2986 0 : my3 = fft_sizes%my3
2987 0 : mz3 = fft_sizes%mz3
2988 0 : CALL fft_alloc(fft_scratch_new%fft_scratch%a1buf, [mx1*my1, n(3)])
2989 0 : CALL fft_alloc(fft_scratch_new%fft_scratch%a2buf, [n(3), mx1*my1])
2990 0 : CALL fft_alloc(fft_scratch_new%fft_scratch%a3buf, [mx2*mz2, n(2)])
2991 0 : CALL fft_alloc(fft_scratch_new%fft_scratch%a4buf, [n(2), mx2*mz2])
2992 0 : CALL fft_alloc(fft_scratch_new%fft_scratch%a5buf, [my3*mz3, n(1)])
2993 0 : CALL fft_alloc(fft_scratch_new%fft_scratch%a6buf, [n(1), my3*mz3])
2994 0 : fft_scratch_new%fft_scratch%group = fft_sizes%rs_group
2995 :
2996 0 : dim = fft_sizes%rs_group%num_pe_cart
2997 0 : pos = fft_sizes%rs_group%mepos_cart
2998 0 : fft_scratch_new%fft_scratch%mip = fft_sizes%rs_group%mepos
2999 0 : fft_scratch_new%fft_scratch%dim = dim
3000 0 : fft_scratch_new%fft_scratch%pos = pos
3001 0 : mcz1 = fft_sizes%mcz1
3002 0 : mcx2 = fft_sizes%mcx2
3003 0 : mcz2 = fft_sizes%mcz2
3004 0 : mcy3 = fft_sizes%mcy3
3005 0 : ALLOCATE (fft_scratch_new%fft_scratch%rbuf1(mx2*my1*mcz2, 0:DIM(2) - 1))
3006 0 : ALLOCATE (fft_scratch_new%fft_scratch%rbuf2(mx1*my1*mcz2, 0:DIM(2) - 1))
3007 0 : ALLOCATE (fft_scratch_new%fft_scratch%rbuf3(mx2*mz3*mcy3, 0:DIM(1) - 1))
3008 0 : ALLOCATE (fft_scratch_new%fft_scratch%rbuf4(mx2*mz2*mcy3, 0:DIM(1) - 1))
3009 :
3010 0 : dims = [.TRUE., .FALSE.]
3011 0 : CALL fft_scratch_new%fft_scratch%cart_sub_comm(1)%from_sub(fft_sizes%rs_group, dims)
3012 0 : dims = [.FALSE., .TRUE.]
3013 0 : CALL fft_scratch_new%fft_scratch%cart_sub_comm(2)%from_sub(fft_sizes%rs_group, dims)
3014 :
3015 : !initialise pgcube
3016 0 : DO i = 0, DIM(1) - 1
3017 0 : coord = [i, pos(2)]
3018 0 : CALL fft_sizes%rs_group%rank_cart(coord, fft_scratch_new%fft_scratch%pgcube(i, 1))
3019 : END DO
3020 0 : DO i = 0, DIM(2) - 1
3021 0 : coord = [pos(1), i]
3022 0 : CALL fft_sizes%rs_group%rank_cart(coord, fft_scratch_new%fft_scratch%pgcube(i, 2))
3023 : END DO
3024 :
3025 : !set up fft plans
3026 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(1), FWFFT, .TRUE., .FALSE., &
3027 : mx1*my1, n(3), n(3), mx1*my1, &
3028 0 : fft_scratch_new%fft_scratch%a1buf, fft_scratch_new%fft_scratch%a2buf)
3029 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(2), FWFFT, .TRUE., .FALSE., &
3030 : mx2*mz2, n(2), n(2), mx2*mz2, &
3031 0 : fft_scratch_new%fft_scratch%a3buf, fft_scratch_new%fft_scratch%a4buf)
3032 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(3), FWFFT, .TRUE., .FALSE., &
3033 : my3*mz3, n(1), n(1), my3*mz3, &
3034 0 : fft_scratch_new%fft_scratch%a5buf, fft_scratch_new%fft_scratch%a6buf)
3035 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(4), BWFFT, .FALSE., .TRUE., &
3036 : n(1), my3*mz3, n(1), my3*mz3, &
3037 0 : fft_scratch_new%fft_scratch%a6buf, fft_scratch_new%fft_scratch%a5buf)
3038 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(5), BWFFT, .FALSE., .TRUE., &
3039 : n(2), mx2*mz2, n(2), mx2*mz2, &
3040 0 : fft_scratch_new%fft_scratch%a4buf, fft_scratch_new%fft_scratch%a3buf)
3041 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(6), BWFFT, .FALSE., .TRUE., &
3042 : n(3), mx1*my1, n(3), mx1*my1, &
3043 0 : fft_scratch_new%fft_scratch%a2buf, fft_scratch_new%fft_scratch%a1buf)
3044 :
3045 : CASE (101) ! fft3d_pb: full cube distribution (dim 1)
3046 8 : CPASSERT(PRESENT(fft_sizes))
3047 8 : mx1 = fft_sizes%mx1
3048 8 : my1 = fft_sizes%my1
3049 8 : mz1 = fft_sizes%mz1
3050 8 : my3 = fft_sizes%my3
3051 8 : mz3 = fft_sizes%mz3
3052 24 : CALL fft_alloc(fft_scratch_new%fft_scratch%a1buf, [mx1*my1, n(3)])
3053 24 : CALL fft_alloc(fft_scratch_new%fft_scratch%a2buf, [n(3), mx1*my1])
3054 8 : fft_scratch_new%fft_scratch%group = fft_sizes%rs_group
3055 24 : CALL fft_alloc(fft_scratch_new%fft_scratch%a3buf, [mx1*mz1, n(2)])
3056 24 : CALL fft_alloc(fft_scratch_new%fft_scratch%a4buf, [n(2), mx1*mz1])
3057 24 : CALL fft_alloc(fft_scratch_new%fft_scratch%a5buf, [my3*mz3, n(1)])
3058 24 : CALL fft_alloc(fft_scratch_new%fft_scratch%a6buf, [n(1), my3*mz3])
3059 :
3060 24 : dim = fft_sizes%rs_group%num_pe_cart
3061 24 : pos = fft_sizes%rs_group%mepos_cart
3062 8 : fft_scratch_new%fft_scratch%mip = fft_sizes%rs_group%mepos
3063 24 : fft_scratch_new%fft_scratch%dim = dim
3064 24 : fft_scratch_new%fft_scratch%pos = pos
3065 8 : mcy3 = fft_sizes%mcy3
3066 32 : ALLOCATE (fft_scratch_new%fft_scratch%rbuf5(mx1*mz3*mcy3, 0:DIM(1) - 1))
3067 32 : ALLOCATE (fft_scratch_new%fft_scratch%rbuf6(mx1*mz1*mcy3, 0:DIM(1) - 1))
3068 :
3069 : !set up fft plans
3070 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(1), FWFFT, .TRUE., .FALSE., &
3071 : mx1*my1, n(3), n(3), mx1*my1, &
3072 8 : fft_scratch_new%fft_scratch%a1buf, fft_scratch_new%fft_scratch%a3buf)
3073 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(2), FWFFT, .TRUE., .FALSE., &
3074 : mx1*mz1, n(2), n(2), mx1*mz1, &
3075 8 : fft_scratch_new%fft_scratch%a3buf, fft_scratch_new%fft_scratch%a4buf)
3076 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(3), FWFFT, .TRUE., .FALSE., &
3077 : my3*mz3, n(1), n(1), my3*mz3, &
3078 8 : fft_scratch_new%fft_scratch%a5buf, fft_scratch_new%fft_scratch%a6buf)
3079 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(4), BWFFT, .FALSE., .TRUE., &
3080 : n(1), my3*mz3, n(1), my3*mz3, &
3081 8 : fft_scratch_new%fft_scratch%a6buf, fft_scratch_new%fft_scratch%a5buf)
3082 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(5), BWFFT, .FALSE., .TRUE., &
3083 : n(2), mx1*mz1, n(2), mx1*mz1, &
3084 8 : fft_scratch_new%fft_scratch%a4buf, fft_scratch_new%fft_scratch%a3buf)
3085 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(6), BWFFT, .FALSE., .TRUE., &
3086 : n(3), mx1*my1, n(3), mx1*my1, &
3087 8 : fft_scratch_new%fft_scratch%a3buf, fft_scratch_new%fft_scratch%a1buf)
3088 :
3089 : CASE (200) ! fft3d_ps: plane distribution
3090 34729 : CPASSERT(PRESENT(fft_sizes))
3091 34729 : nx = fft_sizes%nx
3092 34729 : ny = fft_sizes%ny
3093 34729 : nz = fft_sizes%nz
3094 34729 : mx2 = fft_sizes%mx2
3095 34729 : lmax = fft_sizes%lmax
3096 34729 : mmax = fft_sizes%mmax
3097 34729 : lg = fft_sizes%lg
3098 34729 : mg = fft_sizes%mg
3099 34729 : np = fft_sizes%numtask
3100 34729 : nmray = fft_sizes%nmray
3101 34729 : nyzray = fft_sizes%nyzray
3102 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
3103 : length = INT(2*dp_size*MAX(mmax, 1)*MAX(lmax, 1), KIND=C_SIZE_T)
3104 : ierr = offload_malloc_pinned_mem(cptr_r1buf, length)
3105 : CPASSERT(ierr == 0)
3106 : CALL c_f_pointer(cptr_r1buf, fft_scratch_new%fft_scratch%r1buf, [MAX(mmax, 1), MAX(lmax, 1)])
3107 : length = INT(2*dp_size*MAX(ny, 1)*MAX(nz, 1)*MAX(nx, 1), KIND=C_SIZE_T)
3108 : ierr = offload_malloc_pinned_mem(cptr_tbuf, length)
3109 : CPASSERT(ierr == 0)
3110 : CALL c_f_pointer(cptr_tbuf, fft_scratch_new%fft_scratch%tbuf, [MAX(ny, 1), MAX(nz, 1), MAX(nx, 1)])
3111 : #else
3112 104187 : CALL fft_alloc(fft_scratch_new%fft_scratch%r1buf, [mmax, lmax])
3113 138916 : CALL fft_alloc(fft_scratch_new%fft_scratch%tbuf, [ny, nz, nx])
3114 : #endif
3115 34729 : fft_scratch_new%fft_scratch%group = fft_sizes%rs_group
3116 104187 : CALL fft_alloc(fft_scratch_new%fft_scratch%r2buf, [lg, mg])
3117 34729 : nm = nmray*mx2
3118 34729 : IF (alltoall_sgl) THEN
3119 32 : ALLOCATE (fft_scratch_new%fft_scratch%ss(mmax, lmax))
3120 32 : ALLOCATE (fft_scratch_new%fft_scratch%tt(nm, 0:np - 1))
3121 : ELSE
3122 138884 : ALLOCATE (fft_scratch_new%fft_scratch%rr(nm, 0:np - 1))
3123 : END IF
3124 :
3125 : !set up fft plans
3126 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(1), FWFFT, .TRUE., .FALSE., &
3127 : nx*ny, nz, nz, nx*ny, &
3128 34729 : fft_scratch_new%fft_scratch%tbuf, fft_scratch_new%fft_scratch%r1buf)
3129 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(2), FWFFT, .TRUE., .FALSE., &
3130 : nx*nz, ny, ny, nx*nz, &
3131 34729 : fft_scratch_new%fft_scratch%r1buf, fft_scratch_new%fft_scratch%tbuf)
3132 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(3), FWFFT, .TRUE., .FALSE., &
3133 : nyzray, n(1), n(1), nyzray, &
3134 34729 : fft_scratch_new%fft_scratch%r1buf, fft_scratch_new%fft_scratch%r2buf)
3135 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(4), BWFFT, .FALSE., .TRUE., &
3136 : n(1), nyzray, n(1), nyzray, &
3137 34729 : fft_scratch_new%fft_scratch%r2buf, fft_scratch_new%fft_scratch%r1buf)
3138 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(5), BWFFT, .FALSE., .TRUE., &
3139 : ny, nx*nz, ny, nx*nz, &
3140 34729 : fft_scratch_new%fft_scratch%tbuf, fft_scratch_new%fft_scratch%r1buf)
3141 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(6), BWFFT, .FALSE., .TRUE., &
3142 : nz, nx*ny, nz, nx*ny, &
3143 34729 : fft_scratch_new%fft_scratch%r1buf, fft_scratch_new%fft_scratch%tbuf)
3144 :
3145 : CASE (300) ! fft3d_ps: block distribution
3146 0 : CPASSERT(PRESENT(fft_sizes))
3147 0 : mx1 = fft_sizes%mx1
3148 0 : mx2 = fft_sizes%mx2
3149 0 : my1 = fft_sizes%my1
3150 0 : mz2 = fft_sizes%mz2
3151 0 : mcx2 = fft_sizes%mcx2
3152 0 : lg = fft_sizes%lg
3153 0 : mg = fft_sizes%mg
3154 0 : nmax = fft_sizes%nmax
3155 0 : nmray = fft_sizes%nmray
3156 0 : nyzray = fft_sizes%nyzray
3157 0 : m1 = fft_sizes%r_dim(1)
3158 0 : m2 = fft_sizes%r_dim(2)
3159 0 : nbx = fft_sizes%nbx
3160 0 : nbz = fft_sizes%nbz
3161 0 : CALL fft_alloc(fft_scratch_new%fft_scratch%p1buf, [mx1*my1, n(3)])
3162 0 : CALL fft_alloc(fft_scratch_new%fft_scratch%p6buf, [lg, mg])
3163 : #if defined(__OFFLOAD) && !defined(__NO_OFFLOAD_PW)
3164 : length = INT(2*dp_size*MAX(n(3), 1)*MAX(mx1*my1, 1), KIND=C_SIZE_T)
3165 : ierr = offload_malloc_pinned_mem(cptr_p2buf, length)
3166 : CPASSERT(ierr == 0)
3167 : CALL c_f_pointer(cptr_p2buf, fft_scratch_new%fft_scratch%p2buf, [MAX(n(3), 1), MAX(mx1*my1, 1)])
3168 : length = INT(2*dp_size*MAX(mx2*mz2, 1)*MAX(n(2), 1), KIND=C_SIZE_T)
3169 : ierr = offload_malloc_pinned_mem(cptr_p3buf, length)
3170 : CPASSERT(ierr == 0)
3171 : CALL c_f_pointer(cptr_p3buf, fft_scratch_new%fft_scratch%p3buf, [MAX(mx2*mz2, 1), MAX(n(2), 1)])
3172 : length = INT(2*dp_size*MAX(n(2), 1)*MAX(mx2*mz2, 1), KIND=C_SIZE_T)
3173 : ierr = offload_malloc_pinned_mem(cptr_p4buf, length)
3174 : CPASSERT(ierr == 0)
3175 : CALL c_f_pointer(cptr_p4buf, fft_scratch_new%fft_scratch%p4buf, [MAX(n(2), 1), MAX(mx2*mz2, 1)])
3176 : length = INT(2*dp_size*MAX(nyzray, 1)*MAX(n(1), 1), KIND=C_SIZE_T)
3177 : ierr = offload_malloc_pinned_mem(cptr_p5buf, length)
3178 : CPASSERT(ierr == 0)
3179 : CALL c_f_pointer(cptr_p5buf, fft_scratch_new%fft_scratch%p5buf, [MAX(nyzray, 1), MAX(n(1), 1)])
3180 : length = INT(2*dp_size*MAX(mg, 1)*MAX(lg, 1), KIND=C_SIZE_T)
3181 : ierr = offload_malloc_pinned_mem(cptr_p7buf, length)
3182 : CPASSERT(ierr == 0)
3183 : CALL c_f_pointer(cptr_p7buf, fft_scratch_new%fft_scratch%p7buf, [MAX(mg, 1), MAX(lg, 1)])
3184 : #else
3185 0 : CALL fft_alloc(fft_scratch_new%fft_scratch%p2buf, [n(3), mx1*my1])
3186 0 : CALL fft_alloc(fft_scratch_new%fft_scratch%p3buf, [mx2*mz2, n(2)])
3187 0 : CALL fft_alloc(fft_scratch_new%fft_scratch%p4buf, [n(2), mx2*mz2])
3188 0 : CALL fft_alloc(fft_scratch_new%fft_scratch%p5buf, [nyzray, n(1)])
3189 0 : CALL fft_alloc(fft_scratch_new%fft_scratch%p7buf, [mg, lg])
3190 : #endif
3191 0 : IF (alltoall_sgl) THEN
3192 0 : ALLOCATE (fft_scratch_new%fft_scratch%yzbuf_sgl(mg*lg))
3193 0 : ALLOCATE (fft_scratch_new%fft_scratch%xzbuf_sgl(n(2)*mx2*mz2))
3194 : ELSE
3195 0 : ALLOCATE (fft_scratch_new%fft_scratch%yzbuf(mg*lg))
3196 0 : ALLOCATE (fft_scratch_new%fft_scratch%xzbuf(n(2)*mx2*mz2))
3197 : END IF
3198 0 : ALLOCATE (fft_scratch_new%fft_scratch%pgrid(0:m1 - 1, 0:m2 - 1))
3199 0 : ALLOCATE (fft_scratch_new%fft_scratch%xcor(nbx))
3200 0 : ALLOCATE (fft_scratch_new%fft_scratch%zcor(nbz))
3201 0 : ALLOCATE (fft_scratch_new%fft_scratch%pzcoord(0:np - 1))
3202 : ALLOCATE (fft_scratch_new%fft_scratch%xzcount(0:np - 1), &
3203 0 : fft_scratch_new%fft_scratch%yzcount(0:np - 1))
3204 : ALLOCATE (fft_scratch_new%fft_scratch%xzdispl(0:np - 1), &
3205 0 : fft_scratch_new%fft_scratch%yzdispl(0:np - 1))
3206 0 : fft_scratch_new%fft_scratch%group = fft_sizes%rs_group
3207 :
3208 0 : dim = fft_sizes%rs_group%num_pe_cart
3209 0 : pos = fft_sizes%rs_group%mepos_cart
3210 0 : fft_scratch_new%fft_scratch%mip = fft_sizes%rs_group%mepos
3211 0 : fft_scratch_new%fft_scratch%dim = dim
3212 0 : fft_scratch_new%fft_scratch%pos = pos
3213 0 : mcz1 = fft_sizes%mcz1
3214 0 : mcz2 = fft_sizes%mcz2
3215 0 : ALLOCATE (fft_scratch_new%fft_scratch%rbuf1(mx2*my1*mcz2, 0:DIM(2) - 1))
3216 0 : ALLOCATE (fft_scratch_new%fft_scratch%rbuf2(mx1*my1*mcz2, 0:DIM(2) - 1))
3217 :
3218 0 : dims = [.FALSE., .TRUE.]
3219 0 : CALL fft_scratch_new%fft_scratch%cart_sub_comm(2)%from_sub(fft_sizes%rs_group, dims)
3220 :
3221 : !initialise pgcube
3222 0 : DO i = 0, DIM(2) - 1
3223 0 : coord = [pos(1), i]
3224 0 : CALL fft_sizes%rs_group%rank_cart(coord, fft_scratch_new%fft_scratch%pgcube(i, 2))
3225 : END DO
3226 :
3227 : !initialise pgrid
3228 0 : DO ix = 0, m1 - 1
3229 0 : DO iz = 0, m2 - 1
3230 0 : coord = [ix, iz]
3231 0 : CALL fft_sizes%rs_group%rank_cart(coord, fft_scratch_new%fft_scratch%pgrid(ix, iz))
3232 : END DO
3233 : END DO
3234 :
3235 : !initialise pzcoord
3236 0 : DO i = 0, np - 1
3237 0 : CALL fft_sizes%rs_group%coords(i, pcoord)
3238 0 : fft_scratch_new%fft_scratch%pzcoord(i) = pcoord(2)
3239 : END DO
3240 :
3241 : !set up fft plans
3242 : ! The scratch buffers must cover the geometry declared in each call
3243 : ! below, because the FFTW planner may write the buffers during
3244 : ! planning. The plans pair the buffers such that each of p1-p6 is
3245 : ! checked once; p7 shares the size of p6.
3246 0 : CPASSERT(SIZE(fft_scratch_new%fft_scratch%p1buf) >= mx1*my1*n(3))
3247 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(1), FWFFT, .TRUE., .FALSE., &
3248 : mx1*my1, n(3), n(3), mx1*my1, &
3249 0 : fft_scratch_new%fft_scratch%p1buf, fft_scratch_new%fft_scratch%p2buf)
3250 0 : CPASSERT(SIZE(fft_scratch_new%fft_scratch%p3buf) >= mx2*mz2*n(2))
3251 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(2), FWFFT, .TRUE., .FALSE., &
3252 : mx2*mz2, n(2), n(2), mx2*mz2, &
3253 0 : fft_scratch_new%fft_scratch%p3buf, fft_scratch_new%fft_scratch%p4buf)
3254 0 : CPASSERT(SIZE(fft_scratch_new%fft_scratch%p5buf) >= nyzray*n(1))
3255 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(3), FWFFT, .TRUE., .FALSE., &
3256 : nyzray, n(1), n(1), nyzray, &
3257 0 : fft_scratch_new%fft_scratch%p5buf, fft_scratch_new%fft_scratch%p6buf)
3258 0 : CPASSERT(SIZE(fft_scratch_new%fft_scratch%p6buf) >= n(1)*nyzray)
3259 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(4), BWFFT, .FALSE., .TRUE., &
3260 : n(1), nyzray, n(1), nyzray, &
3261 0 : fft_scratch_new%fft_scratch%p6buf, fft_scratch_new%fft_scratch%p7buf)
3262 0 : CPASSERT(SIZE(fft_scratch_new%fft_scratch%p4buf) >= n(2)*mx2*mz2)
3263 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(5), BWFFT, .FALSE., .TRUE., &
3264 : n(2), mx2*mz2, n(2), mx2*mz2, &
3265 0 : fft_scratch_new%fft_scratch%p4buf, fft_scratch_new%fft_scratch%p3buf)
3266 0 : CPASSERT(SIZE(fft_scratch_new%fft_scratch%p2buf) >= n(3)*mx1*my1)
3267 : CALL fft_create_plan_1d(fft_scratch_new%fft_scratch%fft_plan(6), BWFFT, .FALSE., .TRUE., &
3268 : n(3), mx1*my1, n(3), mx1*my1, &
3269 0 : fft_scratch_new%fft_scratch%p2buf, fft_scratch_new%fft_scratch%p1buf)
3270 :
3271 : CASE (400) ! serial FFT
3272 15871 : np = 0
3273 15871 : CALL fft_alloc(fft_scratch_new%fft_scratch%ziptr, n)
3274 15871 : CALL fft_alloc(fft_scratch_new%fft_scratch%zoptr, n)
3275 :
3276 : !in place plans
3277 : CALL fft_create_plan_3d(fft_scratch_new%fft_scratch%fft_plan(1), .TRUE., FWFFT, n, &
3278 15871 : fft_scratch_new%fft_scratch%ziptr, fft_scratch_new%fft_scratch%ziptr)
3279 : CALL fft_create_plan_3d(fft_scratch_new%fft_scratch%fft_plan(2), .TRUE., BWFFT, n, &
3280 15871 : fft_scratch_new%fft_scratch%ziptr, fft_scratch_new%fft_scratch%ziptr)
3281 : ! out of place plans
3282 : CALL fft_create_plan_3d(fft_scratch_new%fft_scratch%fft_plan(3), .FALSE., FWFFT, n, &
3283 15871 : fft_scratch_new%fft_scratch%ziptr, fft_scratch_new%fft_scratch%zoptr)
3284 : CALL fft_create_plan_3d(fft_scratch_new%fft_scratch%fft_plan(4), .FALSE., BWFFT, n, &
3285 66479 : fft_scratch_new%fft_scratch%ziptr, fft_scratch_new%fft_scratch%zoptr)
3286 :
3287 : END SELECT
3288 :
3289 50608 : NULLIFY (fft_scratch_new%fft_scratch_next)
3290 : fft_scratch_new%fft_scratch%fft_scratch_id = &
3291 50608 : fft_scratch_last%fft_scratch%fft_scratch_id + 1
3292 50608 : fft_scratch_new%fft_scratch%in_use = .TRUE.
3293 202432 : fft_scratch_new%fft_scratch%nfft = n
3294 50608 : fft_scratch_last%fft_scratch_next => fft_scratch_new
3295 50608 : fft_scratch_new%fft_scratch%tf_type = tf_type
3296 50608 : fft_scratch => fft_scratch_new%fft_scratch
3297 50608 : EXIT
3298 :
3299 : END IF
3300 : END DO
3301 :
3302 4208195 : !$OMP ATOMIC READ
3303 : fft_scratch%last_tick = tick_fft_pool
3304 :
3305 4208195 : CALL timestop(handle)
3306 :
3307 4208195 : END SUBROUTINE get_fft_scratch
3308 :
3309 : ! **************************************************************************************************
3310 : !> \brief ...
3311 : !> \param fft_scratch ...
3312 : ! **************************************************************************************************
3313 4208195 : SUBROUTINE release_fft_scratch(fft_scratch)
3314 :
3315 : TYPE(fft_scratch_type), POINTER :: fft_scratch
3316 :
3317 : INTEGER :: scratch_id
3318 : TYPE(fft_scratch_pool_type), POINTER :: fft_scratch_current
3319 :
3320 4208195 : scratch_id = fft_scratch%fft_scratch_id
3321 :
3322 4208195 : fft_scratch_current => fft_scratch_first
3323 8780608 : DO
3324 12988803 : IF (ASSOCIATED(fft_scratch_current)) THEN
3325 12988803 : IF (scratch_id == fft_scratch_current%fft_scratch%fft_scratch_id) THEN
3326 4208195 : fft_scratch%in_use = .FALSE.
3327 4208195 : NULLIFY (fft_scratch)
3328 4208195 : EXIT
3329 : END IF
3330 8780608 : fft_scratch_current => fft_scratch_current%fft_scratch_next
3331 : ELSE
3332 : ! We cannot find the scratch type in this pool
3333 0 : CPABORT("Invalid scratch type.")
3334 0 : EXIT
3335 : END IF
3336 : END DO
3337 :
3338 4208195 : END SUBROUTINE release_fft_scratch
3339 :
3340 : ! **************************************************************************************************
3341 : !> \brief ...
3342 : !> \param rs ...
3343 : !> \param scount ...
3344 : !> \param sdispl ...
3345 : !> \param rq ...
3346 : !> \param rcount ...
3347 : !> \param rdispl ...
3348 : !> \param group ...
3349 : ! **************************************************************************************************
3350 0 : SUBROUTINE sparse_alltoall(rs, scount, sdispl, rq, rcount, rdispl, group)
3351 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: rs
3352 : INTEGER, DIMENSION(:), POINTER :: scount, sdispl
3353 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: rq
3354 : INTEGER, DIMENSION(:), POINTER :: rcount, rdispl
3355 :
3356 : CLASS(mp_comm_type), INTENT(IN) :: group
3357 :
3358 0 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: msgin, msgout
3359 : INTEGER :: ip, n, nr, ns, pos
3360 0 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: rreq, sreq
3361 :
3362 0 : CALL group%sync()
3363 0 : n = group%num_pe
3364 0 : pos = group%mepos
3365 0 : ALLOCATE (sreq(0:n - 1))
3366 0 : ALLOCATE (rreq(0:n - 1))
3367 0 : nr = 0
3368 0 : DO ip = 0, n - 1
3369 0 : IF (rcount(ip) == 0) CYCLE
3370 0 : IF (ip == pos) CYCLE
3371 0 : msgout => rq(rdispl(ip) + 1:rdispl(ip) + rcount(ip))
3372 0 : CALL group%irecv(msgout, ip, rreq(nr))
3373 0 : nr = nr + 1
3374 : END DO
3375 0 : ns = 0
3376 0 : DO ip = 0, n - 1
3377 0 : IF (scount(ip) == 0) CYCLE
3378 0 : IF (ip == pos) CYCLE
3379 0 : msgin => rs(sdispl(ip) + 1:sdispl(ip) + scount(ip))
3380 0 : CALL group%isend(msgin, ip, sreq(ns))
3381 0 : ns = ns + 1
3382 : END DO
3383 0 : IF (rcount(pos) /= 0) THEN
3384 0 : IF (rcount(pos) /= scount(pos)) CPABORT("Invalid count.")
3385 0 : rq(rdispl(pos) + 1:rdispl(pos) + rcount(pos)) = rs(sdispl(pos) + 1:sdispl(pos) + scount(pos))
3386 : END IF
3387 0 : CALL mp_waitall(sreq(0:ns - 1))
3388 0 : CALL mp_waitall(rreq(0:nr - 1))
3389 0 : DEALLOCATE (sreq)
3390 0 : DEALLOCATE (rreq)
3391 0 : CALL group%sync()
3392 :
3393 0 : END SUBROUTINE sparse_alltoall
3394 :
3395 : ! **************************************************************************************************
3396 : !> \brief test data structures for equality. It is assumed that if they are
3397 : !> different for one mpi task they are different for all (??)
3398 : !> \param fft_size_1 ...
3399 : !> \param fft_size_2 ...
3400 : !> \param equal ...
3401 : ! **************************************************************************************************
3402 3495633 : SUBROUTINE is_equal(fft_size_1, fft_size_2, equal)
3403 : TYPE(fft_scratch_sizes) :: fft_size_1, fft_size_2
3404 : LOGICAL :: equal
3405 :
3406 : equal = .TRUE.
3407 :
3408 3495633 : equal = equal .AND. fft_size_1%nx == fft_size_2%nx
3409 3495633 : equal = equal .AND. fft_size_1%ny == fft_size_2%ny
3410 3495633 : equal = equal .AND. fft_size_1%nz == fft_size_2%nz
3411 :
3412 3495633 : equal = equal .AND. fft_size_1%lmax == fft_size_2%lmax
3413 3495633 : equal = equal .AND. fft_size_1%mmax == fft_size_2%mmax
3414 3495633 : equal = equal .AND. fft_size_1%nmax == fft_size_2%nmax
3415 :
3416 3495633 : equal = equal .AND. fft_size_1%mx1 == fft_size_2%mx1
3417 3495633 : equal = equal .AND. fft_size_1%mx2 == fft_size_2%mx2
3418 3495633 : equal = equal .AND. fft_size_1%mx3 == fft_size_2%mx3
3419 :
3420 3495633 : equal = equal .AND. fft_size_1%my1 == fft_size_2%my1
3421 3495633 : equal = equal .AND. fft_size_1%my2 == fft_size_2%my2
3422 3495633 : equal = equal .AND. fft_size_1%my3 == fft_size_2%my3
3423 :
3424 3495633 : equal = equal .AND. fft_size_1%mcz1 == fft_size_2%mcz1
3425 3495633 : equal = equal .AND. fft_size_1%mcx2 == fft_size_2%mcx2
3426 3495633 : equal = equal .AND. fft_size_1%mcz2 == fft_size_2%mcz2
3427 3495633 : equal = equal .AND. fft_size_1%mcy3 == fft_size_2%mcy3
3428 :
3429 3495633 : equal = equal .AND. fft_size_1%lg == fft_size_2%lg
3430 3495633 : equal = equal .AND. fft_size_1%mg == fft_size_2%mg
3431 :
3432 3495633 : equal = equal .AND. fft_size_1%nbx == fft_size_2%nbx
3433 3495633 : equal = equal .AND. fft_size_1%nbz == fft_size_2%nbz
3434 :
3435 3495633 : equal = equal .AND. fft_size_1%nmray == fft_size_2%nmray
3436 3495633 : equal = equal .AND. fft_size_1%nyzray == fft_size_2%nyzray
3437 :
3438 3495633 : equal = equal .AND. fft_size_1%rs_group == fft_size_2%rs_group
3439 :
3440 10486899 : equal = equal .AND. ALL(fft_size_1%g_pos == fft_size_2%g_pos)
3441 10486899 : equal = equal .AND. ALL(fft_size_1%r_pos == fft_size_2%r_pos)
3442 10486899 : equal = equal .AND. ALL(fft_size_1%r_dim == fft_size_2%r_dim)
3443 :
3444 3495633 : equal = equal .AND. fft_size_1%numtask == fft_size_2%numtask
3445 :
3446 3495633 : END SUBROUTINE is_equal
3447 :
3448 0 : END MODULE fft_tools
|
