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