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 5736024 : type(zpad_t) :: zpad
5 : type(DFTI_DESCRIPTOR),pointer :: Desc1,Desc2,Desc3
6 :
7 5738054 : iscale__ = merge(1, 0, isign == -1); if (present(iscale)) iscale__ = iscale
8 :
9 5736024 : nthreads = xomp_get_num_threads(open_parallel=.TRUE.)
10 :
11 5736024 : call zpad%init(nx, ny, nz, ldx, ldy, ldz, mgfft, gbound)
12 :
13 2030 : 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 2030 : status = DftiCreateDescriptor(Desc1, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, nx)
19 2030 : 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 2030 : status = DftiSetValue(Desc1, DFTI_NUMBER_OF_USER_THREADS, nthreads)
26 :
27 2030 : status = DftiCommitDescriptor(Desc1)
28 2030 : DFTI_CHECK(status)
29 :
30 : !$OMP PARALLEL DO PRIVATE(dat,line,jj,kk,sidx,status)
31 400722 : do cnt=0,(zpad%nlinex*ndat)-1
32 398692 : dat = 1 + cnt / zpad%nlinex
33 398692 : line = 1 + MOD(cnt, zpad%nlinex)
34 398692 : jj = zpad%linex2ifft_yz(1,line)
35 398692 : kk = zpad%linex2ifft_yz(2,line)
36 398692 : sidx = 1+ (jj-1)*ldx + (kk-1)*ldx*ldy + (dat-1) * ldx*ldy*ldz
37 398692 : status = DftiComputeBackward(Desc1, ff(sidx:)) ! Pass the base address.
38 400722 : DFTI_CHECK(status)
39 : end do
40 :
41 2030 : status = DftiFreeDescriptor(Desc1)
42 2030 : 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 2030 : status = DftiCreateDescriptor(Desc2, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, ny)
48 2030 : DFTI_CHECK(status)
49 :
50 : !status = DftiSetValue(Desc2, DFTI_THREAD_LIMIT, 1)
51 2030 : status = DftiSetValue(Desc2, DFTI_NUMBER_OF_TRANSFORMS, nx)
52 2030 : status = DftiSetValue(Desc2, DFTI_INPUT_DISTANCE, 1)
53 2030 : arr(1) = 0
54 2030 : arr(2) = ldx
55 2030 : status = DftiSetValue(Desc2, DFTI_INPUT_STRIDES, arr)
56 2030 : status = DftiSetValue(Desc2, DFTI_NUMBER_OF_USER_THREADS, nthreads)
57 :
58 2030 : status = DftiCommitDescriptor(Desc2)
59 2030 : DFTI_CHECK(status)
60 :
61 : !$OMP PARALLEL DO PRIVATE(dat,zplane,kk,sidx,status)
62 31862 : do cnt=0,(zpad%n_zplanes*ndat)-1
63 29832 : dat = 1 + cnt / zpad%n_zplanes
64 29832 : zplane = 1 + MOD(cnt, zpad%n_zplanes)
65 29832 : kk = zpad%zplane(1,zplane)
66 29832 : sidx = 1 + (kk-1)*ldx*ldy + (dat-1) *ldx*ldy*ldz
67 29832 : status = DftiComputeBackward(Desc2, ff(sidx:))
68 31862 : DFTI_CHECK(status)
69 : end do
70 :
71 2030 : status = DftiFreeDescriptor(Desc2)
72 2030 : 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 2030 : 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 2030 : status = DftiCreateDescriptor(Desc3, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, nz)
83 2030 : DFTI_CHECK(status)
84 :
85 2030 : status = DftiSetValue(Desc3, DFTI_NUMBER_OF_TRANSFORMS, ldx*ldy)
86 2030 : status = DftiSetValue(Desc3, DFTI_INPUT_DISTANCE, 1)
87 : arr(1) = 0
88 2030 : arr(2) = ldx*ldy
89 2030 : status = DftiSetValue(Desc3, DFTI_INPUT_STRIDES, arr)
90 :
91 : ! Normalize the transform.
92 2030 : if (iscale__ /= 0) then
93 0 : status = DftiSetValue(Desc3, DFTI_BACKWARD_SCALE, one/(nx*ny*nz))
94 : end if
95 :
96 2030 : status = DftiCommitDescriptor(Desc3)
97 2030 : DFTI_CHECK(status)
98 :
99 4102 : do dat=1,ndat
100 2072 : sidx = 1 + (dat-1) *ldx*ldy*ldz
101 2072 : status = DftiComputeBackward(Desc3, ff(sidx:))
102 4102 : DFTI_CHECK(status)
103 : end do
104 :
105 2030 : status = DftiFreeDescriptor(Desc3)
106 2030 : 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 5733994 : 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 5733994 : status = DftiCreateDescriptor(Desc3, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, nz)
148 5733994 : DFTI_CHECK(status)
149 :
150 5733994 : status = DftiSetValue(Desc3, DFTI_NUMBER_OF_TRANSFORMS, ldx*ldy)
151 5733994 : status = DftiSetValue(Desc3, DFTI_INPUT_DISTANCE, 1)
152 5733994 : arr(1) = 0
153 5733994 : arr(2) = ldx*ldy
154 5733994 : status = DftiSetValue(Desc3, DFTI_INPUT_STRIDES, arr)
155 : !status = DftiSetValue(Desc3, DFTI_THREAD_LIMIT, nthreads)
156 :
157 5733994 : status = DftiCommitDescriptor(Desc3)
158 5733994 : DFTI_CHECK(status)
159 :
160 11468036 : do dat=1,ndat
161 5734042 : sidx = 1 + (dat-1) *ldx*ldy*ldz
162 5734042 : status = DftiComputeForward(Desc3, ff(sidx:))
163 11468036 : DFTI_CHECK(status)
164 : end do
165 :
166 5733994 : status = DftiFreeDescriptor(Desc3)
167 5733994 : 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 5733994 : status = DftiCreateDescriptor(Desc2, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, ny)
199 5733994 : DFTI_CHECK(status)
200 :
201 : !status = DftiSetValue(Desc2, DFTI_THREAD_LIMIT, 1)
202 5733994 : status = DftiSetValue(Desc2, DFTI_NUMBER_OF_TRANSFORMS, nx)
203 5733994 : status = DftiSetValue(Desc2, DFTI_INPUT_DISTANCE, 1)
204 5733994 : arr(1) = 0
205 5733994 : arr(2) = ldx
206 5733994 : status = DftiSetValue(Desc2, DFTI_INPUT_STRIDES, arr)
207 :
208 5733994 : status = DftiSetValue(Desc2, DFTI_NUMBER_OF_USER_THREADS, nthreads)
209 :
210 5733994 : status = DftiCommitDescriptor(Desc2)
211 5733994 : DFTI_CHECK(status)
212 :
213 : !$OMP PARALLEL DO PRIVATE(dat,zplane,kk,sidx,status)
214 35475544 : do cnt=0,(zpad%n_zplanes*ndat)-1
215 29741550 : dat = 1 + cnt / zpad%n_zplanes
216 29741550 : zplane = 1 + MOD(cnt, zpad%n_zplanes)
217 29741550 : kk = zpad%zplane(1,zplane)
218 29741550 : sidx = 1 + ldx*ldy*(kk-1) + (dat-1) *ldx*ldy*ldz
219 29741550 : status = DftiComputeForward(Desc2, ff(sidx:))
220 35475544 : DFTI_CHECK(status)
221 : end do
222 :
223 5733994 : status = DftiFreeDescriptor(Desc2)
224 5733994 : DFTI_CHECK(status)
225 : !
226 : ! 3) Transform along x.
227 : !
228 5733994 : status = DftiCreateDescriptor(Desc1, DEV_DFTI_PRECISION, DFTI_COMPLEX, 1, nx)
229 5733994 : DFTI_CHECK(status)
230 :
231 : !status = DftiSetValue(Desc1, DFTI_THREAD_LIMIT, 1)
232 5733994 : status = DftiSetValue(Desc1, DFTI_NUMBER_OF_USER_THREADS, nthreads)
233 :
234 : ! Normalize the transform.
235 5733994 : if (iscale__ /= 0) then
236 5733994 : status = DftiSetValue(Desc1, DFTI_FORWARD_SCALE, one/(nx*ny*nz))
237 : end if
238 :
239 5733994 : status = DftiCommitDescriptor(Desc1)
240 5733994 : DFTI_CHECK(status)
241 :
242 : !$OMP PARALLEL DO PRIVATE(dat,line,jj,kk,sidx,status)
243 117574079 : do cnt=0,(zpad%nlinex*ndat)-1
244 111840085 : dat = 1 + cnt / zpad%nlinex
245 111840085 : line = 1 + MOD(cnt, zpad%nlinex)
246 111840085 : jj = zpad%linex2ifft_yz(1,line)
247 111840085 : kk = zpad%linex2ifft_yz(2,line)
248 111840085 : sidx = 1+ (jj-1)*ldx + (kk-1)*ldx*ldy + (dat-1) * ldx*ldy*ldz
249 111840085 : status = DftiComputeForward(Desc1, ff(sidx:)) ! Pass the base address.
250 117574079 : DFTI_CHECK(status)
251 : end do
252 :
253 5733994 : status = DftiFreeDescriptor(Desc1)
254 5733994 : DFTI_CHECK(status)
255 :
256 : case default
257 5736024 : ABI_BUG("Wrong isign")
258 : end select
259 :
260 5736024 : 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
|
