source: palm/trunk/SOURCE/singleton_mod.f90 @ 4180

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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