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

Last change on this file since 3761 was 3761, checked in by raasch, 5 years ago

unused variables removed, OpenACC directives re-formatted, statements added to avoid compiler warnings

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