Home

Context Navigation

poisfft.f90 @ 1842

Last change on this file since 1842 was 1805, checked in by maronga, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 45.5 KB

Line
1	!> @file poisfft.f90
2	!--------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the terms
6	! of the GNU General Public License as published by the Free Software Foundation,
7	! either version 3 of the License, or (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with
14	! PALM. If not, see <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2014 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------!
18	!
19	! Current revisions:
20	! -----------------
21	!
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: poisfft.f90 1805 2016-04-05 16:32:34Z raasch $
26	!
27	! 1804 2016-04-05 16:30:18Z maronga
28	! Removed code for parameter file check (__check)
29	!
30	! 1682 2015-10-07 23:56:08Z knoop
31	! Code annotations made doxygen readable
32	!
33	! 1482 2014-10-18 12:34:45Z raasch
34	! use 2d-decomposition, if accelerator boards are used
35	!
36	! 1406 2014-05-16 13:47:01Z raasch
37	! bugfix for pgi 14.4: declare create moved after array declaration
38	!
39	! 1320 2014-03-20 08:40:49Z raasch
40	! ONLY-attribute added to USE-statements,
41	! kind-parameters added to all INTEGER and REAL declaration statements,
42	! kinds are defined in new module kinds,
43	! old module precision_kind is removed,
44	! revision history before 2012 removed,
45	! comment fields (!:) to be used for variable explanations added to
46	! all variable declaration statements
47	!
48	! 1318 2014-03-17 13:35:16Z raasch
49	! module interfaces removed
50	!
51	! 1306 2014-03-13 14:30:59Z raasch
52	! openmp sections removed from the overlap branch,
53	! second argument removed from parameter list
54	!
55	! 1216 2013-08-26 09:31:42Z raasch
56	! resorting of arrays moved to separate routines resort_for_...,
57	! one argument, used as temporary work array, removed from all transpose
58	! routines
59	! overlapping fft / transposition implemented
60	!
61	! 1212 2013-08-15 08:46:27Z raasch
62	! tridia routines moved to seperate module tridia_solver
63	!
64	! 1208 2013-08-13 06:41:49Z raasch
65	! acc-update clauses added for "ar" so that ffts other than cufft can also be
66	! used (although they are not ported and will give a poor performance)
67	!
68	! 1111 2013-03-08 23:54:10Z raasch
69	! further openACC porting of non-parallel (MPI) branch:
70	! tridiagonal routines split into extermal subroutines (instead using CONTAINS),
71	! no distinction between parallel/non-parallel in poisfft and tridia any more,
72	! tridia routines moved to end of file because of probable bug in PGI compiler 12.5
73	! (otherwise "invalid device function" is indicated during runtime),
74	! optimization of tridia routines: constant elements and coefficients of tri are
75	! stored in seperate arrays ddzuw and tric, last dimension of tri reduced from 5
76	! to 2,
77	! poisfft_init is now called internally from poisfft, maketri is called from
78	! poisfft_init,
79	! ibc_p_b = 2 removed
80	!
81	! 1106 2013-03-04 05:31:38Z raasch
82	! routines fftx, ffty, fftxp, fftyp removed, calls replaced by fft_x, fft_y,
83	! in the 1D-decomposition routines fft_x, ffty are replaced by fft_x_1d,
84	! fft_y_1d
85	!
86	! 1103 2013-02-20 02:15:53Z raasch
87	! tri, ar, and ar1 arguments in tridia-routines (2d) are removed because they
88	! sometimes cause segmentation faults with intel 12.1 compiler
89	!
90	! 1092 2013-02-02 11:24:22Z raasch
91	! unused variables removed
92	!
93	! 1036 2012-10-22 13:43:42Z raasch
94	! code put under GPL (PALM 3.9)
95	!
96	! 2012-09-21 07:03:55Z raasch
97	! FLOAT type conversion replaced by REAL
98	!
99	! 1003 2012-09-14 14:35:53Z raasch
100	! indices nxa, nya, etc. replaced by nx, ny, etc.
101	!
102	! 940 2012-07-09 14:31:00Z raasch
103	! special handling of tri-array as an argument in tridia_1dd routines switched
104	! off because it caused segmentation faults with intel 12.1 compiler
105	!
106	! 877 2012-04-03 11:21:44Z suehring
107	! Bugfix: Avoid divisions by zero in case of using a 'neumann' bc for the
108	! pressure at the top of the model domain.
109	!
110	! 809 2012-01-30 13:32:58Z maronga
111	! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
112	!
113	! 807 2012-01-25 11:53:51Z maronga
114	! New cpp directive "__check" implemented which is used by check_namelist_files
115	! (most of the code is unneeded by check_namelist_files).
116	!
117	! Revision 1.1 1997/07/24 11:24:14 raasch
118	! Initial revision
119	!
120	!
121	! Description:
122	! ------------
123	!> Solves the Poisson equation with a 2D spectral method
124	!> d^2 p / dx^2 + d^2 p / dy^2 + d^2 p / dz^2 = s
125	!>
126	!> Input:
127	!> real ar contains (nnz,nny,nnx) elements of the velocity divergence,
128	!> starting from (1,nys,nxl)
129	!>
130	!> Output:
131	!> real ar contains the solution for perturbation pressure p
132	!------------------------------------------------------------------------------!
133	MODULE poisfft_mod
134
135
136	USE fft_xy, &
137	ONLY: fft_init, fft_y, fft_y_1d, fft_y_m, fft_x, fft_x_1d, fft_x_m
138
139	USE indices, &
140	ONLY: nnx, nny, nx, nxl, nxr, ny, nys, nyn, nz
141
142	USE transpose_indices, &
143	ONLY: nxl_y, nxl_z, nxr_y, nxr_z, nys_x, nys_z, nyn_x, nyn_z, nzb_x, &
144	nzb_y, nzt_x, nzt_y
145
146	USE tridia_solver, &
147	ONLY: tridia_1dd, tridia_init, tridia_substi, tridia_substi_overlap
148
149	IMPLICIT NONE
150
151	LOGICAL, SAVE :: poisfft_initialized = .FALSE.
152
153	PRIVATE
154
155	PUBLIC poisfft, poisfft_init
156
157	INTERFACE poisfft
158	MODULE PROCEDURE poisfft
159	END INTERFACE poisfft
160
161	INTERFACE poisfft_init
162	MODULE PROCEDURE poisfft_init
163	END INTERFACE poisfft_init
164
165
166	CONTAINS
167
168	!------------------------------------------------------------------------------!
169	! Description:
170	! ------------
171	!> @todo Missing subroutine description.
172	!------------------------------------------------------------------------------!
173	SUBROUTINE poisfft_init
174
175	USE arrays_3d, &
176	ONLY: ddzu_pres, ddzw
177
178	USE kinds
179
180	IMPLICIT NONE
181
182	INTEGER(iwp) :: k !<
183
184
185	CALL fft_init
186
187	CALL tridia_init
188
189	poisfft_initialized = .TRUE.
190
191	END SUBROUTINE poisfft_init
192
193
194
195	!------------------------------------------------------------------------------!
196	! Description:
197	! ------------
198	!> Two-dimensional Fourier Transformation in x- and y-direction.
199	!------------------------------------------------------------------------------!
200	SUBROUTINE poisfft( ar )
201
202	USE control_parameters, &
203	ONLY: fft_method, transpose_compute_overlap
204
205	USE cpulog, &
206	ONLY: cpu_log, cpu_log_nowait, log_point_s
207
208	USE kinds
209
210	USE pegrid
211
212	IMPLICIT NONE
213
214	INTEGER(iwp) :: ii !<
215	INTEGER(iwp) :: iind !<
216	INTEGER(iwp) :: inew !<
217	INTEGER(iwp) :: jj !<
218	INTEGER(iwp) :: jind !<
219	INTEGER(iwp) :: jnew !<
220	INTEGER(iwp) :: ki !<
221	INTEGER(iwp) :: kk !<
222	INTEGER(iwp) :: knew !<
223	INTEGER(iwp) :: n !<
224	INTEGER(iwp) :: nblk !<
225	INTEGER(iwp) :: nnx_y !<
226	INTEGER(iwp) :: nny_z !<
227	INTEGER(iwp) :: nnz_t !<
228	INTEGER(iwp) :: nnz_x !<
229	INTEGER(iwp) :: nxl_y_bound !<
230	INTEGER(iwp) :: nxr_y_bound !<
231
232	INTEGER(iwp), DIMENSION(4) :: isave !<
233
234	REAL(wp), DIMENSION(1:nz,nys:nyn,nxl:nxr) :: ar !<
235	REAL(wp), DIMENSION(nys:nyn,nxl:nxr,1:nz) :: ar_inv !<
236	!$acc declare create( ar_inv )
237
238	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ar1 !<
239	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_in !<
240	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_inv !<
241	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_out_y !<
242	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_out_z !<
243
244
245	CALL cpu_log( log_point_s(3), 'poisfft', 'start' )
246
247	IF ( .NOT. poisfft_initialized ) CALL poisfft_init
248
249	!
250	!-- Two-dimensional Fourier Transformation in x- and y-direction.
251	IF ( pdims(2) == 1 .AND. pdims(1) > 1 .AND. num_acc_per_node == 0 ) &
252	THEN
253
254	!
255	!-- 1d-domain-decomposition along x:
256	!-- FFT along y and transposition y --> x
257	CALL ffty_tr_yx( ar, ar )
258
259	!
260	!-- FFT along x, solving the tridiagonal system and backward FFT
261	CALL fftx_tri_fftx( ar )
262
263	!
264	!-- Transposition x --> y and backward FFT along y
265	CALL tr_xy_ffty( ar, ar )
266
267	ELSEIF ( pdims(1) == 1 .AND. pdims(2) > 1 .AND. num_acc_per_node == 0 ) &
268	THEN
269
270	!
271	!-- 1d-domain-decomposition along y:
272	!-- FFT along x and transposition x --> y
273	CALL fftx_tr_xy( ar, ar )
274
275	!
276	!-- FFT along y, solving the tridiagonal system and backward FFT
277	CALL ffty_tri_ffty( ar )
278
279	!
280	!-- Transposition y --> x and backward FFT along x
281	CALL tr_yx_fftx( ar, ar )
282
283	ELSEIF ( .NOT. transpose_compute_overlap ) THEN
284
285	!
286	!-- 2d-domain-decomposition or no decomposition (1 PE run)
287	!-- Transposition z --> x
288	CALL cpu_log( log_point_s(5), 'transpo forward', 'start' )
289	CALL resort_for_zx( ar, ar_inv )
290	CALL transpose_zx( ar_inv, ar )
291	CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
292
293	CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
294	IF ( fft_method /= 'system-specific' ) THEN
295	!$acc update host( ar )
296	ENDIF
297	CALL fft_x( ar, 'forward' )
298	IF ( fft_method /= 'system-specific' ) THEN
299	!$acc update device( ar )
300	ENDIF
301	CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
302
303	!
304	!-- Transposition x --> y
305	CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
306	CALL resort_for_xy( ar, ar_inv )
307	CALL transpose_xy( ar_inv, ar )
308	CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
309
310	CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
311	IF ( fft_method /= 'system-specific' ) THEN
312	!$acc update host( ar )
313	ENDIF
314	CALL fft_y( ar, 'forward', ar_tr = ar, &
315	nxl_y_bound = nxl_y, nxr_y_bound = nxr_y, &
316	nxl_y_l = nxl_y, nxr_y_l = nxr_y )
317	IF ( fft_method /= 'system-specific' ) THEN
318	!$acc update device( ar )
319	ENDIF
320	CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
321
322	!
323	!-- Transposition y --> z
324	CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
325	CALL resort_for_yz( ar, ar_inv )
326	CALL transpose_yz( ar_inv, ar )
327	CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
328
329	!
330	!-- Solve the tridiagonal equation system along z
331	CALL cpu_log( log_point_s(6), 'tridia', 'start' )
332	CALL tridia_substi( ar )
333	CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
334
335	!
336	!-- Inverse Fourier Transformation
337	!-- Transposition z --> y
338	CALL cpu_log( log_point_s(8), 'transpo invers', 'start' )
339	CALL transpose_zy( ar, ar_inv )
340	CALL resort_for_zy( ar_inv, ar )
341	CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
342
343	CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
344	IF ( fft_method /= 'system-specific' ) THEN
345	!$acc update host( ar )
346	ENDIF
347	CALL fft_y( ar, 'backward', ar_tr = ar, &
348	nxl_y_bound = nxl_y, nxr_y_bound = nxr_y, &
349	nxl_y_l = nxl_y, nxr_y_l = nxr_y )
350	IF ( fft_method /= 'system-specific' ) THEN
351	!$acc update device( ar )
352	ENDIF
353	CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
354
355	!
356	!-- Transposition y --> x
357	CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
358	CALL transpose_yx( ar, ar_inv )
359	CALL resort_for_yx( ar_inv, ar )
360	CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
361
362	CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
363	IF ( fft_method /= 'system-specific' ) THEN
364	!$acc update host( ar )
365	ENDIF
366	CALL fft_x( ar, 'backward' )
367	IF ( fft_method /= 'system-specific' ) THEN
368	!$acc update device( ar )
369	ENDIF
370	CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
371
372	!
373	!-- Transposition x --> z
374	CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
375	CALL transpose_xz( ar, ar_inv )
376	CALL resort_for_xz( ar_inv, ar )
377	CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
378
379	ELSE
380
381	!
382	!-- 2d-domain-decomposition or no decomposition (1 PE run) with
383	!-- overlapping transposition / fft
384	!-- cputime logging must not use barriers, which would prevent overlapping
385	ALLOCATE( f_out_y(0:ny,nxl_y:nxr_y,nzb_y:nzt_y), &
386	f_out_z(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
387	!
388	!-- Transposition z --> x + subsequent fft along x
389	ALLOCATE( f_inv(nys:nyn,nxl:nxr,1:nz) )
390	CALL resort_for_zx( ar, f_inv )
391	!
392	!-- Save original indices and gridpoint counter
393	isave(1) = nz
394	isave(2) = nzb_x
395	isave(3) = nzt_x
396	isave(4) = sendrecvcount_zx
397	!
398	!-- Set new indices for transformation
399	nblk = nz / pdims(1)
400	nz = pdims(1)
401	nnz_x = 1
402	nzb_x = 1 + myidx * nnz_x
403	nzt_x = ( myidx + 1 ) * nnz_x
404	sendrecvcount_zx = nnx * nny * nnz_x
405
406	ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
407	ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) )
408
409	DO kk = 1, nblk
410
411	IF ( kk == 1 ) THEN
412	CALL cpu_log( log_point_s(5), 'transpo forward', 'start', cpu_log_nowait )
413	ELSE
414	CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
415	ENDIF
416
417	DO knew = 1, nz
418	ki = kk + nblk * ( knew - 1 )
419	f_in(:,:,knew) = f_inv(:,:,ki)
420	ENDDO
421
422	CALL transpose_zx( f_in, ar1(:,:,:))
423	CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
424
425	IF ( kk == 1 ) THEN
426	CALL cpu_log( log_point_s(4), 'fft_x', 'start', cpu_log_nowait )
427	ELSE
428	CALL cpu_log( log_point_s(4), 'fft_x', 'continue', cpu_log_nowait )
429	ENDIF
430
431	n = isave(2) + kk - 1
432	CALL fft_x( ar1(:,:,:), 'forward', ar_2d = f_out_z(:,:,n))
433	CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
434
435	ENDDO
436	!
437	!-- Restore original indices/counters
438	nz = isave(1)
439	nzb_x = isave(2)
440	nzt_x = isave(3)
441	sendrecvcount_zx = isave(4)
442
443	DEALLOCATE( ar1, f_in, f_inv )
444
445	!
446	!-- Transposition x --> y + subsequent fft along y
447	ALLOCATE( f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
448	CALL resort_for_xy( f_out_z, f_inv )
449	!
450	!-- Save original indices and gridpoint counter
451	isave(1) = nx
452	isave(2) = nxl_y
453	isave(3) = nxr_y
454	isave(4) = sendrecvcount_xy
455	!
456	!-- Set new indices for transformation
457	nblk = ( ( nx+1 ) / pdims(2) ) - 1
458	nx = pdims(2)
459	nnx_y = 1
460	nxl_y = myidy * nnx_y
461	nxr_y = ( myidy + 1 ) * nnx_y - 1
462	sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 )
463
464	ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) )
465	ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
466
467	DO ii = 0, nblk
468
469	CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
470
471	DO inew = 0, nx-1
472	iind = ii + ( nblk + 1 ) * inew
473	f_in(:,:,inew) = f_inv(:,:,iind)
474	ENDDO
475
476	CALL transpose_xy( f_in, ar1(:,:,:) )
477
478	CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
479
480	IF ( ii == 1 ) THEN
481	CALL cpu_log( log_point_s(7), 'fft_y', 'start', cpu_log_nowait )
482	ELSE
483	CALL cpu_log( log_point_s(7), 'fft_y', 'continue', cpu_log_nowait )
484	ENDIF
485
486	nxl_y_bound = isave(2)
487	nxr_y_bound = isave(3)
488	n = isave(2) + ii
489	CALL fft_y( ar1(:,:,:), 'forward', ar_tr = f_out_y, &
490	nxl_y_bound = nxl_y_bound, nxr_y_bound = nxr_y_bound, &
491	nxl_y_l = n, nxr_y_l = n )
492
493	CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
494
495	ENDDO
496	!
497	!-- Restore original indices/counters
498	nx = isave(1)
499	nxl_y = isave(2)
500	nxr_y = isave(3)
501	sendrecvcount_xy = isave(4)
502
503	DEALLOCATE( ar1, f_in, f_inv )
504
505	!
506	!-- Transposition y --> z + subsequent tridia + resort for z --> y
507	ALLOCATE( f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) )
508	CALL resort_for_yz( f_out_y, f_inv )
509	!
510	!-- Save original indices and gridpoint counter
511	isave(1) = ny
512	isave(2) = nys_z
513	isave(3) = nyn_z
514	isave(4) = sendrecvcount_yz
515	!
516	!-- Set new indices for transformation
517	nblk = ( ( ny+1 ) / pdims(1) ) - 1
518	ny = pdims(1)
519	nny_z = 1
520	nys_z = myidx * nny_z
521	nyn_z = ( myidx + 1 ) * nny_z - 1
522	sendrecvcount_yz = ( nxr_y-nxl_y+1 ) * nny_z * ( nzt_y-nzb_y+1 )
523
524	ALLOCATE( ar1(nxl_z:nxr_z,nys_z:nyn_z,1:nz) )
525	ALLOCATE( f_in(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) )
526
527	DO jj = 0, nblk
528	!
529	!-- Forward Fourier Transformation
530	!-- Transposition y --> z
531	CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
532
533	DO jnew = 0, ny-1
534	jind = jj + ( nblk + 1 ) * jnew
535	f_in(:,:,jnew) = f_inv(:,:,jind)
536	ENDDO
537
538	CALL transpose_yz( f_in, ar1(:,:,:) )
539
540	IF ( jj == nblk ) THEN
541	CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
542	ELSE
543	CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
544	ENDIF
545
546	!
547	!-- Solve the tridiagonal equation system along z
548	CALL cpu_log( log_point_s(6), 'tridia', 'start', cpu_log_nowait )
549
550	n = isave(2) + jj
551	CALL tridia_substi_overlap( ar1(:,:,:), n )
552
553	CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
554
555	!
556	!-- Inverse Fourier Transformation
557	!-- Transposition z --> y
558	!-- Only one thread should call MPI routines, therefore forward and
559	!-- backward tranpose are in the same section
560	IF ( jj == 0 ) THEN
561	CALL cpu_log( log_point_s(8), 'transpo invers', 'start', cpu_log_nowait )
562	ELSE
563	CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
564	ENDIF
565
566	CALL transpose_zy( ar1(:,:,:), f_in )
567
568	DO jnew = 0, ny-1
569	jind = jj + ( nblk + 1 ) * jnew
570	f_inv(:,:,jind) = f_in(:,:,jnew)
571	ENDDO
572
573	CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
574
575	ENDDO
576	!
577	!-- Restore original indices/counters
578	ny = isave(1)
579	nys_z = isave(2)
580	nyn_z = isave(3)
581	sendrecvcount_yz = isave(4)
582
583	CALL resort_for_zy( f_inv, f_out_y )
584
585	DEALLOCATE( ar1, f_in, f_inv )
586
587	!
588	!-- fft along y backward + subsequent transposition y --> x
589	ALLOCATE( f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
590	!
591	!-- Save original indices and gridpoint counter
592	isave(1) = nx
593	isave(2) = nxl_y
594	isave(3) = nxr_y
595	isave(4) = sendrecvcount_xy
596	!
597	!-- Set new indices for transformation
598	nblk = (( nx+1 ) / pdims(2) ) - 1
599	nx = pdims(2)
600	nnx_y = 1
601	nxl_y = myidy * nnx_y
602	nxr_y = ( myidy + 1 ) * nnx_y - 1
603	sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 )
604
605	ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) )
606	ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
607
608	DO ii = 0, nblk
609
610	CALL cpu_log( log_point_s(7), 'fft_y', 'continue', cpu_log_nowait )
611
612	n = isave(2) + ii
613	nxl_y_bound = isave(2)
614	nxr_y_bound = isave(3)
615
616	CALL fft_y( ar1(:,:,:), 'backward', ar_tr = f_out_y, &
617	nxl_y_bound = nxl_y_bound, nxr_y_bound = nxr_y_bound, &
618	nxl_y_l = n, nxr_y_l = n )
619
620	IF ( ii == nblk ) THEN
621	CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
622	ELSE
623	CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
624	ENDIF
625
626	CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
627
628	CALL transpose_yx( ar1(:,:,:), f_in )
629
630	DO inew = 0, nx-1
631	iind = ii + (nblk+1) * inew
632	f_inv(:,:,iind) = f_in(:,:,inew)
633	ENDDO
634
635	CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
636
637	ENDDO
638	!
639	!-- Restore original indices/counters
640	nx = isave(1)
641	nxl_y = isave(2)
642	nxr_y = isave(3)
643	sendrecvcount_xy = isave(4)
644
645	CALL resort_for_yx( f_inv, f_out_z )
646
647	DEALLOCATE( ar1, f_in, f_inv )
648
649	!
650	!-- fft along x backward + subsequent final transposition x --> z
651	ALLOCATE( f_inv(nys:nyn,nxl:nxr,1:nz) )
652	!
653	!-- Save original indices and gridpoint counter
654	isave(1) = nz
655	isave(2) = nzb_x
656	isave(3) = nzt_x
657	isave(4) = sendrecvcount_zx
658	!
659	!-- Set new indices for transformation
660	nblk = nz / pdims(1)
661	nz = pdims(1)
662	nnz_x = 1
663	nzb_x = 1 + myidx * nnz_x
664	nzt_x = ( myidx + 1 ) * nnz_x
665	sendrecvcount_zx = nnx * nny * nnz_x
666
667	ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
668	ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) )
669
670	DO kk = 1, nblk
671
672	CALL cpu_log( log_point_s(4), 'fft_x', 'continue', cpu_log_nowait )
673
674	n = isave(2) + kk - 1
675	CALL fft_x( ar1(:,:,:), 'backward', f_out_z(:,:,n))
676
677	IF ( kk == nblk ) THEN
678	CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
679	ELSE
680	CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
681	ENDIF
682
683	CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
684
685	CALL transpose_xz( ar1(:,:,:), f_in )
686
687	DO knew = 1, nz
688	ki = kk + nblk * (knew-1)
689	f_inv(:,:,ki) = f_in(:,:,knew)
690	ENDDO
691
692	IF ( kk == nblk ) THEN
693	CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
694	ELSE
695	CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
696	ENDIF
697
698	ENDDO
699	!
700	!-- Restore original indices/counters
701	nz = isave(1)
702	nzb_x = isave(2)
703	nzt_x = isave(3)
704	sendrecvcount_zx = isave(4)
705
706	CALL resort_for_xz( f_inv, ar )
707
708	DEALLOCATE( ar1, f_in, f_inv )
709
710	ENDIF
711
712	CALL cpu_log( log_point_s(3), 'poisfft', 'stop' )
713
714	END SUBROUTINE poisfft
715
716
717	!------------------------------------------------------------------------------!
718	! Description:
719	! ------------
720	!> Fourier-transformation along y with subsequent transposition y --> x for
721	!> a 1d-decomposition along x.
722	!>
723	!> @attention The performance of this routine is much faster on the NEC-SX6,
724	!> if the first index of work_ffty_vec is odd. Otherwise
725	!> memory bank conflicts may occur (especially if the index is a
726	!> multiple of 128). That's why work_ffty_vec is dimensioned as
727	!> 0:ny+1.
728	!> Of course, this will not work if users are using an odd number
729	!> of gridpoints along y.
730	!------------------------------------------------------------------------------!
731	SUBROUTINE ffty_tr_yx( f_in, f_out )
732
733
734	USE control_parameters, &
735	ONLY: host
736
737	USE cpulog, &
738	ONLY: cpu_log, log_point_s
739
740	USE kinds
741
742	USE pegrid
743
744	IMPLICIT NONE
745
746	INTEGER(iwp) :: i !<
747	INTEGER(iwp) :: iend !<
748	INTEGER(iwp) :: iouter !<
749	INTEGER(iwp) :: ir !<
750	INTEGER(iwp) :: j !<
751	INTEGER(iwp) :: k !<
752
753	INTEGER(iwp), PARAMETER :: stridex = 4 !<
754
755	REAL(wp), DIMENSION(0:ny,stridex) :: work_ffty !<
756	#if defined( __nec )
757	REAL(wp), DIMENSION(0:ny+1,1:nz,nxl:nxr) :: work_ffty_vec !<
758	#endif
759	REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr) :: f_in !<
760	REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) :: f_out !<
761	REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny) :: work !<
762
763	!
764	!-- Carry out the FFT along y, where all data are present due to the
765	!-- 1d-decomposition along x. Resort the data in a way that x becomes
766	!-- the first index.
767	CALL cpu_log( log_point_s(7), 'fft_y_1d', 'start' )
768
769	IF ( host(1:3) == 'nec' ) THEN
770	#if defined( __nec )
771	!
772	!-- Code optimized for vector processors
773	!$OMP PARALLEL PRIVATE ( i, j, k )
774	!$OMP DO
775	DO i = nxl, nxr
776
777	DO j = 0, ny
778	DO k = 1, nz
779	work_ffty_vec(j,k,i) = f_in(k,j,i)
780	ENDDO
781	ENDDO
782
783	CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'forward' )
784
785	ENDDO
786
787	!$OMP DO
788	DO k = 1, nz
789	DO j = 0, ny
790	DO i = nxl, nxr
791	work(i,k,j) = work_ffty_vec(j,k,i)
792	ENDDO
793	ENDDO
794	ENDDO
795	!$OMP END PARALLEL
796	#endif
797
798	ELSE
799
800	!
801	!-- Cache optimized code.
802	!-- The i-(x-)direction is split into a strided outer loop and an inner
803	!-- loop for better cache performance
804	!$OMP PARALLEL PRIVATE (i,iend,iouter,ir,j,k,work_ffty)
805	!$OMP DO
806	DO iouter = nxl, nxr, stridex
807
808	iend = MIN( iouter+stridex-1, nxr ) ! Upper bound for inner i loop
809
810	DO k = 1, nz
811
812	DO i = iouter, iend
813
814	ir = i-iouter+1 ! counter within a stride
815	DO j = 0, ny
816	work_ffty(j,ir) = f_in(k,j,i)
817	ENDDO
818	!
819	!-- FFT along y
820	CALL fft_y_1d( work_ffty(:,ir), 'forward' )
821
822	ENDDO
823
824	!
825	!-- Resort
826	DO j = 0, ny
827	DO i = iouter, iend
828	work(i,k,j) = work_ffty(j,i-iouter+1)
829	ENDDO
830	ENDDO
831
832	ENDDO
833
834	ENDDO
835	!$OMP END PARALLEL
836
837	ENDIF
838	CALL cpu_log( log_point_s(7), 'fft_y_1d', 'pause' )
839
840	!
841	!-- Transpose array
842	#if defined( __parallel )
843	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
844	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
845	CALL MPI_ALLTOALL( work(nxl,1,0), sendrecvcount_xy, MPI_REAL, &
846	f_out(1,1,nys_x,1), sendrecvcount_xy, MPI_REAL, &
847	comm1dx, ierr )
848	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
849	#endif
850
851	END SUBROUTINE ffty_tr_yx
852
853
854	!------------------------------------------------------------------------------!
855	! Description:
856	! ------------
857	!> Transposition x --> y with a subsequent backward Fourier transformation for
858	!> a 1d-decomposition along x
859	!------------------------------------------------------------------------------!
860	SUBROUTINE tr_xy_ffty( f_in, f_out )
861
862
863	USE control_parameters, &
864	ONLY: host
865
866	USE cpulog, &
867	ONLY: cpu_log, log_point_s
868
869	USE kinds
870
871	USE pegrid
872
873	IMPLICIT NONE
874
875	INTEGER(iwp) :: i !<
876	INTEGER(iwp) :: iend !<
877	INTEGER(iwp) :: iouter !<
878	INTEGER(iwp) :: ir !<
879	INTEGER(iwp) :: j !<
880	INTEGER(iwp) :: k !<
881
882	INTEGER(iwp), PARAMETER :: stridex = 4 !<
883
884	REAL(wp), DIMENSION(0:ny,stridex) :: work_ffty !<
885	#if defined( __nec )
886	REAL(wp), DIMENSION(0:ny+1,1:nz,nxl:nxr) :: work_ffty_vec !<
887	#endif
888	REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) :: f_in !<
889	REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr) :: f_out !<
890	REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny) :: work !<
891
892	!
893	!-- Transpose array
894	#if defined( __parallel )
895	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
896	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
897	CALL MPI_ALLTOALL( f_in(1,1,nys_x,1), sendrecvcount_xy, MPI_REAL, &
898	work(nxl,1,0), sendrecvcount_xy, MPI_REAL, &
899	comm1dx, ierr )
900	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
901	#endif
902
903	!
904	!-- Resort the data in a way that y becomes the first index and carry out the
905	!-- backward fft along y.
906	CALL cpu_log( log_point_s(7), 'fft_y_1d', 'continue' )
907
908	IF ( host(1:3) == 'nec' ) THEN
909	#if defined( __nec )
910	!
911	!-- Code optimized for vector processors
912	!$OMP PARALLEL PRIVATE ( i, j, k )
913	!$OMP DO
914	DO k = 1, nz
915	DO j = 0, ny
916	DO i = nxl, nxr
917	work_ffty_vec(j,k,i) = work(i,k,j)
918	ENDDO
919	ENDDO
920	ENDDO
921
922	!$OMP DO
923	DO i = nxl, nxr
924
925	CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'backward' )
926
927	DO j = 0, ny
928	DO k = 1, nz
929	f_out(k,j,i) = work_ffty_vec(j,k,i)
930	ENDDO
931	ENDDO
932
933	ENDDO
934	!$OMP END PARALLEL
935	#endif
936
937	ELSE
938
939	!
940	!-- Cache optimized code.
941	!-- The i-(x-)direction is split into a strided outer loop and an inner
942	!-- loop for better cache performance
943	!$OMP PARALLEL PRIVATE ( i, iend, iouter, ir, j, k, work_ffty )
944	!$OMP DO
945	DO iouter = nxl, nxr, stridex
946
947	iend = MIN( iouter+stridex-1, nxr ) ! Upper bound for inner i loop
948
949	DO k = 1, nz
950	!
951	!-- Resort
952	DO j = 0, ny
953	DO i = iouter, iend
954	work_ffty(j,i-iouter+1) = work(i,k,j)
955	ENDDO
956	ENDDO
957
958	DO i = iouter, iend
959
960	!
961	!-- FFT along y
962	ir = i-iouter+1 ! counter within a stride
963	CALL fft_y_1d( work_ffty(:,ir), 'backward' )
964
965	DO j = 0, ny
966	f_out(k,j,i) = work_ffty(j,ir)
967	ENDDO
968	ENDDO
969
970	ENDDO
971
972	ENDDO
973	!$OMP END PARALLEL
974
975	ENDIF
976
977	CALL cpu_log( log_point_s(7), 'fft_y_1d', 'stop' )
978
979	END SUBROUTINE tr_xy_ffty
980
981
982	!------------------------------------------------------------------------------!
983	! Description:
984	! ------------
985	!> FFT along x, solution of the tridiagonal system and backward FFT for
986	!> a 1d-decomposition along x
987	!>
988	!> @warning this subroutine may still not work for hybrid parallelization
989	!> with OpenMP (for possible necessary changes see the original
990	!> routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002)
991	!------------------------------------------------------------------------------!
992	SUBROUTINE fftx_tri_fftx( ar )
993
994
995	USE control_parameters, &
996	ONLY: host
997
998	USE cpulog, &
999	ONLY: cpu_log, log_point_s
1000
1001	USE grid_variables, &
1002	ONLY: ddx2, ddy2
1003
1004	USE kinds
1005
1006	USE pegrid
1007
1008	IMPLICIT NONE
1009
1010	INTEGER(iwp) :: i !<
1011	INTEGER(iwp) :: j !<
1012	INTEGER(iwp) :: k !<
1013	INTEGER(iwp) :: m !<
1014	INTEGER(iwp) :: n !<
1015	INTEGER(iwp) :: omp_get_thread_num !<
1016	INTEGER(iwp) :: tn !<
1017
1018	REAL(wp), DIMENSION(0:nx) :: work_fftx !<
1019	REAL(wp), DIMENSION(0:nx,1:nz) :: work_trix !<
1020	REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) :: ar !<
1021	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tri !<
1022
1023
1024	CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'start' )
1025
1026	ALLOCATE( tri(5,0:nx,0:nz-1,0:threads_per_task-1) )
1027
1028	tn = 0 ! Default thread number in case of one thread
1029	!$OMP PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_fftx, work_trix )
1030	DO j = nys_x, nyn_x
1031
1032	!$ tn = omp_get_thread_num()
1033
1034	IF ( host(1:3) == 'nec' ) THEN
1035	!
1036	!-- Code optimized for vector processors
1037	DO k = 1, nz
1038
1039	m = 0
1040	DO n = 1, pdims(1)
1041	DO i = 1, nnx
1042	work_trix(m,k) = ar(i,k,j,n)
1043	m = m + 1
1044	ENDDO
1045	ENDDO
1046
1047	ENDDO
1048
1049	CALL fft_x_m( work_trix, 'forward' )
1050
1051	ELSE
1052	!
1053	!-- Cache optimized code
1054	DO k = 1, nz
1055
1056	m = 0
1057	DO n = 1, pdims(1)
1058	DO i = 1, nnx
1059	work_fftx(m) = ar(i,k,j,n)
1060	m = m + 1
1061	ENDDO
1062	ENDDO
1063
1064	CALL fft_x_1d( work_fftx, 'forward' )
1065
1066	DO i = 0, nx
1067	work_trix(i,k) = work_fftx(i)
1068	ENDDO
1069
1070	ENDDO
1071
1072	ENDIF
1073
1074	!
1075	!-- Solve the linear equation system
1076	CALL tridia_1dd( ddx2, ddy2, nx, ny, j, work_trix, tri(:,:,:,tn) )
1077
1078	IF ( host(1:3) == 'nec' ) THEN
1079	!
1080	!-- Code optimized for vector processors
1081	CALL fft_x_m( work_trix, 'backward' )
1082
1083	DO k = 1, nz
1084
1085	m = 0
1086	DO n = 1, pdims(1)
1087	DO i = 1, nnx
1088	ar(i,k,j,n) = work_trix(m,k)
1089	m = m + 1
1090	ENDDO
1091	ENDDO
1092
1093	ENDDO
1094
1095	ELSE
1096	!
1097	!-- Cache optimized code
1098	DO k = 1, nz
1099
1100	DO i = 0, nx
1101	work_fftx(i) = work_trix(i,k)
1102	ENDDO
1103
1104	CALL fft_x_1d( work_fftx, 'backward' )
1105
1106	m = 0
1107	DO n = 1, pdims(1)
1108	DO i = 1, nnx
1109	ar(i,k,j,n) = work_fftx(m)
1110	m = m + 1
1111	ENDDO
1112	ENDDO
1113
1114	ENDDO
1115
1116	ENDIF
1117
1118	ENDDO
1119
1120	DEALLOCATE( tri )
1121
1122	CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'stop' )
1123
1124	END SUBROUTINE fftx_tri_fftx
1125
1126
1127	!------------------------------------------------------------------------------!
1128	! Description:
1129	! ------------
1130	!> Fourier-transformation along x with subsequent transposition x --> y for
1131	!> a 1d-decomposition along y.
1132	!>
1133	!> @attention NEC-branch of this routine may significantly profit from
1134	!> further optimizations. So far, performance is much worse than
1135	!> for routine ffty_tr_yx (more than three times slower).
1136	!------------------------------------------------------------------------------!
1137	SUBROUTINE fftx_tr_xy( f_in, f_out )
1138
1139
1140	USE control_parameters, &
1141	ONLY: host
1142
1143	USE cpulog, &
1144	ONLY: cpu_log, log_point_s
1145
1146	USE kinds
1147
1148	USE pegrid
1149
1150	IMPLICIT NONE
1151
1152	INTEGER(iwp) :: i !<
1153	INTEGER(iwp) :: j !<
1154	INTEGER(iwp) :: k !<
1155
1156	REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn) :: work_fftx !<
1157	REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx) :: f_in !<
1158	REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) :: f_out !<
1159	REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx) :: work !<
1160
1161	!
1162	!-- Carry out the FFT along x, where all data are present due to the
1163	!-- 1d-decomposition along y. Resort the data in a way that y becomes
1164	!-- the first index.
1165	CALL cpu_log( log_point_s(4), 'fft_x_1d', 'start' )
1166
1167	IF ( host(1:3) == 'nec' ) THEN
1168	!
1169	!-- Code for vector processors
1170	!$OMP PARALLEL PRIVATE ( i, j, k )
1171	!$OMP DO
1172	DO i = 0, nx
1173
1174	DO j = nys, nyn
1175	DO k = 1, nz
1176	work_fftx(i,k,j) = f_in(k,j,i)
1177	ENDDO
1178	ENDDO
1179
1180	ENDDO
1181
1182	!$OMP DO
1183	DO j = nys, nyn
1184
1185	CALL fft_x_m( work_fftx(:,:,j), 'forward' )
1186
1187	DO k = 1, nz
1188	DO i = 0, nx
1189	work(j,k,i) = work_fftx(i,k,j)
1190	ENDDO
1191	ENDDO
1192
1193	ENDDO
1194	!$OMP END PARALLEL
1195
1196	ELSE
1197
1198	!
1199	!-- Cache optimized code (there might be still a potential for better
1200	!-- optimization).
1201	!$OMP PARALLEL PRIVATE (i,j,k)
1202	!$OMP DO
1203	DO i = 0, nx
1204
1205	DO j = nys, nyn
1206	DO k = 1, nz
1207	work_fftx(i,k,j) = f_in(k,j,i)
1208	ENDDO
1209	ENDDO
1210
1211	ENDDO
1212
1213	!$OMP DO
1214	DO j = nys, nyn
1215	DO k = 1, nz
1216
1217	CALL fft_x_1d( work_fftx(0:nx,k,j), 'forward' )
1218
1219	DO i = 0, nx
1220	work(j,k,i) = work_fftx(i,k,j)
1221	ENDDO
1222	ENDDO
1223
1224	ENDDO
1225	!$OMP END PARALLEL
1226
1227	ENDIF
1228	CALL cpu_log( log_point_s(4), 'fft_x_1d', 'pause' )
1229
1230	!
1231	!-- Transpose array
1232	#if defined( __parallel )
1233	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1234	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1235	CALL MPI_ALLTOALL( work(nys,1,0), sendrecvcount_xy, MPI_REAL, &
1236	f_out(1,1,nxl_y,1), sendrecvcount_xy, MPI_REAL, &
1237	comm1dy, ierr )
1238	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1239	#endif
1240
1241	END SUBROUTINE fftx_tr_xy
1242
1243
1244	!------------------------------------------------------------------------------!
1245	! Description:
1246	! ------------
1247	!> Transposition y --> x with a subsequent backward Fourier transformation for
1248	!> a 1d-decomposition along x.
1249	!------------------------------------------------------------------------------!
1250	SUBROUTINE tr_yx_fftx( f_in, f_out )
1251
1252
1253	USE control_parameters, &
1254	ONLY: host
1255
1256	USE cpulog, &
1257	ONLY: cpu_log, log_point_s
1258
1259	USE kinds
1260
1261	USE pegrid
1262
1263	IMPLICIT NONE
1264
1265	INTEGER(iwp) :: i !<
1266	INTEGER(iwp) :: j !<
1267	INTEGER(iwp) :: k !<
1268
1269	REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn) :: work_fftx !<
1270	REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) :: f_in !<
1271	REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx) :: f_out !<
1272	REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx) :: work !<
1273
1274	!
1275	!-- Transpose array
1276	#if defined( __parallel )
1277	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1278	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1279	CALL MPI_ALLTOALL( f_in(1,1,nxl_y,1), sendrecvcount_xy, MPI_REAL, &
1280	work(nys,1,0), sendrecvcount_xy, MPI_REAL, &
1281	comm1dy, ierr )
1282	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1283	#endif
1284
1285	!
1286	!-- Carry out the FFT along x, where all data are present due to the
1287	!-- 1d-decomposition along y. Resort the data in a way that y becomes
1288	!-- the first index.
1289	CALL cpu_log( log_point_s(4), 'fft_x_1d', 'continue' )
1290
1291	IF ( host(1:3) == 'nec' ) THEN
1292	!
1293	!-- Code optimized for vector processors
1294	!$OMP PARALLEL PRIVATE ( i, j, k )
1295	!$OMP DO
1296	DO j = nys, nyn
1297
1298	DO k = 1, nz
1299	DO i = 0, nx
1300	work_fftx(i,k,j) = work(j,k,i)
1301	ENDDO
1302	ENDDO
1303
1304	CALL fft_x_m( work_fftx(:,:,j), 'backward' )
1305
1306	ENDDO
1307
1308	!$OMP DO
1309	DO i = 0, nx
1310	DO j = nys, nyn
1311	DO k = 1, nz
1312	f_out(k,j,i) = work_fftx(i,k,j)
1313	ENDDO
1314	ENDDO
1315	ENDDO
1316	!$OMP END PARALLEL
1317
1318	ELSE
1319
1320	!
1321	!-- Cache optimized code (there might be still a potential for better
1322	!-- optimization).
1323	!$OMP PARALLEL PRIVATE (i,j,k)
1324	!$OMP DO
1325	DO j = nys, nyn
1326	DO k = 1, nz
1327
1328	DO i = 0, nx
1329	work_fftx(i,k,j) = work(j,k,i)
1330	ENDDO
1331
1332	CALL fft_x_1d( work_fftx(0:nx,k,j), 'backward' )
1333
1334	ENDDO
1335	ENDDO
1336
1337	!$OMP DO
1338	DO i = 0, nx
1339	DO j = nys, nyn
1340	DO k = 1, nz
1341	f_out(k,j,i) = work_fftx(i,k,j)
1342	ENDDO
1343	ENDDO
1344	ENDDO
1345	!$OMP END PARALLEL
1346
1347	ENDIF
1348	CALL cpu_log( log_point_s(4), 'fft_x_1d', 'stop' )
1349
1350	END SUBROUTINE tr_yx_fftx
1351
1352
1353	!------------------------------------------------------------------------------!
1354	! Description:
1355	! ------------
1356	!> FFT along y, solution of the tridiagonal system and backward FFT for
1357	!> a 1d-decomposition along y.
1358	!>
1359	!> @warning this subroutine may still not work for hybrid parallelization
1360	!> with OpenMP (for possible necessary changes see the original
1361	!> routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002)
1362	!------------------------------------------------------------------------------!
1363	SUBROUTINE ffty_tri_ffty( ar )
1364
1365
1366	USE control_parameters, &
1367	ONLY: host
1368
1369	USE cpulog, &
1370	ONLY: cpu_log, log_point_s
1371
1372	USE grid_variables, &
1373	ONLY: ddx2, ddy2
1374
1375	USE kinds
1376
1377	USE pegrid
1378
1379	IMPLICIT NONE
1380
1381	INTEGER(iwp) :: i !<
1382	INTEGER(iwp) :: j !<
1383	INTEGER(iwp) :: k !<
1384	INTEGER(iwp) :: m !<
1385	INTEGER(iwp) :: n !<
1386	INTEGER(iwp) :: omp_get_thread_num !<
1387	INTEGER(iwp) :: tn !<
1388
1389	REAL(wp), DIMENSION(0:ny) :: work_ffty !<
1390	REAL(wp), DIMENSION(0:ny,1:nz) :: work_triy !<
1391	REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) :: ar !<
1392	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tri !<
1393
1394
1395	CALL cpu_log( log_point_s(39), 'fft_y_1d + tridia', 'start' )
1396
1397	ALLOCATE( tri(5,0:ny,0:nz-1,0:threads_per_task-1) )
1398
1399	tn = 0 ! Default thread number in case of one thread
1400	!$OMP PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_ffty, work_triy )
1401	DO i = nxl_y, nxr_y
1402
1403	!$ tn = omp_get_thread_num()
1404
1405	IF ( host(1:3) == 'nec' ) THEN
1406	!
1407	!-- Code optimized for vector processors
1408	DO k = 1, nz
1409
1410	m = 0
1411	DO n = 1, pdims(2)
1412	DO j = 1, nny
1413	work_triy(m,k) = ar(j,k,i,n)
1414	m = m + 1
1415	ENDDO
1416	ENDDO
1417
1418	ENDDO
1419
1420	CALL fft_y_m( work_triy, ny, 'forward' )
1421
1422	ELSE
1423	!
1424	!-- Cache optimized code
1425	DO k = 1, nz
1426
1427	m = 0
1428	DO n = 1, pdims(2)
1429	DO j = 1, nny
1430	work_ffty(m) = ar(j,k,i,n)
1431	m = m + 1
1432	ENDDO
1433	ENDDO
1434
1435	CALL fft_y_1d( work_ffty, 'forward' )
1436
1437	DO j = 0, ny
1438	work_triy(j,k) = work_ffty(j)
1439	ENDDO
1440
1441	ENDDO
1442
1443	ENDIF
1444
1445	!
1446	!-- Solve the linear equation system
1447	CALL tridia_1dd( ddy2, ddx2, ny, nx, i, work_triy, tri(:,:,:,tn) )
1448
1449	IF ( host(1:3) == 'nec' ) THEN
1450	!
1451	!-- Code optimized for vector processors
1452	CALL fft_y_m( work_triy, ny, 'backward' )
1453
1454	DO k = 1, nz
1455
1456	m = 0
1457	DO n = 1, pdims(2)
1458	DO j = 1, nny
1459	ar(j,k,i,n) = work_triy(m,k)
1460	m = m + 1
1461	ENDDO
1462	ENDDO
1463
1464	ENDDO
1465
1466	ELSE
1467	!
1468	!-- Cache optimized code
1469	DO k = 1, nz
1470
1471	DO j = 0, ny
1472	work_ffty(j) = work_triy(j,k)
1473	ENDDO
1474
1475	CALL fft_y_1d( work_ffty, 'backward' )
1476
1477	m = 0
1478	DO n = 1, pdims(2)
1479	DO j = 1, nny
1480	ar(j,k,i,n) = work_ffty(m)
1481	m = m + 1
1482	ENDDO
1483	ENDDO
1484
1485	ENDDO
1486
1487	ENDIF
1488
1489	ENDDO
1490
1491	DEALLOCATE( tri )
1492
1493	CALL cpu_log( log_point_s(39), 'fft_y_1d + tridia', 'stop' )
1494
1495	END SUBROUTINE ffty_tri_ffty
1496
1497	END MODULE poisfft_mod

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: palm/trunk/SOURCE/poisfft.f90 @ 1842

Download in other formats: