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

Last change on this file since 1850 was 1850, checked in by maronga, 5 years ago

added _mod string to several filenames to meet the naming convection for modules

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