Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This file contains procedure implementations of [pm_fftpack](@ref pm_fftpack).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, Wednesday 12:20 PM, September 22, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
28 : #if setFactorFFT_ENABLED && DCA_ENABLED
29 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 :
31 2 : allocate(coef(size(data, 1, IK)))
32 12 : factor = getFactorFFT(data, coef)
33 :
34 : !%%%%%%%%%%%%%%%%%%%
35 : #elif setFactorFFT_ENABLED
36 : !%%%%%%%%%%%%%%%%%%%
37 :
38 : integer(IK) :: lenData, i, j, lenFactor, nl, nq, nr, ntry
39 : #if CK_ENABLED
40 : integer(IK), parameter :: NUM_TRYH(4) = [integer(IK) :: 3, 4, 2, 5]
41 : #elif RK_ENABLED
42 : integer(IK), parameter :: NUM_TRYH(4) = [integer(IK) :: 4, 2, 3, 5]
43 : #else
44 : #error "Unrecognized interface."
45 : #endif
46 : #if DCX_ENABLED
47 : integer(IK) :: ido, ii, is, k1, l1, l2, ld, itmp
48 : real(TKC), parameter :: TWO_TIMES_PI = 2.0_TKC * acos( - 1.0_TKC)
49 : real(TKC) :: arg, argh, argld, fi
50 : #endif
51 60074 : CHECK_ASSERTION(__LINE__, size(data, kind = IK) > 1_IK, SK_"@getFactorFFT(): The condition `size(data) > 1` must hold. size(data) = "//getStr([size(data, kind = IK)]))
52 : lenData = size(data, kind = IK)
53 30037 : allocate(factor(15))
54 : lenFactor = 0_IK
55 : nl = lenData
56 : j = 0_IK
57 : loop1: do
58 426400 : j = j + 1_IK
59 426400 : if (j <= 4_IK) then
60 114192 : ntry = NUM_TRYH(j)
61 : else
62 312208 : ntry = ntry + 2_IK
63 : end if
64 : loop2: do
65 458648 : nq = nl / ntry
66 : nr = nl - ntry * nq
67 458648 : if (nr /= 0_IK) cycle loop1
68 62285 : lenFactor = lenFactor + 1_IK
69 62285 : if (lenFactor > size(factor, kind = IK)) call setResized(factor)
70 62285 : factor(lenFactor) = ntry
71 : nl = nq
72 62285 : if (ntry == 2_IK) then
73 8080 : if (lenFactor /= 1_IK) then
74 8198 : do i = lenFactor, 2_IK, -1_IK
75 8198 : factor(i) = factor(i - 1_IK)
76 : end do
77 3367 : factor(1_IK) = 2_IK
78 : end if
79 : end if
80 62285 : if (nl /= 1_IK) cycle loop2
81 428611 : exit loop1
82 : end do loop2
83 : end do loop1
84 184644 : factor = factor(1:lenFactor)
85 : !factor(2_IK) = lenFactor
86 : !factor(1_IK) = lenData
87 : #if DCX_ENABLED
88 15281 : argh = TWO_TIMES_PI / real(lenData, TKC)
89 : #if CK_ENABLED
90 : i = 1_IK
91 : l1 = 1_IK
92 : itmp = lenFactor
93 46837 : do k1 = 1_IK, itmp
94 : ld = 0_IK
95 31556 : l2 = l1 * factor(k1)
96 31556 : ido = lenData / l2
97 31556 : itmp = ido + ido + 2_IK
98 403857 : do j = 1_IK, factor(k1) - 1_IK
99 : is = i
100 372301 : coef(i) = (1._TKC, 0._TKC)
101 190675 : fi = 0._TKC
102 372301 : ld = ld + l1
103 372301 : argld = ld * argh
104 372301 : do ii = 4_IK, itmp, 2_IK
105 822965 : i = i + 1_IK
106 822965 : fi = fi + 1._TKC
107 822965 : arg = fi * argld
108 822965 : coef(i) = cmplx(cos(arg), sin(arg), TKC)
109 : end do
110 403857 : if (factor(k1) > 5_IK) coef(is) = coef(i)
111 : end do
112 15281 : l1 = l2
113 : end do
114 : #elif RK_ENABLED
115 14756 : if (lenFactor == 1_IK) return
116 : itmp = lenFactor - 1_IK
117 : is = 0_IK
118 : l1 = 1_IK
119 26337 : do k1 = 1_IK, itmp
120 : ld = 0_IK
121 15973 : l2 = l1 * factor(k1)
122 15973 : ido = lenData / l2
123 52798 : do j = 1_IK, factor(k1) - 1_IK
124 : i = is
125 18710 : fi = 0._TKC
126 36825 : ld = ld + l1
127 36825 : argld = ld * argh
128 36825 : do ii = 3_IK, ido, 2_IK
129 220593 : i = i + 2_IK
130 220593 : fi = fi + 1._TKC
131 220593 : arg = fi * argld
132 220593 : coef(i) = sin(arg)
133 220593 : coef(i - 1_IK) = cos(arg)
134 : end do
135 52798 : is = is + ido
136 : end do
137 10364 : l1 = l2
138 : end do
139 : #else
140 : #error "Unrecognized interface."
141 : #endif
142 : #elif !DXX_ENABLED
143 : #error "Unrecognized interface."
144 : #endif
145 :
146 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147 : #elif setFFTF_ENABLED && CK_ENABLED
148 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
149 :
150 : integer(IK) :: lenData, idl1, ido, ip, ix1, ix2, ix3, ix4, k1, l1, l2, na, nac, lenFactor
151 26476 : lenFactor = size(factor, kind = IK)
152 26476 : lenData = size(data, kind = IK)
153 79428 : CHECK_ASSERTION(__LINE__, lenData == size(coef, kind = IK), SK_"@setFFTF(): The condition `size(data) == size(coef)` must hold. size(data), size(coef) = "//getStr([lenData, size(coef, kind = IK)]))
154 79428 : CHECK_ASSERTION(__LINE__, lenData == size(work, kind = IK), SK_"@setFFTF(): The condition `size(data) == size(work)` must hold. size(data), size(work) = "//getStr([lenData, size(work, kind = IK)]))
155 26476 : inwork = logical(1_IK < lenData, LK)
156 26476 : if (.not. inwork) return
157 : na = 0_IK
158 26476 : l1 = 1_IK
159 26476 : ix1 = 1_IK
160 81503 : do k1 = 1_IK, lenFactor
161 55027 : ip = factor(k1)
162 55027 : l2 = ip * l1
163 55027 : ido = lenData / l2
164 55027 : idl1 = ido * l1
165 55027 : if (ip == 4_IK) then
166 7794 : ix2 = ix1 + ido
167 7794 : ix3 = ix2 + ido
168 7794 : if (na /= 0_IK) then
169 3737 : call passf4(ido, l1, ix1, ix2, ix3, coef, work, data)
170 : else
171 4057 : call passf4(ido, l1, ix1, ix2, ix3, coef, data, work)
172 : end if
173 7794 : na = 1_IK - na
174 47233 : elseif (ip == 2_IK) then
175 7646 : if (na /= 0_IK) then
176 0 : call passf2(ido, l1, ix1, coef, work, data)
177 : else
178 7646 : call passf2(ido, l1, ix1, coef, data, work)
179 : end if
180 7646 : na = 1_IK - na
181 39587 : elseif (ip == 3_IK) then
182 12442 : ix2 = ix1 + ido
183 12442 : if (na /= 0_IK) then
184 4894 : call passf3(ido, l1, ix1, ix2, coef, work, data)
185 : else
186 7548 : call passf3(ido, l1, ix1, ix2, coef, data, work)
187 : end if
188 12442 : na = 1_IK - na
189 27145 : elseif (ip /= 5_IK) then
190 20621 : if (na /= 0_IK) then
191 8732 : call passf(nac, ido, ip, l1, idl1, ix1, coef, work, work, work, data, data)
192 : else
193 11889 : call passf(nac, ido, ip, l1, idl1, ix1, coef, data, data, data, work, work)
194 : end if
195 20621 : if (nac /= 0_IK) na = 1_IK - na
196 : else
197 6524 : ix2 = ix1 + ido
198 6524 : ix3 = ix2 + ido
199 6524 : ix4 = ix3 + ido
200 6524 : if (na /= 0_IK) then
201 2132 : call passf5(ido, l1, ix1, ix2, ix3, ix4, coef, work, data)
202 : else
203 4392 : call passf5(ido, l1, ix1, ix2, ix3, ix4, coef, data, work)
204 : end if
205 6524 : na = 1_IK - na
206 : end if
207 55027 : l1 = l2
208 81503 : ix1 = ix1 + (ip - 1_IK) * ido
209 : end do
210 26476 : inwork = logical(na /= 0_IK, LK)
211 : !if (inwork) return
212 : !data(1:lenData) = work(1:lenData) ! \todo these copies can be avoided.
213 :
214 : contains
215 :
216 7646 : pure subroutine passf2(ido, l1, ix1, coef, data, work)
217 : integer(IK) , intent(in) :: ido, l1, ix1
218 : complex(TKC), intent(in) , contiguous :: coef(2:)
219 : complex(TKC), intent(in) :: data(ido, 2, l1)
220 : complex(TKC), intent(out) :: work(ido, l1, 2)
221 : integer(IK) :: i, k
222 7646 : if (ido > 1_IK) then
223 15292 : do k = 1_IK, l1
224 237487 : do i = 1_IK, ido
225 222195 : work(i, k, 1) = data(i, 1, k) + data(i, 2, k)
226 229841 : work(i, k, 2) = (data(i, 1, k) - data(i, 2, k)) * conjg(coef(ix1 + i))
227 : end do
228 : end do
229 : else
230 0 : do k = 1_IK, l1
231 0 : work(1, k, 1) = data(1, 1, k) + data(1, 2, k)
232 0 : work(1, k, 2) = data(1, 1, k) - data(1, 2, k)
233 : end do
234 : end if
235 7646 : end subroutine
236 :
237 12442 : pure subroutine passf3(ido, l1, ix1, ix2, coef, data, work)
238 : integer(IK) , intent(in) :: ido, l1, ix1, ix2
239 : complex(TKC), intent(in) , contiguous :: coef(2:)
240 : complex(TKC), intent(in) :: data(ido, 3, l1)
241 : complex(TKC), intent(out) :: work(ido, l1, 3)
242 : complex(TKC), parameter :: TAU = -cmplx(0.5_TKC, sqrt(3._TKC) / 2._TKC, TKC)
243 : complex(TKC) :: c2, c3, d2, d3, t2
244 : integer(IK) :: i, k
245 12442 : if (ido /= 1_IK) then
246 35925 : do k = 1_IK, l1
247 254678 : do i = 1_IK, ido
248 218753 : t2 = data(i, 2, k) + data(i, 3, k)
249 218753 : work(i, k, 1) = data(i, 1, k) + t2
250 218753 : c2 = data(i, 1, k) + TAU%re * t2
251 218753 : c3 = (data(i, 2, k) - data(i, 3, k)) * TAU%im
252 218753 : d2%re = c2%re - c3%im
253 218753 : d2%im = c2%im + c3%re
254 218753 : d3%re = c2%re + c3%im
255 218753 : d3%im = c2%im - c3%re
256 218753 : work(i, k, 2)%im = coef(ix1 + i)%re * d2%im - coef(ix1 + i)%im * d2%re
257 218753 : work(i, k, 2)%re = coef(ix1 + i)%re * d2%re + coef(ix1 + i)%im * d2%im
258 218753 : work(i, k, 3)%im = coef(ix2 + i)%re * d3%im - coef(ix2 + i)%im * d3%re
259 243657 : work(i, k, 3)%re = coef(ix2 + i)%re * d3%re + coef(ix2 + i)%im * d3%im
260 : end do
261 : end do
262 : else
263 20256 : do k = 1_IK, l1
264 18835 : t2 = data(1, 2, k) + data(1, 3, k)
265 18835 : c2 = data(1, 1, k) + TAU%re * t2
266 18835 : work(1, k, 1) = data(1, 1, k) + t2
267 18835 : c3 = (data(1, 2, k) - data(1, 3, k)) * TAU%im
268 18835 : work(1, k, 2)%re = c2%re - c3%im
269 18835 : work(1, k, 2)%im = c2%im + c3%re
270 18835 : work(1, k, 3)%re = c2%re + c3%im
271 20256 : work(1, k, 3)%im = c2%im - c3%re
272 : end do
273 : end if
274 12442 : end subroutine
275 :
276 7794 : pure subroutine passf4(ido, l1, ix1, ix2, ix3, coef, data, work)
277 : integer(IK) , intent(in) :: ido, l1, ix1, ix2, ix3
278 : complex(TKC), intent(in) , contiguous :: coef(2:)
279 : complex(TKC), intent(in) :: data(ido, 4, l1)
280 : complex(TKC), intent(out) :: work(ido, l1, 4)
281 : complex(TKC) :: c2, c3, c4, t1, t2, t3, t4
282 : integer(IK) :: i, k
283 7794 : if (ido /= 1_IK) then
284 17873 : do k = 1_IK, l1
285 109212 : do i = 1_IK, ido
286 91339 : t1 = data(i, 1, k) - data(i, 3, k)
287 91339 : t2 = data(i, 1, k) + data(i, 3, k)
288 91339 : t3 = data(i, 2, k) + data(i, 4, k)
289 91339 : t4%re = data(i, 2, k)%im - data(i, 4, k)%im
290 91339 : t4%im = data(i, 4, k)%re - data(i, 2, k)%re
291 91339 : work(i, k, 1) = t2 + t3
292 91339 : c3 = t2 - t3
293 91339 : c2 = t1 + t4
294 91339 : c4 = t1 - t4
295 91339 : work(i, k, 2)%re = coef(ix1 + i)%re * c2%re + coef(ix1 + i)%im * c2%im
296 91339 : work(i, k, 2)%im = coef(ix1 + i)%re * c2%im - coef(ix1 + i)%im * c2%re
297 91339 : work(i, k, 3)%re = coef(ix2 + i)%re * c3%re + coef(ix2 + i)%im * c3%im
298 91339 : work(i, k, 3)%im = coef(ix2 + i)%re * c3%im - coef(ix2 + i)%im * c3%re
299 91339 : work(i, k, 4)%re = coef(ix3 + i)%re * c4%re + coef(ix3 + i)%im * c4%im
300 103201 : work(i, k, 4)%im = coef(ix3 + i)%re * c4%im - coef(ix3 + i)%im * c4%re
301 : end do
302 : end do
303 : else
304 23799 : do k = 1_IK, l1
305 22016 : t1 = data(1, 1, k) - data(1, 3, k)
306 22016 : t2 = data(1, 1, k) + data(1, 3, k)
307 22016 : t3 = data(1, 2, k) + data(1, 4, k)
308 22016 : t4%re = data(1, 2, k)%im - data(1, 4, k)%im
309 22016 : t4%im = data(1, 4, k)%re - data(1, 2, k)%re
310 22016 : work(1, k, 1) = t2 + t3
311 22016 : work(1, k, 2) = t1 + t4
312 22016 : work(1, k, 3) = t2 - t3
313 23799 : work(1, k, 4) = t1 - t4
314 : end do
315 : end if
316 7794 : end subroutine
317 :
318 6524 : PURE subroutine passf5(ido, l1, ix1, ix2, ix3, ix4, coef, data, work)
319 : integer(IK) , intent(in) :: ido, l1, ix1, ix2, ix3, ix4
320 : complex(TKC), intent(in) , contiguous :: coef(2:)
321 : complex(TKC), intent(in) :: data(ido, 5, l1)
322 : complex(TKC), intent(out) :: work(ido, l1, 5)
323 : real(TKC) , parameter :: PI = acos( - 1._TKC)
324 : complex(TKC), parameter :: T11 = cmplx(cos(2._TKC * PI / 5._TKC), -sin(2._TKC * PI / 5._TKC), TKC)
325 : complex(TKC), parameter :: T12 = cmplx(cos(4._TKC * PI / 5._TKC), -sin(4._TKC * PI / 5._TKC), TKC)
326 : complex(TKC) :: c2, c3, c4, c5, d2, d3, d4, d5, t2, t3, t4, t5
327 : integer(IK) :: i, k
328 6524 : if (ido /= 1_IK) then
329 7403 : do k = 1_IK, l1
330 43924 : do i = 1_IK, ido
331 36521 : t2 = data(i, 2, k) + data(i, 5, k)
332 36521 : t3 = data(i, 3, k) + data(i, 4, k)
333 36521 : t4 = data(i, 3, k) - data(i, 4, k)
334 36521 : t5 = data(i, 2, k) - data(i, 5, k)
335 36521 : work(i, k, 1) = data(i, 1, k) + t2 + t3
336 36521 : c2%re = data(i, 1, k)%re + T11%re * t2%re + T12%re * t3%re
337 36521 : c2%im = data(i, 1, k)%im + T11%re * t2%im + T12%re * t3%im
338 36521 : c3%re = data(i, 1, k)%re + T12%re * t2%re + T11%re * t3%re
339 36521 : c3%im = data(i, 1, k)%im + T12%re * t2%im + T11%re * t3%im
340 36521 : c4%re = T12%im * t5%re - T11%im * t4%re
341 36521 : c4%im = T12%im * t5%im - T11%im * t4%im
342 36521 : c5%re = T11%im * t5%re + T12%im * t4%re
343 36521 : c5%im = T11%im * t5%im + T12%im * t4%im
344 36521 : d2%re = c2%re - c5%im
345 36521 : d2%im = c2%im + c5%re
346 36521 : d3%re = c3%re - c4%im
347 36521 : d3%im = c3%im + c4%re
348 36521 : d4%re = c3%re + c4%im
349 36521 : d4%im = c3%im - c4%re
350 36521 : d5%re = c2%re + c5%im
351 36521 : d5%im = c2%im - c5%re
352 36521 : work(i, k, 2)%re = coef(ix1 + i)%re * d2%re + coef(ix1 + i)%im * d2%im
353 36521 : work(i, k, 2)%im = coef(ix1 + i)%re * d2%im - coef(ix1 + i)%im * d2%re
354 36521 : work(i, k, 3)%re = coef(ix2 + i)%re * d3%re + coef(ix2 + i)%im * d3%im
355 36521 : work(i, k, 3)%im = coef(ix2 + i)%re * d3%im - coef(ix2 + i)%im * d3%re
356 36521 : work(i, k, 4)%re = coef(ix3 + i)%re * d4%re + coef(ix3 + i)%im * d4%im
357 36521 : work(i, k, 4)%im = coef(ix3 + i)%re * d4%im - coef(ix3 + i)%im * d4%re
358 36521 : work(i, k, 5)%re = coef(ix4 + i)%re * d5%re + coef(ix4 + i)%im * d5%im
359 40904 : work(i, k, 5)%im = coef(ix4 + i)%re * d5%im - coef(ix4 + i)%im * d5%re
360 : end do
361 : end do
362 : else
363 38860 : do k = 1_IK, l1
364 35356 : t2 = data(1, 2, k) + data(1, 5, k)
365 35356 : t3 = data(1, 3, k) + data(1, 4, k)
366 35356 : t4 = data(1, 3, k) - data(1, 4, k)
367 35356 : t5 = data(1, 2, k) - data(1, 5, k)
368 35356 : work(1, k, 1) = data(1, 1, k) + t2 + t3
369 35356 : c2%re = data(1, 1, k)%re + T11%re * t2%re + T12%re * t3%re
370 35356 : c2%im = data(1, 1, k)%im + T11%re * t2%im + T12%re * t3%im
371 35356 : c3%re = data(1, 1, k)%re + T12%re * t2%re + T11%re * t3%re
372 35356 : c3%im = data(1, 1, k)%im + T12%re * t2%im + T11%re * t3%im
373 35356 : c5%re = T11%im * t5%re + T12%im * t4%re
374 35356 : c5%im = T11%im * t5%im + T12%im * t4%im
375 35356 : c4%re = T12%im * t5%re - T11%im * t4%re
376 35356 : c4%im = T12%im * t5%im - T11%im * t4%im
377 35356 : work(1, k, 2)%re = c2%re - c5%im
378 35356 : work(1, k, 2)%im = c2%im + c5%re
379 35356 : work(1, k, 3)%re = c3%re - c4%im
380 35356 : work(1, k, 3)%im = c3%im + c4%re
381 35356 : work(1, k, 4)%re = c3%re + c4%im
382 35356 : work(1, k, 4)%im = c3%im - c4%re
383 35356 : work(1, k, 5)%re = c2%re + c5%im
384 38860 : work(1, k, 5)%im = c2%im - c5%re
385 : !work(1, k, 2) = c2%re - c5%im
386 : !work(2, k, 2) = c2%im + c5%re
387 : !work(1, k, 3) = c3%re - c4%im
388 : !work(2, k, 3) = c3%im + c4%re
389 : !work(1, k, 4) = c3%re + c4%im
390 : !work(2, k, 4) = c3%im - c4%re
391 : !work(1, k, 5) = c2%re + c5%im
392 : !work(2, k, 5) = c2%im - c5%re
393 : end do
394 : end if
395 6524 : end subroutine
396 :
397 20621 : pure subroutine passf(nac, ido, ip, l1, idl1, ix1, coef, data, C1, C2, work, Ch2)
398 : integer(IK) , intent(out) :: nac
399 : integer(IK) , intent(in) :: ido, ip, l1, idl1, ix1
400 : complex(TKC), intent(in) , contiguous :: coef(2:)
401 : complex(TKC), intent(inout) :: C1(ido,l1,ip), C2(idl1,ip)
402 : complex(TKC), intent(inout) :: data(ido,ip,l1), work(ido,l1,ip), Ch2(idl1,ip)
403 : integer(IK) :: i, idij, idj, idl, idlj, idp, jk, inc, ipp2, ipph, j, jc, k, l, lc !, nt, idot
404 : !idot = ido / 2_IK
405 : !nt = ip * idl1 what in the world is this doing here?
406 20621 : idp = ip * ido
407 20621 : ipp2 = ip + 2_IK
408 20621 : ipph = (ip + 1_IK) / 2_IK
409 20621 : if (2_IK * ido < l1) then
410 66912 : do j = 2_IK, ipph
411 57307 : jc = ipp2 - j
412 124219 : do i = 1_IK, ido
413 387786 : do k = 1_IK, l1
414 273172 : work(i, k, j) = data(i, j, k) + data(i, jc, k)
415 330479 : work(i, k, jc) = data(i, j, k) - data(i, jc, k)
416 : end do
417 : end do
418 : end do
419 19210 : do i = 1_IK, ido
420 72058 : do k = 1_IK, l1
421 62453 : work(i, k, 1) = data(i, 1, k)
422 : end do
423 : end do
424 : else
425 236845 : do j = 2_IK, ipph
426 225829 : jc = ipp2 - j
427 505283 : do k = 1_IK, l1
428 784935 : do i = 1_IK, ido
429 290668 : work(i, k, j) = data(i, j, k) + data(i, jc, k)
430 559106 : work(i, k, jc) = data(i, j, k) - data(i, jc, k)
431 : end do
432 : end do
433 : end do
434 24980 : do k = 1_IK, l1
435 46354 : do i = 1, ido
436 35338 : work(i, k, 1) = data(i, 1, k)
437 : end do
438 : end do
439 : end if
440 : inc = 0_IK
441 20621 : idl = 1_IK - ido
442 303757 : do l = 2_IK, ipph
443 283136 : lc = ipp2 - l
444 283136 : idl = idl + ido
445 846976 : do jk = 1_IK, idl1
446 563840 : C2(jk, l) = Ch2(jk, 1) + coef(ix1 + idl)%re * Ch2(jk, 2)
447 846976 : C2(jk, lc) = -coef(ix1 + idl)%im * Ch2(jk, ip)
448 : end do
449 : idlj = idl
450 283136 : inc = inc + ido
451 6591861 : do j = 3_IK, ipph
452 6288104 : jc = ipp2 - j
453 6288104 : idlj = idlj + inc
454 6288104 : if (idlj > idp) idlj = idlj - idp
455 14663082 : do jk = 1_IK, idl1
456 8091842 : C2(jk, l) = C2(jk, l) + coef(ix1 + idlj)%re * Ch2(jk, j)
457 14379946 : C2(jk, lc) = C2(jk, lc) - coef(ix1 + idlj)%im * Ch2(jk, jc)
458 : end do
459 : end do
460 : end do
461 303757 : do j = 2_IK, ipph
462 867597 : do jk = 1_IK, idl1
463 846976 : Ch2(jk, 1) = Ch2(jk, 1) + Ch2(jk, j)
464 : end do
465 : end do
466 303757 : do j = 2_IK, ipph
467 283136 : jc = ipp2 - j
468 867597 : do jk = 1_IK, idl1
469 563840 : Ch2(jk, j)%re = C2(jk, j)%re - C2(jk, jc)%im
470 563840 : Ch2(jk, j)%im = C2(jk, j)%im + C2(jk, jc)%re
471 563840 : Ch2(jk, jc)%re = C2(jk, j)%re + C2(jk, jc)%im
472 846976 : Ch2(jk, jc)%im = C2(jk, j)%im - C2(jk, jc)%re
473 : end do
474 : end do
475 20621 : nac = 1_IK
476 20621 : if (ido == 1_IK) return
477 853 : nac = 0_IK
478 9116 : do jk = 1_IK, idl1
479 9116 : C2(jk, 1) = Ch2(jk, 1)
480 : end do
481 5971 : do j = 2_IK, ip
482 11089 : do k = 1_IK, l1
483 10236 : C1(1, k, j) = work(1, k, j)
484 : end do
485 : end do
486 853 : if (ido > l1) then
487 : idj = 1_IK - ido
488 5971 : do j = 2_IK, ip
489 5118 : idj = idj + ido
490 11089 : do k = 1_IK, l1
491 : idij = idj
492 54696 : do i = 2_IK, ido
493 44460 : idij = idij + 1_IK
494 44460 : C1(i, k, j)%re = coef(ix1 + idij)%re * work(i, k, j)%re + coef(ix1 + idij)%im * work(i, k, j)%im
495 49578 : C1(i, k, j)%im = coef(ix1 + idij)%re * work(i, k, j)%im - coef(ix1 + idij)%im * work(i, k, j)%re
496 : end do
497 : end do
498 : end do
499 : else
500 : idij = 0_IK
501 0 : do j = 2_IK, ip
502 0 : idij = idij + 1_IK
503 0 : do i = 2_IK, ido
504 0 : idij = idij + 2_IK
505 0 : do k = 1_IK, l1
506 0 : C1(i, k, j)%re = coef(ix1 + idij)%re * work(i, k, j)%re + coef(ix1 + idij)%im * work(i, k, j)%im
507 0 : C1(i, k, j)%im = coef(ix1 + idij)%re * work(i, k, j)%im - coef(ix1 + idij)%im * work(i, k, j)%re
508 : end do
509 : end do
510 : end do
511 : end if
512 : end subroutine
513 :
514 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
515 : #elif setFFTF_ENABLED && RK_ENABLED
516 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
517 :
518 : integer(IK) :: lenData, idl1, ido, ip, ix1, ix2, ix3, ix4, k1, l1, l2, na, lenFactor
519 26138 : lenFactor = size(factor, kind = IK)
520 26138 : lenData = size(data, kind = IK)
521 78414 : CHECK_ASSERTION(__LINE__, lenData == size(coef, kind = IK), SK_"@setFFTF(): The condition `size(data) == size(coef)` must hold. size(data), size(coef) = "//getStr([lenData, size(coef, kind = IK)]))
522 78414 : CHECK_ASSERTION(__LINE__, lenData == size(work, kind = IK), SK_"@setFFTF(): The condition `size(data) == size(work)` must hold. size(data), size(work) = "//getStr([lenData, size(work, kind = IK)]))
523 26138 : inwork = logical(1_IK < lenData, LK)
524 26138 : if (.not. inwork) return
525 : na = 1_IK
526 : l2 = lenData
527 26138 : ix1 = lenData
528 80982 : do k1 = 1_IK, lenFactor
529 54844 : ip = factor(lenFactor - k1 + 1_IK)
530 54844 : l1 = l2 / ip
531 54844 : ido = lenData / l2
532 54844 : idl1 = ido * l1
533 54844 : ix1 = ix1 - (ip - 1_IK) * ido
534 54844 : na = 1_IK - na
535 54844 : if (ip == 4) then
536 7592 : ix2 = ix1 + ido
537 7592 : ix3 = ix2 + ido
538 7592 : if (na /= 0_IK) then
539 4551 : call radf4(ido, l1, ix1, ix2, ix3, coef, work, data)
540 : else
541 3041 : call radf4(ido, l1, ix1, ix2, ix3, coef, data, work)
542 : end if
543 47252 : elseif (ip /= 2_IK) then
544 40329 : if (ip == 3_IK) then
545 13560 : ix2 = ix1 + ido
546 13560 : if (na /= 0_IK) then
547 8122 : call radf3(ido, l1, ix1, ix2, coef, work, data)
548 : else
549 5438 : call radf3(ido, l1, ix1, ix2, coef, data, work)
550 : end if
551 26769 : elseif (ip /= 5_IK) then
552 20397 : if (ido == 1_IK) na = 1_IK - na
553 20397 : if (na /= 0_IK) then
554 : na = 0_IK
555 20397 : call radfg(ido, ip, l1, idl1, ix1, coef, work, work, work, data, data)
556 : else
557 0 : call radfg(ido, ip, l1, idl1, ix1, coef, data, data, data, work, work)
558 : na = 1_IK
559 : end if
560 : else
561 6372 : ix2 = ix1 + ido
562 6372 : ix3 = ix2 + ido
563 6372 : ix4 = ix3 + ido
564 6372 : if (na /= 0_IK) then
565 2969 : call radf5(ido, l1, ix1, ix2, ix3, ix4, coef, work, data)
566 : else
567 3403 : call radf5(ido, l1, ix1, ix2, ix3, ix4, coef, data, work)
568 : end if
569 : end if
570 6923 : elseif (na /= 0_IK) then
571 4055 : call radf2(ido, l1, ix1, coef, work, data)
572 : else
573 2868 : call radf2(ido, l1, ix1, coef, data, work)
574 : end if
575 26138 : l2 = l1
576 : end do
577 26138 : inwork = logical(na /= 1_IK, LK)
578 : !if (na == 1_IK) return
579 : !data(1:lenData) = work(1:lenData)
580 :
581 : contains
582 :
583 6923 : pure subroutine radf2(ido, l1, ix1, coef, data, work)
584 : integer(IK) , intent(in) :: ido, l1, ix1
585 : real(TKC) , intent(in), contiguous :: coef(2:)
586 : real(TKC) , intent(in) :: data(ido, l1, 2)
587 : real(TKC) , intent(out) :: work(ido, 2, l1)
588 : integer(IK) :: i, ic, idp2, k
589 : complex(TKC) :: t2
590 13846 : do k = 1_IK, l1
591 6923 : work(1, 1, k) = data(1, k, 1) + data(1, k, 2)
592 13846 : work(ido, 2, k) = data(1, k, 1) - data(1, k, 2)
593 : end do
594 6923 : if (ido < 2_IK) return
595 6923 : if (ido /= 2_IK) then
596 6923 : idp2 = ido + 2_IK
597 13846 : do k = 1_IK, l1
598 13846 : do i = 3_IK, ido, 2_IK
599 99811 : ic = idp2 - i
600 99811 : t2%re = coef(ix1 + i - 2) * data(i - 1, k, 2) + coef(ix1 + i - 1) * data(i, k, 2)
601 99811 : t2%im = coef(ix1 + i - 2) * data(i, k, 2) - coef(ix1 + i - 1) * data(i - 1, k, 2)
602 99811 : work(i, 1, k) = data(i, k, 1) + t2%im
603 99811 : work(ic, 2, k) = t2%im - data(i, k, 1)
604 99811 : work(i - 1, 1, k) = data(i - 1, k, 1) + t2%re
605 99811 : work(ic - 1, 2, k) = data(i - 1, k, 1) - t2%re
606 : end do
607 : end do
608 6923 : if (mod(ido, 2_IK) == 1_IK) return
609 : end if
610 3842 : do k = 1_IK, l1
611 1921 : work(1, 2, k) = -data(ido, k, 2)
612 3842 : work(ido, 1, k) = data(ido, k, 1)
613 : end do
614 : end subroutine
615 :
616 13560 : pure subroutine radf3(ido, l1, ix1, ix2, coef, data, work)
617 : integer(IK) , intent(in) :: ido, l1, ix1, ix2
618 : real(TKC) , intent(in), contiguous :: coef(2:)
619 : real(TKC) , intent(in) :: data(ido, l1, 3)
620 : real(TKC) , intent(out) :: work(ido, 3, l1)
621 : complex(TKC), parameter :: TAU = cmplx(-.5_TKC, sqrt(3._TKC) / 2._TKC, TKC)
622 : complex(TKC) :: c2, d2, d3, t2, t3
623 : integer(IK) :: i, ic, idp2, k
624 82281 : do k = 1_IK, l1
625 68721 : c2%re = data(1, k, 2) + data(1, k, 3)
626 68721 : work(1, 1, k) = data(1, k, 1) + c2%re
627 68721 : work(1, 3, k) = TAU%im * (data(1, k, 3) - data(1, k, 2))
628 82281 : work(ido, 2, k) = data(1, k, 1) + TAU%re * c2%re
629 : end do
630 13560 : if (ido == 1_IK) return
631 11268 : idp2 = ido + 2_IK
632 41592 : do k = 1_IK, l1
633 41592 : do i = 3_IK, ido, 2_IK
634 100202 : ic = idp2 - i
635 100202 : d2%re = coef(ix1 + i - 2) * data(i - 1, k, 2) + coef(ix1 + i - 1) * data(i, k, 2)
636 100202 : d2%im = coef(ix1 + i - 2) * data(i, k, 2) - coef(ix1 + i - 1) * data(i - 1, k, 2)
637 100202 : d3%re = coef(ix2 + i - 2) * data(i - 1, k, 3) + coef(ix2 + i - 1) * data(i, k, 3)
638 100202 : d3%im = coef(ix2 + i - 2) * data(i, k, 3) - coef(ix2 + i - 1) * data(i - 1, k, 3)
639 100202 : c2%re = d2%re + d3%re
640 100202 : c2%im = d2%im + d3%im
641 100202 : work(i - 1, 1, k) = data(i - 1, k, 1) + c2%re
642 100202 : work(i, 1, k) = data(i, k, 1) + c2%im
643 100202 : t2%re = data(i - 1, k, 1) + TAU%re * c2%re
644 100202 : t2%im = data(i, k, 1) + TAU%re * c2%im
645 100202 : t3%re = TAU%im * (d2%im - d3%im)
646 100202 : t3%im = TAU%im * (d3%re - d2%re)
647 100202 : work(i - 1, 3, k) = t2%re + t3%re
648 100202 : work(ic - 1, 2, k) = t2%re - t3%re
649 100202 : work(i, 3, k) = t2%im + t3%im
650 100202 : work(ic, 2, k) = t3%im - t2%im
651 : end do
652 : end do
653 : end subroutine
654 :
655 7592 : pure subroutine radf4(ido, l1, ix1, ix2, ix3, coef, data, work)
656 : integer(IK) , intent(in) :: ido, l1, ix1, ix2, ix3
657 : real(TKC) , intent(in), contiguous :: coef(2:)
658 : real(TKC) , intent(in) :: data(ido, l1, 4)
659 : real(TKC) , intent(out) :: work(ido, 4, l1)
660 : real(TKC) , parameter :: NHSQT2 = -sqrt(2._TKC) / 2._TKC
661 : complex(TKC) :: c2, c3, c4, t1, t2, t3, t4
662 : integer(IK) :: i, ic, idp2, k
663 27699 : do k = 1_IK, l1
664 20107 : t1%re = data(1, k, 2) + data(1, k, 4)
665 20107 : t2%re = data(1, k, 1) + data(1, k, 3)
666 20107 : work(1, 1, k) = t1%re + t2%re
667 20107 : work(ido, 4, k) = t2%re - t1%re
668 20107 : work(ido, 2, k) = data(1, k, 1) - data(1, k, 3)
669 27699 : work(1, 3, k) = data(1, k, 4) - data(1, k, 2)
670 : end do
671 7592 : if (ido < 2_IK) return
672 6823 : if (ido /= 2_IK) then
673 6823 : idp2 = ido + 2_IK
674 18270 : do k = 1_IK, l1
675 18270 : do i = 3_IK, ido, 2_IK
676 43532 : ic = idp2 - i
677 43532 : c2%re = coef(ix1 + i - 2) * data(i - 1, k, 2) + coef(ix1 + i - 1) * data(i, k, 2)
678 43532 : c2%im = coef(ix1 + i - 2) * data(i, k, 2) - coef(ix1 + i - 1) * data(i - 1, k, 2)
679 43532 : c3%re = coef(ix2 + i - 2) * data(i - 1, k, 3) + coef(ix2 + i - 1) * data(i, k, 3)
680 43532 : c3%im = coef(ix2 + i - 2) * data(i, k, 3) - coef(ix2 + i - 1) * data(i - 1, k, 3)
681 43532 : c4%re = coef(ix3 + i - 2) * data(i - 1, k, 4) + coef(ix3 + i - 1) * data(i, k, 4)
682 43532 : c4%im = coef(ix3 + i - 2) * data(i, k, 4) - coef(ix3 + i - 1) * data(i - 1, k, 4)
683 43532 : t1%re = c2%re + c4%re
684 43532 : t4%re = c4%re - c2%re
685 43532 : t1%im = c2%im + c4%im
686 43532 : t4%im = c2%im - c4%im
687 43532 : t2%im = data(i, k, 1) + c3%im
688 43532 : t3%im = data(i, k, 1) - c3%im
689 43532 : t2%re = data(i - 1, k, 1) + c3%re
690 43532 : t3%re = data(i - 1, k, 1) - c3%re
691 43532 : work(i - 1, 1, k) = t1%re + t2%re
692 43532 : work(ic - 1, 4, k) = t2%re - t1%re
693 43532 : work(i, 1, k) = t1%im + t2%im
694 43532 : work(ic, 4, k) = t1%im - t2%im
695 43532 : work(i - 1, 3, k) = t4%im + t3%re
696 43532 : work(ic - 1, 2, k) = t3%re - t4%im
697 43532 : work(i, 3, k) = t4%re + t3%im
698 43532 : work(ic, 2, k) = t4%re - t3%im
699 : end do
700 : end do
701 6823 : if (mod(ido, 2_IK) == 1_IK) return
702 : end if
703 4646 : do k = 1_IK, l1
704 3048 : t1%im = NHSQT2 * (data(ido, k, 2) + data(ido, k, 4))
705 3048 : t1%re = NHSQT2 * (data(ido, k, 4) - data(ido, k, 2))
706 3048 : work(ido, 1, k) = t1%re + data(ido, k, 1)
707 3048 : work(ido, 3, k) = data(ido, k, 1) - t1%re
708 3048 : work(1, 2, k) = t1%im - data(ido, k, 3)
709 4646 : work(1, 4, k) = t1%im + data(ido, k, 3)
710 : end do
711 : end subroutine
712 :
713 6372 : pure subroutine radf5(ido, l1, ix1, ix2, ix3, ix4, coef, data, work)
714 : integer(IK) , intent(in) :: ido, l1, ix1, ix2, ix3, ix4
715 : real(TKC) , intent(in), contiguous :: coef(2:)
716 : real(TKC) , intent(in) :: data(ido, l1, 5)
717 : real(TKC) , intent(out) :: work(ido, 5, l1)
718 : real(TKC) , parameter :: PI = acos(-1._TKC)
719 : complex(TKC), parameter :: T11 = cmplx(cos(2._TKC * PI / 5._TKC), sin(2._TKC * PI / 5._TKC), TKC)
720 : complex(TKC), parameter :: T12 = cmplx(cos(4._TKC * PI / 5._TKC), sin(4._TKC * PI / 5._TKC), TKC)
721 : complex(TKC) :: c2, c3, c4, c5, d2, d3, d4, d5, t2, t3, t4, t5
722 : integer(IK) :: i, ic, idp2, k
723 47476 : do k = 1_IK, l1
724 41104 : c2%re = data(1, k, 5) + data(1, k, 2)
725 41104 : c5%im = data(1, k, 5) - data(1, k, 2)
726 41104 : c3%re = data(1, k, 4) + data(1, k, 3)
727 41104 : c4%im = data(1, k, 4) - data(1, k, 3)
728 41104 : work(1, 1, k) = data(1, k, 1) + c2%re + c3%re
729 41104 : work(ido, 2, k) = data(1, k, 1) + T11%re * c2%re + T12%re * c3%re
730 41104 : work(1, 3, k) = T11%im * c5%im + T12%im * c4%im
731 41104 : work(ido, 4, k) = data(1, k, 1) + T12%re * c2%re + T11%re * c3%re
732 47476 : work(1, 5, k) = T12%im * c5%im - T11%im * c4%im
733 : end do
734 6372 : if (ido == 1_IK) return
735 2969 : idp2 = ido + 2_IK
736 7496 : do k = 1_IK, l1
737 7496 : do i = 3_IK, ido, 2_IK
738 16118 : ic = idp2 - i
739 16118 : d2%re = coef(ix1 + i - 2) * data(i - 1, k, 2) + coef(ix1 + i - 1) * data(i, k, 2)
740 16118 : d2%im = coef(ix1 + i - 2) * data(i, k, 2) - coef(ix1 + i - 1) * data(i - 1, k, 2)
741 16118 : d3%re = coef(ix2 + i - 2) * data(i - 1, k, 3) + coef(ix2 + i - 1) * data(i, k, 3)
742 16118 : d3%im = coef(ix2 + i - 2) * data(i, k, 3) - coef(ix2 + i - 1) * data(i - 1, k, 3)
743 16118 : d4%re = coef(ix3 + i - 2) * data(i - 1, k, 4) + coef(ix3 + i - 1) * data(i, k, 4)
744 16118 : d4%im = coef(ix3 + i - 2) * data(i, k, 4) - coef(ix3 + i - 1) * data(i - 1, k, 4)
745 16118 : d5%re = coef(ix4 + i - 2) * data(i - 1, k, 5) + coef(ix4 + i - 1) * data(i, k, 5)
746 16118 : d5%im = coef(ix4 + i - 2) * data(i, k, 5) - coef(ix4 + i - 1) * data(i - 1, k, 5)
747 16118 : c2%re = d2%re + d5%re
748 16118 : c5%im = d5%re - d2%re
749 16118 : c5%re = d2%im - d5%im
750 16118 : c2%im = d2%im + d5%im
751 16118 : c3%re = d3%re + d4%re
752 16118 : c4%im = d4%re - d3%re
753 16118 : c4%re = d3%im - d4%im
754 16118 : c3%im = d3%im + d4%im
755 16118 : work(i - 1, 1, k) = data(i - 1, k, 1) + c2%re + c3%re
756 16118 : work(i, 1, k) = data(i, k, 1) + c2%im + c3%im
757 16118 : t2%re = data(i - 1, k, 1) + T11%re * c2%re + T12%re * c3%re
758 16118 : t2%im = data(i, k, 1) + T11%re * c2%im + T12%re * c3%im
759 16118 : t3%re = data(i - 1, k, 1) + T12%re * c2%re + T11%re * c3%re
760 16118 : t3%im = data(i, k, 1) + T12%re * c2%im + T11%re * c3%im
761 16118 : t5%re = T11%im * c5%re + T12%im * c4%re
762 16118 : t5%im = T11%im * c5%im + T12%im * c4%im
763 16118 : t4%re = T12%im * c5%re - T11%im * c4%re
764 16118 : t4%im = T12%im * c5%im - T11%im * c4%im
765 16118 : work(i - 1, 3, k) = t2%re + t5%re
766 16118 : work(ic - 1, 2, k) = t2%re - t5%re
767 16118 : work(i, 3, k) = t2%im + t5%im
768 16118 : work(ic, 2, k) = t5%im - t2%im
769 16118 : work(i - 1, 5, k) = t3%re + t4%re
770 16118 : work(ic - 1, 4, k) = t3%re - t4%re
771 16118 : work(i, 5, k) = t3%im + t4%im
772 16118 : work(ic, 4, k) = t4%im - t3%im
773 : end do
774 : end do
775 : end subroutine radf5
776 :
777 20397 : pure subroutine radfg(ido, ip, l1, idl1, ix1, coef, data, C1, C2, work, Ch2)
778 : integer(IK) , intent(in) :: ido, ip, l1, idl1, ix1
779 : real(TKC) , intent(in), contiguous :: coef(2:)
780 : real(TKC) , intent(inout) :: C1(ido,l1,ip), C2(idl1,ip)
781 : real(TKC) , intent(inout) :: data(ido,ip,l1), work(ido,l1,ip), Ch2(idl1,ip)
782 : real(TKC) , parameter :: TWO_PI = 2._TKC * acos(-1._TKC)
783 : real(TKC) :: ar1h, ar2h, arg, dc2, dcp, ds2, dsp
784 : integer(IK) :: i, ic, idij, idp2, jk, ipp2, ipph, is, j, j2, jc, k, l, lc, nbd
785 : complex(TKC) :: a1, a2
786 20397 : arg = TWO_PI / real(ip, TKC)
787 20397 : dcp = cos(arg)
788 20397 : dsp = sin(arg)
789 20397 : ipph = (ip + 1_IK) / 2_IK
790 20397 : ipp2 = ip + 2_IK
791 20397 : idp2 = ido + 2_IK
792 20397 : nbd = (ido - 1_IK) / 2_IK
793 20397 : if (ido == 1_IK) then
794 86542 : C2(1:idl1, 1) = Ch2(1:idl1, 1)
795 : else
796 7589 : Ch2(1:idl1, 1) = C2(1:idl1, 1)
797 9405 : work(1, 1:l1, 2:ip) = C1(1, 1:l1, 2:ip)
798 723 : if (nbd > l1) then
799 723 : is = -ido
800 5061 : do j = 2_IK, ip
801 4338 : is = is + ido
802 9405 : do k = 1_IK, l1
803 : idij = is
804 8682 : do i = 3_IK, ido, 2_IK
805 18426 : idij = idij + 2_IK
806 18426 : work(i - 1, k, j) = coef(ix1 + idij - 1) * C1(i - 1, k, j) + coef(ix1 + idij) * C1(i, k, j)
807 18426 : work(i, k, j) = coef(ix1 + idij - 1) * C1(i, k, j) - coef(ix1 + idij) * C1(i - 1, k, j)
808 : end do
809 : end do
810 : end do
811 : else
812 0 : is = -ido
813 0 : do j = 2_IK, ip
814 0 : is = is + ido
815 : idij = is
816 0 : do i = 3_IK, ido, 2_IK
817 0 : idij = idij + 2_IK
818 0 : do k = 1_IK, l1
819 0 : work(i - 1, k, j) = coef(ix1 + idij - 1) * C1(i - 1, k, j) + coef(ix1 + idij) * C1(i, k, j)
820 0 : work(i, k, j) = coef(ix1 + idij - 1) * C1(i, k, j) - coef(ix1 + idij) * C1(i - 1, k, j)
821 : end do
822 : end do
823 : end do
824 : end if
825 723 : if (nbd < l1) then
826 0 : do j = 2_IK, ipph
827 0 : jc = ipp2 - j
828 0 : do i = 3_IK, ido, 2_IK
829 0 : do k = 1_IK, l1
830 0 : C1(i - 1, k, j) = work(i - 1, k, j) + work(i - 1, k, jc)
831 0 : C1(i - 1, k, jc) = work(i, k, j) - work(i, k, jc)
832 0 : C1(i, k, j) = work(i, k, j) + work(i, k, jc)
833 0 : C1(i, k, jc) = work(i - 1, k, jc) - work(i - 1, k, j)
834 : end do
835 : end do
836 : end do
837 : else
838 2892 : do j = 2_IK, ipph
839 2169 : jc = ipp2 - j
840 5064 : do k = 1_IK, l1
841 4341 : do i = 3_IK, ido, 2_IK
842 9213 : C1(i - 1, k, j) = work(i - 1, k, j) + work(i - 1, k, jc)
843 9213 : C1(i - 1, k, jc) = work(i, k, j) - work(i, k, jc)
844 9213 : C1(i, k, j) = work(i, k, j) + work(i, k, jc)
845 9213 : C1(i, k, jc) = work(i - 1, k, jc) - work(i - 1, k, j)
846 : end do
847 : end do
848 : end do
849 : end if
850 : end if
851 339507 : do j = 2_IK, ipph
852 319110 : jc = ipp2 - j
853 913045 : do k = 1_IK, l1
854 573538 : C1(1, k, j) = work(1, k, j) + work(1, k, jc)
855 892648 : C1(1, k, jc) = work(1, k, jc) - work(1, k, j)
856 : end do
857 : end do
858 15131 : a1%re = 1._TKC
859 15131 : a1%im = 0._TKC
860 339507 : do l = 2_IK, ipph
861 319110 : lc = ipp2 - l
862 319110 : ar1h = dcp * a1%re - dsp * a1%im
863 319110 : a1%im = dcp * a1%im + dsp * a1%re
864 319110 : a1%re = ar1h
865 911074 : do jk = 1, idl1
866 591964 : Ch2(jk, l) = C2(jk, 1) + a1%re * C2(jk, 2)
867 911074 : Ch2(jk, lc) = a1%im * C2(jk, ip)
868 : end do
869 137901 : dc2 = a1%re
870 137901 : ds2 = a1%im
871 212479 : a2%re = a1%re
872 212479 : a2%im = a1%im
873 348726807 : do j = 3_IK, ipph
874 348387300 : jc = ipp2 - j
875 348387300 : ar2h = dc2 * a2%re - ds2 * a2%im
876 348387300 : a2%im = dc2 * a2%im + ds2 * a2%re
877 348387300 : a2%re = ar2h
878 698770582 : do jk = 1_IK, idl1
879 350064172 : Ch2(jk, l) = Ch2(jk, l) + a2%re * C2(jk, j)
880 698451472 : Ch2(jk, lc) = Ch2(jk, lc) + a2%im * C2(jk, jc)
881 : end do
882 : end do
883 : end do
884 339507 : do j = 2_IK, ipph
885 931471 : do jk = 1_IK, idl1
886 911074 : Ch2(jk, 1) = Ch2(jk, 1) + C2(jk, j)
887 : end do
888 : end do
889 20397 : if (ido < l1) then
890 24674 : do i = 1, ido
891 84205 : do k = 1, l1
892 71868 : data(i, 1, k) = work(i, k, 1)
893 : end do
894 : end do
895 : else
896 16121 : do k = 1_IK, l1
897 30324 : do i = 1, ido
898 22264 : data(i, 1, k) = work(i, k, 1)
899 : end do
900 : end do
901 : end if
902 339507 : do j = 2_IK, ipph
903 319110 : jc = ipp2 - j
904 319110 : j2 = j + j
905 913045 : do k = 1, l1
906 573538 : data(ido, j2 - 2, k) = work(1, k, j)
907 892648 : data(1, j2 - 1, k) = work(1, k, jc)
908 : end do
909 : end do
910 20397 : if (ido == 1_IK) return
911 723 : if (nbd < l1) then
912 0 : do j = 2_IK, ipph
913 0 : jc = ipp2 - j
914 0 : j2 = j + j
915 0 : do i = 3_IK, ido, 2_IK
916 0 : ic = idp2 - i
917 0 : do k = 1, l1
918 0 : data(i - 1, j2 - 1, k) = work(i - 1, k, j) + work(i - 1, k, jc)
919 0 : data(ic - 1, j2 - 2, k) = work(i - 1, k, j) - work(i - 1, k, jc)
920 0 : data(i, j2 - 1, k) = work(i, k, j) + work(i, k, jc)
921 0 : data(ic, j2 - 2, k) = work(i, k, jc) - work(i, k, j)
922 : end do
923 : end do
924 : end do
925 : else
926 2892 : do j = 2_IK, ipph
927 2169 : jc = ipp2 - j
928 2169 : j2 = j + j
929 5064 : do k = 1_IK, l1
930 4341 : do i = 3_IK, ido, 2_IK
931 9213 : ic = idp2 - i
932 9213 : data(i - 1, j2 - 1, k) = work(i - 1, k, j) + work(i - 1, k, jc)
933 9213 : data(ic - 1, j2 - 2, k) = work(i - 1, k, j) - work(i - 1, k, jc)
934 9213 : data(i, j2 - 1, k) = work(i, k, j) + work(i, k, jc)
935 9213 : data(ic, j2 - 2, k) = work(i, k, jc) - work(i, k, j)
936 : end do
937 : end do
938 : end do
939 : end if
940 : end subroutine
941 :
942 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
943 : #elif setFFTR_ENABLED && CK_ENABLED
944 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
945 :
946 : integer(IK) :: lenData, idl1, ido, ip, ix1, ix2, ix3, ix4, k1, l1, l2, na, nac, lenFactor
947 14923 : lenFactor = size(factor, kind = IK)
948 14923 : lenData = size(data, kind = IK)
949 44769 : CHECK_ASSERTION(__LINE__, lenData == size(coef, kind = IK), SK_"@setFFTR(): The condition `size(data) == size(coef)` must hold. size(data), size(coef) = "//getStr([lenData, size(coef, kind = IK)]))
950 44769 : CHECK_ASSERTION(__LINE__, lenData == size(work, kind = IK), SK_"@setFFTR(): The condition `size(data) == size(work)` must hold. size(data), size(work) = "//getStr([lenData, size(work, kind = IK)]))
951 14923 : inwork = logical(1_IK < lenData, LK)
952 14923 : if (.not. inwork) return
953 : na = 0_IK
954 14923 : l1 = 1_IK
955 14923 : ix1 = 1_IK
956 45577 : do k1 = 1_IK, lenFactor
957 30654 : ip = factor(k1)
958 30654 : l2 = ip * l1
959 30654 : ido = lenData / l2
960 : !idot = ido + ido
961 30654 : idl1 = ido * l1
962 30654 : if (ip == 4_IK) then
963 4190 : ix2 = ix1 + ido
964 4190 : ix3 = ix2 + ido
965 4190 : if (na /= 0_IK) then
966 2005 : call passb4(ido, l1, ix1, ix2, ix3, coef, work, data)
967 : else
968 2185 : call passb4(ido, l1, ix1, ix2, ix3, coef, data, work)
969 : end if
970 4190 : na = 1_IK - na
971 26464 : elseif (ip == 2_IK) then
972 4108 : if (na /= 0_IK) then
973 0 : call passb2(ido, l1, ix1, coef, work, data)
974 : else
975 4108 : call passb2(ido, l1, ix1, coef, data, work)
976 : end if
977 4108 : na = 1_IK - na
978 22356 : elseif (ip == 3_IK) then
979 6964 : ix2 = ix1 + ido
980 6964 : if (na /= 0_IK) then
981 2696 : call passb3(ido, l1, ix1, ix2, coef, work, data)
982 : else
983 4268 : call passb3(ido, l1, ix1, ix2, coef, data, work)
984 : end if
985 6964 : na = 1_IK - na
986 15392 : elseif (ip /= 5_IK) then
987 11724 : if (na /= 0_IK) then
988 4864 : call passb(nac, ido, ip, l1, idl1, ix1, coef, work, work, work, data, data)
989 : else
990 6860 : call passb(nac, ido, ip, l1, idl1, ix1, coef, data, data, data, work, work)
991 : end if
992 11724 : if (nac /= 0_IK) na = 1_IK - na
993 : else
994 3668 : ix2 = ix1 + ido
995 3668 : ix3 = ix2 + ido
996 3668 : ix4 = ix3 + ido
997 3668 : if (na /= 0_IK) then
998 1196 : call passb5(ido, l1, ix1, ix2, ix3, ix4, coef, work, data)
999 : else
1000 2472 : call passb5(ido, l1, ix1, ix2, ix3, ix4, coef, data, work)
1001 : end if
1002 3668 : na = 1_IK - na
1003 : end if
1004 30654 : l1 = l2
1005 45577 : ix1 = ix1 + (ip - 1_IK) * ido
1006 : end do
1007 14923 : inwork = logical(na /= 0_IK, LK)
1008 : !if (na == 0_IK) return
1009 : !data(1 : lenData) = work(1 : lenData)
1010 :
1011 : contains
1012 :
1013 4108 : pure subroutine passb2(ido, l1, ix1, coef, data, work)
1014 : integer(IK) , intent(in) :: ido, l1, ix1
1015 : complex(TKC), intent(in) , contiguous :: coef(2:)
1016 : complex(TKC), intent(in) :: data(ido, 2, l1)
1017 : complex(TKC), intent(out) :: work(ido, l1, 2)
1018 : complex(TKC) :: t2
1019 : integer(IK) :: i, k
1020 4108 : if (ido > 1_IK) then
1021 8216 : do k = 1_IK, l1
1022 128919 : do i = 1_IK, ido
1023 120703 : work(i, k, 1) = data(i, 1, k) + data(i, 2, k)
1024 120703 : t2 = data(i, 1, k) - data(i, 2, k)
1025 120703 : work(i, k, 2)%im = coef(ix1 + i)%re * t2%im + coef(ix1 + i)%im * t2%re
1026 124811 : work(i, k, 2)%re = coef(ix1 + i)%re * t2%re - coef(ix1 + i)%im * t2%im
1027 : end do
1028 : end do
1029 : else
1030 0 : do k = 1_IK, l1
1031 0 : work(1, k, 1) = data(1, 1, k) + data(1, 2, k)
1032 0 : work(1, k, 2) = data(1, 1, k) - data(1, 2, k)
1033 : end do
1034 : end if
1035 4108 : end subroutine
1036 :
1037 6964 : pure subroutine passb3(ido, l1, ix1, ix2, coef, data, work)
1038 : integer(IK) , intent(in) :: ido, l1, ix1, ix2
1039 : complex(TKC), intent(in) , contiguous :: coef(2:)
1040 : complex(TKC), intent(in) :: data(ido, 3, l1)
1041 : complex(TKC), intent(out) :: work(ido, l1, 3)
1042 : complex(TKC), parameter :: TAU = cmplx(-.5_TKC, sqrt(3._TKC) / 2._TKC, TKC)
1043 : complex(TKC) :: c2, c3, d2, d3, t2
1044 : integer(IK) :: i, k
1045 6964 : if (ido /= 1_IK) then
1046 19963 : do k = 1_IK, l1
1047 141770 : do i = 1_IK, ido
1048 121807 : t2 = data(i, 2, k) + data(i, 3, k)
1049 121807 : c2 = data(i, 1, k) + TAU%re * t2
1050 121807 : work(i, k, 1) = data(i, 1, k) + t2
1051 121807 : c3 = TAU%im * (data(i, 2, k) - data(i, 3, k))
1052 121807 : d2%re = c2%re - c3%im
1053 121807 : d2%im = c2%im + c3%re
1054 121807 : d3%re = c2%re + c3%im
1055 121807 : d3%im = c2%im - c3%re
1056 121807 : work(i, k, 2)%re = coef(ix1 + i)%re * d2%re - coef(ix1 + i)%im * d2%im
1057 121807 : work(i, k, 2)%im = coef(ix1 + i)%re * d2%im + coef(ix1 + i)%im * d2%re
1058 121807 : work(i, k, 3)%re = coef(ix2 + i)%re * d3%re - coef(ix2 + i)%im * d3%im
1059 135621 : work(i, k, 3)%im = coef(ix2 + i)%re * d3%im + coef(ix2 + i)%im * d3%re
1060 : end do
1061 : end do
1062 : else
1063 11298 : do k = 1_IK, l1
1064 10483 : t2 = data(1, 2, k) + data(1, 3, k)
1065 10483 : c2 = data(1, 1, k) + TAU%re * t2
1066 10483 : work(1, k, 1) = data(1, 1, k) + t2
1067 10483 : c3 = TAU%im * (data(1, 2, k) - data(1, 3, k))
1068 10483 : work(1, k, 2)%re = c2%re - c3%im
1069 10483 : work(1, k, 2)%im = c2%im + c3%re
1070 10483 : work(1, k, 3)%re = c2%re + c3%im
1071 11298 : work(1, k, 3)%im = c2%im - c3%re
1072 : end do
1073 : end if
1074 6964 : end subroutine
1075 :
1076 4190 : pure subroutine passb4(ido, l1, ix1, ix2, ix3, coef, data, work)
1077 : integer(IK) , intent(in) :: ido, l1, ix1, ix2, ix3
1078 : complex(TKC), intent(in) , contiguous :: coef(2:)
1079 : complex(TKC), intent(in) :: data(ido, 4, l1)
1080 : complex(TKC), intent(out) :: work(ido, l1, 4)
1081 : complex(TKC) :: c2, c3, c4, t1, t2, t3, t4
1082 : integer(IK) :: i, k
1083 4190 : if (ido /= 1_IK) then
1084 9659 : do k = 1_IK, l1
1085 59488 : do i = 1_IK, ido, 1_IK
1086 49829 : t1 = data(i, 1, k) - data(i, 3, k)
1087 49829 : t2 = data(i, 1, k) + data(i, 3, k)
1088 49829 : t3 = data(i, 2, k) + data(i, 4, k)
1089 49829 : t4%re = data(i, 4, k)%im - data(i, 2, k)%im
1090 49829 : t4%im = data(i, 2, k)%re - data(i, 4, k)%re
1091 49829 : work(i, k, 1) = t2 + t3
1092 49829 : c3 = t2 - t3
1093 49829 : c2 = t1 + t4
1094 49829 : c4 = t1 - t4
1095 49829 : work(i, k, 2)%re = coef(ix1 + i)%re * c2%re - coef(ix1 + i)%im * c2%im
1096 49829 : work(i, k, 2)%im = coef(ix1 + i)%re * c2%im + coef(ix1 + i)%im * c2%re
1097 49829 : work(i, k, 3)%re = coef(ix2 + i)%re * c3%re - coef(ix2 + i)%im * c3%im
1098 49829 : work(i, k, 3)%im = coef(ix2 + i)%re * c3%im + coef(ix2 + i)%im * c3%re
1099 49829 : work(i, k, 4)%re = coef(ix3 + i)%re * c4%re - coef(ix3 + i)%im * c4%im
1100 56241 : work(i, k, 4)%im = coef(ix3 + i)%re * c4%im + coef(ix3 + i)%im * c4%re
1101 : end do
1102 : end do
1103 : else
1104 12767 : do k = 1_IK, l1
1105 11824 : t1 = data(1, 1, k) - data(1, 3, k)
1106 11824 : t2 = data(1, 1, k) + data(1, 3, k)
1107 11824 : t3 = data(1, 2, k) + data(1, 4, k)
1108 11824 : t4%re = data(1, 4, k)%im - data(1, 2, k)%im
1109 11824 : t4%im = data(1, 2, k)%re - data(1, 4, k)%re
1110 11824 : work(1, k, 1) = t2 + t3
1111 11824 : work(1, k, 3) = t2 - t3
1112 11824 : work(1, k, 2) = t1 + t4
1113 12767 : work(1, k, 4) = t1 - t4
1114 : end do
1115 : end if
1116 4190 : end subroutine passb4
1117 :
1118 3668 : pure subroutine passb5(ido, l1, ix1, ix2, ix3, ix4, coef, data, work)
1119 : integer(IK) , intent(in) :: ido, l1, ix1, ix2, ix3, ix4
1120 : complex(TKC), intent(in) , contiguous :: coef(2:)
1121 : complex(TKC), intent(in) :: data(ido, 5, l1)
1122 : complex(TKC), intent(out) :: work(ido, l1, 5)
1123 : complex(TKC) :: c2, c3, c4, c5, d2, d3, d4, d5, t2, t3, t4, t5
1124 : real(TKC) , parameter :: pi = acos( - 1._TKC)
1125 : complex(TKC), parameter :: T11 = cmplx(cos(2._TKC * pi / 5._TKC), sin(2._TKC * pi / 5._TKC), TKC)
1126 : complex(TKC), parameter :: T12 = cmplx(cos(4._TKC * pi / 5._TKC), sin(4._TKC * pi / 5._TKC), TKC)
1127 : integer(IK) :: i, k
1128 3668 : if (ido /= 1_IK) then
1129 4229 : do k = 1_IK, l1
1130 25076 : do i = 1_IK, ido
1131 20847 : t2 = data(i, 2, k) + data(i, 5, k)
1132 20847 : t3 = data(i, 3, k) + data(i, 4, k)
1133 20847 : t4 = data(i, 3, k) - data(i, 4, k)
1134 20847 : t5 = data(i, 2, k) - data(i, 5, k)
1135 20847 : work(i, k, 1) = data(i, 1, k) + t2 + t3
1136 20847 : c2%re = data(i, 1, k)%re + T11%re * t2%re + T12%re * t3%re
1137 20847 : c2%im = data(i, 1, k)%im + T11%re * t2%im + T12%re * t3%im
1138 20847 : c3%re = data(i, 1, k)%re + T12%re * t2%re + T11%re * t3%re
1139 20847 : c3%im = data(i, 1, k)%im + T12%re * t2%im + T11%re * t3%im
1140 20847 : c5%re = T11%im * t5%re + T12%im * t4%re
1141 20847 : c5%im = T11%im * t5%im + T12%im * t4%im
1142 20847 : c4%re = T12%im * t5%re - T11%im * t4%re
1143 20847 : c4%im = T12%im * t5%im - T11%im * t4%im
1144 20847 : d3%re = c3%re - c4%im
1145 20847 : d3%im = c3%im + c4%re
1146 20847 : d4%re = c3%re + c4%im
1147 20847 : d4%im = c3%im - c4%re
1148 20847 : d5%re = c2%re + c5%im
1149 20847 : d5%im = c2%im - c5%re
1150 20847 : d2%re = c2%re - c5%im
1151 20847 : d2%im = c2%im + c5%re
1152 20847 : work(i, k, 2)%re = coef(ix1 + i)%re * d2%re - coef(ix1 + i)%im * d2%im
1153 20847 : work(i, k, 2)%im = coef(ix1 + i)%re * d2%im + coef(ix1 + i)%im * d2%re
1154 20847 : work(i, k, 3)%re = coef(ix2 + i)%re * d3%re - coef(ix2 + i)%im * d3%im
1155 20847 : work(i, k, 3)%im = coef(ix2 + i)%re * d3%im + coef(ix2 + i)%im * d3%re
1156 20847 : work(i, k, 4)%re = coef(ix3 + i)%re * d4%re - coef(ix3 + i)%im * d4%im
1157 20847 : work(i, k, 4)%im = coef(ix3 + i)%re * d4%im + coef(ix3 + i)%im * d4%re
1158 20847 : work(i, k, 5)%re = coef(ix4 + i)%re * d5%re - coef(ix4 + i)%im * d5%im
1159 23348 : work(i, k, 5)%im = coef(ix4 + i)%re * d5%im + coef(ix4 + i)%im * d5%re
1160 : end do
1161 : end do
1162 : else
1163 21432 : do k = 1_IK, l1
1164 19492 : t2 = data(1, 2, k) + data(1, 5, k)
1165 19492 : t3 = data(1, 3, k) + data(1, 4, k)
1166 19492 : t4 = data(1, 3, k) - data(1, 4, k)
1167 19492 : t5 = data(1, 2, k) - data(1, 5, k)
1168 19492 : work(1, k, 1) = data(1, 1, k) + t2 + t3
1169 19492 : c2%re = data(1, 1, k)%re + T11%re * t2%re + T12%re * t3%re
1170 19492 : c2%im = data(1, 1, k)%im + T11%re * t2%im + T12%re * t3%im
1171 19492 : c3%re = data(1, 1, k)%re + T12%re * t2%re + T11%re * t3%re
1172 19492 : c3%im = data(1, 1, k)%im + T12%re * t2%im + T11%re * t3%im
1173 19492 : c5%re = T11%im * t5%re + T12%im * t4%re
1174 19492 : c5%im = T11%im * t5%im + T12%im * t4%im
1175 19492 : c4%re = T12%im * t5%re - T11%im * t4%re
1176 19492 : c4%im = T12%im * t5%im - T11%im * t4%im
1177 19492 : work(1, k, 2)%re = c2%re - c5%im
1178 19492 : work(1, k, 2)%im = c2%im + c5%re
1179 19492 : work(1, k, 3)%re = c3%re - c4%im
1180 19492 : work(1, k, 3)%im = c3%im + c4%re
1181 19492 : work(1, k, 4)%re = c3%re + c4%im
1182 19492 : work(1, k, 4)%im = c3%im - c4%re
1183 19492 : work(1, k, 5)%re = c2%re + c5%im
1184 21432 : work(1, k, 5)%im = c2%im - c5%re
1185 : end do
1186 : end if
1187 3668 : end subroutine passb5
1188 :
1189 11724 : pure subroutine passb(nac, ido, ip, l1, idl1, ix1, coef, data, C1, C2, work, Ch2)
1190 : integer(IK) , intent(out) :: nac
1191 : integer(IK) , intent(in) :: ido, ip, l1, idl1, ix1
1192 : complex(TKC), intent(in) , contiguous :: coef(2:)
1193 : complex(TKC), intent(in) :: data(ido, ip, l1)
1194 : complex(TKC), intent(inout) :: C1(ido, l1, ip)
1195 : complex(TKC), intent(inout) :: work(ido, l1, ip), C2(idl1,ip), Ch2(idl1,ip)
1196 : integer(IK) :: i, idij, idj, idl, idlj, idp, jk, inc, ipp2, ipph, j, jc, k, l, lc!, nt
1197 : !idot = ido / 2_IK
1198 : !nt = ip * idl1 what in the world is this doing here?
1199 11724 : idp = ip * ido
1200 11724 : ipp2 = ip + 2_IK
1201 11724 : ipph = (ip + 1_IK) / 2_IK
1202 11724 : if (2_IK * ido < l1) then
1203 37668 : do j = 2_IK, ipph
1204 32263 : jc = ipp2 - j
1205 69931 : do i = 1_IK, ido
1206 218184 : do k = 1_IK, l1
1207 153658 : work(i, k, j) = data(i, j, k) + data(i, jc, k)
1208 185921 : work(i, k, jc) = data(i, j, k) - data(i, jc, k)
1209 : end do
1210 : end do
1211 : end do
1212 10810 : do i = 1_IK, ido
1213 40494 : do k = 1_IK, l1
1214 35089 : work(i, k, 1) = data(i, 1, k)
1215 : end do
1216 : end do
1217 : else
1218 136605 : do j = 2_IK, ipph
1219 130286 : jc = ipp2 - j
1220 290226 : do k = 1_IK, l1
1221 450578 : do i = 1_IK, ido
1222 166671 : work(i, k, j) = data(i, j, k) + data(i, jc, k)
1223 320292 : work(i, k, jc) = data(i, j, k) - data(i, jc, k)
1224 : end do
1225 : end do
1226 : end do
1227 14228 : do k = 1_IK, l1
1228 26487 : do i = 1_IK, ido
1229 20168 : work(i, k, 1) = data(i, 1, k)
1230 : end do
1231 : end do
1232 : end if
1233 : inc = 0_IK
1234 11724 : idl = 1_IK - ido
1235 174273 : do l = 2_IK, ipph
1236 162549 : lc = ipp2 - l
1237 162549 : idl = idl + ido
1238 482878 : do jk = 1_IK, idl1
1239 320329 : C2(jk, l) = Ch2(jk, 1) + coef(ix1 + idl)%re * Ch2(jk, 2)
1240 482878 : C2(jk, lc) = coef(ix1 + idl)%im * Ch2(jk, ip)
1241 : end do
1242 : idlj = idl
1243 162549 : inc = inc + ido
1244 3817431 : do j = 3_IK, ipph
1245 3643158 : jc = ipp2 - j
1246 3643158 : idlj = idlj + inc
1247 3643158 : if (idlj > idp) idlj = idlj - idp
1248 8458159 : do jk = 1_IK, idl1
1249 4652452 : C2(jk, l) = C2(jk, l) + coef(ix1 + idlj)%re * Ch2(jk, j)
1250 8295610 : C2(jk, lc) = C2(jk, lc) + coef(ix1 + idlj)%im * Ch2(jk, jc)
1251 : end do
1252 : end do
1253 : end do
1254 174273 : do j = 2_IK, ipph
1255 494602 : do jk = 1_IK, idl1
1256 482878 : Ch2(jk, 1) = Ch2(jk, 1) + Ch2(jk, j)
1257 : end do
1258 : end do
1259 174273 : do j = 2_IK, ipph
1260 162549 : jc = ipp2 - j
1261 494602 : do jk = 1_IK, idl1
1262 320329 : Ch2(jk, j)%re = C2(jk, j)%re - C2(jk, jc)%im
1263 320329 : Ch2(jk, j)%im = C2(jk, j)%im + C2(jk, jc)%re
1264 320329 : Ch2(jk, jc)%re = C2(jk, j)%re + C2(jk, jc)%im
1265 482878 : Ch2(jk, jc)%im = C2(jk, j)%im - C2(jk, jc)%re
1266 : end do
1267 : end do
1268 11724 : nac = 1_IK
1269 11724 : if (ido == 1_IK) return
1270 499 : nac = 0_IK
1271 5348 : do jk = 1_IK, idl1
1272 5348 : C2(jk, 1) = Ch2(jk, 1)
1273 : end do
1274 3493 : do j = 2_IK, ip
1275 6487 : do k = 1_IK, l1
1276 5988 : C1(1, k, j) = work(1, k, j)
1277 : end do
1278 : end do
1279 499 : if (ido > l1) then
1280 : idj = 1_IK - ido
1281 3493 : do j = 2_IK, ip
1282 2994 : idj = idj + ido
1283 6487 : do k = 1_IK, l1
1284 : idij = idj
1285 32088 : do i = 2_IK, ido
1286 26100 : idij = idij + 1_IK
1287 26100 : C1(i, k, j)%re = coef(ix1 + idij)%re * work(i, k, j)%re - coef(ix1 + idij)%im * work(i, k, j)%im
1288 29094 : C1(i, k, j)%im = coef(ix1 + idij)%re * work(i, k, j)%im + coef(ix1 + idij)%im * work(i, k, j)%re
1289 : end do
1290 : end do
1291 : end do
1292 : else
1293 : idij = 0_IK
1294 0 : do j = 2_IK, ip
1295 0 : idij = idij + 1_IK
1296 0 : do i = 2_IK, ido
1297 0 : idij = idij + 2_IK
1298 0 : do k = 1_IK, l1
1299 0 : C1(i, k, j)%re = coef(ix1 + idij)%re * work(i, k, j)%re - coef(ix1 + idij)%im * work(i, k, j)%im
1300 0 : C1(i, k, j)%im = coef(ix1 + idij)%re * work(i, k, j)%im + coef(ix1 + idij)%im * work(i, k, j)%re
1301 : end do
1302 : end do
1303 : end do
1304 : end if
1305 : end subroutine passb
1306 :
1307 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1308 : #elif setFFTR_ENABLED && RK_ENABLED
1309 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1310 :
1311 : integer(IK) :: lenData, idl1, ido, ip, ix1, ix2, ix3, ix4, k1, l1, l2, na, lenFactor
1312 14749 : lenFactor = size(factor, kind = IK)
1313 14749 : lenData = size(data, kind = IK)
1314 44247 : CHECK_ASSERTION(__LINE__, lenData == size(coef, kind = IK), SK_"@setFFTR(): The condition `size(data) == size(coef)` must hold. size(data), size(coef) = "//getStr([lenData, size(coef, kind = IK)]))
1315 44247 : CHECK_ASSERTION(__LINE__, lenData == size(work, kind = IK), SK_"@setFFTR(): The condition `size(data) == size(work)` must hold. size(data), size(work) = "//getStr([lenData, size(work, kind = IK)]))
1316 14749 : inwork = logical(1_IK < lenData, LK)
1317 14749 : if (.not. inwork) return
1318 : na = 0_IK
1319 14749 : l1 = 1_IK
1320 14749 : ix1 = 1_IK
1321 45436 : do k1 = 1_IK, lenFactor
1322 30687 : ip = factor(k1)
1323 30687 : l2 = ip * l1
1324 30687 : ido = lenData / l2
1325 30687 : idl1 = ido * l1
1326 30687 : if (ip == 4_IK) then
1327 4074 : ix2 = ix1 + ido
1328 4074 : ix3 = ix2 + ido
1329 4074 : if (na /= 0_IK) then
1330 1579 : call radb4(ido, l1, ix1, ix2, ix3, coef, work, data)
1331 : else
1332 2495 : call radb4(ido, l1, ix1, ix2, ix3, coef, data, work)
1333 : end if
1334 4074 : na = 1_IK - na
1335 26613 : elseif (ip == 2_IK) then
1336 3781 : if (na /= 0_IK) then
1337 0 : call radb2(ido, l1, ix1, coef, work, data)
1338 : else
1339 3781 : call radb2(ido, l1, ix1, coef, data, work)
1340 : end if
1341 3781 : na = 1_IK - na
1342 22832 : elseif (ip == 3_IK) then
1343 7664 : ix2 = ix1 + ido
1344 7664 : if (na /= 0_IK) then
1345 3377 : call radb3(ido, l1, ix1, ix2, coef, work, data)
1346 : else
1347 4287 : call radb3(ido, l1, ix1, ix2, coef, data, work)
1348 : end if
1349 7664 : na = 1_IK - na
1350 15168 : elseif (ip /= 5_IK) then
1351 11570 : if (na /= 0_IK) then
1352 4752 : call radbg(ido, ip, l1, idl1, ix1, coef, work, work, work, data, data)
1353 : else
1354 6818 : call radbg(ido, ip, l1, idl1, ix1, coef, data, data, data, work, work)
1355 : end if
1356 11570 : if (ido == 1_IK) na = 1_IK - na
1357 : else
1358 3598 : ix2 = ix1 + ido
1359 3598 : ix3 = ix2 + ido
1360 3598 : ix4 = ix3 + ido
1361 3598 : if (na /= 0_IK) then
1362 1236 : call radb5(ido, l1, ix1, ix2, ix3, ix4, coef, work, data)
1363 : else
1364 2362 : call radb5(ido, l1, ix1, ix2, ix3, ix4, coef, data, work)
1365 : end if
1366 3598 : na = 1_IK - na
1367 : end if
1368 30687 : l1 = l2
1369 45436 : ix1 = ix1 + (ip - 1_IK) * ido
1370 : end do
1371 14749 : inwork = logical(na /= 0_IK, LK)
1372 : !if (na == 0_IK) return
1373 : !data(1 : lenData) = work(1 : lenData)
1374 :
1375 : contains
1376 :
1377 3781 : pure subroutine radb2(ido, l1, ix1, coef, data, work)
1378 : integer(IK) , intent(in) :: ido, l1, ix1
1379 : real(TKC) , intent(in), contiguous :: coef(2:)
1380 : real(TKC) , intent(in) :: data(ido, 2, l1)
1381 : real(TKC) , intent(out) :: work(ido, l1, 2)
1382 : integer(IK) :: i, ic, idp2, k
1383 : complex(TKC) :: t2
1384 7562 : do k = 1_IK, l1
1385 3781 : work(1, k, 1) = data(1, 1, k) + data(ido, 2, k)
1386 7562 : work(1, k, 2) = data(1, 1, k) - data(ido, 2, k)
1387 : end do
1388 3781 : if (ido < 2_IK) return
1389 3781 : if (ido /= 2_IK) then
1390 3781 : idp2 = ido + 2_IK
1391 7562 : do k = 1_IK, l1
1392 7562 : do i = 3_IK, ido, 2_IK
1393 55179 : ic = idp2 - i
1394 55179 : work(i - 1, k, 1) = data(i - 1, 1, k) + data(ic - 1, 2, k)
1395 55179 : t2%re = data(i - 1, 1, k) - data(ic - 1, 2, k)
1396 55179 : work(i, k, 1) = data(i, 1, k) - data(ic, 2, k)
1397 55179 : t2%im = data(i, 1, k) + data(ic, 2, k)
1398 55179 : work(i - 1, k, 2) = coef(ix1 + i - 2) * t2%re - coef(ix1 + i - 1) * t2%im
1399 55179 : work(i, k, 2) = coef(ix1 + i - 2) * t2%im + coef(ix1 + i - 1) * t2%re
1400 : end do
1401 : end do
1402 3781 : if (mod(ido, 2_IK) == 1_IK) return
1403 : end if
1404 2078 : do k = 1_IK, l1
1405 1039 : work(ido, k, 1) = data(ido, 1, k) + data(ido, 1, k)
1406 2078 : work(ido, k, 2) = -(data(1, 2, k) + data(1, 2, k))
1407 : end do
1408 : end subroutine
1409 :
1410 7664 : pure subroutine radb3(ido, l1, ix1, ix2, coef, data, work)
1411 : integer(IK) , intent(in) :: ido, l1, ix1, ix2
1412 : real(TKC) , intent(in), contiguous :: coef(2:)
1413 : real(TKC) , intent(in) :: data(ido, 3, l1)
1414 : real(TKC) , intent(out) :: work(ido, l1, 3)
1415 : complex(TKC), parameter :: TAU = cmplx(-.5_TKC, sqrt(3._TKC) / 2._TKC, TKC)
1416 : complex(TKC) :: c2, c3, d2, d3, t2
1417 : integer(IK) :: i, ic, idp2, k
1418 45681 : do k = 1_IK, l1
1419 38017 : t2%re = data(ido, 2, k) + data(ido, 2, k)
1420 38017 : c2%re = data(1, 1, k) + TAU%re * t2%re
1421 38017 : work(1, k, 1) = data(1, 1, k) + t2%re
1422 38017 : c3%im = TAU%im * (data(1, 3, k) + data(1, 3, k))
1423 38017 : work(1, k, 2) = c2%re - c3%im
1424 45681 : work(1, k, 3) = c2%re + c3%im
1425 : end do
1426 7664 : if (ido == 1_IK) return
1427 6388 : idp2 = ido + 2_IK
1428 23316 : do k = 1_IK, l1
1429 23316 : do i = 3_IK, ido, 2_IK
1430 56638 : ic = idp2 - i
1431 56638 : t2%re = data(i - 1, 3, k) + data(ic - 1, 2, k)
1432 56638 : c2%re = data(i - 1, 1, k) + TAU%re * t2%re
1433 56638 : work(i - 1, k, 1) = data(i - 1, 1, k) + t2%re
1434 56638 : t2%im = data(i, 3, k) - data(ic, 2, k)
1435 56638 : c2%im = data(i, 1, k) + TAU%re * t2%im
1436 56638 : work(i, k, 1) = data(i, 1, k) + t2%im
1437 56638 : c3%re = TAU%im * (data(i - 1, 3, k) - data(ic - 1, 2, k))
1438 56638 : c3%im = TAU%im * (data(i, 3, k) + data(ic, 2, k))
1439 56638 : d2%re = c2%re - c3%im
1440 56638 : d3%re = c2%re + c3%im
1441 56638 : d2%im = c2%im + c3%re
1442 56638 : d3%im = c2%im - c3%re
1443 56638 : work(i - 1, k, 2) = coef(ix1 + i - 2) * d2%re - coef(ix1 + i - 1) * d2%im
1444 56638 : work(i, k, 2) = coef(ix1 + i - 2) * d2%im + coef(ix1 + i - 1) * d2%re
1445 56638 : work(i - 1, k, 3) = coef(ix2 + i - 2) * d3%re - coef(ix2 + i - 1) * d3%im
1446 56638 : work(i, k, 3) = coef(ix2 + i - 2) * d3%im + coef(ix2 + i - 1) * d3%re
1447 : end do
1448 : end do
1449 : end subroutine
1450 :
1451 4074 : pure subroutine radb4(ido, l1, ix1, ix2, ix3, coef, data, work)
1452 : integer(IK) , intent(in) :: ido, l1, ix1, ix2, ix3
1453 : real(TKC) , intent(in), contiguous :: coef(2:)
1454 : real(TKC) , intent(in) :: data(ido, 4, l1)
1455 : real(TKC) , intent(out) :: work(ido, l1, 4)
1456 : real(TKC) , parameter :: NSQT2 = -sqrt(2._TKC)
1457 : complex(TKC) :: c2, c3, c4, t1, t2, t3, t4
1458 : integer(IK) :: i, ic, idp2, k
1459 14959 : do k = 1_IK, l1
1460 10885 : t1%re = data(1, 1, k) - data(ido, 4, k)
1461 10885 : t2%re = data(1, 1, k) + data(ido, 4, k)
1462 10885 : t3%re = data(ido, 2, k) + data(ido, 2, k)
1463 10885 : t4%re = data(1, 3, k) + data(1, 3, k)
1464 10885 : work(1, k, 1) = t2%re + t3%re
1465 10885 : work(1, k, 2) = t1%re - t4%re
1466 10885 : work(1, k, 3) = t2%re - t3%re
1467 14959 : work(1, k, 4) = t1%re + t4%re
1468 : end do
1469 4074 : if (ido < 2_IK) return
1470 3663 : if (ido /= 2_IK) then
1471 3663 : idp2 = ido + 2_IK
1472 9854 : do k = 1_IK, l1
1473 9854 : do i = 3_IK, ido, 2_IK
1474 23712 : ic = idp2 - i
1475 23712 : t1%im = data(i, 1, k) + data(ic, 4, k)
1476 23712 : t2%im = data(i, 1, k) - data(ic, 4, k)
1477 23712 : t3%im = data(i, 3, k) - data(ic, 2, k)
1478 23712 : t4%re = data(i, 3, k) + data(ic, 2, k)
1479 23712 : t1%re = data(i - 1, 1, k) - data(ic - 1, 4, k)
1480 23712 : t2%re = data(i - 1, 1, k) + data(ic - 1, 4, k)
1481 23712 : t4%im = data(i - 1, 3, k) - data(ic - 1, 2, k)
1482 23712 : t3%re = data(i - 1, 3, k) + data(ic - 1, 2, k)
1483 23712 : work(i - 1, k, 1) = t2%re + t3%re
1484 23712 : c3%re = t2%re - t3%re
1485 23712 : work(i, k, 1) = t2%im + t3%im
1486 23712 : c3%im = t2%im - t3%im
1487 23712 : c2%re = t1%re - t4%re
1488 23712 : c4%re = t1%re + t4%re
1489 23712 : c2%im = t1%im + t4%im
1490 23712 : c4%im = t1%im - t4%im
1491 23712 : work(i - 1, k, 2) = coef(ix1 + i - 2) * c2%re - coef(ix1 + i - 1) * c2%im
1492 23712 : work(i, k, 2) = coef(ix1 + i - 2) * c2%im + coef(ix1 + i - 1) * c2%re
1493 23712 : work(i - 1, k, 3) = coef(ix2 + i - 2) * c3%re - coef(ix2 + i - 1) * c3%im
1494 23712 : work(i, k, 3) = coef(ix2 + i - 2) * c3%im + coef(ix2 + i - 1) * c3%re
1495 23712 : work(i - 1, k, 4) = coef(ix3 + i - 2) * c4%re - coef(ix3 + i - 1) * c4%im
1496 23712 : work(i, k, 4) = coef(ix3 + i - 2) * c4%im + coef(ix3 + i - 1) * c4%re
1497 : end do
1498 : end do
1499 3663 : if (mod(ido, 2_IK) == 1_IK) return
1500 : end if
1501 2528 : do k = 1_IK, l1
1502 1660 : t1%im = data(1, 2, k) + data(1, 4, k)
1503 1660 : t2%im = data(1, 4, k) - data(1, 2, k)
1504 1660 : t1%re = data(ido, 1, k) - data(ido, 3, k)
1505 1660 : t2%re = data(ido, 1, k) + data(ido, 3, k)
1506 1660 : work(ido, k, 1) = t2%re + t2%re
1507 1660 : work(ido, k, 2) = NSQT2 * (t1%im - t1%re)
1508 1660 : work(ido, k, 3) = t2%im + t2%im
1509 2528 : work(ido, k, 4) = NSQT2 * (t1%re + t1%im)
1510 : end do
1511 : end subroutine radb4
1512 :
1513 3598 : pure subroutine radb5(ido, l1, ix1, ix2, ix3, ix4, coef, data, work)
1514 : integer(IK) , intent(in) :: ido, l1, ix1, ix2, ix3, ix4
1515 : real(TKC) , intent(in), contiguous :: coef(2:)
1516 : real(TKC) , intent(in) :: data(ido, 5, l1)
1517 : real(TKC) , intent(out) :: work(ido, l1, 5)
1518 : real(TKC) , parameter :: PI = acos(-1._TKC)
1519 : complex(TKC), parameter :: T11 = cmplx(cos(2._TKC * PI / 5._TKC), sin(2._TKC * PI / 5._TKC), TKC)
1520 : complex(TKC), parameter :: T12 = cmplx(cos(4._TKC * PI / 5._TKC), sin(4._TKC * PI / 5._TKC), TKC)
1521 : complex(TKC) :: c2, c3, c4, c5, d2, d3, d4, d5, t2, t3, t4, t5
1522 : integer(IK) :: i, ic, idp2, k
1523 26656 : do k = 1_IK, l1
1524 23058 : t5%im = data(1, 3, k) + data(1, 3, k)
1525 23058 : t4%im = data(1, 5, k) + data(1, 5, k)
1526 23058 : t2%re = data(ido, 2, k) + data(ido, 2, k)
1527 23058 : t3%re = data(ido, 4, k) + data(ido, 4, k)
1528 23058 : work(1, k, 1) = data(1, 1, k) + t2%re + t3%re
1529 23058 : c2%re = data(1, 1, k) + T11%re * t2%re + T12%re * t3%re
1530 23058 : c3%re = data(1, 1, k) + T12%re * t2%re + T11%re * t3%re
1531 23058 : c5%im = T11%im * t5%im + T12%im * t4%im
1532 23058 : c4%im = T12%im * t5%im - T11%im * t4%im
1533 23058 : work(1, k, 2) = c2%re - c5%im
1534 23058 : work(1, k, 3) = c3%re - c4%im
1535 23058 : work(1, k, 4) = c3%re + c4%im
1536 26656 : work(1, k, 5) = c2%re + c5%im
1537 : end do
1538 3598 : if (ido == 1_IK) return
1539 1693 : idp2 = ido + 2_IK
1540 4292 : do k = 1_IK, l1
1541 4292 : do i = 3_IK, ido, 2_IK
1542 9194 : ic = idp2 - i
1543 9194 : t5%im = data(i, 3, k) + data(ic, 2, k)
1544 9194 : t2%im = data(i, 3, k) - data(ic, 2, k)
1545 9194 : t4%im = data(i, 5, k) + data(ic, 4, k)
1546 9194 : t3%im = data(i, 5, k) - data(ic, 4, k)
1547 9194 : t5%re = data(i - 1, 3, k) - data(ic - 1, 2, k)
1548 9194 : t2%re = data(i - 1, 3, k) + data(ic - 1, 2, k)
1549 9194 : t4%re = data(i - 1, 5, k) - data(ic - 1, 4, k)
1550 9194 : t3%re = data(i - 1, 5, k) + data(ic - 1, 4, k)
1551 9194 : work(i - 1, k, 1) = data(i - 1, 1, k) + t2%re + t3%re
1552 9194 : work(i, k, 1) = data(i, 1, k) + t2%im + t3%im
1553 9194 : c2%re = data(i - 1, 1, k) + T11%re * t2%re + T12%re * t3%re
1554 9194 : c2%im = data(i, 1, k) + T11%re * t2%im + T12%re * t3%im
1555 9194 : c3%re = data(i - 1, 1, k) + T12%re * t2%re + T11%re * t3%re
1556 9194 : c3%im = data(i, 1, k) + T12%re * t2%im + T11%re * t3%im
1557 9194 : c5%re = T11%im * t5%re + T12%im * t4%re
1558 9194 : c5%im = T11%im * t5%im + T12%im * t4%im
1559 9194 : c4%re = T12%im * t5%re - T11%im * t4%re
1560 9194 : c4%im = T12%im * t5%im - T11%im * t4%im
1561 9194 : d3%re = c3%re - c4%im
1562 9194 : d4%re = c3%re + c4%im
1563 9194 : d3%im = c3%im + c4%re
1564 9194 : d4%im = c3%im - c4%re
1565 9194 : d5%re = c2%re + c5%im
1566 9194 : d2%re = c2%re - c5%im
1567 9194 : d5%im = c2%im - c5%re
1568 9194 : d2%im = c2%im + c5%re
1569 9194 : work(i - 1, k, 2) = coef(ix1 + i - 2) * d2%re - coef(ix1 + i - 1) * d2%im
1570 9194 : work(i, k, 2) = coef(ix1 + i - 2) * d2%im + coef(ix1 + i - 1) * d2%re
1571 9194 : work(i - 1, k, 3) = coef(ix2 + i - 2) * d3%re - coef(ix2 + i - 1) * d3%im
1572 9194 : work(i, k, 3) = coef(ix2 + i - 2) * d3%im + coef(ix2 + i - 1) * d3%re
1573 9194 : work(i - 1, k, 4) = coef(ix3 + i - 2) * d4%re - coef(ix3 + i - 1) * d4%im
1574 9194 : work(i, k, 4) = coef(ix3 + i - 2) * d4%im + coef(ix3 + i - 1) * d4%re
1575 9194 : work(i - 1, k, 5) = coef(ix4 + i - 2) * d5%re - coef(ix4 + i - 1) * d5%im
1576 9194 : work(i, k, 5) = coef(ix4 + i - 2) * d5%im + coef(ix4 + i - 1) * d5%re
1577 : end do
1578 : end do
1579 : end subroutine
1580 :
1581 11570 : pure subroutine radbg(ido, ip, l1, idl1, ix1, coef, data, C1, C2, work, Ch2)
1582 : integer(IK) , intent(in) :: ido, ip, l1, idl1, ix1
1583 : real(TKC) , intent(in), contiguous :: coef(2:)
1584 : real(TKC) , intent(in) :: data(ido,ip,l1)
1585 : real(TKC) , intent(inout) :: C1(ido,l1,ip)
1586 : real(TKC) , intent(inout) :: work(ido,l1,ip), C2(idl1,ip), Ch2(idl1,ip)
1587 : real(TKC) , parameter :: TWO_PI = 2._TKC * acos(-1._TKC)
1588 : real(TKC) :: ar1h, ar2h, arg, dc2, dcp, ds2, dsp
1589 : integer(IK) :: i, ic, idij, idp2, jk, ipp2, ipph, is, j, j2, jc, k, l, lc, nbd
1590 : complex(TKC) :: a1, a2
1591 11570 : arg = TWO_PI / ip
1592 11570 : dcp = cos(arg)
1593 11570 : dsp = sin(arg)
1594 11570 : idp2 = ido + 2_IK
1595 11570 : nbd = (ido - 1_IK) / 2_IK
1596 11570 : ipp2 = ip + 2
1597 11570 : ipph = (ip + 1_IK) / 2_IK
1598 11570 : if (ido < l1) then
1599 13782 : do i = 1_IK, ido
1600 47113 : do k = 1_IK, l1
1601 40222 : work(i, k, 1) = data(i, 1, k)
1602 : end do
1603 : end do
1604 : else
1605 9359 : do k = 1_IK, l1
1606 17537 : do i = 1_IK, ido
1607 12858 : work(i, k, 1) = data(i, 1, k)
1608 : end do
1609 : end do
1610 : end if
1611 209731 : do j = 2_IK, ipph
1612 198161 : jc = ipp2 - j
1613 198161 : j2 = j + j
1614 550636 : do k = 1_IK, l1
1615 340905 : work(1, k, j) = data(ido, j2 - 2, k) + data(ido, j2 - 2, k)
1616 539066 : work(1, k, jc) = data(1, j2 - 1, k) + data(1, j2 - 1, k)
1617 : end do
1618 : end do
1619 11570 : if (ido /= 1_IK) then
1620 413 : if (nbd < l1) then
1621 0 : do j = 2_IK, ipph
1622 0 : jc = ipp2 - j
1623 0 : do i = 3_IK, ido, 2_IK
1624 0 : ic = idp2 - i
1625 0 : do k = 1_IK, l1
1626 0 : work(i - 1, k, j) = data(i - 1, 2 * j - 1, k) + data(ic - 1, 2 * j - 2, k)
1627 0 : work(i - 1, k, jc) = data(i - 1, 2 * j - 1, k) - data(ic - 1, 2 * j - 2, k)
1628 0 : work(i, k, j) = data(i, 2 * j - 1, k) - data(ic, 2 * j - 2, k)
1629 0 : work(i, k, jc) = data(i, 2 * j - 1, k) + data(ic, 2 * j - 2, k)
1630 : end do
1631 : end do
1632 : end do
1633 : else
1634 1652 : do j = 2_IK, ipph
1635 1239 : jc = ipp2 - j
1636 2894 : do k = 1_IK, l1
1637 2481 : do i = 3_IK, ido, 2_IK
1638 5247 : ic = idp2 - i
1639 5247 : work(i - 1, k, j) = data(i - 1, 2 * j - 1, k) + data(ic - 1, 2 * j - 2, k)
1640 5247 : work(i - 1, k, jc) = data(i - 1, 2 * j - 1, k) - data(ic - 1, 2 * j - 2, k)
1641 5247 : work(i, k, j) = data(i, 2 * j - 1, k) - data(ic, 2 * j - 2, k)
1642 5247 : work(i, k, jc) = data(i, 2 * j - 1, k) + data(ic, 2 * j - 2, k)
1643 : end do
1644 : end do
1645 : end do
1646 : end if
1647 : end if
1648 8583 : a1%re = 1._TKC
1649 8583 : a1%im = 0._TKC
1650 209731 : do l = 2_IK, ipph
1651 198161 : lc = ipp2 - l
1652 198161 : ar1h = dcp * a1%re - dsp * a1%im
1653 198161 : a1%im = dcp * a1%im + dsp * a1%re
1654 198161 : a1%re = ar1h
1655 549560 : do jk = 1_IK, idl1
1656 351399 : C2(jk, l) = Ch2(jk, 1) + a1%re * Ch2(jk, 2)
1657 549560 : C2(jk, lc) = a1%im * Ch2(jk, ip)
1658 : end do
1659 78794 : dc2 = a1%re
1660 78794 : ds2 = a1%im
1661 121455 : a2%re = a1%re
1662 121455 : a2%im = a1%im
1663 345841881 : do j = 3_IK, ipph
1664 345632150 : jc = ipp2 - j
1665 345632150 : ar2h = dc2 * a2%re - ds2 * a2%im
1666 345632150 : a2%im = dc2 * a2%im + ds2 * a2%re
1667 345632150 : a2%re = ar2h
1668 692405197 : do jk = 1_IK, idl1
1669 346574886 : C2(jk, l) = C2(jk, l) + a2%re * Ch2(jk, j)
1670 692207036 : C2(jk, lc) = C2(jk, lc) + a2%im * Ch2(jk, jc)
1671 : end do
1672 : end do
1673 : end do
1674 209731 : do j = 2_IK, ipph
1675 561130 : do jk = 1_IK, idl1
1676 549560 : Ch2(jk, 1) = Ch2(jk, 1) + Ch2(jk, j)
1677 : end do
1678 : end do
1679 209731 : do j = 2_IK, ipph
1680 198161 : jc = ipp2 - j
1681 550636 : do k = 1_IK, l1
1682 340905 : work(1, k, j) = C1(1, k, j) - C1(1, k, jc)
1683 539066 : work(1, k, jc) = C1(1, k, j) + C1(1, k, jc)
1684 : end do
1685 : end do
1686 11570 : if (ido /= 1_IK) then
1687 413 : if (nbd < l1) then
1688 0 : do j = 2_IK, ipph
1689 0 : jc = ipp2 - j
1690 0 : do i = 3_IK, ido, 2_IK
1691 0 : do k = 1_IK, l1
1692 0 : work(i - 1, k, j) = C1(i - 1, k, j) - C1(i, k, jc)
1693 0 : work(i - 1, k, jc) = C1(i - 1, k, j) + C1(i, k, jc)
1694 0 : work(i, k, j) = C1(i, k, j) + C1(i - 1, k, jc)
1695 0 : work(i, k, jc) = C1(i, k, j) - C1(i - 1, k, jc)
1696 : end do
1697 : end do
1698 : end do
1699 : else
1700 1652 : do j = 2_IK, ipph
1701 1239 : jc = ipp2 - j
1702 2894 : do k = 1_IK, l1
1703 2481 : do i = 3_IK, ido, 2_IK
1704 5247 : work(i - 1, k, j) = C1(i - 1, k, j) - C1(i, k, jc)
1705 5247 : work(i - 1, k, jc) = C1(i - 1, k, j) + C1(i, k, jc)
1706 5247 : work(i, k, j) = C1(i, k, j) + C1(i - 1, k, jc)
1707 5247 : work(i, k, jc) = C1(i, k, j) - C1(i - 1, k, jc)
1708 : end do
1709 : end do
1710 : end do
1711 : end if
1712 : end if
1713 11570 : if (ido == 1_IK) return
1714 4325 : C2(1 : idl1, 1) = Ch2(1 : idl1, 1)
1715 5375 : C1(1, 1 : l1, 2 : ip) = work(1, 1 : l1, 2 : ip)
1716 413 : if (nbd > l1) then
1717 413 : is = -ido
1718 2891 : do j = 2_IK, ip
1719 2478 : is = is + ido
1720 5375 : do k = 1_IK, l1
1721 : idij = is
1722 4962 : do i = 3_IK, ido, 2_IK
1723 10494 : idij = idij + 2_IK
1724 10494 : C1(i - 1, k, j) = coef(ix1 + idij - 1) * work(i - 1, k, j) - coef(ix1 + idij) * work(i, k, j)
1725 10494 : C1(i, k, j) = coef(ix1 + idij - 1) * work(i, k, j) + coef(ix1 + idij) * work(i - 1, k, j)
1726 : end do
1727 : end do
1728 : end do
1729 : else
1730 0 : is = -ido
1731 0 : do j = 2_IK, ip
1732 0 : is = is + ido
1733 : idij = is
1734 0 : do i = 3_IK, ido, 2_IK
1735 0 : idij = idij + 2_IK
1736 0 : do k = 1_IK, l1
1737 0 : C1(i - 1, k, j) = coef(ix1 + idij - 1) * work(i - 1, k, j) - coef(ix1 + idij) * work(i, k, j)
1738 0 : C1(i, k, j) = coef(ix1 + idij - 1) * work(i, k, j) + coef(ix1 + idij) * work(i - 1, k, j)
1739 : end do
1740 : end do
1741 : end do
1742 : end if
1743 : end subroutine
1744 :
1745 : !%%%%%%%%%%%%%%
1746 : #elif setFFTI_ENABLED
1747 : !%%%%%%%%%%%%%%
1748 :
1749 : integer(IK) :: i
1750 : real(TKC) :: invLenData
1751 309 : call setFFTR(factor, coef, data, work, inwork)
1752 309 : invLenData = 1._TKC / real(size(data, 1, IK), TKC)
1753 309 : if (inwork) then
1754 : do concurrent(i = 1 : size(data, 1, IK))
1755 5866 : work(i) = work(i) * invLenData
1756 : end do
1757 : else
1758 : do concurrent(i = 1 : size(data, 1, IK))
1759 4710 : data(i) = data(i) * invLenData
1760 : end do
1761 : end if
1762 :
1763 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1764 : #elif getFFTF_ENABLED || getFFTR_ENABLED || getFFTI_ENABLED
1765 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1766 :
1767 : #if CK_ENABLED
1768 699 : complex(TKC), dimension(size(data, 1, IK)) :: coef, work
1769 : #elif RK_ENABLED
1770 165 : real(TKC), dimension(size(data, 1, IK)) :: coef, work
1771 : #else
1772 : #error "Unrecognized interface."
1773 : #endif
1774 : logical(LK) :: inwork
1775 : integer(IK) , allocatable :: factor(:)
1776 3652 : factor = getFactorFFT(data, coef)
1777 31824 : fft = data
1778 : #if getFFTF_ENABLED
1779 531 : call setFFTF(factor, coef, fft, work, inwork)
1780 : #elif getFFTR_ENABLED
1781 30 : call setFFTR(factor, coef, fft, work, inwork)
1782 : #elif getFFTI_ENABLED
1783 303 : call setFFTI(factor, coef, fft, work, inwork)
1784 : #else
1785 : #error "Unrecognized interface."
1786 : #endif
1787 18872 : if (inwork) fft = work
1788 : ! Normalize the result.
1789 : !#if getFFTF_ENABLED || getFFTR_ENABLED
1790 : ! if (inwork) fft = work
1791 : !#elif getFFTI_ENABLED
1792 : ! block
1793 : ! integer(IK) :: i
1794 : ! real(TKC) :: invLenData
1795 : ! invLenData = 1._TKC / real(size(data, 1, IK), TKC)
1796 : ! if (inwork) then
1797 : ! do concurrent(i = 1 : size(data, 1, IK))
1798 : ! fft(i) = fft(i) * invLenData
1799 : ! end do
1800 : ! else
1801 : ! do concurrent(i = 1 : size(data, 1, IK))
1802 : ! fft(i) = work(i) * invLenData
1803 : ! end do
1804 : ! end if
1805 : ! end block
1806 : !#endif
1807 :
1808 : #else
1809 : !%%%%%%%%%%%%%%%%%%%%%%%%
1810 : #error "Unrecognized interface."
1811 : !%%%%%%%%%%%%%%%%%%%%%%%%
1812 : #endif
|
