source: palm/trunk/SOURCE/poisfft.f90 @ 1727

Last change on this file since 1727 was 1683, checked in by knoop, 9 years ago

last commit documented

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