Branch data 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 : 64013623 : integer,allocatable :: indpw_kin(:,:),indpw_kout(:,:)
19 : 64013623 : complex(MYKIND),allocatable :: pha1(:),pha2(:),pha3(:)
20 : 64013623 : complex(MYKIND),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 : 64013623 : G_TO_R = +1; R_TO_G = -1
29 : 64013623 : R_TO_G_SCALE = DFTI_FORWARD_SCALE; G_TO_R_SCALE = DFTI_BACKWARD_SCALE
30 : 64013623 : if (present(abi_convention)) then
31 [ - + - + : 214039 : 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 [ + + + + : 64013623 : iscale__ = 1; if (present(iscale)) iscale__ = iscale
- + ]
37 : :
38 [ - + - - : 64013623 : 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 [ + + - + : 64013623 : 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 : 64013623 : nx=ngfft(1); ny=ngfft(2); nz=ngfft(3); nfft=nx*ny*nz
52 : 64013623 : fftcache=ngfft(8)
53 : 64013623 : ldxy = ldx*ldy
54 : :
55 [ + + + + : 64013623 : if (option/=3) then
+ - ]
56 [ - + + - : 189620061 : ABI_MALLOC(indpw_kin,(4,npwin))
+ - + - -
+ - + - +
+ - + - +
- - + - +
- + + - +
- + - - +
- + ]
57 : 63206687 : call indfftrisc(gboundin(3:3+2*mgfft+4,1),indpw_kin,kg_kin,mgfft,ngbin,ngfft,npwin)
58 : : end if
59 : :
60 [ + + + + : 64013623 : if (option==2 .or. option==3) then
+ - ]
61 [ - + + - : 160312362 : ABI_MALLOC(indpw_kout,(4,npwout))
+ - + - -
+ - + - +
+ - + - +
- - + - +
- + + - +
- + - - +
- + ]
62 : 53437454 : 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 [ - + + + : 64013623 : if (option==2) then
+ - ]
68 : 52630518 : mgb=MAX(ngbin,ngbout)
69 [ + + + + : 11383105 : else if (option==0 .or. option==1) then
- - ]
70 : 10576169 : mgb=ngbin; ngbout=1
71 [ + - + - : 806936 : else if (option==3) then
- - ]
72 : 806936 : mgb=ngbout; ngbin=1
73 : : end if
74 : :
75 [ + + + + : 64013623 : if (MOD(mgb,2)/=1) mgb=mgb+1
+ + ]
76 : :
77 : : !Initialise openmp, if needed
78 : 64013623 : 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 : 64013623 : fraction=0.25
83 : : !First estimation of lot and nlot
84 : 64013623 : 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 : 64013623 : nlot=(ngbin-1)/lot+1
90 : 64013623 : nlot=nthreads*((nlot-1)/nthreads+1)
91 : 64013623 : lotin=(ngbin-1)/nlot+1
92 : 64013623 : nlot=(ngbout-1)/lot+1
93 : 64013623 : nlot=nthreads*((nlot-1)/nthreads+1)
94 : 64013623 : 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 [ + + + + : 64013623 : if (istwf_k==1) then
+ - ]
104 : :
105 : 62255910 : call dfti_alloc_complex(mgb*nz,cptr_cuboid,cuboid)
106 : :
107 [ + + + + : 62255910 : if (option/=3) then
+ - ]
108 : :
109 : : !$OMP PARALLEL DO
110 [ + + + + : 57957683357 : do cidx=1,mgb*nz
+ + ]
111 : 57957683357 : cuboid(cidx) = MYCZERO
112 : : end do
113 : :
114 : : ! Insert fofgin into the work array
115 : : !$OMP PARALLEL DO PRIVATE(igb,iz,cidx)
116 [ + + + + : 16156172964 : do ig=1,npwin
+ + ]
117 : 16094723990 : igb =indpw_kin(4,ig)
118 : 16094723990 : iz =indpw_kin(3,ig)
119 : 16094723990 : cidx = igb + (iz-1)*mgb
120 : 16156172964 : 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 : 61448974 : status = DftiCreateDescriptor(bw_plan, FFT_PRECISION, DFTI_COMPLEX, 1, nz)
126 : 61448974 : status = DftiSetValue(bw_plan, DFTI_NUMBER_OF_TRANSFORMS, ngbin)
127 : 61448974 : arr(1) = 0
128 : 61448974 : arr(2) = mgb
129 : 61448974 : status = DftiSetValue(bw_plan, DFTI_INPUT_STRIDES, arr)
130 : 61448974 : status = DftiSetValue(bw_plan, DFTI_INPUT_DISTANCE, 1)
131 : 61448974 : status = DftiSetValue(bw_plan, DFTI_NUMBER_OF_USER_THREADS, nt1)
132 [ - + - + : 61448974 : DFTI_CHECK(status)
- + ]
133 : :
134 : 61448974 : status = DftiCommitDescriptor(bw_plan)
135 [ - + - + : 61448974 : DFTI_CHECK(status)
- + ]
136 : :
137 : : !status = DftiComputeBackward(bw_plan, cuboid)
138 : : !DFTI_CHECK(status)
139 [ - + - + : 61448974 : 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 : 61448974 : status = DftiFreeDescriptor(bw_plan)
162 [ - + - + : 61448974 : 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 : 62255910 : g2min_in = gboundin(3,1)
173 : 62255910 : g2max_in = gboundin(4,1)
174 : :
175 : 62255910 : g2min_out = gboundout(3,1)
176 : 62255910 : g2max_out = gboundout(4,1)
177 : :
178 : : ! ng2_pos 1D transforms of f(Gx,Gy,z) along Gx for Gy >= 0
179 : 62255910 : ng2_pos = g2max_in + 1
180 : 62255910 : status = DftiCreateDescriptor(bw_plan_xg2pos , FFT_PRECISION, DFTI_COMPLEX, 1, nx)
181 [ - + - + : 62255910 : DFTI_CHECK(status)
- + ]
182 : :
183 : 62255910 : status = DftiSetValue(bw_plan_xg2pos, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
184 : 62255910 : status = DftiSetValue(bw_plan_xg2pos, DFTI_NUMBER_OF_TRANSFORMS, ng2_pos)
185 : 62255910 : status = DftiSetValue(bw_plan_xg2pos, DFTI_INPUT_DISTANCE, ldx)
186 : 62255910 : arr(1) = 0
187 : 62255910 : arr(2) = 1
188 : 62255910 : status = DftiSetValue(bw_plan_xg2pos, DFTI_INPUT_STRIDES, arr)
189 : 62255910 : status = DftiSetValue(bw_plan_xg2pos, DFTI_OUTPUT_DISTANCE, ldx)
190 : 62255910 : status = DftiSetValue(bw_plan_xg2pos, DFTI_OUTPUT_STRIDES, arr)
191 : 62255910 : status = DftiSetValue(bw_plan_xg2pos, DFTI_NUMBER_OF_USER_THREADS, nthreads)
192 : :
193 : 62255910 : status = DftiCommitDescriptor(bw_plan_xg2pos)
194 [ - + - + : 62255910 : DFTI_CHECK(status)
- + ]
195 : :
196 : : ! ng2_neg 1D transforms of f(Gx,Gy,z) along Gx for Gy < 0
197 : 62255910 : ng2_neg = -g2min_in
198 : 62255910 : status = DftiCreateDescriptor(bw_plan_xg2neg , FFT_PRECISION, DFTI_COMPLEX, 1, nx)
199 [ - + - + : 62255910 : DFTI_CHECK(status)
- + ]
200 : :
201 : 62255910 : status = DftiSetValue(bw_plan_xg2neg, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
202 : 62255910 : status = DftiSetValue(bw_plan_xg2neg, DFTI_NUMBER_OF_TRANSFORMS, ng2_neg)
203 : 62255910 : status = DftiSetValue(bw_plan_xg2neg, DFTI_INPUT_DISTANCE, ldx)
204 : : arr(1) = 0
205 : : arr(2) = 1
206 : 62255910 : status = DftiSetValue(bw_plan_xg2neg, DFTI_INPUT_STRIDES, arr)
207 : 62255910 : status = DftiSetValue(bw_plan_xg2neg, DFTI_OUTPUT_DISTANCE, ldx)
208 : 62255910 : status = DftiSetValue(bw_plan_xg2neg, DFTI_OUTPUT_STRIDES, arr)
209 : 62255910 : status = DftiSetValue(bw_plan_xg2neg, DFTI_NUMBER_OF_USER_THREADS, nthreads)
210 : :
211 : 62255910 : status = DftiCommitDescriptor(bw_plan_xg2neg)
212 [ - + - + : 62255910 : DFTI_CHECK(status)
- + ]
213 : :
214 : : ! nx 1D transforms of f(x,Gy,z) along Gy.
215 : 62255910 : status = DftiCreateDescriptor(bw_plan_y, FFT_PRECISION, DFTI_COMPLEX, 1, ny)
216 [ - + - + : 62255910 : DFTI_CHECK(status)
- + ]
217 : :
218 : 62255910 : status = DftiSetValue(bw_plan_y, DFTI_NUMBER_OF_TRANSFORMS, nx)
219 : 62255910 : status = DftiSetValue(bw_plan_y, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
220 : 62255910 : status = DftiSetValue(bw_plan_y, DFTI_INPUT_DISTANCE, 1)
221 : : arr(1) = 0
222 : 62255910 : arr(2) = ldx
223 : 62255910 : status = DftiSetValue(bw_plan_y, DFTI_INPUT_STRIDES, arr)
224 : 62255910 : status = DftiSetValue(bw_plan_y, DFTI_OUTPUT_DISTANCE, 1)
225 : 62255910 : status = DftiSetValue(bw_plan_y, DFTI_OUTPUT_STRIDES, arr)
226 : 62255910 : status = DftiSetValue(bw_plan_y, DFTI_NUMBER_OF_USER_THREADS, nthreads)
227 [ - + - + : 62255910 : DFTI_CHECK(status)
- + ]
228 : :
229 : 62255910 : status = DftiCommitDescriptor(bw_plan_y)
230 [ - + - + : 62255910 : DFTI_CHECK(status)
- + ]
231 : :
232 [ + + + + : 62255910 : if (option==2 .or. option==3) then
+ - ]
233 : : ! nx 1D transforms of f(x,y,Gz) along y.
234 : 51824241 : status = DftiCreateDescriptor(fw_plan_y , FFT_PRECISION, DFTI_COMPLEX, 1, ny)
235 [ - + - + : 51824241 : DFTI_CHECK(status)
- + ]
236 : :
237 : 51824241 : status = DftiSetValue(fw_plan_y, DFTI_NUMBER_OF_TRANSFORMS, nx)
238 : 51824241 : status = DftiSetValue(fw_plan_y, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
239 : 51824241 : status = DftiSetValue(fw_plan_y, DFTI_INPUT_DISTANCE, 1)
240 : : arr(1) = 0
241 : : arr(2) = ldx
242 : 51824241 : status = DftiSetValue(fw_plan_y, DFTI_INPUT_STRIDES, arr)
243 : 51824241 : status = DftiSetValue(fw_plan_y, DFTI_OUTPUT_DISTANCE, 1)
244 : 51824241 : status = DftiSetValue(fw_plan_y, DFTI_OUTPUT_STRIDES, arr)
245 : 51824241 : status = DftiSetValue(fw_plan_y, DFTI_NUMBER_OF_USER_THREADS, nthreads)
246 [ - + - + : 51824241 : DFTI_CHECK(status)
- + ]
247 : :
248 : 51824241 : status = DftiCommitDescriptor(fw_plan_y)
249 [ - + - + : 51824241 : DFTI_CHECK(status)
- + ]
250 : :
251 : : ! ng2_pos 1D transforms of f(x,y,Gz) along x.
252 : 51824241 : ng2_pos = g2max_out + 1
253 : 51824241 : status = DftiCreateDescriptor(fw_plan_xg2pos , FFT_PRECISION, DFTI_COMPLEX, 1, nx)
254 [ - + - + : 51824241 : DFTI_CHECK(status)
- + ]
255 : :
256 : 51824241 : status = DftiSetValue(fw_plan_xg2pos, DFTI_NUMBER_OF_TRANSFORMS, ng2_pos)
257 : 51824241 : status = DftiSetValue(fw_plan_xg2pos, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
258 : 51824241 : status = DftiSetValue(fw_plan_xg2pos, DFTI_INPUT_DISTANCE, ldx)
259 : : !status = DftiSetValue(fw_plan_xg2pos, DFTI_INPUT_STRIDES, arr)
260 : 51824241 : status = DftiSetValue(fw_plan_xg2pos, DFTI_OUTPUT_DISTANCE, ldx)
261 : 51824241 : status = DftiSetValue(fw_plan_xg2pos, DFTI_NUMBER_OF_USER_THREADS, nthreads)
262 [ - + - + : 51824241 : DFTI_CHECK(status)
- + ]
263 : :
264 : 51824241 : status = DftiCommitDescriptor(fw_plan_xg2pos)
265 [ - + - + : 51824241 : DFTI_CHECK(status)
- + ]
266 : :
267 : : ! ng2_neg 1D transforms of f(x,y,Gz) along x.
268 : 51824241 : ng2_neg = -g2min_out
269 : 51824241 : status = DftiCreateDescriptor(fw_plan_xg2neg , FFT_PRECISION, DFTI_COMPLEX, 1, nx)
270 [ - + - + : 51824241 : DFTI_CHECK(status)
- + ]
271 : :
272 : 51824241 : status = DftiSetValue(fw_plan_xg2neg, DFTI_NUMBER_OF_TRANSFORMS, ng2_neg)
273 : 51824241 : status = DftiSetValue(fw_plan_xg2neg, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
274 : 51824241 : status = DftiSetValue(fw_plan_xg2neg, DFTI_INPUT_DISTANCE, ldx)
275 : : !status = DftiSetValue(fw_plan_xg2neg, DFTI_INPUT_STRIDES, arr)
276 : 51824241 : status = DftiSetValue(fw_plan_xg2neg, DFTI_OUTPUT_DISTANCE, ldx)
277 : 51824241 : status = DftiSetValue(fw_plan_xg2neg, DFTI_NUMBER_OF_USER_THREADS, nthreads)
278 : :
279 : 51824241 : status = DftiCommitDescriptor(fw_plan_xg2neg)
280 [ - + - + : 51824241 : 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 : 62255910 : call dfti_alloc_complex(ldxy,cptr_sa,slice_a)
295 : 62255910 : call dfti_alloc_complex(ldxy,cptr_sb,slice_b)
296 : : #endif
297 : :
298 : : ! Big Loop over z parallelized with OpenMP.
299 : : !$OMP DO
300 [ + + + + : 1107174505 : do iz=1,nz
+ + ]
301 : :
302 [ + + + + : 1044918595 : if (option/=3) then
+ - ]
303 : : ! Zero the values on the current plane
304 [ + + + + : >38177*10^7 : slice_a(:) = MYCZERO
+ + ]
305 : : ! Copy the data in the current plane
306 [ + + + + : 58349866823 : do igb=1,ngbin
+ + ]
307 : 57318878240 : ix = indpw_kin(1,igb)
308 : 57318878240 : iy = indpw_kin(2,igb)
309 : 57318878240 : cidx = igb + (iz-1)*mgb
310 : 57318878240 : sidx = ix + (iy-1)*ldx
311 : 58349866823 : slice_a(sidx) = cuboid(cidx)
312 : : end do
313 : : ! Perform Gx transform, taking into account arrays of zeros
314 [ + - + - : 1030988583 : if (g2min_in+ny >= g2max_in+2) then
+ - ]
315 [ + + + + : 10444632219 : do iy=g2max_in+2,g2min_in+ny
+ + ]
316 : 9413643636 : padx = (iy-1)*ldx
317 [ + + + + : >19635*10^7 : do ix=1,nx
+ + ]
318 : >18590*10^7 : sidx = ix + padx
319 : >19532*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 [ - + - + : 1030988583 : DFTI_CHECK(compute_op(bw_plan_xg2pos, G_TO_R, slice_a, slice_b))
- + ]
327 : :
328 : 1030988583 : pt = 1 + (g2min_in+ny)*ldx
329 : : !status = DftiComputeBackward(bw_plan_xg2neg, slice_a(pt:), slice_b(pt:))
330 : : !DFTI_CHECK(status)
331 [ - + - + : 1030988583 : 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 [ - + - + : 1030988583 : 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 : 130999202 : SELECT CASE (option)
342 : : CASE (0)
343 : : ! Copy the transformed function at the right place and we are done!
344 : 130999202 : xyplane = 1 + (iz-1)*ldx*ldy
345 : : #ifndef HAVE_DFTI_MIXED_PRECISION
346 : 500 : if (MYKIND==dp) then
347 : 130998702 : 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 [ - - + + : 40821256 : if (abs(weight_r-weight_i)<tol12) then
- - ]
364 [ - - + + : 796109190 : do iy=1,ny
- - ]
365 : 755374102 : padx = (iy-1)*ldx
366 [ - - + + : 16805456219 : do ix=1,nx
- - ]
367 : 16009347029 : sidx = ix + padx
368 : 16764721131 : 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 [ - - + + : 859168125 : if (cplex==1) then
+ + ]
384 [ - - + + : 13545370759 : do iy=1,ny
+ + ]
385 : 12817608064 : padx = (iy-1)*ldx
386 [ - - + + : >27409*10^7 : do ix=1,nx
+ + ]
387 : >26055*10^7 : sidx = ix + padx
388 : >27336*10^7 : slice_a(sidx) = denpot(ix,iy,iz)*slice_a(sidx)
389 : : end do
390 : : end do
391 : : else
392 [ - - + + : 2348464692 : do iy=1,ny
+ + ]
393 : 2217059262 : padx = (iy-1)*ldx
394 [ - - + + : 43074291222 : do ix=1,nx
+ + ]
395 : 40725826530 : sidx = ix + padx
396 : 40725826530 : vloc = CMPLX(denpot(2*ix-1,iy,iz), denpot(2*ix,iy,iz), KIND=MYKIND)
397 : 42942885792 : 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 [ + - - + : 1040878893 : 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 [ + + + + : 1107174505 : 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 [ - + - + : 873098137 : 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 [ - + - + : 873098137 : DFTI_CHECK(compute_op(fw_plan_xg2pos, R_TO_G, slice_b, slice_a))
- + ]
432 : :
433 : 873098137 : pt = 1 + (g2min_out+ny)*ldx
434 : : !status = DftiComputeForward(fw_plan_xg2neg, slice_b(pt:), slice_a(pt:))
435 : : !DFTI_CHECK(status)
436 [ - + - + : 873098137 : 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 [ + + + + : 51225032404 : do igb=1,ngbout
+ + ]
440 : 50351934267 : ix=indpw_kout(1,igb)
441 : 50351934267 : iy=indpw_kout(2,igb)
442 : 50351934267 : sidx = ix + (iy-1)*ldx
443 : 50351934267 : cidx = igb + (iz-1)*mgb
444 : 51225032404 : 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 : 62255910 : call dfti_free(cptr_sa)
456 : 62255910 : call dfti_free(cptr_sb)
457 : : #endif
458 : : !$OMP END PARALLEL
459 : :
460 : : ! Free plans
461 : 62255910 : status = DftiFreeDescriptor(bw_plan_xg2pos)
462 : 62255910 : status = DftiFreeDescriptor(bw_plan_xg2neg)
463 : 62255910 : status = DftiFreeDescriptor(bw_plan_y)
464 : :
465 [ + + + + : 62255910 : if (option==2 .or. option==3)then
+ - ]
466 : 51824241 : status = DftiFreeDescriptor(fw_plan_y)
467 : 51824241 : status = DftiFreeDescriptor(fw_plan_xg2pos)
468 : 51824241 : status = DftiFreeDescriptor(fw_plan_xg2neg)
469 : : end if
470 : :
471 [ + + + + : 62255910 : 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 : 51824241 : status = DftiCreateDescriptor(fw_plan, FFT_PRECISION, DFTI_COMPLEX, 1, nz)
476 [ - + - + : 51824241 : DFTI_CHECK(status)
- + ]
477 : :
478 : 51824241 : status = DftiSetValue(fw_plan, DFTI_NUMBER_OF_TRANSFORMS, ngbout)
479 : : arr(1) = 0
480 : 51824241 : arr(2) = mgb
481 : 51824241 : status = DftiSetValue(fw_plan, DFTI_INPUT_STRIDES, arr)
482 : 51824241 : status = DftiSetValue(fw_plan, DFTI_INPUT_DISTANCE, 1)
483 : 51824241 : status = DftiSetValue(fw_plan, DFTI_NUMBER_OF_USER_THREADS, nt1)
484 : 51824241 : status = DftiSetValue(fw_plan, R_TO_G_SCALE, one/DBLE(nfft)) ! normalization
485 : :
486 : 51824241 : status = DftiCommitDescriptor(fw_plan)
487 [ - + - + : 51824241 : DFTI_CHECK(status)
- + ]
488 : :
489 : : !status = DftiComputeForward(fw_plan, cuboid)
490 : : !DFTI_CHECK(status)
491 [ - + - + : 51824241 : 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 : 51824241 : status = DftiFreeDescriptor(fw_plan)
518 : :
519 : : ! Transfer the data in the output array, after normalization
520 : : !$OMP PARALLEL DO PRIVATE(igb,iz,cidx)
521 [ + + + + : 14355888525 : do ig=1,npwout
+ + ]
522 : 14304064284 : igb = indpw_kout(4,ig)
523 : 14304064284 : iz = indpw_kout(3,ig)
524 : 14304064284 : cidx = igb + (iz-1)*mgb
525 : 14304064284 : fofgout(1,ig) = REAL (cuboid(cidx))
526 : 14355888525 : fofgout(2,ig) = AIMAG(cuboid(cidx))
527 : : end do
528 : : end if
529 : :
530 : 62255910 : call dfti_free(cptr_cuboid)
531 : :
532 [ - + - + : 1757713 : else if (istwf_k>=2) then
- - ]
533 : :
534 : : !------------------------------------------------------------------
535 : : !Here, use of time-reversal symmetry
536 : 1757713 : nxhalf1 = nx/2+1
537 : 1757713 : nxhalfm = (nx+1)/2
538 : 1757713 : nyhalf1 = ny/2+1
539 : :
540 : 1757713 : call dfti_alloc_complex(mgb*nz,cptr_cuboid,cuboid)
541 : :
542 [ - + + + : 1757713 : if (istwf_k/=2) then
- - ]
543 : : ! Precompute phases.
544 [ - - - - : 3331881 : ABI_MALLOC(pha1,(nx))
- - - - -
- - - - +
+ - + - +
- - + - +
- - - - -
- - - - -
- - ]
545 [ - - - - : 3331881 : ABI_MALLOC(pha2,(ny))
- - - - -
- - - - +
+ - + - +
- - + - +
- - - - -
- - - - -
- - ]
546 [ - - - - : 3331881 : ABI_MALLOC(pha3,(nz))
- - - - -
- - - - +
+ - + - +
- - + - +
- - - - -
- - - - -
- - ]
547 [ - - + + : 25793130 : do ix=1,nx
- - ]
548 : 24682503 : arg = DBLE(ix-1)*pi/DBLE(nx)
549 : 25793130 : pha1(ix)= MYCMPLX(COS(arg), SIN(arg))
550 : : end do
551 [ - - + + : 25780868 : do iy=1,ny
- - ]
552 : 24670241 : arg = DBLE(iy-1)*pi/DBLE(ny)
553 : 25780868 : pha2(iy)= MYCMPLX(COS(arg), SIN(arg))
554 : : end do
555 [ - - + + : 27158212 : do iz=1,nz
- - ]
556 : 26047585 : arg = DBLE(iz-1)*pi/DBLE(nz)
557 : 27158212 : pha3(iz)= MYCMPLX(COS(arg), SIN(arg))
558 : : end do
559 : : end if
560 : :
561 [ + - + - : 1757713 : if (option/=3) then
- - ]
562 : : ! Zero the components of cuboid
563 : : !$OMP PARALLEL DO
564 [ + + + + : 3641566122 : do cidx=1,mgb*nz
- - ]
565 : 3641566122 : cuboid(cidx) = MYCZERO
566 : : end do
567 : :
568 : : ! Insert fofgin into the work array cuboid
569 : : !$OMP PARALLEL DO PRIVATE(igb,iz,cidx)
570 [ + + + + : 890759180 : do ig=1,npwin
- - ]
571 : 889001467 : igb = indpw_kin(4,ig)
572 : 889001467 : iz = indpw_kin(3,ig)
573 : 889001467 : cidx = igb + (iz-1)*mgb
574 : 890759180 : 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 [ + - + + : 1757713 : if (istwf_k==2) then
- - ]
581 [ + + + + : 10290716 : do iz=nz/2+1,nz
- - ]
582 : 9643630 : iz_inv = nz+2-iz
583 : 9643630 : cidx = 1 + (iz-1)*mgb
584 : 9643630 : cidx_inv = 1 + (iz_inv-1)*mgb
585 : 10290716 : cuboid(cidx) = MYCONJG(cuboid(cidx_inv))
586 : : end do
587 [ - - + + : 1110627 : else if (istwf_k==4)then
- - ]
588 [ - - + + : 1586680 : do iz=nz/2+1,nz
- - ]
589 : 1469689 : iz_inv = nz+1-iz
590 : 1469689 : cidx = 1 + (iz-1)*mgb
591 : 1469689 : cidx_inv = 1 + (iz_inv-1)*mgb
592 : 1586680 : 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 : 1757713 : ixmax=gboundin(6,1)+1
598 : :
599 [ + - + + : 1757713 : if (istwf_k==2) then
- - ]
600 : : !$OMP PARALLEL DO PRIVATE(igb_inv,iz_inv,cidx,cidx_inv)
601 [ + + + + : 7525450 : do igb=2,2*ixmax-1
- - ]
602 : 6878364 : igb_inv=2*ixmax+1-igb
603 [ + + + + : 128042430 : do iz=nz/2+1,nz
- - ]
604 : 120516980 : iz_inv = nz+2-iz
605 : 120516980 : cidx = igb + (iz-1)*mgb
606 : 120516980 : cidx_inv = igb_inv + (iz_inv-1)*mgb
607 : 127395344 : cuboid(cidx) = MYCONJG(cuboid(cidx_inv))
608 : : end do
609 : : end do
610 : :
611 [ - - + + : 1110627 : else if (istwf_k==3) then
- - ]
612 : : !$OMP PARALLEL DO PRIVATE(igb_inv,iz_inv,cidx,cidx_inv)
613 [ - - + + : 3383240 : do igb=1,2*ixmax
- - ]
614 : 3074602 : igb_inv=2*ixmax+1-igb
615 [ - - + + : 48689410 : do iz=nz/2+1,nz
- - ]
616 : 45306170 : iz_inv = nz+2-iz
617 : 45306170 : cidx = igb + (iz-1)*mgb
618 : 45306170 : cidx_inv = igb_inv + (iz_inv-1)*mgb
619 : 48380772 : cuboid(cidx) = MYCONJG(cuboid(cidx_inv))
620 : : end do
621 : : end do
622 : :
623 [ - - + + : 801989 : else if (istwf_k==4) then
- - ]
624 : : !$OMP PARALLEL DO PRIVATE(igb_inv,iz_inv,cidx,cidx_inv)
625 [ - - + + : 951037 : do igb=2,2*ixmax-1
- - ]
626 : 834046 : igb_inv=2*ixmax+1-igb
627 [ - - + + : 11879061 : do iz=nz/2+1,nz
- - ]
628 : 10928024 : iz_inv = nz+1-iz
629 : 10928024 : cidx = igb + (iz-1)*mgb
630 : 10928024 : cidx_inv = igb_inv + (iz_inv-1)*mgb
631 : 11762070 : cuboid(cidx) = MYCONJG(cuboid(cidx_inv))
632 : : end do
633 : : end do
634 : :
635 [ - - + + : 684998 : else if (istwf_k==5) then
- - ]
636 : : !$OMP PARALLEL DO PRIVATE(igb_inv,iz_inv,cidx,cidx_inv)
637 [ - - + + : 905177 : do igb=1,2*ixmax
- - ]
638 : 806358 : igb_inv=2*ixmax+1-igb
639 [ - - + + : 9597393 : do iz=nz/2+1,nz
- - ]
640 : 8692216 : iz_inv = nz+1-iz
641 : 8692216 : cidx = igb + (iz-1)*mgb
642 : 8692216 : cidx_inv = igb_inv + (iz_inv-1)*mgb
643 : 9498574 : cuboid(cidx) = MYCONJG(cuboid(cidx_inv))
644 : : end do
645 : : end do
646 : :
647 : : end if
648 : :
649 : : ! Now, iz==1
650 [ + - + + : 1757713 : if (istwf_k==2) then
- - ]
651 [ + + + + : 4086268 : do igb=2,ixmax
- - ]
652 : 3439182 : igb_inv = 2*ixmax+1-igb
653 : 3439182 : cidx_inv = igb_inv !+ (iz_inv-1)*mgb
654 : 4086268 : cuboid(cidx_inv) = MYCONJG(cuboid(igb))
655 : : end do
656 [ - - + + : 1110627 : else if (istwf_k==3) then
- - ]
657 [ - - + + : 1845939 : do igb=1,ixmax
- - ]
658 : 1537301 : igb_inv = 2*ixmax+1-igb
659 : 1537301 : cidx_inv = igb_inv !+ (iz_inv-1)*mgb
660 : 1845939 : 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 : 1757713 : status = DftiCreateDescriptor(bw_plan, FFT_PRECISION, DFTI_COMPLEX, 1, nz)
671 : 1757713 : status = DftiSetValue(bw_plan, DFTI_NUMBER_OF_TRANSFORMS, ngbin)
672 : 1757713 : arr(1) = 0
673 : 1757713 : arr(2) = mgb
674 : 1757713 : status = DftiSetValue(bw_plan, DFTI_INPUT_STRIDES, arr)
675 : 1757713 : status = DftiSetValue(bw_plan, DFTI_INPUT_DISTANCE, 1)
676 : 1757713 : status = DftiSetValue(bw_plan, DFTI_NUMBER_OF_USER_THREADS, nt1)
677 [ - + - + : 1757713 : DFTI_CHECK(status)
- - ]
678 : :
679 : 1757713 : status = DftiCommitDescriptor(bw_plan)
680 [ - + - + : 1757713 : DFTI_CHECK(status)
- - ]
681 : :
682 : : !status = DftiComputeBackward(bw_plan, cuboid)
683 : : !DFTI_CHECK(status)
684 [ - + - + : 1757713 : DFTI_CHECK(compute_ip(bw_plan, G_TO_R, cuboid))
- - ]
685 : :
686 : 1757713 : status = DftiFreeDescriptor(bw_plan)
687 [ - + - + : 1757713 : 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 [ - + + + : 1757713 : 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 [ - - + + : 10421510 : do iz=1,nz
- - ]
696 : 9949722 : phase = pha3(iz)
697 : 9949722 : padz = (iz-1)*mgb
698 [ - - + + : 424078391 : do igb=1,ngbin
- - ]
699 : 413656881 : cidx = igb + padz
700 : 423606603 : cuboid(cidx) = cuboid(cidx) * phase
701 : : end do
702 : : end do
703 : : end if
704 : :
705 : : end if ! if(option/=3)
706 : :
707 : 1757713 : g2max_in= gboundin(4,1)
708 : 1757713 : 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 : 1757713 : status = DftiCreateDescriptor(bw_plan_xg2pos , FFT_PRECISION, DFTI_COMPLEX, 1, nx)
715 [ - + - + : 1757713 : DFTI_CHECK(status)
- - ]
716 : :
717 : 1757713 : status = DftiSetValue(bw_plan_xg2pos, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
718 : 1757713 : status = DftiSetValue(bw_plan_xg2pos, DFTI_NUMBER_OF_TRANSFORMS, ng2_pos)
719 : 1757713 : status = DftiSetValue(bw_plan_xg2pos, DFTI_INPUT_DISTANCE, ldx)
720 : 1757713 : arr(1) = 0
721 : 1757713 : arr(2) = 1
722 : 1757713 : status = DftiSetValue(bw_plan_xg2pos, DFTI_INPUT_STRIDES, arr)
723 : 1757713 : status = DftiSetValue(bw_plan_xg2pos, DFTI_OUTPUT_DISTANCE, ldx)
724 : 1757713 : status = DftiSetValue(bw_plan_xg2pos, DFTI_OUTPUT_STRIDES, arr)
725 : 1757713 : status = DftiSetValue(bw_plan_xg2pos, DFTI_NUMBER_OF_USER_THREADS, nthreads)
726 : :
727 : 1757713 : status = DftiCommitDescriptor(bw_plan_xg2pos)
728 [ - + - + : 1757713 : 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 : 1757713 : status = DftiCreateDescriptor(bw_plan_y, FFT_PRECISION, DFTI_COMPLEX, 1, ny)
736 [ - + - + : 1757713 : DFTI_CHECK(status)
- - ]
737 : :
738 : 1757713 : status = DftiSetValue(bw_plan_y, DFTI_NUMBER_OF_TRANSFORMS, nxhalf1)
739 : 1757713 : status = DftiSetValue(bw_plan_y, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
740 : 1757713 : status = DftiSetValue(bw_plan_y, DFTI_INPUT_DISTANCE, 1)
741 : : arr(1) = 0
742 : 1757713 : arr(2) = ldx
743 : 1757713 : status = DftiSetValue(bw_plan_y, DFTI_INPUT_STRIDES, arr)
744 : 1757713 : status = DftiSetValue(bw_plan_y, DFTI_OUTPUT_DISTANCE, 1)
745 : 1757713 : status = DftiSetValue(bw_plan_y, DFTI_OUTPUT_STRIDES, arr)
746 : 1757713 : status = DftiSetValue(bw_plan_y, DFTI_NUMBER_OF_USER_THREADS, nthreads)
747 [ - + - + : 1757713 : DFTI_CHECK(status)
- - ]
748 : :
749 : 1757713 : status = DftiCommitDescriptor(bw_plan_y)
750 [ - + - + : 1757713 : DFTI_CHECK(status)
- - ]
751 : :
752 [ - + + + : 1757713 : 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 : 1613213 : status = DftiCreateDescriptor(fw_plan_y , FFT_PRECISION, DFTI_COMPLEX, 1, ny)
759 [ - - - + : 1613213 : DFTI_CHECK(status)
- - ]
760 : :
761 : 1613213 : status = DftiSetValue(fw_plan_y, DFTI_NUMBER_OF_TRANSFORMS, nxhalf1)
762 : 1613213 : status = DftiSetValue(fw_plan_y, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
763 : 1613213 : status = DftiSetValue(fw_plan_y, DFTI_INPUT_DISTANCE, 1)
764 : : arr(1) = 0
765 : : arr(2) = ldx
766 : 1613213 : status = DftiSetValue(fw_plan_y, DFTI_INPUT_STRIDES, arr)
767 : 1613213 : status = DftiSetValue(fw_plan_y, DFTI_OUTPUT_DISTANCE, 1)
768 : 1613213 : status = DftiSetValue(fw_plan_y, DFTI_OUTPUT_STRIDES, arr)
769 : 1613213 : status = DftiSetValue(fw_plan_y, DFTI_NUMBER_OF_USER_THREADS, nthreads)
770 [ - - - + : 1613213 : DFTI_CHECK(status)
- - ]
771 : :
772 : 1613213 : status = DftiCommitDescriptor(fw_plan_y)
773 [ - - - + : 1613213 : 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 : 1613213 : status = DftiCreateDescriptor(fw_plan_xg2pos , FFT_PRECISION, DFTI_COMPLEX, 1, nx)
783 [ - - - + : 1613213 : DFTI_CHECK(status)
- - ]
784 : :
785 : 1613213 : status = DftiSetValue(fw_plan_xg2pos, DFTI_NUMBER_OF_TRANSFORMS, ng2_pos)
786 : 1613213 : status = DftiSetValue(fw_plan_xg2pos, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
787 : 1613213 : 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 : 1613213 : status = DftiSetValue(fw_plan_xg2pos, DFTI_OUTPUT_DISTANCE, 2*nxhalfm)
792 : 1613213 : status = DftiSetValue(fw_plan_xg2pos, DFTI_NUMBER_OF_USER_THREADS, nthreads)
793 [ - - - + : 1613213 : DFTI_CHECK(status)
- - ]
794 : :
795 : 1613213 : status = DftiCommitDescriptor(fw_plan_xg2pos)
796 [ - - - + : 1613213 : 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 : 1757713 : call dfti_alloc_complex(ldxy,cptr_sa,slice_a)
813 : 1757713 : call dfti_alloc_complex(ldxy,cptr_sb,slice_b)
814 : :
815 : 1757713 : call dfti_alloc_complex(2*nxhalfm*ldy,cptr_sc,slice_c)
816 : 1757713 : call dfti_alloc_complex(2*nxhalfm*ldy,cptr_sd,slice_d)
817 : : #endif
818 : :
819 : : !$OMP DO
820 [ + + + + : 47027888 : do iz=1,nz
- - ]
821 : :
822 [ + - + - : 45270175 : 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 [ + + + + : 45236629120 : slice_a = MYCZERO
- - ]
833 : : !
834 : : ! Copy the data in the current plane
835 [ + + + + : 3653802184 : do igb=1,ngbin
- - ]
836 : 3608532009 : ix = indpw_kin(1,igb)
837 : 3608532009 : iy = indpw_kin(2,igb)
838 : 3608532009 : cidx = igb + (iz-1)*mgb
839 : 3608532009 : sidx = ix + (iy-1)*ldx
840 : 3653802184 : 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 [ - + - + : 45270175 : 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 [ - + + + : 45270175 : if (istwf_k==3 .or. istwf_k==5 .or. istwf_k==7 .or. istwf_k==9) then
- - ]
853 [ - - + + : 150512754 : do iy=1,g2max_in+1
- - ]
854 : 130970478 : padx = (iy-1)*ldx
855 [ - - + + : 4274650545 : do ix=1,nx
- - ]
856 : 4124137791 : phase = pha1(ix)
857 : 4124137791 : sidx = ix + padx
858 : 4255108269 : slice_b(sidx) = slice_b(sidx) * phase
859 : : end do
860 : : end do
861 : : end if
862 : :
863 : : ! Compute symmetric and antisymmetric combinations
864 [ + - + + : 45270175 : if (istwf_k>=2 .and. istwf_k<=5) then
- - ]
865 [ + + + + : 494553358 : do ix=1,nxhalf1-1
- - ]
866 : 494553358 : 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 [ - + + + : 31943234 : if ((2*nxhalf1-2)/=nx) then
- - ]
870 : 2831913 : slice_a(nxhalf1)= REAL(slice_b(nx))
871 : : end if
872 : : iiy=2
873 : : else
874 : : iiy=1
875 : : end if
876 : :
877 [ + - + - : 45270175 : if (g2max_in+1 >= iiy) then
- - ]
878 [ + + + + : 333586523 : do iy=iiy,g2max_in+1
- - ]
879 : : !
880 : 288316348 : padx = (iy-1)*ldx
881 [ + + + + : 5086666462 : do ix=1,nxhalf1-1
- - ]
882 : 4798350114 : sidx = ix + padx
883 : 4798350114 : cidx = 2*ix-1 + padx
884 : :
885 : : slice_a(sidx) = MYCMPLX( REAL(slice_b(cidx)) - AIMAG(slice_b(cidx+1)),&
886 : 4798350114 : & AIMAG(slice_b(cidx)) + REAL (slice_b(cidx+1)) )
887 : :
888 : 4798350114 : sidx = ix + (ny+iiy-iy-1)*ldx
889 : : slice_a(sidx) = MYCMPLX( REAL(slice_b(cidx)) + AIMAG(slice_b(cidx+1)),&
890 : 5086666462 : & -AIMAG(slice_b(cidx)) + REAL (slice_b(cidx+1)) )
891 : : end do
892 : : !
893 [ - + + + : 333586523 : if ((2*nxhalf1-2)/=nx) then
- - ]
894 : 23148498 : sidx = nxhalf1 + padx
895 : 23148498 : cidx = nx + padx
896 : 23148498 : slice_a(sidx) = slice_b(cidx)
897 : :
898 : 23148498 : sidx = nxhalf1 + (ny+iiy-iy-1)*ldx
899 : 23148498 : slice_a(sidx) = MYCONJG(slice_b(cidx))
900 : : end if
901 : : end do
902 : : end if
903 : :
904 [ + - + - : 45270175 : if (nyhalf1 >= g2max_in+2) then
- - ]
905 [ + + + + : 415163157 : do iy=g2max_in+2,nyhalf1
- - ]
906 [ + + + + : 6389868932 : do ix=1,nxhalf1-1
- - ]
907 : 6019975950 : sidx = ix + (iy-1)*ldx
908 : 6019975950 : slice_a(sidx) = MYCZERO
909 : :
910 : 6019975950 : sidx = ix + (ny+iiy-iy-1)*ldx
911 : 6389868932 : slice_a(sidx) =zero
912 : : end do
913 [ - + + + : 415163157 : if ((2*nxhalf1-2)/=nx) then
- - ]
914 : 28319825 : sidx = nxhalf1 + (iy-1)*ldx
915 : 28319825 : slice_a(sidx) = MYCZERO
916 : :
917 : 28319825 : sidx = nxhalf1 + (ny+iiy-iy-1)*ldx
918 : 28319825 : 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 [ - + - + : 45270175 : 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 [ - + + + : 45270175 : if (istwf_k>=6 .and. istwf_k<=9) then
- - ]
933 [ - - + + : 355387810 : do iy=1,ny
- - ]
934 : 342060869 : phase = pha2(iy)
935 : 342060869 : padx = (iy-1)*ldx
936 [ - - + + : 5402112237 : do ix=1,nxhalfm
- - ]
937 : 5046724427 : sidx = ix + padx
938 : 5388785296 : 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 [ + - + + : 45270175 : 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 [ + + + + : 659700 : do iy=1,ny
- - ]
952 : 638727 : padx = (iy-1)*ldx
953 [ + + + + : 15136134 : do ix=1,nxhalf1-1
- - ]
954 : 14497407 : sidx = ix + (iy-1)*ldx
955 : 14497407 : fidx = 2*ix-1 + padx + (iz-1)*ldx*ldy
956 : :
957 : 14497407 : fofr(1,fidx) = REAL(slice_b(sidx))
958 : 14497407 : fofr(2,fidx) = zero
959 : :
960 : 14497407 : fofr(1,fidx+1) = AIMAG(slice_b(sidx))
961 : 15136134 : fofr(2,fidx+1) = zero
962 : : end do
963 : : ! If nx odd, must add last data
964 [ - + + + : 659700 : if ((2*nxhalf1-2)/=nx) then
- - ]
965 : 98271 : sidx = nxhalf1 + (iy-1)*ldx
966 : 98271 : fidx = nx + (iy-1) * ldx + (iz-1)*ldx*ldy
967 : 98271 : fofr(1,fidx) = REAL(slice_b(sidx))
968 : 98271 : fofr(2,fidx) = zero
969 : : end if
970 : : end do
971 : : end if
972 : :
973 [ - + + + : 45270175 : if (option==1) then ! Accumulate density
- - ]
974 [ - - + + : 126653078 : do iy=1,ny
- - ]
975 : 122552582 : padx = (iy-1)*ldx
976 [ - - + + : 2188879309 : do ix=1,nxhalf1-1
- - ]
977 : 2066326727 : sidx = ix + padx
978 : 2066326727 : denpot(2*ix-1,iy,iz) = denpot(2*ix-1,iy,iz) + weight_r* REAL(slice_b(sidx))**2
979 : 2188879309 : 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 [ - - + + : 126653078 : 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 [ - + + + : 45270175 : if (option==2) then ! Apply local potential
- - ]
990 [ - - + + : 1211711755 : do iy=1,ny
- - ]
991 : 1170563049 : padx = (iy-1)*ldx
992 [ - - + + : 20445313943 : do ix=1,nxhalf1-1
- - ]
993 : 19274750894 : sidx = ix + padx
994 : : slice_a(sidx) = MYCMPLX( denpot(2*ix-1,iy,iz)* REAL(slice_b(sidx)), &
995 : 20445313943 : denpot(2*ix ,iy,iz)*AIMAG(slice_b(sidx)) )
996 : : end do
997 : : ! If nx odd, must add last data
998 [ - - + + : 1211711755 : 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 [ - + - + : 45270175 : 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 [ - + + + : 47027888 : if (option==2 .or. option==3) then
- - ]
1025 : : !
1026 : : ! Change the phase if $k_y \neq 0$
1027 [ - - + + : 41148706 : if (istwf_k>=6 .and. istwf_k<=9) then
- - ]
1028 [ - - + + : 324817380 : do iy=1,ny
- - ]
1029 : 312538593 : padx = (iy-1)*ldx
1030 : 312538593 : phase = MYCONJG(pha2(iy))
1031 [ - - + + : 4901867863 : do ix=1,nxhalfm
- - ]
1032 : 4577050483 : sidx = ix + padx
1033 : 4889589076 : 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 [ - - - + : 41148706 : DFTI_CHECK(compute_op(fw_plan_y, R_TO_G, slice_a, slice_b))
- - ]
1046 : :
1047 : : ! Decompose symmetric and antisymmetric parts
1048 [ - - + + : 41148706 : if (istwf_k>=2 .and. istwf_k<=5) then
- - ]
1049 [ - - + + : 448867646 : do ix=1,nxhalfm
- - ]
1050 : 419997727 : pidx = 2*ix-1
1051 : 419997727 : slice_c(pidx) = MYCMPLX( REAL(slice_b(ix)), zero)
1052 : 448867646 : 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 [ - - + + : 302094664 : do iy=iiy,g2max_in+1
- - ]
1060 : 260945958 : padx = (iy-1)*ldx
1061 [ - - + + : 4655112356 : do ix=1,nxhalfm
- - ]
1062 : 4353017692 : sidx = ix + padx
1063 : 4353017692 : cidx = ix + (ny+iiy-iy-1)*ldx
1064 : 4353017692 : 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 : 4353017692 : & (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 : 4613963650 : & (-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 [ - - + + : 41148706 : if (istwf_k==3 .or. istwf_k==5 .or. istwf_k==7 .or. istwf_k==9) then
- - ]
1078 [ - - + + : 137269634 : do iy=1,g2max_in+1
- - ]
1079 [ - - + + : 3871449429 : do ix=1,nx
- - ]
1080 : 3734179795 : phase = MYCONJG(pha1(ix))
1081 : 3734179795 : pidx = ix + (iy-1)*2*nxhalfm
1082 : 3853522143 : 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 [ - - - + : 41148706 : 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 [ - - + + : 3304287633 : do igb=1,ngbout
- - ]
1097 : 3263138927 : ix=indpw_kout(1,igb)
1098 : 3263138927 : iy=indpw_kout(2,igb)
1099 : 3263138927 : cidx = igb + (iz-1)*mgb
1100 : 3263138927 : pidx = ix + (iy-1)*2*nxhalfm
1101 : 3304287633 : 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 : 1757713 : call dfti_free(cptr_sa)
1116 : 1757713 : call dfti_free(cptr_sb)
1117 : 1757713 : call dfti_free(cptr_sc)
1118 : 1757713 : call dfti_free(cptr_sd)
1119 : : #endif
1120 : : !$OMP END PARALLEL
1121 : :
1122 : 1757713 : status = DftiFreeDescriptor(bw_plan_xg2pos)
1123 : 1757713 : status = DftiFreeDescriptor(bw_plan_y)
1124 : :
1125 : : !call fftw3_destroy_plan(bw_plan_xg2pos)
1126 : : !call fftw3_destroy_plan(bw_plan_y)
1127 : :
1128 [ - + + + : 1757713 : if (option==2 .or. option==3) then
- - ]
1129 : :
1130 : 1613213 : status = DftiFreeDescriptor(fw_plan_y)
1131 : 1613213 : 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 [ - - + + : 1613213 : 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 [ - - + + : 9594385 : do iz=1,nz
- - ]
1140 : 9154364 : phase = MYCONJG(pha3(iz))
1141 : 9154364 : padc = (iz-1)*mgb
1142 [ - - + + : 388436714 : do igb=1,ngbout
- - ]
1143 : 378842329 : cidx = igb + padc
1144 : 387996693 : 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 : 1613213 : status = DftiCreateDescriptor(fw_plan, FFT_PRECISION, DFTI_COMPLEX, 1, nz)
1155 [ - - - + : 1613213 : DFTI_CHECK(status)
- - ]
1156 : :
1157 : 1613213 : status = DftiSetValue(fw_plan, DFTI_NUMBER_OF_TRANSFORMS, ngbout)
1158 : : arr(1) = 0
1159 : 1613213 : arr(2) = mgb
1160 : 1613213 : status = DftiSetValue(fw_plan, DFTI_INPUT_STRIDES, arr)
1161 : 1613213 : status = DftiSetValue(fw_plan, DFTI_INPUT_DISTANCE, 1)
1162 : 1613213 : status = DftiSetValue(fw_plan, DFTI_NUMBER_OF_USER_THREADS, nt1)
1163 [ - - + - : 1613213 : if (iscale__ == 1) then
- - ]
1164 : 1613213 : status = DftiSetValue(fw_plan, R_TO_G_SCALE, one/DBLE(nfft)) ! normalization
1165 : : end if
1166 : :
1167 : 1613213 : status = DftiCommitDescriptor(fw_plan)
1168 [ - - - + : 1613213 : DFTI_CHECK(status)
- - ]
1169 : :
1170 : : !status = DftiComputeForward(fw_plan, cuboid)
1171 : : !DFTI_CHECK(status)
1172 [ - - - + : 1613213 : DFTI_CHECK(compute_ip(fw_plan, R_TO_G, cuboid))
- - ]
1173 : :
1174 : 1613213 : 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 [ - - + + : 802984612 : do ig=1,npwout
- - ]
1183 : 801371399 : igb = indpw_kout(4,ig)
1184 : 801371399 : iz = indpw_kout(3,ig)
1185 : 801371399 : cidx = igb + (iz-1)*mgb
1186 : 801371399 : fofgout(1,ig)=REAL (cuboid(cidx))
1187 : 802984612 : fofgout(2,ig)=AIMAG(cuboid(cidx))
1188 : : end do
1189 : : end if
1190 : :
1191 : 1757713 : call dfti_free(cptr_cuboid)
1192 : :
1193 [ - + + + : 1757713 : if (istwf_k/=2) then
- - ]
1194 [ - - - + : 1110627 : ABI_FREE(pha1)
- - ]
1195 [ - - - + : 1110627 : ABI_FREE(pha2)
- - ]
1196 [ - - - + : 1110627 : ABI_FREE(pha3)
- - ]
1197 : : end if
1198 : : !------------------------------------------------------------------
1199 : :
1200 : : end if ! End time-reversal symmetry
1201 : :
1202 [ + + + + : 64013623 : if (option/=3) then
+ - ]
1203 [ - + - + : 63206687 : ABI_FREE(indpw_kin)
- + ]
1204 : : end if
1205 : :
1206 [ + + + + : 181464700 : if (option==2 .or. option==3) then
- + - + -
+ - + - +
+ + + + -
+ - + - +
- + - + -
+ + - - +
- + - + -
+ - + ]
1207 [ - + - + : 53437454 : ABI_FREE(indpw_kout)
- + ]
1208 : : end if
1209 : :
1210 : : contains
1211 : :
1212 : 116644141 : 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 : 170081595 : select case (isign)
1218 : : case (-1)
1219 : 53437454 : istat = DftiComputeForward(desc, io_data)
1220 : : case (+1)
1221 : 63206687 : istat = DftiComputeBackward(desc, io_data)
1222 : : case default
1223 : : !ABI_ERROR(sjoin("Invalid isign:", itoa(isign)))
1224 [ + + - ]: 116644141 : istat = DFTI_INVALID_CONFIGURATION
1225 : : end select
1226 : 116644141 : end function compute_ip
1227 : :
1228 : 5885097922 : 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 : 8586689745 : select case (isign)
1235 : : case (-1)
1236 : 2701591823 : istat = DftiComputeForward(desc, in_data, out_data)
1237 : : case (+1)
1238 : 3183506099 : istat = DftiComputeBackward(desc, in_data, out_data)
1239 : : case default
1240 : : !ABI_ERROR(sjoin("Invalid isign:", itoa(isign)))
1241 [ + + - ]: 5885097922 : istat = DFTI_INVALID_CONFIGURATION
1242 : : end select
1243 : 5885097922 : end function compute_op
|
