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

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

fortran bugfix in subroutine cascade concerning dimension os MPI status variable

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