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

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

last commit documented

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