Line data Source code
1 : #ifndef CPP_ESSL
2 : #define CPP_Singleton
3 : #endif
4 : ! ou might want to undefine CPP_singleton to use essl FFT
5 : MODULE m_cfft
6 : IMPLICIT NONE
7 : CONTAINS
8 : ! ***************************************************************
9 : ! multivariate complex fourier transform, computed in place
10 : ! using mixed-radix fast fourier transform algorithm.
11 : ! by r. c. singleton, stanford research institute, oct. 1968
12 : ! arrays a and b originally hold the real and imaginary
13 : ! components of the data, and return the real and
14 : ! imaginary components of the resulting fourier coefficients.
15 : ! multivariate data is indexed according to the fortran
16 : ! array element successor function, without limit
17 : ! on the number of implied multiple subscripts.
18 : ! the subroutine is called once for each variate.
19 : ! the calls for a multivariate transform may be in any order.
20 : ! ntot is the total number of complex data values.
21 : ! n is the dimension of the current variable.
22 : ! nspan/n is the spacing of consucutive data values
23 : ! while indexing the current variable.
24 : ! the sign of isn determines the sign of the complex
25 : ! exponential, and the magnitude of isn is normally one.
26 : ! for a single-variate transform,
27 : ! ntot = n = nspan = (number of complex data values), f.g.
28 : ! call cft(a,b,n,n,n,1)
29 : ! a tri-variate transform with a(n1,n2,n3), b(n1,n2,n3)
30 : ! is computed by
31 : ! call cft(a,b,n1*n2*n3,n1,n1,1)
32 : ! call cft(a,b,n1*n2*n3,n2,n1*n2,1)
33 : ! call cft(a,b,n1*n2*n3,n3,n1*n2*n3,1)
34 : ! the data may alternatively be stored in a single complex
35 : ! array a, then the magnitude of isn changed to two to
36 : ! give the correct indexing increment and the second parameter
37 : ! used to pass the initial address for the sequence of
38 : ! imaginary values, e.g.
39 : ! real s(2)
40 : ! equivalence (a,s)
41 : ! ....
42 : ! ....
43 : ! call cft(a,s(2),ntot,n,nspan,2)
44 : ! arrays at(maxf), ck(maxf), bt(maxf), sk(maxf), and np(maxp)
45 : ! are used for temporary storage. if the available storage
46 : ! is insufficient, the program is terminated by a stop.
47 : ! maxf must be .ge. the maximum prime factor of n.
48 : ! maxp must be .gt. the number of prime factors of n.
49 : ! in addition, if the square-free portion k of n has two or
50 : ! more prime factors, then maxp must be .ge. k-1.
51 : ! array storage in nfac for a maximum of 11 factors of n.
52 : ! if n has more than one square-free factor, the product of the
53 : ! square-free factors must be .le. 210
54 : ! *******************************************************************
55 : ! array storage for maximum prime factor of 199
56 : ! the following two constants should agree with the array dimensions
57 : #ifdef CPP_Singleton
58 14812 : SUBROUTINE cfft(a, b, ntot, n, nspan, isn)
59 : USE m_juDFT
60 : USE m_constants
61 : IMPLICIT NONE
62 : ! .. Scalar Arguments ..
63 : INTEGER :: isn, n, nspan, ntot
64 : ! ..
65 : ! .. Array Arguments ..
66 : ! REAL a(*),b(*)
67 : REAL :: a(ntot), b(ntot)
68 : ! ..
69 : ! .. Local Scalars ..
70 : REAL :: aa, aj, ajm, ajp, ak, akm, akp, bb, bj, bjm, bjp, bk, bkm, bkp, c1, c2, c3, &
71 : c72, cd, rad, radf, s1, s120, s2, s3, s72, sd
72 : INTEGER :: i, ii, inc, j, jc, jf, jj, k, k1, k2, k3, k4, kk, ks, kspan, kspnn, kt, m, &
73 : maxf, maxp, nn, nt, maxnf
74 : ! ..
75 : ! .. Local Arrays ..
76 14812 : REAL, ALLOCATABLE :: at(:), bt(:), ck(:), sk(:)
77 14812 : INTEGER, ALLOCATABLE :: nfac(:), np(:)
78 : ! ..
79 : ! .. Equivalences ..
80 : EQUIVALENCE(i, ii)
81 : ! ..
82 14812 : IF (n < 2) RETURN
83 :
84 14812 : maxf = MAX(299, n)
85 14812 : maxp = MAX(503, n - 1)
86 14812 : maxnf = MAX(17, CEILING(sqrt(1.0*n)))
87 88872 : ALLOCATE (at(maxf), bt(maxf), ck(maxf), sk(maxf))
88 74060 : ALLOCATE (nfac(maxnf), np(maxp))
89 :
90 14812 : inc = isn
91 : ! the following constants are rad = 2.*pi , s72 = sin(0.4*pi) ,
92 : ! c72 = cos(0.4*pi) and s120 = sqrt(0.75)
93 14812 : rad = 6.2831853071796
94 14812 : s72 = 0.95105651629515
95 14812 : c72 = 0.30901699437495
96 14812 : s120 = 0.86602540378444
97 14812 : IF (isn < 0) THEN
98 14572 : s72 = -s72
99 14572 : s120 = -s120
100 14572 : rad = -rad
101 14572 : inc = -inc
102 : ENDIF
103 14812 : nt = inc*ntot
104 14812 : ks = inc*nspan
105 14812 : kspan = ks
106 14812 : nn = nt - inc
107 14812 : jc = ks/n
108 14812 : radf = rad*real(jc)*0.5
109 : i = 0
110 14812 : jf = 0
111 : ! determine the factors of n
112 14812 : m = 0
113 14812 : k = n
114 14831 : DO WHILE (k - (k/16)*16 == 0)
115 19 : m = m + 1
116 19 : nfac(m) = 4
117 14831 : k = k/16
118 : ENDDO
119 : j = 3
120 : jj = 9
121 2495 : GO TO 50
122 2495 : 40 m = m + 1
123 2495 : nfac(m) = j
124 2495 : k = k/jj
125 30085 : 50 IF (mod(k,jj).EQ.0) GO TO 40
126 27590 : j = j + 2
127 27590 : jj = j**2
128 27590 : IF (jj.LE.k) GO TO 50
129 14812 : IF (k.GT.4) GO TO 60
130 65 : kt = m
131 65 : nfac(m+1) = k
132 65 : IF (k.NE.1) m = m + 1
133 14747 : GO TO 100
134 14747 : 60 IF (k- (k/4)*4.NE.0) GO TO 70
135 44 : m = m + 1
136 44 : nfac(m) = 2
137 14747 : k = k/4
138 14747 : 70 kt = m
139 14747 : j = 2
140 79349 : 80 IF (mod(k,j).NE.0) GO TO 90
141 29876 : m = m + 1
142 29876 : nfac(m) = j
143 79349 : k = k/j
144 79349 : 90 j = ((j + 1)/2)*2 + 1
145 79349 : IF (j <= k) GO TO 80
146 14812 : 100 IF (kt == 0) GO TO 120
147 5116 : DO j=kt,1,-1
148 2558 : m = m + 1
149 17370 : nfac(m) = nfac(j)
150 : ENDDO
151 : ! compute fourier transform
152 35057 : 120 sd = radf/real(kspan)
153 35057 : cd = 2.0*sin(sd)**2
154 35057 : sd = sin(sd + sd)
155 35057 : kk = 1
156 35057 : i = i + 1
157 35057 : IF (nfac(i) /= 2) GO TO 170
158 : ! transform for factor of 2 (including rotation factor)
159 2963 : kspan = kspan/2
160 2963 : k1 = kspan + 2
161 483036 : 130 k2 = kk + kspan
162 483036 : ak = a(k2)
163 483036 : bk = b(k2)
164 483036 : a(k2) = a(kk) - ak
165 483036 : b(k2) = b(kk) - bk
166 483036 : a(kk) = a(kk) + ak
167 483036 : b(kk) = b(kk) + bk
168 483036 : kk = k2 + kspan
169 483036 : IF (kk <= nn) GO TO 130
170 52645 : kk = kk - nn
171 52645 : IF (kk <= jc) GO TO 130
172 2963 : IF (kk > kspan) GO TO 360
173 45333 : 140 c1 = 1.0 - cd
174 45333 : s1 = sd
175 1722540 : 150 k2 = kk + kspan
176 1722540 : ak = a(kk) - a(k2)
177 1722540 : bk = b(kk) - b(k2)
178 1722540 : a(kk) = a(kk) + a(k2)
179 1722540 : b(kk) = b(kk) + b(k2)
180 1722540 : a(k2) = c1*ak - s1*bk
181 1722540 : b(k2) = s1*ak + c1*bk
182 1722540 : kk = k2 + kspan
183 1722540 : IF (kk < nt) GO TO 150
184 809610 : k2 = kk - nt
185 809610 : c1 = -c1
186 809610 : kk = k1 - k2
187 809610 : IF (kk > k2) GO TO 150
188 412117 : ak = c1 - (cd*c1 + sd*s1)
189 412117 : s1 = (sd*c1 - cd*s1) + s1
190 : ! the following three statements compensate for truncation
191 : ! error. if rounded arithmetic is used, they may be deleted.
192 : ! c1=0.5/(ak**2+s1**2)+0.5
193 : ! s1=c1*s1
194 : ! c1=c1*ak
195 : ! next statement should be deleted if non-rounded arithmetic is used
196 412117 : c1 = ak
197 412117 : kk = kk + jc
198 412117 : IF (kk < k2) GO TO 150
199 45333 : k1 = k1 + inc + inc
200 45333 : kk = (k1 - kspan)/2 + jc
201 45333 : IF (kk <= jc + jc) GO TO 140
202 7653783 : GO TO 120
203 : ! transform for factor of 3 (optional code)
204 7653783 : 160 k1 = kk + kspan
205 7653783 : k2 = k1 + kspan
206 7653783 : ak = a(kk)
207 7653783 : bk = b(kk)
208 7653783 : aj = a(k1) + a(k2)
209 7653783 : bj = b(k1) + b(k2)
210 7653783 : a(kk) = ak + aj
211 7653783 : b(kk) = bk + bj
212 7653783 : ak = -0.5*aj + ak
213 7653783 : bk = -0.5*bj + bk
214 7653783 : aj = (a(k1) - a(k2))*s120
215 7653783 : bj = (b(k1) - b(k2))*s120
216 7653783 : a(k1) = ak - bj
217 7653783 : b(k1) = bk + aj
218 7653783 : a(k2) = ak + bj
219 7653783 : b(k2) = bk - aj
220 7653783 : kk = k2 + kspan
221 7653783 : IF (kk < nn) GO TO 160
222 1431348 : kk = kk - nn
223 1431348 : IF (kk <= kspan) GO TO 160
224 32094 : GO TO 320
225 : ! transform for factor of 4
226 32094 : 170 IF (nfac(i) /= 4) GO TO 260
227 65 : kspnn = kspan
228 65 : kspan = kspan/4
229 49180 : 180 c1 = 1.0
230 49180 : s1 = 0
231 1248534 : 190 k1 = kk + kspan
232 1248534 : k2 = k1 + kspan
233 1248534 : k3 = k2 + kspan
234 1248534 : akp = a(kk) + a(k2)
235 1248534 : akm = a(kk) - a(k2)
236 1248534 : ajp = a(k1) + a(k3)
237 1248534 : ajm = a(k1) - a(k3)
238 1248534 : a(kk) = akp + ajp
239 1248534 : ajp = akp - ajp
240 1248534 : bkp = b(kk) + b(k2)
241 1248534 : bkm = b(kk) - b(k2)
242 1248534 : bjp = b(k1) + b(k3)
243 1248534 : bjm = b(k1) - b(k3)
244 1248534 : b(kk) = bkp + bjp
245 1248534 : bjp = bkp - bjp
246 1248534 : IF (isn < 0) GO TO 220
247 1248534 : akp = akm - bjm
248 1248534 : akm = akm + bjm
249 1248534 : bkp = bkm + ajm
250 1248534 : bkm = bkm - ajm
251 1248534 : IF (s1 == 0.0) GO TO 230
252 641466 : 200 a(k1) = akp*c1 - bkp*s1
253 641466 : b(k1) = akp*s1 + bkp*c1
254 641466 : a(k2) = ajp*c2 - bjp*s2
255 641466 : b(k2) = ajp*s2 + bjp*c2
256 641466 : a(k3) = akm*c3 - bkm*s3
257 641466 : b(k3) = akm*s3 + bkm*c3
258 641466 : kk = k3 + kspan
259 641466 : IF (kk <= nt) GO TO 190
260 287134 : 210 c2 = c1 - (cd*c1 + sd*s1)
261 287134 : s1 = (sd*c1 - cd*s1) + s1
262 : ! the following three statements compensate for truncation
263 : ! error. if rounded arithmetic is used, they may be deleted.
264 : ! c1=0.5/(c2**2+s1**2)+0.5
265 : ! s1=c1*s1
266 : ! c1=c1*c2
267 : ! next statement should be deleted if non-rounded arithmetic is used
268 287134 : c1 = c2
269 287134 : c2 = c1**2 - s1**2
270 287134 : s2 = 2.0*c1*s1
271 287134 : c3 = c2*c1 - s2*s1
272 287134 : s3 = c2*s1 + s2*c1
273 287134 : kk = kk - nt + jc
274 287134 : IF (kk <= kspan) GO TO 190
275 49180 : kk = kk - kspan + inc
276 49180 : IF (kk <= jc) GO TO 180
277 65 : IF (kspan == jc) GO TO 360
278 0 : GO TO 120
279 0 : 220 akp = akm + bjm
280 0 : akm = akm - bjm
281 0 : bkp = bkm - ajm
282 0 : bkm = bkm + ajm
283 0 : IF (s1 /= 0.0) GO TO 200
284 607068 : 230 a(k1) = akp
285 607068 : b(k1) = bkp
286 607068 : a(k2) = ajp
287 607068 : b(k2) = bjp
288 607068 : a(k3) = akm
289 607068 : b(k3) = bkm
290 607068 : kk = k3 + kspan
291 607068 : IF (kk <= nt) GO TO 190
292 161 : GO TO 210
293 : ! transform for factor of 5 (optional code)
294 161 : 240 c2 = c72**2 - s72**2
295 161 : s2 = 2.0*c72*s72
296 1136529 : 250 k1 = kk + kspan
297 1136529 : k2 = k1 + kspan
298 1136529 : k3 = k2 + kspan
299 1136529 : k4 = k3 + kspan
300 1136529 : akp = a(k1) + a(k4)
301 1136529 : akm = a(k1) - a(k4)
302 1136529 : bkp = b(k1) + b(k4)
303 1136529 : bkm = b(k1) - b(k4)
304 1136529 : ajp = a(k2) + a(k3)
305 1136529 : ajm = a(k2) - a(k3)
306 1136529 : bjp = b(k2) + b(k3)
307 1136529 : bjm = b(k2) - b(k3)
308 1136529 : aa = a(kk)
309 1136529 : bb = b(kk)
310 1136529 : a(kk) = aa + akp + ajp
311 1136529 : b(kk) = bb + bkp + bjp
312 1136529 : ak = akp*c72 + ajp*c2 + aa
313 1136529 : bk = bkp*c72 + bjp*c2 + bb
314 1136529 : aj = akm*s72 + ajm*s2
315 1136529 : bj = bkm*s72 + bjm*s2
316 1136529 : a(k1) = ak - bj
317 1136529 : a(k4) = ak + bj
318 1136529 : b(k1) = bk + aj
319 1136529 : b(k4) = bk - aj
320 1136529 : ak = akp*c2 + ajp*c72 + aa
321 1136529 : bk = bkp*c2 + bjp*c72 + bb
322 1136529 : aj = akm*s2 - ajm*s72
323 1136529 : bj = bkm*s2 - bjm*s72
324 1136529 : a(k2) = ak - bj
325 1136529 : a(k3) = ak + bj
326 1136529 : b(k2) = bk + aj
327 1136529 : b(k3) = bk - aj
328 1136529 : kk = k4 + kspan
329 1136529 : IF (kk < nn) GO TO 250
330 106122 : kk = kk - nn
331 106122 : IF (kk <= kspan) GO TO 250
332 32029 : GO TO 320
333 : ! transform for odd factors
334 32029 : 260 k = nfac(i)
335 32029 : kspnn = kspan
336 32029 : kspan = kspan/k
337 32029 : IF (k == 3) GO TO 160
338 12389 : IF (k == 5) GO TO 240
339 12228 : IF (k == jf) GO TO 280
340 12228 : jf = k
341 12228 : s1 = rad/real(k)
342 12228 : c1 = cos(s1)
343 12228 : s1 = sin(s1)
344 12228 : IF (jf > maxf) GO TO 590
345 12228 : ck(jf) = 1.0
346 12228 : sk(jf) = 0.0
347 12228 : j = 1
348 61922 : 270 ck(j) = ck(k)*c1 + sk(k)*s1
349 61922 : sk(j) = ck(k)*s1 - sk(k)*c1
350 61922 : k = k - 1
351 61922 : ck(k) = ck(j)
352 61922 : sk(k) = -sk(j)
353 61922 : j = j + 1
354 61922 : IF (j < k) GO TO 270
355 391059 : 280 k1 = kk
356 391059 : k2 = kk + kspnn
357 391059 : aa = a(kk)
358 391059 : bb = b(kk)
359 391059 : ak = aa
360 391059 : bk = bb
361 391059 : j = 1
362 391059 : k1 = k1 + kspan
363 2243349 : 290 k2 = k2 - kspan
364 2243349 : j = j + 1
365 2243349 : at(j) = a(k1) + a(k2)
366 2243349 : ak = at(j) + ak
367 2243349 : bt(j) = b(k1) + b(k2)
368 2243349 : bk = bt(j) + bk
369 2243349 : j = j + 1
370 2243349 : at(j) = a(k1) - a(k2)
371 2243349 : bt(j) = b(k1) - b(k2)
372 2243349 : k1 = k1 + kspan
373 2243349 : IF (k1 < k2) GO TO 290
374 391059 : a(kk) = ak
375 391059 : b(kk) = bk
376 391059 : k1 = kk
377 391059 : k2 = kk + kspnn
378 391059 : j = 1
379 2243349 : 300 k1 = k1 + kspan
380 2243349 : k2 = k2 - kspan
381 2243349 : jj = j
382 2243349 : ak = aa
383 2243349 : bk = bb
384 2243349 : aj = 0.0
385 2243349 : bj = 0.0
386 2243349 : k = 1
387 14038461 : 310 k = k + 1
388 14038461 : ak = at(k)*ck(jj) + ak
389 14038461 : bk = bt(k)*ck(jj) + bk
390 14038461 : k = k + 1
391 14038461 : aj = at(k)*sk(jj) + aj
392 14038461 : bj = bt(k)*sk(jj) + bj
393 14038461 : jj = jj + j
394 14038461 : IF (jj > jf) jj = jj - jf
395 14038461 : IF (k < jf) GO TO 310
396 2243349 : k = jf - j
397 2243349 : a(k1) = ak - bj
398 2243349 : b(k1) = bk + aj
399 2243349 : a(k2) = ak + bj
400 2243349 : b(k2) = bk - aj
401 2243349 : j = j + 1
402 2243349 : IF (j < k) GO TO 300
403 391059 : kk = kk + kspnn
404 391059 : IF (kk <= nn) GO TO 280
405 37417 : kk = kk - nn
406 37417 : IF (kk <= kspan) GO TO 280
407 : ! multiply by rotation factor (except for factors of 2 and 4)
408 32029 : 320 IF (i == m) GO TO 360
409 17280 : kk = jc + 1
410 168208 : 330 c2 = 1.0 - cd
411 168208 : s1 = sd
412 1308347 : 340 c1 = c2
413 1308347 : s2 = s1
414 1308347 : kk = kk + kspan
415 11833438 : 350 ak = a(kk)
416 11833438 : a(kk) = c2*ak - s2*b(kk)
417 11833438 : b(kk) = s2*ak + c2*b(kk)
418 11833438 : kk = kk + kspnn
419 11833438 : IF (kk <= nt) GO TO 350
420 2746198 : ak = s1*s2
421 2746198 : s2 = s1*c2 + c1*s2
422 2746198 : c2 = c1*c2 - ak
423 2746198 : kk = kk - nt + kspan
424 2746198 : IF (kk <= kspnn) GO TO 350
425 1308347 : c2 = c1 - (cd*c1 + sd*s1)
426 1308347 : s1 = s1 + (sd*c1 - cd*s1)
427 : ! the following three statements compensate for truncation
428 : ! error. if rounded arithmetic is used, they may
429 : ! be deleted.
430 : ! c1=0.5/(c2**2+s1**2)+0.5
431 : ! s1=c1*s1
432 : ! c2=c1*c2
433 1308347 : kk = kk - kspnn + jc
434 1308347 : IF (kk <= kspan) GO TO 340
435 168208 : kk = kk - kspan + jc + inc
436 168208 : IF (kk <= jc + jc) GO TO 330
437 14812 : GO TO 120
438 : ! permute the results to normal order---done in two stages
439 : ! permutation for square factors of n
440 14812 : 360 np(1) = ks
441 14812 : IF (kt == 0) GO TO 450
442 2558 : k = kt + kt + 1
443 2558 : IF (m < k) k = k - 1
444 2558 : j = 1
445 2558 : np(k + 1) = jc
446 2558 : 370 np(j + 1) = np(j)/nfac(j)
447 2558 : np(k) = np(k + 1)*nfac(j)
448 2558 : j = j + 1
449 2558 : k = k - 1
450 2558 : IF (j < k) GO TO 370
451 2558 : k3 = np(k + 1)
452 2558 : kspan = np(2)
453 2558 : kk = jc + 1
454 2558 : k2 = kspan + 1
455 2558 : j = 1
456 2558 : IF (n /= ntot) GO TO 410
457 : ! permutation for single-variate transform (optional code)
458 44040 : 380 ak = a(kk)
459 44040 : a(kk) = a(k2)
460 44040 : a(k2) = ak
461 44040 : bk = b(kk)
462 44040 : b(kk) = b(k2)
463 44040 : b(k2) = bk
464 44040 : kk = kk + inc
465 44040 : k2 = kspan + k2
466 44040 : IF (k2 < ks) GO TO 380
467 46432 : 390 k2 = k2 - np(j)
468 46432 : j = j + 1
469 46432 : k2 = np(j + 1) + k2
470 46432 : IF (k2 > np(j)) GO TO 390
471 110264 : j = 1
472 110264 : 400 IF (kk < k2) GO TO 380
473 83296 : kk = kk + inc
474 83296 : k2 = kspan + k2
475 83296 : IF (k2 < ks) GO TO 400
476 14680 : IF (kk < ks) GO TO 390
477 : jc = k3
478 1060959 : GO TO 450
479 : ! permutation for multivariate transform
480 1060959 : 410 k = kk + jc
481 3134295 : 420 ak = a(kk)
482 3134295 : a(kk) = a(k2)
483 3134295 : a(k2) = ak
484 3134295 : bk = b(kk)
485 3134295 : b(kk) = b(k2)
486 3134295 : b(k2) = bk
487 3134295 : kk = kk + inc
488 3134295 : k2 = k2 + inc
489 3134295 : IF (kk < k) GO TO 420
490 1060959 : kk = kk + ks - jc
491 1060959 : k2 = k2 + ks - jc
492 1060959 : IF (kk < nt) GO TO 410
493 1989 : k2 = k2 - nt + kspan
494 1989 : kk = kk - nt + jc
495 1989 : IF (k2 < ks) GO TO 410
496 2280 : 430 k2 = k2 - np(j)
497 2280 : j = j + 1
498 2280 : k2 = np(j + 1) + k2
499 2280 : IF (k2 > np(j)) GO TO 430
500 5098 : j = 1
501 5098 : 440 IF (kk < k2) GO TO 410
502 3907 : kk = kk + jc
503 3907 : k2 = kspan + k2
504 3907 : IF (k2 < ks) GO TO 440
505 782 : IF (kk < ks) GO TO 430
506 12254 : jc = k3
507 14812 : 450 IF (2*kt + 1 >= m) GO TO 667
508 14695 : kspnn = np(kt + 1)
509 : ! permutation for square-free factors of n
510 14695 : j = m - kt
511 14695 : nfac(j + 1) = 1
512 29824 : 460 nfac(j) = nfac(j)*nfac(j + 1)
513 29824 : j = j - 1
514 29824 : IF (j /= kt) GO TO 460
515 14695 : kt = kt + 1
516 14695 : nn = nfac(kt) - 1
517 14695 : IF (nn > maxp) GO TO 590
518 : jj = 0
519 : j = 0
520 146848 : GO TO 490
521 146848 : 470 jj = jj - k2
522 146848 : k2 = kk
523 146848 : k = k + 1
524 146848 : kk = nfac(k)
525 574019 : 480 jj = kk + jj
526 574019 : IF (jj >= k2) GO TO 470
527 427171 : np(j) = jj
528 441866 : 490 k2 = nfac(kt)
529 441866 : k = kt + 1
530 441866 : kk = nfac(k)
531 441866 : j = j + 1
532 441866 : IF (j <= nn) GO TO 480
533 : ! determine the permutation cycles of length greater than 1
534 : j = 0
535 312997 : GO TO 510
536 312997 : 500 k = kk
537 312997 : kk = np(k)
538 312997 : np(k) = -kk
539 312997 : IF (kk /= j) GO TO 500
540 325409 : k3 = kk
541 427171 : 510 j = j + 1
542 427171 : kk = np(j)
543 427171 : IF (kk < 0) GO TO 510
544 114174 : IF (kk /= j) GO TO 500
545 27107 : np(j) = -j
546 27107 : IF (j /= nn) GO TO 510
547 14695 : maxf = inc*maxf
548 : ! reorder a and b, following the permutation cycles
549 526142 : GO TO 580
550 511447 : 520 j = j - 1
551 511447 : IF (np(j) < 0) GO TO 520
552 629 : jj = jc
553 242638 : 530 kspan = jj
554 : IF (jj > maxf) kspan = maxf
555 242638 : jj = jj - kspan
556 242638 : k = np(j)
557 242638 : kk = jc*k + ii + jj
558 242638 : k1 = kk + kspan
559 242638 : k2 = 0
560 651747 : 540 k2 = k2 + 1
561 651747 : at(k2) = a(k1)
562 651747 : bt(k2) = b(k1)
563 651747 : k1 = k1 - inc
564 651747 : IF (k1 /= kk) GO TO 540
565 1946642 : 550 k1 = kk + kspan
566 1946642 : k2 = k1 - jc*(k + np(k))
567 1946642 : k = -np(k)
568 5404401 : 560 a(k1) = a(k2)
569 5404401 : b(k1) = b(k2)
570 5404401 : k1 = k1 - inc
571 5404401 : k2 = k2 - inc
572 5404401 : IF (k1 /= kk) GO TO 560
573 1946642 : kk = k2
574 1946642 : IF (k /= j) GO TO 550
575 242638 : k1 = kk + kspan
576 242638 : k2 = 0
577 651747 : 570 k2 = k2 + 1
578 651747 : a(k1) = at(k2)
579 651747 : b(k1) = bt(k2)
580 651747 : k1 = k1 - inc
581 651747 : IF (k1 /= kk) GO TO 570
582 242638 : IF (jj /= 0) GO TO 530
583 242009 : IF (j /= 1) GO TO 520
584 104604 : 580 j = k3 + 1
585 104604 : nt = nt - kspnn
586 104604 : ii = nt - inc + 1
587 104604 : IF (nt >= 0) GO TO 520
588 0 : GOTO 667
589 : ! error finish, insufficient array storage
590 : 590 CONTINUE
591 : ! isn = 0
592 0 : WRITE (oUnit, FMT=8000)
593 117 : CALL juDFT_error('array bounds exceeded', calledby='cfft')
594 : 8000 FORMAT('array bounds exceeded within subroutine cft')
595 : 667 CONTINUE
596 14812 : DEALLOCATE (at, bt, ck, sk, nfac, np)
597 : END SUBROUTINE
598 :
599 : #else
600 : SUBROUTINE cfft(a, b, ntot, n, nspan, isn)
601 :
602 : !-------------------------------------------------------------*
603 : ! driver routine for ccfft subroutine instead of cfft on cray *
604 : ! and dcft, dcft2 and dcft3 essl routines on IBM *
605 : !-------------------------------------------------------------*
606 :
607 : IMPLICIT NONE
608 :
609 : ! ... input variables
610 : INTEGER :: ntot, n, nspan, isn
611 : REAL :: a(ntot), b(ntot)
612 :
613 : ! ... local variables
614 : INTEGER :: i, ld1, ld2, n1, n2, n3, dimfft, idim, s(4)
615 :
616 : LOGICAL :: calc
617 :
618 : REAL, DIMENSION(:), ALLOCATABLE :: table, aux
619 : REAL, DIMENSION(:), ALLOCATABLE :: work, aux1, aux2
620 : COMPLEX, DIMENSION(:), ALLOCATABLE :: x
621 :
622 : INTEGER :: naux, naux1, naux2, lam, la1, la2
623 : REAL, PARAMETER :: scale = 1.0
624 : ! ... save variables
625 : SAVE n1, n2, n3
626 :
627 : ! ... data statements
628 : DATA s/1, 1, 1, 1/
629 :
630 : ! ... now, what do we have to do ?
631 :
632 : IF ((ntot == n) .AND. (n == nspan)) THEN
633 : ! ... 1D-FFT
634 : dimfft = 1
635 : n1 = n
636 : n2 = 2
637 : n3 = 2
638 : calc = .TRUE.
639 : ELSE
640 : IF (n == nspan) THEN
641 : ! ... 2D or 3D first step
642 : n1 = n
643 : n2 = 0
644 : calc = .FALSE.
645 : ELSE
646 : IF (ntot == nspan) THEN
647 : ! ... 2D second step or 3D third step
648 : IF (n2 == 0) THEN
649 : dimfft = 2
650 : n2 = n
651 : n3 = 1
652 : calc = .TRUE.
653 : ELSE
654 : dimfft = 3
655 : n3 = n
656 : calc = .TRUE.
657 : ENDIF
658 : ELSE
659 : ! ... 3D second step.
660 : n2 = n
661 : calc = .FALSE.
662 : ENDIF
663 : ENDIF
664 : ENDIF
665 :
666 : IF (calc) THEN
667 :
668 : ! ... build x from a and b
669 :
670 : ALLOCATE (x(ntot))
671 : x = (0.0, 0.0)
672 : DO i = 1, ntot
673 : x(i) = cmplx(a(i), b(i))
674 : ENDDO
675 : ! ... do the FFT
676 :
677 : #ifdef CPP_AIX
678 :
679 : ld1 = n1
680 : ld2 = n1*n2
681 : IF (dimfft == 1) THEN
682 : naux1 = 20000
683 : IF (n1 > 2048) naux1 = naux1 + CEILING(2.28*n1)
684 : ALLOCATE (aux1(naux1), aux2(naux1))
685 : ELSEIF (dimfft == 2) THEN
686 : naux1 = 40000 + CEILING(2.28*(n1 + n2))
687 : IF (max(n1, n2) <= 2048) naux1 = 40000
688 : naux2 = 20000 + CEILING((2*max(n1, n2) + 256)* &
689 : (2.28 + min(64, n1, n2)))
690 : IF (max(n1, n2) < 256) naux2 = 20000
691 : ALLOCATE (aux1(naux1), aux2(naux2))
692 : ELSE IF (dimfft == 3) THEN
693 : IF (max(n2, n3) < 252) THEN
694 : IF (n1 <= 2048) THEN
695 : naux = 60000
696 : ELSE
697 : naux = 60000 + CEILING(4.56*n1)
698 : ENDIF
699 : ELSE
700 : la1 = CEILING((2*n2 + 256)*(min(64, n1) + 4.56))
701 : la2 = CEILING((2*n3 + 256)*(min(64, n1*n2) + 4.56))
702 : lam = max(la2, la1)
703 : IF ((n2 >= 252) .AND. (n3 < 252)) lam = la1
704 : IF ((n2 < 252) .AND. (n3 >= 252)) lam = la2
705 : IF (n1 <= 2048) THEN
706 : naux = 60000 + lam
707 : ELSE
708 : naux = 60000 + CEILING(4.56*n1) + lam
709 : ENDIF
710 : ENDIF
711 : ALLOCATE (aux(naux))
712 : ENDIF
713 : #else
714 :
715 : ld1 = n1
716 : ld2 = n2
717 : s(1) = dimfft
718 :
719 : #ifndef CPP_MPI
720 : ! t,j90:
721 : idim = 1024*n
722 : ALLOCATE (table(16*n + 100), work(idim))
723 : #else
724 : ! t3e:
725 : idim = 2*n1*max(n2, 1)*max(n3, 1)
726 : ALLOCATE (table(12*(n1 + n2 + n3)), work(idim))
727 : #endif
728 : #endif
729 : IF (dimfft == 1) THEN
730 : #ifdef CPP_AIX
731 : CALL dcft(-1, x, 1, 1, x, 1, 1, n, 1, -isn, 1.0, &
732 : aux1, naux1, aux2, naux1)
733 : CALL dcft(0, x, 1, 1, x, 1, 1, n, 1, -isn, 1.0, &
734 : aux1, naux1, aux2, naux1)
735 : #else
736 : CALL ccfft(0, n, 1.0, x, x, table, work, s)
737 : CALL ccfft(isn, n, 1.0, x, x, table, work, s)
738 : #endif
739 : ENDIF
740 : IF (dimfft == 2) THEN
741 : #ifdef CPP_AIX
742 : CALL dcft2(-1, x, 1, n1, x, 1, n1, n1, n2, -isn, 1.0, &
743 : aux1, naux1, aux2, naux2)
744 : CALL dcft2(0, x, 1, n1, x, 1, n1, n1, n2, -isn, 1.0, &
745 : aux1, naux1, aux2, naux2)
746 : #else
747 : CALL ccfft2d(0, n1, n2, 1.0, x, ld1, x, ld1, table, work, s)
748 : CALL ccfft2d(isn, n1, n2, 1.0, x, ld1, x, ld1, table, work, s)
749 : #endif
750 : ENDIF
751 : IF (dimfft == 3) THEN
752 : #ifdef CPP_AIX
753 : CALL dcft3(x, ld1, ld2, x, ld1, ld2, n1, n2, n3, &
754 : -isn, scale, aux, naux)
755 : #else
756 : CALL ccfft3d(0, n1, n2, n3, 1.0, x, ld1, ld2, x, ld1, ld2, &
757 : table, work, s)
758 : CALL ccfft3d(isn, n1, n2, n3, 1.0, x, ld1, ld2, x, ld1, ld2, &
759 : table, work, s)
760 : #endif
761 : ENDIF
762 :
763 : #ifdef CPP_AIX
764 : IF (dimfft == 3) THEN
765 : DEALLOCATE (aux)
766 : ELSE
767 : DEALLOCATE (aux1, aux2)
768 : ENDIF
769 : #else
770 : DEALLOCATE (table, work)
771 : #endif
772 :
773 : ! ... backup a and b
774 :
775 : DO i = 1, ntot
776 : a(i) = real(x(i))
777 : b(i) = aimag(x(i))
778 : ENDDO
779 :
780 : DEALLOCATE (x)
781 : ENDIF
782 :
783 : END subroutine
784 : #endif
785 : END module
|
