source: palm/trunk/SOURCE/poisfft_hybrid.f90 @ 1207

Last change on this file since 1207 was 1112, checked in by raasch, 12 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 33.4 KB
Line 
1 MODULE poisfft_hybrid_mod
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: poisfft_hybrid.f90 1112 2013-03-09 00:34:37Z witha $
27!
28! 1111 2013-03-08 23:54:10Z raasch
29! poisfft_hybrid_ini is now called internally from poisfft_hybrid,
30! ibc_p_b = 2 removed
31!
32! 1106 2013-03-04 05:31:38Z raasch
33! calls of fft_x, fft_y replaced by fft_x_1d, fft_y_1d
34!
35! 1036 2012-10-22 13:43:42Z raasch
36! code put under GPL (PALM 3.9)
37!
38! 1013 2012-09-21 07:03:55Z raasch
39! FLOAT type conversion replaced by REAL
40!
41! 809 2012-01-30 13:32:58Z maronga
42! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
43!
44! 807 2012-01-25 11:53:51Z maronga
45! New cpp directive "__check" implemented which is used by check_namelist_files
46! (most of the code is unneeded by check_namelist_files).
47!
48! 667 2010-12-23 12:06:00Z suehring/gryschka
49! ddzu replaced by ddzu_pres due to changes in zu(0)
50!
51! 415 2009-12-15 10:26:23Z raasch
52! Dimension of array stat in cascade change to prevent type problems with___
53! mpi2 libraries
54!
55! 274 2009-03-26 15:11:21Z heinze
56! Output of messages replaced by message handling routine.
57!
58! Feb. 2007
59! RCS Log replace by Id keyword, revision history cleaned up
60!
61! Revision 1.11  2004/04/30 12:43:14  raasch
62! Renaming of fft routines, additional argument in calls of fft_y_m
63!
64! Revision 1.2  2002/12/19 16:08:31  raasch
65! Preprocessor directive KKMP introduced (OMP does NOT work),
66! array tri will be a shared array in OpenMP loop, to get better cache
67! utilization, the i index  (x-direction) will be executed in stride
68! "istride" as outer loop and in a shorter inner loop,
69! overlapping of computation and communication realized by new routine
70! poisfft_hybrid_nodes, name of old routine poisfft_hybrid changed to
71! poisfft_hybrid_omp, STOP statement replaced by call of subroutine local_stop
72!
73!
74! Description:
75! ------------
76! Solution of the Poisson equation with a 2D spectral method. 
77! Hybrid version for parallel computers using a 1D domain decomposition,
78! realized with MPI, along x and parallelization with OPEN-MP along y
79! (routine poisfft_hybrid_omp). In a second version (poisfft_hybrid_nodes),
80! optimization is realized by overlapping of computation and communication
81! and by simultaneously executing as many communication calls as switches
82! per logical partition (LPAR) are available. This version comes into
83! effect if more than one node is used and if the environment variable
84! tasks_per_node is set in a way that it can be devided by switch_per_lpar
85! without any rest.
86!
87! WARNING: In case of OpenMP, there are problems with allocating large
88! arrays in parallel regions.
89!
90! Copyright   Klaus Ketelsen / Siegfried Raasch   May 2002
91!------------------------------------------------------------------------------!
92
93    USE fft_xy
94    USE indices
95    USE pegrid
96
97    IMPLICIT NONE
98
99    INTEGER, PARAMETER ::  switch_per_lpar = 2
100
101    INTEGER, SAVE ::  nxl_a, nxr_a, &      ! total   x dimension
102                      nxl_p, nxr_p, &      ! partial x dimension
103                      nys_a, nyn_a, &      ! total   y dimension
104                      nys_p, nyn_p, &      ! partial y dimension
105
106                      npe_s,        &      ! total number of PEs for solver
107                      nwords,       &      ! number of points to be exchanged
108                                           ! with MPI_ALLTOALL
109                      n_omp_threads        ! number of OpenMP threads
110
111!
112!-- Variables for multi node version (cluster version) using routine
113!-- poisfft_hybrid_nodes
114    INTEGER, SAVE :: comm_nodes,          & ! communicater nodes
115                     comm_node_all,       & ! communicater all PEs node version
116                     comm_tasks,          & ! communicater tasks
117                     me, me_node, me_task,& ! identity of this PE
118                     nodes,               & ! number of nodes
119                     tasks_per_logical_node = -1    ! default no cluster
120
121    LOGICAL, SAVE ::  poisfft_initialized = .FALSE.
122
123
124    PRIVATE
125
126
127#if ! defined ( __check )
128    PUBLIC poisfft_hybrid, poisfft_hybrid_ini
129
130
131!
132!-- Public interfaces
133    INTERFACE poisfft_hybrid_ini
134       MODULE PROCEDURE poisfft_hybrid_ini
135    END INTERFACE poisfft_hybrid_ini
136
137    INTERFACE poisfft_hybrid
138       MODULE PROCEDURE poisfft_hybrid
139    END INTERFACE poisfft_hybrid
140
141!
142!-- Private interfaces
143    INTERFACE poisfft_hybrid_omp
144       MODULE PROCEDURE poisfft_hybrid_omp
145    END INTERFACE poisfft_hybrid_omp
146
147    INTERFACE poisfft_hybrid_omp_vec
148       MODULE PROCEDURE poisfft_hybrid_omp_vec
149    END INTERFACE poisfft_hybrid_omp_vec
150
151    INTERFACE poisfft_hybrid_nodes
152       MODULE PROCEDURE poisfft_hybrid_nodes
153    END INTERFACE poisfft_hybrid_nodes
154
155    INTERFACE tridia_hybrid
156       MODULE PROCEDURE tridia_hybrid
157    END INTERFACE tridia_hybrid
158
159    INTERFACE cascade
160       MODULE PROCEDURE cascade
161    END INTERFACE cascade
162#else
163    PUBLIC poisfft_hybrid_ini
164
165!
166!-- Public interfaces
167    INTERFACE poisfft_hybrid_ini
168       MODULE PROCEDURE poisfft_hybrid_ini
169    END INTERFACE poisfft_hybrid_ini
170#endif
171
172 CONTAINS
173 
174
175    SUBROUTINE poisfft_hybrid_ini
176
177       USE control_parameters
178       USE pegrid
179
180       IMPLICIT NONE
181
182       CHARACTER(LEN=8)      ::  cdummy
183       INTEGER               ::  idummy, istat
184       INTEGER, DIMENSION(2) ::  coords, dims
185
186       LOGICAL, DIMENSION(2) ::  period = .false., re_dims
187
188
189!
190!--    Set the internal index values for the hybrid solver
191#if defined( __parallel )
192       npe_s = pdims(1)
193#else
194       npe_s = 1
195#endif
196       nxl_a = 0
197       nxr_a = nx
198       nxl_p = 0
199       nxr_p = ( ( nx+1 ) / npe_s ) - 1
200       nys_a = nys 
201       nyn_a = nyn
202       nys_p = 0
203       nyn_p = ( ( ny+1 ) / npe_s ) - 1
204
205       nwords = ( nxr_p-nxl_p+1 ) * nz * ( nyn_p-nys_p+1 )
206
207#if defined( __KKMP ) && ! defined ( __check )
208       CALL LOCAL_GETENV( 'OMP_NUM_THREADS', 15, cdummy, idummy )
209       READ ( cdummy, '(I8)' )  n_omp_threads
210       IF ( n_omp_threads > 1 )  THEN
211          WRITE( message_string, * ) 'Number of OpenMP threads = ', &
212                                     n_omp_threads
213          CALL message( 'poisfft_hybrid_ini', 'PA0280', 0, 0, 0, 6, 0 ) 
214       ENDIF
215#else
216       n_omp_threads = 1
217#endif
218!
219!--    Initialize the one-dimensional FFT routines
220       CALL fft_init
221
222!
223!--    Setup for multi node version (poisfft_hybrid_nodes)
224       IF ( n_omp_threads == 1  .AND.  &
225            ( host(1:4) == 'ibmh' .OR. host(1:4) == 'ibmb' ) )  THEN
226
227          IF ( tasks_per_node /= -9999 )  THEN
228!
229!--          Multi node version requires that the available number of
230!--          switches per logical partition must be an integral divisor
231!--          of the chosen number of tasks per node
232             IF ( MOD( tasks_per_node, switch_per_lpar ) == 0 )  THEN
233!
234!--             Set the switch which decides about usage of the multi node
235!--             version
236                IF ( tasks_per_node / switch_per_lpar > 1  .AND. &
237                     numprocs > tasks_per_node )  THEN
238                   tasks_per_logical_node = tasks_per_node / switch_per_lpar
239                ENDIF
240
241                IF ( tasks_per_logical_node > -1 )  THEN
242
243                   WRITE( message_string, * ) 'running optimized ',         &
244                                              'multinode version',          &
245                                              '&switch_per_lpar        = ', &
246                                              switch_per_lpar,              &
247                                              '&tasks_per_lpar         = ', &
248                                              tasks_per_node,               &
249                                              'tasks_per_logical_node = ',  &
250                                              tasks_per_logical_node
251                   CALL message( 'poisfft_hybrid_ini', 'PA0281', 0, 0, 0, 6, 0 )
252
253                ENDIF
254
255             ENDIF
256          ENDIF
257       ENDIF
258
259!
260!--    Determine sub-topologies for multi node version
261       IF ( tasks_per_logical_node >= 2 )  THEN
262
263#if defined( __parallel ) && ! defined ( __check )
264          nodes   = ( numprocs + tasks_per_logical_node - 1 ) / &
265                    tasks_per_logical_node
266          dims(1) = nodes
267          dims(2) = tasks_per_logical_node
268
269          CALL MPI_CART_CREATE( comm2d, 2, dims, period, .FALSE., &
270                                comm_node_all, istat )
271          CALL MPI_COMM_RANK( comm_node_all, me, istat )
272
273          re_dims(1) = .TRUE.
274          re_dims(2) = .FALSE.
275          CALL MPI_CART_SUB( comm_node_all, re_dims, comm_nodes, istat )
276          CALL MPI_COMM_RANK( comm_nodes, me_node, istat )
277
278          re_dims(1) = .FALSE.
279          re_dims(2) = .TRUE.
280          CALL MPI_CART_SUB( comm_node_all, re_dims, comm_tasks, istat )
281          CALL MPI_COMM_RANK( comm_tasks, me_task, istat )
282
283!          write(0,*) 'who am i',myid,me,me_node,me_task,nodes,&
284!                     tasks_per_logical_node
285#elif ! defined( __parallel )
286          message_string = 'parallel environment (MPI) required'
287          CALL message( 'poisfft_hybrid_ini', 'PA0282', 1, 2, 0, 6, 0 )
288#endif
289       ENDIF
290
291       poisfft_initialized = .TRUE.
292
293    END SUBROUTINE poisfft_hybrid_ini
294
295#if ! defined ( __check )
296    SUBROUTINE poisfft_hybrid( ar )
297
298       USE control_parameters
299       USE interfaces
300
301       IMPLICIT NONE
302
303       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar
304
305       IF ( .NOT. poisfft_initialized )  CALL poisfft_hybrid_ini
306
307       IF ( host(1:3) == 'nec' )  THEN
308          CALL poisfft_hybrid_omp_vec( ar )
309       ELSE
310          IF ( tasks_per_logical_node == -1 )  THEN
311             CALL poisfft_hybrid_omp( ar )
312          ELSE
313             CALL poisfft_hybrid_nodes( ar )
314          ENDIF
315       ENDIF
316
317    END SUBROUTINE poisfft_hybrid
318
319
320    SUBROUTINE poisfft_hybrid_omp ( ar )
321
322       USE cpulog
323       USE interfaces
324
325       IMPLICIT NONE
326
327       INTEGER, PARAMETER ::  istride = 4  ! stride of i loop
328       INTEGER ::  i, ii, ir, iei, iouter, istat, j, jj, k, m, n, jthread
329
330       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar
331
332       REAL, DIMENSION(0:nx)                 ::  fftx_ar
333       REAL, DIMENSION(0:ny,istride)         ::  ffty_ar
334 
335       REAL, DIMENSION(0:nx,nz)              ::  tri_ar
336
337       REAL, DIMENSION(nxl_p:nxr_p,nz,nys_p:nyn_p,npe_s) ::  work1, work2
338#if defined( __KKMP )
339       INTEGER ::  omp_get_thread_num
340       REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  tri
341       ALLOCATE( tri(5,0:nx,0:nz-1,n_omp_threads ) )
342#else
343       REAL, DIMENSION(5,0:nx,0:nz-1,1)      ::  tri
344#endif
345
346
347       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_omp', 'start' )
348
349       CALL cpu_log( log_point_s(7), 'fft_y_1d', 'start' )
350
351!$OMP  PARALLEL PRIVATE (i,iouter,ii,ir,iei,j,k,m,n,ffty_ar)
352!$OMP  DO
353!
354!--    Store grid points to be transformed on a 1d-array, do the fft
355!--    and sample the results on a 4d-array
356       DO  iouter = nxl_p, nxr_p, istride     ! stride loop, better cache
357          iei = MIN( iouter+istride-1, nxr_p )
358          DO  k = 1, nz
359
360             DO  i = iouter, iei
361                ii = nxl + i
362                ir = i - iouter + 1
363
364                DO  j = nys_a, nyn_a
365                   ffty_ar(j,ir) = ar(k,j,ii)
366                ENDDO
367   
368                CALL fft_y_1d( ffty_ar(:,ir), 'forward' )
369             ENDDO
370
371             m = nys_a
372             DO  n = 1, npe_s
373                DO  j = nys_p, nyn_p
374                   DO  i = iouter, iei
375                      ir = i - iouter + 1
376                      work1(i,k,j,n) = ffty_ar(m,ir)
377                   ENDDO
378                   m = m+1
379                ENDDO
380             ENDDO
381   
382          ENDDO
383       ENDDO
384!$OMP  END PARALLEL
385
386       CALL cpu_log( log_point_s(7), 'fft_y_1d', 'pause' )
387
388#if defined( __parallel )
389       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
390
391       CALL MPI_ALLTOALL( work1(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
392                          work2(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
393                          comm2d, istat )
394
395       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
396#else
397       work2 = work1
398#endif
399
400       CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'start' )
401
402#if defined( __KKMP )
403!$OMP  PARALLEL PRIVATE (i,j,jj,k,m,n,fftx_ar,tri_ar,jthread)
404!$OMP  DO
405       DO  j = nys_p, nyn_p
406          jthread = omp_get_thread_num() + 1
407#else
408       DO  j = nys_p, nyn_p
409          jthread = 1
410#endif
411          DO  k = 1, nz
412
413             m = nxl_a
414             DO  n = 1, npe_s
415                DO  i = nxl_p, nxr_p
416                   fftx_ar(m) = work2(i,k,j,n)
417                   m = m+1
418                ENDDO
419             ENDDO
420
421             CALL fft_x_1d( fftx_ar, 'forward' )
422
423             DO  i = nxl_a, nxr_a
424                tri_ar(i,k) = fftx_ar(i)
425             ENDDO
426
427          ENDDO
428 
429          jj = myid * (nyn_p-nys_p+1) + j
430          CALL tridia_hybrid( jj, tri_ar, tri(:,:,:,jthread))
431
432          DO  k = 1, nz
433             DO  i = nxl_a, nxr_a
434                fftx_ar(i) = tri_ar (i,k)
435             ENDDO
436
437             CALL fft_x_1d( fftx_ar, 'backward' )
438
439             m = nxl_a
440             DO  n = 1, npe_s
441                DO  i = nxl_p, nxr_p
442                   work2(i,k,j,n) = fftx_ar(m)
443                   m = m+1
444                ENDDO
445             ENDDO
446
447          ENDDO
448       ENDDO
449#if defined( __KKMP )
450!$OMP  END PARALLEL
451#endif
452
453       CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'stop' )
454
455#if defined( __parallel )
456       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' ) 
457       nwords = (nxr_p-nxl_p+1) * nz * (nyn_p-nys_p+1)
458
459       CALL MPI_ALLTOALL( work2(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
460                          work1(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
461                          comm2d, istat )
462
463       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
464#else
465       work1 = work2
466#endif
467
468       CALL cpu_log( log_point_s(7), 'fft_y_1d', 'continue' )
469
470!$OMP  PARALLEL PRIVATE (i,iouter,ii,ir,iei,j,k,m,n,ffty_ar)
471!$OMP  DO
472       DO  iouter = nxl_p, nxr_p, istride
473          iei = MIN( iouter+istride-1, nxr_p )
474          DO  k = 1, nz
475
476             m = nys_a
477             DO  n = 1, npe_s
478                DO  j = nys_p, nyn_p
479                   DO  i = iouter, iei
480                      ir = i - iouter + 1
481                      ffty_ar(m,ir) = work1 (i,k,j,n)
482                   ENDDO
483                   m = m+1
484                ENDDO
485             ENDDO
486
487             DO  i = iouter, iei
488                ii = nxl + i
489                ir = i - iouter + 1
490                CALL fft_y_1d( ffty_ar(:,ir), 'backward' )
491
492                DO  j = nys_a, nyn_a
493                   ar(k,j,ii) = ffty_ar(j,ir)
494                ENDDO
495             ENDDO
496
497          ENDDO
498       ENDDO
499!$OMP  END PARALLEL
500
501       CALL cpu_log( log_point_s(7), 'fft_y_1d', 'stop' )
502
503       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_omp', 'stop' )
504
505#if defined( __KKMP )
506       DEALLOCATE( tri )
507#endif
508
509    END SUBROUTINE poisfft_hybrid_omp
510
511
512    SUBROUTINE poisfft_hybrid_omp_vec ( ar )
513
514       USE cpulog
515       USE interfaces
516
517       IMPLICIT NONE
518
519       INTEGER, PARAMETER ::  istride = 4  ! stride of i loop
520       INTEGER ::  i, ii, ir, iei, iouter, istat, j, jj, k, m, n, jthread
521
522       REAL, DIMENSION(0:nx,nz)               ::  tri_ar
523
524       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr)  ::  ar
525
526       REAL, DIMENSION(0:ny+3,nz,nxl_p:nxr_p) ::  ffty_ar3
527       REAL, DIMENSION(0:nx+3,nz,nys_p:nyn_p) ::  fftx_ar3
528 
529       REAL, DIMENSION(nxl_p:nxr_p,nz,nys_p:nyn_p,npe_s) ::  work1, work2
530#if defined( __KKMP )
531       INTEGER ::  omp_get_thread_num
532       REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  tri
533       ALLOCATE( tri(5,0:nx,0:nz-1,n_omp_threads ) )
534#else
535       REAL, DIMENSION(5,0:nx,0:nz-1,1)      ::  tri
536#endif
537
538
539       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_vec', 'start' )
540
541       CALL cpu_log( log_point_s(7), 'fft_y_m', 'start' )
542
543!$OMP  PARALLEL PRIVATE (i,j,k,m,n)
544!$OMP  DO
545!
546!--    Store grid points to be transformed on a 1d-array, do the fft
547!--    and sample the results on a 4d-array
548       DO  i = nxl_p, nxr_p
549
550          DO  j = nys_a, nyn_a
551             DO  k = 1, nz
552                ffty_ar3(j,k,i) = ar(k,j,i+nxl)
553             ENDDO
554          ENDDO
555
556          CALL fft_y_m( ffty_ar3(:,:,i), ny+3, 'forward' )
557       ENDDO
558
559!$OMP  DO
560       DO  k = 1, nz
561          m = nys_a
562          DO  n = 1, npe_s
563             DO  j = nys_p, nyn_p
564                 DO  i = nxl_p, nxr_p
565                     work1(i,k,j,n) = ffty_ar3(m,k,i)
566                 ENDDO
567                 m = m+1
568             ENDDO
569          ENDDO
570       ENDDO
571!$OMP  END PARALLEL
572
573       CALL cpu_log( log_point_s(7), 'fft_y_m', 'pause' )
574
575#if defined( __parallel )
576       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
577       CALL MPI_ALLTOALL( work1(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
578                          work2(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
579                          comm2d, istat )
580       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
581#else
582       work2 = work1
583#endif
584
585       CALL cpu_log( log_point_s(33), 'fft_x_m + tridia', 'start' )
586
587#if defined( __KKMP )
588!$OMP  PARALLEL PRIVATE (i,j,jj,k,m,n,tri_ar,jthread)
589!$OMP  DO
590       DO  j = nys_p, nyn_p
591          jthread = omp_get_thread_num() + 1
592#else
593       DO  j = nys_p, nyn_p
594          jthread = 1
595#endif
596          DO  k = 1, nz
597
598             m = nxl_a
599             DO  n = 1, npe_s
600                DO  i = nxl_p, nxr_p
601                   fftx_ar3(m,k,j) = work2(i,k,j,n)
602                   m = m+1
603                ENDDO
604             ENDDO
605          ENDDO
606
607          CALL fft_x_m( fftx_ar3(:,:,j), 'forward' )
608
609          DO  k = 1, nz
610             DO  i = nxl_a, nxr_a
611                tri_ar(i,k) = fftx_ar3(i,k,j)
612             ENDDO
613          ENDDO
614 
615          jj = myid * (nyn_p-nys_p+1) + j
616          CALL tridia_hybrid( jj, tri_ar, tri(:,:,:,jthread))
617
618          DO  k = 1, nz
619             DO  i = nxl_a, nxr_a
620                fftx_ar3(i,k,j) = tri_ar (i,k)
621             ENDDO
622          ENDDO
623
624          CALL fft_x_m( fftx_ar3(:,:,j), 'backward' )
625
626          DO  k = 1, nz
627             m = nxl_a
628             DO  n = 1, npe_s
629                DO  i = nxl_p, nxr_p
630                   work2(i,k,j,n) = fftx_ar3(m,k,j)
631                   m = m+1
632                ENDDO
633             ENDDO
634          ENDDO
635
636       ENDDO
637#if defined( __KKMP )
638!$OMP  END PARALLEL
639#endif
640
641       CALL cpu_log( log_point_s(33), 'fft_x_m + tridia', 'stop' )
642
643#if defined( __parallel )
644       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' ) 
645       nwords = (nxr_p-nxl_p+1) * nz * (nyn_p-nys_p+1)
646       CALL MPI_ALLTOALL( work2(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
647                          work1(nxl_p,1,nys_p,1), nwords, MPI_REAL, &
648                          comm2d, istat )
649       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
650#else
651       work1 = work2
652#endif
653
654       CALL cpu_log( log_point_s(7), 'fft_y_m', 'continue' )
655
656!$OMP  PARALLEL PRIVATE (i,iouter,ii,ir,iei,j,k,m,n)
657!$OMP  DO
658       DO  k = 1, nz
659          m = nys_a
660          DO  n = 1, npe_s
661             DO  j = nys_p, nyn_p
662                 DO  i = nxl_p, nxr_p
663                     ffty_ar3(m,k,i) = work1(i,k,j,n)
664                 ENDDO
665                 m = m+1
666             ENDDO
667          ENDDO
668       ENDDO 
669
670!$OMP  DO
671       DO  i = nxl_p, nxr_p
672          CALL fft_y_m( ffty_ar3(:,:,i), ny+3, 'backward' )
673          DO  j = nys_a, nyn_a
674             DO  k = 1, nz
675                ar(k,j,i+nxl) = ffty_ar3(j,k,i)
676             ENDDO
677          ENDDO
678       ENDDO
679!$OMP  END PARALLEL
680
681       CALL cpu_log( log_point_s(7), 'fft_y_m', 'stop' )
682
683       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_vec', 'stop' )
684
685#if defined( __KKMP )
686       DEALLOCATE( tri )
687#endif
688
689    END SUBROUTINE poisfft_hybrid_omp_vec
690
691
692    SUBROUTINE poisfft_hybrid_nodes ( ar )
693
694       USE cpulog
695       USE interfaces
696
697       IMPLICIT NONE
698
699       INTEGER, PARAMETER ::  istride = 4  ! stride of i loop
700       INTEGER            ::  i, iei, ii, iouter, ir, istat, j, jj, k, m, &
701                              n, nn, nt, nw1, nw2
702
703       REAL, DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar
704
705       REAL, DIMENSION(0:nx)                 ::  fftx_ar
706       REAL, DIMENSION(0:ny,istride)         ::  ffty_ar
707 
708       REAL, DIMENSION(0:nx,nz)              ::  tri_ar
709
710       REAL, DIMENSION(nxl_p:nxr_p,nz,tasks_per_logical_node, &
711                          nodes,nys_p:nyn_p) ::  work1,work2
712       REAL, DIMENSION(5,0:nx,0:nz-1)        ::  tri
713
714
715       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_nodes', 'start' )
716
717       CALL cpu_log( log_point_s(7), 'fft_y_1d', 'start' )
718
719!
720!--    Store grid points to be transformed on a 1d-array, do the fft
721!--    and sample the results on a 4d-array
722       DO  iouter = nxl_p, nxr_p, istride  ! stride loop, better cache
723          iei = MIN( iouter+istride-1, nxr_p )
724          DO  k = 1, nz
725
726             DO  i = iouter, iei
727                ii = nxl + i
728                ir = i - iouter + 1
729
730                DO  j = nys_a, nyn_a
731                   ffty_ar(j,ir) = ar(k,j,ii)
732                ENDDO
733   
734                CALL fft_y_1d( ffty_ar(:,ir), 'forward' )
735             ENDDO
736
737             m = nys_a
738             DO  nn = 1, nodes
739                DO  nt = 1, tasks_per_logical_node
740                   DO  j = nys_p, nyn_p
741                      DO  i = iouter, iei
742                         ir = i - iouter + 1
743                         work1(i,k,nt,nn,j) = ffty_ar(m,ir)
744                      ENDDO
745                      m = m+1
746                   ENDDO
747                ENDDO
748             ENDDO
749   
750          ENDDO
751       ENDDO
752
753       CALL cpu_log( log_point_s(7), 'fft_y_1d', 'pause' )
754
755       CALL cpu_log( log_point_s(32), 'alltoall_task', 'start' )
756       nw1 = SIZE( work1, 1 ) * SIZE( work1, 2 )
757       DO  nn = 1, nodes
758           DO  j = nys_p, nyn_p
759#if defined( __parallel )
760              CALL MPI_ALLTOALL( work1(nxl_p,1,1,nn,j), nw1, MPI_REAL, &
761                                 work2(nxl_p,1,1,nn,j), nw1, MPI_REAL, &
762                                                     comm_tasks, istat )
763#endif
764           ENDDO
765       ENDDO
766       CALL cpu_log( log_point_s(32), 'alltoall_task', 'stop' )
767
768
769       DO  j = nys_p, nyn_p
770       
771          CALL cascade( 1, j, nys_p, nyn_p )
772          nw2 = nw1 * SIZE( work1, 3 )
773          CALL cpu_log( log_point_s(37), 'alltoall_node', 'start' )
774#if defined( __parallel )
775          CALL MPI_ALLTOALL( work2(nxl_p,1,1,1,j), nw2, MPI_REAL, &
776                             work1(nxl_p,1,1,1,j), nw2, MPI_REAL, &
777                                                  comm_nodes, istat )
778#endif
779          CALL cpu_log( log_point_s(37), 'alltoall_node', 'pause' )
780          CALL cascade( 2, j, nys_p, nyn_p )
781
782          CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'start' )
783          DO  k = 1, nz
784
785             m = nxl_a
786             DO  nn = 1, nodes
787                DO  nt = 1, tasks_per_logical_node
788                   DO  i = nxl_p, nxr_p
789                      fftx_ar(m) = work1(i,k,nt,nn,j)
790                      m = m+1
791                   ENDDO
792                ENDDO
793             ENDDO
794
795             CALL fft_x_1d( fftx_ar, 'forward' )
796
797             DO  i = nxl_a, nxr_a
798                tri_ar(i,k) = fftx_ar(i)
799             ENDDO
800
801          ENDDO
802 
803          jj = myid * (nyn_p-nys_p+1) + j
804          CALL tridia_hybrid( jj, tri_ar, tri(:,:,:) )
805
806          DO  k = 1, nz
807             DO  i = nxl_a, nxr_a
808                fftx_ar(i) = tri_ar(i,k)
809             ENDDO
810
811             CALL fft_x_1d( fftx_ar, 'backward' )
812
813             m = nxl_a
814             DO  nn = 1, nodes
815                DO  nt = 1, tasks_per_logical_node
816                   DO  i = nxl_p, nxr_p
817                      work1(i,k,nt,nn,j) = fftx_ar(m)
818                      m = m+1
819                   ENDDO
820                ENDDO
821             ENDDO
822          ENDDO
823
824          CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'stop' )
825          nw2 = nw1 * SIZE( work1, 3 )
826          CALL cpu_log( log_point_s(37), 'alltoall_node', 'continue' )
827#if defined( __parallel )
828          CALL MPI_ALLTOALL( work1(nxl_p,1,1,1,j), nw2, MPI_REAL, &
829                             work2(nxl_p,1,1,1,j), nw2, MPI_REAL, &
830                                                  comm_nodes, istat )
831#endif
832          CALL cpu_log( log_point_s(37), 'alltoall_node', 'stop' )
833
834       ENDDO
835
836       CALL cpu_log( log_point_s(32), 'alltoall_task', 'start' ) 
837       DO  nn = 1, nodes
838           DO  j = nys_p, nyn_p
839#if defined( __parallel )
840              CALL MPI_ALLTOALL( work2(nxl_p,1,1,nn,j), nw1, MPI_REAL, &
841                                 work1(nxl_p,1,1,nn,j), nw1, MPI_REAL, &
842                                                     comm_tasks, istat )
843#endif
844           ENDDO
845       ENDDO
846       CALL cpu_log( log_point_s(32), 'alltoall_task', 'stop' )
847
848       CALL cpu_log( log_point_s(7), 'fft_y_1d', 'continue' )
849
850       DO  iouter = nxl_p, nxr_p, istride
851          iei = MIN( iouter+istride-1, nxr_p )
852          DO  k = 1, nz
853
854             m = nys_a
855             DO  nn = 1, nodes
856                DO  nt = 1, tasks_per_logical_node
857                   DO  j = nys_p, nyn_p
858                      DO  i = iouter, iei
859                         ir = i - iouter + 1
860                         ffty_ar(m,ir) = work1(i,k,nt,nn,j)
861                      ENDDO
862                      m = m+1
863                   ENDDO
864                ENDDO
865             ENDDO
866
867             DO  i = iouter, iei
868                ii = nxl + i
869                ir = i - iouter + 1
870                CALL fft_y_1d( ffty_ar(:,ir), 'backward' )
871
872                DO  j = nys_a, nyn_a
873                   ar(k,j,ii) = ffty_ar(j,ir)
874                ENDDO
875             ENDDO
876
877          ENDDO
878       ENDDO
879
880       CALL cpu_log( log_point_s(7), 'fft_y_1d', 'stop' )
881
882       CALL cpu_log( log_point_s(30), 'poisfft_hybrid_nodes', 'stop' )
883
884    END SUBROUTINE poisfft_hybrid_nodes
885
886
887
888    SUBROUTINE tridia_hybrid( j, ar, tri )
889
890       USE arrays_3d
891       USE control_parameters
892       USE grid_variables
893
894       IMPLICIT NONE
895
896       INTEGER ::  i, j, k, nnyh
897       REAL, DIMENSION(0:nx,nz)       ::  ar
898       REAL, DIMENSION(0:nx,0:nz-1)   ::  ar1
899       REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
900
901       nnyh = (ny+1) / 2
902
903       tri = 0.0
904!
905!--    Define constant elements of the tridiagonal matrix.
906       DO  k = 0, nz-1
907          DO  i = 0,nx
908             tri(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
909             tri(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
910          ENDDO
911       ENDDO
912
913       IF ( j <= nnyh )  THEN
914          CALL maketri_hybrid( j )
915       ELSE
916          CALL maketri_hybrid( ny+1-j)
917       ENDIF
918       CALL zerleg_hybrid
919       CALL substi_hybrid( ar, tri )
920
921    CONTAINS
922
923       SUBROUTINE maketri_hybrid( j )
924
925!----------------------------------------------------------------------!
926!                         maketri                                      !
927!                                                                      !
928!       computes the i- and j-dependent component of the matrix        !
929!----------------------------------------------------------------------!
930
931          USE constants
932
933          IMPLICIT NONE
934
935          INTEGER ::  i, j, k, nnxh
936          REAL    ::  a, c
937
938          REAL, DIMENSION(0:nx) ::  l
939
940
941          nnxh = (nx+1) / 2
942!
943!--       Provide the tridiagonal matrix for solution of the Poisson equation
944!--       in Fourier space. The coefficients are computed following the method
945!--       of Schmidt et al. (DFVLR-Mitteilung 84-15) --> departs from Stephan
946!--       Siano's original version.
947          DO  i = 0,nx
948             IF ( i >= 0 .AND. i < nnxh ) THEN
949                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
950                                        REAL( nx+1 ) ) ) / ( dx * dx ) + &
951                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
952                                        REAL( ny+1 ) ) ) / ( dy * dy )
953             ELSEIF ( i == nnxh ) THEN
954                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
955                                         REAL( nx+1 ) ) ) / ( dx * dx ) + &
956                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
957                                         REAL(ny+1) ) ) / ( dy * dy )
958             ELSE
959                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
960                                         REAL( nx+1 ) ) ) / ( dx * dx ) + &
961                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
962                                         REAL( ny+1 ) ) ) / ( dy * dy )
963             ENDIF
964          ENDDO
965
966          DO  k = 0,nz-1
967             DO  i = 0, nx
968                a = -1.0 * ddzu_pres(k+2) * ddzw(k+1)
969                c = -1.0 * ddzu_pres(k+1) * ddzw(k+1)
970                tri(1,i,k) = a + c - l(i)
971             ENDDO
972          ENDDO
973          IF ( ibc_p_b == 1 )  THEN
974             DO  i = 0,nx
975                tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
976             ENDDO
977          ENDIF
978          IF ( ibc_p_t == 1 )  THEN
979             DO  i = 0,nx
980                tri(1,i,nz-1) = tri(1,i,nz-1) + tri(3,i,nz-1)
981             ENDDO
982          ENDIF
983
984       END SUBROUTINE maketri_hybrid
985
986
987       SUBROUTINE zerleg_hybrid
988
989!----------------------------------------------------------------------!
990!                          zerleg                                      !
991!                                                                      !
992!       Splitting of the tridiagonal matrix (Thomas algorithm)         !
993!----------------------------------------------------------------------!
994
995          USE indices
996
997          IMPLICIT NONE
998
999          INTEGER ::  i, k
1000
1001!
1002!--       Splitting
1003          DO  i = 0, nx
1004             tri(4,i,0) = tri(1,i,0)
1005          ENDDO
1006          DO  k = 1, nz-1
1007             DO  i = 0,nx
1008                tri(5,i,k) = tri(2,i,k) / tri(4,i,k-1)
1009                tri(4,i,k) = tri(1,i,k) - tri(3,i,k-1) * tri(5,i,k)
1010             ENDDO
1011          ENDDO
1012
1013       END SUBROUTINE zerleg_hybrid
1014
1015       SUBROUTINE substi_hybrid( ar, tri )
1016
1017!----------------------------------------------------------------------!
1018!                          substi                                      !
1019!                                                                      !
1020!       Substitution (Forward and Backward) (Thomas algorithm)         !
1021!----------------------------------------------------------------------!
1022
1023          IMPLICIT NONE
1024
1025          INTEGER ::  i, j, k
1026          REAL, DIMENSION(0:nx,nz)       ::  ar
1027          REAL, DIMENSION(0:nx,0:nz-1)   ::  ar1
1028          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
1029
1030!
1031!--       Forward substitution
1032          DO  i = 0, nx
1033             ar1(i,0) = ar(i,1)
1034          ENDDO
1035          DO  k = 1, nz - 1
1036             DO  i = 0,nx
1037                ar1(i,k) = ar(i,k+1) - tri(5,i,k) * ar1(i,k-1)
1038             ENDDO
1039          ENDDO
1040
1041!
1042!--       Backward substitution
1043          DO  i = 0,nx
1044             ar(i,nz) = ar1(i,nz-1) / tri(4,i,nz-1)
1045          ENDDO
1046          DO  k = nz-2, 0, -1
1047             DO  i = 0,nx
1048                ar(i,k+1) = ( ar1(i,k) - tri(3,i,k) * ar(i,k+2) ) &
1049                            / tri(4,i,k)
1050             ENDDO
1051          ENDDO
1052
1053       END SUBROUTINE substi_hybrid
1054
1055    END SUBROUTINE tridia_hybrid
1056
1057
1058    SUBROUTINE cascade( loca, j, nys_p, nyn_p )
1059
1060       USE cpulog
1061
1062       IMPLICIT NONE
1063
1064       INTEGER ::  ier, j, loca, nyn_p, nys_p, req, reqa(1)
1065       INTEGER, SAVE ::  tag = 10
1066#if defined( __parallel )
1067       INTEGER, DIMENSION(MPI_STATUS_SIZE)   :: stat
1068       INTEGER, DIMENSION(MPI_STATUS_SIZE,1) :: stata
1069#endif
1070
1071       REAL ::  buf, buf1
1072
1073
1074       buf  = 1.0
1075       buf1 = 1.1
1076       IF ( me_node == 0 )  THEN  ! first node only
1077
1078          SELECT CASE ( loca )
1079
1080             CASE ( 1 )  ! before alltoall
1081
1082                IF( me_task > 0 )  THEN  ! first task does not wait
1083#if defined( __parallel )
1084                   CALL MPI_SENDRECV( buf,  1, MPI_REAL, me_task-1, 0, &
1085                                      buf1, 1, MPI_REAL, me_task-1, 0, &
1086                                                 comm_tasks, stat, ierr )
1087#endif
1088                ELSEIF ( j > nys_p )  THEN
1089                   req = 0
1090                   tag = MOD( tag-10, 10 ) + 10
1091#if defined( __parallel )
1092                   CALL MPI_IRECV( buf, 1, MPI_REAL, tasks_per_logical_node-1,&
1093                                   tag, comm_tasks, req, ierr )
1094                   reqa = req
1095                   CALL MPI_WAITALL( 1, reqa, stata, ierr )
1096#endif
1097                ENDIF
1098
1099             CASE ( 2 )  ! after alltoall
1100
1101                IF ( me_task < tasks_per_logical_node-1 )  THEN  ! last task
1102#if defined( __parallel )
1103                   CALL MPI_SENDRECV( buf,  1, MPI_REAL, me_task+1, 0, &
1104                                      buf1, 1, MPI_REAL, me_task+1, 0, &
1105                                                 comm_tasks, stat, ierr)
1106#endif
1107                ELSEIF ( j < nyn_p )  THEN
1108                   req = 0
1109                   tag = MOD( tag-10, 10 ) + 10
1110#if defined( __parallel )
1111                   CALL MPI_ISEND( buf, 1, MPI_REAL, 0, tag, comm_tasks, req, &
1112                                   ierr )
1113#endif
1114                ENDIF
1115
1116          END SELECT
1117
1118       ENDIF
1119
1120    END SUBROUTINE cascade
1121#endif
1122 END MODULE poisfft_hybrid_mod
Note: See TracBrowser for help on using the repository browser.