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

Last change on this file since 1 was 1, checked in by raasch, 15 years ago

Initial repository layout and content

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