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

Last change on this file since 4397 was 4366, checked in by raasch, 4 years ago

code vectorization for NEC Aurora: vectorized version of Temperton FFT, vectorization of Newtor iteration for calculating the Obukhov length

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