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

Last change on this file since 1053 was 1037, checked in by raasch, 11 years ago

last commit documented

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