Line data Source code
1 : !Local variables-------------------------------
2 : !scalars
3 : integer,parameter :: nt1=1
4 : integer :: g2max_in,g2min_in,g2max_out,g2min_out,ix,iy,iz,ig,igb,mgb,nx
5 : integer :: ny,nz,nfft,ngbin,ngbout,nthreads,ldxy,ng2_pos,ng2_neg,pt,status
6 : integer :: lot,nlot,lotin,lotout,fftcache,cidx,sidx,xyplane,padx
7 : integer :: nxhalf1,nxhalfm,nyhalf1,cidx_inv,pidx,fidx,padz,padc
8 : integer :: igb_inv,iiy,ixmax,iz_inv, G_TO_R, R_TO_G, iscale__, G_TO_R_SCALE, R_TO_G_SCALE
9 : integer :: arr(2)
10 : real(dp) :: fraction,arg
11 : complex(MYKIND) :: vloc,phase
12 : character(len=500) :: msg
13 : type(DFTI_DESCRIPTOR),pointer :: bw_plan,fw_plan
14 : type(DFTI_DESCRIPTOR),pointer :: fw_plan_xg2pos,fw_plan_xg2neg,fw_plan_y
15 : type(DFTI_DESCRIPTOR),pointer :: bw_plan_xg2pos,bw_plan_xg2neg,bw_plan_y
16 : type(C_PTR) :: cptr_cuboid,cptr_sa,cptr_sb,cptr_sc,cptr_sd
17 : !arrays
18 62924146 : integer,allocatable :: indpw_kin(:,:),indpw_kout(:,:)
19 62924146 : complex(MYKIND),allocatable :: pha1(:),pha2(:),pha3(:)
20 62924146 : complex(MYKIND),ABI_CONTIGUOUS pointer :: cuboid(:),slice_a(:),slice_b(:),slice_c(:),slice_d(:)
21 : ! *************************************************************************
22 :
23 : ! gcc does not like C-pointers in OMP clauses
24 : #if defined HAVE_OPENMP && defined __GFORTRAN__
25 : #define DEV_USE_ABIMALLOC
26 : #endif
27 :
28 62924146 : G_TO_R = +1; R_TO_G = -1
29 62924146 : R_TO_G_SCALE = DFTI_FORWARD_SCALE; G_TO_R_SCALE = DFTI_BACKWARD_SCALE
30 62924146 : if (present(abi_convention)) then
31 213778 : if (.not. abi_convention) then
32 0 : G_TO_R = -1; R_TO_G = +1
33 0 : R_TO_G_SCALE = DFTI_BACKWARD_SCALE; G_TO_R_SCALE = DFTI_FORWARD_SCALE
34 : end if
35 : end if
36 62924146 : iscale__ = 1; if (present(iscale)) iscale__ = iscale
37 :
38 62924146 : if (istwf_k>2 .and. option==0) then
39 0 : write(msg,'(a,i0)')' option=0 is not allowed with istwf_k=',istwf_k
40 0 : ABI_BUG(msg)
41 : end if
42 :
43 62924146 : if (istwf_k>=2 .and. option==3) then
44 0 : write(msg,'(a,i0)')' option=3 is not allowed with istwf_k=',istwf_k
45 0 : ABI_BUG(msg)
46 : end if
47 :
48 : !For all other tests of validity of inputs, assume that they
49 : !have been done in the calling routine
50 :
51 62924146 : nx=ngfft(1); ny=ngfft(2); nz=ngfft(3); nfft=nx*ny*nz
52 62924146 : fftcache=ngfft(8)
53 62924146 : ldxy = ldx*ldy
54 :
55 62924146 : if (option/=3) then
56 186351630 : ABI_MALLOC(indpw_kin,(4,npwin))
57 62117210 : call indfftrisc(gboundin(3:3+2*mgfft+4,1),indpw_kin,kg_kin,mgfft,ngbin,ngfft,npwin)
58 : end if
59 :
60 62924146 : if (option==2 .or. option==3) then
61 157481895 : ABI_MALLOC(indpw_kout,(4,npwout))
62 52493965 : call indfftrisc(gboundout(3:3+2*mgfft+4,1),indpw_kout,kg_kout,mgfft,ngbout,ngfft,npwout)
63 : end if
64 :
65 : ! Define the dimension of the first work arrays, for 1D transforms along z ,
66 : ! taking into account the need to avoid the cache trashing
67 62924146 : if (option==2) then
68 51687029 : mgb=MAX(ngbin,ngbout)
69 11237117 : else if (option==0 .or. option==1) then
70 10430181 : mgb=ngbin; ngbout=1
71 806936 : else if (option==3) then
72 806936 : mgb=ngbout; ngbin=1
73 : end if
74 :
75 62924146 : if (MOD(mgb,2)/=1) mgb=mgb+1
76 :
77 : !Initialise openmp, if needed
78 62924146 : nthreads = xomp_get_num_threads(open_parallel=.TRUE.)
79 :
80 : !For the treatment of the z transform, one tries to use only a fraction of the cache, since the
81 : !treatment of the array cuboid will not involve contiguous segments
82 62924146 : fraction=0.25
83 : !First estimation of lot and nlot
84 62924146 : lot=(fftcache*fraction*1000)/(nz*8*2)+1
85 : !Select the smallest integer multiple of nthreads, larger
86 : !or equal to nlot. In this way, the cache size is not exhausted,
87 : !and one takes care correctly of the number of processors.
88 : !Treat separately the in and out cases
89 62924146 : nlot=(ngbin-1)/lot+1
90 62924146 : nlot=nthreads*((nlot-1)/nthreads+1)
91 62924146 : lotin=(ngbin-1)/nlot+1
92 62924146 : nlot=(ngbout-1)/lot+1
93 62924146 : nlot=nthreads*((nlot-1)/nthreads+1)
94 62924146 : lotout=(ngbout-1)/nlot+1
95 : !The next line impose only one lot. Usually, comment it.
96 : !lotin=mgb; lotout=mgb
97 : !write(std_out,*)"lotin, lotout",lotin,lotout
98 :
99 : ! ===============================
100 : ! ===== General k-point code ====
101 : ! ===============================
102 :
103 62924146 : if (istwf_k==1) then
104 :
105 61236675 : call dfti_alloc_complex(mgb*nz,cptr_cuboid,cuboid)
106 :
107 61236675 : if (option/=3) then
108 :
109 : !$OMP PARALLEL DO
110 56076538038 : do cidx=1,mgb*nz
111 56076538038 : cuboid(cidx) = MYCZERO
112 : end do
113 :
114 : ! Insert fofgin into the work array
115 : !$OMP PARALLEL DO PRIVATE(igb,iz,cidx)
116 15730067058 : do ig=1,npwin
117 15669637319 : igb =indpw_kin(4,ig)
118 15669637319 : iz =indpw_kin(3,ig)
119 15669637319 : cidx = igb + (iz-1)*mgb
120 15730067058 : cuboid(cidx) = CMPLX(fofgin(1,ig),fofgin(2,ig), KIND=MYKIND)
121 : end do
122 :
123 : #ifndef HAVE_OPENMP
124 : ! ngbin 1D in-place transforms of cuboid(Gx,Gy,Gz) along Gz.
125 60429739 : status = DftiCreateDescriptor(bw_plan, FFT_PRECISION, DFTI_COMPLEX, 1, nz)
126 60429739 : status = DftiSetValue(bw_plan, DFTI_NUMBER_OF_TRANSFORMS, ngbin)
127 60429739 : arr(1) = 0
128 60429739 : arr(2) = mgb
129 60429739 : status = DftiSetValue(bw_plan, DFTI_INPUT_STRIDES, arr)
130 60429739 : status = DftiSetValue(bw_plan, DFTI_INPUT_DISTANCE, 1)
131 60429739 : status = DftiSetValue(bw_plan, DFTI_NUMBER_OF_USER_THREADS, nt1)
132 60429739 : DFTI_CHECK(status)
133 :
134 60429739 : status = DftiCommitDescriptor(bw_plan)
135 60429739 : DFTI_CHECK(status)
136 :
137 : !status = DftiComputeBackward(bw_plan, cuboid)
138 : !DFTI_CHECK(status)
139 60429739 : DFTI_CHECK(compute_ip(bw_plan, G_TO_R, cuboid))
140 :
141 : #else
142 : ! Distribute 1D in-place transforms of cuboid(Gx,Gy,Gz) along Gz among the threads.
143 : status = DftiCreateDescriptor(bw_plan, FFT_PRECISION, DFTI_COMPLEX, 1, nz)
144 : arr(1) = 0
145 : arr(2) = mgb
146 : status = DftiSetValue(bw_plan, DFTI_INPUT_STRIDES, arr)
147 : status = DftiSetValue(bw_plan, DFTI_NUMBER_OF_USER_THREADS, nthreads)
148 :
149 : status = DftiCommitDescriptor(bw_plan)
150 : DFTI_CHECK(status)
151 :
152 : !!$OMP PARALLEL DO PRIVATE(status) SCHEDULE(STATIC, lotin)
153 : !$OMP PARALLEL DO PRIVATE(status)
154 : do igb=1,ngbin
155 : !status = DftiComputeBackward(bw_plan, cuboid(igb:))
156 : !DFTI_CHECK(status)
157 : DFTI_CHECK(compute_ip(bw_plan, G_TO_R, cuboid(igb:)))
158 : end do
159 : #endif
160 :
161 60429739 : status = DftiFreeDescriptor(bw_plan)
162 60429739 : DFTI_CHECK(status)
163 : ! Now we have cuboid(Gx,Gy,z) in the cuboid enclosing the G-sphere
164 : end if
165 : !
166 : ! ================================
167 : ! ==== Generate the FFT plans ====
168 : ! ================================
169 : !
170 : !(1:g2max_in+1,n3) ! Positive g_y.
171 : !(g2min_in+ny+1:ny,n3) ! Negative g_y.
172 61236675 : g2min_in = gboundin(3,1)
173 61236675 : g2max_in = gboundin(4,1)
174 :
175 61236675 : g2min_out = gboundout(3,1)
176 61236675 : g2max_out = gboundout(4,1)
177 :
178 : ! ng2_pos 1D transforms of f(Gx,Gy,z) along Gx for Gy >= 0
179 61236675 : ng2_pos = g2max_in + 1
180 61236675 : status = DftiCreateDescriptor(bw_plan_xg2pos , FFT_PRECISION, DFTI_COMPLEX, 1, nx)
181 61236675 : DFTI_CHECK(status)
182 :
183 61236675 : status = DftiSetValue(bw_plan_xg2pos, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
184 61236675 : status = DftiSetValue(bw_plan_xg2pos, DFTI_NUMBER_OF_TRANSFORMS, ng2_pos)
185 61236675 : status = DftiSetValue(bw_plan_xg2pos, DFTI_INPUT_DISTANCE, ldx)
186 61236675 : arr(1) = 0
187 61236675 : arr(2) = 1
188 61236675 : status = DftiSetValue(bw_plan_xg2pos, DFTI_INPUT_STRIDES, arr)
189 61236675 : status = DftiSetValue(bw_plan_xg2pos, DFTI_OUTPUT_DISTANCE, ldx)
190 61236675 : status = DftiSetValue(bw_plan_xg2pos, DFTI_OUTPUT_STRIDES, arr)
191 61236675 : status = DftiSetValue(bw_plan_xg2pos, DFTI_NUMBER_OF_USER_THREADS, nthreads)
192 :
193 61236675 : status = DftiCommitDescriptor(bw_plan_xg2pos)
194 61236675 : DFTI_CHECK(status)
195 :
196 : ! ng2_neg 1D transforms of f(Gx,Gy,z) along Gx for Gy < 0
197 61236675 : ng2_neg = -g2min_in
198 61236675 : status = DftiCreateDescriptor(bw_plan_xg2neg , FFT_PRECISION, DFTI_COMPLEX, 1, nx)
199 61236675 : DFTI_CHECK(status)
200 :
201 61236675 : status = DftiSetValue(bw_plan_xg2neg, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
202 61236675 : status = DftiSetValue(bw_plan_xg2neg, DFTI_NUMBER_OF_TRANSFORMS, ng2_neg)
203 61236675 : status = DftiSetValue(bw_plan_xg2neg, DFTI_INPUT_DISTANCE, ldx)
204 : arr(1) = 0
205 : arr(2) = 1
206 61236675 : status = DftiSetValue(bw_plan_xg2neg, DFTI_INPUT_STRIDES, arr)
207 61236675 : status = DftiSetValue(bw_plan_xg2neg, DFTI_OUTPUT_DISTANCE, ldx)
208 61236675 : status = DftiSetValue(bw_plan_xg2neg, DFTI_OUTPUT_STRIDES, arr)
209 61236675 : status = DftiSetValue(bw_plan_xg2neg, DFTI_NUMBER_OF_USER_THREADS, nthreads)
210 :
211 61236675 : status = DftiCommitDescriptor(bw_plan_xg2neg)
212 61236675 : DFTI_CHECK(status)
213 :
214 : ! nx 1D transforms of f(x,Gy,z) along Gy.
215 61236675 : status = DftiCreateDescriptor(bw_plan_y, FFT_PRECISION, DFTI_COMPLEX, 1, ny)
216 61236675 : DFTI_CHECK(status)
217 :
218 61236675 : status = DftiSetValue(bw_plan_y, DFTI_NUMBER_OF_TRANSFORMS, nx)
219 61236675 : status = DftiSetValue(bw_plan_y, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
220 61236675 : status = DftiSetValue(bw_plan_y, DFTI_INPUT_DISTANCE, 1)
221 : arr(1) = 0
222 61236675 : arr(2) = ldx
223 61236675 : status = DftiSetValue(bw_plan_y, DFTI_INPUT_STRIDES, arr)
224 61236675 : status = DftiSetValue(bw_plan_y, DFTI_OUTPUT_DISTANCE, 1)
225 61236675 : status = DftiSetValue(bw_plan_y, DFTI_OUTPUT_STRIDES, arr)
226 61236675 : status = DftiSetValue(bw_plan_y, DFTI_NUMBER_OF_USER_THREADS, nthreads)
227 61236675 : DFTI_CHECK(status)
228 :
229 61236675 : status = DftiCommitDescriptor(bw_plan_y)
230 61236675 : DFTI_CHECK(status)
231 :
232 61236675 : if (option==2 .or. option==3) then
233 : ! nx 1D transforms of f(x,y,Gz) along y.
234 50942016 : status = DftiCreateDescriptor(fw_plan_y , FFT_PRECISION, DFTI_COMPLEX, 1, ny)
235 50942016 : DFTI_CHECK(status)
236 :
237 50942016 : status = DftiSetValue(fw_plan_y, DFTI_NUMBER_OF_TRANSFORMS, nx)
238 50942016 : status = DftiSetValue(fw_plan_y, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
239 50942016 : status = DftiSetValue(fw_plan_y, DFTI_INPUT_DISTANCE, 1)
240 : arr(1) = 0
241 : arr(2) = ldx
242 50942016 : status = DftiSetValue(fw_plan_y, DFTI_INPUT_STRIDES, arr)
243 50942016 : status = DftiSetValue(fw_plan_y, DFTI_OUTPUT_DISTANCE, 1)
244 50942016 : status = DftiSetValue(fw_plan_y, DFTI_OUTPUT_STRIDES, arr)
245 50942016 : status = DftiSetValue(fw_plan_y, DFTI_NUMBER_OF_USER_THREADS, nthreads)
246 50942016 : DFTI_CHECK(status)
247 :
248 50942016 : status = DftiCommitDescriptor(fw_plan_y)
249 50942016 : DFTI_CHECK(status)
250 :
251 : ! ng2_pos 1D transforms of f(x,y,Gz) along x.
252 50942016 : ng2_pos = g2max_out + 1
253 50942016 : status = DftiCreateDescriptor(fw_plan_xg2pos , FFT_PRECISION, DFTI_COMPLEX, 1, nx)
254 50942016 : DFTI_CHECK(status)
255 :
256 50942016 : status = DftiSetValue(fw_plan_xg2pos, DFTI_NUMBER_OF_TRANSFORMS, ng2_pos)
257 50942016 : status = DftiSetValue(fw_plan_xg2pos, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
258 50942016 : status = DftiSetValue(fw_plan_xg2pos, DFTI_INPUT_DISTANCE, ldx)
259 : !status = DftiSetValue(fw_plan_xg2pos, DFTI_INPUT_STRIDES, arr)
260 50942016 : status = DftiSetValue(fw_plan_xg2pos, DFTI_OUTPUT_DISTANCE, ldx)
261 50942016 : status = DftiSetValue(fw_plan_xg2pos, DFTI_NUMBER_OF_USER_THREADS, nthreads)
262 50942016 : DFTI_CHECK(status)
263 :
264 50942016 : status = DftiCommitDescriptor(fw_plan_xg2pos)
265 50942016 : DFTI_CHECK(status)
266 :
267 : ! ng2_neg 1D transforms of f(x,y,Gz) along x.
268 50942016 : ng2_neg = -g2min_out
269 50942016 : status = DftiCreateDescriptor(fw_plan_xg2neg , FFT_PRECISION, DFTI_COMPLEX, 1, nx)
270 50942016 : DFTI_CHECK(status)
271 :
272 50942016 : status = DftiSetValue(fw_plan_xg2neg, DFTI_NUMBER_OF_TRANSFORMS, ng2_neg)
273 50942016 : status = DftiSetValue(fw_plan_xg2neg, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
274 50942016 : status = DftiSetValue(fw_plan_xg2neg, DFTI_INPUT_DISTANCE, ldx)
275 : !status = DftiSetValue(fw_plan_xg2neg, DFTI_INPUT_STRIDES, arr)
276 50942016 : status = DftiSetValue(fw_plan_xg2neg, DFTI_OUTPUT_DISTANCE, ldx)
277 50942016 : status = DftiSetValue(fw_plan_xg2neg, DFTI_NUMBER_OF_USER_THREADS, nthreads)
278 :
279 50942016 : status = DftiCommitDescriptor(fw_plan_xg2neg)
280 50942016 : DFTI_CHECK(status)
281 : end if
282 :
283 : ! Open OMP parallel region and allocate two 2-dimensional work arrays for out-of-place transforms.
284 : #ifdef DEV_USE_ABIMALLOC
285 : !$OMP PARALLEL PRIVATE(status,ix,iy,cidx,sidx,padx,xyplane,pt,vloc,slice_a,slice_b)
286 :
287 : ABI_MALLOC(slice_a,(ldxy))
288 : ABI_MALLOC(slice_b,(ldxy))
289 :
290 : #else
291 : !$OMP PARALLEL PRIVATE(status,ix,iy,cidx,sidx,padx,xyplane,pt,vloc,slice_a,slice_b,cptr_sa,cptr_sb)
292 : ! This causes a sigfault in gcc44 due to the presence of C-pointers in the OMP statement.
293 :
294 61236675 : call dfti_alloc_complex(ldxy,cptr_sa,slice_a)
295 61236675 : call dfti_alloc_complex(ldxy,cptr_sb,slice_b)
296 : #endif
297 :
298 : ! Big Loop over z parallelized with OpenMP.
299 : !$OMP DO
300 1083699948 : do iz=1,nz
301 :
302 1022463273 : if (option/=3) then
303 : ! Zero the values on the current plane
304 >36740*10^7 : slice_a(:) = MYCZERO
305 : ! Copy the data in the current plane
306 56460882560 : do igb=1,ngbin
307 55452349299 : ix = indpw_kin(1,igb)
308 55452349299 : iy = indpw_kin(2,igb)
309 55452349299 : cidx = igb + (iz-1)*mgb
310 55452349299 : sidx = ix + (iy-1)*ldx
311 56460882560 : slice_a(sidx) = cuboid(cidx)
312 : end do
313 : ! Perform Gx transform, taking into account arrays of zeros
314 1008533261 : if (g2min_in+ny >= g2max_in+2) then
315 10134787743 : do iy=g2max_in+2,g2min_in+ny
316 9126254482 : padx = (iy-1)*ldx
317 >18862*10^7 : do ix=1,nx
318 >17848*10^7 : sidx = ix + padx
319 >18761*10^7 : slice_b(sidx) = MYCZERO
320 : end do
321 : end do
322 : end if
323 :
324 : !status = DftiComputeBackward(bw_plan_xg2pos, slice_a, slice_b)
325 : !DFTI_CHECK(status)
326 1008533261 : DFTI_CHECK(compute_op(bw_plan_xg2pos, G_TO_R, slice_a, slice_b))
327 :
328 1008533261 : pt = 1 + (g2min_in+ny)*ldx
329 : !status = DftiComputeBackward(bw_plan_xg2neg, slice_a(pt:), slice_b(pt:))
330 : !DFTI_CHECK(status)
331 1008533261 : DFTI_CHECK(compute_op(bw_plan_xg2neg, G_TO_R, slice_a(pt:), slice_b(pt:)))
332 :
333 :
334 : ! Got f(x,Gy,z). Now perform y transform
335 : !status = DftiComputeBackward(bw_plan_y, slice_b, slice_a)
336 : !DFTI_CHECK(status)
337 1008533261 : DFTI_CHECK(compute_op(bw_plan_y, G_TO_R, slice_b, slice_a))
338 : end if
339 : ! The wave function is now in real space, for the current z plane
340 :
341 129443879 : SELECT CASE (option)
342 : CASE (0)
343 : ! Copy the transformed function at the right place and we are done!
344 129443879 : xyplane = 1 + (iz-1)*ldx*ldy
345 : #ifndef HAVE_DFTI_MIXED_PRECISION
346 500 : if (MYKIND==dp) then
347 129443379 : call ZCOPY(ldxy,slice_a,1,fofr(1,xyplane),1)
348 : else if (MYKIND==sp) then
349 500 : call CCOPY(ldxy,slice_a,1,fofr(1,xyplane),1)
350 : else
351 : ABI_ERROR("Wrong FFT precision")
352 : end if
353 : #else
354 0 : do ix=0,ldxy-1
355 0 : fofr(1, xyplane + ix) = real(slice_a(ix+1), kind=dp)
356 0 : fofr(2, xyplane + ix) = aimag(slice_a(ix+1))
357 : end do
358 : #endif
359 :
360 : CASE (1)
361 : ! Accumulate density
362 : !THIS IS FOR TESTING PURPOSE
363 38925136 : if (abs(weight_r-weight_i)<tol12) then
364 751615782 : do iy=1,ny
365 712776814 : padx = (iy-1)*ldx
366 15761030051 : do ix=1,nx
367 15009414269 : sidx = ix + padx
368 15722191083 : denpot(ix,iy,iz)=denpot(ix,iy,iz)+weight_r*(REAL(slice_a(sidx))**2+AIMAG(slice_a(sidx))**2)
369 : end do
370 : end do
371 : else
372 1871160 : do iy=1,ny
373 1784992 : padx = (iy-1)*ldx
374 39231608 : do ix=1,nx
375 37360448 : sidx = ix + padx
376 39145440 : denpot(ix,iy,iz)=denpot(ix,iy,iz)+weight_r*(REAL(slice_a(sidx))**2)+weight_i*(AIMAG(slice_a(sidx))**2)
377 : end do
378 : end do
379 : end if
380 :
381 : CASE (2)
382 : ! Apply local potential
383 840164246 : if (cplex==1) then
384 13185840247 : do iy=1,ny
385 12473542973 : padx = (iy-1)*ldx
386 >26523*10^7 : do ix=1,nx
387 >25204*10^7 : sidx = ix + padx
388 >26451*10^7 : slice_a(sidx) = denpot(ix,iy,iz)*slice_a(sidx)
389 : end do
390 : end do
391 : else
392 2249387868 : do iy=1,ny
393 2121520896 : padx = (iy-1)*ldx
394 40395678516 : do ix=1,nx
395 38146290648 : sidx = ix + padx
396 38146290648 : vloc = CMPLX(denpot(2*ix-1,iy,iz), denpot(2*ix,iy,iz), KIND=MYKIND)
397 40267811544 : slice_a(sidx) = vloc * slice_a(sidx)
398 : end do
399 : end do
400 : end if
401 :
402 : CASE (3)
403 : ! Copy the function to be transformed at the right place
404 13930012 : xyplane = 1 + (iz-1)*ldx*ldy
405 : #ifndef HAVE_DFTI_MIXED_PRECISION
406 1018423571 : if (MYKIND==dp) then
407 13929512 : call ZCOPY(ldxy,fofr(1,xyplane),1,slice_a,1)
408 : else if (MYKIND==sp) then
409 500 : call CCOPY(ldxy,fofr(1,xyplane),1,slice_a,1)
410 : else
411 : ABI_ERROR("Wrong FFT precision")
412 : end if
413 : #else
414 4039702 : do ix=1,ldxy
415 0 : slice_a(ix) = cmplx(fofr(1, xyplane + ix -1), fofr(2, xyplane + ix -1), kind=MYKIND)
416 : end do
417 : #endif
418 : END SELECT
419 :
420 1083699948 : if (option==2 .or. option==3) then
421 : ! Go back to G space.
422 : !
423 : ! 1) Perform y transform
424 : !status = DftiComputeForward(fw_plan_y, slice_a, slice_b)
425 : !DFTI_CHECK(status)
426 854094258 : DFTI_CHECK(compute_op(fw_plan_y, R_TO_G, slice_a, slice_b))
427 :
428 : ! Perform x transform, taking into account arrays of zeros
429 : !status = DftiComputeForward(fw_plan_xg2pos, slice_b, slice_a)
430 : !DFTI_CHECK(status)
431 854094258 : DFTI_CHECK(compute_op(fw_plan_xg2pos, R_TO_G, slice_b, slice_a))
432 :
433 854094258 : pt = 1 + (g2min_out+ny)*ldx
434 : !status = DftiComputeForward(fw_plan_xg2neg, slice_b(pt:), slice_a(pt:))
435 : !DFTI_CHECK(status)
436 854094258 : DFTI_CHECK(compute_op(fw_plan_xg2neg, R_TO_G, slice_b(pt:), slice_a(pt:)))
437 :
438 : ! Copy the data from the current plane to cuboid
439 49663948161 : do igb=1,ngbout
440 48809853903 : ix=indpw_kout(1,igb)
441 48809853903 : iy=indpw_kout(2,igb)
442 48809853903 : sidx = ix + (iy-1)*ldx
443 48809853903 : cidx = igb + (iz-1)*mgb
444 49663948161 : cuboid(cidx)=slice_a(sidx)
445 : end do
446 : end if
447 :
448 : end do ! End loop on planes iz
449 : !$OMP END DO
450 :
451 : #ifdef DEV_USE_ABIMALLOC
452 : ABI_FREE(slice_a)
453 : ABI_FREE(slice_b)
454 : #else
455 61236675 : call dfti_free(cptr_sa)
456 61236675 : call dfti_free(cptr_sb)
457 : #endif
458 : !$OMP END PARALLEL
459 :
460 : ! Free plans
461 61236675 : status = DftiFreeDescriptor(bw_plan_xg2pos)
462 61236675 : status = DftiFreeDescriptor(bw_plan_xg2neg)
463 61236675 : status = DftiFreeDescriptor(bw_plan_y)
464 :
465 61236675 : if (option==2 .or. option==3)then
466 50942016 : status = DftiFreeDescriptor(fw_plan_y)
467 50942016 : status = DftiFreeDescriptor(fw_plan_xg2pos)
468 50942016 : status = DftiFreeDescriptor(fw_plan_xg2neg)
469 : end if
470 :
471 61236675 : if (option==2 .or. option==3) then ! Do 1D FFTs on the z direction
472 :
473 : #ifndef HAVE_OPENMP
474 : ! ngbout 1D transforms of f(x,y,Gz) along Gz.
475 50942016 : status = DftiCreateDescriptor(fw_plan, FFT_PRECISION, DFTI_COMPLEX, 1, nz)
476 50942016 : DFTI_CHECK(status)
477 :
478 50942016 : status = DftiSetValue(fw_plan, DFTI_NUMBER_OF_TRANSFORMS, ngbout)
479 : arr(1) = 0
480 50942016 : arr(2) = mgb
481 50942016 : status = DftiSetValue(fw_plan, DFTI_INPUT_STRIDES, arr)
482 50942016 : status = DftiSetValue(fw_plan, DFTI_INPUT_DISTANCE, 1)
483 50942016 : status = DftiSetValue(fw_plan, DFTI_NUMBER_OF_USER_THREADS, nt1)
484 50942016 : status = DftiSetValue(fw_plan, R_TO_G_SCALE, one/DBLE(nfft)) ! normalization
485 :
486 50942016 : status = DftiCommitDescriptor(fw_plan)
487 50942016 : DFTI_CHECK(status)
488 :
489 : !status = DftiComputeForward(fw_plan, cuboid)
490 : !DFTI_CHECK(status)
491 50942016 : DFTI_CHECK(compute_ip(fw_plan, R_TO_G, cuboid))
492 :
493 : #else
494 : ! Distribute ngbout 1D transforms of f(x,y,Gz) along Gz among the threads.
495 : status = DftiCreateDescriptor(fw_plan, FFT_PRECISION, DFTI_COMPLEX, 1, nz)
496 : DFTI_CHECK(status)
497 :
498 : arr(1) = 0
499 : arr(2) = mgb
500 : status = DftiSetValue(fw_plan, DFTI_INPUT_STRIDES, arr)
501 : status = DftiSetValue(fw_plan, DFTI_NUMBER_OF_USER_THREADS, nthreads)
502 : status = DftiSetValue(fw_plan, R_TO_G_SCALE, one/DBLE(nfft)) ! normalization
503 : DFTI_CHECK(status)
504 :
505 : status = DftiCommitDescriptor(fw_plan)
506 : DFTI_CHECK(status)
507 :
508 : !!$OMP PARALLEL DO PRIVATE(status) SCHEDULE(STATIC, lotout)
509 : !$OMP PARALLEL DO PRIVATE(status)
510 : do igb=1,ngbout
511 : !status = DftiComputeForward(fw_plan, cuboid(igb:))
512 : !DFTI_CHECK(status)
513 : DFTI_CHECK(compute_ip(fw_plan, R_TO_G, cuboid(igb:)))
514 : end do
515 : #endif
516 :
517 50942016 : status = DftiFreeDescriptor(fw_plan)
518 :
519 : ! Transfer the data in the output array, after normalization
520 : !$OMP PARALLEL DO PRIVATE(igb,iz,cidx)
521 13990606331 : do ig=1,npwout
522 13939664315 : igb = indpw_kout(4,ig)
523 13939664315 : iz = indpw_kout(3,ig)
524 13939664315 : cidx = igb + (iz-1)*mgb
525 13939664315 : fofgout(1,ig) = REAL (cuboid(cidx))
526 13990606331 : fofgout(2,ig) = AIMAG(cuboid(cidx))
527 : end do
528 : end if
529 :
530 61236675 : call dfti_free(cptr_cuboid)
531 :
532 1687471 : else if (istwf_k>=2) then
533 :
534 : !------------------------------------------------------------------
535 : !Here, use of time-reversal symmetry
536 1687471 : nxhalf1 = nx/2+1
537 1687471 : nxhalfm = (nx+1)/2
538 1687471 : nyhalf1 = ny/2+1
539 :
540 1687471 : call dfti_alloc_complex(mgb*nz,cptr_cuboid,cuboid)
541 :
542 1687471 : if (istwf_k/=2) then
543 : ! Precompute phases.
544 3196569 : ABI_MALLOC(pha1,(nx))
545 3196569 : ABI_MALLOC(pha2,(ny))
546 3196569 : ABI_MALLOC(pha3,(nz))
547 24379326 : do ix=1,nx
548 23313803 : arg = DBLE(ix-1)*pi/DBLE(nx)
549 24379326 : pha1(ix)= MYCMPLX(COS(arg), SIN(arg))
550 : end do
551 24367064 : do iy=1,ny
552 23301541 : arg = DBLE(iy-1)*pi/DBLE(ny)
553 24367064 : pha2(iy)= MYCMPLX(COS(arg), SIN(arg))
554 : end do
555 25744408 : do iz=1,nz
556 24678885 : arg = DBLE(iz-1)*pi/DBLE(nz)
557 25744408 : pha3(iz)= MYCMPLX(COS(arg), SIN(arg))
558 : end do
559 : end if
560 :
561 1687471 : if (option/=3) then
562 : ! Zero the components of cuboid
563 : !$OMP PARALLEL DO
564 3530938720 : do cidx=1,mgb*nz
565 3530938720 : cuboid(cidx) = MYCZERO
566 : end do
567 :
568 : ! Insert fofgin into the work array cuboid
569 : !$OMP PARALLEL DO PRIVATE(igb,iz,cidx)
570 876041864 : do ig=1,npwin
571 874354393 : igb = indpw_kin(4,ig)
572 874354393 : iz = indpw_kin(3,ig)
573 874354393 : cidx = igb + (iz-1)*mgb
574 876041864 : cuboid(cidx) = MYCMPLX(fofgin(1,ig), fofgin(2,ig))
575 : end do
576 :
577 : ! Must complete the iy=1 plane when $k_y \equiv 0$
578 :
579 : ! Take care of ix=1 when $k_x \equiv 0$
580 1687471 : if (istwf_k==2) then
581 9956307 : do iz=nz/2+1,nz
582 9334359 : iz_inv = nz+2-iz
583 9334359 : cidx = 1 + (iz-1)*mgb
584 9334359 : cidx_inv = 1 + (iz_inv-1)*mgb
585 9956307 : cuboid(cidx) = MYCONJG(cuboid(cidx_inv))
586 : end do
587 1065523 : else if (istwf_k==4)then
588 1583996 : do iz=nz/2+1,nz
589 1467249 : iz_inv = nz+1-iz
590 1467249 : cidx = 1 + (iz-1)*mgb
591 1467249 : cidx_inv = 1 + (iz_inv-1)*mgb
592 1583996 : cuboid(cidx) = MYCONJG(cuboid(cidx_inv))
593 : end do
594 : end if
595 :
596 : ! Now, take care of other ix values, except iz==1 when $k_z \equiv 0$
597 1687471 : ixmax=gboundin(6,1)+1
598 :
599 1687471 : if (istwf_k==2) then
600 : !$OMP PARALLEL DO PRIVATE(igb_inv,iz_inv,cidx,cidx_inv)
601 7320544 : do igb=2,2*ixmax-1
602 6698596 : igb_inv=2*ixmax+1-igb
603 125537148 : do iz=nz/2+1,nz
604 118216604 : iz_inv = nz+2-iz
605 118216604 : cidx = igb + (iz-1)*mgb
606 118216604 : cidx_inv = igb_inv + (iz_inv-1)*mgb
607 124915200 : cuboid(cidx) = MYCONJG(cuboid(cidx_inv))
608 : end do
609 : end do
610 :
611 1065523 : else if (istwf_k==3) then
612 : !$OMP PARALLEL DO PRIVATE(igb_inv,iz_inv,cidx,cidx_inv)
613 3245351 : do igb=1,2*ixmax
614 2952034 : igb_inv=2*ixmax+1-igb
615 46698145 : do iz=nz/2+1,nz
616 43452794 : iz_inv = nz+2-iz
617 43452794 : cidx = igb + (iz-1)*mgb
618 43452794 : cidx_inv = igb_inv + (iz_inv-1)*mgb
619 46404828 : cuboid(cidx) = MYCONJG(cuboid(cidx_inv))
620 : end do
621 : end do
622 :
623 772206 : else if (istwf_k==4) then
624 : !$OMP PARALLEL DO PRIVATE(igb_inv,iz_inv,cidx,cidx_inv)
625 949329 : do igb=2,2*ixmax-1
626 832582 : igb_inv=2*ixmax+1-igb
627 11862713 : do iz=nz/2+1,nz
628 10913384 : iz_inv = nz+1-iz
629 10913384 : cidx = igb + (iz-1)*mgb
630 10913384 : cidx_inv = igb_inv + (iz_inv-1)*mgb
631 11745966 : cuboid(cidx) = MYCONJG(cuboid(cidx_inv))
632 : end do
633 : end do
634 :
635 655459 : else if (istwf_k==5) then
636 : !$OMP PARALLEL DO PRIVATE(igb_inv,iz_inv,cidx,cidx_inv)
637 901514 : do igb=1,2*ixmax
638 803102 : igb_inv=2*ixmax+1-igb
639 9561170 : do iz=nz/2+1,nz
640 8659656 : iz_inv = nz+1-iz
641 8659656 : cidx = igb + (iz-1)*mgb
642 8659656 : cidx_inv = igb_inv + (iz_inv-1)*mgb
643 9462758 : cuboid(cidx) = MYCONJG(cuboid(cidx_inv))
644 : end do
645 : end do
646 :
647 : end if
648 :
649 : ! Now, iz==1
650 1687471 : if (istwf_k==2) then
651 3971246 : do igb=2,ixmax
652 3349298 : igb_inv = 2*ixmax+1-igb
653 3349298 : cidx_inv = igb_inv !+ (iz_inv-1)*mgb
654 3971246 : cuboid(cidx_inv) = MYCONJG(cuboid(igb))
655 : end do
656 1065523 : else if (istwf_k==3) then
657 1769334 : do igb=1,ixmax
658 1476017 : igb_inv = 2*ixmax+1-igb
659 1476017 : cidx_inv = igb_inv !+ (iz_inv-1)*mgb
660 1769334 : cuboid(cidx_inv) = MYCONJG(cuboid(igb))
661 : end do
662 : end if
663 : !
664 : ! Perform in-place 1D FFTs on the z direction
665 : ! ngbin 1D in-place transforms of cuboid(Gx,Gy,Gz) along Gz.
666 : ! bw_plan = fftw3_plan_many_dft(1, (/nz/), ngbin, &
667 : !& cuboid, (/mgb, nz/), mgb, 1, &
668 : !& cuboid, (/mgb, nz/), mgb, 1, ABI_FFTW_BACKWARD, ABI_FFTW_ESTIMATE, nthreads)
669 :
670 1687471 : status = DftiCreateDescriptor(bw_plan, FFT_PRECISION, DFTI_COMPLEX, 1, nz)
671 1687471 : status = DftiSetValue(bw_plan, DFTI_NUMBER_OF_TRANSFORMS, ngbin)
672 1687471 : arr(1) = 0
673 1687471 : arr(2) = mgb
674 1687471 : status = DftiSetValue(bw_plan, DFTI_INPUT_STRIDES, arr)
675 1687471 : status = DftiSetValue(bw_plan, DFTI_INPUT_DISTANCE, 1)
676 1687471 : status = DftiSetValue(bw_plan, DFTI_NUMBER_OF_USER_THREADS, nt1)
677 1687471 : DFTI_CHECK(status)
678 :
679 1687471 : status = DftiCommitDescriptor(bw_plan)
680 1687471 : DFTI_CHECK(status)
681 :
682 : !status = DftiComputeBackward(bw_plan, cuboid)
683 : !DFTI_CHECK(status)
684 1687471 : DFTI_CHECK(compute_ip(bw_plan, G_TO_R, cuboid))
685 :
686 1687471 : status = DftiFreeDescriptor(bw_plan)
687 1687471 : DFTI_CHECK(status)
688 :
689 : !call fftw3_execute_dft(bw_plan, cuboid, cuboid)
690 : !call fftw3_destroy_plan(bw_plan)
691 :
692 : ! Change the phase if $k_z \neq 0$
693 1687471 : if (istwf_k==4 .or. istwf_k==5 .or. istwf_k==8 .or. istwf_k==9) then
694 : !$OMP PARALLEL DO PRIVATE(phase,padz,cidx)
695 9970694 : do iz=1,nz
696 9513022 : phase = pha3(iz)
697 9513022 : padz = (iz-1)*mgb
698 397295815 : do igb=1,ngbin
699 387325121 : cidx = igb + padz
700 396838143 : cuboid(cidx) = cuboid(cidx) * phase
701 : end do
702 : end do
703 : end if
704 :
705 : end if ! if(option/=3)
706 :
707 1687471 : g2max_in= gboundin(4,1)
708 1687471 : ng2_pos = g2max_in + 1
709 :
710 : ! bw_plan_xg2pos = fftw3_plan_many_dft(1, (/nx/), ng2_pos, & ! ng2_pos 1D transforms of f(Gx,Gy,z) along Gx
711 : !& slice_a, (/ldx, ldy/), 1, ldx, & ! for Gy >= 0
712 : !& slice_b, (/ldx, ldy/), 1, ldx, ABI_FFTW_BACKWARD, ABI_FFTW_ESTIMATE, nt1)
713 :
714 1687471 : status = DftiCreateDescriptor(bw_plan_xg2pos , FFT_PRECISION, DFTI_COMPLEX, 1, nx)
715 1687471 : DFTI_CHECK(status)
716 :
717 1687471 : status = DftiSetValue(bw_plan_xg2pos, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
718 1687471 : status = DftiSetValue(bw_plan_xg2pos, DFTI_NUMBER_OF_TRANSFORMS, ng2_pos)
719 1687471 : status = DftiSetValue(bw_plan_xg2pos, DFTI_INPUT_DISTANCE, ldx)
720 1687471 : arr(1) = 0
721 1687471 : arr(2) = 1
722 1687471 : status = DftiSetValue(bw_plan_xg2pos, DFTI_INPUT_STRIDES, arr)
723 1687471 : status = DftiSetValue(bw_plan_xg2pos, DFTI_OUTPUT_DISTANCE, ldx)
724 1687471 : status = DftiSetValue(bw_plan_xg2pos, DFTI_OUTPUT_STRIDES, arr)
725 1687471 : status = DftiSetValue(bw_plan_xg2pos, DFTI_NUMBER_OF_USER_THREADS, nthreads)
726 :
727 1687471 : status = DftiCommitDescriptor(bw_plan_xg2pos)
728 1687471 : DFTI_CHECK(status)
729 :
730 : ! TODO clarify this point. why nxhalfm?
731 : ! bw_plan_y = fftw3_plan_many_dft(1, (/ny/), nxhalf1, & ! nx 1D transforms of f(x,Gy,z) along Gy.
732 : !& slice_a, (/ldx, ldy/), ldx, 1, &
733 : !& slice_b, (/ldx, ldy/), ldx, 1, ABI_FFTW_BACKWARD, ABI_FFTW_ESTIMATE, nt1)
734 :
735 1687471 : status = DftiCreateDescriptor(bw_plan_y, FFT_PRECISION, DFTI_COMPLEX, 1, ny)
736 1687471 : DFTI_CHECK(status)
737 :
738 1687471 : status = DftiSetValue(bw_plan_y, DFTI_NUMBER_OF_TRANSFORMS, nxhalf1)
739 1687471 : status = DftiSetValue(bw_plan_y, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
740 1687471 : status = DftiSetValue(bw_plan_y, DFTI_INPUT_DISTANCE, 1)
741 : arr(1) = 0
742 1687471 : arr(2) = ldx
743 1687471 : status = DftiSetValue(bw_plan_y, DFTI_INPUT_STRIDES, arr)
744 1687471 : status = DftiSetValue(bw_plan_y, DFTI_OUTPUT_DISTANCE, 1)
745 1687471 : status = DftiSetValue(bw_plan_y, DFTI_OUTPUT_STRIDES, arr)
746 1687471 : status = DftiSetValue(bw_plan_y, DFTI_NUMBER_OF_USER_THREADS, nthreads)
747 1687471 : DFTI_CHECK(status)
748 :
749 1687471 : status = DftiCommitDescriptor(bw_plan_y)
750 1687471 : DFTI_CHECK(status)
751 :
752 1687471 : if (option==2 .or. option==3)then
753 : ! TODO clarify this point. why nxhalfm?
754 : ! fw_plan_y = fftw3_plan_many_dft(1, (/ny/), nxhalf1, & ! nx 1D transforms of f(x,y,Gz) along y.
755 : !& slice_a, (/ldx, ldy/), ldx, 1, &
756 : !& slice_b, (/ldx, ldy/), ldx, 1, ABI_FFTW_FORWARD, ABI_FFTW_ESTIMATE, nt1)
757 :
758 1551949 : status = DftiCreateDescriptor(fw_plan_y , FFT_PRECISION, DFTI_COMPLEX, 1, ny)
759 1551949 : DFTI_CHECK(status)
760 :
761 1551949 : status = DftiSetValue(fw_plan_y, DFTI_NUMBER_OF_TRANSFORMS, nxhalf1)
762 1551949 : status = DftiSetValue(fw_plan_y, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
763 1551949 : status = DftiSetValue(fw_plan_y, DFTI_INPUT_DISTANCE, 1)
764 : arr(1) = 0
765 : arr(2) = ldx
766 1551949 : status = DftiSetValue(fw_plan_y, DFTI_INPUT_STRIDES, arr)
767 1551949 : status = DftiSetValue(fw_plan_y, DFTI_OUTPUT_DISTANCE, 1)
768 1551949 : status = DftiSetValue(fw_plan_y, DFTI_OUTPUT_STRIDES, arr)
769 1551949 : status = DftiSetValue(fw_plan_y, DFTI_NUMBER_OF_USER_THREADS, nthreads)
770 1551949 : DFTI_CHECK(status)
771 :
772 1551949 : status = DftiCommitDescriptor(fw_plan_y)
773 1551949 : DFTI_CHECK(status)
774 :
775 : ! Note that here gboundin must equal gboudnout
776 : ! TODO clarify this point. why nxhalfm?
777 : ng2_pos = g2max_in + 1
778 : ! fw_plan_xg2pos = fftw3_plan_many_dft(1, (/nx/), ng2_pos, & ! ng2_pos 1D transforms of f(x,y,Gz) along x.
779 : !& slice_c, (/2*nxhalfm, ldy/), 1, 2*nxhalfm, &
780 : !& slice_d, (/2*nxhalfm, ldy/), 1, 2*nxhalfm, ABI_FFTW_FORWARD, ABI_FFTW_ESTIMATE, nt1)
781 :
782 1551949 : status = DftiCreateDescriptor(fw_plan_xg2pos , FFT_PRECISION, DFTI_COMPLEX, 1, nx)
783 1551949 : DFTI_CHECK(status)
784 :
785 1551949 : status = DftiSetValue(fw_plan_xg2pos, DFTI_NUMBER_OF_TRANSFORMS, ng2_pos)
786 1551949 : status = DftiSetValue(fw_plan_xg2pos, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
787 1551949 : status = DftiSetValue(fw_plan_xg2pos, DFTI_INPUT_DISTANCE, 2*nxhalfm)
788 : !arr(1) = 0
789 : !arr(2) = ldx
790 : !status = DftiSetValue(fw_plan_xg2pos, DFTI_INPUT_STRIDES, arr)
791 1551949 : status = DftiSetValue(fw_plan_xg2pos, DFTI_OUTPUT_DISTANCE, 2*nxhalfm)
792 1551949 : status = DftiSetValue(fw_plan_xg2pos, DFTI_NUMBER_OF_USER_THREADS, nthreads)
793 1551949 : DFTI_CHECK(status)
794 :
795 1551949 : status = DftiCommitDescriptor(fw_plan_xg2pos)
796 1551949 : DFTI_CHECK(status)
797 : end if
798 :
799 : ! Do-loop on the planes stacked in the z direction
800 : #ifdef DEV_USE_ABIMALLOC
801 : !$OMP PARALLEL PRIVATE(ix,iy,iiy,cidx,sidx,pidx,fidx,padx,padz,phase,slice_a,slice_b,slice_c,slice_d)
802 : ! Allocate two 2-dimensional work arrays
803 : ABI_MALLOC(slice_a,(ldxy))
804 : ABI_MALLOC(slice_b,(ldxy))
805 :
806 : ABI_MALLOC(slice_c,(2*nxhalfm*ldy))
807 : ABI_MALLOC(slice_d,(2*nxhalfm*ldy))
808 :
809 : #else
810 : !$OMP PARALLEL &
811 : !$OMP& PRIVATE(ix,iy,iiy,cidx,sidx,pidx,fidx,padx,padz,phase,slice_a,slice_b,slice_c,slice_d,cptr_sa,cptr_sb,cptr_sc,cptr_sd)
812 1687471 : call dfti_alloc_complex(ldxy,cptr_sa,slice_a)
813 1687471 : call dfti_alloc_complex(ldxy,cptr_sb,slice_b)
814 :
815 1687471 : call dfti_alloc_complex(2*nxhalfm*ldy,cptr_sc,slice_c)
816 1687471 : call dfti_alloc_complex(2*nxhalfm*ldy,cptr_sd,slice_d)
817 : #endif
818 :
819 : !$OMP DO
820 44970446 : do iz=1,nz
821 :
822 43282975 : if (option/=3) then
823 : !
824 : ! Zero the values on the current plane: needed only from iy=1 to g2max_in+1
825 : !do iy=1,g2max_in+1
826 : ! padx = (iy-1)*ldx
827 : ! do ix=1,nx
828 : ! sidx = ix + padx
829 : ! slice_a(sidx) = MYCZERO
830 : ! end do
831 : !end do
832 43307500470 : slice_a = MYCZERO
833 : !
834 : ! Copy the data in the current plane
835 3543164644 : do igb=1,ngbin
836 3499881669 : ix = indpw_kin(1,igb)
837 3499881669 : iy = indpw_kin(2,igb)
838 3499881669 : cidx = igb + (iz-1)*mgb
839 3499881669 : sidx = ix + (iy-1)*ldx
840 3543164644 : slice_a(sidx) = cuboid(cidx)
841 : end do
842 :
843 : ! Perform x transform, taking into account arrays of zeros
844 : !call sg_fftx(fftcache,mfac,mg,ldx,ldy,1,g2max_in+1,1,slice_a,slice_b,trig1,aft1,now1,bef1,one,ind1,ic1)
845 : !call fftw3_execute_dft(bw_plan_xg2pos, slice_a, slice_b)
846 :
847 : !status = DftiComputeBackward(bw_plan_xg2pos, slice_a, slice_b)
848 : !DFTI_CHECK(status)
849 43282975 : DFTI_CHECK(compute_op(bw_plan_xg2pos, G_TO_R, slice_a, slice_b))
850 :
851 : ! Change the phase if $k_x \neq 0$
852 43282975 : if (istwf_k==3 .or. istwf_k==5 .or. istwf_k==7 .or. istwf_k==9) then
853 139083650 : do iy=1,g2max_in+1
854 120893094 : padx = (iy-1)*ldx
855 3947189153 : do ix=1,nx
856 3808105503 : phase = pha1(ix)
857 3808105503 : sidx = ix + padx
858 3928998597 : slice_b(sidx) = slice_b(sidx) * phase
859 : end do
860 : end do
861 : end if
862 :
863 : ! Compute symmetric and antisymmetric combinations
864 43282975 : if (istwf_k>=2 .and. istwf_k<=5) then
865 477829226 : do ix=1,nxhalf1-1
866 477829226 : slice_a(ix) = MYCMPLX(REAL(slice_b(2*ix-1)), REAL(slice_b(2*ix)))
867 : end do
868 : ! If nx odd, must add last data
869 30848370 : if ((2*nxhalf1-2)/=nx) then
870 2831283 : slice_a(nxhalf1)= REAL(slice_b(nx))
871 : end if
872 : iiy=2
873 : else
874 : iiy=1
875 : end if
876 :
877 43282975 : if (g2max_in+1 >= iiy) then
878 318355531 : do iy=iiy,g2max_in+1
879 : !
880 275072556 : padx = (iy-1)*ldx
881 4869804446 : do ix=1,nxhalf1-1
882 4594731890 : sidx = ix + padx
883 4594731890 : cidx = 2*ix-1 + padx
884 :
885 : slice_a(sidx) = MYCMPLX( REAL(slice_b(cidx)) - AIMAG(slice_b(cidx+1)),&
886 4594731890 : & AIMAG(slice_b(cidx)) + REAL (slice_b(cidx+1)) )
887 :
888 4594731890 : sidx = ix + (ny+iiy-iy-1)*ldx
889 : slice_a(sidx) = MYCMPLX( REAL(slice_b(cidx)) + AIMAG(slice_b(cidx+1)),&
890 4869804446 : & -AIMAG(slice_b(cidx)) + REAL (slice_b(cidx+1)) )
891 : end do
892 : !
893 318355531 : if ((2*nxhalf1-2)/=nx) then
894 23145978 : sidx = nxhalf1 + padx
895 23145978 : cidx = nx + padx
896 23145978 : slice_a(sidx) = slice_b(cidx)
897 :
898 23145978 : sidx = nxhalf1 + (ny+iiy-iy-1)*ldx
899 23145978 : slice_a(sidx) = MYCONJG(slice_b(cidx))
900 : end if
901 : end do
902 : end if
903 :
904 43282975 : if (nyhalf1 >= g2max_in+2) then
905 396019649 : do iy=g2max_in+2,nyhalf1
906 6110654956 : do ix=1,nxhalf1-1
907 5757918282 : sidx = ix + (iy-1)*ldx
908 5757918282 : slice_a(sidx) = MYCZERO
909 :
910 5757918282 : sidx = ix + (ny+iiy-iy-1)*ldx
911 6110654956 : slice_a(sidx) =zero
912 : end do
913 396019649 : if ((2*nxhalf1-2)/=nx) then
914 28317935 : sidx = nxhalf1 + (iy-1)*ldx
915 28317935 : slice_a(sidx) = MYCZERO
916 :
917 28317935 : sidx = nxhalf1 + (ny+iiy-iy-1)*ldx
918 28317935 : slice_a(sidx) = MYCZERO
919 : end if
920 : end do
921 : end if
922 :
923 : !nxi=1
924 : !call sg_ffty(fftcache,mfac,mg,ldx,ldy,1,nxi,nxhalfm,1,1,slice_a,slice_b,trig2,aft2,now2,bef2,one,ind2,ic2)
925 : !call fftw3_execute_dft(bw_plan_y, slice_a, slice_b)
926 :
927 : !status = DftiComputeBackward(bw_plan_y, slice_a, slice_b)
928 : !DFTI_CHECK(status)
929 43282975 : DFTI_CHECK(compute_op(bw_plan_y, G_TO_R, slice_a, slice_b))
930 :
931 : ! Change the phase if $k_y \neq 0$
932 43282975 : if (istwf_k>=6 .and. istwf_k<=9) then
933 326738482 : do iy=1,ny
934 314303877 : phase = pha2(iy)
935 314303877 : padx = (iy-1)*ldx
936 4937328637 : do ix=1,nxhalfm
937 4610590155 : sidx = ix + padx
938 4924894032 : slice_b(sidx) = slice_b(sidx) * phase
939 : end do
940 : end do
941 : end if
942 :
943 : end if ! option/=3
944 :
945 : ! The wave function is now in real space, for the current plane,
946 : ! represented by REAL numbers, although packed in the complex array slice_b
947 :
948 43282975 : if (option==0) then
949 : ! This option is only permitted for istwf_k==2 (Gamma point)
950 : ! Copy the transformed function at the right place
951 615810 : do iy=1,ny
952 596697 : padx = (iy-1)*ldx
953 14578134 : do ix=1,nxhalf1-1
954 13981437 : sidx = ix + (iy-1)*ldx
955 13981437 : fidx = 2*ix-1 + padx + (iz-1)*ldx*ldy
956 :
957 13981437 : fofr(1,fidx) = REAL(slice_b(sidx))
958 13981437 : fofr(2,fidx) = zero
959 :
960 13981437 : fofr(1,fidx+1) = AIMAG(slice_b(sidx))
961 14578134 : fofr(2,fidx+1) = zero
962 : end do
963 : ! If nx odd, must add last data
964 615810 : if ((2*nxhalf1-2)/=nx) then
965 88821 : sidx = nxhalf1 + (iy-1)*ldx
966 88821 : fidx = nx + (iy-1) * ldx + (iz-1)*ldx*ldy
967 88821 : fofr(1,fidx) = REAL(slice_b(sidx))
968 88821 : fofr(2,fidx) = zero
969 : end if
970 : end do
971 : end if
972 :
973 43282975 : if (option==1) then ! Accumulate density
974 118104138 : do iy=1,ny
975 114271702 : padx = (iy-1)*ldx
976 2050839229 : do ix=1,nxhalf1-1
977 1936567527 : sidx = ix + padx
978 1936567527 : denpot(2*ix-1,iy,iz) = denpot(2*ix-1,iy,iz) + weight_r* REAL(slice_b(sidx))**2
979 2050839229 : denpot(2*ix ,iy,iz) = denpot(2*ix ,iy,iz) + weight_i*AIMAG(slice_b(sidx))**2
980 : end do
981 : ! If nx odd, must add last data
982 118104138 : if ((2*nxhalf1-2)/=nx) then
983 10865538 : sidx = nxhalf1 + padx
984 10865538 : denpot(nx,iy,iz) = denpot(nx,iy,iz)+weight_r*REAL(slice_b(sidx))**2
985 : end if
986 : end do
987 : end if
988 :
989 43282975 : if (option==2) then ! Apply local potential
990 1159301227 : do iy=1,ny
991 1119869801 : padx = (iy-1)*ldx
992 19621296663 : do ix=1,nxhalf1-1
993 18501426862 : sidx = ix + padx
994 : slice_a(sidx) = MYCMPLX( denpot(2*ix-1,iy,iz)* REAL(slice_b(sidx)), &
995 19621296663 : denpot(2*ix ,iy,iz)*AIMAG(slice_b(sidx)) )
996 : end do
997 : ! If nx odd, must add last data
998 1159301227 : if ((2*nxhalf1-2)/=nx) then
999 93680901 : sidx = nxhalf1 + padx
1000 93680901 : slice_a(sidx) = MYCMPLX( denpot(nx,iy,iz)*REAL(slice_b(sidx)), zero)
1001 : end if
1002 : end do
1003 : end if
1004 :
1005 43282975 : if (option==3) then
1006 : ! This option is only permitted for istwf_k==2 (Gamma point)
1007 : ! Copy the transformed function at the right place
1008 0 : do iy=1,ny
1009 0 : padx = (iy-1)*ldx
1010 0 : do ix=1,nxhalf1-1
1011 0 : sidx = ix + padx
1012 0 : fidx = 2*ix-1 + padx + (iz-1)*ldx*ldy
1013 0 : slice_b(sidx) = MYCMPLX( fofr(1,fidx), fofr(1,fidx+1) )
1014 : end do
1015 : ! If nx odd, must add last data
1016 0 : if ((2*nxhalf1-2)/=nx) then
1017 0 : sidx = nxhalf1 + padx
1018 0 : fidx = nx + padx + (iz-1)*ldx*ldy
1019 0 : slice_b(sidx) = MYCMPLX( fofr(1,fidx), zero)
1020 : end if
1021 : end do
1022 : end if
1023 :
1024 44970446 : if (option==2 .or. option==3) then
1025 : !
1026 : ! Change the phase if $k_y \neq 0$
1027 39431426 : if (istwf_k>=6 .and. istwf_k<=9) then
1028 300327672 : do iy=1,ny
1029 288813361 : padx = (iy-1)*ldx
1030 288813361 : phase = MYCONJG(pha2(iy))
1031 4505154443 : do ix=1,nxhalfm
1032 4204826771 : sidx = ix + padx
1033 4493640132 : slice_a(sidx) = slice_a(sidx) * phase
1034 : end do
1035 : end do
1036 : end if
1037 : !
1038 : ! Perform y transform
1039 : !nxi=1
1040 : !call sg_ffty(fftcache,mfac,mg,ldx,ldy,1,nxi,nxhalfm,1,1,slice_a,slice_b,trig5,aft5,now5,bef5,-one,ind5,ic5)
1041 : !call fftw3_execute_dft(fw_plan_y, slice_a, slice_b)
1042 :
1043 : !status = DftiComputeForward(fw_plan_y, slice_a, slice_b)
1044 : !DFTI_CHECK(status)
1045 39431426 : DFTI_CHECK(compute_op(fw_plan_y, R_TO_G, slice_a, slice_b))
1046 :
1047 : ! Decompose symmetric and antisymmetric parts
1048 39431426 : if (istwf_k>=2 .and. istwf_k<=5) then
1049 434430834 : do ix=1,nxhalfm
1050 406513719 : pidx = 2*ix-1
1051 406513719 : slice_c(pidx) = MYCMPLX( REAL(slice_b(ix)), zero)
1052 434430834 : slice_c(pidx+1) = MYCMPLX( AIMAG(slice_b(ix)), zero)
1053 : end do
1054 : iiy=2
1055 : else
1056 : iiy=1
1057 : end if
1058 :
1059 289009612 : do iy=iiy,g2max_in+1
1060 249578186 : padx = (iy-1)*ldx
1061 4467824080 : do ix=1,nxhalfm
1062 4178814468 : sidx = ix + padx
1063 4178814468 : cidx = ix + (ny+iiy-iy-1)*ldx
1064 4178814468 : pidx = 2*ix-1 + (iy-1)*2*nxhalfm
1065 :
1066 : slice_c(pidx) = MYCMPLX( &
1067 : & ( REAL(slice_b(sidx)) + REAL(slice_b(cidx)))*0.5d0, &
1068 4178814468 : & (AIMAG(slice_b(sidx)) - AIMAG(slice_b(cidx)))*0.5d0 )
1069 :
1070 : slice_c(pidx+1) = MYCMPLX( &
1071 : & (AIMAG(slice_b(sidx)) + AIMAG(slice_b(cidx)))*0.5d0, &
1072 4428392654 : & (-REAL(slice_b(sidx)) + REAL(slice_b(cidx)))*0.5d0 )
1073 : end do
1074 : end do
1075 : !
1076 : ! Change the phase if $k_x \neq 0$
1077 39431426 : if (istwf_k==3 .or. istwf_k==5 .or. istwf_k==7 .or. istwf_k==9) then
1078 127497730 : do iy=1,g2max_in+1
1079 3591861077 : do ix=1,nx
1080 3464363347 : phase = MYCONJG(pha1(ix))
1081 3464363347 : pidx = ix + (iy-1)*2*nxhalfm
1082 3575091231 : slice_c(pidx)= slice_c(pidx) * phase
1083 : end do
1084 : end do
1085 : end if
1086 : !
1087 : ! Perform x transform: for y=1 to g2max_in+1, to benefit from zeros
1088 : !call sg_fftx(fftcache,mfac,mg,2*nxhalfm,ldy,1,g2max_in+1,1,slice_c,slice_d,trig4,aft4,now4,bef4,-one,ind4,ic4)
1089 : !call fftw3_execute_dft(fw_plan_xg2pos, slice_c, slice_d)
1090 :
1091 : !status = DftiComputeForward(fw_plan_xg2pos, slice_c, slice_d)
1092 : !DFTI_CHECK(status)
1093 39431426 : DFTI_CHECK(compute_op(fw_plan_xg2pos, R_TO_G, slice_c, slice_d))
1094 :
1095 : ! Copy the data from the current plane to cuboid.
1096 3209393273 : do igb=1,ngbout
1097 3169961847 : ix=indpw_kout(1,igb)
1098 3169961847 : iy=indpw_kout(2,igb)
1099 3169961847 : cidx = igb + (iz-1)*mgb
1100 3169961847 : pidx = ix + (iy-1)*2*nxhalfm
1101 3209393273 : cuboid(cidx) = slice_d(pidx)
1102 : end do
1103 :
1104 : end if ! option==2 or 3
1105 :
1106 : end do ! End loop on planes
1107 : !$OMP END DO
1108 :
1109 : #ifdef DEV_USE_ABIMALLOC
1110 : ABI_FREE(slice_a)
1111 : ABI_FREE(slice_b)
1112 : ABI_FREE(slice_c)
1113 : ABI_FREE(slice_d)
1114 : #else
1115 1687471 : call dfti_free(cptr_sa)
1116 1687471 : call dfti_free(cptr_sb)
1117 1687471 : call dfti_free(cptr_sc)
1118 1687471 : call dfti_free(cptr_sd)
1119 : #endif
1120 : !$OMP END PARALLEL
1121 :
1122 1687471 : status = DftiFreeDescriptor(bw_plan_xg2pos)
1123 1687471 : status = DftiFreeDescriptor(bw_plan_y)
1124 :
1125 : !call fftw3_destroy_plan(bw_plan_xg2pos)
1126 : !call fftw3_destroy_plan(bw_plan_y)
1127 :
1128 1687471 : if (option==2 .or. option==3) then
1129 :
1130 1551949 : status = DftiFreeDescriptor(fw_plan_y)
1131 1551949 : status = DftiFreeDescriptor(fw_plan_xg2pos)
1132 :
1133 : !call fftw3_destroy_plan(fw_plan_y)
1134 : !call fftw3_destroy_plan(fw_plan_xg2pos)
1135 :
1136 : ! Change the phase if $k_z \neq 0$
1137 1551949 : if (istwf_k==4 .or. istwf_k==5 .or. istwf_k==8 .or. istwf_k==9) then
1138 : !$OMP PARALLEL DO PRIVATE(phase,cidx,padc)
1139 9206929 : do iz=1,nz
1140 8779104 : phase = MYCONJG(pha3(iz))
1141 8779104 : padc = (iz-1)*mgb
1142 365526778 : do igb=1,ngbout
1143 356319849 : cidx = igb + padc
1144 365098953 : cuboid(cidx) = cuboid(cidx) * phase
1145 : end do
1146 : end do
1147 : end if
1148 :
1149 : ! 1D FFTs on the z direction
1150 : ! fw_plan = fftw3_plan_many_dft(1, (/nz/), ngbout, & ! ngbout 1D transforms of f(x,y,Gz) along Gz.
1151 : !& cuboid, (/mgb, nz, 1/), mgb, 1, &
1152 : !& cuboid, (/mgb, nz, 1/), mgb, 1, ABI_FFTW_FORWARD, ABI_FFTW_ESTIMATE, nthreads)
1153 :
1154 1551949 : status = DftiCreateDescriptor(fw_plan, FFT_PRECISION, DFTI_COMPLEX, 1, nz)
1155 1551949 : DFTI_CHECK(status)
1156 :
1157 1551949 : status = DftiSetValue(fw_plan, DFTI_NUMBER_OF_TRANSFORMS, ngbout)
1158 : arr(1) = 0
1159 1551949 : arr(2) = mgb
1160 1551949 : status = DftiSetValue(fw_plan, DFTI_INPUT_STRIDES, arr)
1161 1551949 : status = DftiSetValue(fw_plan, DFTI_INPUT_DISTANCE, 1)
1162 1551949 : status = DftiSetValue(fw_plan, DFTI_NUMBER_OF_USER_THREADS, nt1)
1163 1551949 : if (iscale__ == 1) then
1164 1551949 : status = DftiSetValue(fw_plan, R_TO_G_SCALE, one/DBLE(nfft)) ! normalization
1165 : end if
1166 :
1167 1551949 : status = DftiCommitDescriptor(fw_plan)
1168 1551949 : DFTI_CHECK(status)
1169 :
1170 : !status = DftiComputeForward(fw_plan, cuboid)
1171 : !DFTI_CHECK(status)
1172 1551949 : DFTI_CHECK(compute_ip(fw_plan, R_TO_G, cuboid))
1173 :
1174 1551949 : status = DftiFreeDescriptor(fw_plan)
1175 :
1176 : !call fftw3_execute_dft(fw_plan, cuboid, cuboid)
1177 : !call fftw3_destroy_plan(fw_plan)
1178 :
1179 : ! Transfer the data in the output array.
1180 :
1181 : !$OMP PARALLEL DO PRIVATE(igb,iz,cidx)
1182 790313959 : do ig=1,npwout
1183 788762010 : igb = indpw_kout(4,ig)
1184 788762010 : iz = indpw_kout(3,ig)
1185 788762010 : cidx = igb + (iz-1)*mgb
1186 788762010 : fofgout(1,ig)=REAL (cuboid(cidx))
1187 790313959 : fofgout(2,ig)=AIMAG(cuboid(cidx))
1188 : end do
1189 : end if
1190 :
1191 1687471 : call dfti_free(cptr_cuboid)
1192 :
1193 1687471 : if (istwf_k/=2) then
1194 1065523 : ABI_FREE(pha1)
1195 1065523 : ABI_FREE(pha2)
1196 1065523 : ABI_FREE(pha3)
1197 : end if
1198 : !------------------------------------------------------------------
1199 :
1200 : end if ! End time-reversal symmetry
1201 :
1202 62924146 : if (option/=3) then
1203 62117210 : ABI_FREE(indpw_kin)
1204 : end if
1205 :
1206 178342257 : if (option==2 .or. option==3) then
1207 52493965 : ABI_FREE(indpw_kout)
1208 : end if
1209 :
1210 : contains
1211 :
1212 114611175 : integer function compute_ip(desc, isign, io_data) result(istat)
1213 : type(DFTI_DESCRIPTOR),pointer :: desc
1214 : integer,intent(in) :: isign
1215 : complex(MYKIND),intent(inout) :: io_data(*)
1216 :
1217 167105140 : select case (isign)
1218 : case (-1)
1219 52493965 : istat = DftiComputeForward(desc, io_data)
1220 : case (+1)
1221 62117210 : istat = DftiComputeBackward(desc, io_data)
1222 : case default
1223 : !ABI_ERROR(sjoin("Invalid isign:", itoa(isign)))
1224 114611175 : istat = DFTI_INVALID_CONFIGURATION
1225 : end select
1226 114611175 : end function compute_ip
1227 :
1228 5753311359 : integer function compute_op(desc, isign, in_data, out_data) result(istat)
1229 : type(DFTI_DESCRIPTOR),pointer :: desc
1230 : integer,intent(in) :: isign
1231 : complex(MYKIND),intent(in) :: in_data(*)
1232 : complex(MYKIND),intent(out) :: out_data(*)
1233 :
1234 8394456985 : select case (isign)
1235 : case (-1)
1236 2641145626 : istat = DftiComputeForward(desc, in_data, out_data)
1237 : case (+1)
1238 3112165733 : istat = DftiComputeBackward(desc, in_data, out_data)
1239 : case default
1240 : !ABI_ERROR(sjoin("Invalid isign:", itoa(isign)))
1241 5753311359 : istat = DFTI_INVALID_CONFIGURATION
1242 : end select
1243 5753311359 : end function compute_op
|
