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

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

Initial repository layout and content

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