Branch data 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 : 71562 : 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 : 71562 : real(dp), allocatable :: zt(:,:,:) ! work arrays for transpositions
131 : 71562 : real(dp), allocatable :: zmpi1(:,:,:,:),zmpi2(:,:,:,:) ! work arrays for MPI
132 : 71562 : real(dp), allocatable :: zw(:,:,:) ! cache work array
133 : : real(dp) :: tsec(2)
134 : : ! FFT work arrays
135 : 71562 : real(dp), allocatable, dimension(:,:) :: trig1,trig2,trig3
136 : 71562 : integer, allocatable, dimension(:) :: after1,now1,before1,after2,now2,before2,after3,now3,before3
137 : :
138 : : ! *************************************************************************
139 : :
140 : 71562 : 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 : 71562 : ncache=4*max(n1,n2,n3,1024)
144 : :
145 [ - + ]: 71562 : 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 [ + - + - : 71562 : 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 : 71562 : n1eff=n1; n2eff=n2; n1zt=n1
160 [ + + ]: 71562 : if (cplex==1) then
161 : 92 : n1eff=(n1+1)/2 ; n2eff=n2/2+1 ; n1zt=2*(n1/2+1)
162 : : end if
163 : :
164 : 71562 : lzt=n2eff
165 [ + + ]: 71562 : if (mod(n2eff,2) == 0) lzt=lzt+1
166 [ + + ]: 71562 : if (mod(n2eff,4) == 0) lzt=lzt+1
167 : :
168 : : ! maximal number of big box 3rd dim slices for all procs
169 : 71562 : nnd3=nd3proc*nproc_fft
170 : :
171 [ - + + - : 214686 : ABI_MALLOC(trig1,(2,n1))
+ - + - -
+ - + ]
172 [ - + ]: 71562 : ABI_MALLOC(after1,(mdata))
173 [ - + ]: 71562 : ABI_MALLOC(now1,(mdata))
174 [ - + ]: 71562 : ABI_MALLOC(before1,(mdata))
175 [ - + + - : 214686 : ABI_MALLOC(trig2,(2,n2))
+ - + - -
+ - + ]
176 [ - + ]: 71562 : ABI_MALLOC(after2,(mdata))
177 [ - + ]: 71562 : ABI_MALLOC(now2,(mdata))
178 [ - + ]: 71562 : ABI_MALLOC(before2,(mdata))
179 [ - + + - : 214686 : ABI_MALLOC(trig3,(2,n3))
+ - + - -
+ - + ]
180 [ - + ]: 71562 : ABI_MALLOC(after3,(mdata))
181 [ - + ]: 71562 : ABI_MALLOC(now3,(mdata))
182 [ - + ]: 71562 : ABI_MALLOC(before3,(mdata))
183 [ - + + - : 286248 : ABI_MALLOC(zw,(2,ncache/4,2))
+ - + - -
+ - + ]
184 [ - + + - : 286248 : ABI_MALLOC(zt,(2,lzt,n1zt))
- + + - +
- + - - +
- + ]
185 [ - + + - : 357810 : ABI_MALLOC(zmpi2,(2,n1,nd2proc,nnd3))
- + + - -
+ + - + -
+ - - + -
+ ]
186 [ + + ]: 71562 : if (nproc_fft>1) then
187 [ - + + - : 58752 : ABI_MALLOC(zmpi1,(2,n1,nd2proc,nnd3))
- + + - -
+ + - - +
- + ]
188 : : end if
189 : :
190 : 71562 : call ctrig(n3,trig3,after3,before3,now3,1,ic3)
191 : 71562 : call ctrig(n1,trig1,after1,before1,now1,1,ic1)
192 : 71562 : 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 [ + + ]: 143132 : do idat=1,ndat
204 : : ! transform along z axis
205 : : ! input: I1,I3,J2,(Jp2)
206 : 71570 : lot=ncache/(4*n3)
207 : :
208 [ + + ]: 1036989 : do j2=1,nd2proc
209 [ + + ]: 1036989 : if (me_fft*nd2proc+j2 <= n2eff) then
210 : :
211 [ - + - - : 2572904 : do i1=1,n1,lot
- + + + ]
212 : 1607490 : ma=i1
213 : 1607490 : mb=min(i1+(lot-1),n1)
214 : 1607490 : n1dfft=mb-ma+1
215 : :
216 : : ! input: G1,G3,G2,(Gp2)
217 : 1607490 : call fill(nd1,nd3,lot,n1dfft,n3,zf(1,i1,1,j2,idat),zw(1,1,1))
218 : :
219 : 1607490 : inzee=1
220 [ + + ]: 5665398 : do i=1,ic3
221 : : call fftstp(lot,n1dfft,n3,lot,n3,zw(1,1,inzee),zw(1,1,3-inzee), &
222 : 4057908 : & trig3,after3(i),now3(i),before3(i),1)
223 : 5665398 : inzee=3-inzee
224 : : end do
225 : :
226 : : ! input: G1,R3,G2,(Gp2)
227 : : ! output: G1,G2,R3,(Gp2)
228 : 2572904 : 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 [ + + ]: 71570 : if (nproc_fft>1) then
237 : 14688 : call timab(543,1,tsec)
238 : : call xmpi_alltoall(zmpi2,2*n1*nd2proc*nd3proc, &
239 : 14688 : & zmpi1,2*n1*nd2proc*nd3proc,comm_fft,ierr)
240 : 14688 : call timab(543,2,tsec)
241 : : end if
242 : :
243 [ + + ]: 1114315 : do j3=1,nd3proc
244 [ + - ]: 1042753 : if (me_fft*nd3proc+j3 <= n3) then
245 : 971183 : Jp2st=1
246 : 971183 : J2st=1
247 : :
248 : : ! transform along x axis
249 : 971183 : lot=ncache/(4*n1)
250 : :
251 [ - + - - : 2580674 : do j=1,n2eff,lot
- + + + ]
252 : 1609491 : ma=j
253 : 1609491 : mb=min(j+(lot-1),n2eff)
254 : 1609491 : n1dfft=mb-ma+1
255 : :
256 : : ! input: G1,G2,R3,Gp2,(Rp3)
257 : : ! output: G2,G1,R3,Jp2,(Rp3)
258 [ + + ]: 1609491 : 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 : 744016 : 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 : 1609491 : inzee=1
267 [ + + ]: 4041908 : do i=1,ic1-1
268 : : call fftstp(lot,n1dfft,n1,lot,n1,zw(1,1,inzee),zw(1,1,3-inzee), &
269 : 2432417 : & trig1,after1(i),now1(i),before1(i),1)
270 : 4041908 : inzee=3-inzee
271 : : end do
272 : :
273 : 1609491 : i=ic1
274 : : call fftstp(lot,n1dfft,n1,lzt,n1zt,zw(1,1,inzee),zt(1,j,1), &
275 : 2580674 : & trig1,after1(i),now1(i),before1(i),1)
276 : : end do
277 : :
278 : : ! transform along y axis
279 : 971183 : lot=ncache/(4*n2)
280 : :
281 [ - + - - : 2580674 : do j=1,n1eff,lot
- + + + ]
282 : 1609491 : ma=j
283 : 1609491 : mb=min(j+(lot-1),n1eff)
284 : 1609491 : n1dfft=mb-ma+1
285 : 1609491 : includelast=1
286 : :
287 [ + + ]: 1609491 : 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 [ + + ]: 1609491 : if (cplex==2) then
296 : 1603011 : 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 : 1609491 : inzee=1
302 [ + + ]: 4041668 : do i=1,ic2-1
303 : : call fftstp(lot,n1dfft,n2,lot,n2,zw(1,1,inzee),zw(1,1,3-inzee), &
304 : 2432177 : & trig2,after2(i),now2(i),before2(i),1)
305 : 4041668 : inzee=3-inzee
306 : : end do
307 : :
308 : 1609491 : i=ic2
309 : : call fftstp(lot,n1dfft,n2,nd1eff,nd2,zw(1,1,inzee),zr(1,j,1,j3,idat), &
310 : 2580674 : & 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 [ - + ]: 71562 : ABI_FREE(trig1)
320 [ - + ]: 71562 : ABI_FREE(after1)
321 [ - + ]: 71562 : ABI_FREE(now1)
322 [ - + ]: 71562 : ABI_FREE(before1)
323 [ - + ]: 71562 : ABI_FREE(trig2)
324 [ - + ]: 71562 : ABI_FREE(after2)
325 [ - + ]: 71562 : ABI_FREE(now2)
326 [ - + ]: 71562 : ABI_FREE(before2)
327 [ - + ]: 71562 : ABI_FREE(trig3)
328 [ - + ]: 71562 : ABI_FREE(after3)
329 [ - + ]: 71562 : ABI_FREE(now3)
330 [ - + ]: 71562 : ABI_FREE(before3)
331 [ - + ]: 71562 : ABI_FREE(zmpi2)
332 [ - + ]: 71562 : ABI_FREE(zw)
333 [ - + ]: 71562 : ABI_FREE(zt)
334 [ + + ]: 71562 : if (nproc_fft>1) then
335 [ - + ]: 14688 : ABI_FREE(zmpi1)
336 : : end if
337 : :
338 [ - + ]: 71562 : 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 : 71059 : 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 : 71059 : real(dp), allocatable :: zt(:,:,:) ! work arrays for transpositions
410 : 71059 : real(dp), allocatable :: zmpi1(:,:,:,:),zmpi2(:,:,:,:) ! work arrays for MPI
411 : 71059 : real(dp), allocatable :: zw(:,:,:) ! cache work array
412 : : real(dp) :: tsec(2)
413 : : ! FFT work arrays
414 : 71059 : real(dp), allocatable, dimension(:,:) :: trig1,trig2,trig3
415 : 71059 : integer, allocatable, dimension(:) :: after1,now1,before1,after2,now2,before2,after3,now3,before3
416 : :
417 : : ! *************************************************************************
418 : :
419 : 71059 : 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 : 71059 : ncache=4*max(n1,n2,n3,1024)
423 [ - + ]: 71059 : 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 [ + - + - : 71059 : 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 : 71059 : n1eff=n1; n2eff=n2; n1zt=n1
438 [ + + ]: 71059 : if (cplex==1) then
439 : 80 : n1eff=(n1+1)/2; n2eff=n2/2+1; n1zt=2*(n1/2+1)
440 : : end if
441 : :
442 : 71059 : lzt=n2eff
443 [ + + ]: 71059 : if (mod(n2eff,2) == 0) lzt=lzt+1
444 [ + + ]: 71059 : if (mod(n2eff,4) == 0) lzt=lzt+1
445 : :
446 : : ! maximal number of big box 3rd dim slices for all procs
447 : 71059 : nnd3=nd3proc*nproc_fft
448 : :
449 [ - + + - : 213177 : ABI_MALLOC(trig1,(2,n1))
+ - + - -
+ - + ]
450 [ - + ]: 71059 : ABI_MALLOC(after1,(mdata))
451 [ - + ]: 71059 : ABI_MALLOC(now1,(mdata))
452 [ - + ]: 71059 : ABI_MALLOC(before1,(mdata))
453 [ - + + - : 213177 : ABI_MALLOC(trig2,(2,n2))
+ - + - -
+ - + ]
454 [ - + ]: 71059 : ABI_MALLOC(after2,(mdata))
455 [ - + ]: 71059 : ABI_MALLOC(now2,(mdata))
456 [ - + ]: 71059 : ABI_MALLOC(before2,(mdata))
457 [ - + + - : 213177 : ABI_MALLOC(trig3,(2,n3))
+ - + - -
+ - + ]
458 [ - + ]: 71059 : ABI_MALLOC(after3,(mdata))
459 [ - + ]: 71059 : ABI_MALLOC(now3,(mdata))
460 [ - + ]: 71059 : ABI_MALLOC(before3,(mdata))
461 [ - + + - : 284236 : ABI_MALLOC(zw,(2,ncache/4,2))
+ - + - -
+ - + ]
462 [ - + + - : 284236 : ABI_MALLOC(zt,(2,lzt,n1zt))
- + + - +
- + - - +
- + ]
463 [ - + + - : 355295 : ABI_MALLOC(zmpi2,(2,n1,nd2proc,nnd3))
- + + - -
+ + - + -
+ - - + -
+ ]
464 [ + + ]: 71059 : if (nproc_fft>1) then
465 [ - + + - : 58560 : ABI_MALLOC(zmpi1,(2,n1,nd2proc,nnd3))
- + + - -
+ + - - +
- + ]
466 : : end if
467 : :
468 : 71059 : call ctrig(n2,trig2,after2,before2,now2,-1,ic2)
469 : 71059 : call ctrig(n1,trig1,after1,before1,now1,-1,ic1)
470 : 71059 : call ctrig(n3,trig3,after3,before3,now3,-1,ic3)
471 : :
472 [ + + ]: 142126 : do idat=1,ndat
473 [ + + ]: 1063870 : do j3=1,nd3proc
474 [ + + ]: 1063870 : if (me_fft*(nd3proc)+j3 <= n3) then
475 : 991030 : Jp2st=1; J2st=1
476 : :
477 : : ! transform along y axis
478 : : ! input: R1,R2,R3,(Rp3)
479 : 991030 : lot=ncache/(4*n2)
480 : :
481 [ - + - - : 2766492 : do j=1,n1eff,lot
- + + + ]
482 : 1775462 : ma=j
483 : 1775462 : mb=min(j+(lot-1),n1eff)
484 : 1775462 : n1dfft=mb-ma+1
485 : 1775462 : i=1
486 : : call fftstp(nd1eff,n1dfft,nd2,lot,n2,zr(1,j,1,j3,idat),zw(1,1,1), &
487 : 1775462 : & trig2,after2(i),now2(i),before2(i),-1)
488 : :
489 : 1775462 : inzee=1
490 [ + + ]: 4578983 : do i=2,ic2
491 : : call fftstp(lot,n1dfft,n2,lot,n2,zw(1,1,inzee),zw(1,1,3-inzee), &
492 : 2803521 : & trig2,after2(i),now2(i),before2(i),-1)
493 : 4578983 : inzee=3-inzee
494 : : end do
495 : :
496 : : ! input: R1,G2,R3,(Rp3)
497 : : ! output: G2,R1,R3,(Rp3)
498 [ + + ]: 2766492 : if(cplex==2)then
499 : 1769540 : 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 : 991030 : lot=ncache/(4*n1)
508 : :
509 [ - + - - : 2766492 : do j=1,n2eff,lot
- + + + ]
510 : 1775462 : ma=j
511 : 1775462 : mb=min(j+(lot-1),n2eff)
512 : 1775462 : n1dfft=mb-ma+1
513 : :
514 : 1775462 : i=1
515 : : call fftstp(lzt,n1dfft,n1zt,lot,n1,zt(1,j,1),zw(1,1,1), &
516 : 1775462 : & trig1,after1(i),now1(i),before1(i),-1)
517 : :
518 : 1775462 : inzee=1
519 [ + + ]: 4579223 : do i=2,ic1
520 : : call fftstp(lot,n1dfft,n1,lot,n1,zw(1,1,inzee),zw(1,1,3-inzee), &
521 : 2803761 : & trig1,after1(i),now1(i),before1(i),-1)
522 : 4579223 : 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 [ + + ]: 2766492 : 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 : 923844 : 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 [ + + ]: 71067 : if (nproc_fft>1) then
543 : 14640 : call timab(544,1,tsec)
544 : : call xmpi_alltoall(zmpi1,2*n1*nd2proc*nd3proc, &
545 : 14640 : & zmpi2,2*n1*nd2proc*nd3proc,comm_fft,ierr)
546 : 14640 : call timab(544,2,tsec)
547 : : end if
548 : :
549 : : ! transform along z axis
550 : : ! input: G1,G2,R3,(Gp2)
551 : 71067 : lot=ncache/(4*n3)
552 : :
553 [ + + ]: 1128162 : do j2=1,nd2proc
554 [ + + ]: 1057103 : if (me_fft*(nd2proc)+j2 <= n2eff) then
555 [ - + - - : 2759032 : do i1=1,n1,lot
- + + + ]
556 : 1773001 : ma=i1
557 : 1773001 : mb=min(i1+(lot-1),n1)
558 : 1773001 : n1dfft=mb-ma+1
559 : :
560 : : ! input: G1,G2,R3,(Gp2)
561 : : ! output: G1,R3,G2,(Gp2)
562 : 1773001 : call unscramble(i1,j2,lot,n1dfft,n1,n3,nd2proc,nd3,zmpi2,zw(1,1,1))
563 : :
564 : 1773001 : inzee=1
565 [ + + ]: 6365019 : do i=1,ic3
566 : : call fftstp(lot,n1dfft,n3,lot,n3,zw(1,1,inzee),zw(1,1,3-inzee), &
567 : 4592018 : & trig3,after3(i),now3(i),before3(i),-1)
568 : 6365019 : inzee=3-inzee
569 : : end do
570 : :
571 : 2759032 : 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 [ - + ]: 71059 : ABI_FREE(trig1)
580 [ - + ]: 71059 : ABI_FREE(after1)
581 [ - + ]: 71059 : ABI_FREE(now1)
582 [ - + ]: 71059 : ABI_FREE(before1)
583 [ - + ]: 71059 : ABI_FREE(trig2)
584 [ - + ]: 71059 : ABI_FREE(after2)
585 [ - + ]: 71059 : ABI_FREE(now2)
586 [ - + ]: 71059 : ABI_FREE(before2)
587 [ - + ]: 71059 : ABI_FREE(trig3)
588 [ - + ]: 71059 : ABI_FREE(after3)
589 [ - + ]: 71059 : ABI_FREE(now3)
590 [ - + ]: 71059 : ABI_FREE(before3)
591 [ - + ]: 71059 : ABI_FREE(zmpi2)
592 [ - + ]: 71059 : ABI_FREE(zw)
593 [ - + ]: 71059 : ABI_FREE(zt)
594 [ + + ]: 71059 : if (nproc_fft>1) then
595 [ - + ]: 14640 : ABI_FREE(zmpi1)
596 : : end if
597 : :
598 [ - + ]: 71059 : 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 : 800864 : subroutine sg2002_mpiback_wf(icplexwf,ndat,n1,n2,n3,nd1,nd2,nd3proc,&
662 : 800864 : & 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 : 800864 : real(dp),allocatable :: zt(:,:,:) ! work arrays for transpositions
677 : 800864 : real(dp),allocatable :: zmpi1(:,:,:,:,:),zmpi2(:,:,:,:,:) ! work arrays for MPI
678 : 800864 : real(dp),allocatable :: zw(:,:,:) ! cache work array
679 : : ! FFT work arrays
680 : 800864 : real(dp),allocatable :: trig1(:,:),trig2(:,:),trig3(:,:)
681 : 800864 : integer,allocatable :: after1(:),now1(:),before1(:),after2(:)
682 : 800864 : 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 : 800864 : ioption = 1
691 : : !if (paral_kgb==1) ioption=1
692 : :
693 : 800864 : 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 : 800864 : ncache=4*max(n1,n2,n3,1024)
697 [ - + ]: 800864 : 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 : 800864 : n1eff=n1; m2eff=m2; m1zt=n1
707 [ + + ]: 800864 : if (icplexwf==1) then
708 : 51 : n1eff=(n1+1)/2; m2eff=m2/2+1; m1zt=2*(n1/2+1)
709 : : end if
710 : :
711 : 800864 : lzt=m2eff
712 [ + + ]: 800864 : if (mod(m2eff,2)==0) lzt=lzt+1
713 [ + + ]: 703803 : if (mod(m2eff,4)==0) lzt=lzt+1
714 : :
715 : : ! maximal number of big box 3rd dim slices for all procs
716 : 800864 : nnd3=nd3proc*nproc_fft
717 : :
718 [ - + + - : 2402592 : ABI_MALLOC(trig1,(2,n1))
+ - + - -
+ - + ]
719 [ - + ]: 800864 : ABI_MALLOC(after1,(mdata))
720 [ - + ]: 800864 : ABI_MALLOC(now1,(mdata))
721 [ - + ]: 800864 : ABI_MALLOC(before1,(mdata))
722 [ - + + - : 2402592 : ABI_MALLOC(trig2,(2,n2))
+ - + - -
+ - + ]
723 [ - + ]: 800864 : ABI_MALLOC(after2,(mdata))
724 [ - + ]: 800864 : ABI_MALLOC(now2,(mdata))
725 [ - + ]: 800864 : ABI_MALLOC(before2,(mdata))
726 [ - + + - : 2402592 : ABI_MALLOC(trig3,(2,n3))
+ - + - -
+ - + ]
727 [ - + ]: 800864 : ABI_MALLOC(after3,(mdata))
728 [ - + ]: 800864 : ABI_MALLOC(now3,(mdata))
729 [ - + ]: 800864 : ABI_MALLOC(before3,(mdata))
730 : :
731 : : ! Allocate cache work array and work arrays for MPI transpositions.
732 [ - + + - : 3203456 : ABI_MALLOC(zw,(2,ncache/4,2))
+ - + - -
+ - + ]
733 [ - + + - : 3203456 : ABI_MALLOC(zt,(2,lzt,m1zt))
- + + - +
- + - - +
- + ]
734 [ - + + - : 4805184 : ABI_MALLOC(zmpi2,(2,md1,md2proc,nnd3,ndat))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
735 [ + + ]: 800864 : if (nproc_fft>1) then
736 [ - + + - : 3495210 : ABI_MALLOC(zmpi1,(2,md1,md2proc,nnd3,ndat))
- + + - -
+ + - - +
+ - - + -
+ ]
737 : : end if
738 : :
739 : : ! Compute twiddle coefficients.
740 : 800864 : call ctrig(n3,trig3,after3,before3,now3,1,ic3)
741 : 800864 : call ctrig(n1,trig1,after1,before1,now1,1,ic1)
742 : 800864 : 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 [ + + ]: 1655281 : do idat=1,ndat
754 : :
755 : : ! transform along z axis
756 : : ! input: G1,G3,G2,(Gp2)
757 : 854417 : lot=ncache/(4*n3)
758 : :
759 [ + + + + : 5252101299 : zw(:,:,:)=zero
+ + ]
760 [ + + + + : 675738053 : zt(:,:,:)=zero
+ + ]
761 : :
762 : : ! Loop over the y planes treated by this node and trasform n1ddft G_z lines.
763 [ + + ]: 5563180 : do j2=1,md2proc
764 : :
765 : : ! if (me_fft*md2proc+j2<=m2eff) then !a faire plus tard
766 : :
767 [ - + - - : 10639876 : do i1=1,m1,lot
- + + + ]
768 : 5076696 : ma=i1
769 : 5076696 : mb=min(i1+(lot-1),m1)
770 : 5076696 : 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 : 5076696 : 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 : 5076696 : inzee=1
779 [ + + ]: 15832187 : do i=1,ic3
780 : : call fftstp(lot,n1dfft,n3,lot,n3,zw(1,1,inzee),zw(1,1,3-inzee), &
781 : 10755491 : & trig3,after3(i),now3(i),before3(i),1)
782 : 15832187 : inzee=3-inzee
783 : : end do
784 : :
785 : : ! Local rotation.
786 : : ! input: G1,R3,G2,(Gp2)
787 : : ! output: G1,G2,R3,(Gp2)
788 : 9785459 : 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 [ + + ]: 854417 : if (nproc_fft>1) then
797 : 720698 : call timab(543,1,tsec)
798 : : call xmpi_alltoall(zmpi2(:,:,:,:,idat),2*md1*md2proc*nd3proc, &
799 : 720698 : & zmpi1(:,:,:,:,idat),2*md1*md2proc*nd3proc,comm_fft,ierr)
800 : 720698 : call timab(543,2,tsec)
801 : : end if
802 : :
803 : : ! Loop over the z treated by this node.
804 [ + + ]: 11428834 : do j3=1,nd3proc
805 : : !j3glob = j3 + me_fft*nd3proc
806 [ + - ]: 10627970 : if (me_fft*nd3proc+j3 <= n3) then
807 : 9773553 : Jp2st=1; J2st=1
808 : :
809 : 9773553 : lot=ncache/(4*n1)
810 : :
811 : : ! Loop over G_y in the small box.
812 [ - + - - : 20288834 : do j=1,m2eff,lot
- + + + ]
813 : 10515281 : ma=j
814 : 10515281 : mb=min(j+(lot-1),m2eff)
815 : 10515281 : n1dfft=mb-ma+1
816 : :
817 : : ! Zero-pad input.
818 : : ! input: G1,G2,R3,JG2,(Rp3)
819 : : ! output: G2,G1,R3,JG2,(Rp3)
820 [ + + ]: 10515281 : 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 : 7886350 : & 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 : 10515281 : inzee=1
832 [ + + ]: 22494564 : do i=1,ic1-1
833 : : call fftstp(lot,n1dfft,n1,lot,n1,zw(1,1,inzee),zw(1,1,3-inzee), &
834 : 11979283 : & trig1,after1(i),now1(i),before1(i),1)
835 : 22494564 : inzee=3-inzee
836 : : end do
837 : :
838 : 10515281 : i=ic1
839 : : call fftstp(lot,n1dfft,n1,lzt,m1zt,zw(1,1,inzee),zt(1,j,1), &
840 : 20288834 : & 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 : 9773553 : lot=ncache/(4*n2)
846 : :
847 [ - + - - : 21247926 : do j=1,n1eff,lot
- + + + ]
848 : 11474373 : ma=j
849 : 11474373 : mb=min(j+(lot-1),n1eff)
850 : 11474373 : n1dfft=mb-ma+1
851 : 11474373 : includelast=1
852 [ + + ]: 11474373 : 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 [ + + ]: 11474373 : if (icplexwf==2) then
861 : 11459291 : 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 : 11474373 : inzee=1
869 [ + + ]: 25130130 : do i=1,ic2-1
870 : : call fftstp(lot,n1dfft,n2,lot,n2,zw(1,1,inzee),zw(1,1,3-inzee), &
871 : 13655757 : & trig2,after2(i),now2(i),before2(i),1)
872 : 25130130 : inzee=3-inzee
873 : : end do
874 : :
875 : 11474373 : i=ic2
876 : :
877 : : call fftstp(lot,n1dfft,n2,nd1,nd2,zw(1,1,inzee),zr(1,j,1,j3,idat), &
878 : 21247926 : & trig2,after2(i),now2(i),before2(i),1)
879 : :
880 : :
881 : : end do
882 : :
883 : : ! Treat real wavefunctions.
884 [ + + ]: 9773553 : 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 [ - + ]: 800864 : ABI_FREE(trig1)
909 [ - + ]: 800864 : ABI_FREE(after1)
910 [ - + ]: 800864 : ABI_FREE(now1)
911 [ - + ]: 800864 : ABI_FREE(before1)
912 [ - + ]: 800864 : ABI_FREE(trig2)
913 [ - + ]: 800864 : ABI_FREE(after2)
914 [ - + ]: 800864 : ABI_FREE(now2)
915 [ - + ]: 800864 : ABI_FREE(before2)
916 [ - + ]: 800864 : ABI_FREE(trig3)
917 [ - + ]: 800864 : ABI_FREE(after3)
918 [ - + ]: 800864 : ABI_FREE(now3)
919 [ - + ]: 800864 : ABI_FREE(before3)
920 [ - + ]: 800864 : ABI_FREE(zmpi2)
921 [ - + ]: 800864 : ABI_FREE(zw)
922 [ - + ]: 800864 : ABI_FREE(zt)
923 [ + + ]: 800864 : if (nproc_fft>1) then
924 [ - + ]: 699042 : ABI_FREE(zmpi1)
925 : : end if
926 : :
927 : : !call timab(541,2,tsec)
928 : :
929 [ - + ]: 800864 : 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 : 737556 : subroutine sg2002_mpiforw_wf(icplexwf,ndat,n1,n2,n3,nd1,nd2,nd3proc,&
994 : 737556 : & 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 : 737556 : real(dp), allocatable :: zt(:,:,:) ! work arrays for transpositions
1012 : 737556 : real(dp), allocatable :: zmpi1(:,:,:,:,:),zmpi2(:,:,:,:,:) ! work arrays for MPI
1013 : 737556 : real(dp), allocatable :: zw(:,:,:) ! cache work array
1014 : : ! FFT work arrays
1015 : 737556 : real(dp), allocatable :: trig1(:,:),trig2(:,:),trig3(:,:)
1016 : 737556 : 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 : 737556 : ioption = 1
1026 : : !if (paral_kgb==1) ioption=1
1027 : :
1028 : 737556 : 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 : 737556 : ncache=4*max(n1,n2,n3,1024)
1032 [ - + ]: 737556 : 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 : 737556 : n1eff=n1; m2eff=m2; m1zt=n1
1042 [ + + ]: 737556 : if (icplexwf==1) then
1043 : 35 : n1eff=(n1+1)/2; m2eff=m2/2+1; m1zt=2*(n1/2+1)
1044 : : end if
1045 : :
1046 : 737556 : lzt=m2eff
1047 [ + + ]: 737556 : if (mod(m2eff,2)==0) lzt=lzt+1
1048 [ + + ]: 654763 : if (mod(m2eff,4)==0) lzt=lzt+1
1049 : :
1050 : : ! maximal number of big box 3rd dim slices for all procs
1051 : 737556 : nnd3=nd3proc*nproc_fft
1052 : :
1053 [ - + + - : 2212668 : ABI_MALLOC(trig1,(2,n1))
+ - + - -
+ - + ]
1054 [ - + ]: 737556 : ABI_MALLOC(after1,(mdata))
1055 [ - + ]: 737556 : ABI_MALLOC(now1,(mdata))
1056 [ - + ]: 737556 : ABI_MALLOC(before1,(mdata))
1057 [ - + + - : 2212668 : ABI_MALLOC(trig2,(2,n2))
+ - + - -
+ - + ]
1058 [ - + ]: 737556 : ABI_MALLOC(after2,(mdata))
1059 [ - + ]: 737556 : ABI_MALLOC(now2,(mdata))
1060 [ - + ]: 737556 : ABI_MALLOC(before2,(mdata))
1061 [ - + + - : 2212668 : ABI_MALLOC(trig3,(2,n3))
+ - + - -
+ - + ]
1062 [ - + ]: 737556 : ABI_MALLOC(after3,(mdata))
1063 [ - + ]: 737556 : ABI_MALLOC(now3,(mdata))
1064 [ - + ]: 737556 : ABI_MALLOC(before3,(mdata))
1065 [ - + + - : 2950224 : ABI_MALLOC(zw,(2,ncache/4,2))
+ - + - -
+ - + ]
1066 [ - + + - : 2950224 : ABI_MALLOC(zt,(2,lzt,m1zt))
- + + - +
- + - - +
- + ]
1067 [ - + + - : 4425336 : ABI_MALLOC(zmpi2,(2,md1,md2proc,nnd3,ndat))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
1068 [ + + ]: 737556 : if (nproc_fft>1) then
1069 [ - + + - : 3255040 : ABI_MALLOC(zmpi1,(2,md1,md2proc,nnd3,ndat))
- + + - -
+ + - - +
+ - - + -
+ ]
1070 : : end if
1071 : :
1072 : 737556 : call ctrig(n2,trig2,after2,before2,now2,-1,ic2)
1073 : 737556 : call ctrig(n1,trig1,after1,before1,now1,-1,ic1)
1074 : 737556 : 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 [ + + ]: 1528655 : do idat=1,ndat
1086 : : ! Loop over the z-planes treated by this node
1087 [ + + ]: 9799417 : do j3=1,nd3proc
1088 : :
1089 [ + - ]: 9799417 : if (me_fft*nd3proc+j3 <= n3) then
1090 : 9008318 : Jp2st=1
1091 : 9008318 : J2st=1
1092 : :
1093 : : ! Treat real wavefunctions.
1094 [ + + ]: 9008318 : 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 : 9008318 : lot=ncache/(4*n2)
1115 : :
1116 [ - + - - : 19563916 : do j=1,n1eff,lot
- + + + ]
1117 : 10555598 : ma=j
1118 : 10555598 : mb=min(j+(lot-1),n1eff)
1119 : 10555598 : n1dfft=mb-ma+1
1120 : 10555598 : i=1
1121 : : call fftstp(nd1,n1dfft,nd2,lot,n2,zr(1,j,1,j3,idat),zw(1,1,1), &
1122 : 10555598 : & trig2,after2(i),now2(i),before2(i),-1)
1123 : :
1124 : 10555598 : inzee=1
1125 [ + + ]: 23092280 : do i=2,ic2
1126 : : call fftstp(lot,n1dfft,n2,lot,n2,zw(1,1,inzee),zw(1,1,3-inzee), &
1127 : 12536682 : & trig2,after2(i),now2(i),before2(i),-1)
1128 : 23092280 : inzee=3-inzee
1129 : : end do
1130 : :
1131 : : ! input: R1,G2,R3,(Rp3)
1132 : : ! output: G2,R1,R3,(Rp3)
1133 [ + + ]: 19563916 : if(icplexwf==2)then
1134 : 10544372 : 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 : 9008318 : lot=ncache/(4*n1)
1143 : :
1144 [ - + - - : 18688644 : do j=1,m2eff,lot
- + + + ]
1145 : 9680326 : ma=j
1146 : 9680326 : mb=min(j+(lot-1),m2eff)
1147 : 9680326 : n1dfft=mb-ma+1
1148 : 9680326 : i=1
1149 : : call fftstp(lzt,n1dfft,m1zt,lot,n1,zt(1,j,1),zw(1,1,1), &
1150 : 9680326 : & trig1,after1(i),now1(i),before1(i),-1)
1151 : :
1152 : 9680326 : inzee=1
1153 [ + + ]: 20701844 : do i=2,ic1
1154 : : call fftstp(lot,n1dfft,n1,lot,n1,zw(1,1,inzee),zw(1,1,3-inzee), &
1155 : 11021518 : & trig1,after1(i),now1(i),before1(i),-1)
1156 : 20701844 : 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 [ + + ]: 18688644 : 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 : 7376560 : & 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 [ + + ]: 791099 : if (nproc_fft>1) then
1178 : 672664 : call timab(544,1,tsec)
1179 : : call xmpi_alltoall(zmpi1(:,:,:,:,idat),2*md1*md2proc*nd3proc, &
1180 : 672664 : & zmpi2(:,:,:,:,idat),2*md1*md2proc*nd3proc,comm_fft,ierr)
1181 : :
1182 : 672664 : call timab(544,2,tsec)
1183 : : end if
1184 : :
1185 : : ! transform along z axis
1186 : : ! input: G1,G2,R3,(Gp2)
1187 : 791099 : lot=ncache/(4*n3)
1188 : :
1189 [ + + ]: 5135480 : do j2=1,md2proc
1190 [ + + ]: 5135480 : 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 [ - + - - : 9018376 : do i1=1,m1,lot
- + + + ]
1193 : 4674620 : ma=i1
1194 : 4674620 : mb=min(i1+(lot-1),m1)
1195 : 4674620 : n1dfft=mb-ma+1
1196 : :
1197 : : ! input: G1,G2,R3,(Gp2)
1198 : : ! output: G1,R3,G2,(Gp2)
1199 : 4674620 : call unscramble(i1,j2,lot,n1dfft,md1,n3,md2proc,nnd3,zmpi2(:,:,:,:,idat),zw(1,1,1))
1200 : :
1201 : 4674620 : inzee=1
1202 [ + + ]: 14566053 : do i=1,ic3
1203 : : call fftstp(lot,n1dfft,n3,lot,n3,zw(1,1,inzee),zw(1,1,3-inzee), &
1204 : 9891433 : & trig3,after3(i),now3(i),before3(i),-1)
1205 : 14566053 : inzee=3-inzee
1206 : : end do
1207 : :
1208 : 9018376 : 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 [ + + ]: 1528655 : 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 [ - + ]: 737556 : ABI_FREE(trig1)
1239 [ - + ]: 737556 : ABI_FREE(after1)
1240 [ - + ]: 737556 : ABI_FREE(now1)
1241 [ - + ]: 737556 : ABI_FREE(before1)
1242 [ - + ]: 737556 : ABI_FREE(trig2)
1243 [ - + ]: 737556 : ABI_FREE(after2)
1244 [ - + ]: 737556 : ABI_FREE(now2)
1245 [ - + ]: 737556 : ABI_FREE(before2)
1246 [ - + ]: 737556 : ABI_FREE(trig3)
1247 [ - + ]: 737556 : ABI_FREE(after3)
1248 [ - + ]: 737556 : ABI_FREE(now3)
1249 [ - + ]: 737556 : ABI_FREE(before3)
1250 [ - + ]: 737556 : ABI_FREE(zmpi2)
1251 [ - + ]: 737556 : ABI_FREE(zw)
1252 [ - + ]: 737556 : ABI_FREE(zt)
1253 [ + + ]: 737556 : if (nproc_fft>1) then
1254 [ - + ]: 651008 : ABI_FREE(zmpi1)
1255 : : end if
1256 : :
1257 : : !call timab(542,2,tsec)
1258 : :
1259 [ - + ]: 737556 : 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 : 137428 : subroutine sg2002_mpifourdp(cplex,nfft,ngfft,ndat,isign,&
1299 : 137428 : & 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 : 137428 : real(dp),allocatable :: workf(:,:,:,:,:),workr(:,:,:,:,:)
1315 : :
1316 : : ! *************************************************************************
1317 : :
1318 : : ! Note the only c2c is supported in parallel.
1319 : 137428 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
1320 : 137428 : n4=ngfft(4); n5=ngfft(5); n6=ngfft(6)
1321 : 137428 : me_fft=ngfft(11); nproc_fft=ngfft(10)
1322 : :
1323 : 137428 : nd2proc=((n2-1)/nproc_fft) +1
1324 : 137428 : nd3proc=((n6-1)/nproc_fft) +1
1325 [ - + + - : 824568 : ABI_MALLOC(workr,(2,n4,n5,nd3proc,ndat))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
1326 [ - + + - : 824568 : ABI_MALLOC(workf,(2,n4,n6,nd2proc,ndat))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
1327 : :
1328 : : ! Complex to Complex
1329 : 206227 : select case (isign)
1330 : : case (1)
1331 : : ! G --> R
1332 : 68799 : call mpifft_fg2dbox(nfft,ndat,fofg,n1,n2,n3,n4,nd2proc,n6,fftn2_distrib,ffti2_local,me_fft,workf)
1333 : :
1334 : 68799 : call sg2002_back(2,ndat,n1,n2,n3,n4,n5,n6,n4,nd2proc,nd3proc,2,workf,workr,comm_fft)
1335 : :
1336 : 68799 : 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 : 68629 : call mpifft_fr2dbox(cplex,nfft,ndat,fofr,n1,n2,n3,n4,n5,nd3proc,fftn3_distrib,ffti3_local,me_fft,workr)
1341 : :
1342 : 68629 : 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 : 68629 : call mpifft_dbox2fg(n1,n2,n3,n4,nd2proc,n6,ndat,fftn2_distrib,ffti2_local,me_fft,workf,nfft,fofg)
1346 : :
1347 : : case default
1348 [ + + - ]: 137428 : ABI_BUG("Wrong isign")
1349 : : end select
1350 : :
1351 [ - + ]: 137428 : ABI_FREE(workr)
1352 [ - + ]: 137428 : ABI_FREE(workf)
1353 : :
1354 : 137428 : 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 : 5046474 : 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 [ + - ]: 59425342 : do i=1,ndata
2513 [ + + ]: 59425342 : if (n.eq.ifftdata(1,i)) then
2514 : 5046474 : ic=0
2515 [ - + ]: 15237119 : do j=1,(mdata-1)
2516 : 15237119 : itt=ifftdata(1+j,i)
2517 [ + + ]: 15237119 : if (itt.gt.1) then
2518 : 10190645 : ic=ic+1
2519 : 10190645 : 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 : 5046474 : ABI_ERROR("Aborting now")
2532 : :
2533 : : 1000 continue
2534 : 5046474 : after(1)=1
2535 : 5046474 : before(ic)=1
2536 [ + + ]: 10190645 : do i=2,ic
2537 : 5144171 : after(i)=after(i-1)*now(i-1)
2538 : 10190645 : before(ic-i+1)=before(ic-i+2)*now(ic-i+2)
2539 : : end do
2540 : :
2541 : 5046474 : angle=isign*two_pi/n
2542 [ + + ]: 5046474 : if (mod(n,2).eq.0) then
2543 : 4990975 : nh=n/2
2544 : 4990975 : trig(1,1)=one
2545 : 4990975 : trig(2,1)=zero
2546 : 4990975 : trig(1,nh+1)=-one
2547 : 4990975 : trig(2,nh+1)=zero
2548 [ + + ]: 50884149 : do i=1,nh-1
2549 : 45893174 : trigc=cos(i*angle)
2550 : 45893174 : trigs=sin(i*angle)
2551 : 45893174 : trig(1,i+1)=trigc
2552 : 45893174 : trig(2,i+1)=trigs
2553 : 45893174 : trig(1,n-i+1)=trigc
2554 : 50884149 : 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 : 5046474 : 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 : 139623218 : 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 : 139623218 : atn=after*now
2613 : 139623218 : atb=after*before
2614 : :
2615 : : ! sqrt(.5d0)
2616 [ + + + + : 139623218 : 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 [ + + ]: 34170863 : if (isign.eq.1) then
2795 : 17581258 : ia=1
2796 : 17581258 : nin1=ia-after
2797 : 17581258 : nout1=ia-atn
2798 [ + + ]: 70219554 : do ib=1,before
2799 : 52638296 : nin1=nin1+after
2800 : 52638296 : nin2=nin1+atb
2801 : 52638296 : nin3=nin2+atb
2802 : 52638296 : nin4=nin3+atb
2803 : 52638296 : nout1=nout1+atn
2804 : 52638296 : nout2=nout1+after
2805 : 52638296 : nout3=nout2+after
2806 : 52638296 : nout4=nout3+after
2807 [ + + ]: 785407455 : do j=1,n1dfft
2808 : 715187901 : r1=zin(1,j,nin1)
2809 : 715187901 : s1=zin(2,j,nin1)
2810 : 715187901 : r2=zin(1,j,nin2)
2811 : 715187901 : s2=zin(2,j,nin2)
2812 : 715187901 : r3=zin(1,j,nin3)
2813 : 715187901 : s3=zin(2,j,nin3)
2814 : 715187901 : r4=zin(1,j,nin4)
2815 : 715187901 : s4=zin(2,j,nin4)
2816 : 715187901 : r=r1 + r3
2817 : 715187901 : s=r2 + r4
2818 : 715187901 : zout(1,j,nout1) = r + s
2819 : 715187901 : zout(1,j,nout3) = r - s
2820 : 715187901 : r=r1 - r3
2821 : 715187901 : s=s2 - s4
2822 : 715187901 : zout(1,j,nout2) = r - s
2823 : 715187901 : zout(1,j,nout4) = r + s
2824 : 715187901 : r=s1 + s3
2825 : 715187901 : s=s2 + s4
2826 : 715187901 : zout(2,j,nout1) = r + s
2827 : 715187901 : zout(2,j,nout3) = r - s
2828 : 715187901 : r=s1 - s3
2829 : 715187901 : s=r2 - r4
2830 : 715187901 : zout(2,j,nout2) = r + s
2831 : 767826197 : zout(2,j,nout4) = r - s
2832 : : enddo
2833 : : enddo
2834 [ + + ]: 90154048 : do 4000,ia=2,after
2835 : 72572790 : ias=ia-1
2836 [ + + ]: 72572790 : 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 : 67734356 : itt=ias*before
2881 : 67734356 : itrig=itt+1
2882 : 67734356 : cr2=trig(1,itrig)
2883 : 67734356 : ci2=trig(2,itrig)
2884 : 67734356 : itrig=itrig+itt
2885 : 67734356 : cr3=trig(1,itrig)
2886 : 67734356 : ci3=trig(2,itrig)
2887 : 67734356 : itrig=itrig+itt
2888 : 67734356 : cr4=trig(1,itrig)
2889 : 67734356 : ci4=trig(2,itrig)
2890 : 67734356 : nin1=ia-after
2891 : 67734356 : nout1=ia-atn
2892 [ + + ]: 141780240 : do ib=1,before
2893 : 74045884 : nin1=nin1+after
2894 : 74045884 : nin2=nin1+atb
2895 : 74045884 : nin3=nin2+atb
2896 : 74045884 : nin4=nin3+atb
2897 : 74045884 : nout1=nout1+atn
2898 : 74045884 : nout2=nout1+after
2899 : 74045884 : nout3=nout2+after
2900 : 74045884 : nout4=nout3+after
2901 [ + + ]: 1028281902 : do j=1,n1dfft
2902 : 886501662 : r1=zin(1,j,nin1)
2903 : 886501662 : s1=zin(2,j,nin1)
2904 : 886501662 : r=zin(1,j,nin2)
2905 : 886501662 : s=zin(2,j,nin2)
2906 : 886501662 : r2=r*cr2 - s*ci2
2907 : 886501662 : s2=r*ci2 + s*cr2
2908 : 886501662 : r=zin(1,j,nin3)
2909 : 886501662 : s=zin(2,j,nin3)
2910 : 886501662 : r3=r*cr3 - s*ci3
2911 : 886501662 : s3=r*ci3 + s*cr3
2912 : 886501662 : r=zin(1,j,nin4)
2913 : 886501662 : s=zin(2,j,nin4)
2914 : 886501662 : r4=r*cr4 - s*ci4
2915 : 886501662 : s4=r*ci4 + s*cr4
2916 : 886501662 : r=r1 + r3
2917 : 886501662 : s=r2 + r4
2918 : 886501662 : zout(1,j,nout1) = r + s
2919 : 886501662 : zout(1,j,nout3) = r - s
2920 : 886501662 : r=r1 - r3
2921 : 886501662 : s=s2 - s4
2922 : 886501662 : zout(1,j,nout2) = r - s
2923 : 886501662 : zout(1,j,nout4) = r + s
2924 : 886501662 : r=s1 + s3
2925 : 886501662 : s=s2 + s4
2926 : 886501662 : zout(2,j,nout1) = r + s
2927 : 886501662 : zout(2,j,nout3) = r - s
2928 : 886501662 : r=s1 - s3
2929 : 886501662 : s=r2 - r4
2930 : 886501662 : zout(2,j,nout2) = r + s
2931 : 960547546 : zout(2,j,nout4) = r - s
2932 : : enddo
2933 : : enddo
2934 : : end if
2935 : 17581258 : 4000 continue
2936 : : else
2937 : 16589605 : ia=1
2938 : 16589605 : nin1=ia-after
2939 : 16589605 : nout1=ia-atn
2940 [ + + ]: 67415150 : do ib=1,before
2941 : 50825545 : nin1=nin1+after
2942 : 50825545 : nin2=nin1+atb
2943 : 50825545 : nin3=nin2+atb
2944 : 50825545 : nin4=nin3+atb
2945 : 50825545 : nout1=nout1+atn
2946 : 50825545 : nout2=nout1+after
2947 : 50825545 : nout3=nout2+after
2948 : 50825545 : nout4=nout3+after
2949 [ + + ]: 771329250 : do j=1,n1dfft
2950 : 703914100 : r1=zin(1,j,nin1)
2951 : 703914100 : s1=zin(2,j,nin1)
2952 : 703914100 : r2=zin(1,j,nin2)
2953 : 703914100 : s2=zin(2,j,nin2)
2954 : 703914100 : r3=zin(1,j,nin3)
2955 : 703914100 : s3=zin(2,j,nin3)
2956 : 703914100 : r4=zin(1,j,nin4)
2957 : 703914100 : s4=zin(2,j,nin4)
2958 : 703914100 : r=r1 + r3
2959 : 703914100 : s=r2 + r4
2960 : 703914100 : zout(1,j,nout1) = r + s
2961 : 703914100 : zout(1,j,nout3) = r - s
2962 : 703914100 : r=r1 - r3
2963 : 703914100 : s=s2 - s4
2964 : 703914100 : zout(1,j,nout2) = r + s
2965 : 703914100 : zout(1,j,nout4) = r - s
2966 : 703914100 : r=s1 + s3
2967 : 703914100 : s=s2 + s4
2968 : 703914100 : zout(2,j,nout1) = r + s
2969 : 703914100 : zout(2,j,nout3) = r - s
2970 : 703914100 : r=s1 - s3
2971 : 703914100 : s=r2 - r4
2972 : 703914100 : zout(2,j,nout2) = r - s
2973 : 754739645 : zout(2,j,nout4) = r + s
2974 : : enddo
2975 : : enddo
2976 [ + + ]: 95980785 : do 4100,ia=2,after
2977 : 79391180 : ias=ia-1
2978 [ + + ]: 79391180 : 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 : 75146200 : itt=ias*before
3023 : 75146200 : itrig=itt+1
3024 : 75146200 : cr2=trig(1,itrig)
3025 : 75146200 : ci2=trig(2,itrig)
3026 : 75146200 : itrig=itrig+itt
3027 : 75146200 : cr3=trig(1,itrig)
3028 : 75146200 : ci3=trig(2,itrig)
3029 : 75146200 : itrig=itrig+itt
3030 : 75146200 : cr4=trig(1,itrig)
3031 : 75146200 : ci4=trig(2,itrig)
3032 : 75146200 : nin1=ia-after
3033 : 75146200 : nout1=ia-atn
3034 [ + + ]: 155564352 : do ib=1,before
3035 : 80418152 : nin1=nin1+after
3036 : 80418152 : nin2=nin1+atb
3037 : 80418152 : nin3=nin2+atb
3038 : 80418152 : nin4=nin3+atb
3039 : 80418152 : nout1=nout1+atn
3040 : 80418152 : nout2=nout1+after
3041 : 80418152 : nout3=nout2+after
3042 : 80418152 : nout4=nout3+after
3043 [ + + ]: 1089478268 : do j=1,n1dfft
3044 : 933913916 : r1=zin(1,j,nin1)
3045 : 933913916 : s1=zin(2,j,nin1)
3046 : 933913916 : r=zin(1,j,nin2)
3047 : 933913916 : s=zin(2,j,nin2)
3048 : 933913916 : r2=r*cr2 - s*ci2
3049 : 933913916 : s2=r*ci2 + s*cr2
3050 : 933913916 : r=zin(1,j,nin3)
3051 : 933913916 : s=zin(2,j,nin3)
3052 : 933913916 : r3=r*cr3 - s*ci3
3053 : 933913916 : s3=r*ci3 + s*cr3
3054 : 933913916 : r=zin(1,j,nin4)
3055 : 933913916 : s=zin(2,j,nin4)
3056 : 933913916 : r4=r*cr4 - s*ci4
3057 : 933913916 : s4=r*ci4 + s*cr4
3058 : 933913916 : r=r1 + r3
3059 : 933913916 : s=r2 + r4
3060 : 933913916 : zout(1,j,nout1) = r + s
3061 : 933913916 : zout(1,j,nout3) = r - s
3062 : 933913916 : r=r1 - r3
3063 : 933913916 : s=s2 - s4
3064 : 933913916 : zout(1,j,nout2) = r + s
3065 : 933913916 : zout(1,j,nout4) = r - s
3066 : 933913916 : r=s1 + s3
3067 : 933913916 : s=s2 + s4
3068 : 933913916 : zout(2,j,nout1) = r + s
3069 : 933913916 : zout(2,j,nout3) = r - s
3070 : 933913916 : r=s1 - s3
3071 : 933913916 : s=r2 - r4
3072 : 933913916 : zout(2,j,nout2) = r - s
3073 : 1014332068 : zout(2,j,nout4) = r + s
3074 : : enddo
3075 : : enddo
3076 : : end if
3077 : 16589605 : 4100 continue
3078 : : end if
3079 : : else if (now.eq.8) then
3080 [ + + ]: 19455422 : if (isign.eq.-1) then
3081 : 9322472 : ia=1
3082 : 9322472 : nin1=ia-after
3083 : 9322472 : nout1=ia-atn
3084 [ + + ]: 52978320 : do ib=1,before
3085 : 43655848 : nin1=nin1+after
3086 : 43655848 : nin2=nin1+atb
3087 : 43655848 : nin3=nin2+atb
3088 : 43655848 : nin4=nin3+atb
3089 : 43655848 : nin5=nin4+atb
3090 : 43655848 : nin6=nin5+atb
3091 : 43655848 : nin7=nin6+atb
3092 : 43655848 : nin8=nin7+atb
3093 : 43655848 : nout1=nout1+atn
3094 : 43655848 : nout2=nout1+after
3095 : 43655848 : nout3=nout2+after
3096 : 43655848 : nout4=nout3+after
3097 : 43655848 : nout5=nout4+after
3098 : 43655848 : nout6=nout5+after
3099 : 43655848 : nout7=nout6+after
3100 : 43655848 : nout8=nout7+after
3101 [ + + ]: 666578480 : do j=1,n1dfft
3102 : 613600160 : r1=zin(1,j,nin1)
3103 : 613600160 : s1=zin(2,j,nin1)
3104 : 613600160 : r2=zin(1,j,nin2)
3105 : 613600160 : s2=zin(2,j,nin2)
3106 : 613600160 : r3=zin(1,j,nin3)
3107 : 613600160 : s3=zin(2,j,nin3)
3108 : 613600160 : r4=zin(1,j,nin4)
3109 : 613600160 : s4=zin(2,j,nin4)
3110 : 613600160 : r5=zin(1,j,nin5)
3111 : 613600160 : s5=zin(2,j,nin5)
3112 : 613600160 : r6=zin(1,j,nin6)
3113 : 613600160 : s6=zin(2,j,nin6)
3114 : 613600160 : r7=zin(1,j,nin7)
3115 : 613600160 : s7=zin(2,j,nin7)
3116 : 613600160 : r8=zin(1,j,nin8)
3117 : 613600160 : s8=zin(2,j,nin8)
3118 : 613600160 : r=r1 + r5
3119 : 613600160 : s=r3 + r7
3120 : 613600160 : ap=r + s
3121 : 613600160 : am=r - s
3122 : 613600160 : r=r2 + r6
3123 : 613600160 : s=r4 + r8
3124 : 613600160 : bp=r + s
3125 : 613600160 : bm=r - s
3126 : 613600160 : r=s1 + s5
3127 : 613600160 : s=s3 + s7
3128 : 613600160 : cp=r + s
3129 : 613600160 : cm=r - s
3130 : 613600160 : r=s2 + s6
3131 : 613600160 : s=s4 + s8
3132 : 613600160 : dpp=r + s
3133 : 613600160 : dm=r - s
3134 : 613600160 : zout(1,j,nout1) = ap + bp
3135 : 613600160 : zout(2,j,nout1) = cp + dpp
3136 : 613600160 : zout(1,j,nout5) = ap - bp
3137 : 613600160 : zout(2,j,nout5) = cp - dpp
3138 : 613600160 : zout(1,j,nout3) = am + dm
3139 : 613600160 : zout(2,j,nout3) = cm - bm
3140 : 613600160 : zout(1,j,nout7) = am - dm
3141 : 613600160 : zout(2,j,nout7) = cm + bm
3142 : 613600160 : r=r1 - r5
3143 : 613600160 : s=s3 - s7
3144 : 613600160 : ap=r + s
3145 : 613600160 : am=r - s
3146 : 613600160 : r=s1 - s5
3147 : 613600160 : s=r3 - r7
3148 : 613600160 : bp=r + s
3149 : 613600160 : bm=r - s
3150 : 613600160 : r=s4 - s8
3151 : 613600160 : s=r2 - r6
3152 : 613600160 : cp=r + s
3153 : 613600160 : cm=r - s
3154 : 613600160 : r=s2 - s6
3155 : 613600160 : s=r4 - r8
3156 : 613600160 : dpp=r + s
3157 : 613600160 : dm=r - s
3158 : 613600160 : r = ( cp + dm)*rt2i
3159 : 613600160 : s = ( dm - cp)*rt2i
3160 : 613600160 : cp= ( cm + dpp)*rt2i
3161 : 613600160 : dpp = ( cm - dpp)*rt2i
3162 : 613600160 : zout(1,j,nout2) = ap + r
3163 : 613600160 : zout(2,j,nout2) = bm + s
3164 : 613600160 : zout(1,j,nout6) = ap - r
3165 : 613600160 : zout(2,j,nout6) = bm - s
3166 : 613600160 : zout(1,j,nout4) = am + cp
3167 : 613600160 : zout(2,j,nout4) = bp + dpp
3168 : 613600160 : zout(1,j,nout8) = am - cp
3169 : 657256008 : zout(2,j,nout8) = bp - dpp
3170 : : enddo
3171 : : enddo
3172 [ + - ]: 9322472 : 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 : 9322472 : 8000 continue
3301 : :
3302 : : else
3303 : 10132950 : ia=1
3304 : 10132950 : nin1=ia-after
3305 : 10132950 : nout1=ia-atn
3306 [ + + ]: 58943728 : do ib=1,before
3307 : 48810778 : nin1=nin1+after
3308 : 48810778 : nin2=nin1+atb
3309 : 48810778 : nin3=nin2+atb
3310 : 48810778 : nin4=nin3+atb
3311 : 48810778 : nin5=nin4+atb
3312 : 48810778 : nin6=nin5+atb
3313 : 48810778 : nin7=nin6+atb
3314 : 48810778 : nin8=nin7+atb
3315 : 48810778 : nout1=nout1+atn
3316 : 48810778 : nout2=nout1+after
3317 : 48810778 : nout3=nout2+after
3318 : 48810778 : nout4=nout3+after
3319 : 48810778 : nout5=nout4+after
3320 : 48810778 : nout6=nout5+after
3321 : 48810778 : nout7=nout6+after
3322 : 48810778 : nout8=nout7+after
3323 [ + + ]: 737778982 : do j=1,n1dfft
3324 : 678835254 : r1=zin(1,j,nin1)
3325 : 678835254 : s1=zin(2,j,nin1)
3326 : 678835254 : r2=zin(1,j,nin2)
3327 : 678835254 : s2=zin(2,j,nin2)
3328 : 678835254 : r3=zin(1,j,nin3)
3329 : 678835254 : s3=zin(2,j,nin3)
3330 : 678835254 : r4=zin(1,j,nin4)
3331 : 678835254 : s4=zin(2,j,nin4)
3332 : 678835254 : r5=zin(1,j,nin5)
3333 : 678835254 : s5=zin(2,j,nin5)
3334 : 678835254 : r6=zin(1,j,nin6)
3335 : 678835254 : s6=zin(2,j,nin6)
3336 : 678835254 : r7=zin(1,j,nin7)
3337 : 678835254 : s7=zin(2,j,nin7)
3338 : 678835254 : r8=zin(1,j,nin8)
3339 : 678835254 : s8=zin(2,j,nin8)
3340 : 678835254 : r=r1 + r5
3341 : 678835254 : s=r3 + r7
3342 : 678835254 : ap=r + s
3343 : 678835254 : am=r - s
3344 : 678835254 : r=r2 + r6
3345 : 678835254 : s=r4 + r8
3346 : 678835254 : bp=r + s
3347 : 678835254 : bm=r - s
3348 : 678835254 : r=s1 + s5
3349 : 678835254 : s=s3 + s7
3350 : 678835254 : cp=r + s
3351 : 678835254 : cm=r - s
3352 : 678835254 : r=s2 + s6
3353 : 678835254 : s=s4 + s8
3354 : 678835254 : dpp=r + s
3355 : 678835254 : dm=r - s
3356 : 678835254 : zout(1,j,nout1) = ap + bp
3357 : 678835254 : zout(2,j,nout1) = cp + dpp
3358 : 678835254 : zout(1,j,nout5) = ap - bp
3359 : 678835254 : zout(2,j,nout5) = cp - dpp
3360 : 678835254 : zout(1,j,nout3) = am - dm
3361 : 678835254 : zout(2,j,nout3) = cm + bm
3362 : 678835254 : zout(1,j,nout7) = am + dm
3363 : 678835254 : zout(2,j,nout7) = cm - bm
3364 : 678835254 : r= r1 - r5
3365 : 678835254 : s=-s3 + s7
3366 : 678835254 : ap=r + s
3367 : 678835254 : am=r - s
3368 : 678835254 : r=s1 - s5
3369 : 678835254 : s=r7 - r3
3370 : 678835254 : bp=r + s
3371 : 678835254 : bm=r - s
3372 : 678835254 : r=-s4 + s8
3373 : 678835254 : s= r2 - r6
3374 : 678835254 : cp=r + s
3375 : 678835254 : cm=r - s
3376 : 678835254 : r=-s2 + s6
3377 : 678835254 : s= r4 - r8
3378 : 678835254 : dpp=r + s
3379 : 678835254 : dm=r - s
3380 : 678835254 : r = ( cp + dm)*rt2i
3381 : 678835254 : s = ( cp - dm)*rt2i
3382 : 678835254 : cp= ( cm + dpp)*rt2i
3383 : 678835254 : dpp= ( dpp - cm)*rt2i
3384 : 678835254 : zout(1,j,nout2) = ap + r
3385 : 678835254 : zout(2,j,nout2) = bm + s
3386 : 678835254 : zout(1,j,nout6) = ap - r
3387 : 678835254 : zout(2,j,nout6) = bm - s
3388 : 678835254 : zout(1,j,nout4) = am + cp
3389 : 678835254 : zout(2,j,nout4) = bp + dpp
3390 : 678835254 : zout(1,j,nout8) = am - cp
3391 : 727646032 : zout(2,j,nout8) = bp - dpp
3392 : : enddo
3393 : : enddo
3394 : :
3395 [ + - ]: 10132950 : 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 : 10132950 : 8001 continue
3524 : :
3525 : : end if
3526 : : else if (now.eq.3) then
3527 : : ! .5d0*sqrt(3.d0)
3528 : 50736432 : bb=isign*0.8660254037844387d0
3529 : 50736432 : ia=1
3530 : 50736432 : nin1=ia-after
3531 : 50736432 : nout1=ia-atn
3532 [ + + ]: 120271422 : do ib=1,before
3533 : 69534990 : nin1=nin1+after
3534 : 69534990 : nin2=nin1+atb
3535 : 69534990 : nin3=nin2+atb
3536 : 69534990 : nout1=nout1+atn
3537 : 69534990 : nout2=nout1+after
3538 : 69534990 : nout3=nout2+after
3539 [ + + ]: 1088650822 : do j=1,n1dfft
3540 : 968379400 : r1=zin(1,j,nin1)
3541 : 968379400 : s1=zin(2,j,nin1)
3542 : 968379400 : r2=zin(1,j,nin2)
3543 : 968379400 : s2=zin(2,j,nin2)
3544 : 968379400 : r3=zin(1,j,nin3)
3545 : 968379400 : s3=zin(2,j,nin3)
3546 : 968379400 : r=r2 + r3
3547 : 968379400 : s=s2 + s3
3548 : 968379400 : zout(1,j,nout1) = r + r1
3549 : 968379400 : zout(2,j,nout1) = s + s1
3550 : 968379400 : r1=r1 - .5d0*r
3551 : 968379400 : s1=s1 - .5d0*s
3552 : 968379400 : r2=bb*(r2-r3)
3553 : 968379400 : s2=bb*(s2-s3)
3554 : 968379400 : zout(1,j,nout2) = r1 - s2
3555 : 968379400 : zout(2,j,nout2) = s1 + r2
3556 : 968379400 : zout(1,j,nout3) = r1 + s2
3557 : 1037914390 : zout(2,j,nout3) = s1 - r2
3558 : : enddo
3559 : : enddo
3560 [ + + ]: 462635212 : do 3000,ia=2,after
3561 : 411898780 : ias=ia-1
3562 [ + + ]: 411898780 : if (4*ias.eq.3*after) then
3563 [ + + ]: 32549600 : if (isign.eq.1) then
3564 : 16618994 : nin1=ia-after
3565 : 16618994 : nout1=ia-atn
3566 [ + + ]: 39794536 : do ib=1,before
3567 : 23175542 : nin1=nin1+after
3568 : 23175542 : nin2=nin1+atb
3569 : 23175542 : nin3=nin2+atb
3570 : 23175542 : nout1=nout1+atn
3571 : 23175542 : nout2=nout1+after
3572 : 23175542 : nout3=nout2+after
3573 [ + + ]: 365634038 : do j=1,n1dfft
3574 : 325839502 : r1=zin(1,j,nin1)
3575 : 325839502 : s1=zin(2,j,nin1)
3576 : 325839502 : r2=zin(2,j,nin2)
3577 : 325839502 : s2=zin(1,j,nin2)
3578 : 325839502 : r3=zin(1,j,nin3)
3579 : 325839502 : s3=zin(2,j,nin3)
3580 : 325839502 : r=r3 + r2
3581 : 325839502 : s=s2 - s3
3582 : 325839502 : zout(1,j,nout1) = r1 - r
3583 : 325839502 : zout(2,j,nout1) = s + s1
3584 : 325839502 : r1=r1 + .5d0*r
3585 : 325839502 : s1=s1 - .5d0*s
3586 : 325839502 : r2=bb*(r2-r3)
3587 : 325839502 : s2=bb*(s2+s3)
3588 : 325839502 : zout(1,j,nout2) = r1 - s2
3589 : 325839502 : zout(2,j,nout2) = s1 - r2
3590 : 325839502 : zout(1,j,nout3) = r1 + s2
3591 : 349015044 : zout(2,j,nout3) = s1 + r2
3592 : : enddo
3593 : : enddo
3594 : : else
3595 : 15930606 : nin1=ia-after
3596 : 15930606 : nout1=ia-atn
3597 [ + + ]: 38336082 : do ib=1,before
3598 : 22405476 : nin1=nin1+after
3599 : 22405476 : nin2=nin1+atb
3600 : 22405476 : nin3=nin2+atb
3601 : 22405476 : nout1=nout1+atn
3602 : 22405476 : nout2=nout1+after
3603 : 22405476 : nout3=nout2+after
3604 [ + + ]: 360229810 : do j=1,n1dfft
3605 : 321893728 : r1=zin(1,j,nin1)
3606 : 321893728 : s1=zin(2,j,nin1)
3607 : 321893728 : r2=zin(2,j,nin2)
3608 : 321893728 : s2=zin(1,j,nin2)
3609 : 321893728 : r3=zin(1,j,nin3)
3610 : 321893728 : s3=zin(2,j,nin3)
3611 : 321893728 : r=r2 - r3
3612 : 321893728 : s=s2 + s3
3613 : 321893728 : zout(1,j,nout1) = r + r1
3614 : 321893728 : zout(2,j,nout1) = s1 - s
3615 : 321893728 : r1=r1 - .5d0*r
3616 : 321893728 : s1=s1 + .5d0*s
3617 : 321893728 : r2=bb*(r2+r3)
3618 : 321893728 : s2=bb*(s2-s3)
3619 : 321893728 : zout(1,j,nout2) = r1 + s2
3620 : 321893728 : zout(2,j,nout2) = s1 + r2
3621 : 321893728 : zout(1,j,nout3) = r1 - s2
3622 : 344299204 : zout(2,j,nout3) = s1 - r2
3623 : : enddo
3624 : : enddo
3625 : : end if
3626 [ + + ]: 379349180 : else if (8*ias.eq.3*after) then
3627 [ + + ]: 25373874 : if (isign.eq.1) then
3628 : 13323682 : nin1=ia-after
3629 : 13323682 : nout1=ia-atn
3630 [ + + ]: 31772852 : do ib=1,before
3631 : 18449170 : nin1=nin1+after
3632 : 18449170 : nin2=nin1+atb
3633 : 18449170 : nin3=nin2+atb
3634 : 18449170 : nout1=nout1+atn
3635 : 18449170 : nout2=nout1+after
3636 : 18449170 : nout3=nout2+after
3637 [ + + ]: 283814886 : do j=1,n1dfft
3638 : 252042034 : r1=zin(1,j,nin1)
3639 : 252042034 : s1=zin(2,j,nin1)
3640 : 252042034 : r=zin(1,j,nin2)
3641 : 252042034 : s=zin(2,j,nin2)
3642 : 252042034 : r2=(r - s)*rt2i
3643 : 252042034 : s2=(r + s)*rt2i
3644 : 252042034 : r3=zin(2,j,nin3)
3645 : 252042034 : s3=zin(1,j,nin3)
3646 : 252042034 : r=r2 - r3
3647 : 252042034 : s=s2 + s3
3648 : 252042034 : zout(1,j,nout1) = r + r1
3649 : 252042034 : zout(2,j,nout1) = s + s1
3650 : 252042034 : r1=r1 - .5d0*r
3651 : 252042034 : s1=s1 - .5d0*s
3652 : 252042034 : r2=bb*(r2+r3)
3653 : 252042034 : s2=bb*(s2-s3)
3654 : 252042034 : zout(1,j,nout2) = r1 - s2
3655 : 252042034 : zout(2,j,nout2) = s1 + r2
3656 : 252042034 : zout(1,j,nout3) = r1 + s2
3657 : 270491204 : zout(2,j,nout3) = s1 - r2
3658 : : enddo
3659 : : enddo
3660 : : else
3661 : 12050192 : nin1=ia-after
3662 : 12050192 : nout1=ia-atn
3663 [ + + ]: 28533304 : do ib=1,before
3664 : 16483112 : nin1=nin1+after
3665 : 16483112 : nin2=nin1+atb
3666 : 16483112 : nin3=nin2+atb
3667 : 16483112 : nout1=nout1+atn
3668 : 16483112 : nout2=nout1+after
3669 : 16483112 : nout3=nout2+after
3670 [ + + ]: 254932688 : do j=1,n1dfft
3671 : 226399384 : r1=zin(1,j,nin1)
3672 : 226399384 : s1=zin(2,j,nin1)
3673 : 226399384 : r=zin(1,j,nin2)
3674 : 226399384 : s=zin(2,j,nin2)
3675 : 226399384 : r2=(r + s)*rt2i
3676 : 226399384 : s2=(s - r)*rt2i
3677 : 226399384 : r3=zin(2,j,nin3)
3678 : 226399384 : s3=zin(1,j,nin3)
3679 : 226399384 : r=r2 + r3
3680 : 226399384 : s=s2 - s3
3681 : 226399384 : zout(1,j,nout1) = r + r1
3682 : 226399384 : zout(2,j,nout1) = s + s1
3683 : 226399384 : r1=r1 - .5d0*r
3684 : 226399384 : s1=s1 - .5d0*s
3685 : 226399384 : r2=bb*(r2-r3)
3686 : 226399384 : s2=bb*(s2+s3)
3687 : 226399384 : zout(1,j,nout2) = r1 - s2
3688 : 226399384 : zout(2,j,nout2) = s1 + r2
3689 : 226399384 : zout(1,j,nout3) = r1 + s2
3690 : 242882496 : zout(2,j,nout3) = s1 - r2
3691 : : enddo
3692 : : enddo
3693 : : end if
3694 : : else
3695 : 353975306 : itt=ias*before
3696 : 353975306 : itrig=itt+1
3697 : 353975306 : cr2=trig(1,itrig)
3698 : 353975306 : ci2=trig(2,itrig)
3699 : 353975306 : itrig=itrig+itt
3700 : 353975306 : cr3=trig(1,itrig)
3701 : 353975306 : ci3=trig(2,itrig)
3702 : 353975306 : nin1=ia-after
3703 : 353975306 : nout1=ia-atn
3704 [ + + ]: 764995528 : do ib=1,before
3705 : 411020222 : nin1=nin1+after
3706 : 411020222 : nin2=nin1+atb
3707 : 411020222 : nin3=nin2+atb
3708 : 411020222 : nout1=nout1+atn
3709 : 411020222 : nout2=nout1+after
3710 : 411020222 : nout3=nout2+after
3711 [ + + ]: 6245276626 : do j=1,n1dfft
3712 : 5480281098 : r1=zin(1,j,nin1)
3713 : 5480281098 : s1=zin(2,j,nin1)
3714 : 5480281098 : r=zin(1,j,nin2)
3715 : 5480281098 : s=zin(2,j,nin2)
3716 : 5480281098 : r2=r*cr2 - s*ci2
3717 : 5480281098 : s2=r*ci2 + s*cr2
3718 : 5480281098 : r=zin(1,j,nin3)
3719 : 5480281098 : s=zin(2,j,nin3)
3720 : 5480281098 : r3=r*cr3 - s*ci3
3721 : 5480281098 : s3=r*ci3 + s*cr3
3722 : 5480281098 : r=r2 + r3
3723 : 5480281098 : s=s2 + s3
3724 : 5480281098 : zout(1,j,nout1) = r + r1
3725 : 5480281098 : zout(2,j,nout1) = s + s1
3726 : 5480281098 : r1=r1 - .5d0*r
3727 : 5480281098 : s1=s1 - .5d0*s
3728 : 5480281098 : r2=bb*(r2-r3)
3729 : 5480281098 : s2=bb*(s2-s3)
3730 : 5480281098 : zout(1,j,nout2) = r1 - s2
3731 : 5480281098 : zout(2,j,nout2) = s1 + r2
3732 : 5480281098 : zout(1,j,nout3) = r1 + s2
3733 : 5891301320 : zout(2,j,nout3) = s1 - r2
3734 : : enddo
3735 : : enddo
3736 : : end if
3737 : 50736432 : 3000 continue
3738 : : else if (now==5) then
3739 : : ! cos(2.d0*pi/5.d0)
3740 : 18921133 : cos2=0.3090169943749474d0
3741 : : ! cos(4.d0*pi/5.d0)
3742 : 18921133 : cos4=-0.8090169943749474d0
3743 : : ! sin(2.d0*pi/5.d0)
3744 : 18921133 : sin2=isign*0.9510565162951536d0
3745 : : ! sin(4.d0*pi/5.d0)
3746 : 18921133 : sin4=isign*0.5877852522924731d0
3747 : 18921133 : ia=1
3748 : 18921133 : nin1=ia-after
3749 : 18921133 : nout1=ia-atn
3750 [ + + ]: 152188684 : do ib=1,before
3751 : 133267551 : nin1=nin1+after
3752 : 133267551 : nin2=nin1+atb
3753 : 133267551 : nin3=nin2+atb
3754 : 133267551 : nin4=nin3+atb
3755 : 133267551 : nin5=nin4+atb
3756 : 133267551 : nout1=nout1+atn
3757 : 133267551 : nout2=nout1+after
3758 : 133267551 : nout3=nout2+after
3759 : 133267551 : nout4=nout3+after
3760 : 133267551 : nout5=nout4+after
3761 [ + + ]: 1623901152 : do j=1,n1dfft
3762 : 1471712468 : r1=zin(1,j,nin1)
3763 : 1471712468 : s1=zin(2,j,nin1)
3764 : 1471712468 : r2=zin(1,j,nin2)
3765 : 1471712468 : s2=zin(2,j,nin2)
3766 : 1471712468 : r3=zin(1,j,nin3)
3767 : 1471712468 : s3=zin(2,j,nin3)
3768 : 1471712468 : r4=zin(1,j,nin4)
3769 : 1471712468 : s4=zin(2,j,nin4)
3770 : 1471712468 : r5=zin(1,j,nin5)
3771 : 1471712468 : s5=zin(2,j,nin5)
3772 : 1471712468 : r25 = r2 + r5
3773 : 1471712468 : r34 = r3 + r4
3774 : 1471712468 : s25 = s2 - s5
3775 : 1471712468 : s34 = s3 - s4
3776 : 1471712468 : zout(1,j,nout1) = r1 + r25 + r34
3777 : 1471712468 : r = r1 + cos2*r25 + cos4*r34
3778 : 1471712468 : s = sin2*s25 + sin4*s34
3779 : 1471712468 : zout(1,j,nout2) = r - s
3780 : 1471712468 : zout(1,j,nout5) = r + s
3781 : 1471712468 : r = r1 + cos4*r25 + cos2*r34
3782 : 1471712468 : s = sin4*s25 - sin2*s34
3783 : 1471712468 : zout(1,j,nout3) = r - s
3784 : 1471712468 : zout(1,j,nout4) = r + s
3785 : 1471712468 : r25 = r2 - r5
3786 : 1471712468 : r34 = r3 - r4
3787 : 1471712468 : s25 = s2 + s5
3788 : 1471712468 : s34 = s3 + s4
3789 : 1471712468 : zout(2,j,nout1) = s1 + s25 + s34
3790 : 1471712468 : r = s1 + cos2*s25 + cos4*s34
3791 : 1471712468 : s = sin2*r25 + sin4*r34
3792 : 1471712468 : zout(2,j,nout2) = r + s
3793 : 1471712468 : zout(2,j,nout5) = r - s
3794 : 1471712468 : r = s1 + cos4*s25 + cos2*s34
3795 : 1471712468 : s = sin4*r25 - sin2*r34
3796 : 1471712468 : zout(2,j,nout3) = r + s
3797 : 1604980019 : zout(2,j,nout4) = r - s
3798 : : enddo
3799 : : enddo
3800 [ + + ]: 41495105 : do 5000,ia=2,after
3801 : 22573972 : ias=ia-1
3802 [ + + ]: 22573972 : 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 : 22213780 : ias=ia-1
3920 : 22213780 : itt=ias*before
3921 : 22213780 : itrig=itt+1
3922 : 22213780 : cr2=trig(1,itrig)
3923 : 22213780 : ci2=trig(2,itrig)
3924 : 22213780 : itrig=itrig+itt
3925 : 22213780 : cr3=trig(1,itrig)
3926 : 22213780 : ci3=trig(2,itrig)
3927 : 22213780 : itrig=itrig+itt
3928 : 22213780 : cr4=trig(1,itrig)
3929 : 22213780 : ci4=trig(2,itrig)
3930 : 22213780 : itrig=itrig+itt
3931 : 22213780 : cr5=trig(1,itrig)
3932 : 22213780 : ci5=trig(2,itrig)
3933 : 22213780 : nin1=ia-after
3934 : 22213780 : nout1=ia-atn
3935 [ + + ]: 94210152 : do ib=1,before
3936 : 71996372 : nin1=nin1+after
3937 : 71996372 : nin2=nin1+atb
3938 : 71996372 : nin3=nin2+atb
3939 : 71996372 : nin4=nin3+atb
3940 : 71996372 : nin5=nin4+atb
3941 : 71996372 : nout1=nout1+atn
3942 : 71996372 : nout2=nout1+after
3943 : 71996372 : nout3=nout2+after
3944 : 71996372 : nout4=nout3+after
3945 : 71996372 : nout5=nout4+after
3946 [ + + ]: 848473019 : do j=1,n1dfft
3947 : 754262867 : r1=zin(1,j,nin1)
3948 : 754262867 : s1=zin(2,j,nin1)
3949 : 754262867 : r=zin(1,j,nin2)
3950 : 754262867 : s=zin(2,j,nin2)
3951 : 754262867 : r2=r*cr2 - s*ci2
3952 : 754262867 : s2=r*ci2 + s*cr2
3953 : 754262867 : r=zin(1,j,nin3)
3954 : 754262867 : s=zin(2,j,nin3)
3955 : 754262867 : r3=r*cr3 - s*ci3
3956 : 754262867 : s3=r*ci3 + s*cr3
3957 : 754262867 : r=zin(1,j,nin4)
3958 : 754262867 : s=zin(2,j,nin4)
3959 : 754262867 : r4=r*cr4 - s*ci4
3960 : 754262867 : s4=r*ci4 + s*cr4
3961 : 754262867 : r=zin(1,j,nin5)
3962 : 754262867 : s=zin(2,j,nin5)
3963 : 754262867 : r5=r*cr5 - s*ci5
3964 : 754262867 : s5=r*ci5 + s*cr5
3965 : 754262867 : r25 = r2 + r5
3966 : 754262867 : r34 = r3 + r4
3967 : 754262867 : s25 = s2 - s5
3968 : 754262867 : s34 = s3 - s4
3969 : 754262867 : zout(1,j,nout1) = r1 + r25 + r34
3970 : 754262867 : r = r1 + cos2*r25 + cos4*r34
3971 : 754262867 : s = sin2*s25 + sin4*s34
3972 : 754262867 : zout(1,j,nout2) = r - s
3973 : 754262867 : zout(1,j,nout5) = r + s
3974 : 754262867 : r = r1 + cos4*r25 + cos2*r34
3975 : 754262867 : s = sin4*s25 - sin2*s34
3976 : 754262867 : zout(1,j,nout3) = r - s
3977 : 754262867 : zout(1,j,nout4) = r + s
3978 : 754262867 : r25 = r2 - r5
3979 : 754262867 : r34 = r3 - r4
3980 : 754262867 : s25 = s2 + s5
3981 : 754262867 : s34 = s3 + s4
3982 : 754262867 : zout(2,j,nout1) = s1 + s25 + s34
3983 : 754262867 : r = s1 + cos2*s25 + cos4*s34
3984 : 754262867 : s = sin2*r25 + sin4*r34
3985 : 754262867 : zout(2,j,nout2) = r + s
3986 : 754262867 : zout(2,j,nout5) = r - s
3987 : 754262867 : r = s1 + cos4*s25 + cos2*s34
3988 : 754262867 : s = sin4*r25 - sin2*r34
3989 : 754262867 : zout(2,j,nout3) = r + s
3990 : 826259239 : zout(2,j,nout4) = r - s
3991 : : enddo
3992 : : enddo
3993 : : end if
3994 : 18921133 : 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 : 139623218 : end subroutine fftstp
4074 : : !!***
4075 : :
4076 : : end module m_sg2002
4077 : : !!***
|
