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

Last change on this file since 1107 was 1107, checked in by raasch, 9 years ago

last commit documented

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