source: palm/trunk/SOURCE/poisfft_mod.f90 @ 4719

Last change on this file since 4719 was 4717, checked in by pavelkrc, 4 years ago

Fixes and optimizations of OpenMP parallelization, formatting of OpenMP directives

  • Property svn:keywords set to Id
File size: 45.7 KB
Line 
1!> @file poisfft_mod.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2020 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: poisfft_mod.f90 4717 2020-09-30 22:27:40Z pavelkrc $
27! Formatting of OpenMP directives (J. Resler)
28!
29! 4671 2020-09-09 20:27:58Z pavelkrc
30! OMP bugfix
31!
32! 4649 2020-08-25 12:11:17Z raasch
33! File re-formatted to follow the PALM coding standard
34!
35! 4429 2020-02-27 15:24:30Z raasch
36! Statements added to avoid compile errors due to unused dummy arguments in serial mode
37!
38! 4366 2020-01-09 08:12:43Z raasch
39! Modification concerning NEC vectorizatio
40!
41! 4360 2020-01-07 11:25:50Z suehring
42! Corrected "Former revisions" section
43!
44! 3690 2019-01-22 22:56:42Z knoop
45! OpenACC port for SPEC
46!
47! Revision 1.1  1997/07/24 11:24:14  raasch
48! Initial revision
49!
50!
51!--------------------------------------------------------------------------------------------------!
52! Description:
53! ------------
54!> Solves the Poisson equation with a 2D spectral method
55!>        d^2 p / dx^2 + d^2 p / dy^2 + d^2 p / dz^2 = s
56!>
57!> Input:
58!> real   ar   contains (nnz,nny,nnx) elements of the velocity divergence, starting from (1,nys,nxl)
59!>
60!> Output:
61!> real   ar   contains the solution for perturbation pressure p
62!--------------------------------------------------------------------------------------------------!
63 MODULE poisfft_mod
64
65
66    USE fft_xy,                                                                                    &
67        ONLY:  fft_init,                                                                           &
68               fft_y,                                                                              &
69               fft_y_1d,                                                                           &
70               fft_y_m,                                                                            &
71               fft_x,                                                                              &
72               fft_x_1d,                                                                           &
73               fft_x_m,                                                                            &
74               temperton_fft_vec
75
76    USE indices,                                                                                   &
77        ONLY:  nnx,                                                                                &
78               nny,                                                                                &
79               nx,                                                                                 &
80               nxl,                                                                                &
81               nxr,                                                                                &
82               ny,                                                                                 &
83               nys,                                                                                &
84               nyn,                                                                                &
85               nz
86
87    USE transpose_indices,                                                                         &
88        ONLY:  nxl_y,                                                                              &
89               nxl_z,                                                                              &
90               nxr_y,                                                                              &
91               nxr_z,                                                                              &
92               nys_x,                                                                              &
93               nys_z,                                                                              &
94               nyn_x,                                                                              &
95               nyn_z,                                                                              &
96               nzb_x,                                                                              &
97               nzb_y,                                                                              &
98               nzt_x,                                                                              &
99               nzt_y
100
101    USE tridia_solver,                                                                             &
102        ONLY:  tridia_1dd,                                                                         &
103               tridia_init,                                                                        &
104               tridia_substi,                                                                      &
105               tridia_substi_overlap
106
107    IMPLICIT NONE
108
109    LOGICAL, SAVE ::  poisfft_initialized = .FALSE.  !<
110
111    PRIVATE
112
113    PUBLIC  poisfft, poisfft_init
114
115    INTERFACE poisfft
116       MODULE PROCEDURE poisfft
117    END INTERFACE poisfft
118
119    INTERFACE poisfft_init
120       MODULE PROCEDURE poisfft_init
121    END INTERFACE poisfft_init
122
123
124 CONTAINS
125
126!--------------------------------------------------------------------------------------------------!
127! Description:
128! ------------
129!> Setup coefficients for FFT and the tridiagonal solver
130!--------------------------------------------------------------------------------------------------!
131 SUBROUTINE poisfft_init
132
133    IMPLICIT NONE
134
135
136    CALL fft_init
137
138    CALL tridia_init
139
140    poisfft_initialized = .TRUE.
141
142 END SUBROUTINE poisfft_init
143
144
145
146!--------------------------------------------------------------------------------------------------!
147! Description:
148! ------------
149!> Two-dimensional Fourier Transformation in x- and y-direction.
150!--------------------------------------------------------------------------------------------------!
151 SUBROUTINE poisfft( ar )
152
153    USE control_parameters,                                                                        &
154        ONLY:  transpose_compute_overlap
155
156    USE cpulog,                                                                                    &
157        ONLY:  cpu_log,                                                                            &
158               cpu_log_nowait,                                                                     &
159               log_point_s
160
161    USE kinds
162
163    USE pegrid
164
165    IMPLICIT NONE
166
167    INTEGER(iwp) ::  ii           !<
168    INTEGER(iwp) ::  iind         !<
169    INTEGER(iwp) ::  inew         !<
170    INTEGER(iwp) ::  jj           !<
171    INTEGER(iwp) ::  jind         !<
172    INTEGER(iwp) ::  jnew         !<
173    INTEGER(iwp) ::  ki           !<
174    INTEGER(iwp) ::  kk           !<
175    INTEGER(iwp) ::  knew         !<
176    INTEGER(iwp) ::  n            !<
177    INTEGER(iwp) ::  nblk         !<
178    INTEGER(iwp) ::  nnx_y        !<
179    INTEGER(iwp) ::  nny_z        !<
180    INTEGER(iwp) ::  nnz_x        !<
181    INTEGER(iwp) ::  nxl_y_bound  !<
182    INTEGER(iwp) ::  nxr_y_bound  !<
183
184    INTEGER(iwp), DIMENSION(4) ::  isave  !<
185
186    REAL(wp), DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar      !<
187    REAL(wp), DIMENSION(nys:nyn,nxl:nxr,1:nz) ::  ar_inv  !<
188
189#define __acc_fft_device ( defined( _OPENACC ) && ( defined ( __cuda_fft ) ) )
190#if __acc_fft_device
191    !$ACC DECLARE CREATE(ar_inv)
192#endif
193
194    REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  ar1      !<
195    REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  f_in     !<
196    REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  f_inv    !<
197    REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  f_out_y  !<
198    REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  f_out_z  !<
199
200
201    CALL cpu_log( log_point_s(3), 'poisfft', 'start' )
202
203    IF ( .NOT. poisfft_initialized )  CALL poisfft_init
204
205#if !__acc_fft_device
206    !$ACC UPDATE HOST(ar)
207#endif
208
209#ifndef _OPENACC
210!
211!-- Two-dimensional Fourier Transformation in x- and y-direction.
212    IF ( pdims(2) == 1  .AND.  pdims(1) > 1 )  THEN
213
214!
215!--    1d-domain-decomposition along x:
216!--    FFT along y and transposition y --> x
217       CALL ffty_tr_yx( ar, ar )
218
219!
220!--    FFT along x, solving the tridiagonal system and backward FFT
221       CALL fftx_tri_fftx( ar )
222
223!
224!--    Transposition x --> y and backward FFT along y
225       CALL tr_xy_ffty( ar, ar )
226
227    ELSEIF ( pdims(1) == 1 .AND. pdims(2) > 1 )  THEN
228
229!
230!--    1d-domain-decomposition along y:
231!--    FFT along x and transposition x --> y
232       CALL fftx_tr_xy( ar, ar )
233
234!
235!--    FFT along y, solving the tridiagonal system and backward FFT
236       CALL ffty_tri_ffty( ar )
237
238!
239!--    Transposition y --> x and backward FFT along x
240       CALL tr_yx_fftx( ar, ar )
241
242    ELSEIF ( .NOT. transpose_compute_overlap )  THEN
243#endif
244
245!
246!--    2d-domain-decomposition or no decomposition (1 PE run)
247!--    Transposition z --> x
248       CALL cpu_log( log_point_s(5), 'transpo forward', 'start' )
249       CALL resort_for_zx( ar, ar_inv )
250       CALL transpose_zx( ar_inv, ar )
251       CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
252
253       CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
254       IF ( temperton_fft_vec )  THEN
255!
256!--       Vector version outputs a transformed array ar_inv that does not require resorting
257!--       (which is done for ar further below)
258          CALL fft_x( ar, 'forward',  ar_inv=ar_inv)
259       ELSE
260          CALL fft_x( ar, 'forward')
261       ENDIF
262       CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
263
264!
265!--    Transposition x --> y
266       CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
267       IF( .NOT. temperton_fft_vec )  CALL resort_for_xy( ar, ar_inv )
268       CALL transpose_xy( ar_inv, ar )
269       CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
270
271       CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
272       IF ( temperton_fft_vec )  THEN
273!
274!--       Input array ar_inv from fft_x can be directly used here.
275!--       The output (also in array ar_inv) does not require resorting below.
276          CALL fft_y( ar, 'forward', ar_inv = ar_inv, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,    &
277                      nxl_y_l = nxl_y, nxr_y_l = nxr_y )
278       ELSE
279          CALL fft_y( ar, 'forward', ar_tr = ar, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,         &
280                      nxl_y_l = nxl_y, nxr_y_l = nxr_y )
281       ENDIF
282       CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
283
284!
285!--    Transposition y --> z
286       CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
287       IF ( .NOT. temperton_fft_vec )  CALL resort_for_yz( ar, ar_inv )
288       CALL transpose_yz( ar_inv, ar )
289       CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
290
291!
292!--    Solve the tridiagonal equation system along z
293       CALL cpu_log( log_point_s(6), 'tridia', 'start' )
294       CALL tridia_substi( ar )
295       CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
296
297!
298!--    Inverse Fourier Transformation
299!--    Transposition z --> y
300       CALL cpu_log( log_point_s(8), 'transpo invers', 'start' )
301       CALL transpose_zy( ar, ar_inv )
302!
303!--    The fft_y below (vector branch) can directly process ar_inv (i.e. does not require a
304!--    resorting)
305       IF ( .NOT. temperton_fft_vec )  CALL resort_for_zy( ar_inv, ar )
306       CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
307
308       CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
309       IF ( temperton_fft_vec )  THEN
310!
311!--       Output array ar_inv can be used as input to the below fft_x routine without resorting
312          CALL fft_y( ar, 'backward', ar_inv = ar_inv, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,   &
313                      nxl_y_l = nxl_y, nxr_y_l = nxr_y )
314       ELSE
315          CALL fft_y( ar, 'backward', ar_tr = ar, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,        &
316                      nxl_y_l = nxl_y, nxr_y_l = nxr_y )
317       ENDIF
318
319       CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
320
321!
322!--    Transposition y --> x
323       CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
324       CALL transpose_yx( ar, ar_inv )
325       IF ( .NOT. temperton_fft_vec )  CALL resort_for_yx( ar_inv, ar )
326       CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
327
328       CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
329       IF ( temperton_fft_vec )  THEN
330          CALL fft_x( ar, 'backward',  ar_inv=ar_inv )
331       ELSE
332          CALL fft_x( ar, 'backward' )
333       ENDIF
334       CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
335
336!
337!--    Transposition x --> z
338       CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
339       CALL transpose_xz( ar, ar_inv )
340       CALL resort_for_xz( ar_inv, ar )
341       CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
342
343#ifndef _OPENACC
344    ELSE
345
346!
347!--    2d-domain-decomposition or no decomposition (1 PE run) with overlapping transposition / fft
348!--    cputime logging must not use barriers, which would prevent overlapping
349       ALLOCATE( f_out_y(0:ny,nxl_y:nxr_y,nzb_y:nzt_y), &
350                 f_out_z(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
351!
352!--    Transposition z --> x + subsequent fft along x
353       ALLOCATE( f_inv(nys:nyn,nxl:nxr,1:nz) )
354       CALL resort_for_zx( ar, f_inv )
355!
356!--    Save original indices and gridpoint counter
357       isave(1) = nz
358       isave(2) = nzb_x
359       isave(3) = nzt_x
360       isave(4) = sendrecvcount_zx
361!
362!--    Set new indices for transformation
363       nblk  = nz / pdims(1)
364       nz    = pdims(1)
365       nnz_x = 1
366       nzb_x = 1 + myidx * nnz_x
367       nzt_x = ( myidx + 1 ) * nnz_x
368       sendrecvcount_zx = nnx * nny * nnz_x
369
370       ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
371       ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) )
372
373       DO  kk = 1, nblk
374
375          IF ( kk == 1 )  THEN
376             CALL cpu_log( log_point_s(5), 'transpo forward', 'start', cpu_log_nowait )
377          ELSE
378             CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
379          ENDIF
380
381          DO  knew = 1, nz
382             ki = kk + nblk * ( knew - 1 )
383             f_in(:,:,knew) = f_inv(:,:,ki)
384          ENDDO
385
386          CALL transpose_zx( f_in, ar1(:,:,:))
387          CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
388
389          IF ( kk == 1 )  THEN
390             CALL cpu_log( log_point_s(4), 'fft_x', 'start', cpu_log_nowait )
391          ELSE
392             CALL cpu_log( log_point_s(4), 'fft_x', 'continue', cpu_log_nowait )
393          ENDIF
394
395          n = isave(2) + kk - 1
396          CALL fft_x( ar1(:,:,:), 'forward',  ar_2d = f_out_z(:,:,n) )
397          CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
398
399       ENDDO
400!
401!--    Restore original indices/counters
402       nz               = isave(1)
403       nzb_x            = isave(2)
404       nzt_x            = isave(3)
405       sendrecvcount_zx = isave(4)
406
407       DEALLOCATE( ar1, f_in, f_inv )
408
409!
410!--    Transposition x --> y + subsequent fft along y
411       ALLOCATE( f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
412       CALL resort_for_xy( f_out_z, f_inv )
413!
414!--    Save original indices and gridpoint counter
415       isave(1) = nx
416       isave(2) = nxl_y
417       isave(3) = nxr_y
418       isave(4) = sendrecvcount_xy
419!
420!--    Set new indices for transformation
421       nblk  = ( ( nx+1 ) / pdims(2) ) - 1
422       nx    = pdims(2)
423       nnx_y = 1
424       nxl_y = myidy * nnx_y
425       nxr_y = ( myidy + 1 ) * nnx_y - 1
426       sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 )
427
428       ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) )
429       ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
430
431       DO  ii = 0, nblk
432
433          CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
434
435          DO  inew = 0, nx-1
436             iind = ii + ( nblk + 1 ) * inew
437             f_in(:,:,inew) = f_inv(:,:,iind)
438          ENDDO
439
440          CALL transpose_xy( f_in, ar1(:,:,:) )
441
442          CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
443
444          IF ( ii == 1 )  THEN
445             CALL cpu_log( log_point_s(7), 'fft_y', 'start', cpu_log_nowait )
446          ELSE
447             CALL cpu_log( log_point_s(7), 'fft_y', 'continue', cpu_log_nowait )
448          ENDIF
449
450          nxl_y_bound = isave(2)
451          nxr_y_bound = isave(3)
452          n           = isave(2) + ii
453          CALL fft_y( ar1(:,:,:), 'forward', ar_tr = f_out_y, nxl_y_bound = nxl_y_bound,           &
454                      nxr_y_bound = nxr_y_bound, nxl_y_l = n, nxr_y_l = n )
455
456          CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
457
458       ENDDO
459!
460!--    Restore original indices/counters
461       nx               = isave(1)
462       nxl_y            = isave(2)
463       nxr_y            = isave(3)
464       sendrecvcount_xy = isave(4)
465
466       DEALLOCATE( ar1, f_in, f_inv )
467
468!
469!--    Transposition y --> z + subsequent tridia + resort for z --> y
470       ALLOCATE( f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) )
471       CALL resort_for_yz( f_out_y, f_inv )
472!
473!--    Save original indices and gridpoint counter
474       isave(1) = ny
475       isave(2) = nys_z
476       isave(3) = nyn_z
477       isave(4) = sendrecvcount_yz
478!
479!--    Set new indices for transformation
480       nblk             = ( ( ny+1 ) / pdims(1) ) - 1
481       ny               = pdims(1)
482       nny_z            = 1
483       nys_z            = myidx * nny_z
484       nyn_z            = ( myidx + 1 ) * nny_z - 1
485       sendrecvcount_yz = ( nxr_y-nxl_y+1 ) * nny_z * ( nzt_y-nzb_y+1 )
486
487       ALLOCATE( ar1(nxl_z:nxr_z,nys_z:nyn_z,1:nz) )
488       ALLOCATE( f_in(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) )
489
490       DO  jj = 0, nblk
491!
492!--       Forward Fourier Transformation
493!--       Transposition y --> z
494          CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
495
496          DO  jnew = 0, ny-1
497             jind = jj + ( nblk + 1 ) * jnew
498             f_in(:,:,jnew) = f_inv(:,:,jind)
499          ENDDO
500
501          CALL transpose_yz( f_in, ar1(:,:,:) )
502
503          IF ( jj == nblk )  THEN
504             CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
505          ELSE
506             CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
507          ENDIF
508
509!
510!--       Solve the tridiagonal equation system along z
511          CALL cpu_log( log_point_s(6), 'tridia', 'start', cpu_log_nowait )
512
513          n = isave(2) + jj
514          CALL tridia_substi_overlap( ar1(:,:,:), n )
515
516          CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
517
518!
519!--       Inverse Fourier Transformation
520!--       Transposition z --> y
521!--       Only one thread should call MPI routines, therefore forward and backward tranpose are in
522!--       the same section
523          IF ( jj == 0 )  THEN
524             CALL cpu_log( log_point_s(8), 'transpo invers', 'start', cpu_log_nowait )
525          ELSE
526             CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
527          ENDIF
528
529          CALL transpose_zy( ar1(:,:,:), f_in )
530
531          DO  jnew = 0, ny-1
532             jind = jj + ( nblk + 1 ) * jnew
533             f_inv(:,:,jind) = f_in(:,:,jnew)
534          ENDDO
535
536          CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
537
538       ENDDO
539!
540!--    Restore original indices/counters
541       ny               = isave(1)
542       nys_z            = isave(2)
543       nyn_z            = isave(3)
544       sendrecvcount_yz = isave(4)
545
546       CALL resort_for_zy( f_inv, f_out_y )
547
548       DEALLOCATE( ar1, f_in, f_inv )
549
550!
551!--    fft along y backward + subsequent transposition y --> x
552       ALLOCATE( f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
553!
554!--    Save original indices and gridpoint counter
555       isave(1) = nx
556       isave(2) = nxl_y
557       isave(3) = nxr_y
558       isave(4) = sendrecvcount_xy
559!
560!--    Set new indices for transformation
561       nblk             = ( ( nx+1 ) / pdims(2) ) - 1
562       nx               = pdims(2)
563       nnx_y            = 1
564       nxl_y            = myidy * nnx_y
565       nxr_y            = ( myidy + 1 ) * nnx_y - 1
566       sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 )
567
568       ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) )
569       ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
570
571       DO  ii = 0, nblk
572
573          CALL cpu_log( log_point_s(7), 'fft_y', 'continue', cpu_log_nowait )
574
575          n = isave(2) + ii
576          nxl_y_bound = isave(2)
577          nxr_y_bound = isave(3)
578
579          CALL fft_y( ar1(:,:,:), 'backward', ar_tr = f_out_y, nxl_y_bound = nxl_y_bound,          &
580                      nxr_y_bound = nxr_y_bound, nxl_y_l = n, nxr_y_l = n )
581
582          IF ( ii == nblk )  THEN
583             CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
584          ELSE
585             CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
586          ENDIF
587
588          CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
589
590          CALL transpose_yx( ar1(:,:,:), f_in )
591
592          DO  inew = 0, nx-1
593             iind = ii + (nblk+1) * inew
594             f_inv(:,:,iind) = f_in(:,:,inew)
595          ENDDO
596
597          CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
598
599       ENDDO
600!
601!--    Restore original indices/counters
602       nx               = isave(1)
603       nxl_y            = isave(2)
604       nxr_y            = isave(3)
605       sendrecvcount_xy = isave(4)
606
607       CALL resort_for_yx( f_inv, f_out_z )
608
609       DEALLOCATE( ar1, f_in, f_inv )
610
611!
612!--    fft along x backward + subsequent final transposition x --> z
613       ALLOCATE( f_inv(nys:nyn,nxl:nxr,1:nz) )
614!
615!--    Save original indices and gridpoint counter
616       isave(1) = nz
617       isave(2) = nzb_x
618       isave(3) = nzt_x
619       isave(4) = sendrecvcount_zx
620!
621!--    Set new indices for transformation
622       nblk             = nz / pdims(1)
623       nz               = pdims(1)
624       nnz_x            = 1
625       nzb_x            = 1 + myidx * nnz_x
626       nzt_x            = ( myidx + 1 ) * nnz_x
627       sendrecvcount_zx = nnx * nny * nnz_x
628
629       ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
630       ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) )
631
632       DO  kk = 1, nblk
633
634          CALL cpu_log( log_point_s(4), 'fft_x', 'continue', cpu_log_nowait )
635
636          n = isave(2) + kk - 1
637          CALL fft_x( ar1(:,:,:), 'backward', f_out_z(:,:,n) )
638
639          IF ( kk == nblk )  THEN
640             CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
641          ELSE
642             CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
643          ENDIF
644
645          CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
646
647          CALL transpose_xz( ar1(:,:,:), f_in )
648
649          DO  knew = 1, nz
650             ki = kk + nblk * ( knew - 1 )
651             f_inv(:,:,ki) = f_in(:,:,knew)
652          ENDDO
653
654          IF ( kk == nblk )  THEN
655             CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
656          ELSE
657             CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
658          ENDIF
659
660       ENDDO
661!
662!--    Restore original indices/counters
663       nz               = isave(1)
664       nzb_x            = isave(2)
665       nzt_x            = isave(3)
666       sendrecvcount_zx = isave(4)
667
668       CALL resort_for_xz( f_inv, ar )
669
670       DEALLOCATE( ar1, f_in, f_inv )
671
672    ENDIF
673#endif
674
675#if !__acc_fft_device
676    !$ACC UPDATE DEVICE(ar)
677#endif
678
679    CALL cpu_log( log_point_s(3), 'poisfft', 'stop' )
680
681 END SUBROUTINE poisfft
682
683
684!--------------------------------------------------------------------------------------------------!
685! Description:
686! ------------
687!> Fourier-transformation along y with subsequent transposition y --> x for a 1d-decomposition
688!> along x.
689!>
690!> @attention The performance of this routine is much faster on the NEC-SX6, if the first index of
691!>            work_ffty_vec is odd. Otherwise memory bank conflicts may occur (especially if the
692!>            index is a multiple of 128). That's why work_ffty_vec is dimensioned as 0:ny+1.
693!>            Of course, this will not work if users are using an odd number of gridpoints along y.
694!--------------------------------------------------------------------------------------------------!
695 SUBROUTINE ffty_tr_yx( f_in, f_out )
696
697    USE control_parameters,                                                                        &
698        ONLY:  loop_optimization
699
700    USE cpulog,                                                                                    &
701        ONLY:  cpu_log,                                                                            &
702               log_point_s
703
704    USE kinds
705
706    USE pegrid
707
708    IMPLICIT NONE
709
710    INTEGER(iwp) ::  i       !<
711    INTEGER(iwp) ::  iend    !<
712    INTEGER(iwp) ::  iouter  !<
713    INTEGER(iwp) ::  ir      !<
714    INTEGER(iwp) ::  j       !<
715    INTEGER(iwp) ::  k       !<
716
717    INTEGER(iwp), PARAMETER ::  stridex = 4  !<
718
719    REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr)             ::  f_in   !<
720    REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) ::  f_out  !<
721    REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny)             ::  work   !<
722
723    REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  work_ffty      !<
724    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  work_ffty_vec  !<
725
726!
727!-- Carry out the FFT along y, where all data are present due to the 1d-decomposition along x.
728!-- Resort the data in a way that x becomes the first index.
729    CALL cpu_log( log_point_s(7), 'fft_y_1d', 'start' )
730
731    IF ( loop_optimization == 'vector' )  THEN
732
733       ALLOCATE( work_ffty_vec(0:ny+1,1:nz,nxl:nxr) )
734!
735!--    Code optimized for vector processors
736       !$OMP PARALLEL PRIVATE ( i, j, k )
737       !$OMP DO
738       DO  i = nxl, nxr
739
740          DO  j = 0, ny
741             DO  k = 1, nz
742                work_ffty_vec(j,k,i) = f_in(k,j,i)
743             ENDDO
744          ENDDO
745
746          CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'forward' )
747
748       ENDDO
749
750       !$OMP DO
751       DO  k = 1, nz
752          DO  j = 0, ny
753             DO  i = nxl, nxr
754                work(i,k,j) = work_ffty_vec(j,k,i)
755             ENDDO
756          ENDDO
757       ENDDO
758       !$OMP END PARALLEL
759
760       DEALLOCATE( work_ffty_vec )
761
762    ELSE
763!
764!--    Cache optimized code.
765       ALLOCATE( work_ffty(0:ny,stridex) )
766!
767!--    The i-(x-)direction is split into a strided outer loop and an inner loop for better cache
768!--    performance
769       !$OMP PARALLEL PRIVATE (i,iend,iouter,ir,j,k,work_ffty)
770       !$OMP DO
771       DO  iouter = nxl, nxr, stridex
772
773          iend = MIN( iouter+stridex-1, nxr )  ! Upper bound for inner i loop
774
775          DO  k = 1, nz
776
777             DO  i = iouter, iend
778
779                ir = i-iouter+1  ! Counter within a stride
780                DO  j = 0, ny
781                   work_ffty(j,ir) = f_in(k,j,i)
782                ENDDO
783!
784!--             FFT along y
785                CALL fft_y_1d( work_ffty(:,ir), 'forward' )
786
787             ENDDO
788
789!
790!--          Resort
791             DO  j = 0, ny
792                DO  i = iouter, iend
793                   work(i,k,j) = work_ffty(j,i-iouter+1)
794                ENDDO
795             ENDDO
796
797          ENDDO
798
799       ENDDO
800       !$OMP END PARALLEL
801
802       DEALLOCATE( work_ffty )
803
804    ENDIF
805
806    CALL cpu_log( log_point_s(7), 'fft_y_1d', 'pause' )
807
808!
809!-- Transpose array
810#if defined( __parallel )
811    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
812    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
813    CALL MPI_ALLTOALL( work(nxl,1,0), sendrecvcount_xy, MPI_REAL, f_out(1,1,nys_x,1),              &
814                       sendrecvcount_xy, MPI_REAL, comm1dx, ierr )
815    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
816#else
817!
818!-- Next line required to avoid compile error about unused dummy argument in serial mode
819    i = SIZE( f_out )
820#endif
821
822 END SUBROUTINE ffty_tr_yx
823
824
825!--------------------------------------------------------------------------------------------------!
826! Description:
827! ------------
828!> Transposition x --> y with a subsequent backward Fourier transformation for a 1d-decomposition
829!> along x
830!--------------------------------------------------------------------------------------------------!
831 SUBROUTINE tr_xy_ffty( f_in, f_out )
832
833    USE control_parameters,                                                                        &
834        ONLY:  loop_optimization
835
836    USE cpulog,                                                                                    &
837        ONLY:  cpu_log,                                                                            &
838               log_point_s
839
840    USE kinds
841
842    USE pegrid
843
844    IMPLICIT NONE
845
846    INTEGER(iwp) ::  i       !<
847    INTEGER(iwp) ::  iend    !<
848    INTEGER(iwp) ::  iouter  !<
849    INTEGER(iwp) ::  ir      !<
850    INTEGER(iwp) ::  j       !<
851    INTEGER(iwp) ::  k       !<
852
853    INTEGER(iwp), PARAMETER ::  stridex = 4  !<
854
855    REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) ::  f_in   !<
856    REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr)             ::  f_out  !<
857    REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny)             ::  work   !<
858
859    REAL(wp), DIMENSION(:,:),   ALLOCATABLE ::  work_ffty      !<
860    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  work_ffty_vec  !<
861
862!
863!-- Transpose array
864#if defined( __parallel )
865    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
866    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
867    CALL MPI_ALLTOALL( f_in(1,1,nys_x,1), sendrecvcount_xy, MPI_REAL, work(nxl,1,0),               &
868                       sendrecvcount_xy, MPI_REAL, comm1dx, ierr )
869    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
870#else
871!
872!-- Next line required to avoid compile error about unused dummy argument in serial mode
873    i = SIZE( f_in )
874#endif
875
876!
877!-- Resort the data in a way that y becomes the first index and carry out the backward fft along y.
878    CALL cpu_log( log_point_s(7), 'fft_y_1d', 'continue' )
879
880    IF ( loop_optimization == 'vector' )  THEN
881
882       ALLOCATE( work_ffty_vec(0:ny+1,1:nz,nxl:nxr) )
883!
884!--    Code optimized for vector processors
885       !$OMP PARALLEL PRIVATE ( i, j, k )
886       !$OMP DO
887       DO  k = 1, nz
888          DO  j = 0, ny
889             DO  i = nxl, nxr
890                work_ffty_vec(j,k,i) = work(i,k,j)
891             ENDDO
892          ENDDO
893       ENDDO
894
895       !$OMP DO
896       DO  i = nxl, nxr
897
898          CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'backward' )
899
900          DO  j = 0, ny
901             DO  k = 1, nz
902                f_out(k,j,i) = work_ffty_vec(j,k,i)
903             ENDDO
904          ENDDO
905
906       ENDDO
907       !$OMP END PARALLEL
908
909       DEALLOCATE( work_ffty_vec )
910
911    ELSE
912!
913!--    Cache optimized code.
914       ALLOCATE( work_ffty(0:ny,stridex) )
915!
916!--    The i-(x-)direction is split into a strided outer loop and an inner loop for better cache
917!--    performance
918       !$OMP PARALLEL PRIVATE ( i, iend, iouter, ir, j, k, work_ffty )
919       !$OMP DO
920       DO  iouter = nxl, nxr, stridex
921
922          iend = MIN( iouter+stridex-1, nxr )  ! Upper bound for inner i loop
923
924          DO  k = 1, nz
925!
926!--          Resort
927             DO  j = 0, ny
928                DO  i = iouter, iend
929                   work_ffty(j,i-iouter+1) = work(i,k,j)
930                ENDDO
931             ENDDO
932
933             DO  i = iouter, iend
934
935!
936!--             FFT along y
937                ir = i-iouter+1  ! Counter within a stride
938                CALL fft_y_1d( work_ffty(:,ir), 'backward' )
939
940                DO  j = 0, ny
941                   f_out(k,j,i) = work_ffty(j,ir)
942                ENDDO
943             ENDDO
944
945          ENDDO
946
947       ENDDO
948       !$OMP END PARALLEL
949
950       DEALLOCATE( work_ffty )
951
952    ENDIF
953
954    CALL cpu_log( log_point_s(7), 'fft_y_1d', 'stop' )
955
956 END SUBROUTINE tr_xy_ffty
957
958
959!--------------------------------------------------------------------------------------------------!
960! Description:
961! ------------
962!> FFT along x, solution of the tridiagonal system and backward FFT for a 1d-decomposition along x
963!>
964!> @warning This subroutine may still not work for hybrid parallelization with OpenMP (for possible
965!>          necessary changes see the original routine poisfft_hybrid, developed by Klaus Ketelsen,
966!>          May 2002)
967!--------------------------------------------------------------------------------------------------!
968 SUBROUTINE fftx_tri_fftx( ar )
969
970    USE control_parameters,                                                                        &
971        ONLY:  loop_optimization
972
973    USE cpulog,                                                                                    &
974        ONLY:  cpu_log,                                                                            &
975               log_point_s
976
977    USE grid_variables,                                                                            &
978        ONLY:  ddx2,                                                                               &
979               ddy2
980
981    USE kinds
982
983    USE pegrid
984
985    IMPLICIT NONE
986
987    INTEGER(iwp) ::  i                   !<
988    INTEGER(iwp) ::  j                   !<
989    INTEGER(iwp) ::  k                   !<
990    INTEGER(iwp) ::  m                   !<
991    INTEGER(iwp) ::  n                   !<
992!$  INTEGER(iwp) ::  omp_get_thread_num  !<
993    INTEGER(iwp) ::  tn                  !<
994
995    REAL(wp), DIMENSION(0:nx)                          ::  work_fftx  !<
996    REAL(wp), DIMENSION(0:nx,1:nz)                     ::  work_trix  !<
997    REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) ::  ar         !<
998    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE          ::  tri        !<
999
1000
1001    CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'start' )
1002
1003    ALLOCATE( tri(5,0:nx,0:nz-1,0:threads_per_task-1) )
1004
1005    tn = 0              ! Default thread number in case of one thread
1006    !$OMP  PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_fftx, work_trix )
1007    DO  j = nys_x, nyn_x
1008
1009!$     tn = omp_get_thread_num()
1010
1011       IF ( loop_optimization == 'vector' )  THEN
1012!
1013!--       Code optimized for vector processors
1014          DO  k = 1, nz
1015
1016             m = 0
1017             DO  n = 1, pdims(1)
1018                DO  i = 1, nnx
1019                   work_trix(m,k) = ar(i,k,j,n)
1020                   m = m + 1
1021                ENDDO
1022             ENDDO
1023
1024          ENDDO
1025
1026          CALL fft_x_m( work_trix, 'forward' )
1027
1028       ELSE
1029!
1030!--       Cache optimized code
1031          DO  k = 1, nz
1032
1033             m = 0
1034             DO  n = 1, pdims(1)
1035                DO  i = 1, nnx
1036                   work_fftx(m) = ar(i,k,j,n)
1037                   m = m + 1
1038                ENDDO
1039             ENDDO
1040
1041             CALL fft_x_1d( work_fftx, 'forward' )
1042
1043             DO  i = 0, nx
1044                work_trix(i,k) = work_fftx(i)
1045             ENDDO
1046
1047          ENDDO
1048
1049       ENDIF
1050
1051!
1052!--    Solve the linear equation system
1053       CALL tridia_1dd( ddx2, ddy2, nx, ny, j, work_trix, tri(:,:,:,tn) )
1054
1055       IF ( loop_optimization == 'vector' )  THEN
1056!
1057!--       Code optimized for vector processors
1058          CALL fft_x_m( work_trix, 'backward' )
1059
1060          DO  k = 1, nz
1061
1062             m = 0
1063             DO  n = 1, pdims(1)
1064                DO  i = 1, nnx
1065                   ar(i,k,j,n) = work_trix(m,k)
1066                   m = m + 1
1067                ENDDO
1068             ENDDO
1069
1070          ENDDO
1071
1072       ELSE
1073!
1074!--       Cache optimized code
1075          DO  k = 1, nz
1076
1077             DO  i = 0, nx
1078                work_fftx(i) = work_trix(i,k)
1079             ENDDO
1080
1081             CALL fft_x_1d( work_fftx, 'backward' )
1082
1083             m = 0
1084             DO  n = 1, pdims(1)
1085                DO  i = 1, nnx
1086                   ar(i,k,j,n) = work_fftx(m)
1087                   m = m + 1
1088                ENDDO
1089             ENDDO
1090
1091          ENDDO
1092
1093       ENDIF
1094
1095    ENDDO
1096
1097    DEALLOCATE( tri )
1098
1099    CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'stop' )
1100
1101 END SUBROUTINE fftx_tri_fftx
1102
1103
1104!--------------------------------------------------------------------------------------------------!
1105! Description:
1106! ------------
1107!> Fourier-transformation along x with subsequent transposition x --> y for a 1d-decomposition
1108!> along y.
1109!>
1110!> @attention NEC-branch of this routine may significantly profit from further optimizations. So
1111!>            far, performance is much worse than for routine ffty_tr_yx (more than three times
1112!>            slower).
1113!--------------------------------------------------------------------------------------------------!
1114 SUBROUTINE fftx_tr_xy( f_in, f_out )
1115
1116
1117    USE control_parameters,                                                                        &
1118        ONLY:  loop_optimization
1119
1120    USE cpulog,                                                                                    &
1121        ONLY:  cpu_log,                                                                            &
1122               log_point_s
1123
1124    USE kinds
1125
1126    USE pegrid
1127
1128    IMPLICIT NONE
1129
1130    INTEGER(iwp) ::  i  !<
1131    INTEGER(iwp) ::  j  !<
1132    INTEGER(iwp) ::  k  !<
1133
1134    REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn)             ::  work_fftx  !<
1135    REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx)             ::  f_in       !<
1136    REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) ::  f_out      !<
1137    REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx)             ::  work       !<
1138
1139!
1140!--  Carry out the FFT along x, where all data are present due to the 1d-decomposition along y.
1141!--  Resort the data in a way that y becomes the first index.
1142     CALL cpu_log( log_point_s(4), 'fft_x_1d', 'start' )
1143
1144     IF ( loop_optimization == 'vector' )  THEN
1145!
1146!--    Code for vector processors
1147       !$OMP  PARALLEL PRIVATE ( i, j, k )
1148       !$OMP  DO
1149       DO  i = 0, nx
1150
1151          DO  j = nys, nyn
1152             DO  k = 1, nz
1153                work_fftx(i,k,j) = f_in(k,j,i)
1154             ENDDO
1155          ENDDO
1156
1157       ENDDO
1158
1159       !$OMP  DO
1160       DO  j = nys, nyn
1161
1162          CALL fft_x_m( work_fftx(:,:,j), 'forward' )
1163
1164          DO  k = 1, nz
1165             DO  i = 0, nx
1166                work(j,k,i) = work_fftx(i,k,j)
1167             ENDDO
1168          ENDDO
1169
1170       ENDDO
1171       !$OMP  END PARALLEL
1172
1173    ELSE
1174
1175!
1176!--    Cache optimized code (there might still be a potential for better optimization).
1177       !$OMP  PARALLEL PRIVATE (i,j,k)
1178       !$OMP  DO
1179       DO  i = 0, nx
1180
1181          DO  j = nys, nyn
1182             DO  k = 1, nz
1183                work_fftx(i,k,j) = f_in(k,j,i)
1184             ENDDO
1185          ENDDO
1186
1187       ENDDO
1188
1189       !$OMP  DO
1190       DO  j = nys, nyn
1191          DO  k = 1, nz
1192
1193             CALL fft_x_1d( work_fftx(0:nx,k,j), 'forward' )
1194
1195             DO  i = 0, nx
1196                work(j,k,i) = work_fftx(i,k,j)
1197             ENDDO
1198          ENDDO
1199
1200       ENDDO
1201       !$OMP  END PARALLEL
1202
1203    ENDIF
1204    CALL cpu_log( log_point_s(4), 'fft_x_1d', 'pause' )
1205
1206!
1207!-- Transpose array
1208#if defined( __parallel )
1209    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1210    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1211    CALL MPI_ALLTOALL( work(nys,1,0), sendrecvcount_xy, MPI_REAL, f_out(1,1,nxl_y,1),              &
1212                       sendrecvcount_xy, MPI_REAL, comm1dy, ierr )
1213    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1214#else
1215!
1216!-- Next line required to avoid compile error about unused dummy argument in serial mode
1217    i = SIZE( f_out )
1218#endif
1219
1220 END SUBROUTINE fftx_tr_xy
1221
1222
1223!--------------------------------------------------------------------------------------------------!
1224! Description:
1225! ------------
1226!> Transposition y --> x with a subsequent backward Fourier transformation for a 1d-decomposition
1227!> along x.
1228!--------------------------------------------------------------------------------------------------!
1229 SUBROUTINE tr_yx_fftx( f_in, f_out )
1230
1231
1232    USE control_parameters,                                                                        &
1233        ONLY:  loop_optimization
1234
1235    USE cpulog,                                                                                    &
1236        ONLY:  cpu_log,                                                                            &
1237               log_point_s
1238
1239    USE kinds
1240
1241    USE pegrid
1242
1243    IMPLICIT NONE
1244
1245    INTEGER(iwp) ::  i  !<
1246    INTEGER(iwp) ::  j  !<
1247    INTEGER(iwp) ::  k  !<
1248
1249    REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) ::  f_in       !<
1250    REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx)             ::  f_out      !<
1251    REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx)             ::  work       !<
1252    REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn)             ::  work_fftx  !<
1253
1254
1255!
1256!-- Transpose array
1257#if defined( __parallel )
1258    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1259    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1260    CALL MPI_ALLTOALL( f_in(1,1,nxl_y,1), sendrecvcount_xy, MPI_REAL, work(nys,1,0),               &
1261                       sendrecvcount_xy, MPI_REAL, comm1dy, ierr )
1262    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1263#else
1264!
1265!-- Next line required to avoid compile error about unused dummy argument in serial mode
1266    i = SIZE( f_in )
1267#endif
1268
1269!
1270!-- Carry out the FFT along x, where all data are present due to the 1d-decomposition along y.
1271!-- Resort the data in a way that y becomes the first index.
1272    CALL cpu_log( log_point_s(4), 'fft_x_1d', 'continue' )
1273
1274    IF ( loop_optimization == 'vector' )  THEN
1275!
1276!--    Code optimized for vector processors
1277!$OMP  PARALLEL PRIVATE ( i, j, k )
1278!$OMP  DO
1279       DO  j = nys, nyn
1280
1281          DO  k = 1, nz
1282             DO  i = 0, nx
1283                work_fftx(i,k,j) = work(j,k,i)
1284             ENDDO
1285          ENDDO
1286
1287          CALL fft_x_m( work_fftx(:,:,j), 'backward' )
1288
1289       ENDDO
1290
1291!$OMP  DO
1292       DO  i = 0, nx
1293          DO  j = nys, nyn
1294             DO  k = 1, nz
1295                f_out(k,j,i) = work_fftx(i,k,j)
1296             ENDDO
1297          ENDDO
1298       ENDDO
1299!$OMP  END PARALLEL
1300
1301    ELSE
1302
1303!
1304!--    Cache optimized code (there might be still a potential for better optimization).
1305!$OMP  PARALLEL PRIVATE (i,j,k)
1306!$OMP  DO
1307       DO  j = nys, nyn
1308          DO  k = 1, nz
1309
1310             DO  i = 0, nx
1311                work_fftx(i,k,j) = work(j,k,i)
1312             ENDDO
1313
1314             CALL fft_x_1d( work_fftx(0:nx,k,j), 'backward' )
1315
1316          ENDDO
1317       ENDDO
1318
1319!$OMP  DO
1320       DO  i = 0, nx
1321          DO  j = nys, nyn
1322             DO  k = 1, nz
1323                f_out(k,j,i) = work_fftx(i,k,j)
1324             ENDDO
1325          ENDDO
1326       ENDDO
1327!$OMP  END PARALLEL
1328
1329    ENDIF
1330    CALL cpu_log( log_point_s(4), 'fft_x_1d', 'stop' )
1331
1332 END SUBROUTINE tr_yx_fftx
1333
1334
1335!--------------------------------------------------------------------------------------------------!
1336! Description:
1337! ------------
1338!> FFT along y, solution of the tridiagonal system and backward FFT for a 1d-decomposition along y.
1339!>
1340!> @warning This subroutine may still not work for hybrid parallelization with OpenMP (for possible
1341!>          necessary changes see the original routine poisfft_hybrid, developed by Klaus Ketelsen,
1342!>          May 2002)
1343!--------------------------------------------------------------------------------------------------!
1344 SUBROUTINE ffty_tri_ffty( ar )
1345
1346
1347    USE control_parameters,                                                                        &
1348        ONLY:  loop_optimization
1349
1350    USE cpulog,                                                                                    &
1351        ONLY:  cpu_log,                                                                            &
1352               log_point_s
1353
1354    USE grid_variables,                                                                            &
1355        ONLY:  ddx2,                                                                               &
1356               ddy2
1357
1358    USE kinds
1359
1360    USE pegrid
1361
1362    IMPLICIT NONE
1363
1364    INTEGER(iwp) ::  i                   !<
1365    INTEGER(iwp) ::  j                   !<
1366    INTEGER(iwp) ::  k                   !<
1367    INTEGER(iwp) ::  m                   !<
1368    INTEGER(iwp) ::  n                   !<
1369!$  INTEGER(iwp) ::  omp_get_thread_num  !<
1370    INTEGER(iwp) ::  tn                  !<
1371
1372    REAL(wp), DIMENSION(0:ny)                          ::  work_ffty  !<
1373    REAL(wp), DIMENSION(0:ny,1:nz)                     ::  work_triy  !<
1374    REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) ::  ar         !<
1375    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE          ::  tri        !<
1376
1377
1378    CALL cpu_log( log_point_s(39), 'fft_y_1d + tridia', 'start' )
1379
1380    ALLOCATE( tri(5,0:ny,0:nz-1,0:threads_per_task-1) )
1381
1382    tn = 0           ! Default thread number in case of one thread
1383!$OMP PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_ffty, work_triy )
1384    DO  i = nxl_y, nxr_y
1385
1386!$     tn = omp_get_thread_num()
1387
1388       IF ( loop_optimization == 'vector' )  THEN
1389!
1390!--       Code optimized for vector processors
1391          DO  k = 1, nz
1392
1393             m = 0
1394             DO  n = 1, pdims(2)
1395                DO  j = 1, nny
1396                   work_triy(m,k) = ar(j,k,i,n)
1397                   m = m + 1
1398                ENDDO
1399             ENDDO
1400
1401          ENDDO
1402
1403          CALL fft_y_m( work_triy, ny, 'forward' )
1404
1405       ELSE
1406!
1407!--       Cache optimized code
1408          DO  k = 1, nz
1409
1410             m = 0
1411             DO  n = 1, pdims(2)
1412                DO  j = 1, nny
1413                   work_ffty(m) = ar(j,k,i,n)
1414                   m = m + 1
1415                ENDDO
1416             ENDDO
1417
1418             CALL fft_y_1d( work_ffty, 'forward' )
1419
1420             DO  j = 0, ny
1421                work_triy(j,k) = work_ffty(j)
1422             ENDDO
1423
1424          ENDDO
1425
1426       ENDIF
1427
1428!
1429!--    Solve the linear equation system
1430       CALL tridia_1dd( ddy2, ddx2, ny, nx, i, work_triy, tri(:,:,:,tn) )
1431
1432       IF ( loop_optimization == 'vector' )  THEN
1433!
1434!--       Code optimized for vector processors
1435          CALL fft_y_m( work_triy, ny, 'backward' )
1436
1437          DO  k = 1, nz
1438
1439             m = 0
1440             DO  n = 1, pdims(2)
1441                DO  j = 1, nny
1442                   ar(j,k,i,n) = work_triy(m,k)
1443                   m = m + 1
1444                ENDDO
1445             ENDDO
1446
1447          ENDDO
1448
1449       ELSE
1450!
1451!--       Cache optimized code
1452          DO  k = 1, nz
1453
1454             DO  j = 0, ny
1455                work_ffty(j) = work_triy(j,k)
1456             ENDDO
1457
1458             CALL fft_y_1d( work_ffty, 'backward' )
1459
1460             m = 0
1461             DO  n = 1, pdims(2)
1462                DO  j = 1, nny
1463                   ar(j,k,i,n) = work_ffty(m)
1464                   m = m + 1
1465                ENDDO
1466             ENDDO
1467
1468          ENDDO
1469
1470       ENDIF
1471
1472    ENDDO
1473
1474    DEALLOCATE( tri )
1475
1476    CALL cpu_log( log_point_s(39), 'fft_y_1d + tridia', 'stop' )
1477
1478 END SUBROUTINE ffty_tri_ffty
1479
1480 END MODULE poisfft_mod
Note: See TracBrowser for help on using the repository browser.