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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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