1 | !> @file singleton_mod.f90 |
---|
2 | !------------------------------------------------------------------------------! |
---|
3 | ! Current revisions: |
---|
4 | ! ----------------- |
---|
5 | ! |
---|
6 | ! |
---|
7 | ! Former revisions: |
---|
8 | ! ----------------- |
---|
9 | ! $Id: singleton_mod.f90 4180 2019-08-21 14:37:54Z scharf $ |
---|
10 | ! statement added to prevent compiler warning about unused variables |
---|
11 | ! |
---|
12 | ! |
---|
13 | ! Description: |
---|
14 | ! ------------ |
---|
15 | !> Multivariate Fast Fourier Transform |
---|
16 | !> |
---|
17 | !> Fortran 90 Implementation of Singleton's mixed-radix algorithm, |
---|
18 | !> RC Singleton, Stanford Research Institute, Sept. 1968. |
---|
19 | !> |
---|
20 | !> Adapted from fftn.c, translated from Fortran 66 to C by Mark Olesen and |
---|
21 | !> John Beale. |
---|
22 | !> |
---|
23 | !> Fourier transforms can be computed either in place, using assumed size |
---|
24 | !> arguments, or by generic function, using assumed shape arguments. |
---|
25 | !> |
---|
26 | !> |
---|
27 | !> Public: |
---|
28 | !> |
---|
29 | !> fftkind kind parameter of complex arguments |
---|
30 | !> and function results. |
---|
31 | !> |
---|
32 | !> fft(array, dim, inv, stat) generic transform function |
---|
33 | !> COMPLEX(fftkind), DIMENSION(:,...,:), INTENT(IN) :: array |
---|
34 | !> INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
35 | !> LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
36 | !> INTEGER, INTENT(OUT), OPTIONAL:: stat |
---|
37 | !> |
---|
38 | !> fftn(array, shape, dim, inv, stat) in place transform subroutine |
---|
39 | !> COMPLEX(fftkind), DIMENSION(*), INTENT(INOUT) :: array |
---|
40 | !> INTEGER, DIMENSION(:), INTENT(IN) :: shape |
---|
41 | !> INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
42 | !> LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
43 | !> INTEGER, INTENT(OUT), OPTIONAL:: stat |
---|
44 | !> |
---|
45 | !> |
---|
46 | !> Formal Parameters: |
---|
47 | !> |
---|
48 | !> array The complex array to be transformed. array can be of arbitrary |
---|
49 | !> rank (i.e. up to seven). |
---|
50 | !> |
---|
51 | !> shape With subroutine fftn, the shape of the array to be transformed |
---|
52 | !> has to be passed separately, since fftradix - the internal trans- |
---|
53 | !> formation routine - will treat array always as one dimensional. |
---|
54 | !> The product of elements in shape must be the number of |
---|
55 | !> elements in array. |
---|
56 | !> Although passing array with assumed shape would have been nicer, |
---|
57 | !> I prefered assumed size in order to prevent the compiler from |
---|
58 | !> using a copy-in-copy-out mechanism. That would generally be |
---|
59 | !> necessary with fftn passing array to fftradix and with fftn |
---|
60 | !> being prepared for accepting non consecutive array sections. |
---|
61 | !> Using assumed size, it's up to the user to pass an array argu- |
---|
62 | !> ment, that can be addressed as continous one dimensional array |
---|
63 | !> without copying. Otherwise, transformation will not really be |
---|
64 | !> performed in place. |
---|
65 | !> On the other hand, since the rank of array and the size of |
---|
66 | !> shape needn't match, fftn is appropriate for handling more than |
---|
67 | !> seven dimensions. |
---|
68 | !> As far as function fft is concerned all this doesn't matter, |
---|
69 | !> because the argument will be copied anyway. Thus no extra |
---|
70 | !> shape argument is needed for fft. |
---|
71 | !> |
---|
72 | !> Optional Parameters: |
---|
73 | !> |
---|
74 | !> dim One dimensional integer array, containing the dimensions to be |
---|
75 | !> transformed. Default is (/1,...,N/) with N being the rank of |
---|
76 | !> array, i.e. complete transform. dim can restrict transformation |
---|
77 | !> to a subset of available dimensions. Its size must not exceed the |
---|
78 | !> rank of array or the size of shape respectivly. |
---|
79 | !> |
---|
80 | !> inv If .true., inverse transformation will be performed. Default is |
---|
81 | !> .false., i.e. forward transformation. |
---|
82 | !> |
---|
83 | !> stat If present, a system dependent nonzero status value will be |
---|
84 | !> returned in stat, if allocation of temporary storage failed. |
---|
85 | !> |
---|
86 | !> |
---|
87 | !> Scaling: |
---|
88 | !> |
---|
89 | !> Transformation results will always be scaled by the square root of the |
---|
90 | !> product of sizes of each dimension in dim. (See examples below) |
---|
91 | !> |
---|
92 | !> |
---|
93 | !> Examples: |
---|
94 | !> |
---|
95 | !> Let A be a L*M*N three dimensional complex array. Then |
---|
96 | !> |
---|
97 | !> result = fft(A) |
---|
98 | !> |
---|
99 | !> will produce a three dimensional transform, scaled by sqrt(L*M*N), while |
---|
100 | !> |
---|
101 | !> call fftn(A, SHAPE(A)) |
---|
102 | !> |
---|
103 | !> will do the same in place. |
---|
104 | !> |
---|
105 | !> result = fft(A, dim=(/1,3/)) |
---|
106 | !> |
---|
107 | !> will transform with respect to the first and the third dimension, scaled |
---|
108 | !> by sqrt(L*N). |
---|
109 | !> |
---|
110 | !> result = fft(fft(A), inv=.true.) |
---|
111 | !> |
---|
112 | !> should (approximately) reproduce A. |
---|
113 | !> With B having the same shape as A |
---|
114 | !> |
---|
115 | !> result = fft(fft(A) * CONJG(fft(B)), inv=.true.) |
---|
116 | !> |
---|
117 | !> will correlate A and B. |
---|
118 | !> |
---|
119 | !> |
---|
120 | !> Remarks: |
---|
121 | !> |
---|
122 | !> Following changes have been introduced with respect to fftn.c: |
---|
123 | !> - complex arguments and results are of type complex, rather than |
---|
124 | !> real an imaginary part separately. |
---|
125 | !> - increment parameter (magnitude of isign) has been dropped, |
---|
126 | !> inc is always one, direction of transform is given by inv. |
---|
127 | !> - maxf and maxp have been dropped. The amount of temporary storage |
---|
128 | !> needed is determined by the fftradix routine. Both fftn and fft |
---|
129 | !> can handle any size of array. (Maybe they take a lot of time and |
---|
130 | !> memory, but they will do it) |
---|
131 | !> |
---|
132 | !> Redesigning fftradix in a way, that it handles assumed shape arrays |
---|
133 | !> would have been desirable. However, I found it rather hard to do this |
---|
134 | !> in an efficient way. Problems were: |
---|
135 | !> - to prevent stride multiplications when indexing arrays. At least our |
---|
136 | !> compiler was not clever enough to discover that in fact additions |
---|
137 | !> would do the job as well. On the other hand, I haven't been clever |
---|
138 | !> enough to find an implementation using array operations. |
---|
139 | !> - fftradix is rather large and different versions would be necessaray |
---|
140 | !> for each possible rank of array. |
---|
141 | !> Consequently, in place transformation still needs the argument stored |
---|
142 | !> in a consecutive bunch of memory and can't be performed on array |
---|
143 | !> sections like A(100:199:-3, 50:1020). Calling fftn with such sections |
---|
144 | !> will most probably imply copy-in-copy-out. However, the function fft |
---|
145 | !> works with everything it gets and should be convenient to use. |
---|
146 | !> |
---|
147 | !> Michael Steffens, 09.12.96, <Michael.Steffens@mbox.muk.uni-hannover.de> |
---|
148 | !> Restructured fftradix for better optimization. M. Steffens, 4 June 1997 |
---|
149 | !------------------------------------------------------------------------------! |
---|
150 | MODULE singleton |
---|
151 | |
---|
152 | |
---|
153 | USE kinds |
---|
154 | |
---|
155 | IMPLICIT NONE |
---|
156 | |
---|
157 | PRIVATE |
---|
158 | PUBLIC:: fft, fftn |
---|
159 | |
---|
160 | REAL(wp), PARAMETER:: sin60 = 0.86602540378443865_wp |
---|
161 | REAL(wp), PARAMETER:: cos72 = 0.30901699437494742_wp |
---|
162 | REAL(wp), PARAMETER:: sin72 = 0.95105651629515357_wp |
---|
163 | REAL(wp), PARAMETER:: pi = 3.14159265358979323_wp |
---|
164 | |
---|
165 | INTERFACE fft |
---|
166 | MODULE PROCEDURE fft1d |
---|
167 | MODULE PROCEDURE fft2d |
---|
168 | MODULE PROCEDURE fft3d |
---|
169 | MODULE PROCEDURE fft4d |
---|
170 | MODULE PROCEDURE fft5d |
---|
171 | MODULE PROCEDURE fft6d |
---|
172 | MODULE PROCEDURE fft7d |
---|
173 | END INTERFACE |
---|
174 | |
---|
175 | |
---|
176 | CONTAINS |
---|
177 | |
---|
178 | |
---|
179 | !------------------------------------------------------------------------------! |
---|
180 | ! Description: |
---|
181 | ! ------------ |
---|
182 | !> @todo Missing function description. |
---|
183 | !------------------------------------------------------------------------------! |
---|
184 | FUNCTION fft1d(array, dim, inv, stat) RESULT(ft) |
---|
185 | ! |
---|
186 | !-- Formal parameters |
---|
187 | COMPLEX(wp), DIMENSION(:), INTENT(IN) :: array |
---|
188 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
189 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
190 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
191 | ! |
---|
192 | !-- Function result |
---|
193 | COMPLEX(wp), DIMENSION(SIZE(array, 1)):: ft |
---|
194 | |
---|
195 | INTEGER(iwp):: idum |
---|
196 | INTEGER(iwp):: ishape(1) |
---|
197 | |
---|
198 | ! |
---|
199 | !-- Intrinsics used |
---|
200 | INTRINSIC SIZE, SHAPE |
---|
201 | |
---|
202 | ft = array |
---|
203 | ishape = SHAPE( array ) |
---|
204 | CALL fftn(ft, ishape, inv = inv, stat = stat) |
---|
205 | ! |
---|
206 | !-- Next statement to prevent compiler warning about unused variable |
---|
207 | IF ( PRESENT( dim ) ) idum = 1 |
---|
208 | |
---|
209 | END FUNCTION fft1d |
---|
210 | |
---|
211 | |
---|
212 | !------------------------------------------------------------------------------! |
---|
213 | ! Description: |
---|
214 | ! ------------ |
---|
215 | !> @todo Missing function description. |
---|
216 | !------------------------------------------------------------------------------! |
---|
217 | FUNCTION fft2d(array, dim, inv, stat) RESULT(ft) |
---|
218 | ! |
---|
219 | !-- Formal parameters |
---|
220 | COMPLEX(wp), DIMENSION(:,:), INTENT(IN) :: array |
---|
221 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
222 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
223 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
224 | ! |
---|
225 | !-- Function result |
---|
226 | COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2)):: ft |
---|
227 | |
---|
228 | INTEGER(iwp) :: ishape(2) |
---|
229 | ! |
---|
230 | !-- Intrinsics used |
---|
231 | INTRINSIC SIZE, SHAPE |
---|
232 | |
---|
233 | ft = array |
---|
234 | ishape = SHAPE( array ) |
---|
235 | CALL fftn(ft, ishape, dim, inv, stat) |
---|
236 | |
---|
237 | END FUNCTION fft2d |
---|
238 | |
---|
239 | |
---|
240 | !------------------------------------------------------------------------------! |
---|
241 | ! Description: |
---|
242 | ! ------------ |
---|
243 | !> @todo Missing function description. |
---|
244 | !------------------------------------------------------------------------------! |
---|
245 | FUNCTION fft3d(array, dim, inv, stat) RESULT(ft) |
---|
246 | ! |
---|
247 | !-- Formal parameters |
---|
248 | COMPLEX(wp), DIMENSION(:,:,:), INTENT(IN) :: array |
---|
249 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
250 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
251 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
252 | ! |
---|
253 | !-- Function result |
---|
254 | COMPLEX(wp), & |
---|
255 | DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3)):: ft |
---|
256 | |
---|
257 | INTEGER(iwp) :: ishape(3) |
---|
258 | |
---|
259 | ! |
---|
260 | !-- Intrinsics used |
---|
261 | INTRINSIC SIZE, SHAPE |
---|
262 | |
---|
263 | ft = array |
---|
264 | ishape = SHAPE( array) |
---|
265 | CALL fftn(ft, ishape, dim, inv, stat) |
---|
266 | |
---|
267 | END FUNCTION fft3d |
---|
268 | |
---|
269 | |
---|
270 | !------------------------------------------------------------------------------! |
---|
271 | ! Description: |
---|
272 | ! ------------ |
---|
273 | !> @todo Missing function description. |
---|
274 | !------------------------------------------------------------------------------! |
---|
275 | FUNCTION fft4d(array, dim, inv, stat) RESULT(ft) |
---|
276 | ! |
---|
277 | !-- Formal parameters |
---|
278 | COMPLEX(wp), DIMENSION(:,:,:,:), INTENT(IN) :: array |
---|
279 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
280 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
281 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
282 | ! |
---|
283 | !-- Function result |
---|
284 | COMPLEX(wp), DIMENSION( & |
---|
285 | SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4)):: ft |
---|
286 | |
---|
287 | INTEGER(iwp) :: ishape(4) |
---|
288 | ! |
---|
289 | !-- Intrinsics used |
---|
290 | INTRINSIC SIZE, SHAPE |
---|
291 | |
---|
292 | ft = array |
---|
293 | ishape = SHAPE( array ) |
---|
294 | CALL fftn(ft, ishape, dim, inv, stat) |
---|
295 | |
---|
296 | END FUNCTION fft4d |
---|
297 | |
---|
298 | |
---|
299 | !------------------------------------------------------------------------------! |
---|
300 | ! Description: |
---|
301 | ! ------------ |
---|
302 | !> @todo Missing function description. |
---|
303 | !------------------------------------------------------------------------------! |
---|
304 | FUNCTION fft5d(array, dim, inv, stat) RESULT(ft) |
---|
305 | ! |
---|
306 | !-- Formal parameters |
---|
307 | COMPLEX(wp), DIMENSION(:,:,:,:,:), INTENT(IN) :: array |
---|
308 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
309 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
310 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
311 | ! |
---|
312 | !-- Function result |
---|
313 | COMPLEX(wp), DIMENSION( & |
---|
314 | SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), & |
---|
315 | SIZE(array, 5)):: ft |
---|
316 | |
---|
317 | INTEGER(iwp) :: ishape(5) |
---|
318 | |
---|
319 | ! |
---|
320 | !-- Intrinsics used |
---|
321 | INTRINSIC SIZE, SHAPE |
---|
322 | |
---|
323 | ft = array |
---|
324 | ishape = SHAPE( array ) |
---|
325 | CALL fftn(ft, ishape, dim, inv, stat) |
---|
326 | |
---|
327 | END FUNCTION fft5d |
---|
328 | |
---|
329 | |
---|
330 | !------------------------------------------------------------------------------! |
---|
331 | ! Description: |
---|
332 | ! ------------ |
---|
333 | !> @todo Missing function description. |
---|
334 | !------------------------------------------------------------------------------! |
---|
335 | FUNCTION fft6d(array, dim, inv, stat) RESULT(ft) |
---|
336 | ! |
---|
337 | !-- Formal parameters |
---|
338 | COMPLEX(wp), DIMENSION(:,:,:,:,:,:), INTENT(IN) :: array |
---|
339 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
340 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
341 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
342 | ! |
---|
343 | !-- Function result |
---|
344 | COMPLEX(wp), DIMENSION( & |
---|
345 | SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), & |
---|
346 | SIZE(array, 5), SIZE(array, 6)):: ft |
---|
347 | |
---|
348 | INTEGER(iwp) :: ishape(6) |
---|
349 | |
---|
350 | ! |
---|
351 | !-- Intrinsics used |
---|
352 | INTRINSIC SIZE, SHAPE |
---|
353 | |
---|
354 | ft = array |
---|
355 | ishape = SHAPE( array ) |
---|
356 | CALL fftn(ft, ishape, dim, inv, stat) |
---|
357 | |
---|
358 | END FUNCTION fft6d |
---|
359 | |
---|
360 | |
---|
361 | !------------------------------------------------------------------------------! |
---|
362 | ! Description: |
---|
363 | ! ------------ |
---|
364 | !> @todo Missing function description. |
---|
365 | !------------------------------------------------------------------------------! |
---|
366 | FUNCTION fft7d(array, dim, inv, stat) RESULT(ft) |
---|
367 | ! |
---|
368 | !-- Formal parameters |
---|
369 | COMPLEX(wp), DIMENSION(:,:,:,:,:,:,:), INTENT(IN) :: array |
---|
370 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
371 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
372 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
373 | ! |
---|
374 | !-- Function result |
---|
375 | COMPLEX(wp), DIMENSION( & |
---|
376 | SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), & |
---|
377 | SIZE(array, 5), SIZE(array, 6), SIZE(array, 7)):: ft |
---|
378 | |
---|
379 | INTEGER(iwp) :: ishape(7) |
---|
380 | |
---|
381 | ! |
---|
382 | !-- Intrinsics used |
---|
383 | INTRINSIC SIZE, SHAPE |
---|
384 | |
---|
385 | ft = array |
---|
386 | ishape = SHAPE( array ) |
---|
387 | CALL fftn(ft, ishape, dim, inv, stat) |
---|
388 | |
---|
389 | END FUNCTION fft7d |
---|
390 | |
---|
391 | |
---|
392 | !------------------------------------------------------------------------------! |
---|
393 | ! Description: |
---|
394 | ! ------------ |
---|
395 | !> @todo Missing subroutine description. |
---|
396 | !------------------------------------------------------------------------------! |
---|
397 | SUBROUTINE fftn(array, shape, dim, inv, stat) |
---|
398 | ! |
---|
399 | !-- Formal parameters |
---|
400 | COMPLEX(wp), DIMENSION(*), INTENT(INOUT) :: array |
---|
401 | INTEGER(iwp), DIMENSION(:), INTENT(IN) :: shape |
---|
402 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
403 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
404 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
405 | ! |
---|
406 | !-- Local arrays |
---|
407 | INTEGER(iwp), DIMENSION(SIZE(shape)):: d |
---|
408 | ! |
---|
409 | !-- Local scalars |
---|
410 | LOGICAL :: inverse |
---|
411 | INTEGER(iwp) :: i, ndim, ntotal |
---|
412 | REAL(wp):: scale |
---|
413 | ! |
---|
414 | !-- Intrinsics used |
---|
415 | INTRINSIC PRESENT, MIN, PRODUCT, SIZE, SQRT |
---|
416 | |
---|
417 | ! |
---|
418 | !-- Optional parameter settings |
---|
419 | IF (PRESENT(inv)) THEN |
---|
420 | inverse = inv |
---|
421 | ELSE |
---|
422 | inverse = .FALSE. |
---|
423 | END IF |
---|
424 | IF (PRESENT(dim)) THEN |
---|
425 | ndim = MIN(SIZE(dim), SIZE(d)) |
---|
426 | d(1:ndim) = DIM(1:ndim) |
---|
427 | ELSE |
---|
428 | ndim = SIZE(d) |
---|
429 | d = (/(i, i = 1, SIZE(d))/) |
---|
430 | END IF |
---|
431 | |
---|
432 | ntotal = PRODUCT(shape) |
---|
433 | scale = SQRT(1.0_wp / PRODUCT(shape(d(1:ndim)))) |
---|
434 | DO i = 1, ntotal |
---|
435 | array(i) = CMPLX(REAL(array(i)) * scale, AIMAG(array(i)) * scale, & |
---|
436 | KIND=wp) |
---|
437 | END DO |
---|
438 | |
---|
439 | DO i = 1, ndim |
---|
440 | CALL fftradix(array, ntotal, shape(d(i)), PRODUCT(shape(1:d(i))), & |
---|
441 | inverse, stat) |
---|
442 | IF (PRESENT(stat)) THEN |
---|
443 | IF (stat /=0) RETURN |
---|
444 | END IF |
---|
445 | END DO |
---|
446 | |
---|
447 | END SUBROUTINE fftn |
---|
448 | |
---|
449 | |
---|
450 | !------------------------------------------------------------------------------! |
---|
451 | ! Description: |
---|
452 | ! ------------ |
---|
453 | !> @todo Missing subroutine description. |
---|
454 | !------------------------------------------------------------------------------! |
---|
455 | SUBROUTINE fftradix(array, ntotal, npass, nspan, inv, stat) |
---|
456 | ! |
---|
457 | !-- Formal parameters |
---|
458 | COMPLEX(wp), DIMENSION(*), INTENT(INOUT) :: array |
---|
459 | INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan |
---|
460 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
461 | LOGICAL, INTENT(IN) :: inv |
---|
462 | ! |
---|
463 | !-- Local arrays |
---|
464 | COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: ctmp |
---|
465 | INTEGER(iwp), DIMENSION(BIT_SIZE(0)) :: factor |
---|
466 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: perm |
---|
467 | REAL(wp), DIMENSION(:), ALLOCATABLE :: sine, cosine |
---|
468 | ! |
---|
469 | !-- Local scalars |
---|
470 | INTEGER(iwp) :: maxfactor, nfactor, nsquare, nperm |
---|
471 | ! |
---|
472 | !-- Intrinsics used |
---|
473 | INTRINSIC MAXVAL, MOD, PRESENT, ISHFT, BIT_SIZE, SIN, COS, & |
---|
474 | CMPLX, REAL, AIMAG |
---|
475 | |
---|
476 | IF (npass <= 1) RETURN |
---|
477 | |
---|
478 | CALL factorize(npass, factor, nfactor, nsquare) |
---|
479 | |
---|
480 | maxfactor = MAXVAL(factor(:nfactor)) |
---|
481 | IF (nfactor - ISHFT(nsquare, 1) > 0) THEN |
---|
482 | nperm = MAX(nfactor + 1, PRODUCT(factor(nsquare+1: nfactor-nsquare)) - 1) |
---|
483 | ELSE |
---|
484 | nperm = nfactor + 1 |
---|
485 | END IF |
---|
486 | |
---|
487 | IF (PRESENT(stat)) THEN |
---|
488 | ALLOCATE(ctmp(maxfactor), sine(maxfactor), cosine(maxfactor), STAT=stat) |
---|
489 | IF (stat /= 0) RETURN |
---|
490 | CALL transform(array, ntotal, npass, nspan, & |
---|
491 | factor, nfactor, ctmp, sine, cosine, inv) |
---|
492 | DEALLOCATE(sine, cosine, STAT=stat) |
---|
493 | IF (stat /= 0) RETURN |
---|
494 | ALLOCATE(perm(nperm), STAT=stat) |
---|
495 | IF (stat /= 0) RETURN |
---|
496 | CALL permute(array, ntotal, npass, nspan, & |
---|
497 | factor, nfactor, nsquare, maxfactor, & |
---|
498 | ctmp, perm) |
---|
499 | DEALLOCATE(perm, ctmp, STAT=stat) |
---|
500 | IF (stat /= 0) RETURN |
---|
501 | ELSE |
---|
502 | ALLOCATE(ctmp(maxfactor), sine(maxfactor), cosine(maxfactor)) |
---|
503 | CALL transform(array, ntotal, npass, nspan, & |
---|
504 | factor, nfactor, ctmp, sine, cosine, inv) |
---|
505 | DEALLOCATE(sine, cosine) |
---|
506 | ALLOCATE(perm(nperm)) |
---|
507 | CALL permute(array, ntotal, npass, nspan, & |
---|
508 | factor, nfactor, nsquare, maxfactor, & |
---|
509 | ctmp, perm) |
---|
510 | DEALLOCATE(perm, ctmp) |
---|
511 | END IF |
---|
512 | |
---|
513 | |
---|
514 | CONTAINS |
---|
515 | |
---|
516 | |
---|
517 | !------------------------------------------------------------------------------! |
---|
518 | ! Description: |
---|
519 | ! ------------ |
---|
520 | !> @todo Missing subroutine description. |
---|
521 | !------------------------------------------------------------------------------! |
---|
522 | SUBROUTINE factorize(npass, factor, nfactor, nsquare) |
---|
523 | ! |
---|
524 | !-- Formal parameters |
---|
525 | INTEGER(iwp), INTENT(IN) :: npass |
---|
526 | INTEGER(iwp), DIMENSION(*), INTENT(OUT):: factor |
---|
527 | INTEGER(iwp), INTENT(OUT):: nfactor, nsquare |
---|
528 | ! |
---|
529 | !-- Local scalars |
---|
530 | INTEGER(iwp):: j, jj, k |
---|
531 | |
---|
532 | nfactor = 0 |
---|
533 | k = npass |
---|
534 | DO WHILE (MOD(k, 16) == 0) |
---|
535 | nfactor = nfactor + 1 |
---|
536 | factor(nfactor) = 4 |
---|
537 | k = k / 16 |
---|
538 | END DO |
---|
539 | j = 3 |
---|
540 | jj = 9 |
---|
541 | DO |
---|
542 | DO WHILE (MOD(k, jj) == 0) |
---|
543 | nfactor = nfactor + 1 |
---|
544 | factor(nfactor) = j |
---|
545 | k = k / jj |
---|
546 | END DO |
---|
547 | j = j + 2 |
---|
548 | jj = j * j |
---|
549 | IF (jj > k) EXIT |
---|
550 | END DO |
---|
551 | IF (k <= 4) THEN |
---|
552 | nsquare = nfactor |
---|
553 | factor(nfactor + 1) = k |
---|
554 | IF (k /= 1) nfactor = nfactor + 1 |
---|
555 | ELSE |
---|
556 | IF (k - ISHFT(k / 4, 2) == 0) THEN |
---|
557 | nfactor = nfactor + 1 |
---|
558 | factor(nfactor) = 2 |
---|
559 | k = k / 4 |
---|
560 | END IF |
---|
561 | nsquare = nfactor |
---|
562 | j = 2 |
---|
563 | DO |
---|
564 | IF (MOD(k, j) == 0) THEN |
---|
565 | nfactor = nfactor + 1 |
---|
566 | factor(nfactor) = j |
---|
567 | k = k / j |
---|
568 | END IF |
---|
569 | j = ISHFT((j + 1) / 2, 1) + 1 |
---|
570 | IF (j > k) EXIT |
---|
571 | END DO |
---|
572 | END IF |
---|
573 | IF (nsquare > 0) THEN |
---|
574 | j = nsquare |
---|
575 | DO |
---|
576 | nfactor = nfactor + 1 |
---|
577 | factor(nfactor) = factor(j) |
---|
578 | j = j - 1 |
---|
579 | IF (j==0) EXIT |
---|
580 | END DO |
---|
581 | END IF |
---|
582 | |
---|
583 | END SUBROUTINE factorize |
---|
584 | |
---|
585 | |
---|
586 | !------------------------------------------------------------------------------! |
---|
587 | ! Description: |
---|
588 | ! ------------ |
---|
589 | !> @todo Missing subroutine description. |
---|
590 | !------------------------------------------------------------------------------! |
---|
591 | SUBROUTINE transform(array, ntotal, npass, nspan, & |
---|
592 | factor, nfactor, ctmp, sine, cosine, inv) !-- compute fourier transform |
---|
593 | ! |
---|
594 | !-- Formal parameters |
---|
595 | COMPLEX(wp), DIMENSION(*), INTENT(IN OUT):: array |
---|
596 | COMPLEX(wp), DIMENSION(*), INTENT(OUT) :: ctmp |
---|
597 | INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan |
---|
598 | INTEGER(iwp), DIMENSION(*), INTENT(IN) :: factor |
---|
599 | INTEGER(iwp), INTENT(IN) :: nfactor |
---|
600 | LOGICAL, INTENT(IN) :: inv |
---|
601 | REAL(wp), DIMENSION(*), INTENT(OUT) :: sine, cosine |
---|
602 | ! |
---|
603 | !-- Local scalars |
---|
604 | INTEGER(iwp):: ii, ispan |
---|
605 | INTEGER(iwp):: j, jc, jf, jj |
---|
606 | INTEGER(iwp):: k, kk, kspan, k1, k2, k3, k4 |
---|
607 | INTEGER(iwp):: nn, nt |
---|
608 | REAL(wp) :: s60, c72, s72, pi2, radf |
---|
609 | REAL(wp) :: c1, s1, c2, s2, c3, s3, cd, sd, ak |
---|
610 | COMPLEX(wp) :: cc, cj, ck, cjp, cjm, ckp, ckm |
---|
611 | |
---|
612 | c72 = cos72 |
---|
613 | IF (inv) THEN |
---|
614 | s72 = sin72 |
---|
615 | s60 = sin60 |
---|
616 | pi2 = pi |
---|
617 | ELSE |
---|
618 | s72 = -sin72 |
---|
619 | s60 = -sin60 |
---|
620 | pi2 = -pi |
---|
621 | END IF |
---|
622 | |
---|
623 | nt = ntotal |
---|
624 | nn = nt - 1 |
---|
625 | kspan = nspan |
---|
626 | jc = nspan / npass |
---|
627 | radf = pi2 * jc |
---|
628 | pi2 = pi2 * 2.0_wp !-- use 2 PI from here on |
---|
629 | |
---|
630 | ii = 0 |
---|
631 | jf = 0 |
---|
632 | DO |
---|
633 | sd = radf / kspan |
---|
634 | cd = SIN(sd) |
---|
635 | cd = 2.0_wp * cd * cd |
---|
636 | sd = SIN(sd + sd) |
---|
637 | kk = 1 |
---|
638 | ii = ii + 1 |
---|
639 | |
---|
640 | SELECT CASE (factor(ii)) |
---|
641 | CASE (2) |
---|
642 | ! |
---|
643 | !-- Transform for factor of 2 (including rotation factor) |
---|
644 | kspan = kspan / 2 |
---|
645 | k1 = kspan + 2 |
---|
646 | DO |
---|
647 | DO |
---|
648 | k2 = kk + kspan |
---|
649 | ck = array(k2) |
---|
650 | array(k2) = array(kk)-ck |
---|
651 | array(kk) = array(kk) + ck |
---|
652 | kk = k2 + kspan |
---|
653 | IF (kk > nn) EXIT |
---|
654 | END DO |
---|
655 | kk = kk - nn |
---|
656 | IF (kk > jc) EXIT |
---|
657 | END DO |
---|
658 | IF (kk > kspan) RETURN |
---|
659 | DO |
---|
660 | c1 = 1.0_wp - cd |
---|
661 | s1 = sd |
---|
662 | DO |
---|
663 | DO |
---|
664 | DO |
---|
665 | k2 = kk + kspan |
---|
666 | ck = array(kk) - array(k2) |
---|
667 | array(kk) = array(kk) + array(k2) |
---|
668 | array(k2) = ck * CMPLX(c1, s1, KIND=wp) |
---|
669 | kk = k2 + kspan |
---|
670 | IF (kk >= nt) EXIT |
---|
671 | END DO |
---|
672 | k2 = kk - nt |
---|
673 | c1 = -c1 |
---|
674 | kk = k1 - k2 |
---|
675 | IF (kk <= k2) EXIT |
---|
676 | END DO |
---|
677 | ak = c1 - (cd * c1 + sd * s1) |
---|
678 | s1 = sd * c1 - cd * s1 + s1 |
---|
679 | c1 = 2.0_wp - (ak * ak + s1 * s1) |
---|
680 | s1 = s1 * c1 |
---|
681 | c1 = c1 * ak |
---|
682 | kk = kk + jc |
---|
683 | IF (kk >= k2) EXIT |
---|
684 | END DO |
---|
685 | k1 = k1 + 1 + 1 |
---|
686 | kk = (k1 - kspan) / 2 + jc |
---|
687 | IF (kk > jc + jc) EXIT |
---|
688 | END DO |
---|
689 | |
---|
690 | CASE (4) !-- transform for factor of 4 |
---|
691 | ispan = kspan |
---|
692 | kspan = kspan / 4 |
---|
693 | |
---|
694 | DO |
---|
695 | c1 = 1.0_wp |
---|
696 | s1 = 0.0_wp |
---|
697 | DO |
---|
698 | DO |
---|
699 | k1 = kk + kspan |
---|
700 | k2 = k1 + kspan |
---|
701 | k3 = k2 + kspan |
---|
702 | ckp = array(kk) + array(k2) |
---|
703 | ckm = array(kk) - array(k2) |
---|
704 | cjp = array(k1) + array(k3) |
---|
705 | cjm = array(k1) - array(k3) |
---|
706 | array(kk) = ckp + cjp |
---|
707 | cjp = ckp - cjp |
---|
708 | IF (inv) THEN |
---|
709 | ckp = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND=wp) |
---|
710 | ckm = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND=wp) |
---|
711 | ELSE |
---|
712 | ckp = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND=wp) |
---|
713 | ckm = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND=wp) |
---|
714 | END IF |
---|
715 | ! |
---|
716 | !-- Avoid useless multiplies |
---|
717 | IF (s1 == 0.0_wp) THEN |
---|
718 | array(k1) = ckp |
---|
719 | array(k2) = cjp |
---|
720 | array(k3) = ckm |
---|
721 | ELSE |
---|
722 | array(k1) = ckp * CMPLX(c1, s1, KIND=wp) |
---|
723 | array(k2) = cjp * CMPLX(c2, s2, KIND=wp) |
---|
724 | array(k3) = ckm * CMPLX(c3, s3, KIND=wp) |
---|
725 | END IF |
---|
726 | kk = k3 + kspan |
---|
727 | IF (kk > nt) EXIT |
---|
728 | END DO |
---|
729 | |
---|
730 | c2 = c1 - (cd * c1 + sd * s1) |
---|
731 | s1 = sd * c1 - cd * s1 + s1 |
---|
732 | c1 = 2.0_wp - (c2 * c2 + s1 * s1) |
---|
733 | s1 = s1 * c1 |
---|
734 | c1 = c1 * c2 |
---|
735 | ! |
---|
736 | !-- Values of c2, c3, s2, s3 that will get used next time |
---|
737 | c2 = c1 * c1 - s1 * s1 |
---|
738 | s2 = 2.0_wp * c1 * s1 |
---|
739 | c3 = c2 * c1 - s2 * s1 |
---|
740 | s3 = c2 * s1 + s2 * c1 |
---|
741 | kk = kk - nt + jc |
---|
742 | IF (kk > kspan) EXIT |
---|
743 | END DO |
---|
744 | kk = kk - kspan + 1 |
---|
745 | IF (kk > jc) EXIT |
---|
746 | END DO |
---|
747 | IF (kspan == jc) RETURN |
---|
748 | |
---|
749 | CASE default |
---|
750 | ! |
---|
751 | !-- Transform for odd factors |
---|
752 | k = factor(ii) |
---|
753 | ispan = kspan |
---|
754 | kspan = kspan / k |
---|
755 | |
---|
756 | SELECT CASE (k) |
---|
757 | CASE (3) !-- transform for factor of 3 (optional code) |
---|
758 | DO |
---|
759 | DO |
---|
760 | k1 = kk + kspan |
---|
761 | k2 = k1 + kspan |
---|
762 | ck = array(kk) |
---|
763 | cj = array(k1) + array(k2) |
---|
764 | array(kk) = ck + cj |
---|
765 | ck = ck - CMPLX( & |
---|
766 | 0.5_wp * REAL (cj), & |
---|
767 | 0.5_wp * AIMAG(cj), & |
---|
768 | KIND=wp) |
---|
769 | cj = CMPLX( & |
---|
770 | (REAL (array(k1)) - REAL (array(k2))) * s60, & |
---|
771 | (AIMAG(array(k1)) - AIMAG(array(k2))) * s60, & |
---|
772 | KIND=wp) |
---|
773 | array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp) |
---|
774 | array(k2) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp) |
---|
775 | kk = k2 + kspan |
---|
776 | IF (kk >= nn) EXIT |
---|
777 | END DO |
---|
778 | kk = kk - nn |
---|
779 | IF (kk > kspan) EXIT |
---|
780 | END DO |
---|
781 | |
---|
782 | CASE (5) !-- transform for factor of 5 (optional code) |
---|
783 | c2 = c72 * c72 - s72 * s72 |
---|
784 | s2 = 2.0_wp * c72 * s72 |
---|
785 | DO |
---|
786 | DO |
---|
787 | k1 = kk + kspan |
---|
788 | k2 = k1 + kspan |
---|
789 | k3 = k2 + kspan |
---|
790 | k4 = k3 + kspan |
---|
791 | ckp = array(k1) + array(k4) |
---|
792 | ckm = array(k1) - array(k4) |
---|
793 | cjp = array(k2) + array(k3) |
---|
794 | cjm = array(k2) - array(k3) |
---|
795 | cc = array(kk) |
---|
796 | array(kk) = cc + ckp + cjp |
---|
797 | ck = CMPLX(REAL(ckp) * c72, AIMAG(ckp) * c72, & |
---|
798 | KIND=wp) + & |
---|
799 | CMPLX(REAL(cjp) * c2, AIMAG(cjp) * c2, & |
---|
800 | KIND=wp) + cc |
---|
801 | cj = CMPLX(REAL(ckm) * s72, AIMAG(ckm) * s72, & |
---|
802 | KIND=wp) + & |
---|
803 | CMPLX(REAL(cjm) * s2, AIMAG(cjm) * s2, & |
---|
804 | KIND=wp) |
---|
805 | array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp) |
---|
806 | array(k4) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp) |
---|
807 | ck = CMPLX(REAL(ckp) * c2, AIMAG(ckp) * c2, & |
---|
808 | KIND=wp) + & |
---|
809 | CMPLX(REAL(cjp) * c72, AIMAG(cjp) * c72, & |
---|
810 | KIND=wp) + cc |
---|
811 | cj = CMPLX(REAL(ckm) * s2, AIMAG(ckm) * s2, & |
---|
812 | KIND=wp) - & |
---|
813 | CMPLX(REAL(cjm) * s72, AIMAG(cjm) * s72, & |
---|
814 | KIND=wp) |
---|
815 | array(k2) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp) |
---|
816 | array(k3) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp) |
---|
817 | kk = k4 + kspan |
---|
818 | IF (kk >= nn) EXIT |
---|
819 | END DO |
---|
820 | kk = kk - nn |
---|
821 | IF (kk > kspan) EXIT |
---|
822 | END DO |
---|
823 | |
---|
824 | CASE default |
---|
825 | IF (k /= jf) THEN |
---|
826 | jf = k |
---|
827 | s1 = pi2 / k |
---|
828 | c1 = COS(s1) |
---|
829 | s1 = SIN(s1) |
---|
830 | cosine (jf) = 1.0_wp |
---|
831 | sine (jf) = 0.0_wp |
---|
832 | j = 1 |
---|
833 | DO |
---|
834 | cosine (j) = cosine (k) * c1 + sine (k) * s1 |
---|
835 | sine (j) = cosine (k) * s1 - sine (k) * c1 |
---|
836 | k = k-1 |
---|
837 | cosine (k) = cosine (j) |
---|
838 | sine (k) = -sine (j) |
---|
839 | j = j + 1 |
---|
840 | IF (j >= k) EXIT |
---|
841 | END DO |
---|
842 | END IF |
---|
843 | DO |
---|
844 | DO |
---|
845 | k1 = kk |
---|
846 | k2 = kk + ispan |
---|
847 | cc = array(kk) |
---|
848 | ck = cc |
---|
849 | j = 1 |
---|
850 | k1 = k1 + kspan |
---|
851 | DO |
---|
852 | k2 = k2 - kspan |
---|
853 | j = j + 1 |
---|
854 | ctmp(j) = array(k1) + array(k2) |
---|
855 | ck = ck + ctmp(j) |
---|
856 | j = j + 1 |
---|
857 | ctmp(j) = array(k1) - array(k2) |
---|
858 | k1 = k1 + kspan |
---|
859 | IF (k1 >= k2) EXIT |
---|
860 | END DO |
---|
861 | array(kk) = ck |
---|
862 | k1 = kk |
---|
863 | k2 = kk + ispan |
---|
864 | j = 1 |
---|
865 | DO |
---|
866 | k1 = k1 + kspan |
---|
867 | k2 = k2 - kspan |
---|
868 | jj = j |
---|
869 | ck = cc |
---|
870 | cj = (0.0_wp, 0.0_wp) |
---|
871 | k = 1 |
---|
872 | DO |
---|
873 | k = k + 1 |
---|
874 | ck = ck + CMPLX( & |
---|
875 | REAL (ctmp(k)) * cosine(jj), & |
---|
876 | AIMAG(ctmp(k)) * cosine(jj), KIND=wp) |
---|
877 | k = k + 1 |
---|
878 | cj = cj + CMPLX( & |
---|
879 | REAL (ctmp(k)) * sine(jj), & |
---|
880 | AIMAG(ctmp(k)) * sine(jj), KIND=wp) |
---|
881 | jj = jj + j |
---|
882 | IF (jj > jf) jj = jj - jf |
---|
883 | IF (k >= jf) EXIT |
---|
884 | END DO |
---|
885 | k = jf - j |
---|
886 | array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), & |
---|
887 | KIND=wp) |
---|
888 | array(k2) = ck + CMPLX(AIMAG(cj), -REAL(cj), & |
---|
889 | KIND=wp) |
---|
890 | j = j + 1 |
---|
891 | IF (j >= k) EXIT |
---|
892 | END DO |
---|
893 | kk = kk + ispan |
---|
894 | IF (kk > nn) EXIT |
---|
895 | END DO |
---|
896 | kk = kk - nn |
---|
897 | IF (kk > kspan) EXIT |
---|
898 | END DO |
---|
899 | |
---|
900 | END SELECT |
---|
901 | ! |
---|
902 | !-- Multiply by rotation factor (except for factors of 2 and 4) |
---|
903 | IF (ii == nfactor) RETURN |
---|
904 | kk = jc + 1 |
---|
905 | DO |
---|
906 | c2 = 1.0_wp - cd |
---|
907 | s1 = sd |
---|
908 | DO |
---|
909 | c1 = c2 |
---|
910 | s2 = s1 |
---|
911 | kk = kk + kspan |
---|
912 | DO |
---|
913 | DO |
---|
914 | array(kk) = CMPLX(c2, s2, KIND=wp) * array(kk) |
---|
915 | kk = kk + ispan |
---|
916 | IF (kk > nt) EXIT |
---|
917 | END DO |
---|
918 | ak = s1 * s2 |
---|
919 | s2 = s1 * c2 + c1 * s2 |
---|
920 | c2 = c1 * c2 - ak |
---|
921 | kk = kk - nt + kspan |
---|
922 | IF (kk > ispan) EXIT |
---|
923 | END DO |
---|
924 | c2 = c1 - (cd * c1 + sd * s1) |
---|
925 | s1 = s1 + sd * c1 - cd * s1 |
---|
926 | c1 = 2.0_wp - (c2 * c2 + s1 * s1) |
---|
927 | s1 = s1 * c1 |
---|
928 | c2 = c2 * c1 |
---|
929 | kk = kk - ispan + jc |
---|
930 | IF (kk > kspan) EXIT |
---|
931 | END DO |
---|
932 | kk = kk - kspan + jc + 1 |
---|
933 | IF (kk > jc + jc) EXIT |
---|
934 | END DO |
---|
935 | |
---|
936 | END SELECT |
---|
937 | END DO |
---|
938 | END SUBROUTINE transform |
---|
939 | |
---|
940 | |
---|
941 | !------------------------------------------------------------------------------! |
---|
942 | ! Description: |
---|
943 | ! ------------ |
---|
944 | !> @todo Missing subroutine description. |
---|
945 | !------------------------------------------------------------------------------! |
---|
946 | SUBROUTINE permute(array, ntotal, npass, nspan, & |
---|
947 | factor, nfactor, nsquare, maxfactor, & |
---|
948 | ctmp, perm) |
---|
949 | ! |
---|
950 | !-- Formal parameters |
---|
951 | COMPLEX(wp), DIMENSION(*), INTENT(IN OUT):: array |
---|
952 | COMPLEX(wp), DIMENSION(*), INTENT(OUT) :: ctmp |
---|
953 | INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan |
---|
954 | INTEGER(iwp), DIMENSION(*), INTENT(IN OUT):: factor |
---|
955 | INTEGER(iwp), INTENT(IN) :: nfactor, nsquare |
---|
956 | INTEGER(iwp), INTENT(IN) :: maxfactor |
---|
957 | INTEGER(iwp), DIMENSION(*), INTENT(OUT) :: perm |
---|
958 | ! |
---|
959 | !-- Local scalars |
---|
960 | COMPLEX(wp) :: ck |
---|
961 | INTEGER(iwp):: ii, ispan |
---|
962 | INTEGER(iwp):: j, jc, jj |
---|
963 | INTEGER(iwp):: k, kk, kspan, kt, k1, k2, k3 |
---|
964 | INTEGER(iwp):: nn, nt |
---|
965 | ! |
---|
966 | !-- Permute the results to normal order---done in two stages |
---|
967 | !-- Permutation for square factors of n |
---|
968 | |
---|
969 | nt = ntotal |
---|
970 | nn = nt - 1 |
---|
971 | kt = nsquare |
---|
972 | kspan = nspan |
---|
973 | jc = nspan / npass |
---|
974 | |
---|
975 | perm (1) = nspan |
---|
976 | IF (kt > 0) THEN |
---|
977 | k = kt + kt + 1 |
---|
978 | IF (nfactor < k) k = k - 1 |
---|
979 | j = 1 |
---|
980 | perm (k + 1) = jc |
---|
981 | DO |
---|
982 | perm (j + 1) = perm (j) / factor(j) |
---|
983 | perm (k) = perm (k + 1) * factor(j) |
---|
984 | j = j + 1 |
---|
985 | k = k - 1 |
---|
986 | IF (j >= k) EXIT |
---|
987 | END DO |
---|
988 | k3 = perm (k + 1) |
---|
989 | kspan = perm (2) |
---|
990 | kk = jc + 1 |
---|
991 | k2 = kspan + 1 |
---|
992 | j = 1 |
---|
993 | |
---|
994 | IF (npass /= ntotal) THEN |
---|
995 | permute_multi: DO |
---|
996 | DO |
---|
997 | DO |
---|
998 | k = kk + jc |
---|
999 | DO |
---|
1000 | ! |
---|
1001 | !-- Swap array(kk) <> array(k2) |
---|
1002 | ck = array(kk) |
---|
1003 | array(kk) = array(k2) |
---|
1004 | array(k2) = ck |
---|
1005 | kk = kk + 1 |
---|
1006 | k2 = k2 + 1 |
---|
1007 | IF (kk >= k) EXIT |
---|
1008 | END DO |
---|
1009 | kk = kk + nspan - jc |
---|
1010 | k2 = k2 + nspan - jc |
---|
1011 | IF (kk >= nt) EXIT |
---|
1012 | END DO |
---|
1013 | kk = kk - nt + jc |
---|
1014 | k2 = k2 - nt + kspan |
---|
1015 | IF (k2 >= nspan) EXIT |
---|
1016 | END DO |
---|
1017 | DO |
---|
1018 | DO |
---|
1019 | k2 = k2 - perm (j) |
---|
1020 | j = j + 1 |
---|
1021 | k2 = perm (j + 1) + k2 |
---|
1022 | IF (k2 <= perm (j)) EXIT |
---|
1023 | END DO |
---|
1024 | j = 1 |
---|
1025 | DO |
---|
1026 | IF (kk < k2) CYCLE permute_multi |
---|
1027 | kk = kk + jc |
---|
1028 | k2 = k2 + kspan |
---|
1029 | IF (k2 >= nspan) EXIT |
---|
1030 | END DO |
---|
1031 | IF (kk >= nspan) EXIT |
---|
1032 | END DO |
---|
1033 | EXIT |
---|
1034 | END DO permute_multi |
---|
1035 | ELSE |
---|
1036 | permute_single: DO |
---|
1037 | DO |
---|
1038 | ! |
---|
1039 | !-- Swap array(kk) <> array(k2) |
---|
1040 | ck = array(kk) |
---|
1041 | array(kk) = array(k2) |
---|
1042 | array(k2) = ck |
---|
1043 | kk = kk + 1 |
---|
1044 | k2 = k2 + kspan |
---|
1045 | IF (k2 >= nspan) EXIT |
---|
1046 | END DO |
---|
1047 | DO |
---|
1048 | DO |
---|
1049 | k2 = k2 - perm (j) |
---|
1050 | j = j + 1 |
---|
1051 | k2 = perm (j + 1) + k2 |
---|
1052 | IF (k2 <= perm (j)) EXIT |
---|
1053 | END DO |
---|
1054 | j = 1 |
---|
1055 | DO |
---|
1056 | IF (kk < k2) CYCLE permute_single |
---|
1057 | kk = kk + 1 |
---|
1058 | k2 = k2 + kspan |
---|
1059 | IF (k2 >= nspan) EXIT |
---|
1060 | END DO |
---|
1061 | IF (kk >= nspan) EXIT |
---|
1062 | END DO |
---|
1063 | EXIT |
---|
1064 | END DO permute_single |
---|
1065 | END IF |
---|
1066 | jc = k3 |
---|
1067 | END IF |
---|
1068 | |
---|
1069 | IF (ISHFT(kt, 1) + 1 >= nfactor) RETURN |
---|
1070 | |
---|
1071 | ispan = perm (kt + 1) |
---|
1072 | ! |
---|
1073 | !-- Permutation for square-free factors of n |
---|
1074 | j = nfactor - kt |
---|
1075 | factor(j + 1) = 1 |
---|
1076 | DO |
---|
1077 | factor(j) = factor(j) * factor(j+1) |
---|
1078 | j = j - 1 |
---|
1079 | IF (j == kt) EXIT |
---|
1080 | END DO |
---|
1081 | kt = kt + 1 |
---|
1082 | nn = factor(kt) - 1 |
---|
1083 | j = 0 |
---|
1084 | jj = 0 |
---|
1085 | DO |
---|
1086 | k = kt + 1 |
---|
1087 | k2 = factor(kt) |
---|
1088 | kk = factor(k) |
---|
1089 | j = j + 1 |
---|
1090 | IF (j > nn) EXIT !-- exit infinite loop |
---|
1091 | jj = jj + kk |
---|
1092 | DO WHILE (jj >= k2) |
---|
1093 | jj = jj - k2 |
---|
1094 | k2 = kk |
---|
1095 | k = k + 1 |
---|
1096 | kk = factor(k) |
---|
1097 | jj = jj + kk |
---|
1098 | END DO |
---|
1099 | perm (j) = jj |
---|
1100 | END DO |
---|
1101 | ! |
---|
1102 | !-- Determine the permutation cycles of length greater than 1 |
---|
1103 | j = 0 |
---|
1104 | DO |
---|
1105 | DO |
---|
1106 | j = j + 1 |
---|
1107 | kk = perm(j) |
---|
1108 | IF (kk >= 0) EXIT |
---|
1109 | END DO |
---|
1110 | IF (kk /= j) THEN |
---|
1111 | DO |
---|
1112 | k = kk |
---|
1113 | kk = perm (k) |
---|
1114 | perm (k) = -kk |
---|
1115 | IF (kk == j) EXIT |
---|
1116 | END DO |
---|
1117 | k3 = kk |
---|
1118 | ELSE |
---|
1119 | perm (j) = -j |
---|
1120 | IF (j == nn) EXIT !-- exit infinite loop |
---|
1121 | END IF |
---|
1122 | END DO |
---|
1123 | ! |
---|
1124 | !-- Reorder a and b, following the permutation cycles |
---|
1125 | DO |
---|
1126 | j = k3 + 1 |
---|
1127 | nt = nt - ispan |
---|
1128 | ii = nt - 1 + 1 |
---|
1129 | IF (nt < 0) EXIT !-- exit infinite loop |
---|
1130 | DO |
---|
1131 | DO |
---|
1132 | j = j-1 |
---|
1133 | IF (perm(j) >= 0) EXIT |
---|
1134 | END DO |
---|
1135 | jj = jc |
---|
1136 | DO |
---|
1137 | kspan = jj |
---|
1138 | IF (jj > maxfactor) kspan = maxfactor |
---|
1139 | jj = jj - kspan |
---|
1140 | k = perm(j) |
---|
1141 | kk = jc * k + ii + jj |
---|
1142 | k1 = kk + kspan |
---|
1143 | k2 = 0 |
---|
1144 | DO |
---|
1145 | k2 = k2 + 1 |
---|
1146 | ctmp(k2) = array(k1) |
---|
1147 | k1 = k1 - 1 |
---|
1148 | IF (k1 == kk) EXIT |
---|
1149 | END DO |
---|
1150 | DO |
---|
1151 | k1 = kk + kspan |
---|
1152 | k2 = k1 - jc * (k + perm(k)) |
---|
1153 | k = -perm(k) |
---|
1154 | DO |
---|
1155 | array(k1) = array(k2) |
---|
1156 | k1 = k1 - 1 |
---|
1157 | k2 = k2 - 1 |
---|
1158 | IF (k1 == kk) EXIT |
---|
1159 | END DO |
---|
1160 | kk = k2 |
---|
1161 | IF (k == j) EXIT |
---|
1162 | END DO |
---|
1163 | k1 = kk + kspan |
---|
1164 | k2 = 0 |
---|
1165 | DO |
---|
1166 | k2 = k2 + 1 |
---|
1167 | array(k1) = ctmp(k2) |
---|
1168 | k1 = k1 - 1 |
---|
1169 | IF (k1 == kk) EXIT |
---|
1170 | END DO |
---|
1171 | IF (jj == 0) EXIT |
---|
1172 | END DO |
---|
1173 | IF (j == 1) EXIT |
---|
1174 | END DO |
---|
1175 | END DO |
---|
1176 | |
---|
1177 | END SUBROUTINE permute |
---|
1178 | |
---|
1179 | END SUBROUTINE fftradix |
---|
1180 | |
---|
1181 | END MODULE singleton |
---|