MODULE singleton !----------------------------------------------------------------------------- ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: singleton.f90 1321 2014-03-20 09:40:40Z witha $ ! ! 1320 2014-03-20 08:40:49Z raasch ! kind-parameters added to all INTEGER and REAL declaration statements, ! kinds are defined in new module kinds, ! revision history before 2012 removed, ! ! 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 trans- ! formation routine - will treat array always 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 argu- ! ment, 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 !----------------------------------------------------------------------------- USE kinds IMPLICIT NONE PRIVATE PUBLIC:: fft, fftn REAL(wp), PARAMETER:: sin60 = 0.86602540378443865_wp REAL(wp), PARAMETER:: cos72 = 0.30901699437494742_wp REAL(wp), PARAMETER:: sin72 = 0.95105651629515357_wp REAL(wp), PARAMETER:: pi = 3.14159265358979323_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 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):: ishape(1) ! !-- Intrinsics used INTRINSIC SIZE, SHAPE ft = array ishape = SHAPE( array ) CALL fftn(ft, ishape, inv = inv, stat = stat) END FUNCTION fft1d 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 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 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 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 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 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 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 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 SUBROUTINE factorize(npass, factor, nfactor, nsquare) ! !-- Formal parameters INTEGER(iwp), INTENT(IN) :: npass INTEGER(iwp), DIMENSION(*), INTENT(OUT):: factor INTEGER(iwp), INTENT(OUT):: nfactor, nsquare ! !-- 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 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 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 COMPLEX(wp) :: cc, cj, ck, cjp, cjm, ckp, ckm 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 CASE (4) !-- transform for factor of 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) CASE (3) !-- transform for factor of 3 (optional code) 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 CASE (5) !-- transform for factor of 5 (optional code) 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 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), DIMENSION(*), INTENT(IN OUT):: factor INTEGER(iwp), INTENT(IN) :: nfactor, nsquare INTEGER(iwp), INTENT(IN) :: maxfactor 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