source: palm/trunk/SOURCE/singleton.f90 @ 1016

Last change on this file since 1016 was 484, checked in by raasch, 15 years ago

typo in file headers removed

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