source: palm/trunk/SOURCE/transpose.f90 @ 3694

Last change on this file since 3694 was 3694, checked in by knoop, 2 years ago

Bugfix: moved OpenACC specific cpp #endif , so cpu_log call is not affected.

  • Property svn:keywords set to Id
File size: 33.9 KB
Line 
1!> @file transpose.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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: transpose.f90 3694 2019-01-23 17:01:49Z knoop $
27! OpenACC port for SPEC
28!
29! 3241 2018-09-12 15:02:00Z raasch
30! unused variables removed
31!
32! 2718 2018-01-02 08:49:38Z maronga
33! Corrected "Former revisions" section
34!
35! 2696 2017-12-14 17:12:51Z kanani
36! Change in file header (GPL part)
37!
38! 2119 2017-01-17 16:51:50Z raasch
39!
40! 2118 2017-01-17 16:38:49Z raasch
41! OpenACC directives removed
42!
43! 2000 2016-08-20 18:09:15Z knoop
44! Forced header and separation lines into 80 columns
45!
46! 1682 2015-10-07 23:56:08Z knoop
47! Code annotations made doxygen readable
48!
49! 1324 2014-03-21 09:13:16Z suehring
50! Bugfix: ONLY statement for module pegrid removed
51!
52! 1320 2014-03-20 08:40:49Z raasch
53! ONLY-attribute added to USE-statements,
54! kind-parameters added to all INTEGER and REAL declaration statements,
55! kinds are defined in new module kinds,
56! old module precision_kind is removed,
57! revision history before 2012 removed,
58! comment fields (!:) to be used for variable explanations added to
59! all variable declaration statements
60!
61! 1318 2014-03-17 13:35:16Z raasch
62! cpu_log_nowait parameter added to cpu measurements of the transpositions
63! required for solving the Poisson equation (poisfft),
64! module interfaces removed
65!
66! 1257 2013-11-08 15:18:40Z raasch
67! openacc loop and loop vector clauses removed
68!
69! 1216 2013-08-26 09:31:42Z raasch
70! re-sorting of the transposed / to be transposed arrays moved to separate
71! routines resort_for_...
72!
73! 1111 2013-03-08 23:54:10Z raasch
74! openACC directives added,
75! resorting data from/to work changed, work got 4 dimensions instead of 1
76!
77! 1106 2013-03-04 05:31:38Z raasch
78! preprocessor lines rearranged so that routines can also be used in serial
79! (non-parallel) mode
80!
81! 1092 2013-02-02 11:24:22Z raasch
82! unused variables removed
83!
84! 1036 2012-10-22 13:43:42Z raasch
85! code put under GPL (PALM 3.9)
86!
87! 1003 2012-09-14 14:35:53Z raasch
88! indices nxa, nya, etc. replaced by nx, ny, etc.
89!
90! Revision 1.1  1997/07/24 11:25:18  raasch
91! Initial revision
92!
93
94#define __acc_fft_device ( defined( _OPENACC ) && ( defined ( __cuda_fft ) ) )
95
96!------------------------------------------------------------------------------!
97! Description:
98! ------------
99!> Resorting data for the transposition from x to y. The transposition itself
100!> is carried out in transpose_xy
101!------------------------------------------------------------------------------!
102 SUBROUTINE resort_for_xy( f_in, f_inv )
103 
104
105     USE indices,                                                              &
106         ONLY:  nx
107
108     USE kinds
109
110     USE transpose_indices,                                                    &
111         ONLY:  nyn_x, nys_x, nzb_x, nzt_x
112
113     IMPLICIT NONE
114
115     REAL(wp) ::  f_in(0:nx,nys_x:nyn_x,nzb_x:nzt_x)  !<
116     REAL(wp) ::  f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) !<
117
118
119     INTEGER(iwp) ::  i !<
120     INTEGER(iwp) ::  j !<
121     INTEGER(iwp) ::  k !<
122!
123!-- Rearrange indices of input array in order to make data to be send
124!-- by MPI contiguous
125    !$OMP  PARALLEL PRIVATE ( i, j, k )
126    !$OMP  DO
127#if __acc_fft_device
128     !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
129     !$ACC PRESENT(f_inv, f_in)
130#endif
131     DO  i = 0, nx
132         DO  k = nzb_x, nzt_x
133             DO  j = nys_x, nyn_x
134                 f_inv(j,k,i) = f_in(i,j,k)
135             ENDDO
136         ENDDO
137     ENDDO
138     !$OMP  END PARALLEL
139
140 END SUBROUTINE resort_for_xy
141
142
143!------------------------------------------------------------------------------!
144! Description:
145! ------------
146!> Transposition of input array (f_in) from x to y. For the input array, all
147!> elements along x reside on the same PE, while after transposition, all
148!> elements along y reside on the same PE.
149!------------------------------------------------------------------------------!
150 SUBROUTINE transpose_xy( f_inv, f_out )
151
152
153    USE cpulog,                                                                &
154        ONLY:  cpu_log, cpu_log_nowait, log_point_s
155
156    USE indices,                                                               &
157        ONLY:  nx, ny
158       
159    USE kinds
160
161    USE pegrid
162
163    USE transpose_indices,                                                     &
164        ONLY:  nxl_y, nxr_y, nyn_x, nys_x, nzb_x, nzb_y, nzt_x, nzt_y
165
166    IMPLICIT NONE
167
168    INTEGER(iwp) ::  i  !<
169    INTEGER(iwp) ::  j  !<
170    INTEGER(iwp) ::  k  !<
171    INTEGER(iwp) ::  l  !<
172    INTEGER(iwp) ::  ys !<
173 
174    REAL(wp) ::  f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) !<
175    REAL(wp) ::  f_out(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) !<
176
177    REAL(wp), DIMENSION(nyn_x-nys_x+1,nzb_y:nzt_y,nxl_y:nxr_y,0:pdims(2)-1) ::  work !<
178#if __acc_fft_device
179    !$ACC DECLARE CREATE(work)
180#endif
181
182
183    IF ( numprocs /= 1 )  THEN
184
185#if defined( __parallel )
186!
187!--    Transpose array
188       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start', cpu_log_nowait )
189
190#if __acc_fft_device
191#ifndef __cuda_aware_mpi
192       !$ACC UPDATE HOST(f_inv)
193#else
194       !$ACC HOST_DATA USE_DEVICE(work, f_inv)
195#endif
196#endif
197
198       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
199       CALL MPI_ALLTOALL( f_inv(nys_x,nzb_x,0),  sendrecvcount_xy, MPI_REAL, &
200                          work(1,nzb_y,nxl_y,0), sendrecvcount_xy, MPI_REAL, &
201                          comm1dy, ierr )
202
203#if __acc_fft_device
204#ifndef __cuda_aware_mpi
205       !$ACC UPDATE DEVICE(work)
206#else
207       !$ACC END HOST_DATA
208#endif
209#endif
210
211       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
212
213!
214!--    Reorder transposed array
215!$OMP  PARALLEL PRIVATE ( i, j, k, l, ys )
216!$OMP  DO
217       DO  l = 0, pdims(2) - 1
218          ys = 0 + l * ( nyn_x - nys_x + 1 )
219#if __acc_fft_device
220          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
221          !$ACC PRESENT(f_out, work)
222#endif
223          DO  i = nxl_y, nxr_y
224             DO  k = nzb_y, nzt_y
225                DO  j = ys, ys + nyn_x - nys_x
226                   f_out(j,i,k) = work(j-ys+1,k,i,l)
227                ENDDO
228             ENDDO
229          ENDDO
230       ENDDO
231!$OMP  END PARALLEL
232#endif
233
234    ELSE
235
236!
237!--    Reorder transposed array
238!$OMP  PARALLEL PRIVATE ( i, j, k )
239!$OMP  DO
240#if __acc_fft_device
241       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
242       !$ACC PRESENT(f_out, f_inv)
243#endif
244       DO  k = nzb_y, nzt_y
245          DO  i = nxl_y, nxr_y
246             DO  j = 0, ny
247                f_out(j,i,k) = f_inv(j,k,i)
248             ENDDO
249          ENDDO
250       ENDDO
251!$OMP  END PARALLEL
252
253    ENDIF
254
255 END SUBROUTINE transpose_xy
256
257
258!------------------------------------------------------------------------------!
259! Description:
260! ------------
261!> Resorting data after the transposition from x to z. The transposition itself
262!> is carried out in transpose_xz
263!------------------------------------------------------------------------------!
264 SUBROUTINE resort_for_xz( f_inv, f_out )
265
266
267     USE indices,                                                              &
268         ONLY:  nxl, nxr, nyn, nys, nz
269
270     USE kinds
271
272     IMPLICIT NONE
273
274     REAL(wp) ::  f_inv(nys:nyn,nxl:nxr,1:nz) !<
275     REAL(wp) ::  f_out(1:nz,nys:nyn,nxl:nxr) !<
276
277     INTEGER(iwp) ::  i !<
278     INTEGER(iwp) ::  j !<
279     INTEGER(iwp) ::  k !<
280!
281!-- Rearrange indices of input array in order to make data to be send
282!-- by MPI contiguous.
283!-- In case of parallel fft/transposition, scattered store is faster in
284!-- backward direction!!!
285    !$OMP  PARALLEL PRIVATE ( i, j, k )
286    !$OMP  DO
287#if __acc_fft_device
288     !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
289     !$ACC PRESENT(f_out, f_inv)
290#endif
291     DO  k = 1, nz
292         DO  i = nxl, nxr
293             DO  j = nys, nyn
294                 f_out(k,j,i) = f_inv(j,i,k)
295             ENDDO
296         ENDDO
297     ENDDO
298     !$OMP  END PARALLEL
299
300 END SUBROUTINE resort_for_xz
301
302
303!------------------------------------------------------------------------------!
304! Description:
305! ------------
306!> Transposition of input array (f_in) from x to z. For the input array, all
307!> elements along x reside on the same PE, while after transposition, all
308!> elements along z reside on the same PE.
309!------------------------------------------------------------------------------!
310 SUBROUTINE transpose_xz( f_in, f_inv )
311
312
313    USE cpulog,                                                                &
314        ONLY:  cpu_log, cpu_log_nowait, log_point_s
315
316    USE indices,                                                               &
317        ONLY:  nnx, nx, nxl, nxr, nyn, nys, nz
318
319    USE kinds
320
321    USE pegrid
322
323    USE transpose_indices,                                                     &
324        ONLY:  nyn_x, nys_x, nzb_x, nzt_x
325
326    IMPLICIT NONE
327
328    INTEGER(iwp) ::  i  !<
329    INTEGER(iwp) ::  j  !<
330    INTEGER(iwp) ::  k  !<
331    INTEGER(iwp) ::  l  !<
332    INTEGER(iwp) ::  xs !<
333
334    REAL(wp) ::  f_in(0:nx,nys_x:nyn_x,nzb_x:nzt_x) !<
335    REAL(wp) ::  f_inv(nys:nyn,nxl:nxr,1:nz) !<
336
337    REAL(wp), DIMENSION(nys_x:nyn_x,nnx,nzb_x:nzt_x,0:pdims(1)-1) ::  work !<
338#if __acc_fft_device
339    !$ACC DECLARE CREATE(work)
340#endif
341
342
343!
344!-- If the PE grid is one-dimensional along y, the array has only to be
345!-- reordered locally and therefore no transposition has to be done.
346    IF ( pdims(1) /= 1 )  THEN
347
348#if defined( __parallel )
349!
350!--    Reorder input array for transposition
351!$OMP  PARALLEL PRIVATE ( i, j, k, l, xs )
352!$OMP  DO
353       DO  l = 0, pdims(1) - 1
354          xs = 0 + l * nnx
355#if __acc_fft_device
356          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
357          !$ACC PRESENT(work, f_in)
358#endif
359          DO  k = nzb_x, nzt_x
360             DO  i = xs, xs + nnx - 1
361                DO  j = nys_x, nyn_x
362                   work(j,i-xs+1,k,l) = f_in(i,j,k)
363                ENDDO
364             ENDDO
365          ENDDO
366       ENDDO
367!$OMP  END PARALLEL
368
369!
370!--    Transpose array
371       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start', cpu_log_nowait )
372
373#if __acc_fft_device
374#ifndef __cuda_aware_mpi
375       !$ACC UPDATE HOST(work)
376#else
377       !$ACC HOST_DATA USE_DEVICE(work, f_inv)
378#endif
379#endif
380
381       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
382       CALL MPI_ALLTOALL( work(nys_x,1,nzb_x,0), sendrecvcount_zx, MPI_REAL, &
383                          f_inv(nys,nxl,1),      sendrecvcount_zx, MPI_REAL, &
384                          comm1dx, ierr )
385
386#if __acc_fft_device
387#ifndef __cuda_aware_mpi
388       !$ACC UPDATE DEVICE(f_inv)
389#else
390       !$ACC END HOST_DATA
391#endif
392#endif
393
394       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
395#endif
396
397    ELSE
398
399!
400!--    Reorder the array in a way that the z index is in first position
401!$OMP  PARALLEL PRIVATE ( i, j, k )
402!$OMP  DO
403#if __acc_fft_device
404       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
405       !$ACC PRESENT(f_inv, f_in)
406#endif
407       DO  i = nxl, nxr
408          DO  j = nys, nyn
409             DO  k = 1, nz
410                f_inv(j,i,k) = f_in(i,j,k)
411             ENDDO
412          ENDDO
413       ENDDO
414!$OMP  END PARALLEL
415
416    ENDIF
417
418 END SUBROUTINE transpose_xz
419
420
421!------------------------------------------------------------------------------!
422! Description:
423! ------------
424!> Resorting data after the transposition from y to x. The transposition itself
425!> is carried out in transpose_yx
426!------------------------------------------------------------------------------!
427 SUBROUTINE resort_for_yx( f_inv, f_out )
428
429
430     USE indices,                                                              &
431         ONLY:  nx
432
433     USE kinds
434
435     USE transpose_indices,                                                    &
436         ONLY:  nyn_x, nys_x, nzb_x, nzt_x
437
438     IMPLICIT NONE
439
440     REAL(wp) ::  f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) !<
441     REAL(wp) ::  f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x) !<
442
443
444     INTEGER(iwp) ::  i !<
445     INTEGER(iwp) ::  j !<
446     INTEGER(iwp) ::  k !<
447!
448!-- Rearrange indices of input array in order to make data to be send
449!-- by MPI contiguous
450    !$OMP  PARALLEL PRIVATE ( i, j, k )
451    !$OMP  DO
452#if __acc_fft_device
453     !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
454     !$ACC PRESENT(f_out, f_inv)
455#endif
456     DO  i = 0, nx
457         DO  k = nzb_x, nzt_x
458             DO  j = nys_x, nyn_x
459                 f_out(i,j,k) = f_inv(j,k,i)
460             ENDDO
461         ENDDO
462     ENDDO
463     !$OMP  END PARALLEL
464
465 END SUBROUTINE resort_for_yx
466
467
468!------------------------------------------------------------------------------!
469! Description:
470! ------------
471!> Transposition of input array (f_in) from y to x. For the input array, all
472!> elements along y reside on the same PE, while after transposition, all
473!> elements along x reside on the same PE.
474!------------------------------------------------------------------------------!
475 SUBROUTINE transpose_yx( f_in, f_inv )
476
477
478    USE cpulog,                                                                &
479        ONLY:  cpu_log, cpu_log_nowait, log_point_s
480
481    USE indices,                                                               &
482        ONLY:  nx, ny
483
484    USE kinds
485
486    USE pegrid
487
488    USE transpose_indices,                                                     &
489        ONLY:  nxl_y, nxr_y, nyn_x, nys_x, nzb_x, nzb_y, nzt_x, nzt_y
490
491    IMPLICIT NONE
492
493    INTEGER(iwp) ::  i  !<
494    INTEGER(iwp) ::  j  !<
495    INTEGER(iwp) ::  k  !<
496    INTEGER(iwp) ::  l  !<
497    INTEGER(iwp) ::  ys !<
498
499    REAL(wp) ::  f_in(0:ny,nxl_y:nxr_y,nzb_y:nzt_y)  !<
500    REAL(wp) ::  f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) !<
501
502    REAL(wp), DIMENSION(nyn_x-nys_x+1,nzb_y:nzt_y,nxl_y:nxr_y,0:pdims(2)-1) ::  work !<
503#if __acc_fft_device
504    !$ACC DECLARE CREATE(work)
505#endif
506
507
508    IF ( numprocs /= 1 )  THEN
509
510#if defined( __parallel )
511!
512!--    Reorder input array for transposition
513!$OMP  PARALLEL PRIVATE ( i, j, k, l, ys )
514!$OMP  DO
515       DO  l = 0, pdims(2) - 1
516          ys = 0 + l * ( nyn_x - nys_x + 1 )
517#if __acc_fft_device
518          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
519          !$ACC PRESENT(work, f_in)
520#endif
521          DO  i = nxl_y, nxr_y
522             DO  k = nzb_y, nzt_y
523                DO  j = ys, ys + nyn_x - nys_x
524                   work(j-ys+1,k,i,l) = f_in(j,i,k)
525                ENDDO
526             ENDDO
527          ENDDO
528       ENDDO
529!$OMP  END PARALLEL
530
531!
532!--    Transpose array
533       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start', cpu_log_nowait )
534
535#if __acc_fft_device
536#ifndef __cuda_aware_mpi
537       !$ACC UPDATE HOST(work)
538#else
539       !$ACC HOST_DATA USE_DEVICE(work, f_inv)
540#endif
541#endif
542
543       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
544       CALL MPI_ALLTOALL( work(1,nzb_y,nxl_y,0), sendrecvcount_xy, MPI_REAL, &
545                          f_inv(nys_x,nzb_x,0),  sendrecvcount_xy, MPI_REAL, &
546                          comm1dy, ierr )
547
548#if __acc_fft_device
549#ifndef __cuda_aware_mpi
550       !$ACC UPDATE DEVICE(f_inv)
551#else
552       !$ACC END HOST_DATA
553#endif
554#endif
555
556       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
557#endif
558
559    ELSE
560
561!
562!--    Reorder array f_in the same way as ALLTOALL did it
563!$OMP  PARALLEL PRIVATE ( i, j, k )
564!$OMP  DO
565#if __acc_fft_device
566       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
567       !$ACC PRESENT(f_inv, f_in)
568#endif
569       DO  i = nxl_y, nxr_y
570          DO  k = nzb_y, nzt_y
571             DO  j = 0, ny
572                f_inv(j,k,i) = f_in(j,i,k)
573             ENDDO
574          ENDDO
575       ENDDO
576!$OMP  END PARALLEL
577
578    ENDIF
579
580 END SUBROUTINE transpose_yx
581
582
583!------------------------------------------------------------------------------!
584! Description:
585! ------------
586!> Transposition of input array (f_in) from y to x. For the input array, all
587!> elements along y reside on the same PE, while after transposition, all
588!> elements along x reside on the same PE.
589!> This is a direct transposition for arrays with indices in regular order
590!> (k,j,i) (cf. transpose_yx).
591!------------------------------------------------------------------------------!
592 SUBROUTINE transpose_yxd( f_in, f_out )
593
594
595    USE cpulog,                                                                &
596        ONLY:  cpu_log, log_point_s
597
598    USE indices,                                                               &
599        ONLY:  nnx, nny, nnz, nx, nxl, nxr, nyn, nys, nz
600
601    USE kinds
602
603    USE pegrid
604
605    USE transpose_indices,                                                     &
606        ONLY:  nyn_x, nys_x, nzb_x, nzt_x
607
608    IMPLICIT NONE
609
610    INTEGER(iwp) ::  i  !<
611    INTEGER(iwp) ::  j  !<
612    INTEGER(iwp) ::  k  !<
613    INTEGER(iwp) ::  l  !<
614    INTEGER(iwp) ::  m  !<
615    INTEGER(iwp) ::  xs !<
616
617    REAL(wp) ::  f_in(1:nz,nys:nyn,nxl:nxr)          !<
618    REAL(wp) ::  f_inv(nxl:nxr,1:nz,nys:nyn)         !<
619    REAL(wp) ::  f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x) !<
620    REAL(wp) ::  work(nnx*nny*nnz)                   !<
621#if defined( __parallel )
622
623!
624!-- Rearrange indices of input array in order to make data to be send
625!-- by MPI contiguous
626    DO  k = 1, nz
627       DO  j = nys, nyn
628          DO  i = nxl, nxr
629             f_inv(i,k,j) = f_in(k,j,i)
630          ENDDO
631       ENDDO
632    ENDDO
633
634!
635!-- Transpose array
636    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
637    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
638    CALL MPI_ALLTOALL( f_inv(nxl,1,nys), sendrecvcount_xy, MPI_REAL, &
639                       work(1),          sendrecvcount_xy, MPI_REAL, &
640                       comm1dx, ierr )
641    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
642
643!
644!-- Reorder transposed array
645    m = 0
646    DO  l = 0, pdims(1) - 1
647       xs = 0 + l * nnx
648       DO  j = nys_x, nyn_x
649          DO  k = 1, nz
650             DO  i = xs, xs + nnx - 1
651                m = m + 1
652                f_out(i,j,k) = work(m)
653             ENDDO
654          ENDDO
655       ENDDO
656    ENDDO
657
658#endif
659
660 END SUBROUTINE transpose_yxd
661
662
663!------------------------------------------------------------------------------!
664! Description:
665! ------------
666!> Resorting data for the transposition from y to z. The transposition itself
667!> is carried out in transpose_yz
668!------------------------------------------------------------------------------!
669 SUBROUTINE resort_for_yz( f_in, f_inv )
670
671
672     USE indices,                                                              &
673         ONLY:  ny
674
675     USE kinds
676
677     USE transpose_indices,                                                    &
678         ONLY:  nxl_y, nxr_y, nzb_y, nzt_y
679
680     IMPLICIT NONE
681
682     REAL(wp) ::  f_in(0:ny,nxl_y:nxr_y,nzb_y:nzt_y)  !<
683     REAL(wp) ::  f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) !<
684
685     INTEGER(iwp) ::  i !<
686     INTEGER(iwp) ::  j !<
687     INTEGER(iwp) ::  k !<
688
689!
690!-- Rearrange indices of input array in order to make data to be send
691!-- by MPI contiguous
692    !$OMP  PARALLEL PRIVATE ( i, j, k )
693    !$OMP  DO
694#if __acc_fft_device
695     !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
696     !$ACC PRESENT(f_inv, f_in)
697#endif
698     DO  j = 0, ny
699         DO  k = nzb_y, nzt_y
700             DO  i = nxl_y, nxr_y
701                 f_inv(i,k,j) = f_in(j,i,k)
702             ENDDO
703         ENDDO
704     ENDDO
705     !$OMP  END PARALLEL
706
707 END SUBROUTINE resort_for_yz
708
709
710!------------------------------------------------------------------------------!
711! Description:
712! ------------
713!> Transposition of input array (f_in) from y to z. For the input array, all
714!> elements along y reside on the same PE, while after transposition, all
715!> elements along z reside on the same PE.
716!------------------------------------------------------------------------------!
717 SUBROUTINE transpose_yz( f_inv, f_out )
718
719
720    USE cpulog,                                                                &
721        ONLY:  cpu_log, cpu_log_nowait, log_point_s
722
723    USE indices,                                                               &
724        ONLY:  ny, nz
725
726    USE kinds
727
728    USE pegrid
729
730    USE transpose_indices,                                                     &
731        ONLY:  nxl_y, nxl_z, nxr_y, nxr_z, nyn_z, nys_z, nzb_y, nzt_y
732
733    IMPLICIT NONE
734
735    INTEGER(iwp) ::  i  !<
736    INTEGER(iwp) ::  j  !<
737    INTEGER(iwp) ::  k  !<
738    INTEGER(iwp) ::  l  !<
739    INTEGER(iwp) ::  zs !<
740
741    REAL(wp) ::  f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) !<
742    REAL(wp) ::  f_out(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !<
743
744    REAL(wp), DIMENSION(nxl_z:nxr_z,nzt_y-nzb_y+1,nys_z:nyn_z,0:pdims(1)-1) ::  work !<
745#if __acc_fft_device
746    !$ACC DECLARE CREATE(work)
747#endif
748
749
750!
751!-- If the PE grid is one-dimensional along y, only local reordering
752!-- of the data is necessary and no transposition has to be done.
753    IF ( pdims(1) == 1 )  THEN
754
755!$OMP  PARALLEL PRIVATE ( i, j, k )
756!$OMP  DO
757#if __acc_fft_device
758       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
759       !$ACC PRESENT(f_out, f_inv)
760#endif
761       DO  j = 0, ny
762          DO  k = nzb_y, nzt_y
763             DO  i = nxl_y, nxr_y
764                f_out(i,j,k) = f_inv(i,k,j)
765             ENDDO
766          ENDDO
767       ENDDO
768!$OMP  END PARALLEL
769
770    ELSE
771
772#if defined( __parallel )
773!
774!--    Transpose array
775       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start', cpu_log_nowait )
776
777#if __acc_fft_device
778#ifndef __cuda_aware_mpi
779       !$ACC UPDATE HOST(f_inv)
780#else
781       !$ACC HOST_DATA USE_DEVICE(work, f_inv)
782#endif
783#endif
784
785       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
786       CALL MPI_ALLTOALL( f_inv(nxl_y,nzb_y,0),  sendrecvcount_yz, MPI_REAL, &
787                          work(nxl_z,1,nys_z,0), sendrecvcount_yz, MPI_REAL, &
788                          comm1dx, ierr )
789
790#if __acc_fft_device
791#ifndef __cuda_aware_mpi
792       !$ACC UPDATE DEVICE(work)
793#else
794       !$ACC END HOST_DATA
795#endif
796#endif
797
798       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
799
800!
801!--    Reorder transposed array
802!$OMP  PARALLEL PRIVATE ( i, j, k, l, zs )
803!$OMP  DO
804       DO  l = 0, pdims(1) - 1
805          zs = 1 + l * ( nzt_y - nzb_y + 1 )
806#if __acc_fft_device
807          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
808          !$ACC PRESENT(f_out, work)
809#endif
810          DO  j = nys_z, nyn_z
811             DO  k = zs, zs + nzt_y - nzb_y
812                DO  i = nxl_z, nxr_z
813                   f_out(i,j,k) = work(i,k-zs+1,j,l)
814                ENDDO
815             ENDDO
816          ENDDO
817       ENDDO
818!$OMP  END PARALLEL
819#endif
820
821   ENDIF
822
823 END SUBROUTINE transpose_yz
824
825
826!------------------------------------------------------------------------------!
827! Description:
828! ------------
829!> Resorting data for the transposition from z to x. The transposition itself
830!> is carried out in transpose_zx
831!------------------------------------------------------------------------------!
832 SUBROUTINE resort_for_zx( f_in, f_inv )
833
834
835     USE indices,                                                              &
836         ONLY:  nxl, nxr, nyn, nys, nz
837
838     USE kinds
839
840     IMPLICIT NONE
841
842     REAL(wp) ::  f_in(1:nz,nys:nyn,nxl:nxr)  !<
843     REAL(wp) ::  f_inv(nys:nyn,nxl:nxr,1:nz) !<
844
845     INTEGER(iwp) ::  i !<
846     INTEGER(iwp) ::  j !<
847     INTEGER(iwp) ::  k !<
848
849!
850!-- Rearrange indices of input array in order to make data to be send
851!-- by MPI contiguous
852    !$OMP  PARALLEL PRIVATE ( i, j, k )
853    !$OMP  DO
854#if __acc_fft_device
855    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
856    !$ACC PRESENT(f_in, f_inv)
857#endif
858     DO  k = 1,nz
859         DO  i = nxl, nxr
860             DO  j = nys, nyn
861                 f_inv(j,i,k) = f_in(k,j,i)
862             ENDDO
863         ENDDO
864     ENDDO
865     !$OMP  END PARALLEL
866
867 END SUBROUTINE resort_for_zx
868
869
870!------------------------------------------------------------------------------!
871! Description:
872! ------------
873!> Transposition of input array (f_in) from z to x. For the input array, all
874!> elements along z reside on the same PE, while after transposition, all
875!> elements along x reside on the same PE.
876!------------------------------------------------------------------------------!
877 SUBROUTINE transpose_zx( f_inv, f_out )
878
879
880    USE cpulog,                                                                &
881        ONLY:  cpu_log, cpu_log_nowait, log_point_s
882
883    USE indices,                                                               &
884        ONLY:  nnx, nx, nxl, nxr, nyn, nys, nz
885
886    USE kinds
887
888    USE pegrid
889
890    USE transpose_indices,                                                     &
891        ONLY:  nyn_x, nys_x, nzb_x, nzt_x
892
893    IMPLICIT NONE
894
895    INTEGER(iwp) ::  i  !<
896    INTEGER(iwp) ::  j  !<
897    INTEGER(iwp) ::  k  !<
898    INTEGER(iwp) ::  l  !<
899    INTEGER(iwp) ::  xs !<
900
901    REAL(wp) ::  f_inv(nys:nyn,nxl:nxr,1:nz)         !<
902    REAL(wp) ::  f_out(0:nx,nys_x:nyn_x,nzb_x:nzt_x) !<
903
904    REAL(wp), DIMENSION(nys_x:nyn_x,nnx,nzb_x:nzt_x,0:pdims(1)-1) ::  work !<
905#if __acc_fft_device
906    !$ACC DECLARE CREATE(work)
907#endif
908
909
910!
911!-- If the PE grid is one-dimensional along y, only local reordering
912!-- of the data is necessary and no transposition has to be done.
913    IF ( pdims(1) == 1 )  THEN
914
915!$OMP  PARALLEL PRIVATE ( i, j, k )
916!$OMP  DO
917#if __acc_fft_device
918       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
919       !$ACC PRESENT(f_out, f_inv)
920#endif
921       DO  k = 1, nz
922          DO  i = nxl, nxr
923             DO  j = nys, nyn
924                f_out(i,j,k) = f_inv(j,i,k)
925             ENDDO
926          ENDDO
927       ENDDO
928!$OMP  END PARALLEL
929
930    ELSE
931
932#if defined( __parallel )
933!
934!--    Transpose array
935       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start', cpu_log_nowait )
936
937#if __acc_fft_device
938#ifndef __cuda_aware_mpi
939       !$ACC UPDATE HOST(f_inv)
940#else
941       !$ACC HOST_DATA USE_DEVICE(work, f_inv)
942#endif
943#endif
944
945       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
946       CALL MPI_ALLTOALL( f_inv(nys,nxl,1),      sendrecvcount_zx, MPI_REAL, &
947                          work(nys_x,1,nzb_x,0), sendrecvcount_zx, MPI_REAL, &
948                          comm1dx, ierr )
949
950#if __acc_fft_device
951#ifndef __cuda_aware_mpi
952       !$ACC UPDATE DEVICE(work)
953#else
954       !$ACC END HOST_DATA
955#endif
956#endif
957
958       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
959
960!
961!--    Reorder transposed array
962!$OMP  PARALLEL PRIVATE ( i, j, k, l, xs )
963!$OMP  DO
964       DO  l = 0, pdims(1) - 1
965          xs = 0 + l * nnx
966#if __acc_fft_device
967          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
968          !$ACC PRESENT(f_out, work)
969#endif
970          DO  k = nzb_x, nzt_x
971             DO  i = xs, xs + nnx - 1
972                DO  j = nys_x, nyn_x
973                   f_out(i,j,k) = work(j,i-xs+1,k,l)
974                ENDDO
975             ENDDO
976          ENDDO
977       ENDDO
978!$OMP  END PARALLEL
979#endif
980
981    ENDIF
982
983 END SUBROUTINE transpose_zx
984
985
986!------------------------------------------------------------------------------!
987! Description:
988! ------------
989!> Resorting data after the transposition from z to y. The transposition itself
990!> is carried out in transpose_zy
991!------------------------------------------------------------------------------!
992 SUBROUTINE resort_for_zy( f_inv, f_out )
993
994
995     USE indices,                                                              &
996         ONLY:  ny
997
998     USE kinds
999
1000     USE transpose_indices,                                                    &
1001         ONLY:  nxl_y, nxr_y, nzb_y, nzt_y
1002
1003     IMPLICIT NONE
1004
1005     REAL(wp) ::  f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) !<
1006     REAL(wp) ::  f_out(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) !<
1007
1008
1009     INTEGER(iwp) ::  i !<
1010     INTEGER(iwp) ::  j !<
1011     INTEGER(iwp) ::  k !<
1012
1013!
1014!-- Rearrange indices of input array in order to make data to be send
1015!-- by MPI contiguous
1016    !$OMP  PARALLEL PRIVATE ( i, j, k )
1017    !$OMP  DO
1018#if __acc_fft_device
1019    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
1020    !$ACC PRESENT(f_out, f_inv)
1021#endif
1022     DO  k = nzb_y, nzt_y
1023         DO  j = 0, ny
1024             DO  i = nxl_y, nxr_y
1025                 f_out(j,i,k) = f_inv(i,k,j)
1026             ENDDO
1027         ENDDO
1028     ENDDO
1029     !$OMP  END PARALLEL
1030
1031 END SUBROUTINE resort_for_zy
1032
1033
1034!------------------------------------------------------------------------------!
1035! Description:cpu_log_nowait
1036! ------------
1037!> Transposition of input array (f_in) from z to y. For the input array, all
1038!> elements along z reside on the same PE, while after transposition, all
1039!> elements along y reside on the same PE.
1040!------------------------------------------------------------------------------!
1041 SUBROUTINE transpose_zy( f_in, f_inv )
1042
1043
1044    USE cpulog,                                                                &
1045        ONLY:  cpu_log, cpu_log_nowait, log_point_s
1046
1047    USE indices,                                                               &
1048        ONLY:  ny, nz
1049
1050    USE kinds
1051
1052    USE pegrid
1053
1054    USE transpose_indices,                                                     &
1055        ONLY:  nxl_y, nxl_z, nxr_y, nxr_z, nyn_z, nys_z, nzb_y, nzt_y
1056
1057    IMPLICIT NONE
1058
1059    INTEGER(iwp) ::  i  !<
1060    INTEGER(iwp) ::  j  !<
1061    INTEGER(iwp) ::  k  !<
1062    INTEGER(iwp) ::  l  !<
1063    INTEGER(iwp) ::  zs !<
1064
1065    REAL(wp) ::  f_in(nxl_z:nxr_z,nys_z:nyn_z,1:nz)  !<
1066    REAL(wp) ::  f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) !<
1067
1068    REAL(wp), DIMENSION(nxl_z:nxr_z,nzt_y-nzb_y+1,nys_z:nyn_z,0:pdims(1)-1) ::  work !<
1069#if __acc_fft_device
1070    !$ACC DECLARE CREATE(work)
1071#endif
1072
1073!
1074!-- If the PE grid is one-dimensional along y, the array has only to be
1075!-- reordered locally and therefore no transposition has to be done.
1076    IF ( pdims(1) /= 1 )  THEN
1077
1078#if defined( __parallel )
1079!
1080!--    Reorder input array for transposition
1081!$OMP  PARALLEL PRIVATE ( i, j, k, l, zs )
1082!$OMP  DO
1083       DO  l = 0, pdims(1) - 1
1084          zs = 1 + l * ( nzt_y - nzb_y + 1 )
1085#if __acc_fft_device
1086          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
1087          !$ACC PRESENT(work, f_in)
1088#endif
1089          DO  j = nys_z, nyn_z
1090             DO  k = zs, zs + nzt_y - nzb_y
1091                DO  i = nxl_z, nxr_z
1092                   work(i,k-zs+1,j,l) = f_in(i,j,k)
1093                ENDDO
1094             ENDDO
1095          ENDDO
1096       ENDDO
1097!$OMP  END PARALLEL
1098
1099!
1100!--    Transpose array
1101       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start', cpu_log_nowait )
1102
1103#if __acc_fft_device
1104#ifndef __cuda_aware_mpi
1105       !$ACC UPDATE HOST(work)
1106#else
1107       !$ACC HOST_DATA USE_DEVICE(work, f_inv)
1108#endif
1109#endif
1110
1111       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1112       CALL MPI_ALLTOALL( work(nxl_z,1,nys_z,0), sendrecvcount_yz, MPI_REAL, &
1113                          f_inv(nxl_y,nzb_y,0),  sendrecvcount_yz, MPI_REAL, &
1114                          comm1dx, ierr )
1115
1116#if __acc_fft_device
1117#ifndef __cuda_aware_mpi
1118       !$ACC UPDATE DEVICE(f_inv)
1119#else
1120       !$ACC END HOST_DATA
1121#endif
1122#endif
1123
1124       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1125#endif
1126
1127    ELSE
1128!
1129!--    Reorder the array in the same way like ALLTOALL did it
1130!$OMP  PARALLEL PRIVATE ( i, j, k )
1131!$OMP  DO
1132#if __acc_fft_device
1133       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
1134       !$ACC PRESENT(f_inv, f_in)
1135#endif
1136       DO  k = nzb_y, nzt_y
1137          DO  j = 0, ny
1138             DO  i = nxl_y, nxr_y
1139                f_inv(i,k,j) = f_in(i,j,k)
1140             ENDDO
1141          ENDDO
1142       ENDDO
1143!$OMP  END PARALLEL
1144
1145    ENDIF
1146
1147 END SUBROUTINE transpose_zy
1148
1149
1150!------------------------------------------------------------------------------!
1151! Description:
1152! ------------
1153!> Transposition of input array (f_in) from z to y. For the input array, all
1154!> elements along z reside on the same PE, while after transposition, all
1155!> elements along y reside on the same PE.
1156!> This is a direct transposition for arrays with indices in regular order
1157!> (k,j,i) (cf. transpose_zy).
1158!------------------------------------------------------------------------------!
1159 SUBROUTINE transpose_zyd( f_in, f_out )
1160
1161
1162    USE cpulog,                                                                &
1163        ONLY:  cpu_log, log_point_s
1164
1165    USE indices,                                                               &
1166        ONLY:  nnx, nny, nnz, nxl, nxr, nyn, nys, ny, nz
1167
1168    USE kinds
1169
1170    USE pegrid
1171
1172    USE transpose_indices,                                                     &
1173        ONLY:  nxl_yd, nxr_yd, nzb_yd, nzt_yd
1174
1175    IMPLICIT NONE
1176
1177    INTEGER(iwp) ::  i  !<
1178    INTEGER(iwp) ::  j  !<
1179    INTEGER(iwp) ::  k  !<
1180    INTEGER(iwp) ::  l  !<
1181    INTEGER(iwp) ::  m  !<
1182    INTEGER(iwp) ::  ys !<
1183
1184    REAL(wp) ::  f_in(1:nz,nys:nyn,nxl:nxr)              !<
1185    REAL(wp) ::  f_inv(nys:nyn,nxl:nxr,1:nz)             !<
1186    REAL(wp) ::  f_out(0:ny,nxl_yd:nxr_yd,nzb_yd:nzt_yd) !<
1187    REAL(wp) ::  work(nnx*nny*nnz)                       !<
1188
1189#if defined( __parallel )
1190
1191!
1192!-- Rearrange indices of input array in order to make data to be send
1193!-- by MPI contiguous
1194    DO  i = nxl, nxr
1195       DO  j = nys, nyn
1196          DO  k = 1, nz
1197             f_inv(j,i,k) = f_in(k,j,i)
1198          ENDDO
1199       ENDDO
1200    ENDDO
1201
1202!
1203!-- Move data to different array, because memory location of work1 is
1204!-- needed further below (work1 = work2).
1205!-- If the PE grid is one-dimensional along x, only local reordering
1206!-- of the data is necessary and no transposition has to be done.
1207    IF ( pdims(2) == 1 )  THEN
1208       DO  k = 1, nz
1209          DO  i = nxl, nxr
1210             DO  j = nys, nyn
1211                f_out(j,i,k) = f_inv(j,i,k)
1212             ENDDO
1213          ENDDO
1214       ENDDO
1215       RETURN
1216    ENDIF
1217
1218!
1219!-- Transpose array
1220    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1221    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1222    CALL MPI_ALLTOALL( f_inv(nys,nxl,1), sendrecvcount_zyd, MPI_REAL, &
1223                       work(1),          sendrecvcount_zyd, MPI_REAL, &
1224                       comm1dy, ierr )
1225    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1226
1227!
1228!-- Reorder transposed array
1229    m = 0
1230    DO  l = 0, pdims(2) - 1
1231       ys = 0 + l * nny
1232       DO  k = nzb_yd, nzt_yd
1233          DO  i = nxl_yd, nxr_yd
1234             DO  j = ys, ys + nny - 1
1235                m = m + 1
1236                f_out(j,i,k) = work(m)
1237             ENDDO
1238          ENDDO
1239       ENDDO
1240    ENDDO
1241
1242#endif
1243
1244 END SUBROUTINE transpose_zyd
Note: See TracBrowser for help on using the repository browser.