source: palm/trunk/SOURCE/singleton_mod.f90

Last change on this file was 4828, checked in by Giersch, 3 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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