Line data Source code
1 : !!****m* ABINIT/m_sg2002
2 : !! NAME
3 : !! m_sg2002
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2002-2007 Stefan Goedecker, CEA Grenoble
9 : !! Copyright (C) 2014-2026 ABINIT group (XG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! NOTES
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : module m_sg2002
25 :
26 : use defs_basis
27 : use defs_fftdata
28 : use m_abicore
29 : use m_errors
30 : use m_xmpi
31 :
32 : use m_time, only : timab
33 : use m_fstrings, only : itoa
34 : use m_fftcore, only : sphere_fft1, fill, scramble, switchreal, switch, mpiswitch,&
35 : & unfill, unscramble, unswitchreal, unswitch, unmpiswitch,&
36 : & fill_cent, switch_cent, switchreal_cent, mpiswitch_cent, multpot, addrho,&
37 : & unfill_cent, unswitchreal_cent, unswitch_cent, unmpiswitch_cent, unscramble,&
38 : & mpifft_fg2dbox, mpifft_dbox2fr, mpifft_fr2dbox, mpifft_dbox2fg
39 :
40 : implicit none
41 :
42 : private
43 :
44 : ! Public API:
45 : !public :: sg2002_seqfourdp ! seq-FFT of densities and potentials.
46 : public :: sg2002_mpifourdp ! MPI-FFT of densities and potentials.
47 : !public :: mpi_fourwf
48 : !public :: sg2002_seqfourwf ! seq-FFT of wavefunctions.
49 : !public :: sg2002_mpifourwf ! MPI-FFT of wavefunctions.
50 :
51 : ! Low-level tools.
52 : ! These procedure shouls be accessed via a wrapper that selected the library via fftalg
53 : public :: sg2002_back ! G --> R for densities and potentials
54 : public :: sg2002_forw ! R --> G for densities and potentials
55 : public :: sg2002_mpiback_wf ! G --> R for wavefunctions
56 : public :: sg2002_mpiforw_wf ! R --> G for wavefunctions
57 : public :: sg2002_applypot ! Compute <G|vloc|u> where u is given in reciprocal space.
58 : public :: sg2002_applypot_many ! Compute <G|vloc|u> where u is given in reciprocal space.
59 : public :: sg2002_accrho ! Compute rho = weigth_r*Re(u(r))**2 + weigth_i*Im(u(r))**2
60 :
61 : contains
62 : !!***
63 :
64 : !!****f* m_sg2002/sg2002_back
65 : !! NAME
66 : !! sg2002_back
67 : !!
68 : !! FUNCTION
69 : !! CALCULATES THE DISCRETE FOURIER TRANSFORM in parallel using MPI/OpenMP
70 : !!
71 : !! ZR(I1,I2,I3)= \sum_(j1,j2,j3) EXP(isign*i*2*pi*(j1*i1/n1+j2*i2/n2+j3*i3/n3)) ZF(j1,j3,j2)
72 : !!
73 : !! Adopt standard convention that isign=1 for backward transform
74 : !!
75 : !! INPUTS:
76 : !! cplex=1 for real --> complex, 2 for complex --> complex
77 : !! ZF: input array in G-space (note the switch of i2 and i3)
78 : !!
79 : !! real(F(i1,i3,i2,idat))=ZF(1,i1,i3,i2,idat)
80 : !! imag(F(i1,i3,i2,idat))=ZF(2,i1,i3,i2,idat)
81 : !!
82 : !! i1=1,n1 , i2=1,n2 , i3=1,n3 , idat=1,ndat
83 : !! OUTPUTS:
84 : !! ZR: output array in R space.
85 : !!
86 : !! ZR(1,i1,i2,i3,idat)=real(R(i1,i2,i3,idat))
87 : !! ZR(2,i1,i2,i3,idat)=imag(R(i1,i2,i3,idat))
88 : !!
89 : !! i1=1,n1 , i2=1,n2 , i3=1,n3 , idat=1,ndat
90 : !!
91 : !! nproc_fft: number of processors used as returned by MPI_COMM_SIZE
92 : !! me_fft: [0:nproc_fft-1] number of processor as returned by MPI_COMM_RANK
93 : !! n1,n2,n3: logical dimension of the transform. As transform lengths
94 : !! most products of the prime factors 2,3,5 are allowed.
95 : !! The detailed table with allowed transform lengths can
96 : !! be found in subroutine CTRIG
97 : !! nd1,nd2,nd3: Dimension of ZF and ZR
98 : !! nd2proc=((nd2-1)/nproc_fft)+1 maximal number of 2nd dim slices
99 : !! nd3proc=((nd3-1)/nproc_fft)+1 maximal number of 3rd dim slices
100 : !!
101 : !! NOTES:
102 : !! The maximum number of processors that can reasonably be used is max(n2,n3)
103 : !! It is very important to find the optimal
104 : !! value of NCACHE. NCACHE determines the size of the work array ZW, that
105 : !! has to fit into cache. It has therefore to be chosen to equal roughly
106 : !! half the size of the physical cache in units of real*8 numbers.
107 : !! The optimal value of ncache can easily be determined by numerical
108 : !! experimentation. A too large value of ncache leads to a dramatic
109 : !! and sudden decrease of performance, a too small value to a to a
110 : !! slow and less dramatic decrease of performance. If NCACHE is set
111 : !! to a value so small, that not even a single one dimensional transform
112 : !! can be done in the workarray zw, the program stops with an error message.
113 : !!
114 : !! SOURCE
115 :
116 70446 : subroutine sg2002_back(cplex,ndat,n1,n2,n3,nd1,nd2,nd3,nd1eff,nd2proc,nd3proc,option,zf,zr,comm_fft)
117 :
118 : !Arguments ------------------------------------
119 : ! real space input
120 : integer,intent(in) :: cplex,ndat,n1,n2,n3,nd1,nd2,nd3,nd1eff,nd2proc,nd3proc,option,comm_fft
121 : real(dp),intent(in) :: zf(2,nd1,nd3,nd2proc,ndat)
122 : real(dp),intent(out) :: zr(2,nd1eff,nd2,nd3proc,ndat)
123 :
124 : !Local variables-------------------------------
125 : !scalars
126 : integer :: i,j,i1,ic1,ic2,ic3,idat,ierr,includelast,inzee,j2,j2st,j3,jeff,jp2st,lot,lzt
127 : integer :: ma,mb,n1dfft,n1eff,n2eff,n1zt,ncache,nnd3,nproc_fft,me_fft
128 : character(len=500) :: msg
129 : !arrays
130 70446 : real(dp), allocatable :: zt(:,:,:) ! work arrays for transpositions
131 70446 : real(dp), allocatable :: zmpi1(:,:,:,:),zmpi2(:,:,:,:) ! work arrays for MPI
132 70446 : real(dp), allocatable :: zw(:,:,:) ! cache work array
133 : real(dp) :: tsec(2)
134 : ! FFT work arrays
135 70446 : real(dp), allocatable, dimension(:,:) :: trig1,trig2,trig3
136 70446 : integer, allocatable, dimension(:) :: after1,now1,before1,after2,now2,before2,after3,now3,before3
137 :
138 : ! *************************************************************************
139 :
140 70446 : nproc_fft = xmpi_comm_size(comm_fft); me_fft = xmpi_comm_rank(comm_fft)
141 :
142 : ! find cache size that gives optimal performance on machine
143 70446 : ncache=4*max(n1,n2,n3,1024)
144 :
145 70446 : if (ncache/(4*max(n1,n2,n3))<1) then
146 : write(msg,'(5a)') &
147 0 : & 'ncache has to be enlarged to be able to hold at',ch10, &
148 0 : & 'least one 1-d FFT of each size even though this will',ch10,&
149 0 : & 'reduce the performance for shorter transform lengths'
150 0 : ABI_ERROR(msg)
151 : end if
152 :
153 : ! check input
154 70446 : if (nd1<n1 .or. nd2<n2 .or. nd3<n3) then
155 0 : ABI_ERROR("nd1<n1 .or. nd2<n2 .or. nd3<n3")
156 : end if
157 :
158 : ! Effective n1 and n2 (complex-to-complex or real-to-complex)
159 70446 : n1eff=n1; n2eff=n2; n1zt=n1
160 70446 : if (cplex==1) then
161 92 : n1eff=(n1+1)/2 ; n2eff=n2/2+1 ; n1zt=2*(n1/2+1)
162 : end if
163 :
164 70446 : lzt=n2eff
165 70446 : if (mod(n2eff,2) == 0) lzt=lzt+1
166 70446 : if (mod(n2eff,4) == 0) lzt=lzt+1
167 :
168 : ! maximal number of big box 3rd dim slices for all procs
169 70446 : nnd3=nd3proc*nproc_fft
170 :
171 211338 : ABI_MALLOC(trig1,(2,n1))
172 70446 : ABI_MALLOC(after1,(mdata))
173 70446 : ABI_MALLOC(now1,(mdata))
174 70446 : ABI_MALLOC(before1,(mdata))
175 211338 : ABI_MALLOC(trig2,(2,n2))
176 70446 : ABI_MALLOC(after2,(mdata))
177 70446 : ABI_MALLOC(now2,(mdata))
178 70446 : ABI_MALLOC(before2,(mdata))
179 211338 : ABI_MALLOC(trig3,(2,n3))
180 70446 : ABI_MALLOC(after3,(mdata))
181 70446 : ABI_MALLOC(now3,(mdata))
182 70446 : ABI_MALLOC(before3,(mdata))
183 281784 : ABI_MALLOC(zw,(2,ncache/4,2))
184 281784 : ABI_MALLOC(zt,(2,lzt,n1zt))
185 352230 : ABI_MALLOC(zmpi2,(2,n1,nd2proc,nnd3))
186 70446 : if (nproc_fft>1) then
187 54288 : ABI_MALLOC(zmpi1,(2,n1,nd2proc,nnd3))
188 : end if
189 :
190 70446 : call ctrig(n3,trig3,after3,before3,now3,1,ic3)
191 70446 : call ctrig(n1,trig1,after1,before1,now1,1,ic1)
192 70446 : call ctrig(n2,trig2,after2,before2,now2,1,ic2)
193 :
194 : !DEBUG
195 : ! write(std_out,'(a,3i4)' )'sg2002_back,zf n1,n2,n3',n1,n2,n3
196 : ! write(std_out,'(a,3i4)' )'nd1,nd2,nd3proc',nd1,nd2,nd3proc
197 : ! write(std_out,'(a,3i4)' )'m1,m2,m3',m1,m2,m3
198 : ! write(std_out,'(a,3i4)' )'max1,max2,max3',max1,max2,max3
199 : ! write(std_out,'(a,3i4)' )'md1,md2proc,md3',md1,md2proc,md3
200 : ! write(std_out,'(a,3i4)' )'n1eff,m2eff,m1zt',n1eff,m2eff,m1zt
201 : !ENDDEBUG
202 :
203 140900 : do idat=1,ndat
204 : ! transform along z axis
205 : ! input: I1,I3,J2,(Jp2)
206 70454 : lot=ncache/(4*n3)
207 :
208 985673 : do j2=1,nd2proc
209 985673 : if (me_fft*nd2proc+j2 <= n2eff) then
210 :
211 2078304 : do i1=1,n1,lot
212 1163090 : ma=i1
213 1163090 : mb=min(i1+(lot-1),n1)
214 1163090 : n1dfft=mb-ma+1
215 :
216 : ! input: G1,G3,G2,(Gp2)
217 1163090 : call fill(nd1,nd3,lot,n1dfft,n3,zf(1,i1,1,j2,idat),zw(1,1,1))
218 :
219 1163090 : inzee=1
220 3887798 : do i=1,ic3
221 : call fftstp(lot,n1dfft,n3,lot,n3,zw(1,1,inzee),zw(1,1,3-inzee), &
222 2724708 : & trig3,after3(i),now3(i),before3(i),1)
223 3887798 : inzee=3-inzee
224 : end do
225 :
226 : ! input: G1,R3,G2,(Gp2)
227 : ! output: G1,G2,R3,(Gp2)
228 2078304 : call scramble(i1,j2,lot,n1dfft,n1,n3,nd2proc,nd3,zw(1,1,inzee),zmpi2)
229 : end do
230 : end if
231 : end do
232 :
233 : ! Interprocessor data transposition
234 : ! input: G1,G2,R3,Rp3,(Gp2)
235 : ! output: G1,G2,G3,Gp2,(Rp3)
236 70454 : if (nproc_fft>1) then
237 13572 : call timab(543,1,tsec)
238 : call xmpi_alltoall(zmpi2,2*n1*nd2proc*nd3proc, &
239 13572 : & zmpi1,2*n1*nd2proc*nd3proc,comm_fft,ierr)
240 13572 : call timab(543,2,tsec)
241 : end if
242 :
243 1061883 : do j3=1,nd3proc
244 991437 : if (me_fft*nd3proc+j3 <= n3) then
245 920983 : Jp2st=1
246 920983 : J2st=1
247 :
248 : ! transform along x axis
249 920983 : lot=ncache/(4*n1)
250 :
251 2086074 : do j=1,n2eff,lot
252 1165091 : ma=j
253 1165091 : mb=min(j+(lot-1),n2eff)
254 1165091 : n1dfft=mb-ma+1
255 :
256 : ! input: G1,G2,R3,Gp2,(Rp3)
257 : ! output: G2,G1,R3,Jp2,(Rp3)
258 1165091 : if (nproc_fft == 1) then
259 865475 : call mpiswitch(j3,n1dfft,Jp2st,J2st,lot,n1,nd2proc,nd3proc,nproc_fft,option,zmpi2,zw(1,1,1))
260 : else
261 299616 : call mpiswitch(j3,n1dfft,Jp2st,J2st,lot,n1,nd2proc,nd3proc,nproc_fft,option,zmpi1,zw(1,1,1))
262 : end if
263 :
264 : ! input: G2,G1,R3,(Rp3)
265 : ! output: G2,R1,R3,(Rp3)
266 1165091 : inzee=1
267 2708708 : do i=1,ic1-1
268 : call fftstp(lot,n1dfft,n1,lot,n1,zw(1,1,inzee),zw(1,1,3-inzee), &
269 1543617 : & trig1,after1(i),now1(i),before1(i),1)
270 2708708 : inzee=3-inzee
271 : end do
272 :
273 1165091 : i=ic1
274 : call fftstp(lot,n1dfft,n1,lzt,n1zt,zw(1,1,inzee),zt(1,j,1), &
275 2086074 : & trig1,after1(i),now1(i),before1(i),1)
276 : end do
277 :
278 : ! transform along y axis
279 920983 : lot=ncache/(4*n2)
280 :
281 2086074 : do j=1,n1eff,lot
282 1165091 : ma=j
283 1165091 : mb=min(j+(lot-1),n1eff)
284 1165091 : n1dfft=mb-ma+1
285 1165091 : includelast=1
286 :
287 1165091 : if (cplex==1) then
288 6480 : jeff=2*j-1
289 : includelast=1
290 6480 : if (mb==n1eff .and. n1eff*2/=n1) includelast=0
291 : end if
292 :
293 : ! input: G2,R1,R3,(Rp3)
294 : ! output: R1,G2,R3,(Rp3)
295 1165091 : if (cplex==2) then
296 1158611 : call switch(n1dfft,n2,lot,n1,lzt,zt(1,1,j),zw(1,1,1))
297 : else
298 6480 : call switchreal(includelast,n1dfft,n2,n2eff,lot,n1zt,lzt,zt(1,1,jeff),zw(1,1,1))
299 : end if
300 :
301 1165091 : inzee=1
302 2708468 : do i=1,ic2-1
303 : call fftstp(lot,n1dfft,n2,lot,n2,zw(1,1,inzee),zw(1,1,3-inzee), &
304 1543377 : & trig2,after2(i),now2(i),before2(i),1)
305 2708468 : inzee=3-inzee
306 : end do
307 :
308 1165091 : i=ic2
309 : call fftstp(lot,n1dfft,n2,nd1eff,nd2,zw(1,1,inzee),zr(1,j,1,j3,idat), &
310 2086074 : & trig2,after2(i),now2(i),before2(i),1)
311 :
312 : end do
313 : ! output: R1,R2,R3,(Rp3)
314 :
315 : end if
316 : end do
317 : end do ! idat
318 :
319 70446 : ABI_FREE(trig1)
320 70446 : ABI_FREE(after1)
321 70446 : ABI_FREE(now1)
322 70446 : ABI_FREE(before1)
323 70446 : ABI_FREE(trig2)
324 70446 : ABI_FREE(after2)
325 70446 : ABI_FREE(now2)
326 70446 : ABI_FREE(before2)
327 70446 : ABI_FREE(trig3)
328 70446 : ABI_FREE(after3)
329 70446 : ABI_FREE(now3)
330 70446 : ABI_FREE(before3)
331 70446 : ABI_FREE(zmpi2)
332 70446 : ABI_FREE(zw)
333 70446 : ABI_FREE(zt)
334 70446 : if (nproc_fft>1) then
335 13572 : ABI_FREE(zmpi1)
336 : end if
337 :
338 70446 : end subroutine sg2002_back
339 : !!***
340 :
341 : !----------------------------------------------------------------------
342 :
343 : !!****f* m_sg2002/sg2002_forw
344 : !! NAME
345 : !! sg2002_forw
346 : !!
347 : !! FUNCTION
348 : !! Adopt standard convention that isign=-1 for forward transform
349 : !! CALCULATES THE DISCRETE FOURIERTRANSFORM ZF(I1,I3,I2)=
350 : !! S_(j1,j2,j3) EXP(isign*i*2*pi*(j1*i1/n1+j2*i2/n2+j3*i3/n3)) ZR(j1,j2,j3)
351 : !! in parallel using MPI/OpenMP and BLAS library calls.
352 : !!
353 : !! INPUTS
354 : !! ZR: input array
355 : !! ZR(1,i1,i2,i3,idat)=real(R(i1,i2,i3,idat))
356 : !! ZR(2,i1,i2,i3,idat)=imag(R(i1,i2,i3,idat))
357 : !! i1=1,n1 , i2=1,n2 , i3=1,n3 , idat=1,ndat
358 : !! OUTPUTS
359 : !! ZF: output array (note the switch of i2 and i3)
360 : !! real(F(i1,i3,i2,idat))=ZF(1,i1,i3,i2,idat)
361 : !! imag(F(i1,i3,i2,idat))=ZF(2,i1,i3,i2,idat)
362 : !! i1=1,n1 , i2=1,n2 , i3=1,n3 , idat=1,ndat
363 : !! nproc_fft: number of processors used as returned by MPI_COMM_SIZE
364 : !! me_fft: [0:nproc_fft-1] number of processor as returned by MPI_COMM_RANK
365 : !! n1,n2,n3: logical dimension of the transform. As transform lengths
366 : !! most products of the prime factors 2,3,5 are allowed.
367 : !! The detailed table with allowed transform lengths can
368 : !! be found in subroutine CTRIG
369 : !! nd1,nd2,nd3: Dimension of ZR and ZF
370 : !! nd2proc=((nd2-1)/nproc_fft)+1 maximal number of 2nd dim slices
371 : !! nd3proc=((nd3-1)/nproc_fft)+1 maximal number of 3rd dim slices
372 : !!
373 : !! NOTES
374 : !! SHOULD describe nd1eff
375 : !! SHOULD put cplex and nd1eff in OMP declarations
376 : !! SHOULD describe the change of value of nd2prod
377 : !!
378 : !! The maximum number of processors that can reasonably be used is max(n2,n3)
379 : !!
380 : !! It is very important to find the optimal
381 : !! value of NCACHE. NCACHE determines the size of the work array ZW, that
382 : !! has to fit into cache. It has therefore to be chosen to equal roughly
383 : !! half the size of the physical cache in units of real*8 numbers.
384 : !! The optimal value of ncache can easily be determined by numerical
385 : !! experimentation. A too large value of ncache leads to a dramatic
386 : !! and sudden decrease of performance, a too small value to a to a
387 : !! slow and less dramatic decrease of performance. If NCACHE is set
388 : !! to a value so small, that not even a single one dimensional transform
389 : !! can be done in the workarray zw, the program stops with an error message.
390 : !!
391 : !! SOURCE
392 :
393 69867 : subroutine sg2002_forw(cplex,ndat,n1,n2,n3,nd1,nd2,nd3,nd1eff,nd2proc,nd3proc,option,zr,zf,comm_fft)
394 :
395 : !Arguments ------------------------------------
396 : !scalars
397 : integer,intent(in) :: cplex,comm_fft
398 : integer,intent(in) :: ndat,n1,n2,n3,nd1,nd2,nd3,nd1eff,nd2proc,nd3proc,option
399 : !arrays
400 : real(dp),intent(in) :: zr(2,nd1eff,nd2,nd3proc,ndat)
401 : real(dp),intent(out) :: zf(2,nd1,nd3,nd2proc,ndat)
402 :
403 : !Local variables-------------------------------
404 : !scalars
405 : integer :: i,j,i1,ic1,ic2,ic3,idat,ierr,inzee,j2,j2st,j3,jp2st,lot,lzt
406 : integer :: ma,mb,n1dfft,n1eff,n2eff,n1zt,ncache,nnd3,nproc_fft,me_fft
407 : character(len=500) :: msg
408 : !arrays
409 69867 : real(dp), allocatable :: zt(:,:,:) ! work arrays for transpositions
410 69867 : real(dp), allocatable :: zmpi1(:,:,:,:),zmpi2(:,:,:,:) ! work arrays for MPI
411 69867 : real(dp), allocatable :: zw(:,:,:) ! cache work array
412 : real(dp) :: tsec(2)
413 : ! FFT work arrays
414 69867 : real(dp), allocatable, dimension(:,:) :: trig1,trig2,trig3
415 69867 : integer, allocatable, dimension(:) :: after1,now1,before1,after2,now2,before2,after3,now3,before3
416 :
417 : ! *************************************************************************
418 :
419 69867 : nproc_fft = xmpi_comm_size(comm_fft); me_fft = xmpi_comm_rank(comm_fft)
420 :
421 : ! find cache size that gives optimal performance on machine
422 69867 : ncache=4*max(n1,n2,n3,1024)
423 69867 : if (ncache/(4*max(n1,n2,n3))<1) then
424 : write(msg,'(5a)')&
425 0 : & 'ncache has to be enlarged to be able to hold at',ch10, &
426 0 : & 'least one 1-d FFT of each size even though this will',ch10,&
427 0 : & 'reduce the performance for shorter transform lengths'
428 0 : ABI_ERROR(msg)
429 : end if
430 :
431 : ! check input
432 69867 : if (nd1<n1 .or. nd2<n2 .or. nd3<n3) then
433 0 : ABI_ERROR("nd1<n1 .or. nd2<n2 .or. nd3<n3")
434 : end if
435 :
436 : !Effective n1 and n2 (complex-to-complex or real-to-complex)
437 69867 : n1eff=n1; n2eff=n2; n1zt=n1
438 69867 : if (cplex==1) then
439 80 : n1eff=(n1+1)/2; n2eff=n2/2+1; n1zt=2*(n1/2+1)
440 : end if
441 :
442 69867 : lzt=n2eff
443 69867 : if (mod(n2eff,2) == 0) lzt=lzt+1
444 69867 : if (mod(n2eff,4) == 0) lzt=lzt+1
445 :
446 : ! maximal number of big box 3rd dim slices for all procs
447 69867 : nnd3=nd3proc*nproc_fft
448 :
449 209601 : ABI_MALLOC(trig1,(2,n1))
450 69867 : ABI_MALLOC(after1,(mdata))
451 69867 : ABI_MALLOC(now1,(mdata))
452 69867 : ABI_MALLOC(before1,(mdata))
453 209601 : ABI_MALLOC(trig2,(2,n2))
454 69867 : ABI_MALLOC(after2,(mdata))
455 69867 : ABI_MALLOC(now2,(mdata))
456 69867 : ABI_MALLOC(before2,(mdata))
457 209601 : ABI_MALLOC(trig3,(2,n3))
458 69867 : ABI_MALLOC(after3,(mdata))
459 69867 : ABI_MALLOC(now3,(mdata))
460 69867 : ABI_MALLOC(before3,(mdata))
461 279468 : ABI_MALLOC(zw,(2,ncache/4,2))
462 279468 : ABI_MALLOC(zt,(2,lzt,n1zt))
463 349335 : ABI_MALLOC(zmpi2,(2,n1,nd2proc,nnd3))
464 69867 : if (nproc_fft>1) then
465 53792 : ABI_MALLOC(zmpi1,(2,n1,nd2proc,nnd3))
466 : end if
467 :
468 69867 : call ctrig(n2,trig2,after2,before2,now2,-1,ic2)
469 69867 : call ctrig(n1,trig1,after1,before1,now1,-1,ic1)
470 69867 : call ctrig(n3,trig3,after3,before3,now3,-1,ic3)
471 :
472 139742 : do idat=1,ndat
473 1004478 : do j3=1,nd3proc
474 1004478 : if (me_fft*(nd3proc)+j3 <= n3) then
475 932830 : Jp2st=1; J2st=1
476 :
477 : ! transform along y axis
478 : ! input: R1,R2,R3,(Rp3)
479 932830 : lot=ncache/(4*n2)
480 :
481 2140692 : do j=1,n1eff,lot
482 1207862 : ma=j
483 1207862 : mb=min(j+(lot-1),n1eff)
484 1207862 : n1dfft=mb-ma+1
485 1207862 : i=1
486 : call fftstp(nd1eff,n1dfft,nd2,lot,n2,zr(1,j,1,j3,idat),zw(1,1,1), &
487 1207862 : & trig2,after2(i),now2(i),before2(i),-1)
488 :
489 1207862 : inzee=1
490 2876183 : do i=2,ic2
491 : call fftstp(lot,n1dfft,n2,lot,n2,zw(1,1,inzee),zw(1,1,3-inzee), &
492 1668321 : & trig2,after2(i),now2(i),before2(i),-1)
493 2876183 : inzee=3-inzee
494 : end do
495 :
496 : ! input: R1,G2,R3,(Rp3)
497 : ! output: G2,R1,R3,(Rp3)
498 2140692 : if(cplex==2)then
499 1201940 : call unswitch(n1dfft,n2,lot,n1zt,lzt,zw(1,1,inzee),zt(1,1,j))
500 : else
501 5922 : call unswitchreal(n1dfft,n2,n2eff,lot,n1zt,lzt,zw(1,1,inzee),zt(1,1,2*j-1))
502 : end if
503 : end do
504 :
505 : ! transform along x axis
506 : ! input: G2,R1,R3,(Rp3)
507 932830 : lot=ncache/(4*n1)
508 :
509 2140692 : do j=1,n2eff,lot
510 1207862 : ma=j
511 1207862 : mb=min(j+(lot-1),n2eff)
512 1207862 : n1dfft=mb-ma+1
513 :
514 1207862 : i=1
515 : call fftstp(lzt,n1dfft,n1zt,lot,n1,zt(1,j,1),zw(1,1,1), &
516 1207862 : & trig1,after1(i),now1(i),before1(i),-1)
517 :
518 1207862 : inzee=1
519 2876423 : do i=2,ic1
520 : call fftstp(lot,n1dfft,n1,lot,n1,zw(1,1,inzee),zw(1,1,3-inzee), &
521 1668561 : & trig1,after1(i),now1(i),before1(i),-1)
522 2876423 : inzee=3-inzee
523 : end do
524 : ! output: G2,G1,R3,(Rp3)
525 :
526 : ! input: G2,G1,R3,Gp2,(Rp3)
527 : ! output: G1,G2,R3,Gp2,(Rp3)
528 : ! write(std_out,*) 'J2st,Jp2st',J2st,Jp2st
529 2140692 : if (nproc_fft == 1) then
530 851618 : call unmpiswitch(j3,n1dfft,Jp2st,J2st,lot,n1,nd2proc,nd3proc,nproc_fft,option,zw(1,1,inzee),zmpi2)
531 : else
532 356244 : call unmpiswitch(j3,n1dfft,Jp2st,J2st,lot,n1,nd2proc,nd3proc,nproc_fft,option,zw(1,1,inzee),zmpi1)
533 : end if
534 : end do
535 :
536 : end if
537 : end do ! j3
538 :
539 : ! Interprocessor data transposition
540 : ! input: G1,G2,R3,Gp2,(Rp3)
541 : ! output: G1,G2,R3,Rp3,(Gp2)
542 69875 : if (nproc_fft>1) then
543 13448 : call timab(544,1,tsec)
544 : call xmpi_alltoall(zmpi1,2*n1*nd2proc*nd3proc, &
545 13448 : & zmpi2,2*n1*nd2proc*nd3proc,comm_fft,ierr)
546 13448 : call timab(544,2,tsec)
547 : end if
548 :
549 : ! transform along z axis
550 : ! input: G1,G2,R3,(Gp2)
551 69875 : lot=ncache/(4*n3)
552 :
553 1067578 : do j2=1,nd2proc
554 997711 : if (me_fft*(nd2proc)+j2 <= n2eff) then
555 2133232 : do i1=1,n1,lot
556 1205401 : ma=i1
557 1205401 : mb=min(i1+(lot-1),n1)
558 1205401 : n1dfft=mb-ma+1
559 :
560 : ! input: G1,G2,R3,(Gp2)
561 : ! output: G1,R3,G2,(Gp2)
562 1205401 : call unscramble(i1,j2,lot,n1dfft,n1,n3,nd2proc,nd3,zmpi2,zw(1,1,1))
563 :
564 1205401 : inzee=1
565 4094619 : do i=1,ic3
566 : call fftstp(lot,n1dfft,n3,lot,n3,zw(1,1,inzee),zw(1,1,3-inzee), &
567 2889218 : & trig3,after3(i),now3(i),before3(i),-1)
568 4094619 : inzee=3-inzee
569 : end do
570 :
571 2133232 : call unfill(nd1,nd3,lot,n1dfft,n3,zw(1,1,inzee),zf(1,i1,1,j2,idat))
572 : ! output: G1,G3,G2,(Gp2)
573 : end do
574 : end if
575 : end do
576 :
577 : end do ! idat
578 :
579 69867 : ABI_FREE(trig1)
580 69867 : ABI_FREE(after1)
581 69867 : ABI_FREE(now1)
582 69867 : ABI_FREE(before1)
583 69867 : ABI_FREE(trig2)
584 69867 : ABI_FREE(after2)
585 69867 : ABI_FREE(now2)
586 69867 : ABI_FREE(before2)
587 69867 : ABI_FREE(trig3)
588 69867 : ABI_FREE(after3)
589 69867 : ABI_FREE(now3)
590 69867 : ABI_FREE(before3)
591 69867 : ABI_FREE(zmpi2)
592 69867 : ABI_FREE(zw)
593 69867 : ABI_FREE(zt)
594 69867 : if (nproc_fft>1) then
595 13448 : ABI_FREE(zmpi1)
596 : end if
597 :
598 69867 : end subroutine sg2002_forw
599 : !!***
600 :
601 : !----------------------------------------------------------------------
602 :
603 : !!****f* m_sg2002/sg2002_mpiback_wf
604 : !! NAME
605 : !! sg2002_mpiback_wf
606 : !!
607 : !! FUNCTION
608 : !! Does multiple 3-dim backward FFTs from Fourier into real space
609 : !! Adopt standard convention that isign=1 for backward transform
610 : !!
611 : !! CALCULATES THE DISCRETE FOURIER TRANSFORM ZF(I1,I2,I3)=
612 : !!
613 : !! S_(j1,j2,j3) EXP(isign*i*2*pi*(j1*i1/n1+j2*i2/n2+j3*i3/n3)) ZF(j1,j3,j2)
614 : !!
615 : !! in parallel using MPI/OpenMP.
616 : !!
617 : !! INPUTS:
618 : !! icplexwf=1 if wavefunction is real, 2 if complex
619 : !! ndat=Number of wavefunctions to transform.
620 : !! n1,n2,n3: logical dimension of the transform. As transform lengths
621 : !! most products of the prime factors 2,3,5 are allowed.
622 : !! The detailed table with allowed transform lengths can be found in subroutine CTRIG
623 : !! nd1,nd2,nd3: Leading Dimension of ZR
624 : !! nd3proc=((nd3-1)/nproc_fft)+1 maximal number of big box 3rd dim slices for one proc
625 : !! max1 is positive or zero; m1 >=max1+1
626 : !! i1= 1... max1+1 corresponds to positive and zero wavevectors 0 ... max1
627 : !! then, if m1 > max1+1, one has min1=max1-m1+1 and
628 : !! i1= max1+2 ... m1 corresponds to negative wavevectors min1 ... -1
629 : !! max2 and max3 have a similar definition of range
630 : !! m1,m2,m3=Size of the box enclosing the G-sphere.
631 : !! md1,md2,md3: Dimension of ZF given on the **small** FFT box.
632 : !! md2proc=((md2-1)/nproc_fft)+1 maximal number of small box 2nd dim slices for one proc
633 : !! nproc_fft: number of processors used as returned by MPI_COMM_SIZE
634 : !! comm_fft=MPI communicator for the FFT.
635 : !! ZF: input array (note the switch of i2 and i3)
636 : !! real(F(i1,i3,i2,idat))=ZF(1,i1,i3,i2,idat)
637 : !! imag(F(i1,i3,i2,idat))=ZF(2,i1,i3,i2,idat)
638 : !!
639 : !! OUTPUTS
640 : !! ZR: output array
641 : !! ZR(1,i1,i2,i3,idat)=real(R(i1,i2,i3,idat))
642 : !! ZR(2,i1,i2,i3,idat)=imag(R(i1,i2,i3,idat))
643 : !! i1=1,n1 , i2=1,n2 , i3=1,n3 , idat=1,ndat
644 : !!
645 : !! NOTES
646 : !! The maximum number of processors that can reasonably be used is max(n2/2,n3/2)
647 : !!
648 : !! It is very important to find the optimal
649 : !! value of NCACHE. NCACHE determines the size of the work array ZW, that
650 : !! has to fit into cache. It has therefore to be chosen to equal roughly
651 : !! half the size of the physical cache in units of real*8 numbers.
652 : !! The optimal value of ncache can easily be determined by numerical
653 : !! experimentation. A too large value of ncache leads to a dramatic
654 : !! and sudden decrease of performance, a too small value to a to a
655 : !! slow and less dramatic decrease of performance. If NCACHE is set
656 : !! to a value so small, that not even a single one dimensional transform
657 : !! can be done in the workarray zw, the program stops with an error message.
658 : !!
659 : !! SOURCE
660 :
661 798252 : subroutine sg2002_mpiback_wf(icplexwf,ndat,n1,n2,n3,nd1,nd2,nd3proc,&
662 798252 : & max1,max2,max3,m1,m2,m3,md1,md2proc,md3,zf,zr,comm_fft)
663 :
664 : !Arguments ------------------------------------
665 : integer,intent(in) :: icplexwf,ndat,n1,n2,n3,nd1,nd2,nd3proc
666 : integer,intent(in) :: max1,max2,max3,m1,m2,m3,md1,md2proc,md3,comm_fft
667 : real(dp),intent(in) :: zf(2,md1,md3,md2proc,ndat)
668 : real(dp),intent(out) :: zr(2,nd1,nd2,nd3proc,ndat)
669 :
670 : !Local variables-------------------------------
671 : integer :: i,j,i1,i2,ic1,ic2,ic3,idat,ierr,inzee,includelast
672 : integer :: ioption,j2,j3,j2st,jp2st,jeff,lot,lzt,m1zt,ma,mb,n1dfft,nnd3
673 : integer :: m2eff,ncache,n1eff,n1half,nproc_fft,me_fft
674 : character(len=500) :: msg
675 : !arrays
676 798252 : real(dp),allocatable :: zt(:,:,:) ! work arrays for transpositions
677 798252 : real(dp),allocatable :: zmpi1(:,:,:,:,:),zmpi2(:,:,:,:,:) ! work arrays for MPI
678 798252 : real(dp),allocatable :: zw(:,:,:) ! cache work array
679 : ! FFT work arrays
680 798252 : real(dp),allocatable :: trig1(:,:),trig2(:,:),trig3(:,:)
681 798252 : integer,allocatable :: after1(:),now1(:),before1(:),after2(:)
682 798252 : integer,allocatable :: now2(:),before2(:),after3(:),now3(:),before3(:)
683 : real(dp) :: tsec(2)
684 :
685 : ! *************************************************************************
686 :
687 : ! call timab(541,1,tsec)
688 : ! FIXME must provide a default value but which one?
689 : ! ioption = 0
690 798252 : ioption = 1
691 : !if (paral_kgb==1) ioption=1
692 :
693 798252 : nproc_fft = xmpi_comm_size(comm_fft); me_fft = xmpi_comm_rank(comm_fft)
694 :
695 : ! Find cache size that gives optimal performance on machine
696 798252 : ncache=4*max(n1,n2,n3,1024)
697 798252 : if (ncache/(4*max(n1,n2,n3))<1) then
698 : write(msg,"(5a)") &
699 0 : & 'ncache has to be enlarged to be able to hold at',ch10, &
700 0 : & 'least one 1-d FFT of each size even though this will',ch10,&
701 0 : & 'reduce the performance for shorter transform lengths'
702 0 : ABI_ERROR(msg)
703 : end if
704 :
705 : ! Effective m1 and m2 (complex-to-complex or real-to-complex)
706 798252 : n1eff=n1; m2eff=m2; m1zt=n1
707 798252 : if (icplexwf==1) then
708 51 : n1eff=(n1+1)/2; m2eff=m2/2+1; m1zt=2*(n1/2+1)
709 : end if
710 :
711 798252 : lzt=m2eff
712 798252 : if (mod(m2eff,2)==0) lzt=lzt+1
713 701191 : if (mod(m2eff,4)==0) lzt=lzt+1
714 :
715 : ! maximal number of big box 3rd dim slices for all procs
716 798252 : nnd3=nd3proc*nproc_fft
717 :
718 2394756 : ABI_MALLOC(trig1,(2,n1))
719 798252 : ABI_MALLOC(after1,(mdata))
720 798252 : ABI_MALLOC(now1,(mdata))
721 798252 : ABI_MALLOC(before1,(mdata))
722 2394756 : ABI_MALLOC(trig2,(2,n2))
723 798252 : ABI_MALLOC(after2,(mdata))
724 798252 : ABI_MALLOC(now2,(mdata))
725 798252 : ABI_MALLOC(before2,(mdata))
726 2394756 : ABI_MALLOC(trig3,(2,n3))
727 798252 : ABI_MALLOC(after3,(mdata))
728 798252 : ABI_MALLOC(now3,(mdata))
729 798252 : ABI_MALLOC(before3,(mdata))
730 :
731 : ! Allocate cache work array and work arrays for MPI transpositions.
732 3193008 : ABI_MALLOC(zw,(2,ncache/4,2))
733 3193008 : ABI_MALLOC(zt,(2,lzt,m1zt))
734 4789512 : ABI_MALLOC(zmpi2,(2,md1,md2proc,nnd3,ndat))
735 798252 : if (nproc_fft>1) then
736 3482150 : ABI_MALLOC(zmpi1,(2,md1,md2proc,nnd3,ndat))
737 : end if
738 :
739 : ! Compute twiddle coefficients.
740 798252 : call ctrig(n3,trig3,after3,before3,now3,1,ic3)
741 798252 : call ctrig(n1,trig1,after1,before1,now1,1,ic1)
742 798252 : call ctrig(n2,trig2,after2,before2,now2,1,ic2)
743 :
744 : !DEBUG
745 : ! write(std_out,'(2a,3i4)' )itoa(me_fft),': sg2002_mpiback_wf,zf n1,n2,n3',n1,n2,n3
746 : ! write(std_out,'(2a,3i4)' )itoa(me_fft),': nd1,nd2,nd3proc',nd1,nd2,nd3proc
747 : ! write(std_out,'(2a,3i4)' )itoa(me_fft),': m1,m2,m3',m1,m2,m3
748 : ! write(std_out,'(2a,3i4)' )itoa(me_fft),': max1,max2,max3',max1,max2,max3
749 : ! write(std_out,'(2a,3i4)' )itoa(me_fft),': md1,md2proc,md3',md1,md2proc,md3
750 : ! write(std_out,'(2a,3i4)' )itoa(me_fft),'n1eff,m2eff,m1zt',n1eff,m2eff,m1zt
751 : !ENDDEBUG
752 :
753 1647045 : do idat=1,ndat
754 :
755 : ! transform along z axis
756 : ! input: G1,G3,G2,(Gp2)
757 848793 : lot=ncache/(4*n3)
758 :
759 5217530571 : zw(:,:,:)=zero
760 629165709 : zt(:,:,:)=zero
761 :
762 : ! Loop over the y planes treated by this node and trasform n1ddft G_z lines.
763 5456324 : do j2=1,md2proc
764 :
765 : ! if (me_fft*md2proc+j2<=m2eff) then !a faire plus tard
766 :
767 10229324 : do i1=1,m1,lot
768 4773000 : ma=i1
769 4773000 : mb=min(i1+(lot-1),m1)
770 4773000 : n1dfft=mb-ma+1
771 :
772 : ! zero-pad n1dfft G_z lines
773 : ! input: G1,G3,G2,(Gp2)
774 : ! output: G1,R3,G2,(Gp2)
775 4773000 : call fill_cent(md1,md3,lot,n1dfft,max3,m3,n3,zf(1,i1,1,j2,idat),zw(1,1,1))
776 :
777 : ! Transform along z.
778 4773000 : inzee=1
779 14617403 : do i=1,ic3
780 : call fftstp(lot,n1dfft,n3,lot,n3,zw(1,1,inzee),zw(1,1,3-inzee), &
781 9844403 : & trig3,after3(i),now3(i),before3(i),1)
782 14617403 : inzee=3-inzee
783 : end do
784 :
785 : ! Local rotation.
786 : ! input: G1,R3,G2,(Gp2)
787 : ! output: G1,G2,R3,(Gp2)
788 9380531 : call scramble(i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3,zw(1,1,inzee),zmpi2(:,:,:,:,idat))
789 : end do
790 : !
791 : end do ! j2
792 :
793 : ! Interprocessor data transposition
794 : ! input: G1,G2,R3,Rp3,(Gp2)
795 : ! output: G1,G2,R3,Gp2,(Rp3)
796 848793 : if (nproc_fft>1) then
797 715074 : call timab(543,1,tsec)
798 : call xmpi_alltoall(zmpi2(:,:,:,:,idat),2*md1*md2proc*nd3proc, &
799 715074 : & zmpi1(:,:,:,:,idat),2*md1*md2proc*nd3proc,comm_fft,ierr)
800 715074 : call timab(543,2,tsec)
801 : end if
802 :
803 : ! Loop over the z treated by this node.
804 11218134 : do j3=1,nd3proc
805 : !j3glob = j3 + me_fft*nd3proc
806 10419882 : if (me_fft*nd3proc+j3 <= n3) then
807 9571089 : Jp2st=1; J2st=1
808 :
809 9571089 : lot=ncache/(4*n1)
810 :
811 : ! Loop over G_y in the small box.
812 19478978 : do j=1,m2eff,lot
813 9907889 : ma=j
814 9907889 : mb=min(j+(lot-1),m2eff)
815 9907889 : n1dfft=mb-ma+1
816 :
817 : ! Zero-pad input.
818 : ! input: G1,G2,R3,JG2,(Rp3)
819 : ! output: G2,G1,R3,JG2,(Rp3)
820 9907889 : if (nproc_fft==1) then
821 : call mpiswitch_cent(j3,n1dfft,Jp2st,J2st,lot,max1,md1,m1,n1,&
822 2628931 : & md2proc,nd3proc,nproc_fft,ioption,zmpi2(:,:,:,:,idat),zw(1,1,1),max2,m2,n2)
823 : else
824 : call mpiswitch_cent(j3,n1dfft,Jp2st,J2st,lot,max1,md1,m1,n1,&
825 7278958 : & md2proc,nd3proc,nproc_fft,ioption,zmpi1(:,:,:,:,idat),zw(1,1,1),max2,m2,n2)
826 : end if
827 :
828 : ! Transform along x
829 : ! input: G2,G1,R3,(Rp3)
830 : ! output: G2,R1,R3,(Rp3)
831 9907889 : inzee=1
832 20672388 : do i=1,ic1-1
833 : call fftstp(lot,n1dfft,n1,lot,n1,zw(1,1,inzee),zw(1,1,3-inzee), &
834 10764499 : & trig1,after1(i),now1(i),before1(i),1)
835 20672388 : inzee=3-inzee
836 : end do
837 :
838 9907889 : i=ic1
839 : call fftstp(lot,n1dfft,n1,lzt,m1zt,zw(1,1,inzee),zt(1,j,1), &
840 19478978 : & trig1,after1(i),now1(i),before1(i),1)
841 : end do
842 :
843 : ! Transform along y axis (take into account c2c or c2r case).
844 : ! Must loop over the full box.
845 9571089 : lot=ncache/(4*n2)
846 :
847 19830678 : do j=1,n1eff,lot
848 10259589 : ma=j
849 10259589 : mb=min(j+(lot-1),n1eff)
850 10259589 : n1dfft=mb-ma+1
851 10259589 : includelast=1
852 10259589 : if (icplexwf==1) then
853 15082 : jeff=2*j-1
854 15082 : if (mb==n1eff .and. n1eff*2/=n1) includelast=0
855 : end if
856 :
857 : ! Zero-pad the input.
858 : ! input: G2,R1,R3,(Rp3)
859 : ! output: R1,G2,R3,(Rp3)
860 10259589 : if (icplexwf==2) then
861 10244507 : call switch_cent(n1dfft,max2,m2,n2,lot,n1,lzt,zt(1,1,j),zw(1,1,1))
862 : else
863 15082 : call switchreal_cent(includelast,n1dfft,max2,n2,lot,m1zt,lzt,zt(1,1,jeff),zw(1,1,1))
864 : end if
865 :
866 : ! input: R1,G2,R3,(Rp3)
867 : ! output: R1,R2,R3,(Rp3)
868 10259589 : inzee=1
869 21485778 : do i=1,ic2-1
870 : call fftstp(lot,n1dfft,n2,lot,n2,zw(1,1,inzee),zw(1,1,3-inzee), &
871 11226189 : & trig2,after2(i),now2(i),before2(i),1)
872 21485778 : inzee=3-inzee
873 : end do
874 :
875 10259589 : i=ic2
876 :
877 : call fftstp(lot,n1dfft,n2,nd1,nd2,zw(1,1,inzee),zr(1,j,1,j3,idat), &
878 19830678 : & trig2,after2(i),now2(i),before2(i),1)
879 :
880 :
881 : end do
882 :
883 : ! Treat real wavefunctions.
884 9571089 : if (icplexwf==1) then
885 3882 : n1half=n1/2
886 : ! If odd
887 3882 : if (n1half*2/=n1) then
888 116800 : do i2=1,n2
889 115200 : zr(1,n1,i2,j3,idat)=zr(1,n1eff,i2,j3,idat)
890 116800 : zr(2,n1,i2,j3,idat)=zero
891 : end do
892 : end if
893 325850 : do i2=1,n2
894 14675946 : do i1=n1half,1,-1
895 14350096 : zr(1,2*i1-1,i2,j3,idat)=zr(1,i1,i2,j3,idat)
896 14350096 : zr(1,2*i1 ,i2,j3,idat)=zr(2,i1,i2,j3,idat)
897 14350096 : zr(2,2*i1-1,i2,j3,idat)=zero
898 14672064 : zr(2,2*i1 ,i2,j3,idat)=zero
899 : end do
900 : end do
901 : end if
902 :
903 : end if
904 :
905 : end do ! j3
906 : end do ! idat
907 :
908 798252 : ABI_FREE(trig1)
909 798252 : ABI_FREE(after1)
910 798252 : ABI_FREE(now1)
911 798252 : ABI_FREE(before1)
912 798252 : ABI_FREE(trig2)
913 798252 : ABI_FREE(after2)
914 798252 : ABI_FREE(now2)
915 798252 : ABI_FREE(before2)
916 798252 : ABI_FREE(trig3)
917 798252 : ABI_FREE(after3)
918 798252 : ABI_FREE(now3)
919 798252 : ABI_FREE(before3)
920 798252 : ABI_FREE(zmpi2)
921 798252 : ABI_FREE(zw)
922 798252 : ABI_FREE(zt)
923 798252 : if (nproc_fft>1) then
924 696430 : ABI_FREE(zmpi1)
925 : end if
926 :
927 : !call timab(541,2,tsec)
928 :
929 798252 : end subroutine sg2002_mpiback_wf
930 : !!***
931 :
932 : !----------------------------------------------------------------------
933 :
934 : !!****f* m_sg2002/sg2002_mpiforw_wf
935 : !! NAME
936 : !! sg2002_mpiforw_wf
937 : !!
938 : !! FUNCTION
939 : !! Does multiple 3-dim backward FFTs from real into Fourier space
940 : !! Adopt standard convention that isign=-1 for forward transform
941 : !! CALCULATES THE DISCRETE FOURIERTRANSFORM
942 : !!
943 : !! ZF(I1,I3,I2)=S_(j1,j2,j3) EXP(isign*i*2*pi*(j1*i1/n1+j2*i2/n2+j3*i3/n3)) ZR(j1,j2,j3)
944 : !!
945 : !! in parallel using MPI/OpenMP.
946 : !!
947 : !! INPUT:
948 : !! ZR: input array
949 : !! ZR(1,i1,i2,i3,idat)=real(R(i1,i2,i3,idat))
950 : !! ZR(2,i1,i2,i3,idat)=imag(R(i1,i2,i3,idat))
951 : !! i1=1,n1 , i2=1,n2 , i3=1,n3 , idat=1,ndat
952 : !! NOTE that ZR is changed by the routine
953 : !!
954 : !! n1,n2,n3: logical dimension of the transform. As transform lengths
955 : !! most products of the prime factors 2,3,5 are allowed.
956 : !! The detailed table with allowed transform lengths can
957 : !! be found in subroutine CTRIG
958 : !! nd1,nd2,nd3: Dimension of ZR
959 : !! nd3proc=((nd3-1)/nproc_fft)+1 maximal number of big box 3rd dim slices for one proc
960 : !!
961 : !! OUTPUT:
962 : !! ZF: output array (note the switch of i2 and i3)
963 : !! real(F(i1,i3,i2,idat))=ZF(1,i1,i3,i2,idat)
964 : !! imag(F(i1,i3,i2,idat))=ZF(2,i1,i3,i2,idat)
965 : !! max1 is positive or zero ; m1 >=max1+1
966 : !! i1= 1... max1+1 corresponds to positive and zero wavevectors 0 ... max1
967 : !! then, if m1 > max1+1, one has min1=max1-m1+1 and
968 : !! i1= max1+2 ... m1 corresponds to negative wavevectors min1 ... -1
969 : !! i2 and i3 have a similar definition of range
970 : !! idat=1,ndat
971 : !! md1,md2,md3: Dimension of ZF
972 : !! md2proc=((md2-1)/nproc_fft)+1 maximal number of small box 2nd dim slices for one proc
973 : !! nproc_fft: number of processors used as returned by MPI_COMM_SIZE
974 : !! me_fft: [0:nproc-1] rank of the processor in the FFT communicator.
975 : !! comm_fft=MPI communicator for parallel FFT.
976 : !!
977 : !! NOTES
978 : !! The maximum number of processors that can reasonably be used is max(n2/2,n3/2)
979 : !!
980 : !! It is very important to find the optimal
981 : !! value of NCACHE. NCACHE determines the size of the work array ZW, that
982 : !! has to fit into cache. It has therefore to be chosen to equal roughly
983 : !! half the size of the physical cache in units of real*8 numbers.
984 : !! The optimal value of ncache can easily be determined by numerical
985 : !! experimentation. A too large value of ncache leads to a dramatic
986 : !! and sudden decrease of performance, a too small value to a to a
987 : !! slow and less dramatic decrease of performance. If NCACHE is set
988 : !! to a value so small, that not even a single one dimensional transform
989 : !! can be done in the workarray zw, the program stops with an error message.
990 : !!
991 : !! SOURCE
992 :
993 735144 : subroutine sg2002_mpiforw_wf(icplexwf,ndat,n1,n2,n3,nd1,nd2,nd3proc,&
994 735144 : & max1,max2,max3,m1,m2,m3,md1,md2proc,md3,zr,zf,comm_fft)
995 :
996 : !Arguments ------------------------------------
997 : !scalars
998 : integer,intent(in) :: icplexwf,ndat,n1,n2,n3,nd1,nd2,nd3proc
999 : integer,intent(in) :: max1,max2,max3,m1,m2,m3,md1,md2proc,md3,comm_fft
1000 : !arrays
1001 : real(dp),intent(inout) :: zr(2,nd1,nd2,nd3proc,ndat)
1002 : real(dp),intent(out) :: zf(2,md1,md3,md2proc,ndat)
1003 :
1004 : !Local variables-------------------------------
1005 : !scalars
1006 : integer :: i,j,i1,i2,i3,ic1,ic2,ic3,idat,ierr,inzee,nproc_fft,me_fft
1007 : integer :: ioption,j2,j3,j2st,jp2st,lot,lzt,m1zt,ma,mb,n1dfft,nnd3
1008 : integer :: m2eff,ncache,n1eff,n1half,i1inv,i2inv,i3inv
1009 : character(len=500) :: msg
1010 : !arrays
1011 735144 : real(dp), allocatable :: zt(:,:,:) ! work arrays for transpositions
1012 735144 : real(dp), allocatable :: zmpi1(:,:,:,:,:),zmpi2(:,:,:,:,:) ! work arrays for MPI
1013 735144 : real(dp), allocatable :: zw(:,:,:) ! cache work array
1014 : ! FFT work arrays
1015 735144 : real(dp), allocatable :: trig1(:,:),trig2(:,:),trig3(:,:)
1016 735144 : integer, allocatable :: after1(:),now1(:),before1(:),after2(:),now2(:),before2(:),after3(:),now3(:),before3(:)
1017 : real(dp) :: tsec(2)
1018 :
1019 : ! *************************************************************************
1020 :
1021 : ! call timab(542,1,tsec)
1022 :
1023 : ! FIXME must provide a default value but which one?
1024 : !ioption = 0
1025 735144 : ioption = 1
1026 : !if (paral_kgb==1) ioption=1
1027 :
1028 735144 : nproc_fft = xmpi_comm_size(comm_fft); me_fft = xmpi_comm_rank(comm_fft)
1029 :
1030 : ! find cache size that gives optimal performance on machine
1031 735144 : ncache=4*max(n1,n2,n3,1024)
1032 735144 : if (ncache/(4*max(n1,n2,n3))<1) then
1033 : write(msg,'(5a)') &
1034 0 : & 'ncache has to be enlarged to be able to hold at',ch10, &
1035 0 : & 'least one 1-d FFT of each size even though this will',ch10,&
1036 0 : & 'reduce the performance for shorter transform lengths'
1037 0 : ABI_ERROR(msg)
1038 : end if
1039 :
1040 : ! Effective m1 and m2 (complex-to-complex or real-to-complex)
1041 735144 : n1eff=n1; m2eff=m2; m1zt=n1
1042 735144 : if (icplexwf==1) then
1043 35 : n1eff=(n1+1)/2; m2eff=m2/2+1; m1zt=2*(n1/2+1)
1044 : end if
1045 :
1046 735144 : lzt=m2eff
1047 735144 : if (mod(m2eff,2)==0) lzt=lzt+1
1048 652351 : if (mod(m2eff,4)==0) lzt=lzt+1
1049 :
1050 : ! maximal number of big box 3rd dim slices for all procs
1051 735144 : nnd3=nd3proc*nproc_fft
1052 :
1053 2205432 : ABI_MALLOC(trig1,(2,n1))
1054 735144 : ABI_MALLOC(after1,(mdata))
1055 735144 : ABI_MALLOC(now1,(mdata))
1056 735144 : ABI_MALLOC(before1,(mdata))
1057 2205432 : ABI_MALLOC(trig2,(2,n2))
1058 735144 : ABI_MALLOC(after2,(mdata))
1059 735144 : ABI_MALLOC(now2,(mdata))
1060 735144 : ABI_MALLOC(before2,(mdata))
1061 2205432 : ABI_MALLOC(trig3,(2,n3))
1062 735144 : ABI_MALLOC(after3,(mdata))
1063 735144 : ABI_MALLOC(now3,(mdata))
1064 735144 : ABI_MALLOC(before3,(mdata))
1065 2940576 : ABI_MALLOC(zw,(2,ncache/4,2))
1066 2940576 : ABI_MALLOC(zt,(2,lzt,m1zt))
1067 4410864 : ABI_MALLOC(zmpi2,(2,md1,md2proc,nnd3,ndat))
1068 735144 : if (nproc_fft>1) then
1069 3242980 : ABI_MALLOC(zmpi1,(2,md1,md2proc,nnd3,ndat))
1070 : end if
1071 :
1072 735144 : call ctrig(n2,trig2,after2,before2,now2,-1,ic2)
1073 735144 : call ctrig(n1,trig1,after1,before1,now1,-1,ic1)
1074 735144 : call ctrig(n3,trig3,after3,before3,now3,-1,ic3)
1075 :
1076 : !DEBUG
1077 : ! write(std_out,'(2a,3i4)' )itoa(me_fft),'sg2002_mpiforw_wf, enter', i1,i2,i3,zr,n1,n2,n3',n1,n2,n3
1078 : ! write(std_out,'(2a,3i4)' )itoa(me_fft),'nd1,nd2,nd3proc',nd1,nd2,nd3proc
1079 : ! write(std_out,'(2a,3i4)' )itoa(me_fft),'m1,m2,m3',m1,m2,m3
1080 : ! write(std_out,'(2a,3i4)' )itoa(me_fft),'max1,max2,max3',max1,max2,max3
1081 : ! write(std_out,'(2a,3i4)' )itoa(me_fft),'md1,md2proc,md3',md1,md2proc,md3
1082 : ! write(std_out,'(2a,3i4)' )itoa(me_fft),'n1eff,m2eff,m1zt',n1eff,m2eff,m1zt
1083 : !ENDDEBUG
1084 :
1085 1520819 : do idat=1,ndat
1086 : ! Loop over the z-planes treated by this node
1087 9598729 : do j3=1,nd3proc
1088 :
1089 9598729 : if (me_fft*nd3proc+j3 <= n3) then
1090 8813054 : Jp2st=1
1091 8813054 : J2st=1
1092 :
1093 : ! Treat real wavefunctions.
1094 8813054 : if (icplexwf==1) then
1095 2826 : n1half=n1/2
1096 242250 : do i2=1,n2
1097 10977498 : do i1=1,n1half
1098 10735248 : zr(1,i1,i2,j3,idat)=zr(1,2*i1-1,i2,j3,idat)
1099 10974672 : zr(2,i1,i2,j3,idat)=zr(1,2*i1 ,i2,j3,idat)
1100 : end do
1101 : end do
1102 : ! If odd
1103 2826 : if(n1half*2/=n1)then
1104 87600 : do i2=1,n2
1105 86400 : zr(1,n1eff,i2,j3,idat)=zr(1,n1,i2,j3,idat)
1106 87600 : zr(2,n1eff,i2,j3,idat)=zero
1107 : end do
1108 : end if
1109 : end if
1110 :
1111 : ! transform along y axis
1112 : ! input: R1,R2,R3,(Rp3)
1113 : ! input: R1,G2,R3,(Rp3)
1114 8813054 : lot=ncache/(4*n2)
1115 :
1116 18197068 : do j=1,n1eff,lot
1117 9384014 : ma=j
1118 9384014 : mb=min(j+(lot-1),n1eff)
1119 9384014 : n1dfft=mb-ma+1
1120 9384014 : i=1
1121 : call fftstp(nd1,n1dfft,nd2,lot,n2,zr(1,j,1,j3,idat),zw(1,1,1), &
1122 9384014 : & trig2,after2(i),now2(i),before2(i),-1)
1123 :
1124 9384014 : inzee=1
1125 19577528 : do i=2,ic2
1126 : call fftstp(lot,n1dfft,n2,lot,n2,zw(1,1,inzee),zw(1,1,3-inzee), &
1127 10193514 : & trig2,after2(i),now2(i),before2(i),-1)
1128 19577528 : inzee=3-inzee
1129 : end do
1130 :
1131 : ! input: R1,G2,R3,(Rp3)
1132 : ! output: G2,R1,R3,(Rp3)
1133 18197068 : if(icplexwf==2)then
1134 9372788 : call unswitch_cent(n1dfft,max2,m2,n2,lot,n1,lzt,zw(1,1,inzee),zt(1,1,j))
1135 : else
1136 11226 : call unswitchreal_cent(n1dfft,max2,n2,lot,n1,lzt,zw(1,1,inzee),zt(1,1,2*j-1))
1137 : end if
1138 : end do
1139 :
1140 : ! transform along x axis
1141 : ! input: G2,R1,R3,(Rp3)
1142 8813054 : lot=ncache/(4*n1)
1143 :
1144 17907588 : do j=1,m2eff,lot
1145 9094534 : ma=j
1146 9094534 : mb=min(j+(lot-1),m2eff)
1147 9094534 : n1dfft=mb-ma+1
1148 9094534 : i=1
1149 : call fftstp(lzt,n1dfft,m1zt,lot,n1,zt(1,j,1),zw(1,1,1), &
1150 9094534 : & trig1,after1(i),now1(i),before1(i),-1)
1151 :
1152 9094534 : inzee=1
1153 18944468 : do i=2,ic1
1154 : call fftstp(lot,n1dfft,n1,lot,n1,zw(1,1,inzee),zw(1,1,3-inzee), &
1155 9849934 : & trig1,after1(i),now1(i),before1(i),-1)
1156 18944468 : inzee=3-inzee
1157 : end do
1158 : ! output: G2,G1,R3,(Rp3)
1159 :
1160 : ! input: G2,G1,R3,Gp2,(Rp3)
1161 : ! output: G1,G2,R3,Gp2,(Rp3)
1162 17907588 : if (nproc_fft==1) then
1163 : call unmpiswitch_cent(j3,n1dfft,Jp2st,J2st,lot,max1,md1,m1,n1,&
1164 2303766 : & md2proc,nd3proc,nproc_fft,ioption,zw(1,1,inzee),zmpi2(:,:,:,:,idat))
1165 : else
1166 : call unmpiswitch_cent(j3,n1dfft,Jp2st,J2st,lot,max1,md1,m1,n1,&
1167 6790768 : & md2proc,nd3proc,nproc_fft,ioption,zw(1,1,inzee),zmpi1(:,:,:,:,idat))
1168 : end if
1169 : end do
1170 :
1171 : end if
1172 : end do ! j3
1173 :
1174 : ! Interprocessor data transposition
1175 : ! input: G1,G2,R3,Gp2,(Rp3)
1176 : ! output: G1,G2,R3,Rp3,(Gp2)
1177 785675 : if (nproc_fft>1) then
1178 667240 : call timab(544,1,tsec)
1179 : call xmpi_alltoall(zmpi1(:,:,:,:,idat),2*md1*md2proc*nd3proc, &
1180 667240 : & zmpi2(:,:,:,:,idat),2*md1*md2proc*nd3proc,comm_fft,ierr)
1181 :
1182 667240 : call timab(544,2,tsec)
1183 : end if
1184 :
1185 : ! transform along z axis
1186 : ! input: G1,G2,R3,(Gp2)
1187 785675 : lot=ncache/(4*n3)
1188 :
1189 5032424 : do j2=1,md2proc
1190 5032424 : if (me_fft*md2proc+j2 <= m2eff) then
1191 : ! write(std_out,*)' forwf_wf : before unscramble, j2,md2proc,me_fft,m2=',j2,md2proc,me_fft,m2
1192 8627848 : do i1=1,m1,lot
1193 4381724 : ma=i1
1194 4381724 : mb=min(i1+(lot-1),m1)
1195 4381724 : n1dfft=mb-ma+1
1196 :
1197 : ! input: G1,G2,R3,(Gp2)
1198 : ! output: G1,R3,G2,(Gp2)
1199 4381724 : call unscramble(i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3,zmpi2(:,:,:,:,idat),zw(1,1,1))
1200 :
1201 4381724 : inzee=1
1202 13394469 : do i=1,ic3
1203 : call fftstp(lot,n1dfft,n3,lot,n3,zw(1,1,inzee),zw(1,1,3-inzee), &
1204 9012745 : & trig3,after3(i),now3(i),before3(i),-1)
1205 13394469 : inzee=3-inzee
1206 : end do
1207 :
1208 8627848 : call unfill_cent(md1,md3,lot,n1dfft,max3,m3,n3,zw(1,1,inzee),zf(1,i1,1,j2,idat))
1209 : ! output: G1,G3,G2,(Gp2)
1210 : end do
1211 : end if
1212 : end do
1213 :
1214 1520819 : if (icplexwf==1) then
1215 : ! Complete missing values with complex conjugate
1216 : ! Inverse of ix is located at nx+2-ix , except for ix=1, for which it is 1.
1217 1414 : do i3=1,m3
1218 1379 : i3inv=m3+2-i3
1219 1379 : if(i3==1)i3inv=1
1220 :
1221 1414 : if (m2eff>1) then
1222 28674 : do i2=2,m2eff
1223 27295 : i2inv=m2+2-i2
1224 27295 : zf(1,1,i3inv,i2inv,idat)= zf(1,1,i3,i2,idat)
1225 27295 : zf(2,1,i3inv,i2inv,idat)=-zf(2,1,i3,i2,idat)
1226 1215604 : do i1=2,m1
1227 1186930 : i1inv=m1+2-i1
1228 1186930 : zf(1,i1inv,i3inv,i2inv,idat)= zf(1,i1,i3,i2,idat)
1229 1214225 : zf(2,i1inv,i3inv,i2inv,idat)=-zf(2,i1,i3,i2,idat)
1230 : end do
1231 : end do
1232 : end if
1233 : end do
1234 : end if
1235 :
1236 : end do ! idat
1237 :
1238 735144 : ABI_FREE(trig1)
1239 735144 : ABI_FREE(after1)
1240 735144 : ABI_FREE(now1)
1241 735144 : ABI_FREE(before1)
1242 735144 : ABI_FREE(trig2)
1243 735144 : ABI_FREE(after2)
1244 735144 : ABI_FREE(now2)
1245 735144 : ABI_FREE(before2)
1246 735144 : ABI_FREE(trig3)
1247 735144 : ABI_FREE(after3)
1248 735144 : ABI_FREE(now3)
1249 735144 : ABI_FREE(before3)
1250 735144 : ABI_FREE(zmpi2)
1251 735144 : ABI_FREE(zw)
1252 735144 : ABI_FREE(zt)
1253 735144 : if (nproc_fft>1) then
1254 648596 : ABI_FREE(zmpi1)
1255 : end if
1256 :
1257 : !call timab(542,2,tsec)
1258 :
1259 735144 : end subroutine sg2002_mpiforw_wf
1260 : !!***
1261 :
1262 : !----------------------------------------------------------------------
1263 :
1264 : !!****f* m_sg2002/sg2002_mpifourdp
1265 : !! NAME
1266 : !! sg2002_mpifourdp
1267 : !!
1268 : !! FUNCTION
1269 : !! Conduct Fourier transform of REAL or COMPLEX function f(r)=fofr defined on
1270 : !! fft grid in real space, to create complex f(G)=fofg defined on full fft grid
1271 : !! in reciprocal space, in full storage mode, or the reverse operation.
1272 : !! For the reverse operation, the final data is divided by nfftot.
1273 : !! REAL case when cplex=1, COMPLEX case when cplex=2
1274 : !! Usually used for density and potentials.
1275 : !!
1276 : !! INPUTS
1277 : !! cplex=1 if fofr is real, 2 if fofr is complex
1278 : !! nfft=(effective) number of FFT grid points (for this processor)
1279 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1280 : !! ndat=Numbre of FFT transforms
1281 : !! isign=sign of Fourier transform exponent: current convention uses
1282 : !! +1 for transforming from G to r
1283 : !! -1 for transforming from r to G.
1284 : !! fftn2_distrib(2),ffti2_local(2)
1285 : !! fftn3_distrib(3),ffti3_local(3)
1286 : !! comm_fft=MPI communicator
1287 : !!
1288 : !! SIDE EFFECTS
1289 : !! Input/Output
1290 : !! fofg(2,nfft)=f(G), complex.
1291 : !! fofr(cplex*nfft)=input function f(r) (real or complex)
1292 : !!
1293 : !! TODO
1294 : !! Write simplified API for sequential version.
1295 : !!
1296 : !! SOURCE
1297 :
1298 135120 : subroutine sg2002_mpifourdp(cplex,nfft,ngfft,ndat,isign,&
1299 135120 : & fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local,fofg,fofr,comm_fft)
1300 :
1301 : !Arguments ------------------------------------
1302 : !scalars
1303 : integer,intent(in) :: cplex,isign,nfft,ndat,comm_fft
1304 : !arrays
1305 : integer,intent(in) :: ngfft(18)
1306 : integer,intent(in) :: fftn2_distrib(ngfft(2)),ffti2_local(ngfft(2))
1307 : integer,intent(in) :: fftn3_distrib(ngfft(3)),ffti3_local(ngfft(3))
1308 : real(dp),intent(inout) :: fofg(2,nfft*ndat),fofr(cplex*nfft*ndat)
1309 :
1310 : !Local variables-------------------------------
1311 : !scalars
1312 : integer :: n1,n2,n3,n4,n5,n6,nd2proc,nd3proc,nproc_fft,me_fft
1313 : !arrays
1314 135120 : real(dp),allocatable :: workf(:,:,:,:,:),workr(:,:,:,:,:)
1315 :
1316 : ! *************************************************************************
1317 :
1318 : ! Note the only c2c is supported in parallel.
1319 135120 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
1320 135120 : n4=ngfft(4); n5=ngfft(5); n6=ngfft(6)
1321 135120 : me_fft=ngfft(11); nproc_fft=ngfft(10)
1322 :
1323 135120 : nd2proc=((n2-1)/nproc_fft) +1
1324 135120 : nd3proc=((n6-1)/nproc_fft) +1
1325 810720 : ABI_MALLOC(workr,(2,n4,n5,nd3proc,ndat))
1326 810720 : ABI_MALLOC(workf,(2,n4,n6,nd2proc,ndat))
1327 :
1328 : ! Complex to Complex
1329 202803 : select case (isign)
1330 : case (1)
1331 : ! G --> R
1332 67683 : call mpifft_fg2dbox(nfft,ndat,fofg,n1,n2,n3,n4,nd2proc,n6,fftn2_distrib,ffti2_local,me_fft,workf)
1333 :
1334 67683 : call sg2002_back(2,ndat,n1,n2,n3,n4,n5,n6,n4,nd2proc,nd3proc,2,workf,workr,comm_fft)
1335 :
1336 67683 : call mpifft_dbox2fr(n1,n2,n3,n4,n5,nd3proc,ndat,fftn3_distrib,ffti3_local,me_fft,workr,cplex,nfft,fofr)
1337 :
1338 : case (-1)
1339 : ! R --> G
1340 67437 : call mpifft_fr2dbox(cplex,nfft,ndat,fofr,n1,n2,n3,n4,n5,nd3proc,fftn3_distrib,ffti3_local,me_fft,workr)
1341 :
1342 67437 : call sg2002_forw(2,ndat,n1,n2,n3,n4,n5,n6,n4,nd2proc,nd3proc,2,workr,workf,comm_fft)
1343 :
1344 : ! Transfer FFT output to the original fft box.
1345 67437 : call mpifft_dbox2fg(n1,n2,n3,n4,nd2proc,n6,ndat,fftn2_distrib,ffti2_local,me_fft,workf,nfft,fofg)
1346 :
1347 : case default
1348 135120 : ABI_BUG("Wrong isign")
1349 : end select
1350 :
1351 135120 : ABI_FREE(workr)
1352 135120 : ABI_FREE(workf)
1353 :
1354 135120 : end subroutine sg2002_mpifourdp
1355 : !!***
1356 :
1357 : !----------------------------------------------------------------------
1358 :
1359 : !!****f* m_sg2002/sg2002_applypot
1360 : !! NAME
1361 : !! sg2002_applypot
1362 : !!
1363 : !! FUNCTION
1364 : !! Applies the local real space potential to multiple wavefunctions in Fourier space
1365 : !!
1366 : !! INPUTS
1367 : !! ZF: Wavefunction (input/output) (note the switch of i2 and i3)
1368 : !! real(F(i1,i3,i2,idat))=ZF(1,i1,i3,i2,idat)
1369 : !! imag(F(i1,i3,i2,idat))=ZF(2,i1,i3,i2,idat)
1370 : !! max1 is positive or zero ; m1 >=max1+1
1371 : !! i1= 1... max1+1 corresponds to positive and zero wavevectors 0 ... max1
1372 : !! then, if m1 > max1+1, one has min1=max1-m1+1 and
1373 : !! i1= max1+2 ... m1 corresponds to negative wavevectors min1 ... -1
1374 : !! i2 and i3 have a similar definition of range
1375 : !! idat=1,ndat
1376 : !! md1,md2,md3: Dimension of ZF (input as well as output), distributed on different procs
1377 : !! md2proc=((md2-1)/nproc_fft)+1 maximal number of small box 2nd dim slices for one proc
1378 : !!
1379 : !! POT: Potential
1380 : !! POT(cplex*i1,i2,i3)
1381 : !! cplex=1 or 2 , i1=1,n1 , i2=1,n2 , i3=1,n3
1382 : !! nd1,nd2,nd3: dimension of pot
1383 : !! comm_fft: MPI communicator
1384 : !! nproc_fft: number of processors used as returned by MPI_COMM_SIZE
1385 : !! me_fft: [0:nproc_fft-1] number of processor as returned by MPI_COMM_RANK
1386 : !! n1,n2,n3: logical dimension of the transform. As transform lengths
1387 : !! most products of the prime factors 2,3,5 are allowed.
1388 : !! The detailed table with allowed transform lengths can
1389 : !! be found in subroutine CTRIG
1390 : !!
1391 : !! NOTES:
1392 : !! PERFORMANCE CONSIDERATIONS:
1393 : !! The maximum number of processors that can reasonably be used is max(n2/2,n3/2)
1394 : !!
1395 : !! It is very important to find the optimal
1396 : !! value of NCACHE. NCACHE determines the size of the work array ZW, that
1397 : !! has to fit into cache. It has therefore to be chosen to equal roughly
1398 : !! half the size of the physical cache in units of real*8 numbers.
1399 : !! The optimal value of ncache can easily be determined by numerical
1400 : !! experimentation. A too large value of ncache leads to a dramatic
1401 : !! and sudden decrease of performance, a too small value to a to a
1402 : !! slow and less dramatic decrease of performance. If NCACHE is set
1403 : !! to a value so small, that not even a single one dimensional transform
1404 : !! can be done in the workarray zw, the program stops with an error message.
1405 : !!
1406 : !! SOURCE
1407 :
1408 :
1409 223 : subroutine sg2002_applypot(icplexwf,cplex,ndat,n1,n2,n3,nd1,nd2,nd3,nd3proc,&
1410 : & max1i,max2i,max3i,m1i,m2i,m3i,md1,md2proc,md3,&
1411 223 : & max1o,max2o,max3o,m1o,m2o,m3o,comm_fft,nproc_fft,me_fft,pot,zf)
1412 :
1413 : !Arguments ------------------------------------
1414 : integer,intent(in) :: icplexwf,cplex,ndat,n1,n2,n3,nd1,nd2,nd3,nd3proc
1415 : integer,intent(in) :: max1i,max2i,max3i,m1i,m2i,m3i,md1,md2proc,md3
1416 : integer,intent(in) :: max1o,max2o,max3o,m1o,m2o,m3o,comm_fft,nproc_fft,me_fft
1417 : real(dp),intent(in) :: pot(cplex*nd1,nd2,nd3)
1418 : real(dp),intent(inout) :: zf(2,md1,md3,md2proc,ndat)
1419 :
1420 : !Local variables-------------------------------
1421 : !scalars
1422 : integer,parameter :: unused0=0
1423 : integer :: i,j,i1,i2,i3,ic1,ic2,ic3,idat,ierr,inzee,j3glob
1424 : integer :: ioption,j2,j3,lot,lzt,m1zt,ma,mb,n1dfft,nnd3
1425 : integer :: m2eff,ncache,n1eff,i1inv,i2inv,i3inv,jeff,includelast,j2stb
1426 : integer :: jx,j2stf,Jp2stb,Jp2stf,m2ieff,m2oeff
1427 : !arrays
1428 : real(dp) :: tsec(2)
1429 223 : real(dp), allocatable :: zt(:,:,:) ! work arrays for transpositions
1430 223 : real(dp), allocatable :: zmpi1(:,:,:,:),zmpi2(:,:,:,:) ! work arrays for MPI
1431 223 : real(dp), allocatable :: zw(:,:,:) ! cache work array
1432 : ! FFT work arrays
1433 223 : real(dp), allocatable, dimension(:,:) :: btrig1,btrig2,btrig3
1434 223 : real(dp), allocatable, dimension(:,:) :: ftrig1,ftrig2,ftrig3
1435 223 : integer, allocatable, dimension(:) :: after1,now1,before1,after2,now2,before2,after3,now3,before3
1436 :
1437 : ! *************************************************************************
1438 :
1439 : !ioption=0 ! This was in the old version.
1440 223 : ioption=1 ! This one is needed to be compatible with paral_kgb
1441 :
1442 223 : ncache=4*max(n1,n2,n3,1024)
1443 223 : if (ncache/(4*max(n1,n2,n3)) < 1) then
1444 : write(std_out,*) &
1445 0 : & 'ncache has to be enlarged to be able to hold at', &
1446 0 : & 'least one 1-d FFT of each size even though this will', &
1447 0 : & 'reduce the performance for shorter transform lengths'
1448 0 : ABI_ERROR("Aborting now")
1449 : end if
1450 :
1451 : ! Effective m1 and m2 (complex-to-complex or real-to-complex)
1452 223 : n1eff=n1; m2ieff=m2i; m2oeff=m2o; m1zt=n1
1453 223 : if (icplexwf==1) then
1454 15 : n1eff=(n1+1)/2; m2ieff=m2i/2+1; m2oeff=m2o/2+1; m1zt=2*(n1/2+1)
1455 : end if
1456 :
1457 223 : m2eff=max(m2ieff,m2oeff)
1458 223 : lzt=m2eff
1459 223 : if (mod(m2eff,2) == 0) lzt=lzt+1
1460 39 : if (mod(m2eff,4) == 0) lzt=lzt+1
1461 :
1462 : ! maximal number of big box 3rd dim slices for all procs
1463 223 : nnd3=nd3proc*nproc_fft
1464 :
1465 669 : ABI_MALLOC(btrig1,(2,n1))
1466 446 : ABI_MALLOC(ftrig1,(2,n1))
1467 223 : ABI_MALLOC(after1,(mdata))
1468 223 : ABI_MALLOC(now1,(mdata))
1469 223 : ABI_MALLOC(before1,(mdata))
1470 669 : ABI_MALLOC(btrig2,(2,n2))
1471 446 : ABI_MALLOC(ftrig2,(2,n2))
1472 223 : ABI_MALLOC(after2,(mdata))
1473 223 : ABI_MALLOC(now2,(mdata))
1474 223 : ABI_MALLOC(before2,(mdata))
1475 669 : ABI_MALLOC(btrig3,(2,n3))
1476 446 : ABI_MALLOC(ftrig3,(2,n3))
1477 223 : ABI_MALLOC(after3,(mdata))
1478 223 : ABI_MALLOC(now3,(mdata))
1479 223 : ABI_MALLOC(before3,(mdata))
1480 :
1481 892 : ABI_MALLOC(zw,(2,ncache/4,2))
1482 892 : ABI_MALLOC(zt,(2,lzt,m1zt))
1483 1115 : ABI_MALLOC(zmpi2,(2,md1,md2proc,nnd3))
1484 223 : if (nproc_fft > 1) then
1485 0 : ABI_MALLOC(zmpi1,(2,md1,md2proc,nnd3))
1486 : end if
1487 :
1488 223 : call ctrig(n3,btrig3,after3,before3,now3,1,ic3)
1489 223 : call ctrig(n1,btrig1,after1,before1,now1,1,ic1)
1490 223 : call ctrig(n2,btrig2,after2,before2,now2,1,ic2)
1491 :
1492 13345 : do j=1,n1
1493 13122 : ftrig1(1,j)= btrig1(1,j)
1494 13345 : ftrig1(2,j)=-btrig1(2,j)
1495 : end do
1496 14701 : do j=1,n2
1497 14478 : ftrig2(1,j)= btrig2(1,j)
1498 14701 : ftrig2(2,j)=-btrig2(2,j)
1499 : end do
1500 10569 : do j=1,n3
1501 10346 : ftrig3(1,j)= btrig3(1,j)
1502 10569 : ftrig3(2,j)=-btrig3(2,j)
1503 : end do
1504 :
1505 446 : do idat=1,ndat
1506 : !
1507 : ! transform along z axis
1508 : ! input: G1,G3,G2,(Gp2)
1509 223 : lot=ncache/(4*n3)
1510 6924 : do j2=1,md2proc
1511 6924 : if (me_fft*md2proc+j2 <= m2ieff) then
1512 20118 : do i1=1,m1i,lot
1513 13642 : ma=i1
1514 13642 : mb=min(i1+(lot-1),m1i)
1515 13642 : n1dfft=mb-ma+1
1516 :
1517 : ! zero-pad n1dfft G_z lines
1518 : ! input: G1,G3,G2,(Gp2)
1519 13642 : call fill_cent(md1,md3,lot,n1dfft,max3i,m3i,n3,zf(1,i1,1,j2,idat),zw(1,1,1))
1520 :
1521 13642 : inzee=1
1522 54287 : do i=1,ic3
1523 : call fftstp(lot,n1dfft,n3,lot,n3,zw(1,1,inzee),zw(1,1,3-inzee), &
1524 40645 : & btrig3,after3(i),now3(i),before3(i),1)
1525 54287 : inzee=3-inzee
1526 : end do
1527 :
1528 : ! Local rotation.
1529 : ! input: G1,R3,G2,(Gp2)
1530 : ! output: G1,G2,R3,(Gp2)
1531 20118 : call scramble(i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3,zw(1,1,inzee),zmpi2)
1532 : end do
1533 : end if
1534 : end do
1535 :
1536 : ! Interprocessor data transposition
1537 : ! input: G1,G2,R3,Rp3,(Gp2)
1538 : ! output: G1,G2,R3,Gp2,(Rp3)
1539 223 : if (nproc_fft > 1) then
1540 0 : call timab(543,1,tsec)
1541 : call xmpi_alltoall(zmpi2,2*md1*md2proc*nd3proc,&
1542 0 : & zmpi1,2*md1*md2proc*nd3proc,comm_fft,ierr)
1543 0 : call timab(543,2,tsec)
1544 : end if
1545 :
1546 10569 : do j3=1,nd3proc
1547 10346 : j3glob = j3 + me_fft*nd3proc
1548 :
1549 10569 : if (me_fft*nd3proc+j3 <= n3) then
1550 10346 : Jp2stb=1; J2stb=1
1551 10346 : Jp2stf=1; J2stf=1
1552 :
1553 : ! transform along x axis
1554 10346 : lot=ncache/(4*n1)
1555 :
1556 39516 : do j=1,m2ieff,lot
1557 29170 : ma=j
1558 29170 : mb=min(j+(lot-1),m2ieff)
1559 29170 : n1dfft=mb-ma+1
1560 :
1561 : ! Zero-pad input.
1562 : ! input: G1,G2,R3,G2,(Rp3)
1563 : ! output: G2,G1,R3,G2,(Rp3)
1564 29170 : if (nproc_fft == 1) then
1565 : call mpiswitch_cent(j3,n1dfft,Jp2stb,J2stb,lot,max1i,md1,m1i,n1,&
1566 29170 : & md2proc,nd3proc,nproc_fft,ioption,zmpi2,zw(1,1,1), unused0, unused0, unused0)
1567 : else
1568 : call mpiswitch_cent(j3,n1dfft,Jp2stb,J2stb,lot,max1i,md1,m1i,n1,&
1569 0 : & md2proc,nd3proc,nproc_fft,ioption,zmpi1,zw(1,1,1), unused0, unused0, unused0)
1570 : end if
1571 :
1572 : ! Transform along x
1573 : ! input: G2,G1,R3,(Rp3)
1574 : ! output: G2,R1,R3,(Rp3)
1575 29170 : inzee=1
1576 86808 : do i=1,ic1-1
1577 : call fftstp(lot,n1dfft,n1,lot,n1,zw(1,1,inzee),zw(1,1,3-inzee), &
1578 57638 : & btrig1,after1(i),now1(i),before1(i),1)
1579 86808 : inzee=3-inzee
1580 : end do
1581 :
1582 29170 : i=ic1
1583 : call fftstp(lot,n1dfft,n1,lzt,m1zt,zw(1,1,inzee),zt(1,j,1), &
1584 39516 : & btrig1,after1(i),now1(i),before1(i),1)
1585 : end do
1586 :
1587 : ! Transform along y axis (take into account c2c or c2r case).
1588 : ! Must loop over the full box.
1589 10346 : lot=ncache/(4*n2)
1590 :
1591 10346 : if (icplexwf==1) then
1592 1026 : if(mod(lot,2).ne.0)lot=lot-1 ! needed to introduce jeff
1593 : end if
1594 :
1595 71068 : do j=1,n1eff,lot
1596 60722 : ma=j
1597 60722 : mb=min(j+(lot-1),n1eff)
1598 60722 : n1dfft=mb-ma+1
1599 60722 : jeff=j
1600 60722 : includelast=1
1601 :
1602 60722 : if (icplexwf==1) then
1603 3826 : jeff=2*j-1
1604 : includelast=1
1605 3826 : if (mb==n1eff .and. n1eff*2/=n1) includelast=0
1606 : end if
1607 :
1608 : ! Zero-pad the input.
1609 : ! input: G2,R1,R3,(Rp3)
1610 : ! output: R1,G2,R3,(Rp3)
1611 60722 : if (icplexwf==2) then
1612 56896 : call switch_cent(n1dfft,max2i,m2i,n2,lot,n1,lzt,zt(1,1,jeff),zw(1,1,1))
1613 : else
1614 3826 : call switchreal_cent(includelast,n1dfft,max2i,n2,lot,m1zt,lzt,zt(1,1,jeff),zw(1,1,1))
1615 : end if
1616 :
1617 : ! input: R1,G2,R3,(Rp3)
1618 : ! output: R1,R2,R3,(Rp3)
1619 60722 : inzee=1
1620 242186 : do i=1,ic2
1621 : call fftstp(lot,n1dfft,n2,lot,n2,zw(1,1,inzee),zw(1,1,3-inzee), &
1622 181464 : & btrig2,after2(i),now2(i),before2(i),1)
1623 242186 : inzee=3-inzee
1624 : end do
1625 : ! output: R1,R2,R3,(Rp3)
1626 :
1627 : ! Multiply with potential in real space
1628 60722 : jx=cplex*(jeff-1)+1
1629 60722 : call multpot(icplexwf,cplex,includelast,nd1,nd2,n2,lot,n1dfft,pot(jx,1,j3glob),zw(1,1,inzee))
1630 :
1631 : ! TRANSFORM BACK IN FOURIER SPACE
1632 : ! transform along y axis
1633 : ! input: R1,R2,R3,(Rp3)
1634 242186 : do i=1,ic2
1635 : call fftstp(lot,n1dfft,n2,lot,n2,zw(1,1,inzee),zw(1,1,3-inzee), &
1636 181464 : & ftrig2,after2(i),now2(i),before2(i),-1)
1637 242186 : inzee=3-inzee
1638 : end do
1639 :
1640 : ! input: R1,G2,R3,(Rp3)
1641 : ! output: G2,R1,R3,(Rp3)
1642 71068 : if (icplexwf==2) then
1643 56896 : call unswitch_cent(n1dfft,max2o,m2o,n2,lot,n1,lzt,zw(1,1,inzee),zt(1,1,jeff))
1644 : else
1645 3826 : call unswitchreal_cent(n1dfft,max2o,n2,lot,n1,lzt,zw(1,1,inzee),zt(1,1,jeff))
1646 : end if
1647 :
1648 : end do ! j
1649 :
1650 : ! transform along x axis
1651 : ! input: R2,R1,R3,(Rp3)
1652 : ! output: R2,G1,R3,(Rp3)
1653 10346 : lot=ncache/(4*n1)
1654 :
1655 39516 : do j=1,m2oeff,lot
1656 29170 : ma=j
1657 29170 : mb=min(j+(lot-1),m2oeff)
1658 29170 : n1dfft=mb-ma+1
1659 29170 : i=1
1660 : call fftstp(lzt,n1dfft,m1zt,lot,n1,zt(1,j,1),zw(1,1,1), &
1661 29170 : & ftrig1,after1(i),now1(i),before1(i),-1)
1662 :
1663 29170 : inzee=1
1664 86808 : do i=2,ic1
1665 : call fftstp(lot,n1dfft,n1,lot,n1,zw(1,1,inzee),zw(1,1,3-inzee), &
1666 57638 : & ftrig1,after1(i),now1(i),before1(i),-1)
1667 86808 : inzee=3-inzee
1668 : end do
1669 :
1670 : ! input: G2,G1,R3,Gp2,(Rp3)
1671 : ! output: G1,G2,R3,Gp2,(Rp3)
1672 39516 : if (nproc_fft == 1) then
1673 : call unmpiswitch_cent(j3,n1dfft,Jp2stf,J2stf,lot,max1o,md1,m1o,n1,&
1674 29170 : & md2proc,nd3proc,nproc_fft,ioption,zw(1,1,inzee),zmpi2)
1675 : else
1676 : call unmpiswitch_cent(j3,n1dfft,Jp2stf,J2stf,lot,max1o,md1,m1o,n1,&
1677 0 : & md2proc,nd3proc,nproc_fft,ioption,zw(1,1,inzee),zmpi1)
1678 : end if
1679 : end do ! j
1680 : end if
1681 : end do
1682 :
1683 : ! Interprocessor data transposition
1684 : ! input: G1,G2,R3,Gp2,(Rp3)
1685 : ! output: G1,G2,R3,Rp3,(Gp2)
1686 223 : if (nproc_fft > 1) then
1687 0 : call timab(544,1,tsec)
1688 : call xmpi_alltoall(zmpi1,2*md1*md2proc*nd3proc, &
1689 0 : & zmpi2,2*md1*md2proc*nd3proc,comm_fft,ierr)
1690 0 : call timab(544,2,tsec)
1691 : end if
1692 :
1693 : ! transform along z axis
1694 : ! input: G1,G2,R3,(Gp2)
1695 223 : lot=ncache/(4*n3)
1696 6924 : do j2=1,md2proc
1697 6924 : if (me_fft*md2proc+j2 <= m2oeff) then
1698 20118 : do i1=1,m1o,lot
1699 13642 : ma=i1
1700 13642 : mb=min(i1+(lot-1),m1o)
1701 13642 : n1dfft=mb-ma+1
1702 :
1703 : ! input: G1,G2,R3,(Gp2)
1704 : ! output: G1,R3,G2,(Gp2)
1705 13642 : call unscramble(i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3,zmpi2,zw(1,1,1))
1706 :
1707 13642 : inzee=1
1708 54287 : do i=1,ic3
1709 : call fftstp(lot,n1dfft,n3,lot,n3,zw(1,1,inzee),zw(1,1,3-inzee), &
1710 40645 : & ftrig3,after3(i),now3(i),before3(i),-1)
1711 54287 : inzee=3-inzee
1712 : end do
1713 :
1714 20118 : call unfill_cent(md1,md3,lot,n1dfft,max3o,m3o,n3,zw(1,1,inzee),zf(1,i1,1,j2,idat))
1715 : ! output: G1,G3,G2,(Gp2)
1716 : end do
1717 : end if
1718 : end do
1719 :
1720 : ! Complete missing values with complex conjugate
1721 : ! Inverse of ix is located at nx+2-ix , except for ix=1, for which it is 1.
1722 446 : if (icplexwf==1) then
1723 514 : do i3=1,m3o
1724 499 : i3inv=m3o+2-i3
1725 499 : if (i3==1) i3inv=1
1726 514 : if (m2oeff>1)then
1727 9794 : do i2=2,m2oeff
1728 9295 : i2inv=m2o+2-i2
1729 9295 : zf(1,1,i3inv,i2inv,idat)= zf(1,1,i3,i2,idat)
1730 9295 : zf(2,1,i3inv,i2inv,idat)=-zf(2,1,i3,i2,idat)
1731 407604 : do i1=2,m1o
1732 397810 : i1inv=m1o+2-i1
1733 397810 : zf(1,i1inv,i3inv,i2inv,idat)= zf(1,i1,i3,i2,idat)
1734 407105 : zf(2,i1inv,i3inv,i2inv,idat)=-zf(2,i1,i3,i2,idat)
1735 : end do
1736 : end do
1737 : end if
1738 : end do
1739 : end if
1740 :
1741 : end do ! idat
1742 :
1743 223 : ABI_FREE(btrig1)
1744 223 : ABI_FREE(ftrig1)
1745 223 : ABI_FREE(after1)
1746 223 : ABI_FREE(now1)
1747 223 : ABI_FREE(before1)
1748 223 : ABI_FREE(btrig2)
1749 223 : ABI_FREE(ftrig2)
1750 223 : ABI_FREE(after2)
1751 223 : ABI_FREE(now2)
1752 223 : ABI_FREE(before2)
1753 223 : ABI_FREE(btrig3)
1754 223 : ABI_FREE(ftrig3)
1755 223 : ABI_FREE(after3)
1756 223 : ABI_FREE(now3)
1757 223 : ABI_FREE(before3)
1758 :
1759 223 : ABI_FREE(zmpi2)
1760 223 : ABI_FREE(zw)
1761 223 : ABI_FREE(zt)
1762 223 : if (nproc_fft > 1) then
1763 0 : ABI_FREE(zmpi1)
1764 : end if
1765 :
1766 223 : end subroutine sg2002_applypot
1767 : !!***
1768 :
1769 : !----------------------------------------------------------------------
1770 :
1771 : !!****f* m_sg2002/sg2002_applypot_many
1772 : !! NAME
1773 : !! sg2002_applypot_many
1774 : !!
1775 : !! FUNCTION
1776 : !! Applies the local real space potential to multiple wavefunctions in Fourier space
1777 : !!
1778 : !! INPUTS
1779 : !! ZF: Wavefunction (input/output) (note the switch of i2 and i3)
1780 : !! real(F(i1,i3,i2,idat))=ZF(1,i1,i3,i2,idat)
1781 : !! imag(F(i1,i3,i2,idat))=ZF(2,i1,i3,i2,idat)
1782 : !! max1 is positive or zero ; m1 >=max1+1
1783 : !! i1= 1... max1+1 corresponds to positive and zero wavevectors 0 ... max1
1784 : !! then, if m1 > max1+1, one has min1=max1-m1+1 and
1785 : !! i1= max1+2 ... m1 corresponds to negative wavevectors min1 ... -1
1786 : !! i2 and i3 have a similar definition of range
1787 : !! idat=1,ndat
1788 : !! md1,md2,md3: Dimension of ZF (input as well as output), distributed on different procs
1789 : !! md2proc=((md2-1)/nproc_fft)+1 maximal number of small box 2nd dim slices for one proc
1790 : !!
1791 : !! POT: Potential
1792 : !! POT(cplex*i1,i2,i3)
1793 : !! cplex=1 or 2 , i1=1,n1 , i2=1,n2 , i3=1,n3
1794 : !! nd1,nd2,nd3: dimension of pot
1795 : !! comm_fft: MPI communicator
1796 : !! nproc_fft: number of processors used as returned by MPI_COMM_SIZE
1797 : !! me_fft: [0:nproc_fft-1] number of processor as returned by MPI_COMM_RANK
1798 : !! n1,n2,n3: logical dimension of the transform. As transform lengths
1799 : !! most products of the prime factors 2,3,5 are allowed.
1800 : !! The detailed table with allowed transform lengths can
1801 : !! be found in subroutine CTRIG
1802 : !!
1803 : !! NOTES:
1804 : !! PERFORMANCE CONSIDERATIONS:
1805 : !! The maximum number of processors that can reasonably be used is max(n2/2,n3/2)
1806 : !!
1807 : !! It is very important to find the optimal
1808 : !! value of NCACHE. NCACHE determines the size of the work array ZW, that
1809 : !! has to fit into cache. It has therefore to be chosen to equal roughly
1810 : !! half the size of the physical cache in units of real*8 numbers.
1811 : !! The optimal value of ncache can easily be determined by numerical
1812 : !! experimentation. A too large value of ncache leads to a dramatic
1813 : !! and sudden decrease of performance, a too small value to a to a
1814 : !! slow and less dramatic decrease of performance. If NCACHE is set
1815 : !! to a value so small, that not even a single one dimensional transform
1816 : !! can be done in the workarray zw, the program stops with an error message.
1817 : !!
1818 : !! SOURCE
1819 :
1820 :
1821 680 : subroutine sg2002_applypot_many(icplexwf,cplex,ndat,n1,n2,n3,nd1,nd2,nd3,nd3proc,&
1822 : & max1i,max2i,max3i,m1i,m2i,m3i,md1,md2proc,md3,&
1823 680 : & max1o,max2o,max3o,m1o,m2o,m3o,comm_fft,nproc_fft,me_fft,pot,zf)
1824 :
1825 : !Arguments ------------------------------------
1826 : integer,intent(in) :: icplexwf,cplex,ndat,n1,n2,n3,nd1,nd2,nd3,nd3proc
1827 : integer,intent(in) :: max1i,max2i,max3i,m1i,m2i,m3i,md1,md2proc,md3
1828 : integer,intent(in) :: max1o,max2o,max3o,m1o,m2o,m3o,comm_fft,nproc_fft,me_fft
1829 : real(dp),intent(in) :: pot(cplex*nd1,nd2,nd3)
1830 : real(dp),intent(inout) :: zf(2,md1,md3,md2proc,ndat)
1831 :
1832 : !Local variables-------------------------------
1833 : !scalars
1834 : integer,parameter :: unused0=0
1835 : integer :: i,j,i1,i2,i3,ic1,ic2,ic3,idat,ierr,inzee,j3glob
1836 : integer :: ioption,j2,j3,lot,lzt,m1zt,ma,mb,n1dfft,nnd3
1837 : integer :: m2eff,ncache,n1eff,i1inv,i2inv,i3inv,jeff,includelast,j2stb
1838 : integer :: jx,j2stf,Jp2stb,Jp2stf,m2ieff,m2oeff
1839 : !arrays
1840 1360 : integer :: requests(ndat)
1841 : real(dp) :: tsec(2)
1842 680 : real(dp), allocatable :: zt(:,:,:) ! work arrays for transpositions
1843 680 : real(dp), allocatable :: zmpi1(:,:,:,:,:),zmpi2(:,:,:,:,:) ! work arrays for MPI
1844 680 : real(dp), allocatable :: zw(:,:,:) ! cache work array
1845 : ! FFT work arrays
1846 680 : real(dp), allocatable, dimension(:,:) :: btrig1,btrig2,btrig3
1847 680 : real(dp), allocatable, dimension(:,:) :: ftrig1,ftrig2,ftrig3
1848 680 : integer, allocatable, dimension(:) :: after1,now1,before1,after2,now2,before2,after3,now3,before3
1849 :
1850 : ! *************************************************************************
1851 :
1852 : !ioption=0 ! This was in the old version.
1853 680 : ioption=1 ! This one is needed to be compatible with paral_kgb
1854 :
1855 : ! call timab(541,1,tsec)
1856 680 : ncache=4*max(n1,n2,n3,1024)
1857 680 : if (ncache/(4*max(n1,n2,n3)) < 1) then
1858 : write(std_out,*) &
1859 0 : & 'ncache has to be enlarged to be able to hold at', &
1860 0 : & 'least one 1-d FFT of each size even though this will', &
1861 0 : & 'reduce the performance for shorter transform lengths'
1862 0 : ABI_ERROR("Aborting now")
1863 : end if
1864 :
1865 : ! Effective m1 and m2 (complex-to-complex or real-to-complex)
1866 680 : n1eff=n1; m2ieff=m2i; m2oeff=m2o; m1zt=n1
1867 680 : if (icplexwf==1) then
1868 0 : n1eff=(n1+1)/2; m2ieff=m2i/2+1; m2oeff=m2o/2+1; m1zt=2*(n1/2+1)
1869 : end if
1870 :
1871 680 : m2eff=max(m2ieff,m2oeff)
1872 680 : lzt=m2eff
1873 680 : if (mod(m2eff,2) == 0) lzt=lzt+1
1874 336 : if (mod(m2eff,4) == 0) lzt=lzt+1
1875 :
1876 : ! maximal number of big box 3rd dim slices for all procs
1877 680 : nnd3=nd3proc*nproc_fft
1878 :
1879 2040 : ABI_MALLOC(btrig1,(2,n1))
1880 1360 : ABI_MALLOC(ftrig1,(2,n1))
1881 680 : ABI_MALLOC(after1,(mdata))
1882 680 : ABI_MALLOC(now1,(mdata))
1883 680 : ABI_MALLOC(before1,(mdata))
1884 2040 : ABI_MALLOC(btrig2,(2,n2))
1885 1360 : ABI_MALLOC(ftrig2,(2,n2))
1886 680 : ABI_MALLOC(after2,(mdata))
1887 680 : ABI_MALLOC(now2,(mdata))
1888 680 : ABI_MALLOC(before2,(mdata))
1889 2040 : ABI_MALLOC(btrig3,(2,n3))
1890 1360 : ABI_MALLOC(ftrig3,(2,n3))
1891 680 : ABI_MALLOC(after3,(mdata))
1892 680 : ABI_MALLOC(now3,(mdata))
1893 680 : ABI_MALLOC(before3,(mdata))
1894 :
1895 2720 : ABI_MALLOC(zw,(2,ncache/4,2))
1896 2720 : ABI_MALLOC(zt,(2,lzt,m1zt))
1897 4080 : ABI_MALLOC(zmpi2,(2,md1,md2proc,nnd3,ndat))
1898 680 : if (nproc_fft > 1) then
1899 1680 : ABI_MALLOC(zmpi1,(2,md1,md2proc,nnd3,ndat))
1900 : end if
1901 :
1902 680 : call ctrig(n3,btrig3,after3,before3,now3,1,ic3)
1903 680 : call ctrig(n1,btrig1,after1,before1,now1,1,ic1)
1904 680 : call ctrig(n2,btrig2,after2,before2,now2,1,ic2)
1905 :
1906 13656 : do j=1,n1
1907 12976 : ftrig1(1,j)= btrig1(1,j)
1908 13656 : ftrig1(2,j)=-btrig1(2,j)
1909 : end do
1910 13656 : do j=1,n2
1911 12976 : ftrig2(1,j)= btrig2(1,j)
1912 13656 : ftrig2(2,j)=-btrig2(2,j)
1913 : end do
1914 13656 : do j=1,n3
1915 12976 : ftrig3(1,j)= btrig3(1,j)
1916 13656 : ftrig3(2,j)=-btrig3(2,j)
1917 : end do
1918 :
1919 : ! Here we take advantage of non-blocking IALLTOALL:
1920 : ! Perform the first step of MPI-FFT for ndat wavefunctions.
1921 3732 : do idat=1,ndat
1922 :
1923 : !
1924 : ! transform along z axis
1925 : ! input: G1,G3,G2,(Gp2)
1926 3052 : lot=ncache/(4*n3)
1927 23672 : do j2=1,md2proc
1928 23672 : if (me_fft*md2proc+j2 <= m2ieff) then
1929 49164 : do i1=1,m1i,lot
1930 28544 : ma=i1
1931 28544 : mb=min(i1+(lot-1),m1i)
1932 28544 : n1dfft=mb-ma+1
1933 :
1934 : ! zero-pad n1dfft G_z lines
1935 : ! input: G1,G3,G2,(Gp2)
1936 28544 : call fill_cent(md1,md3,lot,n1dfft,max3i,m3i,n3,zf(1,i1,1,j2,idat),zw(1,1,1))
1937 :
1938 28544 : inzee=1
1939 95024 : do i=1,ic3
1940 : call fftstp(lot,n1dfft,n3,lot,n3,zw(1,1,inzee),zw(1,1,3-inzee), &
1941 66480 : & btrig3,after3(i),now3(i),before3(i),1)
1942 95024 : inzee=3-inzee
1943 : end do
1944 :
1945 : ! Local rotation.
1946 : ! input: G1,R3,G2,(Gp2)
1947 : ! output: G1,G2,R3,(Gp2)
1948 49164 : call scramble(i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3,zw(1,1,inzee),zmpi2(:,:,:,:,idat))
1949 : end do
1950 : end if
1951 : end do
1952 :
1953 : ! Interprocessor data transposition
1954 : ! input: G1,G2,R3,Rp3,(Gp2)
1955 : ! output: G1,G2,R3,Gp2,(Rp3)
1956 3732 : if (nproc_fft > 1) then
1957 2016 : call timab(543,1,tsec)
1958 : call xmpi_ialltoall(zmpi2(:,:,:,:,idat),2*md1*md2proc*nd3proc,&
1959 2016 : & zmpi1(:,:,:,:,idat),2*md1*md2proc*nd3proc,comm_fft,requests(idat))
1960 2016 : call timab(543,2,tsec)
1961 : end if
1962 : end do ! idat
1963 :
1964 : ! The second step of MPI-FFT
1965 3732 : do idat=1,ndat
1966 : ! Make sure communication is completed.
1967 3052 : if (nproc_fft>1) call xmpi_wait(requests(idat),ierr)
1968 :
1969 42380 : do j3=1,nd3proc
1970 39328 : j3glob = j3 + me_fft*nd3proc
1971 :
1972 42380 : if (me_fft*nd3proc+j3 <= n3) then
1973 39328 : Jp2stb=1; J2stb=1
1974 39328 : Jp2stf=1; J2stf=1
1975 :
1976 : ! transform along x axis
1977 39328 : lot=ncache/(4*n1)
1978 :
1979 95136 : do j=1,m2ieff,lot
1980 55808 : ma=j
1981 55808 : mb=min(j+(lot-1),m2ieff)
1982 55808 : n1dfft=mb-ma+1
1983 :
1984 : ! Zero-pad input.
1985 : ! input: G1,G2,R3,G2,(Rp3)
1986 : ! output: G2,G1,R3,G2,(Rp3)
1987 55808 : if (nproc_fft == 1) then
1988 : call mpiswitch_cent(j3,n1dfft,Jp2stb,J2stb,lot,max1i,md1,m1i,n1,&
1989 37664 : & md2proc,nd3proc,nproc_fft,ioption,zmpi2(:,:,:,:,idat),zw(1,1,1), unused0, unused0, unused0)
1990 : else
1991 : call mpiswitch_cent(j3,n1dfft,Jp2stb,J2stb,lot,max1i,md1,m1i,n1,&
1992 18144 : & md2proc,nd3proc,nproc_fft,ioption,zmpi1(:,:,:,:,idat),zw(1,1,1), unused0, unused0, unused0)
1993 : end if
1994 :
1995 : ! Transform along x
1996 : ! input: G2,G1,R3,(Rp3)
1997 : ! output: G2,R1,R3,(Rp3)
1998 55808 : inzee=1
1999 131136 : do i=1,ic1-1
2000 : call fftstp(lot,n1dfft,n1,lot,n1,zw(1,1,inzee),zw(1,1,3-inzee), &
2001 75328 : & btrig1,after1(i),now1(i),before1(i),1)
2002 131136 : inzee=3-inzee
2003 : end do
2004 :
2005 55808 : i=ic1
2006 : call fftstp(lot,n1dfft,n1,lzt,m1zt,zw(1,1,inzee),zt(1,j,1), &
2007 95136 : & btrig1,after1(i),now1(i),before1(i),1)
2008 : end do
2009 :
2010 : ! Transform along y axis (take into account c2c or c2r case).
2011 : ! Must loop over the full box.
2012 39328 : lot=ncache/(4*n2)
2013 :
2014 39328 : if (icplexwf==1) then
2015 0 : if(mod(lot,2).ne.0)lot=lot-1 ! needed to introduce jeff
2016 : end if
2017 :
2018 113216 : do j=1,n1eff,lot
2019 73888 : ma=j
2020 73888 : mb=min(j+(lot-1),n1eff)
2021 73888 : n1dfft=mb-ma+1
2022 73888 : jeff=j
2023 73888 : includelast=1
2024 :
2025 73888 : if (icplexwf==1) then
2026 0 : jeff=2*j-1
2027 : includelast=1
2028 0 : if (mb==n1eff .and. n1eff*2/=n1) includelast=0
2029 : end if
2030 :
2031 : ! Zero-pad the input.
2032 : ! input: G2,R1,R3,(Rp3)
2033 : ! output: R1,G2,R3,(Rp3)
2034 73888 : if (icplexwf==2) then
2035 73888 : call switch_cent(n1dfft,max2i,m2i,n2,lot,n1,lzt,zt(1,1,jeff),zw(1,1,1))
2036 : else
2037 0 : call switchreal_cent(includelast,n1dfft,max2i,n2,lot,m1zt,lzt,zt(1,1,jeff),zw(1,1,1))
2038 : end if
2039 :
2040 : ! input: R1,G2,R3,(Rp3)
2041 : ! output: R1,R2,R3,(Rp3)
2042 73888 : inzee=1
2043 259264 : do i=1,ic2
2044 : call fftstp(lot,n1dfft,n2,lot,n2,zw(1,1,inzee),zw(1,1,3-inzee), &
2045 185376 : & btrig2,after2(i),now2(i),before2(i),1)
2046 259264 : inzee=3-inzee
2047 : end do
2048 : ! output: R1,R2,R3,(Rp3)
2049 :
2050 : ! Multiply with potential in real space
2051 73888 : jx=cplex*(jeff-1)+1
2052 73888 : call multpot(icplexwf,cplex,includelast,nd1,nd2,n2,lot,n1dfft,pot(jx,1,j3glob),zw(1,1,inzee))
2053 :
2054 : ! TRANSFORM BACK IN FOURIER SPACE
2055 : ! transform along y axis
2056 : ! input: R1,R2,R3,(Rp3)
2057 259264 : do i=1,ic2
2058 : call fftstp(lot,n1dfft,n2,lot,n2,zw(1,1,inzee),zw(1,1,3-inzee), &
2059 185376 : & ftrig2,after2(i),now2(i),before2(i),-1)
2060 259264 : inzee=3-inzee
2061 : end do
2062 :
2063 : ! input: R1,G2,R3,(Rp3)
2064 : ! output: G2,R1,R3,(Rp3)
2065 113216 : if (icplexwf==2) then
2066 73888 : call unswitch_cent(n1dfft,max2o,m2o,n2,lot,n1,lzt,zw(1,1,inzee),zt(1,1,jeff))
2067 : else
2068 0 : call unswitchreal_cent(n1dfft,max2o,n2,lot,n1,lzt,zw(1,1,inzee),zt(1,1,jeff))
2069 : end if
2070 :
2071 : end do ! j
2072 :
2073 : ! transform along x axis
2074 : ! input: R2,R1,R3,(Rp3)
2075 : ! output: R2,G1,R3,(Rp3)
2076 39328 : lot=ncache/(4*n1)
2077 :
2078 95136 : do j=1,m2oeff,lot
2079 55808 : ma=j
2080 55808 : mb=min(j+(lot-1),m2oeff)
2081 55808 : n1dfft=mb-ma+1
2082 55808 : i=1
2083 : call fftstp(lzt,n1dfft,m1zt,lot,n1,zt(1,j,1),zw(1,1,1), &
2084 55808 : & ftrig1,after1(i),now1(i),before1(i),-1)
2085 :
2086 55808 : inzee=1
2087 131136 : do i=2,ic1
2088 : call fftstp(lot,n1dfft,n1,lot,n1,zw(1,1,inzee),zw(1,1,3-inzee), &
2089 75328 : & ftrig1,after1(i),now1(i),before1(i),-1)
2090 131136 : inzee=3-inzee
2091 : end do
2092 :
2093 : ! input: G2,G1,R3,Gp2,(Rp3)
2094 : ! output: G1,G2,R3,Gp2,(Rp3)
2095 95136 : if (nproc_fft == 1) then
2096 : call unmpiswitch_cent(j3,n1dfft,Jp2stf,J2stf,lot,max1o,md1,m1o,n1,&
2097 37664 : & md2proc,nd3proc,nproc_fft,ioption,zw(1,1,inzee),zmpi2(:,:,:,:,idat))
2098 : else
2099 : call unmpiswitch_cent(j3,n1dfft,Jp2stf,J2stf,lot,max1o,md1,m1o,n1,&
2100 18144 : & md2proc,nd3proc,nproc_fft,ioption,zw(1,1,inzee),zmpi1(:,:,:,:,idat))
2101 : end if
2102 : end do ! j
2103 : end if
2104 : end do
2105 :
2106 : ! Interprocessor data transposition
2107 : ! input: G1,G2,R3,Gp2,(Rp3)
2108 : ! output: G1,G2,R3,Rp3,(Gp2)
2109 3732 : if (nproc_fft > 1) then
2110 2016 : call timab(544,1,tsec)
2111 : call xmpi_ialltoall(zmpi1(:,:,:,:,idat),2*md1*md2proc*nd3proc, &
2112 2016 : & zmpi2(:,:,:,:,idat),2*md1*md2proc*nd3proc,comm_fft,requests(idat))
2113 2016 : call timab(544,2,tsec)
2114 : end if
2115 : end do ! idat
2116 :
2117 3732 : do idat=1,ndat
2118 3052 : if (nproc_fft>1) call xmpi_wait(requests(idat),ierr)
2119 :
2120 : ! transform along z axis
2121 : ! input: G1,G2,R3,(Gp2)
2122 3052 : lot=ncache/(4*n3)
2123 23672 : do j2=1,md2proc
2124 23672 : if (me_fft*md2proc+j2 <= m2oeff) then
2125 49164 : do i1=1,m1o,lot
2126 28544 : ma=i1
2127 28544 : mb=min(i1+(lot-1),m1o)
2128 28544 : n1dfft=mb-ma+1
2129 :
2130 : ! input: G1,G2,R3,(Gp2)
2131 : ! output: G1,R3,G2,(Gp2)
2132 28544 : call unscramble(i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3,zmpi2(:,:,:,:,idat),zw(1,1,1))
2133 :
2134 28544 : inzee=1
2135 95024 : do i=1,ic3
2136 : call fftstp(lot,n1dfft,n3,lot,n3,zw(1,1,inzee),zw(1,1,3-inzee), &
2137 66480 : & ftrig3,after3(i),now3(i),before3(i),-1)
2138 95024 : inzee=3-inzee
2139 : end do
2140 :
2141 49164 : call unfill_cent(md1,md3,lot,n1dfft,max3o,m3o,n3,zw(1,1,inzee),zf(1,i1,1,j2,idat))
2142 : ! output: G1,G3,G2,(Gp2)
2143 : end do
2144 : end if
2145 : end do
2146 :
2147 : ! Complete missing values with complex conjugate
2148 : ! Inverse of ix is located at nx+2-ix , except for ix=1, for which it is 1.
2149 3732 : if (icplexwf==1) then
2150 0 : do i3=1,m3o
2151 0 : i3inv=m3o+2-i3
2152 0 : if (i3==1) i3inv=1
2153 0 : if (m2oeff>1)then
2154 0 : do i2=2,m2oeff
2155 0 : i2inv=m2o+2-i2
2156 0 : zf(1,1,i3inv,i2inv,idat)= zf(1,1,i3,i2,idat)
2157 0 : zf(2,1,i3inv,i2inv,idat)=-zf(2,1,i3,i2,idat)
2158 0 : do i1=2,m1o
2159 0 : i1inv=m1o+2-i1
2160 0 : zf(1,i1inv,i3inv,i2inv,idat)= zf(1,i1,i3,i2,idat)
2161 0 : zf(2,i1inv,i3inv,i2inv,idat)=-zf(2,i1,i3,i2,idat)
2162 : end do
2163 : end do
2164 : end if
2165 : end do
2166 : end if
2167 :
2168 : end do ! idat
2169 :
2170 680 : ABI_FREE(btrig1)
2171 680 : ABI_FREE(ftrig1)
2172 680 : ABI_FREE(after1)
2173 680 : ABI_FREE(now1)
2174 680 : ABI_FREE(before1)
2175 680 : ABI_FREE(btrig2)
2176 680 : ABI_FREE(ftrig2)
2177 680 : ABI_FREE(after2)
2178 680 : ABI_FREE(now2)
2179 680 : ABI_FREE(before2)
2180 680 : ABI_FREE(btrig3)
2181 680 : ABI_FREE(ftrig3)
2182 680 : ABI_FREE(after3)
2183 680 : ABI_FREE(now3)
2184 680 : ABI_FREE(before3)
2185 :
2186 680 : ABI_FREE(zmpi2)
2187 680 : ABI_FREE(zw)
2188 680 : ABI_FREE(zt)
2189 680 : if (nproc_fft > 1) then
2190 336 : ABI_FREE(zmpi1)
2191 : end if
2192 :
2193 680 : end subroutine sg2002_applypot_many
2194 : !!***
2195 :
2196 : !----------------------------------------------------------------------
2197 :
2198 : !!****f* m_sg2002/sg2002_accrho
2199 : !! NAME
2200 : !! sg2002_accrho
2201 : !!
2202 : !! FUNCTION
2203 : !! Accumulates the real space density rho from the ndat wavefunctions zf
2204 : !! by transforming zf into real space and adding all the amplitudes squared
2205 : !!
2206 : !! INPUTS:
2207 : !! ZF: input array (note the switch of i2 and i3)
2208 : !! real(F(i1,i3,i2,idat))=ZF(1,i1,i3,i2,idat)
2209 : !! imag(F(i1,i3,i2,idat))=ZF(2,i1,i3,i2,idat)
2210 : !! max1 is positive or zero ; m1 >=max1+1
2211 : !! i1= 1... max1+1 corresponds to positive and zero wavevectors 0 ... max1
2212 : !! then, if m1 > max1+1, one has min1=max1-m1+1 and
2213 : !! i1= max1+2 ... m1 corresponds to negative wavevectors min1 ... -1
2214 : !! i2 and i3 have a similar definition of range
2215 : !! idat=1,ndat
2216 : !! md1,md2,md3: Dimension of ZF
2217 : !! md2proc=((md2-1)/nproc_fft)+1 ! maximal number of small box 2nd dim slices for one proc
2218 : !! weight(ndat)= weight for the density accumulation
2219 : !!
2220 : !! OUTPUTS:
2221 : !! RHOoutput(i1,i2,i3) = RHOinput(i1,i2,i3) + sum on idat of (Re(FFT(ZF))**2 *weight_r + weight_i*Im(FFT(ZF))**2
2222 : !! i1=1,n1 , i2=1,n2 , i3=1,n3
2223 : !! comm_fft: MPI communicator
2224 : !! nproc_fft: number of processors used as returned by MPI_COMM_SIZE
2225 : !! me_fft: [0:nproc_fft-1] number of processor as returned by MPI_COMM_RANK
2226 : !! n1,n2,n3: logical dimension of the transform. As transform lengths
2227 : !! most products of the prime factors 2,3,5 are allowed.
2228 : !! The detailed table with allowed transform lengths can
2229 : !! be found in subroutine CTRIG
2230 : !! nd1,nd2,nd3: Dimension of RHO
2231 : !! nd3proc=((nd3-1)/nproc_fft)+1 ! maximal number of big box 3rd dim slices for one proc
2232 : !!
2233 : !! NOTES:
2234 : !! PERFORMANCE CONSIDERATIONS:
2235 : !! The maximum number of processors that can reasonably be used is max(n2/2,n3/2)
2236 : !!
2237 : !! It is very important to find the optimal
2238 : !! value of NCACHE. NCACHE determines the size of the work array ZW, that
2239 : !! has to fit into cache. It has therefore to be chosen to equal roughly
2240 : !! half the size of the physical cache in units of real*8 numbers.
2241 : !! The optimal value of ncache can easily be determined by numerical
2242 : !! experimentation. A too large value of ncache leads to a dramatic
2243 : !! and sudden decrease of performance, a too small value to a to a
2244 : !! slow and less dramatic decrease of performance. If NCACHE is set
2245 : !! to a value so small, that not even a single one dimensional transform
2246 : !! can be done in the workarray zw, the program stops with an error message.
2247 : !!
2248 : !! SOURCE
2249 :
2250 214 : subroutine sg2002_accrho(icplexwf,ndat,n1,n2,n3,nd1,nd2,nd3,nd3proc,&
2251 214 : & max1,max2,max3,m1,m2,m3,md1,md2proc,md3,comm_fft,nproc_fft,me_fft,zf,rho,weight_r,weight_i)
2252 :
2253 : !Arguments ------------------------------------
2254 : integer,intent(in) :: icplexwf,ndat,n1,n2,n3,nd1,nd2,nd3,nd3proc
2255 : integer,intent(in) :: max1,max2,max3,m1,m2,m3,md1,md2proc,md3,comm_fft,nproc_fft,me_fft
2256 : real(dp),intent(in) :: zf(2,md1,md3,md2proc,ndat)
2257 : real(dp),intent(in) :: weight_r(ndat), weight_i(ndat)
2258 : real(dp),intent(inout) :: rho(nd1,nd2,nd3)
2259 :
2260 : !Local variables-------------------------------
2261 : !scalars
2262 : integer,parameter :: unused0=0
2263 : integer :: i,j,i1,ic1,ic2,ic3,idat,ierr,inzee,j3glob
2264 : integer :: ioption,j2,j3,j2st,jp2st,lot,lzt,m1zt,ma,mb,n1dfft,nnd3
2265 : integer :: m2eff,ncache,n1eff,jeff,includelast
2266 : !arrays
2267 214 : real(dp), allocatable :: zmpi1(:,:,:,:),zmpi2(:,:,:,:) ! work arrays for MPI
2268 214 : real(dp), allocatable :: zt(:,:,:) ! work arrays for transpositions
2269 214 : real(dp), allocatable :: zw(:,:,:) ! cache work array
2270 : real(dp) :: tsec(2)
2271 : ! FFT work arrays
2272 214 : real(dp), allocatable, dimension(:,:) :: trig1,trig2,trig3
2273 214 : integer, allocatable, dimension(:) :: after1,now1,before1, after2,now2,before2,after3,now3,before3
2274 :
2275 : ! *************************************************************************
2276 :
2277 : !ioption=0 ! This was in the old version.
2278 214 : ioption=1 ! This one is needed to be compatible with paral_kgb
2279 :
2280 : !nproc_fft = xmpi_comm_size(comm_fft); me_fft = xmpi_comm_rank(comm_fft)
2281 :
2282 : ! find cache size that gives optimal performance on machine
2283 214 : ncache=4*max(n1,n2,n3,1024)
2284 214 : if (ncache/(4*max(n1,n2,n3)) < 1) then
2285 : write(std_out,*) &
2286 0 : & 'ncache has to be enlarged to be able to hold at', &
2287 0 : & 'least one 1-d FFT of each size even though this will', &
2288 0 : & 'reduce the performance for shorter transform lengths'
2289 0 : ABI_ERROR("Aborting now")
2290 : end if
2291 :
2292 : !Effective m1 and m2 (complex-to-complex or real-to-complex)
2293 214 : n1eff=n1; m2eff=m2 ; m1zt=n1
2294 214 : if (icplexwf==1) then
2295 16 : n1eff=(n1+1)/2; m2eff=m2/2+1; m1zt=2*(n1/2+1)
2296 : end if
2297 :
2298 214 : lzt=m2eff
2299 214 : if (mod(m2eff,2) == 0) lzt=lzt+1
2300 106 : if (mod(m2eff,4) == 0) lzt=lzt+1
2301 :
2302 : ! maximal number of big box 3rd dim slices for all procs
2303 214 : nnd3=nd3proc*nproc_fft
2304 :
2305 642 : ABI_MALLOC(trig1,(2,n1))
2306 214 : ABI_MALLOC(after1,(mdata))
2307 214 : ABI_MALLOC(now1,(mdata))
2308 214 : ABI_MALLOC(before1,(mdata))
2309 642 : ABI_MALLOC(trig2,(2,n2))
2310 214 : ABI_MALLOC(after2,(mdata))
2311 214 : ABI_MALLOC(now2,(mdata))
2312 214 : ABI_MALLOC(before2,(mdata))
2313 642 : ABI_MALLOC(trig3,(2,n3))
2314 214 : ABI_MALLOC(after3,(mdata))
2315 214 : ABI_MALLOC(now3,(mdata))
2316 214 : ABI_MALLOC(before3,(mdata))
2317 :
2318 856 : ABI_MALLOC(zw,(2,ncache/4,2))
2319 856 : ABI_MALLOC(zt,(2,lzt,m1zt))
2320 1070 : ABI_MALLOC(zmpi2,(2,md1,md2proc,nnd3))
2321 214 : if (nproc_fft > 1) then
2322 256 : ABI_MALLOC(zmpi1,(2,md1,md2proc,nnd3))
2323 : end if
2324 :
2325 214 : call ctrig(n3,trig3,after3,before3,now3,1,ic3)
2326 214 : call ctrig(n1,trig1,after1,before1,now1,1,ic1)
2327 214 : call ctrig(n2,trig2,after2,before2,now2,1,ic2)
2328 :
2329 438 : do idat=1,ndat
2330 : ! transform along z axis
2331 : ! input: I1,I3,J2,(Jp2)
2332 224 : lot=ncache/(4*n3)
2333 :
2334 : ! Loop over the y planes treated by this node and trasform n1ddft G_z lines.
2335 4548 : do j2=1,md2proc
2336 4548 : if (me_fft*md2proc+j2 <= m2eff) then ! MG REMOVED TO BE COSISTENT WITH BACK_WF
2337 17538 : do i1=1,m1,lot
2338 13444 : ma=i1
2339 13444 : mb=min(i1+(lot-1),m1)
2340 13444 : n1dfft=mb-ma+1
2341 :
2342 : ! zero-pad n1dfft G_z lines
2343 : ! input: G1,G3,G2,(Gp2)
2344 : ! output: G1,R3,G2,(Gp2)
2345 13444 : call fill_cent(md1,md3,lot,n1dfft,max3,m3,n3,zf(1,i1,1,j2,idat),zw(1,1,1))
2346 :
2347 : ! Transform along z.
2348 13444 : inzee=1
2349 52886 : do i=1,ic3
2350 : call fftstp(lot,n1dfft,n3,lot,n3,zw(1,1,inzee),zw(1,1,3-inzee), &
2351 39442 : & trig3,after3(i),now3(i),before3(i),1)
2352 52886 : inzee=3-inzee
2353 : end do
2354 :
2355 : ! Local rotation.
2356 : ! input: G1,R3,G2,(Gp2)
2357 : ! output: G1,G2,R3,(Gp2)
2358 17538 : call scramble(i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3,zw(1,1,inzee),zmpi2)
2359 : end do
2360 : end if
2361 : end do
2362 :
2363 : ! Interprocessor data transposition
2364 : ! input: G1,G2,R3,Rp3,(Gp2)
2365 : ! output: G1,G2,R3,Gp2,(Rp3)
2366 224 : if (nproc_fft > 1) then
2367 64 : call timab(543,1,tsec)
2368 : call xmpi_alltoall(zmpi2,2*md1*md2proc*nd3proc, &
2369 64 : & zmpi1,2*md1*md2proc*nd3proc,comm_fft,ierr)
2370 64 : call timab(543,2,tsec)
2371 : end if
2372 :
2373 : ! Loop over the z treated by this node.
2374 9074 : do j3=1,nd3proc
2375 8636 : j3glob = j3 + me_fft*nd3proc
2376 : !ABI_CHECK(j3glob <= n3, "j3glob")
2377 :
2378 8860 : if (me_fft*nd3proc+j3 <= n3) then
2379 8636 : Jp2st=1; J2st=1
2380 :
2381 8636 : lot=ncache/(4*n1)
2382 :
2383 : ! Loop over G_y in the small box.
2384 36984 : do j=1,m2eff,lot
2385 28348 : ma=j
2386 28348 : mb=min(j+(lot-1),m2eff)
2387 28348 : n1dfft=mb-ma+1
2388 :
2389 : ! Zero-pad input.
2390 : ! input: G1,G2,R3,JG2,(Rp3)
2391 : ! output: G2,G1,R3,JG2,(Rp3)
2392 :
2393 28348 : if (nproc_fft == 1) then
2394 : call mpiswitch_cent(j3,n1dfft,Jp2st,J2st,lot,max1,md1,m1,n1,&
2395 27772 : & md2proc,nd3proc,nproc_fft,ioption,zmpi2,zw(1,1,1),unused0, unused0,unused0)
2396 : else
2397 : call mpiswitch_cent(j3,n1dfft,Jp2st,J2st,lot,max1,md1,m1,n1,&
2398 576 : & md2proc,nd3proc,nproc_fft,ioption,zmpi1,zw(1,1,1), unused0,unused0,unused0)
2399 : end if
2400 :
2401 : ! Transform along x
2402 : ! input: G2,G1,R3,(Rp3)
2403 : ! output: G2,R1,R3,(Rp3)
2404 28348 : inzee=1
2405 83160 : do i=1,ic1-1
2406 : call fftstp(lot,n1dfft,n1,lot,n1,zw(1,1,inzee),zw(1,1,3-inzee), &
2407 54812 : & trig1,after1(i),now1(i),before1(i),1)
2408 83160 : inzee=3-inzee
2409 : end do
2410 :
2411 28348 : i=ic1
2412 : call fftstp(lot,n1dfft,n1,lzt,m1zt,zw(1,1,inzee),zt(1,j,1), &
2413 36984 : & trig1,after1(i),now1(i),before1(i),1)
2414 : end do
2415 :
2416 : ! Transform along y axis (take into account c2c or c2r case).
2417 : ! Must loop over the full box.
2418 8636 : lot=ncache/(4*n2)
2419 8636 : if (icplexwf==1) then
2420 1056 : if (mod(lot,2).ne.0) lot=lot-1 ! needed to introduce jeff
2421 : end if
2422 :
2423 62380 : do j=1,n1eff,lot
2424 53744 : ma=j
2425 53744 : mb=min(j+(lot-1),n1eff)
2426 53744 : n1dfft=mb-ma+1
2427 53744 : jeff=j
2428 53744 : includelast=1
2429 :
2430 53744 : if (icplexwf==1) then
2431 3856 : jeff=2*j-1
2432 : includelast=1
2433 3856 : if (mb==n1eff .and. n1eff*2/=n1) includelast=0
2434 : end if
2435 :
2436 : ! Zero-pad the input.
2437 : ! input: G2,R1,R3,(Rp3)
2438 : ! output: R1,G2,R3,(Rp3)
2439 53744 : if (icplexwf==2) then
2440 49888 : call switch_cent(n1dfft,max2,m2,n2,lot,n1,lzt,zt(1,1,j),zw(1,1,1))
2441 : else
2442 3856 : call switchreal_cent(includelast,n1dfft,max2,n2,lot,m1zt,lzt,zt(1,1,jeff),zw(1,1,1))
2443 : end if
2444 :
2445 53744 : inzee=1
2446 213092 : do i=1,ic2
2447 : call fftstp(lot,n1dfft,n2,lot,n2,zw(1,1,inzee),zw(1,1,3-inzee), &
2448 159348 : & trig2,after2(i),now2(i),before2(i),1)
2449 213092 : inzee=3-inzee
2450 : end do
2451 :
2452 : ! Accumulate
2453 : call addrho(icplexwf,includelast,nd1,nd2,n2,lot,n1dfft,&
2454 62380 : & zw(1,1,inzee),rho(jeff,1,j3glob),weight_r(idat),weight_i(idat))
2455 : end do
2456 : ! output: i1,i2,j3,(jp3)
2457 :
2458 : end if
2459 : end do ! j3
2460 : end do ! idat
2461 :
2462 214 : ABI_FREE(trig1)
2463 214 : ABI_FREE(after1)
2464 214 : ABI_FREE(now1)
2465 214 : ABI_FREE(before1)
2466 214 : ABI_FREE(trig2)
2467 214 : ABI_FREE(after2)
2468 214 : ABI_FREE(now2)
2469 214 : ABI_FREE(before2)
2470 214 : ABI_FREE(trig3)
2471 214 : ABI_FREE(after3)
2472 214 : ABI_FREE(now3)
2473 214 : ABI_FREE(before3)
2474 :
2475 214 : ABI_FREE(zmpi2)
2476 214 : ABI_FREE(zw)
2477 214 : ABI_FREE(zt)
2478 214 : if (nproc_fft > 1) then
2479 64 : ABI_FREE(zmpi1)
2480 : end if
2481 :
2482 214 : end subroutine sg2002_accrho
2483 : !!***
2484 :
2485 : !!****f* m_sg2002/ctrig
2486 : !! NAME
2487 : !! ctrig
2488 : !!
2489 : !! FUNCTION
2490 : !!
2491 : !! INPUTS
2492 : !!
2493 : !! OUTPUT
2494 : !!
2495 : !! SOURCE
2496 :
2497 5024478 : subroutine ctrig(n,trig,after,before,now,isign,ic)
2498 :
2499 : !Arguments ------------------------------------
2500 : integer,intent(in) :: n,isign
2501 : integer,intent(inout) :: ic
2502 : integer,intent(inout) :: after(mdata),before(mdata),now(mdata)
2503 : real(dp),intent(inout) :: trig(2,n)
2504 :
2505 : !Local variables-------------------------------
2506 : !scalars
2507 : integer :: i,itt,j,nh
2508 : real(dp) :: angle,trigc,trigs
2509 :
2510 : ! *************************************************************************
2511 :
2512 58820902 : do i=1,ndata
2513 58820902 : if (n.eq.ifftdata(1,i)) then
2514 5024478 : ic=0
2515 15149135 : do j=1,(mdata-1)
2516 15149135 : itt=ifftdata(1+j,i)
2517 15149135 : if (itt.gt.1) then
2518 10124657 : ic=ic+1
2519 10124657 : now(j)=ifftdata(1+j,i)
2520 : else
2521 : goto 1000
2522 : end if
2523 : end do
2524 : goto 1000
2525 : end if
2526 : end do
2527 :
2528 0 : write(std_out,*) 'VALUE OF',n,'NOT ALLOWED FOR FFT, ALLOWED VALUES ARE:'
2529 : 37 format(15(i5))
2530 0 : write(std_out,37) (ifftdata(1,i),i=1,ndata)
2531 5024478 : ABI_ERROR("Aborting now")
2532 :
2533 : 1000 continue
2534 5024478 : after(1)=1
2535 5024478 : before(ic)=1
2536 10124657 : do i=2,ic
2537 5100179 : after(i)=after(i-1)*now(i-1)
2538 10124657 : before(ic-i+1)=before(ic-i+2)*now(ic-i+2)
2539 : end do
2540 :
2541 5024478 : angle=isign*two_pi/n
2542 5024478 : if (mod(n,2).eq.0) then
2543 4968979 : nh=n/2
2544 4968979 : trig(1,1)=one
2545 4968979 : trig(2,1)=zero
2546 4968979 : trig(1,nh+1)=-one
2547 4968979 : trig(2,nh+1)=zero
2548 50016357 : do i=1,nh-1
2549 45047378 : trigc=cos(i*angle)
2550 45047378 : trigs=sin(i*angle)
2551 45047378 : trig(1,i+1)=trigc
2552 45047378 : trig(2,i+1)=trigs
2553 45047378 : trig(1,n-i+1)=trigc
2554 50016357 : trig(2,n-i+1)=-trigs
2555 : end do
2556 : else
2557 55499 : nh=(n-1)/2
2558 55499 : trig(1,1)=one
2559 55499 : trig(2,1)=zero
2560 667156 : do i=1,nh
2561 611657 : trigc=cos(i*angle)
2562 611657 : trigs=sin(i*angle)
2563 611657 : trig(1,i+1)=trigc
2564 611657 : trig(2,i+1)=trigs
2565 611657 : trig(1,n-i+1)=trigc
2566 667156 : trig(2,n-i+1)=-trigs
2567 : end do
2568 : end if
2569 :
2570 5024478 : end subroutine ctrig
2571 : !!***
2572 :
2573 : !!****f* m_sg2002/fftstp
2574 : !! NAME
2575 : !! fftstp
2576 : !!
2577 : !! FUNCTION
2578 : !!
2579 : !! INPUTS
2580 : !! mm
2581 : !! n1dfft
2582 : !! m
2583 : !! nn
2584 : !! n
2585 : !! zin
2586 : !! trig
2587 : !! after
2588 : !! now
2589 : !! before
2590 : !! isign
2591 : !!
2592 : !! OUTPUT
2593 : !! zout
2594 : !!
2595 : !! SOURCE
2596 :
2597 117986786 : subroutine fftstp(mm,n1dfft,m,nn,n,zin,zout,trig,after,now,before,isign)
2598 :
2599 : !Arguments ------------------------------------
2600 : integer,intent(in) :: after,before,mm,n1dfft,m,nn,n,now,isign
2601 : real(dp),intent(in) :: trig(2,n),zin(2,mm,m)
2602 : real(dp),intent(inout) :: zout(2,nn,n)
2603 :
2604 : !Local variables-------------------------------
2605 : integer :: atn,atb,ia,ias,ib,itrig,itt,j,nin1,nin2,nin3,nin4,nin5,nin6,nin7,nin8
2606 : integer :: nout1,nout2,nout3,nout4,nout5,nout6,nout7,nout8
2607 : real(dp) :: am,ap,bm,bp,ci3,ci4,ci5,ci6,ci7,ci8,cm,cos2,cos4,cp,cr2,cr3,cr4,cr5,cr6,cr7,cr8
2608 : real(dp) :: dm,bb,ci2,dpp,r,r2,r25,r3,r34,r4,r5,r6,r7,r8,rt2i,s,r1,s1,s2,s3,s25,s34,s4,s5,s6,s7,s8
2609 : real(dp) :: sin2,ui1,ui2,ui3,ur1,ur2,ur3,sin4,vi1,vi2,vi3,vr1,vr2,vr3
2610 :
2611 : ! *************************************************************************
2612 117986786 : atn=after*now
2613 117986786 : atb=after*before
2614 :
2615 : ! sqrt(.5d0)
2616 117986786 : rt2i=half_sqrt2
2617 : if (now.eq.2) then
2618 496368 : ia=1
2619 496368 : nin1=ia-after
2620 496368 : nout1=ia-atn
2621 992736 : do ib=1,before
2622 496368 : nin1=nin1+after
2623 496368 : nin2=nin1+atb
2624 496368 : nout1=nout1+atn
2625 496368 : nout2=nout1+after
2626 6517536 : do j=1,n1dfft
2627 5524800 : r1=zin(1,j,nin1)
2628 5524800 : s1=zin(2,j,nin1)
2629 5524800 : r2=zin(1,j,nin2)
2630 5524800 : s2=zin(2,j,nin2)
2631 5524800 : zout(1,j,nout1)= r2 + r1
2632 5524800 : zout(2,j,nout1)= s2 + s1
2633 5524800 : zout(1,j,nout2)= r1 - r2
2634 6021168 : zout(2,j,nout2)= s1 - s2
2635 : enddo
2636 : enddo
2637 11864880 : do 2000,ia=2,after
2638 11368512 : ias=ia-1
2639 11368512 : if (2*ias.eq.after) then
2640 0 : if (isign.eq.1) then
2641 0 : nin1=ia-after
2642 0 : nout1=ia-atn
2643 0 : do ib=1,before
2644 0 : nin1=nin1+after
2645 0 : nin2=nin1+atb
2646 0 : nout1=nout1+atn
2647 0 : nout2=nout1+after
2648 0 : do j=1,n1dfft
2649 0 : r1=zin(1,j,nin1)
2650 0 : s1=zin(2,j,nin1)
2651 0 : r2=zin(2,j,nin2)
2652 0 : s2=zin(1,j,nin2)
2653 0 : zout(1,j,nout1)= r1 - r2
2654 0 : zout(2,j,nout1)= s2 + s1
2655 0 : zout(1,j,nout2)= r2 + r1
2656 0 : zout(2,j,nout2)= s1 - s2
2657 : enddo
2658 : enddo
2659 : else
2660 0 : nin1=ia-after
2661 0 : nout1=ia-atn
2662 0 : do ib=1,before
2663 0 : nin1=nin1+after
2664 0 : nin2=nin1+atb
2665 0 : nout1=nout1+atn
2666 0 : nout2=nout1+after
2667 0 : do j=1,n1dfft
2668 0 : r1=zin(1,j,nin1)
2669 0 : s1=zin(2,j,nin1)
2670 0 : r2=zin(2,j,nin2)
2671 0 : s2=zin(1,j,nin2)
2672 0 : zout(1,j,nout1)= r2 + r1
2673 0 : zout(2,j,nout1)= s1 - s2
2674 0 : zout(1,j,nout2)= r1 - r2
2675 0 : zout(2,j,nout2)= s2 + s1
2676 : enddo
2677 : enddo
2678 : end if
2679 11368512 : else if (4*ias.eq.after) then
2680 0 : if (isign.eq.1) then
2681 0 : nin1=ia-after
2682 0 : nout1=ia-atn
2683 0 : do ib=1,before
2684 0 : nin1=nin1+after
2685 0 : nin2=nin1+atb
2686 0 : nout1=nout1+atn
2687 0 : nout2=nout1+after
2688 0 : do j=1,n1dfft
2689 0 : r1=zin(1,j,nin1)
2690 0 : s1=zin(2,j,nin1)
2691 0 : r=zin(1,j,nin2)
2692 0 : s=zin(2,j,nin2)
2693 0 : r2=(r - s)*rt2i
2694 0 : s2=(r + s)*rt2i
2695 0 : zout(1,j,nout1)= r2 + r1
2696 0 : zout(2,j,nout1)= s2 + s1
2697 0 : zout(1,j,nout2)= r1 - r2
2698 0 : zout(2,j,nout2)= s1 - s2
2699 : enddo
2700 : enddo
2701 : else
2702 0 : nin1=ia-after
2703 0 : nout1=ia-atn
2704 0 : do ib=1,before
2705 0 : nin1=nin1+after
2706 0 : nin2=nin1+atb
2707 0 : nout1=nout1+atn
2708 0 : nout2=nout1+after
2709 0 : do j=1,n1dfft
2710 0 : r1=zin(1,j,nin1)
2711 0 : s1=zin(2,j,nin1)
2712 0 : r=zin(1,j,nin2)
2713 0 : s=zin(2,j,nin2)
2714 0 : r2=(r + s)*rt2i
2715 0 : s2=(s - r)*rt2i
2716 0 : zout(1,j,nout1)= r2 + r1
2717 0 : zout(2,j,nout1)= s2 + s1
2718 0 : zout(1,j,nout2)= r1 - r2
2719 0 : zout(2,j,nout2)= s1 - s2
2720 : enddo
2721 : enddo
2722 : end if
2723 11368512 : else if (4*ias.eq.3*after) then
2724 0 : if (isign.eq.1) then
2725 0 : nin1=ia-after
2726 0 : nout1=ia-atn
2727 0 : do ib=1,before
2728 0 : nin1=nin1+after
2729 0 : nin2=nin1+atb
2730 0 : nout1=nout1+atn
2731 0 : nout2=nout1+after
2732 0 : do j=1,n1dfft
2733 0 : r1=zin(1,j,nin1)
2734 0 : s1=zin(2,j,nin1)
2735 0 : r=zin(1,j,nin2)
2736 0 : s=zin(2,j,nin2)
2737 0 : r2=(r + s)*rt2i
2738 0 : s2=(r - s)*rt2i
2739 0 : zout(1,j,nout1)= r1 - r2
2740 0 : zout(2,j,nout1)= s2 + s1
2741 0 : zout(1,j,nout2)= r2 + r1
2742 0 : zout(2,j,nout2)= s1 - s2
2743 : enddo
2744 : enddo
2745 : else
2746 0 : nin1=ia-after
2747 0 : nout1=ia-atn
2748 0 : do ib=1,before
2749 0 : nin1=nin1+after
2750 0 : nin2=nin1+atb
2751 0 : nout1=nout1+atn
2752 0 : nout2=nout1+after
2753 0 : do j=1,n1dfft
2754 0 : r1=zin(1,j,nin1)
2755 0 : s1=zin(2,j,nin1)
2756 0 : r=zin(1,j,nin2)
2757 0 : s=zin(2,j,nin2)
2758 0 : r2=(s - r)*rt2i
2759 0 : s2=(r + s)*rt2i
2760 0 : zout(1,j,nout1)= r2 + r1
2761 0 : zout(2,j,nout1)= s1 - s2
2762 0 : zout(1,j,nout2)= r1 - r2
2763 0 : zout(2,j,nout2)= s2 + s1
2764 : enddo
2765 : enddo
2766 : end if
2767 : else
2768 11368512 : itrig=ias*before+1
2769 11368512 : cr2=trig(1,itrig)
2770 11368512 : ci2=trig(2,itrig)
2771 11368512 : nin1=ia-after
2772 11368512 : nout1=ia-atn
2773 22737024 : do ib=1,before
2774 11368512 : nin1=nin1+after
2775 11368512 : nin2=nin1+atb
2776 11368512 : nout1=nout1+atn
2777 11368512 : nout2=nout1+after
2778 150977664 : do j=1,n1dfft
2779 128240640 : r1=zin(1,j,nin1)
2780 128240640 : s1=zin(2,j,nin1)
2781 128240640 : r=zin(1,j,nin2)
2782 128240640 : s=zin(2,j,nin2)
2783 128240640 : r2=r*cr2 - s*ci2
2784 128240640 : s2=r*ci2 + s*cr2
2785 128240640 : zout(1,j,nout1)= r2 + r1
2786 128240640 : zout(2,j,nout1)= s2 + s1
2787 128240640 : zout(1,j,nout2)= r1 - r2
2788 139609152 : zout(2,j,nout2)= s1 - s2
2789 : enddo
2790 : enddo
2791 : end if
2792 496368 : 2000 continue
2793 : else if (now.eq.4) then
2794 31458863 : if (isign.eq.1) then
2795 16507258 : ia=1
2796 16507258 : nin1=ia-after
2797 16507258 : nout1=ia-atn
2798 68071554 : do ib=1,before
2799 51564296 : nin1=nin1+after
2800 51564296 : nin2=nin1+atb
2801 51564296 : nin3=nin2+atb
2802 51564296 : nin4=nin3+atb
2803 51564296 : nout1=nout1+atn
2804 51564296 : nout2=nout1+after
2805 51564296 : nout3=nout2+after
2806 51564296 : nout4=nout3+after
2807 772519455 : do j=1,n1dfft
2808 704447901 : r1=zin(1,j,nin1)
2809 704447901 : s1=zin(2,j,nin1)
2810 704447901 : r2=zin(1,j,nin2)
2811 704447901 : s2=zin(2,j,nin2)
2812 704447901 : r3=zin(1,j,nin3)
2813 704447901 : s3=zin(2,j,nin3)
2814 704447901 : r4=zin(1,j,nin4)
2815 704447901 : s4=zin(2,j,nin4)
2816 704447901 : r=r1 + r3
2817 704447901 : s=r2 + r4
2818 704447901 : zout(1,j,nout1) = r + s
2819 704447901 : zout(1,j,nout3) = r - s
2820 704447901 : r=r1 - r3
2821 704447901 : s=s2 - s4
2822 704447901 : zout(1,j,nout2) = r - s
2823 704447901 : zout(1,j,nout4) = r + s
2824 704447901 : r=s1 + s3
2825 704447901 : s=s2 + s4
2826 704447901 : zout(2,j,nout1) = r + s
2827 704447901 : zout(2,j,nout3) = r - s
2828 704447901 : r=s1 - s3
2829 704447901 : s=r2 - r4
2830 704447901 : zout(2,j,nout2) = r + s
2831 756012197 : zout(2,j,nout4) = r - s
2832 : enddo
2833 : enddo
2834 63304048 : do 4000,ia=2,after
2835 46796790 : ias=ia-1
2836 46796790 : if (2*ias.eq.after) then
2837 4838434 : nin1=ia-after
2838 4838434 : nout1=ia-atn
2839 12276516 : do ib=1,before
2840 7438082 : nin1=nin1+after
2841 7438082 : nin2=nin1+atb
2842 7438082 : nin3=nin2+atb
2843 7438082 : nin4=nin3+atb
2844 7438082 : nout1=nout1+atn
2845 7438082 : nout2=nout1+after
2846 7438082 : nout3=nout2+after
2847 7438082 : nout4=nout3+after
2848 114174315 : do j=1,n1dfft
2849 101897799 : r1=zin(1,j,nin1)
2850 101897799 : s1=zin(2,j,nin1)
2851 101897799 : r=zin(1,j,nin2)
2852 101897799 : s=zin(2,j,nin2)
2853 101897799 : r2=(r-s)*rt2i
2854 101897799 : s2=(r+s)*rt2i
2855 101897799 : r3=zin(2,j,nin3)
2856 101897799 : s3=zin(1,j,nin3)
2857 101897799 : r=zin(1,j,nin4)
2858 101897799 : s=zin(2,j,nin4)
2859 101897799 : r4=(r + s)*rt2i
2860 101897799 : s4=(r - s)*rt2i
2861 101897799 : r=r1 - r3
2862 101897799 : s=r2 - r4
2863 101897799 : zout(1,j,nout1) = r + s
2864 101897799 : zout(1,j,nout3) = r - s
2865 101897799 : r=r1 + r3
2866 101897799 : s=s2 - s4
2867 101897799 : zout(1,j,nout2) = r - s
2868 101897799 : zout(1,j,nout4) = r + s
2869 101897799 : r=s1 + s3
2870 101897799 : s=s2 + s4
2871 101897799 : zout(2,j,nout1) = r + s
2872 101897799 : zout(2,j,nout3) = r - s
2873 101897799 : r=s1 - s3
2874 101897799 : s=r2 + r4
2875 101897799 : zout(2,j,nout2) = r + s
2876 109335881 : zout(2,j,nout4) = r - s
2877 : enddo
2878 : enddo
2879 : else
2880 41958356 : itt=ias*before
2881 41958356 : itrig=itt+1
2882 41958356 : cr2=trig(1,itrig)
2883 41958356 : ci2=trig(2,itrig)
2884 41958356 : itrig=itrig+itt
2885 41958356 : cr3=trig(1,itrig)
2886 41958356 : ci3=trig(2,itrig)
2887 41958356 : itrig=itrig+itt
2888 41958356 : cr4=trig(1,itrig)
2889 41958356 : ci4=trig(2,itrig)
2890 41958356 : nin1=ia-after
2891 41958356 : nout1=ia-atn
2892 90228240 : do ib=1,before
2893 48269884 : nin1=nin1+after
2894 48269884 : nin2=nin1+atb
2895 48269884 : nin3=nin2+atb
2896 48269884 : nin4=nin3+atb
2897 48269884 : nout1=nout1+atn
2898 48269884 : nout2=nout1+after
2899 48269884 : nout3=nout2+after
2900 48269884 : nout4=nout3+after
2901 718969902 : do j=1,n1dfft
2902 628741662 : r1=zin(1,j,nin1)
2903 628741662 : s1=zin(2,j,nin1)
2904 628741662 : r=zin(1,j,nin2)
2905 628741662 : s=zin(2,j,nin2)
2906 628741662 : r2=r*cr2 - s*ci2
2907 628741662 : s2=r*ci2 + s*cr2
2908 628741662 : r=zin(1,j,nin3)
2909 628741662 : s=zin(2,j,nin3)
2910 628741662 : r3=r*cr3 - s*ci3
2911 628741662 : s3=r*ci3 + s*cr3
2912 628741662 : r=zin(1,j,nin4)
2913 628741662 : s=zin(2,j,nin4)
2914 628741662 : r4=r*cr4 - s*ci4
2915 628741662 : s4=r*ci4 + s*cr4
2916 628741662 : r=r1 + r3
2917 628741662 : s=r2 + r4
2918 628741662 : zout(1,j,nout1) = r + s
2919 628741662 : zout(1,j,nout3) = r - s
2920 628741662 : r=r1 - r3
2921 628741662 : s=s2 - s4
2922 628741662 : zout(1,j,nout2) = r - s
2923 628741662 : zout(1,j,nout4) = r + s
2924 628741662 : r=s1 + s3
2925 628741662 : s=s2 + s4
2926 628741662 : zout(2,j,nout1) = r + s
2927 628741662 : zout(2,j,nout3) = r - s
2928 628741662 : r=s1 - s3
2929 628741662 : s=r2 - r4
2930 628741662 : zout(2,j,nout2) = r + s
2931 677011546 : zout(2,j,nout4) = r - s
2932 : enddo
2933 : enddo
2934 : end if
2935 16507258 : 4000 continue
2936 : else
2937 14951605 : ia=1
2938 14951605 : nin1=ia-after
2939 14951605 : nout1=ia-atn
2940 64139150 : do ib=1,before
2941 49187545 : nin1=nin1+after
2942 49187545 : nin2=nin1+atb
2943 49187545 : nin3=nin2+atb
2944 49187545 : nin4=nin3+atb
2945 49187545 : nout1=nout1+atn
2946 49187545 : nout2=nout1+after
2947 49187545 : nout3=nout2+after
2948 49187545 : nout4=nout3+after
2949 751673250 : do j=1,n1dfft
2950 687534100 : r1=zin(1,j,nin1)
2951 687534100 : s1=zin(2,j,nin1)
2952 687534100 : r2=zin(1,j,nin2)
2953 687534100 : s2=zin(2,j,nin2)
2954 687534100 : r3=zin(1,j,nin3)
2955 687534100 : s3=zin(2,j,nin3)
2956 687534100 : r4=zin(1,j,nin4)
2957 687534100 : s4=zin(2,j,nin4)
2958 687534100 : r=r1 + r3
2959 687534100 : s=r2 + r4
2960 687534100 : zout(1,j,nout1) = r + s
2961 687534100 : zout(1,j,nout3) = r - s
2962 687534100 : r=r1 - r3
2963 687534100 : s=s2 - s4
2964 687534100 : zout(1,j,nout2) = r + s
2965 687534100 : zout(1,j,nout4) = r - s
2966 687534100 : r=s1 + s3
2967 687534100 : s=s2 + s4
2968 687534100 : zout(2,j,nout1) = r + s
2969 687534100 : zout(2,j,nout3) = r - s
2970 687534100 : r=s1 - s3
2971 687534100 : s=r2 - r4
2972 687534100 : zout(2,j,nout2) = r - s
2973 736721645 : zout(2,j,nout4) = r + s
2974 : enddo
2975 : enddo
2976 55030785 : do 4100,ia=2,after
2977 40079180 : ias=ia-1
2978 40079180 : if (2*ias.eq.after) then
2979 4244980 : nin1=ia-after
2980 4244980 : nout1=ia-atn
2981 10782712 : do ib=1,before
2982 6537732 : nin1=nin1+after
2983 6537732 : nin2=nin1+atb
2984 6537732 : nin3=nin2+atb
2985 6537732 : nin4=nin3+atb
2986 6537732 : nout1=nout1+atn
2987 6537732 : nout2=nout1+after
2988 6537732 : nout3=nout2+after
2989 6537732 : nout4=nout3+after
2990 100812550 : do j=1,n1dfft
2991 90029838 : r1=zin(1,j,nin1)
2992 90029838 : s1=zin(2,j,nin1)
2993 90029838 : r=zin(1,j,nin2)
2994 90029838 : s=zin(2,j,nin2)
2995 90029838 : r2=(r + s)*rt2i
2996 90029838 : s2=(s - r)*rt2i
2997 90029838 : r3=zin(2,j,nin3)
2998 90029838 : s3=zin(1,j,nin3)
2999 90029838 : r=zin(1,j,nin4)
3000 90029838 : s=zin(2,j,nin4)
3001 90029838 : r4=(s - r)*rt2i
3002 90029838 : s4=(r + s)*rt2i
3003 90029838 : r=r1 + r3
3004 90029838 : s=r2 + r4
3005 90029838 : zout(1,j,nout1) = r + s
3006 90029838 : zout(1,j,nout3) = r - s
3007 90029838 : r=r1 - r3
3008 90029838 : s=s2 + s4
3009 90029838 : zout(1,j,nout2) = r + s
3010 90029838 : zout(1,j,nout4) = r - s
3011 90029838 : r=s1 - s3
3012 90029838 : s=s2 - s4
3013 90029838 : zout(2,j,nout1) = r + s
3014 90029838 : zout(2,j,nout3) = r - s
3015 90029838 : r=s1 + s3
3016 90029838 : s=r2 - r4
3017 90029838 : zout(2,j,nout2) = r - s
3018 96567570 : zout(2,j,nout4) = r + s
3019 : enddo
3020 : enddo
3021 : else
3022 35834200 : itt=ias*before
3023 35834200 : itrig=itt+1
3024 35834200 : cr2=trig(1,itrig)
3025 35834200 : ci2=trig(2,itrig)
3026 35834200 : itrig=itrig+itt
3027 35834200 : cr3=trig(1,itrig)
3028 35834200 : ci3=trig(2,itrig)
3029 35834200 : itrig=itrig+itt
3030 35834200 : cr4=trig(1,itrig)
3031 35834200 : ci4=trig(2,itrig)
3032 35834200 : nin1=ia-after
3033 35834200 : nout1=ia-atn
3034 76940352 : do ib=1,before
3035 41106152 : nin1=nin1+after
3036 41106152 : nin2=nin1+atb
3037 41106152 : nin3=nin2+atb
3038 41106152 : nin4=nin3+atb
3039 41106152 : nout1=nout1+atn
3040 41106152 : nout2=nout1+after
3041 41106152 : nout3=nout2+after
3042 41106152 : nout4=nout3+after
3043 617734268 : do j=1,n1dfft
3044 540793916 : r1=zin(1,j,nin1)
3045 540793916 : s1=zin(2,j,nin1)
3046 540793916 : r=zin(1,j,nin2)
3047 540793916 : s=zin(2,j,nin2)
3048 540793916 : r2=r*cr2 - s*ci2
3049 540793916 : s2=r*ci2 + s*cr2
3050 540793916 : r=zin(1,j,nin3)
3051 540793916 : s=zin(2,j,nin3)
3052 540793916 : r3=r*cr3 - s*ci3
3053 540793916 : s3=r*ci3 + s*cr3
3054 540793916 : r=zin(1,j,nin4)
3055 540793916 : s=zin(2,j,nin4)
3056 540793916 : r4=r*cr4 - s*ci4
3057 540793916 : s4=r*ci4 + s*cr4
3058 540793916 : r=r1 + r3
3059 540793916 : s=r2 + r4
3060 540793916 : zout(1,j,nout1) = r + s
3061 540793916 : zout(1,j,nout3) = r - s
3062 540793916 : r=r1 - r3
3063 540793916 : s=s2 - s4
3064 540793916 : zout(1,j,nout2) = r + s
3065 540793916 : zout(1,j,nout4) = r - s
3066 540793916 : r=s1 + s3
3067 540793916 : s=s2 + s4
3068 540793916 : zout(2,j,nout1) = r + s
3069 540793916 : zout(2,j,nout3) = r - s
3070 540793916 : r=s1 - s3
3071 540793916 : s=r2 - r4
3072 540793916 : zout(2,j,nout2) = r - s
3073 581900068 : zout(2,j,nout4) = r + s
3074 : enddo
3075 : enddo
3076 : end if
3077 14951605 : 4100 continue
3078 : end if
3079 : else if (now.eq.8) then
3080 14955278 : if (isign.eq.-1) then
3081 7207400 : ia=1
3082 7207400 : nin1=ia-after
3083 7207400 : nout1=ia-atn
3084 31827600 : do ib=1,before
3085 24620200 : nin1=nin1+after
3086 24620200 : nin2=nin1+atb
3087 24620200 : nin3=nin2+atb
3088 24620200 : nin4=nin3+atb
3089 24620200 : nin5=nin4+atb
3090 24620200 : nin6=nin5+atb
3091 24620200 : nin7=nin6+atb
3092 24620200 : nin8=nin7+atb
3093 24620200 : nout1=nout1+atn
3094 24620200 : nout2=nout1+after
3095 24620200 : nout3=nout2+after
3096 24620200 : nout4=nout3+after
3097 24620200 : nout5=nout4+after
3098 24620200 : nout6=nout5+after
3099 24620200 : nout7=nout6+after
3100 24620200 : nout8=nout7+after
3101 417878672 : do j=1,n1dfft
3102 386051072 : r1=zin(1,j,nin1)
3103 386051072 : s1=zin(2,j,nin1)
3104 386051072 : r2=zin(1,j,nin2)
3105 386051072 : s2=zin(2,j,nin2)
3106 386051072 : r3=zin(1,j,nin3)
3107 386051072 : s3=zin(2,j,nin3)
3108 386051072 : r4=zin(1,j,nin4)
3109 386051072 : s4=zin(2,j,nin4)
3110 386051072 : r5=zin(1,j,nin5)
3111 386051072 : s5=zin(2,j,nin5)
3112 386051072 : r6=zin(1,j,nin6)
3113 386051072 : s6=zin(2,j,nin6)
3114 386051072 : r7=zin(1,j,nin7)
3115 386051072 : s7=zin(2,j,nin7)
3116 386051072 : r8=zin(1,j,nin8)
3117 386051072 : s8=zin(2,j,nin8)
3118 386051072 : r=r1 + r5
3119 386051072 : s=r3 + r7
3120 386051072 : ap=r + s
3121 386051072 : am=r - s
3122 386051072 : r=r2 + r6
3123 386051072 : s=r4 + r8
3124 386051072 : bp=r + s
3125 386051072 : bm=r - s
3126 386051072 : r=s1 + s5
3127 386051072 : s=s3 + s7
3128 386051072 : cp=r + s
3129 386051072 : cm=r - s
3130 386051072 : r=s2 + s6
3131 386051072 : s=s4 + s8
3132 386051072 : dpp=r + s
3133 386051072 : dm=r - s
3134 386051072 : zout(1,j,nout1) = ap + bp
3135 386051072 : zout(2,j,nout1) = cp + dpp
3136 386051072 : zout(1,j,nout5) = ap - bp
3137 386051072 : zout(2,j,nout5) = cp - dpp
3138 386051072 : zout(1,j,nout3) = am + dm
3139 386051072 : zout(2,j,nout3) = cm - bm
3140 386051072 : zout(1,j,nout7) = am - dm
3141 386051072 : zout(2,j,nout7) = cm + bm
3142 386051072 : r=r1 - r5
3143 386051072 : s=s3 - s7
3144 386051072 : ap=r + s
3145 386051072 : am=r - s
3146 386051072 : r=s1 - s5
3147 386051072 : s=r3 - r7
3148 386051072 : bp=r + s
3149 386051072 : bm=r - s
3150 386051072 : r=s4 - s8
3151 386051072 : s=r2 - r6
3152 386051072 : cp=r + s
3153 386051072 : cm=r - s
3154 386051072 : r=s2 - s6
3155 386051072 : s=r4 - r8
3156 386051072 : dpp=r + s
3157 386051072 : dm=r - s
3158 386051072 : r = ( cp + dm)*rt2i
3159 386051072 : s = ( dm - cp)*rt2i
3160 386051072 : cp= ( cm + dpp)*rt2i
3161 386051072 : dpp = ( cm - dpp)*rt2i
3162 386051072 : zout(1,j,nout2) = ap + r
3163 386051072 : zout(2,j,nout2) = bm + s
3164 386051072 : zout(1,j,nout6) = ap - r
3165 386051072 : zout(2,j,nout6) = bm - s
3166 386051072 : zout(1,j,nout4) = am + cp
3167 386051072 : zout(2,j,nout4) = bp + dpp
3168 386051072 : zout(1,j,nout8) = am - cp
3169 410671272 : zout(2,j,nout8) = bp - dpp
3170 : enddo
3171 : enddo
3172 7207400 : do 8000,ia=2,after
3173 0 : ias=ia-1
3174 0 : itt=ias*before
3175 0 : itrig=itt+1
3176 0 : cr2=trig(1,itrig)
3177 0 : ci2=trig(2,itrig)
3178 0 : itrig=itrig+itt
3179 0 : cr3=trig(1,itrig)
3180 0 : ci3=trig(2,itrig)
3181 0 : itrig=itrig+itt
3182 0 : cr4=trig(1,itrig)
3183 0 : ci4=trig(2,itrig)
3184 0 : itrig=itrig+itt
3185 0 : cr5=trig(1,itrig)
3186 0 : ci5=trig(2,itrig)
3187 0 : itrig=itrig+itt
3188 0 : cr6=trig(1,itrig)
3189 0 : ci6=trig(2,itrig)
3190 0 : itrig=itrig+itt
3191 0 : cr7=trig(1,itrig)
3192 0 : ci7=trig(2,itrig)
3193 0 : itrig=itrig+itt
3194 0 : cr8=trig(1,itrig)
3195 0 : ci8=trig(2,itrig)
3196 0 : nin1=ia-after
3197 0 : nout1=ia-atn
3198 0 : do ib=1,before
3199 0 : nin1=nin1+after
3200 0 : nin2=nin1+atb
3201 0 : nin3=nin2+atb
3202 0 : nin4=nin3+atb
3203 0 : nin5=nin4+atb
3204 0 : nin6=nin5+atb
3205 0 : nin7=nin6+atb
3206 0 : nin8=nin7+atb
3207 0 : nout1=nout1+atn
3208 0 : nout2=nout1+after
3209 0 : nout3=nout2+after
3210 0 : nout4=nout3+after
3211 0 : nout5=nout4+after
3212 0 : nout6=nout5+after
3213 0 : nout7=nout6+after
3214 0 : nout8=nout7+after
3215 0 : do j=1,n1dfft
3216 0 : r1=zin(1,j,nin1)
3217 0 : s1=zin(2,j,nin1)
3218 0 : r=zin(1,j,nin2)
3219 0 : s=zin(2,j,nin2)
3220 0 : r2=r*cr2 - s*ci2
3221 0 : s2=r*ci2 + s*cr2
3222 0 : r=zin(1,j,nin3)
3223 0 : s=zin(2,j,nin3)
3224 0 : r3=r*cr3 - s*ci3
3225 0 : s3=r*ci3 + s*cr3
3226 0 : r=zin(1,j,nin4)
3227 0 : s=zin(2,j,nin4)
3228 0 : r4=r*cr4 - s*ci4
3229 0 : s4=r*ci4 + s*cr4
3230 0 : r=zin(1,j,nin5)
3231 0 : s=zin(2,j,nin5)
3232 0 : r5=r*cr5 - s*ci5
3233 0 : s5=r*ci5 + s*cr5
3234 0 : r=zin(1,j,nin6)
3235 0 : s=zin(2,j,nin6)
3236 0 : r6=r*cr6 - s*ci6
3237 0 : s6=r*ci6 + s*cr6
3238 0 : r=zin(1,j,nin7)
3239 0 : s=zin(2,j,nin7)
3240 0 : r7=r*cr7 - s*ci7
3241 0 : s7=r*ci7 + s*cr7
3242 0 : r=zin(1,j,nin8)
3243 0 : s=zin(2,j,nin8)
3244 0 : r8=r*cr8 - s*ci8
3245 0 : s8=r*ci8 + s*cr8
3246 0 : r=r1 + r5
3247 0 : s=r3 + r7
3248 0 : ap=r + s
3249 0 : am=r - s
3250 0 : r=r2 + r6
3251 0 : s=r4 + r8
3252 0 : bp=r + s
3253 0 : bm=r - s
3254 0 : r=s1 + s5
3255 0 : s=s3 + s7
3256 0 : cp=r + s
3257 0 : cm=r - s
3258 0 : r=s2 + s6
3259 0 : s=s4 + s8
3260 0 : dpp=r + s
3261 0 : dm=r - s
3262 0 : zout(1,j,nout1) = ap + bp
3263 0 : zout(2,j,nout1) = cp + dpp
3264 0 : zout(1,j,nout5) = ap - bp
3265 0 : zout(2,j,nout5) = cp - dpp
3266 0 : zout(1,j,nout3) = am + dm
3267 0 : zout(2,j,nout3) = cm - bm
3268 0 : zout(1,j,nout7) = am - dm
3269 0 : zout(2,j,nout7) = cm + bm
3270 0 : r=r1 - r5
3271 0 : s=s3 - s7
3272 0 : ap=r + s
3273 0 : am=r - s
3274 0 : r=s1 - s5
3275 0 : s=r3 - r7
3276 0 : bp=r + s
3277 0 : bm=r - s
3278 0 : r=s4 - s8
3279 0 : s=r2 - r6
3280 0 : cp=r + s
3281 0 : cm=r - s
3282 0 : r=s2 - s6
3283 0 : s=r4 - r8
3284 0 : dpp=r + s
3285 0 : dm=r - s
3286 0 : r = ( cp + dm)*rt2i
3287 0 : s = ( dm - cp)*rt2i
3288 0 : cp= ( cm + dpp)*rt2i
3289 0 : dpp = ( cm - dpp)*rt2i
3290 0 : zout(1,j,nout2) = ap + r
3291 0 : zout(2,j,nout2) = bm + s
3292 0 : zout(1,j,nout6) = ap - r
3293 0 : zout(2,j,nout6) = bm - s
3294 0 : zout(1,j,nout4) = am + cp
3295 0 : zout(2,j,nout4) = bp + dpp
3296 0 : zout(1,j,nout8) = am - cp
3297 0 : zout(2,j,nout8) = bp - dpp
3298 : enddo
3299 : enddo
3300 7207400 : 8000 continue
3301 :
3302 : else
3303 7747878 : ia=1
3304 7747878 : nin1=ia-after
3305 7747878 : nout1=ia-atn
3306 35093008 : do ib=1,before
3307 27345130 : nin1=nin1+after
3308 27345130 : nin2=nin1+atb
3309 27345130 : nin3=nin2+atb
3310 27345130 : nin4=nin3+atb
3311 27345130 : nin5=nin4+atb
3312 27345130 : nin6=nin5+atb
3313 27345130 : nin7=nin6+atb
3314 27345130 : nin8=nin7+atb
3315 27345130 : nout1=nout1+atn
3316 27345130 : nout2=nout1+after
3317 27345130 : nout3=nout2+after
3318 27345130 : nout4=nout3+after
3319 27345130 : nout5=nout4+after
3320 27345130 : nout6=nout5+after
3321 27345130 : nout7=nout6+after
3322 27345130 : nout8=nout7+after
3323 457251574 : do j=1,n1dfft
3324 422158566 : r1=zin(1,j,nin1)
3325 422158566 : s1=zin(2,j,nin1)
3326 422158566 : r2=zin(1,j,nin2)
3327 422158566 : s2=zin(2,j,nin2)
3328 422158566 : r3=zin(1,j,nin3)
3329 422158566 : s3=zin(2,j,nin3)
3330 422158566 : r4=zin(1,j,nin4)
3331 422158566 : s4=zin(2,j,nin4)
3332 422158566 : r5=zin(1,j,nin5)
3333 422158566 : s5=zin(2,j,nin5)
3334 422158566 : r6=zin(1,j,nin6)
3335 422158566 : s6=zin(2,j,nin6)
3336 422158566 : r7=zin(1,j,nin7)
3337 422158566 : s7=zin(2,j,nin7)
3338 422158566 : r8=zin(1,j,nin8)
3339 422158566 : s8=zin(2,j,nin8)
3340 422158566 : r=r1 + r5
3341 422158566 : s=r3 + r7
3342 422158566 : ap=r + s
3343 422158566 : am=r - s
3344 422158566 : r=r2 + r6
3345 422158566 : s=r4 + r8
3346 422158566 : bp=r + s
3347 422158566 : bm=r - s
3348 422158566 : r=s1 + s5
3349 422158566 : s=s3 + s7
3350 422158566 : cp=r + s
3351 422158566 : cm=r - s
3352 422158566 : r=s2 + s6
3353 422158566 : s=s4 + s8
3354 422158566 : dpp=r + s
3355 422158566 : dm=r - s
3356 422158566 : zout(1,j,nout1) = ap + bp
3357 422158566 : zout(2,j,nout1) = cp + dpp
3358 422158566 : zout(1,j,nout5) = ap - bp
3359 422158566 : zout(2,j,nout5) = cp - dpp
3360 422158566 : zout(1,j,nout3) = am - dm
3361 422158566 : zout(2,j,nout3) = cm + bm
3362 422158566 : zout(1,j,nout7) = am + dm
3363 422158566 : zout(2,j,nout7) = cm - bm
3364 422158566 : r= r1 - r5
3365 422158566 : s=-s3 + s7
3366 422158566 : ap=r + s
3367 422158566 : am=r - s
3368 422158566 : r=s1 - s5
3369 422158566 : s=r7 - r3
3370 422158566 : bp=r + s
3371 422158566 : bm=r - s
3372 422158566 : r=-s4 + s8
3373 422158566 : s= r2 - r6
3374 422158566 : cp=r + s
3375 422158566 : cm=r - s
3376 422158566 : r=-s2 + s6
3377 422158566 : s= r4 - r8
3378 422158566 : dpp=r + s
3379 422158566 : dm=r - s
3380 422158566 : r = ( cp + dm)*rt2i
3381 422158566 : s = ( cp - dm)*rt2i
3382 422158566 : cp= ( cm + dpp)*rt2i
3383 422158566 : dpp= ( dpp - cm)*rt2i
3384 422158566 : zout(1,j,nout2) = ap + r
3385 422158566 : zout(2,j,nout2) = bm + s
3386 422158566 : zout(1,j,nout6) = ap - r
3387 422158566 : zout(2,j,nout6) = bm - s
3388 422158566 : zout(1,j,nout4) = am + cp
3389 422158566 : zout(2,j,nout4) = bp + dpp
3390 422158566 : zout(1,j,nout8) = am - cp
3391 449503696 : zout(2,j,nout8) = bp - dpp
3392 : enddo
3393 : enddo
3394 :
3395 7747878 : do 8001,ia=2,after
3396 0 : ias=ia-1
3397 0 : itt=ias*before
3398 0 : itrig=itt+1
3399 0 : cr2=trig(1,itrig)
3400 0 : ci2=trig(2,itrig)
3401 0 : itrig=itrig+itt
3402 0 : cr3=trig(1,itrig)
3403 0 : ci3=trig(2,itrig)
3404 0 : itrig=itrig+itt
3405 0 : cr4=trig(1,itrig)
3406 0 : ci4=trig(2,itrig)
3407 0 : itrig=itrig+itt
3408 0 : cr5=trig(1,itrig)
3409 0 : ci5=trig(2,itrig)
3410 0 : itrig=itrig+itt
3411 0 : cr6=trig(1,itrig)
3412 0 : ci6=trig(2,itrig)
3413 0 : itrig=itrig+itt
3414 0 : cr7=trig(1,itrig)
3415 0 : ci7=trig(2,itrig)
3416 0 : itrig=itrig+itt
3417 0 : cr8=trig(1,itrig)
3418 0 : ci8=trig(2,itrig)
3419 0 : nin1=ia-after
3420 0 : nout1=ia-atn
3421 0 : do ib=1,before
3422 0 : nin1=nin1+after
3423 0 : nin2=nin1+atb
3424 0 : nin3=nin2+atb
3425 0 : nin4=nin3+atb
3426 0 : nin5=nin4+atb
3427 0 : nin6=nin5+atb
3428 0 : nin7=nin6+atb
3429 0 : nin8=nin7+atb
3430 0 : nout1=nout1+atn
3431 0 : nout2=nout1+after
3432 0 : nout3=nout2+after
3433 0 : nout4=nout3+after
3434 0 : nout5=nout4+after
3435 0 : nout6=nout5+after
3436 0 : nout7=nout6+after
3437 0 : nout8=nout7+after
3438 0 : do j=1,n1dfft
3439 0 : r1=zin(1,j,nin1)
3440 0 : s1=zin(2,j,nin1)
3441 0 : r=zin(1,j,nin2)
3442 0 : s=zin(2,j,nin2)
3443 0 : r2=r*cr2 - s*ci2
3444 0 : s2=r*ci2 + s*cr2
3445 0 : r=zin(1,j,nin3)
3446 0 : s=zin(2,j,nin3)
3447 0 : r3=r*cr3 - s*ci3
3448 0 : s3=r*ci3 + s*cr3
3449 0 : r=zin(1,j,nin4)
3450 0 : s=zin(2,j,nin4)
3451 0 : r4=r*cr4 - s*ci4
3452 0 : s4=r*ci4 + s*cr4
3453 0 : r=zin(1,j,nin5)
3454 0 : s=zin(2,j,nin5)
3455 0 : r5=r*cr5 - s*ci5
3456 0 : s5=r*ci5 + s*cr5
3457 0 : r=zin(1,j,nin6)
3458 0 : s=zin(2,j,nin6)
3459 0 : r6=r*cr6 - s*ci6
3460 0 : s6=r*ci6 + s*cr6
3461 0 : r=zin(1,j,nin7)
3462 0 : s=zin(2,j,nin7)
3463 0 : r7=r*cr7 - s*ci7
3464 0 : s7=r*ci7 + s*cr7
3465 0 : r=zin(1,j,nin8)
3466 0 : s=zin(2,j,nin8)
3467 0 : r8=r*cr8 - s*ci8
3468 0 : s8=r*ci8 + s*cr8
3469 0 : r=r1 + r5
3470 0 : s=r3 + r7
3471 0 : ap=r + s
3472 0 : am=r - s
3473 0 : r=r2 + r6
3474 0 : s=r4 + r8
3475 0 : bp=r + s
3476 0 : bm=r - s
3477 0 : r=s1 + s5
3478 0 : s=s3 + s7
3479 0 : cp=r + s
3480 0 : cm=r - s
3481 0 : r=s2 + s6
3482 0 : s=s4 + s8
3483 0 : dpp=r + s
3484 0 : dm=r - s
3485 0 : zout(1,j,nout1) = ap + bp
3486 0 : zout(2,j,nout1) = cp + dpp
3487 0 : zout(1,j,nout5) = ap - bp
3488 0 : zout(2,j,nout5) = cp - dpp
3489 0 : zout(1,j,nout3) = am - dm
3490 0 : zout(2,j,nout3) = cm + bm
3491 0 : zout(1,j,nout7) = am + dm
3492 0 : zout(2,j,nout7) = cm - bm
3493 0 : r= r1 - r5
3494 0 : s=-s3 + s7
3495 0 : ap=r + s
3496 0 : am=r - s
3497 0 : r=s1 - s5
3498 0 : s=r7 - r3
3499 0 : bp=r + s
3500 0 : bm=r - s
3501 0 : r=-s4 + s8
3502 0 : s= r2 - r6
3503 0 : cp=r + s
3504 0 : cm=r - s
3505 0 : r=-s2 + s6
3506 0 : s= r4 - r8
3507 0 : dpp=r + s
3508 0 : dm=r - s
3509 0 : r = ( cp + dm)*rt2i
3510 0 : s = ( cp - dm)*rt2i
3511 0 : cp= ( cm + dpp)*rt2i
3512 0 : dpp= ( dpp - cm)*rt2i
3513 0 : zout(1,j,nout2) = ap + r
3514 0 : zout(2,j,nout2) = bm + s
3515 0 : zout(1,j,nout6) = ap - r
3516 0 : zout(2,j,nout6) = bm - s
3517 0 : zout(1,j,nout4) = am + cp
3518 0 : zout(2,j,nout4) = bp + dpp
3519 0 : zout(1,j,nout8) = am - cp
3520 0 : zout(2,j,nout8) = bp - dpp
3521 : enddo
3522 : enddo
3523 7747878 : 8001 continue
3524 :
3525 : end if
3526 : else if (now.eq.3) then
3527 : ! .5d0*sqrt(3.d0)
3528 41736144 : bb=isign*0.8660254037844387d0
3529 41736144 : ia=1
3530 41736144 : nin1=ia-after
3531 41736144 : nout1=ia-atn
3532 93270558 : do ib=1,before
3533 51534414 : nin1=nin1+after
3534 51534414 : nin2=nin1+atb
3535 51534414 : nin3=nin2+atb
3536 51534414 : nout1=nout1+atn
3537 51534414 : nout2=nout1+after
3538 51534414 : nout3=nout2+after
3539 846438502 : do j=1,n1dfft
3540 753167944 : r1=zin(1,j,nin1)
3541 753167944 : s1=zin(2,j,nin1)
3542 753167944 : r2=zin(1,j,nin2)
3543 753167944 : s2=zin(2,j,nin2)
3544 753167944 : r3=zin(1,j,nin3)
3545 753167944 : s3=zin(2,j,nin3)
3546 753167944 : r=r2 + r3
3547 753167944 : s=s2 + s3
3548 753167944 : zout(1,j,nout1) = r + r1
3549 753167944 : zout(2,j,nout1) = s + s1
3550 753167944 : r1=r1 - .5d0*r
3551 753167944 : s1=s1 - .5d0*s
3552 753167944 : r2=bb*(r2-r3)
3553 753167944 : s2=bb*(s2-s3)
3554 753167944 : zout(1,j,nout2) = r1 - s2
3555 753167944 : zout(2,j,nout2) = s1 + r2
3556 753167944 : zout(1,j,nout3) = r1 + s2
3557 804702358 : zout(2,j,nout3) = s1 - r2
3558 : enddo
3559 : enddo
3560 318630604 : do 3000,ia=2,after
3561 276894460 : ias=ia-1
3562 276894460 : if (4*ias.eq.3*after) then
3563 23549312 : if (isign.eq.1) then
3564 11848850 : nin1=ia-after
3565 11848850 : nout1=ia-atn
3566 25484104 : do ib=1,before
3567 13635254 : nin1=nin1+after
3568 13635254 : nin2=nin1+atb
3569 13635254 : nin3=nin2+atb
3570 13635254 : nout1=nout1+atn
3571 13635254 : nout2=nout1+after
3572 13635254 : nout3=nout2+after
3573 237245078 : do j=1,n1dfft
3574 211760974 : r1=zin(1,j,nin1)
3575 211760974 : s1=zin(2,j,nin1)
3576 211760974 : r2=zin(2,j,nin2)
3577 211760974 : s2=zin(1,j,nin2)
3578 211760974 : r3=zin(1,j,nin3)
3579 211760974 : s3=zin(2,j,nin3)
3580 211760974 : r=r3 + r2
3581 211760974 : s=s2 - s3
3582 211760974 : zout(1,j,nout1) = r1 - r
3583 211760974 : zout(2,j,nout1) = s + s1
3584 211760974 : r1=r1 + .5d0*r
3585 211760974 : s1=s1 - .5d0*s
3586 211760974 : r2=bb*(r2-r3)
3587 211760974 : s2=bb*(s2+s3)
3588 211760974 : zout(1,j,nout2) = r1 - s2
3589 211760974 : zout(2,j,nout2) = s1 - r2
3590 211760974 : zout(1,j,nout3) = r1 + s2
3591 225396228 : zout(2,j,nout3) = s1 + r2
3592 : enddo
3593 : enddo
3594 : else
3595 11700462 : nin1=ia-after
3596 11700462 : nout1=ia-atn
3597 25645650 : do ib=1,before
3598 13945188 : nin1=nin1+after
3599 13945188 : nin2=nin1+atb
3600 13945188 : nin3=nin2+atb
3601 13945188 : nout1=nout1+atn
3602 13945188 : nout2=nout1+after
3603 13945188 : nout3=nout2+after
3604 246406450 : do j=1,n1dfft
3605 220760800 : r1=zin(1,j,nin1)
3606 220760800 : s1=zin(2,j,nin1)
3607 220760800 : r2=zin(2,j,nin2)
3608 220760800 : s2=zin(1,j,nin2)
3609 220760800 : r3=zin(1,j,nin3)
3610 220760800 : s3=zin(2,j,nin3)
3611 220760800 : r=r2 - r3
3612 220760800 : s=s2 + s3
3613 220760800 : zout(1,j,nout1) = r + r1
3614 220760800 : zout(2,j,nout1) = s1 - s
3615 220760800 : r1=r1 - .5d0*r
3616 220760800 : s1=s1 + .5d0*s
3617 220760800 : r2=bb*(r2+r3)
3618 220760800 : s2=bb*(s2-s3)
3619 220760800 : zout(1,j,nout2) = r1 + s2
3620 220760800 : zout(2,j,nout2) = s1 + r2
3621 220760800 : zout(1,j,nout3) = r1 - s2
3622 234705988 : zout(2,j,nout3) = s1 - r2
3623 : enddo
3624 : enddo
3625 : end if
3626 253345148 : else if (8*ias.eq.3*after) then
3627 16373586 : if (isign.eq.1) then
3628 8553538 : nin1=ia-after
3629 8553538 : nout1=ia-atn
3630 17462420 : do ib=1,before
3631 8908882 : nin1=nin1+after
3632 8908882 : nin2=nin1+atb
3633 8908882 : nin3=nin2+atb
3634 8908882 : nout1=nout1+atn
3635 8908882 : nout2=nout1+after
3636 8908882 : nout3=nout2+after
3637 155425926 : do j=1,n1dfft
3638 137963506 : r1=zin(1,j,nin1)
3639 137963506 : s1=zin(2,j,nin1)
3640 137963506 : r=zin(1,j,nin2)
3641 137963506 : s=zin(2,j,nin2)
3642 137963506 : r2=(r - s)*rt2i
3643 137963506 : s2=(r + s)*rt2i
3644 137963506 : r3=zin(2,j,nin3)
3645 137963506 : s3=zin(1,j,nin3)
3646 137963506 : r=r2 - r3
3647 137963506 : s=s2 + s3
3648 137963506 : zout(1,j,nout1) = r + r1
3649 137963506 : zout(2,j,nout1) = s + s1
3650 137963506 : r1=r1 - .5d0*r
3651 137963506 : s1=s1 - .5d0*s
3652 137963506 : r2=bb*(r2+r3)
3653 137963506 : s2=bb*(s2-s3)
3654 137963506 : zout(1,j,nout2) = r1 - s2
3655 137963506 : zout(2,j,nout2) = s1 + r2
3656 137963506 : zout(1,j,nout3) = r1 + s2
3657 146872388 : zout(2,j,nout3) = s1 - r2
3658 : enddo
3659 : enddo
3660 : else
3661 7820048 : nin1=ia-after
3662 7820048 : nout1=ia-atn
3663 15842872 : do ib=1,before
3664 8022824 : nin1=nin1+after
3665 8022824 : nin2=nin1+atb
3666 8022824 : nin3=nin2+atb
3667 8022824 : nout1=nout1+atn
3668 8022824 : nout2=nout1+after
3669 8022824 : nout3=nout2+after
3670 141109328 : do j=1,n1dfft
3671 125266456 : r1=zin(1,j,nin1)
3672 125266456 : s1=zin(2,j,nin1)
3673 125266456 : r=zin(1,j,nin2)
3674 125266456 : s=zin(2,j,nin2)
3675 125266456 : r2=(r + s)*rt2i
3676 125266456 : s2=(s - r)*rt2i
3677 125266456 : r3=zin(2,j,nin3)
3678 125266456 : s3=zin(1,j,nin3)
3679 125266456 : r=r2 + r3
3680 125266456 : s=s2 - s3
3681 125266456 : zout(1,j,nout1) = r + r1
3682 125266456 : zout(2,j,nout1) = s + s1
3683 125266456 : r1=r1 - .5d0*r
3684 125266456 : s1=s1 - .5d0*s
3685 125266456 : r2=bb*(r2-r3)
3686 125266456 : s2=bb*(s2+s3)
3687 125266456 : zout(1,j,nout2) = r1 - s2
3688 125266456 : zout(2,j,nout2) = s1 + r2
3689 125266456 : zout(1,j,nout3) = r1 + s2
3690 133289280 : zout(2,j,nout3) = s1 - r2
3691 : enddo
3692 : enddo
3693 : end if
3694 : else
3695 236971562 : itt=ias*before
3696 236971562 : itrig=itt+1
3697 236971562 : cr2=trig(1,itrig)
3698 236971562 : ci2=trig(2,itrig)
3699 236971562 : itrig=itrig+itt
3700 236971562 : cr3=trig(1,itrig)
3701 236971562 : ci3=trig(2,itrig)
3702 236971562 : nin1=ia-after
3703 236971562 : nout1=ia-atn
3704 485986600 : do ib=1,before
3705 249015038 : nin1=nin1+after
3706 249015038 : nin2=nin1+atb
3707 249015038 : nin3=nin2+atb
3708 249015038 : nout1=nout1+atn
3709 249015038 : nout2=nout1+after
3710 249015038 : nout3=nout2+after
3711 4029364594 : do j=1,n1dfft
3712 3543377994 : r1=zin(1,j,nin1)
3713 3543377994 : s1=zin(2,j,nin1)
3714 3543377994 : r=zin(1,j,nin2)
3715 3543377994 : s=zin(2,j,nin2)
3716 3543377994 : r2=r*cr2 - s*ci2
3717 3543377994 : s2=r*ci2 + s*cr2
3718 3543377994 : r=zin(1,j,nin3)
3719 3543377994 : s=zin(2,j,nin3)
3720 3543377994 : r3=r*cr3 - s*ci3
3721 3543377994 : s3=r*ci3 + s*cr3
3722 3543377994 : r=r2 + r3
3723 3543377994 : s=s2 + s3
3724 3543377994 : zout(1,j,nout1) = r + r1
3725 3543377994 : zout(2,j,nout1) = s + s1
3726 3543377994 : r1=r1 - .5d0*r
3727 3543377994 : s1=s1 - .5d0*s
3728 3543377994 : r2=bb*(r2-r3)
3729 3543377994 : s2=bb*(s2-s3)
3730 3543377994 : zout(1,j,nout2) = r1 - s2
3731 3543377994 : zout(2,j,nout2) = s1 + r2
3732 3543377994 : zout(1,j,nout3) = r1 + s2
3733 3792393032 : zout(2,j,nout3) = s1 - r2
3734 : enddo
3735 : enddo
3736 : end if
3737 41736144 : 3000 continue
3738 : else if (now==5) then
3739 : ! cos(2.d0*pi/5.d0)
3740 13497133 : cos2=0.3090169943749474d0
3741 : ! cos(4.d0*pi/5.d0)
3742 13497133 : cos4=-0.8090169943749474d0
3743 : ! sin(2.d0*pi/5.d0)
3744 13497133 : sin2=isign*0.9510565162951536d0
3745 : ! sin(4.d0*pi/5.d0)
3746 13497133 : sin4=isign*0.5877852522924731d0
3747 13497133 : ia=1
3748 13497133 : nin1=ia-after
3749 13497133 : nout1=ia-atn
3750 81676684 : do ib=1,before
3751 68179551 : nin1=nin1+after
3752 68179551 : nin2=nin1+atb
3753 68179551 : nin3=nin2+atb
3754 68179551 : nin4=nin3+atb
3755 68179551 : nin5=nin4+atb
3756 68179551 : nout1=nout1+atn
3757 68179551 : nout2=nout1+after
3758 68179551 : nout3=nout2+after
3759 68179551 : nout4=nout3+after
3760 68179551 : nout5=nout4+after
3761 902509152 : do j=1,n1dfft
3762 820832468 : r1=zin(1,j,nin1)
3763 820832468 : s1=zin(2,j,nin1)
3764 820832468 : r2=zin(1,j,nin2)
3765 820832468 : s2=zin(2,j,nin2)
3766 820832468 : r3=zin(1,j,nin3)
3767 820832468 : s3=zin(2,j,nin3)
3768 820832468 : r4=zin(1,j,nin4)
3769 820832468 : s4=zin(2,j,nin4)
3770 820832468 : r5=zin(1,j,nin5)
3771 820832468 : s5=zin(2,j,nin5)
3772 820832468 : r25 = r2 + r5
3773 820832468 : r34 = r3 + r4
3774 820832468 : s25 = s2 - s5
3775 820832468 : s34 = s3 - s4
3776 820832468 : zout(1,j,nout1) = r1 + r25 + r34
3777 820832468 : r = r1 + cos2*r25 + cos4*r34
3778 820832468 : s = sin2*s25 + sin4*s34
3779 820832468 : zout(1,j,nout2) = r - s
3780 820832468 : zout(1,j,nout5) = r + s
3781 820832468 : r = r1 + cos4*r25 + cos2*r34
3782 820832468 : s = sin4*s25 - sin2*s34
3783 820832468 : zout(1,j,nout3) = r - s
3784 820832468 : zout(1,j,nout4) = r + s
3785 820832468 : r25 = r2 - r5
3786 820832468 : r34 = r3 - r4
3787 820832468 : s25 = s2 + s5
3788 820832468 : s34 = s3 + s4
3789 820832468 : zout(2,j,nout1) = s1 + s25 + s34
3790 820832468 : r = s1 + cos2*s25 + cos4*s34
3791 820832468 : s = sin2*r25 + sin4*r34
3792 820832468 : zout(2,j,nout2) = r + s
3793 820832468 : zout(2,j,nout5) = r - s
3794 820832468 : r = s1 + cos4*s25 + cos2*s34
3795 820832468 : s = sin4*r25 - sin2*r34
3796 820832468 : zout(2,j,nout3) = r + s
3797 889012019 : zout(2,j,nout4) = r - s
3798 : enddo
3799 : enddo
3800 25223105 : do 5000,ia=2,after
3801 11725972 : ias=ia-1
3802 11725972 : if (8*ias.eq.5*after) then
3803 360192 : if (isign.eq.1) then
3804 205824 : nin1=ia-after
3805 205824 : nout1=ia-atn
3806 823296 : do ib=1,before
3807 617472 : nin1=nin1+after
3808 617472 : nin2=nin1+atb
3809 617472 : nin3=nin2+atb
3810 617472 : nin4=nin3+atb
3811 617472 : nin5=nin4+atb
3812 617472 : nout1=nout1+atn
3813 617472 : nout2=nout1+after
3814 617472 : nout3=nout2+after
3815 617472 : nout4=nout3+after
3816 617472 : nout5=nout4+after
3817 5525184 : do j=1,n1dfft
3818 4701888 : r1=zin(1,j,nin1)
3819 4701888 : s1=zin(2,j,nin1)
3820 4701888 : r=zin(1,j,nin2)
3821 4701888 : s=zin(2,j,nin2)
3822 4701888 : r2=(r - s)*rt2i
3823 4701888 : s2=(r + s)*rt2i
3824 4701888 : r3=zin(2,j,nin3)
3825 4701888 : s3=zin(1,j,nin3)
3826 4701888 : r=zin(1,j,nin4)
3827 4701888 : s=zin(2,j,nin4)
3828 4701888 : r4=(r + s)*rt2i
3829 4701888 : s4=(r - s)*rt2i
3830 4701888 : r5=zin(1,j,nin5)
3831 4701888 : s5=zin(2,j,nin5)
3832 4701888 : r25 = r2 - r5
3833 4701888 : r34 = r3 + r4
3834 4701888 : s25 = s2 + s5
3835 4701888 : s34 = s3 - s4
3836 4701888 : zout(1,j,nout1) = r1 + r25 - r34
3837 4701888 : r = r1 + cos2*r25 - cos4*r34
3838 4701888 : s = sin2*s25 + sin4*s34
3839 4701888 : zout(1,j,nout2) = r - s
3840 4701888 : zout(1,j,nout5) = r + s
3841 4701888 : r = r1 + cos4*r25 - cos2*r34
3842 4701888 : s = sin4*s25 - sin2*s34
3843 4701888 : zout(1,j,nout3) = r - s
3844 4701888 : zout(1,j,nout4) = r + s
3845 4701888 : r25 = r2 + r5
3846 4701888 : r34 = r4 - r3
3847 4701888 : s25 = s2 - s5
3848 4701888 : s34 = s3 + s4
3849 4701888 : zout(2,j,nout1) = s1 + s25 + s34
3850 4701888 : r = s1 + cos2*s25 + cos4*s34
3851 4701888 : s = sin2*r25 + sin4*r34
3852 4701888 : zout(2,j,nout2) = r + s
3853 4701888 : zout(2,j,nout5) = r - s
3854 4701888 : r = s1 + cos4*s25 + cos2*s34
3855 4701888 : s = sin4*r25 - sin2*r34
3856 4701888 : zout(2,j,nout3) = r + s
3857 5319360 : zout(2,j,nout4) = r - s
3858 : enddo
3859 : enddo
3860 : else
3861 154368 : nin1=ia-after
3862 154368 : nout1=ia-atn
3863 617472 : do ib=1,before
3864 463104 : nin1=nin1+after
3865 463104 : nin2=nin1+atb
3866 463104 : nin3=nin2+atb
3867 463104 : nin4=nin3+atb
3868 463104 : nin5=nin4+atb
3869 463104 : nout1=nout1+atn
3870 463104 : nout2=nout1+after
3871 463104 : nout3=nout2+after
3872 463104 : nout4=nout3+after
3873 463104 : nout5=nout4+after
3874 4143888 : do j=1,n1dfft
3875 3526416 : r1=zin(1,j,nin1)
3876 3526416 : s1=zin(2,j,nin1)
3877 3526416 : r=zin(1,j,nin2)
3878 3526416 : s=zin(2,j,nin2)
3879 3526416 : r2=(r + s)*rt2i
3880 3526416 : s2=(s - r)*rt2i
3881 3526416 : r3=zin(2,j,nin3)
3882 3526416 : s3=zin(1,j,nin3)
3883 3526416 : r=zin(1,j,nin4)
3884 3526416 : s=zin(2,j,nin4)
3885 3526416 : r4=(s - r)*rt2i
3886 3526416 : s4=(r + s)*rt2i
3887 3526416 : r5=zin(1,j,nin5)
3888 3526416 : s5=zin(2,j,nin5)
3889 3526416 : r25 = r2 - r5
3890 3526416 : r34 = r3 + r4
3891 3526416 : s25 = s2 + s5
3892 3526416 : s34 = s4 - s3
3893 3526416 : zout(1,j,nout1) = r1 + r25 + r34
3894 3526416 : r = r1 + cos2*r25 + cos4*r34
3895 3526416 : s = sin2*s25 + sin4*s34
3896 3526416 : zout(1,j,nout2) = r - s
3897 3526416 : zout(1,j,nout5) = r + s
3898 3526416 : r = r1 + cos4*r25 + cos2*r34
3899 3526416 : s = sin4*s25 - sin2*s34
3900 3526416 : zout(1,j,nout3) = r - s
3901 3526416 : zout(1,j,nout4) = r + s
3902 3526416 : r25 = r2 + r5
3903 3526416 : r34 = r3 - r4
3904 3526416 : s25 = s2 - s5
3905 3526416 : s34 = s3 + s4
3906 3526416 : zout(2,j,nout1) = s1 + s25 - s34
3907 3526416 : r = s1 + cos2*s25 - cos4*s34
3908 3526416 : s = sin2*r25 + sin4*r34
3909 3526416 : zout(2,j,nout2) = r + s
3910 3526416 : zout(2,j,nout5) = r - s
3911 3526416 : r = s1 + cos4*s25 - cos2*s34
3912 3526416 : s = sin4*r25 - sin2*r34
3913 3526416 : zout(2,j,nout3) = r + s
3914 3989520 : zout(2,j,nout4) = r - s
3915 : enddo
3916 : enddo
3917 : end if
3918 : else
3919 11365780 : ias=ia-1
3920 11365780 : itt=ias*before
3921 11365780 : itrig=itt+1
3922 11365780 : cr2=trig(1,itrig)
3923 11365780 : ci2=trig(2,itrig)
3924 11365780 : itrig=itrig+itt
3925 11365780 : cr3=trig(1,itrig)
3926 11365780 : ci3=trig(2,itrig)
3927 11365780 : itrig=itrig+itt
3928 11365780 : cr4=trig(1,itrig)
3929 11365780 : ci4=trig(2,itrig)
3930 11365780 : itrig=itrig+itt
3931 11365780 : cr5=trig(1,itrig)
3932 11365780 : ci5=trig(2,itrig)
3933 11365780 : nin1=ia-after
3934 11365780 : nout1=ia-atn
3935 39970152 : do ib=1,before
3936 28604372 : nin1=nin1+after
3937 28604372 : nin2=nin1+atb
3938 28604372 : nin3=nin2+atb
3939 28604372 : nin4=nin3+atb
3940 28604372 : nin5=nin4+atb
3941 28604372 : nout1=nout1+atn
3942 28604372 : nout2=nout1+after
3943 28604372 : nout3=nout2+after
3944 28604372 : nout4=nout3+after
3945 28604372 : nout5=nout4+after
3946 360313019 : do j=1,n1dfft
3947 320342867 : r1=zin(1,j,nin1)
3948 320342867 : s1=zin(2,j,nin1)
3949 320342867 : r=zin(1,j,nin2)
3950 320342867 : s=zin(2,j,nin2)
3951 320342867 : r2=r*cr2 - s*ci2
3952 320342867 : s2=r*ci2 + s*cr2
3953 320342867 : r=zin(1,j,nin3)
3954 320342867 : s=zin(2,j,nin3)
3955 320342867 : r3=r*cr3 - s*ci3
3956 320342867 : s3=r*ci3 + s*cr3
3957 320342867 : r=zin(1,j,nin4)
3958 320342867 : s=zin(2,j,nin4)
3959 320342867 : r4=r*cr4 - s*ci4
3960 320342867 : s4=r*ci4 + s*cr4
3961 320342867 : r=zin(1,j,nin5)
3962 320342867 : s=zin(2,j,nin5)
3963 320342867 : r5=r*cr5 - s*ci5
3964 320342867 : s5=r*ci5 + s*cr5
3965 320342867 : r25 = r2 + r5
3966 320342867 : r34 = r3 + r4
3967 320342867 : s25 = s2 - s5
3968 320342867 : s34 = s3 - s4
3969 320342867 : zout(1,j,nout1) = r1 + r25 + r34
3970 320342867 : r = r1 + cos2*r25 + cos4*r34
3971 320342867 : s = sin2*s25 + sin4*s34
3972 320342867 : zout(1,j,nout2) = r - s
3973 320342867 : zout(1,j,nout5) = r + s
3974 320342867 : r = r1 + cos4*r25 + cos2*r34
3975 320342867 : s = sin4*s25 - sin2*s34
3976 320342867 : zout(1,j,nout3) = r - s
3977 320342867 : zout(1,j,nout4) = r + s
3978 320342867 : r25 = r2 - r5
3979 320342867 : r34 = r3 - r4
3980 320342867 : s25 = s2 + s5
3981 320342867 : s34 = s3 + s4
3982 320342867 : zout(2,j,nout1) = s1 + s25 + s34
3983 320342867 : r = s1 + cos2*s25 + cos4*s34
3984 320342867 : s = sin2*r25 + sin4*r34
3985 320342867 : zout(2,j,nout2) = r + s
3986 320342867 : zout(2,j,nout5) = r - s
3987 320342867 : r = s1 + cos4*s25 + cos2*s34
3988 320342867 : s = sin4*r25 - sin2*r34
3989 320342867 : zout(2,j,nout3) = r + s
3990 348947239 : zout(2,j,nout4) = r - s
3991 : enddo
3992 : enddo
3993 : end if
3994 13497133 : 5000 continue
3995 : else if (now.eq.6) then
3996 : ! .5d0*sqrt(3.d0)
3997 15843000 : bb=isign*0.8660254037844387d0
3998 :
3999 15843000 : ia=1
4000 15843000 : nin1=ia-after
4001 15843000 : nout1=ia-atn
4002 65168744 : do ib=1,before
4003 49325744 : nin1=nin1+after
4004 49325744 : nin2=nin1+atb
4005 49325744 : nin3=nin2+atb
4006 49325744 : nin4=nin3+atb
4007 49325744 : nin5=nin4+atb
4008 49325744 : nin6=nin5+atb
4009 49325744 : nout1=nout1+atn
4010 49325744 : nout2=nout1+after
4011 49325744 : nout3=nout2+after
4012 49325744 : nout4=nout3+after
4013 49325744 : nout5=nout4+after
4014 49325744 : nout6=nout5+after
4015 737408051 : do j=1,n1dfft
4016 672239307 : r2=zin(1,j,nin3)
4017 672239307 : s2=zin(2,j,nin3)
4018 672239307 : r3=zin(1,j,nin5)
4019 672239307 : s3=zin(2,j,nin5)
4020 672239307 : r=r2 + r3
4021 672239307 : s=s2 + s3
4022 672239307 : r1=zin(1,j,nin1)
4023 672239307 : s1=zin(2,j,nin1)
4024 672239307 : ur1 = r + r1
4025 672239307 : ui1 = s + s1
4026 672239307 : r1=r1 - .5d0*r
4027 672239307 : s1=s1 - .5d0*s
4028 672239307 : r=r2-r3
4029 672239307 : s=s2-s3
4030 672239307 : ur2 = r1 - s*bb
4031 672239307 : ui2 = s1 + r*bb
4032 672239307 : ur3 = r1 + s*bb
4033 672239307 : ui3 = s1 - r*bb
4034 :
4035 672239307 : r2=zin(1,j,nin6)
4036 672239307 : s2=zin(2,j,nin6)
4037 672239307 : r3=zin(1,j,nin2)
4038 672239307 : s3=zin(2,j,nin2)
4039 672239307 : r=r2 + r3
4040 672239307 : s=s2 + s3
4041 672239307 : r1=zin(1,j,nin4)
4042 672239307 : s1=zin(2,j,nin4)
4043 672239307 : vr1 = r + r1
4044 672239307 : vi1 = s + s1
4045 672239307 : r1=r1 - .5d0*r
4046 672239307 : s1=s1 - .5d0*s
4047 672239307 : r=r2-r3
4048 672239307 : s=s2-s3
4049 672239307 : vr2 = r1 - s*bb
4050 672239307 : vi2 = s1 + r*bb
4051 672239307 : vr3 = r1 + s*bb
4052 672239307 : vi3 = s1 - r*bb
4053 :
4054 672239307 : zout(1,j,nout1)=ur1+vr1
4055 672239307 : zout(2,j,nout1)=ui1+vi1
4056 672239307 : zout(1,j,nout5)=ur2+vr2
4057 672239307 : zout(2,j,nout5)=ui2+vi2
4058 672239307 : zout(1,j,nout3)=ur3+vr3
4059 672239307 : zout(2,j,nout3)=ui3+vi3
4060 672239307 : zout(1,j,nout4)=ur1-vr1
4061 672239307 : zout(2,j,nout4)=ui1-vi1
4062 672239307 : zout(1,j,nout2)=ur2-vr2
4063 672239307 : zout(2,j,nout2)=ui2-vi2
4064 672239307 : zout(1,j,nout6)=ur3-vr3
4065 721565051 : zout(2,j,nout6)=ui3-vi3
4066 : enddo
4067 : enddo
4068 :
4069 : else
4070 0 : ABI_ERROR('error fftstp')
4071 : end if
4072 :
4073 117986786 : end subroutine fftstp
4074 : !!***
4075 :
4076 : end module m_sg2002
4077 : !!***
|
