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

Last change on this file since 1745 was 1683, checked in by knoop, 9 years ago

last commit documented

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