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

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

Forced header and separation lines into 80 columns

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