Branch data Line data Source code
1 : : #ifdef HAVE_DFTI
2 : : integer :: kk,ii,jj,sidx,cnt,dat,line,zplane,status,nthreads, iscale__,arr(2)
3 : : character(len=500) :: msg
4 : 5739117 : type(zpad_t) :: zpad
5 : : type(DFTI_DESCRIPTOR),pointer :: Desc1,Desc2,Desc3
6 : :
7 [ + + + + : 5741321 : iscale__ = merge(1, 0, isign == -1); if (present(iscale)) iscale__ = iscale
+ + + + ]
8 : :
9 : 5739117 : nthreads = xomp_get_num_threads(open_parallel=.TRUE.)
10 : :
11 : 5739117 : call zpad%init(nx, ny, nz, ldx, ldy, ldz, mgfft, gbound)
12 : :
13 : 2204 : select case (isign)
14 : : case (+1) ! G --> R
15 : : ! 1) Transform along x.
16 : : ! 1D transform of f(Gx,Gy,Gz) along Gx for each Gx-line intersecting the G-sphere
17 : :
18 : 2204 : status = DftiCreateDescriptor(Desc1, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, nx)
19 [ - + - + ]: 2204 : DFTI_CHECK(status)
20 : : !
21 : : ! TODO
22 : : ! this value was added in mkl version 11.0.
23 : : ! It's the only way to avoid side effect in mkl_set_num_threads thus making this routine thread safe
24 : : !status = DftiSetValue(Desc1, DFTI_THREAD_LIMIT, 1)
25 : 2204 : status = DftiSetValue(Desc1, DFTI_NUMBER_OF_USER_THREADS, nthreads)
26 : :
27 : 2204 : status = DftiCommitDescriptor(Desc1)
28 [ - + - + ]: 2204 : DFTI_CHECK(status)
29 : :
30 : : !$OMP PARALLEL DO PRIVATE(dat,line,jj,kk,sidx,status)
31 [ + + + + ]: 411771 : do cnt=0,(zpad%nlinex*ndat)-1
32 : 409567 : dat = 1 + cnt / zpad%nlinex
33 : 409567 : line = 1 + MOD(cnt, zpad%nlinex)
34 : 409567 : jj = zpad%linex2ifft_yz(1,line)
35 : 409567 : kk = zpad%linex2ifft_yz(2,line)
36 : 409567 : sidx = 1+ (jj-1)*ldx + (kk-1)*ldx*ldy + (dat-1) * ldx*ldy*ldz
37 : 409567 : status = DftiComputeBackward(Desc1, ff(sidx:)) ! Pass the base address.
38 [ - + - + ]: 411771 : DFTI_CHECK(status)
39 : : end do
40 : :
41 : 2204 : status = DftiFreeDescriptor(Desc1)
42 [ - + - + ]: 2204 : DFTI_CHECK(status)
43 : : !
44 : : ! 2) Transform along y:
45 : : ! nx 1D transforms of f(x,Gy,Gz) along Gy for each Gz plane intersecting the sphere.
46 : :
47 : 2204 : status = DftiCreateDescriptor(Desc2, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, ny)
48 [ - + - + ]: 2204 : DFTI_CHECK(status)
49 : :
50 : : !status = DftiSetValue(Desc2, DFTI_THREAD_LIMIT, 1)
51 : 2204 : status = DftiSetValue(Desc2, DFTI_NUMBER_OF_TRANSFORMS, nx)
52 : 2204 : status = DftiSetValue(Desc2, DFTI_INPUT_DISTANCE, 1)
53 : 2204 : arr(1) = 0
54 : 2204 : arr(2) = ldx
55 : 2204 : status = DftiSetValue(Desc2, DFTI_INPUT_STRIDES, arr)
56 : 2204 : status = DftiSetValue(Desc2, DFTI_NUMBER_OF_USER_THREADS, nthreads)
57 : :
58 : 2204 : status = DftiCommitDescriptor(Desc2)
59 [ - + - + ]: 2204 : DFTI_CHECK(status)
60 : :
61 : : !$OMP PARALLEL DO PRIVATE(dat,zplane,kk,sidx,status)
62 [ + + + + ]: 33602 : do cnt=0,(zpad%n_zplanes*ndat)-1
63 : 31398 : dat = 1 + cnt / zpad%n_zplanes
64 : 31398 : zplane = 1 + MOD(cnt, zpad%n_zplanes)
65 : 31398 : kk = zpad%zplane(1,zplane)
66 : 31398 : sidx = 1 + (kk-1)*ldx*ldy + (dat-1) *ldx*ldy*ldz
67 : 31398 : status = DftiComputeBackward(Desc2, ff(sidx:))
68 [ - + - + ]: 33602 : DFTI_CHECK(status)
69 : : end do
70 : :
71 : 2204 : status = DftiFreeDescriptor(Desc2)
72 [ - + - + ]: 2204 : DFTI_CHECK(status)
73 : : !
74 : : ! 3) Transform along z.
75 : : ! ldx*ldy 1D transforms of f(x,y,Gz) along Gz.
76 : : ! Note that we have to visit the entire augmented x-y plane!
77 : :
78 [ + + - + : 2204 : if (ndat==1 .or. nthreads==1 .or. MOD(ndat,nthreads)/=0) then
- - + + -
+ - - ]
79 : : !
80 : : ! Use MKL internal threading if single FFT or ndat is not divisible by nthreads
81 : : !status = DftiSetValue(Desc3, DFTI_THREAD_LIMIT, nthreads)
82 : 2204 : status = DftiCreateDescriptor(Desc3, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, nz)
83 [ - + - + ]: 2204 : DFTI_CHECK(status)
84 : :
85 : 2204 : status = DftiSetValue(Desc3, DFTI_NUMBER_OF_TRANSFORMS, ldx*ldy)
86 : 2204 : status = DftiSetValue(Desc3, DFTI_INPUT_DISTANCE, 1)
87 : : arr(1) = 0
88 : 2204 : arr(2) = ldx*ldy
89 : 2204 : status = DftiSetValue(Desc3, DFTI_INPUT_STRIDES, arr)
90 : :
91 : : ! Normalize the transform.
92 [ - + - + ]: 2204 : if (iscale__ /= 0) then
93 : 0 : status = DftiSetValue(Desc3, DFTI_BACKWARD_SCALE, one/(nx*ny*nz))
94 : : end if
95 : :
96 : 2204 : status = DftiCommitDescriptor(Desc3)
97 [ - + - + ]: 2204 : DFTI_CHECK(status)
98 : :
99 [ + + + + ]: 4450 : do dat=1,ndat
100 : 2246 : sidx = 1 + (dat-1) *ldx*ldy*ldz
101 : 2246 : status = DftiComputeBackward(Desc3, ff(sidx:))
102 [ - + - + ]: 4450 : DFTI_CHECK(status)
103 : : end do
104 : :
105 : 2204 : status = DftiFreeDescriptor(Desc3)
106 [ - + - + ]: 2204 : DFTI_CHECK(status)
107 : :
108 : : else
109 : : ! Split ndat transforms among the threads.
110 : : !call mkl_set_num_threads(1)
111 : 0 : status = DftiCreateDescriptor(Desc3, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, nz)
112 [ # # # # ]: 0 : DFTI_CHECK(status)
113 : :
114 : 0 : status = DftiSetValue(Desc3, DFTI_NUMBER_OF_TRANSFORMS, ldx*ldy)
115 : 0 : status = DftiSetValue(Desc3, DFTI_INPUT_DISTANCE, 1)
116 : : arr(1) = 0
117 : 0 : arr(2) = ldx*ldy
118 : 0 : status = DftiSetValue(Desc3, DFTI_INPUT_STRIDES, arr)
119 : 0 : status = DftiSetValue(Desc3, DFTI_NUMBER_OF_USER_THREADS, nthreads)
120 : : !status = DftiSetValue(Desc2, DFTI_THREAD_LIMIT, 1)
121 : :
122 : : ! Normalize the transform.
123 [ # # # # ]: 0 : if (iscale__ /= 0) then
124 : 0 : status = DftiSetValue(Desc3, DFTI_BACKWARD_SCALE, one/(nx*ny*nz))
125 : : end if
126 : :
127 : 0 : status = DftiCommitDescriptor(Desc3)
128 [ # # # # ]: 0 : DFTI_CHECK(status)
129 : :
130 : : !$OMP PARALLEL DO PRIVATE(sidx,status)
131 [ # # # # ]: 0 : do dat=1,ndat
132 : 0 : sidx = 1 + (dat-1) *ldx*ldy*ldz
133 : 0 : status = DftiComputeBackward(Desc3, ff(sidx:))
134 [ # # # # ]: 0 : DFTI_CHECK(status)
135 : : end do
136 : :
137 : 0 : status = DftiFreeDescriptor(Desc3)
138 [ # # # # ]: 0 : DFTI_CHECK(status)
139 : : end if
140 : :
141 : : case (-1) ! R --> G. Inverse of the algoritm used for G --> R
142 : : !
143 : : ! 1) Transform along z. Note We have to visit the entire augmented x-y plane!
144 [ + + - + : 5736913 : if (ndat==1 .or. nthreads==1 .or. MOD(ndat,nthreads)/=0) then
- - + + -
+ - - ]
145 : : !
146 : : ! Use MKL internal threading if single FFT or ndat is not divisible by nthreads
147 : 5736913 : status = DftiCreateDescriptor(Desc3, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, nz)
148 [ - + - + ]: 5736913 : DFTI_CHECK(status)
149 : :
150 : 5736913 : status = DftiSetValue(Desc3, DFTI_NUMBER_OF_TRANSFORMS, ldx*ldy)
151 : 5736913 : status = DftiSetValue(Desc3, DFTI_INPUT_DISTANCE, 1)
152 : 5736913 : arr(1) = 0
153 : 5736913 : arr(2) = ldx*ldy
154 : 5736913 : status = DftiSetValue(Desc3, DFTI_INPUT_STRIDES, arr)
155 : : !status = DftiSetValue(Desc3, DFTI_THREAD_LIMIT, nthreads)
156 : :
157 : 5736913 : status = DftiCommitDescriptor(Desc3)
158 [ - + - + ]: 5736913 : DFTI_CHECK(status)
159 : :
160 [ + + + + ]: 11473874 : do dat=1,ndat
161 : 5736961 : sidx = 1 + (dat-1) *ldx*ldy*ldz
162 : 5736961 : status = DftiComputeForward(Desc3, ff(sidx:))
163 [ - + - + ]: 11473874 : DFTI_CHECK(status)
164 : : end do
165 : :
166 : 5736913 : status = DftiFreeDescriptor(Desc3)
167 [ - + - + ]: 5736913 : DFTI_CHECK(status)
168 : : else
169 : : ! Split ndat transforms among the threads.
170 : 0 : status = DftiCreateDescriptor(Desc3, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, nz)
171 [ # # # # ]: 0 : DFTI_CHECK(status)
172 : :
173 : 0 : status = DftiSetValue(Desc3, DFTI_NUMBER_OF_TRANSFORMS, ldx*ldy)
174 : 0 : status = DftiSetValue(Desc3, DFTI_INPUT_DISTANCE, 1)
175 : 0 : arr(1) = 0
176 : 0 : arr(2) = ldx*ldy
177 : 0 : status = DftiSetValue(Desc3, DFTI_INPUT_STRIDES, arr)
178 : 0 : status = DftiSetValue(Desc3, DFTI_NUMBER_OF_USER_THREADS, nthreads)
179 : : !status = DftiSetValue(Desc3, DFTI_THREAD_LIMIT, 1)
180 : :
181 : 0 : status = DftiCommitDescriptor(Desc3)
182 [ # # # # ]: 0 : DFTI_CHECK(status)
183 : :
184 : : !$OMP PARALLEL DO PRIVATE(sidx,status)
185 [ # # # # ]: 0 : do dat=1,ndat
186 : 0 : sidx = 1 + (dat-1) *ldx*ldy*ldz
187 : 0 : status = DftiComputeForward(Desc3, ff(sidx:))
188 [ # # # # ]: 0 : DFTI_CHECK(status)
189 : : end do
190 : :
191 : 0 : status = DftiFreeDescriptor(Desc3)
192 [ # # # # ]: 0 : DFTI_CHECK(status)
193 : : end if
194 : : !
195 : : ! 2) Transform along y.
196 : :
197 : : !call mkl_set_num_threads(1)
198 : 5736913 : status = DftiCreateDescriptor(Desc2, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, ny)
199 [ - + - + ]: 5736913 : DFTI_CHECK(status)
200 : :
201 : : !status = DftiSetValue(Desc2, DFTI_THREAD_LIMIT, 1)
202 : 5736913 : status = DftiSetValue(Desc2, DFTI_NUMBER_OF_TRANSFORMS, nx)
203 : 5736913 : status = DftiSetValue(Desc2, DFTI_INPUT_DISTANCE, 1)
204 : 5736913 : arr(1) = 0
205 : 5736913 : arr(2) = ldx
206 : 5736913 : status = DftiSetValue(Desc2, DFTI_INPUT_STRIDES, arr)
207 : :
208 : 5736913 : status = DftiSetValue(Desc2, DFTI_NUMBER_OF_USER_THREADS, nthreads)
209 : :
210 : 5736913 : status = DftiCommitDescriptor(Desc2)
211 [ - + - + ]: 5736913 : DFTI_CHECK(status)
212 : :
213 : : !$OMP PARALLEL DO PRIVATE(dat,zplane,kk,sidx,status)
214 [ + + + + ]: 35492260 : do cnt=0,(zpad%n_zplanes*ndat)-1
215 : 29755347 : dat = 1 + cnt / zpad%n_zplanes
216 : 29755347 : zplane = 1 + MOD(cnt, zpad%n_zplanes)
217 : 29755347 : kk = zpad%zplane(1,zplane)
218 : 29755347 : sidx = 1 + ldx*ldy*(kk-1) + (dat-1) *ldx*ldy*ldz
219 : 29755347 : status = DftiComputeForward(Desc2, ff(sidx:))
220 [ - + - + ]: 35492260 : DFTI_CHECK(status)
221 : : end do
222 : :
223 : 5736913 : status = DftiFreeDescriptor(Desc2)
224 [ - + - + ]: 5736913 : DFTI_CHECK(status)
225 : : !
226 : : ! 3) Transform along x.
227 : : !
228 : 5736913 : status = DftiCreateDescriptor(Desc1, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, nx)
229 [ - + - + ]: 5736913 : DFTI_CHECK(status)
230 : :
231 : : !status = DftiSetValue(Desc1, DFTI_THREAD_LIMIT, 1)
232 : 5736913 : status = DftiSetValue(Desc1, DFTI_NUMBER_OF_USER_THREADS, nthreads)
233 : :
234 : : ! Normalize the transform.
235 [ + - + - ]: 5736913 : if (iscale__ /= 0) then
236 : 5736913 : status = DftiSetValue(Desc1, DFTI_FORWARD_SCALE, one/(nx*ny*nz))
237 : : end if
238 : :
239 : 5736913 : status = DftiCommitDescriptor(Desc1)
240 [ - + - + ]: 5736913 : DFTI_CHECK(status)
241 : :
242 : : !$OMP PARALLEL DO PRIVATE(dat,line,jj,kk,sidx,status)
243 [ + + + + ]: 117622298 : do cnt=0,(zpad%nlinex*ndat)-1
244 : 111885385 : dat = 1 + cnt / zpad%nlinex
245 : 111885385 : line = 1 + MOD(cnt, zpad%nlinex)
246 : 111885385 : jj = zpad%linex2ifft_yz(1,line)
247 : 111885385 : kk = zpad%linex2ifft_yz(2,line)
248 : 111885385 : sidx = 1+ (jj-1)*ldx + (kk-1)*ldx*ldy + (dat-1) * ldx*ldy*ldz
249 : 111885385 : status = DftiComputeForward(Desc1, ff(sidx:)) ! Pass the base address.
250 [ - + - + ]: 117622298 : DFTI_CHECK(status)
251 : : end do
252 : :
253 : 5736913 : status = DftiFreeDescriptor(Desc1)
254 [ - + - + ]: 5736913 : DFTI_CHECK(status)
255 : :
256 : : case default
257 [ + + - + : 5739117 : ABI_BUG("Wrong isign")
+ - ]
258 : : end select
259 : :
260 : 5739117 : call zpad%free()
261 : :
262 : : #else
263 : : ABI_ERROR("FFT_DFTI support not activated")
264 : : ABI_UNUSED((/nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,isign,nthreads/))
265 : : ABI_UNUSED(gbound(1,1))
266 : : ABI_UNUSED(ff(1))
267 : : #endif
|
