!> @file singleton_mod.f90
!--------------------------------------------------------------------------------------------------!
! This file is part of the PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
! Public License as published by the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
! Public License for more details.
!
! You should have received a copy of the GNU General Public License along with PALM. If not, see
! .
!
! Copyright 1997-2021 Leibniz Universitaet Hannover
!--------------------------------------------------------------------------------------------------!
!
!
! Current revisions:
! -----------------
!
!
! Former revisions:
! -----------------
! $Id: singleton_mod.f90 4828 2021-01-05 11:21:41Z moh.hefny $
! File re-formatted to follow the PALM coding standard
!
!
! 4182 2019-08-22 15:20:23Z scharf
! Corrected "Former revisions" section
!
! 3761 2019-02-25 15:31:42Z raasch
! Statement added to prevent compiler warning about unused variables
!
! Revision 1.1 2002/05/02 18:56:59 raasch
! Initial revision
!
!
!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> Multivariate Fast Fourier Transform
!>
!> Fortran 90 Implementation of Singleton's mixed-radix algorithm, RC Singleton, Stanford Research
!> Institute, Sept. 1968.
!>
!> Adapted from fftn.c, translated from Fortran 66 to C by Mark Olesen and John Beale.
!>
!> Fourier transforms can be computed either in place, using assumed size arguments, or by generic
!> function, using assumed shape arguments.
!>
!> Public:
!>
!> fftkind kind parameter of complex arguments and function results.
!>
!> fft(array, dim, inv, stat) generic transform function
!> COMPLEX(fftkind), DIMENSION(:,...,:), INTENT(IN) :: array
!> INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim
!> LOGICAL, INTENT(IN), OPTIONAL:: inv
!> INTEGER, INTENT(OUT), OPTIONAL:: stat
!>
!> fftn(array, shape, dim, inv, stat) in place transform subroutine
!> COMPLEX(fftkind), DIMENSION(*), INTENT(INOUT) :: array
!> INTEGER, DIMENSION(:), INTENT(IN) :: shape
!> INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim
!> LOGICAL, INTENT(IN), OPTIONAL:: inv
!> INTEGER, INTENT(OUT), OPTIONAL:: stat
!>
!>
!> Formal Parameters:
!>
!> array The complex array to be transformed. Array can be of arbitrary rank (i.e. up to seven).
!>
!> shape With subroutine fftn, the shape of the array to be transformed has to be passed
!> separately, since fftradix - the internal transformation routine - will always treat
!> array as one dimensional. The product of elements in shape must be the number of
!> elements in array.
!> Although passing array with assumed shape would have been nicer, I prefered assumed
!> size in order to prevent the compiler from using a copy-in-copy-out mechanism. That
!> would generally be necessary with fftn passing array to fftradix and with fftn being
!> prepared for accepting non consecutive array sections. Using assumed size, it's up to
!> the user to pass an array argument, that can be addressed as continous one dimensional
!> array without copying. Otherwise, transformation will not really be performed in place.
!> On the other hand, since the rank of array and the size of shape needn't match, fftn
!> is appropriate for handling more than seven dimensions. As far as function fft is
!> concerned all this doesn't matter, because the argument will be copied anyway. Thus no
!> extra shape argument is needed for fft.
!>
!> Optional Parameters:
!>
!> dim One dimensional integer array, containing the dimensions to be transformed. Default
!> is (/1,...,N/) with N being the rank of array, i.e. complete transform. dim can
!> restrict transformation to a subset of available dimensions. Its size must not exceed
!> the rank of array or the size of shape respectivly.
!>
!> inv If .true., inverse transformation will be performed. Default is .false., i.e. forward
!> transformation.
!>
!> stat If present, a system dependent nonzero status value will be returned in stat, if
!> allocation of temporary storage failed.
!>
!>
!> Scaling:
!>
!> Transformation results will always be scaled by the square root of the product of sizes of each
!> dimension in dim. (See examples below)
!>
!>
!> Examples:
!>
!> Let A be a L*M*N three dimensional complex array. Then
!>
!> result = fft(A)
!>
!> will produce a three dimensional transform, scaled by sqrt(L*M*N), while
!>
!> call fftn(A, SHAPE(A))
!>
!> will do the same in place.
!>
!> result = fft(A, dim=(/1,3/))
!>
!> will transform with respect to the first and the third dimension, scaled by sqrt(L*N).
!>
!> result = fft(fft(A), inv=.true.)
!>
!> should (approximately) reproduce A.
!> With B having the same shape as A
!>
!> result = fft(fft(A) * CONJG(fft(B)), inv=.true.)
!>
!> will correlate A and B.
!>
!>
!> Remarks:
!>
!> Following changes have been introduced with respect to fftn.c:
!> - Complex arguments and results are of type complex, rather than real an imaginary part
!> separately.
!> - Increment parameter (magnitude of isign) has been dropped, inc is always one, direction of
!> transform is given by inv.
!> - maxf and maxp have been dropped. The amount of temporary storage needed is determined by the
!> fftradix routine. Both fftn and fft can handle any size of array. (Maybe they take a lot of
!> time and memory, but they will do it)
!>
!> Redesigning fftradix in a way, that it handles assumed shape arrays would have been desirable.
!> However, I found it rather hard to do this in an efficient way. Problems were:
!> - To prevent stride multiplications when indexing arrays. At least our compiler was not clever
!> enough to discover that in fact additions would do the job as well. On the other hand, I
!> haven't been clever enough to find an implementation using array operations.
!> - fftradix is rather large and different versions would be necessaray for each possible rank of
!> array.
!> Consequently, in place transformation still needs the argument stored in a consecutive bunch of
!> memory and can't be performed on array sections like A(100:199:-3, 50:1020). Calling fftn with
!> such sections will most probably imply copy-in-copy-out. However, the function fft works with
!> everything it gets and should be convenient to use.
!>
!> Michael Steffens, 09.12.96,
!> Restructured fftradix for better optimization. M. Steffens, 4 June 1997
!--------------------------------------------------------------------------------------------------!
MODULE singleton
USE kinds
IMPLICIT NONE
PRIVATE
PUBLIC :: fft !<
PUBLIC :: fftn !<
REAL(wp), PARAMETER :: cos72 = 0.30901699437494742_wp !<
REAL(wp), PARAMETER :: pi = 3.14159265358979323_wp !<
REAL(wp), PARAMETER :: sin60 = 0.86602540378443865_wp !<
REAL(wp), PARAMETER :: sin72 = 0.95105651629515357_wp !<
INTERFACE fft
MODULE PROCEDURE fft1d
MODULE PROCEDURE fft2d
MODULE PROCEDURE fft3d
MODULE PROCEDURE fft4d
MODULE PROCEDURE fft5d
MODULE PROCEDURE fft6d
MODULE PROCEDURE fft7d
END INTERFACE
CONTAINS
!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> @todo Missing function description.
!--------------------------------------------------------------------------------------------------!
FUNCTION fft1d( array, dim, inv, stat ) RESULT( ft )
!
!-- Formal parameters
COMPLEX(wp), DIMENSION(:), INTENT(IN) :: array !<
INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL :: dim !<
INTEGER(iwp), INTENT(OUT), OPTIONAL :: stat !<
LOGICAL, INTENT(IN), OPTIONAL :: inv !<
!
!-- Function result
COMPLEX(wp), DIMENSION(SIZE(array, 1)) :: ft !<
INTEGER(iwp) :: idum !<
INTEGER(iwp) :: ishape(1) !<
!
!-- Intrinsics used
INTRINSIC SIZE, SHAPE
ft = array
ishape = SHAPE( array )
CALL fftn( ft, ishape, inv = inv, stat = stat )
!
!-- Next statement to prevent compiler warning about unused variable
IF ( PRESENT( dim ) ) idum = 1
END FUNCTION fft1d
!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> @todo Missing function description.
!--------------------------------------------------------------------------------------------------!
FUNCTION fft2d( array, dim, inv, stat ) RESULT( ft )
!
!-- Formal parameters
COMPLEX(wp), DIMENSION(:,:), INTENT(IN) :: array !<
INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL :: dim !<
INTEGER(iwp), INTENT(OUT), OPTIONAL :: stat !<
LOGICAL, INTENT(IN), OPTIONAL :: inv !<
!
!-- Function result
COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2)) :: ft !<
INTEGER(iwp) :: ishape(2) !<
!
!-- Intrinsics used
INTRINSIC SIZE, SHAPE
ft = array
ishape = SHAPE( array )
CALL fftn( ft, ishape, dim, inv, stat )
END FUNCTION fft2d
!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> @todo Missing function description.
!--------------------------------------------------------------------------------------------------!
FUNCTION fft3d( array, dim, inv, stat ) RESULT( ft )
!
!-- Formal parameters
COMPLEX(wp), DIMENSION(:,:,:), INTENT(IN) :: array !<
INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL :: dim !<
INTEGER(iwp), INTENT(OUT), OPTIONAL :: stat !<
LOGICAL, INTENT(IN), OPTIONAL :: inv !<
!
!-- Function result
COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3)) :: ft !<
INTEGER(iwp) :: ishape(3) !<
!
!-- Intrinsics used
INTRINSIC SIZE, SHAPE
ft = array
ishape = SHAPE( array)
CALL fftn(ft, ishape, dim, inv, stat)
END FUNCTION fft3d
!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> @todo Missing function description.
!--------------------------------------------------------------------------------------------------!
FUNCTION fft4d( array, dim, inv, stat ) RESULT( ft )
!
!-- Formal parameters
COMPLEX(wp), DIMENSION(:,:,:,:), INTENT(IN) :: array !<
INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL :: dim !<
INTEGER(iwp), INTENT(OUT), OPTIONAL :: stat !<
LOGICAL, INTENT(IN), OPTIONAL :: inv !<
!
!-- Function result
COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4)) :: ft !<
INTEGER(iwp) :: ishape(4) !<
!
!-- Intrinsics used
INTRINSIC SIZE, SHAPE
ft = array
ishape = SHAPE( array )
CALL fftn(ft, ishape, dim, inv, stat)
END FUNCTION fft4d
!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> @todo Missing function description.
!--------------------------------------------------------------------------------------------------!
FUNCTION fft5d( array, dim, inv, stat ) RESULT( ft )
!
!-- Formal parameters
COMPLEX(wp), DIMENSION(:,:,:,:,:), INTENT(IN) :: array !<
INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL :: dim !<
INTEGER(iwp), INTENT(OUT), OPTIONAL :: stat !<
LOGICAL, INTENT(IN), OPTIONAL :: inv !<
!
!-- Function result
COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), SIZE(array, 5)) :: ft !<
INTEGER(iwp) :: ishape(5) !<
!
!-- Intrinsics used
INTRINSIC SIZE, SHAPE
ft = array
ishape = SHAPE( array )
CALL fftn(ft, ishape, dim, inv, stat)
END FUNCTION fft5d
!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> @todo Missing function description.
!--------------------------------------------------------------------------------------------------!
FUNCTION fft6d( array, dim, inv, stat ) RESULT( ft )
!
!-- Formal parameters
COMPLEX(wp), DIMENSION(:,:,:,:,:,:), INTENT(IN) :: array !<
INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL :: dim !<
INTEGER(iwp), INTENT(OUT), OPTIONAL :: stat !<
LOGICAL, INTENT(IN), OPTIONAL :: inv !<
!
!-- Function result
COMPLEX(wp), DIMENSION( SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), &
SIZE(array, 5), SIZE(array, 6)) :: ft !<
INTEGER(iwp) :: ishape(6) !<
!
!-- Intrinsics used
INTRINSIC SIZE, SHAPE
ft = array
ishape = SHAPE( array )
CALL fftn(ft, ishape, dim, inv, stat)
END FUNCTION fft6d
!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> @todo Missing function description.
!--------------------------------------------------------------------------------------------------!
FUNCTION fft7d( array, dim, inv, stat ) RESULT( ft )
!
!-- Formal parameters
COMPLEX(wp), DIMENSION(:,:,:,:,:,:,:), INTENT(IN) :: array !<
INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL :: dim !<
INTEGER(iwp), INTENT(OUT), OPTIONAL :: stat !<
LOGICAL, INTENT(IN), OPTIONAL :: inv !<
!
!-- Function result
COMPLEX(wp), DIMENSION( SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), &
SIZE(array, 5), SIZE(array, 6), SIZE(array, 7)) :: ft !<
INTEGER(iwp) :: ishape(7) !<
!
!-- Intrinsics used
INTRINSIC SIZE, SHAPE
ft = array
ishape = SHAPE( array )
CALL fftn(ft, ishape, dim, inv, stat)
END FUNCTION fft7d
!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> @todo Missing subroutine description.
!--------------------------------------------------------------------------------------------------!
SUBROUTINE fftn( array, shape, dim, inv, stat )
!
!-- Formal parameters
COMPLEX(wp), DIMENSION(*), INTENT(INOUT) :: array !<
INTEGER(iwp), DIMENSION(:), INTENT(IN) :: shape !<
INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL :: dim !<
INTEGER(iwp), INTENT(OUT), OPTIONAL :: stat !<
LOGICAL, INTENT(IN), OPTIONAL :: inv !<
!
!-- Local arrays
INTEGER(iwp), DIMENSION(SIZE(shape)) :: d !<
!
!-- Local scalars
LOGICAL :: inverse !<
INTEGER(iwp) :: i, ndim, ntotal !<
REAL(wp) :: scale !<
!
!-- Intrinsics used
INTRINSIC PRESENT, MIN, PRODUCT, SIZE, SQRT
!
!-- Optional parameter settings
IF ( PRESENT( inv ) ) THEN
inverse = inv
ELSE
inverse = .FALSE.
END IF
IF ( PRESENT( dim ) ) THEN
ndim = MIN( SIZE( dim ), SIZE( d ) )
d(1:ndim) = DIM( 1:ndim )
ELSE
ndim = SIZE( d )
d = (/( i, i = 1, SIZE( d ) )/)
END IF
ntotal = PRODUCT( shape )
scale = SQRT( 1.0_wp / PRODUCT( shape( d(1:ndim) ) ) )
DO i = 1, ntotal
array(i) = CMPLX( REAL( array(i) ) * scale, AIMAG( array(i) ) * scale, KIND = wp )
END DO
DO i = 1, ndim
CALL fftradix( array, ntotal, shape( d(i) ), PRODUCT( shape( 1:d(i) ) ), inverse, stat )
IF ( PRESENT( stat ) ) THEN
IF ( stat /= 0 ) RETURN
END IF
END DO
END SUBROUTINE fftn
!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> @todo Missing subroutine description.
!--------------------------------------------------------------------------------------------------!
SUBROUTINE fftradix( array, ntotal, npass, nspan, inv, stat )
!
!-- Formal parameters
COMPLEX(wp), DIMENSION(*), INTENT(INOUT) :: array !<
INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan !<
INTEGER(iwp), INTENT(OUT), OPTIONAL :: stat !<
LOGICAL, INTENT(IN) :: inv !<
!
!-- Local arrays
COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: ctmp !<
INTEGER(iwp), DIMENSION(BIT_SIZE(0)) :: factor !<
INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: perm !<
REAL(wp), DIMENSION(:), ALLOCATABLE :: sine, cosine !<
!
!-- Local scalars
INTEGER(iwp) :: maxfactor, nfactor, nsquare, nperm !<
!
!-- Intrinsics used
INTRINSIC MAXVAL, MOD, PRESENT, ISHFT, BIT_SIZE, SIN, COS, CMPLX, REAL, AIMAG
IF ( npass <= 1 ) RETURN
CALL factorize( npass, factor, nfactor, nsquare )
maxfactor = MAXVAL( factor(:nfactor) )
IF ( nfactor - ISHFT( nsquare, 1 ) > 0 ) THEN
nperm = MAX( nfactor + 1, PRODUCT( factor(nsquare+1: nfactor-nsquare) ) - 1 )
ELSE
nperm = nfactor + 1
END IF
IF ( PRESENT( stat ) ) THEN
ALLOCATE( ctmp(maxfactor), sine(maxfactor), cosine(maxfactor), STAT = stat )
IF ( stat /= 0 ) RETURN
CALL transform( array, ntotal, npass, nspan, factor, nfactor, ctmp, sine, cosine, inv )
DEALLOCATE( sine, cosine, STAT = stat )
IF ( stat /= 0 ) RETURN
ALLOCATE( perm(nperm), STAT = stat )
IF ( stat /= 0 ) RETURN
CALL permute( array, ntotal, npass, nspan, factor, nfactor, nsquare, maxfactor, ctmp, perm )
DEALLOCATE( perm, ctmp, STAT = stat )
IF ( stat /= 0 ) RETURN
ELSE
ALLOCATE( ctmp(maxfactor), sine(maxfactor), cosine(maxfactor) )
CALL transform( array, ntotal, npass, nspan, factor, nfactor, ctmp, sine, cosine, inv )
DEALLOCATE( sine, cosine )
ALLOCATE( perm(nperm) )
CALL permute( array, ntotal, npass, nspan, factor, nfactor, nsquare, maxfactor, ctmp, perm )
DEALLOCATE( perm, ctmp )
END IF
CONTAINS
!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> @todo Missing subroutine description.
!--------------------------------------------------------------------------------------------------!
SUBROUTINE factorize( npass, factor, nfactor, nsquare )
!
!-- Formal parameters
INTEGER(iwp), INTENT(IN) :: npass !<
INTEGER(iwp), INTENT(OUT) :: nfactor, nsquare !<
INTEGER(iwp), DIMENSION(*), INTENT(OUT) :: factor !<
!
!-- Local scalars
INTEGER(iwp) :: j, jj, k !<
nfactor = 0
k = npass
DO WHILE ( MOD( k, 16 ) == 0 )
nfactor = nfactor + 1
factor(nfactor) = 4
k = k / 16
END DO
j = 3
jj = 9
DO
DO WHILE ( MOD( k, jj ) == 0 )
nfactor = nfactor + 1
factor(nfactor) = j
k = k / jj
END DO
j = j + 2
jj = j * j
IF ( jj > k ) EXIT
END DO
IF ( k <= 4 ) THEN
nsquare = nfactor
factor(nfactor + 1) = k
IF ( k /= 1 ) nfactor = nfactor + 1
ELSE
IF ( k - ISHFT( k / 4, 2 ) == 0 ) THEN
nfactor = nfactor + 1
factor(nfactor) = 2
k = k / 4
END IF
nsquare = nfactor
j = 2
DO
IF ( MOD(k, j) == 0 ) THEN
nfactor = nfactor + 1
factor(nfactor) = j
k = k / j
END IF
j = ISHFT( (j + 1) / 2, 1 ) + 1
IF ( j > k ) EXIT
END DO
END IF
IF ( nsquare > 0 ) THEN
j = nsquare
DO
nfactor = nfactor + 1
factor(nfactor) = factor(j)
j = j - 1
IF ( j == 0 ) EXIT
END DO
END IF
END SUBROUTINE factorize
!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> @todo Missing subroutine description.
!--------------------------------------------------------------------------------------------------!
SUBROUTINE transform( array, ntotal, npass, nspan, factor, nfactor, ctmp, sine, cosine, inv )
!-- Compute fourier transform
!
!-- Formal parameters
COMPLEX(wp), DIMENSION(*), INTENT(IN OUT) :: array !<
COMPLEX(wp), DIMENSION(*), INTENT(OUT) :: ctmp !<
INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan !<
INTEGER(iwp), DIMENSION(*), INTENT(IN) :: factor !<
INTEGER(iwp), INTENT(IN) :: nfactor !<
LOGICAL, INTENT(IN) :: inv !<
REAL(wp), DIMENSION(*), INTENT(OUT) :: sine, cosine !<
!
!-- Local scalars
COMPLEX(wp) :: cc, cj, ck, cjp, cjm, ckp, ckm !<
INTEGER(iwp) :: ii, ispan !<
INTEGER(iwp) :: j, jc, jf, jj !<
INTEGER(iwp) :: k, kk, kspan, k1, k2, k3, k4 !<
INTEGER(iwp) :: nn, nt !<
REAL(wp) :: s60, c72, s72, pi2, radf !<
REAL(wp) :: c1, s1, c2, s2, c3, s3, cd, sd, ak !<
c72 = cos72
IF ( inv ) THEN
s72 = sin72
s60 = sin60
pi2 = pi
ELSE
s72 = - sin72
s60 = - sin60
pi2 = - pi
END IF
nt = ntotal
nn = nt - 1
kspan = nspan
jc = nspan / npass
radf = pi2 * jc
pi2 = pi2 * 2.0_wp !-- Use 2 PI from here on
ii = 0
jf = 0
DO
sd = radf / kspan
cd = SIN( sd )
cd = 2.0_wp * cd * cd
sd = SIN( sd + sd )
kk = 1
ii = ii + 1
SELECT CASE ( factor(ii) )
CASE ( 2 )
!
!-- Transform for factor of 2 (including rotation factor)
kspan = kspan / 2
k1 = kspan + 2
DO
DO
k2 = kk + kspan
ck = array(k2)
array(k2) = array(kk) - ck
array(kk) = array(kk) + ck
kk = k2 + kspan
IF ( kk > nn ) EXIT
END DO
kk = kk - nn
IF ( kk > jc ) EXIT
END DO
IF ( kk > kspan ) RETURN
DO
c1 = 1.0_wp - cd
s1 = sd
DO
DO
DO
k2 = kk + kspan
ck = array(kk) - array(k2)
array(kk) = array(kk) + array(k2)
array(k2) = ck * CMPLX( c1, s1, KIND = wp )
kk = k2 + kspan
IF ( kk >= nt ) EXIT
END DO
k2 = kk - nt
c1 = - c1
kk = k1 - k2
IF ( kk <= k2 ) EXIT
END DO
ak = c1 - (cd * c1 + sd * s1)
s1 = sd * c1 - cd * s1 + s1
c1 = 2.0_wp - ( ak * ak + s1 * s1 )
s1 = s1 * c1
c1 = c1 * ak
kk = kk + jc
IF ( kk >= k2 ) EXIT
END DO
k1 = k1 + 1 + 1
kk = ( k1 - kspan ) / 2 + jc
IF ( kk > jc + jc ) EXIT
END DO
!
!-- Transform for factor of 4
CASE ( 4 )
ispan = kspan
kspan = kspan / 4
DO
c1 = 1.0_wp
s1 = 0.0_wp
DO
DO
k1 = kk + kspan
k2 = k1 + kspan
k3 = k2 + kspan
ckp = array(kk) + array(k2)
ckm = array(kk) - array(k2)
cjp = array(k1) + array(k3)
cjm = array(k1) - array(k3)
array(kk) = ckp + cjp
cjp = ckp - cjp
IF ( inv ) THEN
ckp = ckm + CMPLX( - AIMAG( cjm ), REAL( cjm ), KIND = wp )
ckm = ckm + CMPLX( AIMAG( cjm ), - REAL( cjm ), KIND = wp )
ELSE
ckp = ckm + CMPLX( AIMAG( cjm ), - REAL( cjm ), KIND = wp )
ckm = ckm + CMPLX( - AIMAG( cjm ), REAL( cjm ), KIND = wp )
END IF
!
!-- Avoid useless multiplies
IF ( s1 == 0.0_wp ) THEN
array(k1) = ckp
array(k2) = cjp
array(k3) = ckm
ELSE
array(k1) = ckp * CMPLX( c1, s1, KIND = wp )
array(k2) = cjp * CMPLX( c2, s2, KIND = wp )
array(k3) = ckm * CMPLX( c3, s3, KIND = wp )
END IF
kk = k3 + kspan
IF ( kk > nt ) EXIT
END DO
c2 = c1 - ( cd * c1 + sd * s1 )
s1 = sd * c1 - cd * s1 + s1
c1 = 2.0_wp - ( c2 * c2 + s1 * s1 )
s1 = s1 * c1
c1 = c1 * c2
!
!-- Values of c2, c3, s2, s3 that will get used next time
c2 = c1 * c1 - s1 * s1
s2 = 2.0_wp * c1 * s1
c3 = c2 * c1 - s2 * s1
s3 = c2 * s1 + s2 * c1
kk = kk - nt + jc
IF ( kk > kspan ) EXIT
END DO
kk = kk - kspan + 1
IF ( kk > jc ) EXIT
END DO
IF ( kspan == jc ) RETURN
CASE default
!
!-- Transform for odd factors
k = factor(ii)
ispan = kspan
kspan = kspan / k
SELECT CASE ( k )
!
!-- Transform for factor of 3 (optional code)
CASE ( 3 )
DO
DO
k1 = kk + kspan
k2 = k1 + kspan
ck = array(kk)
cj = array(k1) + array(k2)
array(kk) = ck + cj
ck = ck - CMPLX( 0.5_wp * REAL( cj ), 0.5_wp * AIMAG( cj ), KIND = wp )
cj = CMPLX( ( REAL( array(k1) ) - REAL( array(k2) ) ) * s60, &
( AIMAG( array(k1) ) - AIMAG( array(k2) ) ) * s60, KIND = wp )
array(k1) = ck + CMPLX( - AIMAG( cj ), REAL( cj ), KIND = wp )
array(k2) = ck + CMPLX( AIMAG( cj ), - REAL( cj ), KIND = wp )
kk = k2 + kspan
IF ( kk >= nn ) EXIT
END DO
kk = kk - nn
IF ( kk > kspan ) EXIT
END DO
!
!-- Transform for factor of 5 (optional code)
CASE ( 5 )
c2 = c72 * c72 - s72 * s72
s2 = 2.0_wp * c72 * s72
DO
DO
k1 = kk + kspan
k2 = k1 + kspan
k3 = k2 + kspan
k4 = k3 + kspan
ckp = array(k1) + array(k4)
ckm = array(k1) - array(k4)
cjp = array(k2) + array(k3)
cjm = array(k2) - array(k3)
cc = array(kk)
array(kk) = cc + ckp + cjp
ck = CMPLX( REAL( ckp ) * c72, AIMAG( ckp ) * c72, KIND = wp ) + &
CMPLX( REAL( cjp ) * c2, AIMAG( cjp ) * c2, KIND = wp ) + cc
cj = CMPLX( REAL( ckm ) * s72, AIMAG( ckm ) * s72, KIND = wp) + &
CMPLX( REAL( cjm ) * s2, AIMAG( cjm ) * s2, KIND = wp )
array(k1) = ck + CMPLX( - AIMAG( cj ), REAL( cj ), KIND = wp )
array(k4) = ck + CMPLX( AIMAG( cj ), - REAL( cj ), KIND = wp )
ck = CMPLX( REAL( ckp ) * c2, AIMAG( ckp ) * c2, KIND = wp ) + &
CMPLX( REAL( cjp ) * c72, AIMAG( cjp ) * c72, KIND = wp ) + cc
cj = CMPLX( REAL( ckm ) * s2, AIMAG( ckm ) * s2, KIND = wp ) - &
CMPLX( REAL( cjm ) * s72, AIMAG( cjm ) * s72, KIND = wp )
array(k2) = ck + CMPLX( -AIMAG( cj ), REAL( cj ), KIND = wp )
array(k3) = ck + CMPLX( AIMAG( cj ), - REAL( cj ), KIND = wp )
kk = k4 + kspan
IF ( kk >= nn ) EXIT
END DO
kk = kk - nn
IF ( kk > kspan ) EXIT
END DO
CASE default
IF ( k /= jf ) THEN
jf = k
s1 = pi2 / k
c1 = COS( s1 )
s1 = SIN( s1 )
cosine(jf) = 1.0_wp
sine(jf) = 0.0_wp
j = 1
DO
cosine(j) = cosine(k) * c1 + sine(k) * s1
sine(j) = cosine(k) * s1 - sine(k) * c1
k = k - 1
cosine(k) = cosine(j)
sine(k) = - sine(j)
j = j + 1
IF ( j >= k ) EXIT
END DO
END IF
DO
DO
k1 = kk
k2 = kk + ispan
cc = array(kk)
ck = cc
j = 1
k1 = k1 + kspan
DO
k2 = k2 - kspan
j = j + 1
ctmp(j) = array(k1) + array(k2)
ck = ck + ctmp(j)
j = j + 1
ctmp(j) = array(k1) - array(k2)
k1 = k1 + kspan
IF ( k1 >= k2 ) EXIT
END DO
array(kk) = ck
k1 = kk
k2 = kk + ispan
j = 1
DO
k1 = k1 + kspan
k2 = k2 - kspan
jj = j
ck = cc
cj = ( 0.0_wp, 0.0_wp )
k = 1
DO
k = k + 1
ck = ck + CMPLX( REAL( ctmp(k) ) * cosine(jj), AIMAG( ctmp(k) ) * &
cosine(jj), KIND = wp )
k = k + 1
cj = cj + CMPLX( REAL( ctmp(k) ) * sine(jj), AIMAG( ctmp(k) ) * sine(jj), &
KIND = wp )
jj = jj + j
IF ( jj > jf ) jj = jj - jf
IF ( k >= jf ) EXIT
END DO
k = jf - j
array(k1) = ck + CMPLX( - AIMAG( cj ), REAL( cj ), KIND = wp )
array(k2) = ck + CMPLX( AIMAG( cj ), -REAL( cj ), KIND = wp )
j = j + 1
IF ( j >= k ) EXIT
END DO
kk = kk + ispan
IF ( kk > nn ) EXIT
END DO
kk = kk - nn
IF ( kk > kspan ) EXIT
END DO
END SELECT
!
!-- Multiply by rotation factor (except for factors of 2 and 4)
IF ( ii == nfactor ) RETURN
kk = jc + 1
DO
c2 = 1.0_wp - cd
s1 = sd
DO
c1 = c2
s2 = s1
kk = kk + kspan
DO
DO
array(kk) = CMPLX( c2, s2, KIND = wp ) * array(kk)
kk = kk + ispan
IF ( kk > nt ) EXIT
END DO
ak = s1 * s2
s2 = s1 * c2 + c1 * s2
c2 = c1 * c2 - ak
kk = kk - nt + kspan
IF ( kk > ispan ) EXIT
END DO
c2 = c1 - ( cd * c1 + sd * s1 )
s1 = s1 + sd * c1 - cd * s1
c1 = 2.0_wp - ( c2 * c2 + s1 * s1 )
s1 = s1 * c1
c2 = c2 * c1
kk = kk - ispan + jc
IF ( kk > kspan ) EXIT
END DO
kk = kk - kspan + jc + 1
IF ( kk > jc + jc ) EXIT
END DO
END SELECT
END DO
END SUBROUTINE transform
!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> @todo Missing subroutine description.
!--------------------------------------------------------------------------------------------------!
SUBROUTINE permute( array, ntotal, npass, nspan, factor, nfactor, nsquare, maxfactor, ctmp, perm )
!
!-- Formal parameters
COMPLEX(wp), DIMENSION(*), INTENT(IN OUT) :: array !<
COMPLEX(wp), DIMENSION(*), INTENT(OUT) :: ctmp !<
INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan !<
INTEGER(iwp), INTENT(IN) :: nfactor, nsquare !<
INTEGER(iwp), INTENT(IN) :: maxfactor !<
INTEGER(iwp), DIMENSION(*), INTENT(IN OUT) :: factor !<
INTEGER(iwp), DIMENSION(*), INTENT(OUT) :: perm !<
!
!-- Local scalars
COMPLEX(wp) :: ck !<
INTEGER(iwp) :: ii, ispan !<
INTEGER(iwp) :: j, jc, jj !<
INTEGER(iwp) :: k, kk, kspan, kt, k1, k2, k3 !<
INTEGER(iwp) :: nn, nt !<
!
!-- Permute the results to normal order---done in two stages
!-- Permutation for square factors of n
nt = ntotal
nn = nt - 1
kt = nsquare
kspan = nspan
jc = nspan / npass
perm (1) = nspan
IF ( kt > 0 ) THEN
k = kt + kt + 1
IF ( nfactor < k ) k = k - 1
j = 1
perm(k + 1) = jc
DO
perm(j + 1) = perm(j) / factor(j)
perm(k) = perm(k + 1) * factor(j)
j = j + 1
k = k - 1
IF ( j >= k ) EXIT
END DO
k3 = perm(k + 1)
kspan = perm(2)
kk = jc + 1
k2 = kspan + 1
j = 1
IF ( npass /= ntotal ) THEN
permute_multi: DO
DO
DO
k = kk + jc
DO
!
!-- Swap array(kk) <> array(k2)
ck = array(kk)
array(kk) = array(k2)
array(k2) = ck
kk = kk + 1
k2 = k2 + 1
IF ( kk >= k ) EXIT
END DO
kk = kk + nspan - jc
k2 = k2 + nspan - jc
IF ( kk >= nt ) EXIT
END DO
kk = kk - nt + jc
k2 = k2 - nt + kspan
IF ( k2 >= nspan ) EXIT
END DO
DO
DO
k2 = k2 - perm(j)
j = j + 1
k2 = perm(j + 1) + k2
IF ( k2 <= perm(j) ) EXIT
END DO
j = 1
DO
IF ( kk < k2 ) CYCLE permute_multi
kk = kk + jc
k2 = k2 + kspan
IF ( k2 >= nspan ) EXIT
END DO
IF ( kk >= nspan ) EXIT
END DO
EXIT
END DO permute_multi
ELSE
permute_single: DO
DO
!
!-- Swap array(kk) <> array(k2)
ck = array(kk)
array(kk) = array(k2)
array(k2) = ck
kk = kk + 1
k2 = k2 + kspan
IF ( k2 >= nspan ) EXIT
END DO
DO
DO
k2 = k2 - perm(j)
j = j + 1
k2 = perm(j + 1) + k2
IF ( k2 <= perm(j) ) EXIT
END DO
j = 1
DO
IF ( kk < k2 ) CYCLE permute_single
kk = kk + 1
k2 = k2 + kspan
IF ( k2 >= nspan ) EXIT
END DO
IF ( kk >= nspan ) EXIT
END DO
EXIT
END DO permute_single
END IF
jc = k3
END IF
IF ( ISHFT( kt, 1 ) + 1 >= nfactor ) RETURN
ispan = perm(kt + 1)
!
!-- Permutation for square-free factors of n
j = nfactor - kt
factor( j + 1 ) = 1
DO
factor(j) = factor(j) * factor(j+1)
j = j - 1
IF ( j == kt ) EXIT
END DO
kt = kt + 1
nn = factor( kt ) - 1
j = 0
jj = 0
DO
k = kt + 1
k2 = factor(kt)
kk = factor(k)
j = j + 1
IF ( j > nn ) EXIT !-- Exit infinite loop
jj = jj + kk
DO WHILE ( jj >= k2 )
jj = jj - k2
k2 = kk
k = k + 1
kk = factor(k)
jj = jj + kk
END DO
perm(j) = jj
END DO
!
!-- Determine the permutation cycles of length greater than 1
j = 0
DO
DO
j = j + 1
kk = perm(j)
IF ( kk >= 0 ) EXIT
END DO
IF ( kk /= j ) THEN
DO
k = kk
kk = perm(k)
perm(k) = - kk
IF ( kk == j ) EXIT
END DO
k3 = kk
ELSE
perm(j) = - j
IF ( j == nn ) EXIT !-- Exit infinite loop
END IF
END DO
!
!-- Reorder a and b, following the permutation cycles
DO
j = k3 + 1
nt = nt - ispan
ii = nt - 1 + 1
IF ( nt < 0 ) EXIT !-- Exit infinite loop
DO
DO
j = j - 1
IF ( perm(j) >= 0 ) EXIT
END DO
jj = jc
DO
kspan = jj
IF ( jj > maxfactor ) kspan = maxfactor
jj = jj - kspan
k = perm(j)
kk = jc * k + ii + jj
k1 = kk + kspan
k2 = 0
DO
k2 = k2 + 1
ctmp(k2) = array(k1)
k1 = k1 - 1
IF ( k1 == kk ) EXIT
END DO
DO
k1 = kk + kspan
k2 = k1 - jc * ( k + perm(k) )
k = - perm(k)
DO
array(k1) = array(k2)
k1 = k1 - 1
k2 = k2 - 1
IF ( k1 == kk ) EXIT
END DO
kk = k2
IF ( k == j ) EXIT
END DO
k1 = kk + kspan
k2 = 0
DO
k2 = k2 + 1
array(k1) = ctmp(k2)
k1 = k1 - 1
IF ( k1 == kk ) EXIT
END DO
IF ( jj == 0 ) EXIT
END DO
IF ( j == 1 ) EXIT
END DO
END DO
END SUBROUTINE permute
END SUBROUTINE fftradix
END MODULE singleton