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

Last change on this file since 3046 was 2001, checked in by knoop, 5 years ago

last commit documented

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