Skip to content

GitLab

fleur
fleur
Commit 5a74619b authored by Matthias Redies's avatar Matthias Redies
Browse files
Options

update rfft

parent 765db73a
Hide whitespace changes
Inline Side-by-side
Showing with 1900 additions and 1914 deletions
math/CMakeLists.txt
View file @ 5a74619b
......@@ -16,13 +16,13 @@ math/outint.f
math/points.f
math/qranf.f
math/qsf.f
math/rfft.F
math/sphbes.f
math/sphpts.f
math/util.F
math/difcub.f
)
set(fleur_F90 ${fleur_F90}
math/rfft.f90
math/differentiate.f90
math/fft2d.F90
math/fft3d.f90
......
math/rfft.F deleted 100644 → 0
View file @ 765db73a
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_rfft
use m_juDFT
PRIVATE
!Module also contains routines vrffti,vrfftf&vrfftb as private routines below
PUBLIC rfft
CONTAINS
SUBROUTINE rfft(
> isn,n1d,n2d,n3d,n1,n2,n3,
> nw1,nw2,nw3,wsave,b,
X a)
c **********************************************************************
c
c isn = +1 : FFT a real, 3dim array "a" of dimensions n1*n2*n3 to a complex
c array of size n1*n2*(n3+1)/2 [if n3 is odd] or n1*n2*(n3/2+1)
c [n3 is even].
c isn = -1 : back-FFT of a complex array "a" of size n1*n2*(n3+1)/2 [odd]
c or n1*n2*(n3/2+1) [even] to a real, 3dim array
c
c the actual array is assumed to be located between
c 1 ... nw1+1 & n1-nw1+1 ... n1
c 1 ... nw2+1 & n2-nw2+1 ... n2
c 1 ... nw3+1 & n3-nw3+1 ... n3
c and padded with zeros in between.
c if nw1 >= (n1-1)/2, no padding is assumed (-> aliasing errors!)
c
c G.Bihlmayer (UniWien)
c
c **********************************************************************
USE m_cfft
IMPLICIT NONE
INTEGER n1d,n2d,n3d,n1,n2,n3,nw1,nw2,nw3,isn
REAL a(n1d,n2d,0:n3d),b(n1d,n2d,n3d),wsave(n3d+15)
INTEGER i1,i2,i3,nup
REAL factor
LOGICAL l_nopad
c
c a ... array for FFT
c b ... work array
c wsave ... stores tables for r-2-c FFT
c n1,n2,n3 ... dimensions to be transformed
c nw1,nw2,nw3 ... actual dimensions of a before FFT
c n1d,n2d,n3d ... dimensions of a,b
c
c
c check for input errors
c
IF ((isn/=-1).AND.(isn.NE.1)) CALL juDFT_error("choose isn=+/- 1"
+ ,calledby ="rfft")
IF ((n1d.lt.n1).OR.(n2d.lt.n2).OR.(n3d.lt.n3)) THEN
WRITE (6,*) 'n1d,n2d,n3d =',n1d,n2d,n3d
WRITE (6,*) 'n1 ,n2 ,n3 =',n1 ,n2 ,n3
CALL juDFT_error("n(i) > n(i)d",calledby ="rfft")
ENDIF
IF ((n1.le.2*nw1+1).OR.
+ (n2.le.2*nw2+1).OR.
+ (n3.le.2*nw3+1)) THEN
c WRITE (6,*) 'n1 ,n2 ,n3 =',n1 ,n2 ,n3
c WRITE (6,*) 'nw1,nw2,nw3 =',nw1,nw2,nw3
l_nopad=.true.
ELSE
l_nopad=.false.
ENDIF
c
c ******************** f o r w a r d - t r a n s f o r m *******************
c
IF (isn.eq.1) THEN
c
c array a is assumed to be zero from (1,1,0) to (n1d,n2d,0) and the array
c to be FFT'd starts at (1,1,1) as n1*n2*n3 real numbers.
c first transform n1*n2 real sequences of lenghth n3 to n3/2 complex values
c
CALL vrffti(n3,wsave)
IF (l_nopad) THEN
CALL vrfftf(n1*n2,n3,a(1,1,1),b,n1d*n2d,wsave)
ELSE
DO i2=1,nw2+1
CALL vrfftf(nw1+1,n3,a(1,i2,1),b,n1d*n2d,wsave)
CALL vrfftf(nw1 ,n3,a(n1-nw1+1,i2,1),b,n1d*n2d,wsave)
ENDDO
DO i2=n2-nw2+1,n2
CALL vrfftf(nw1+1,n3,a(1,i2,1),b,n1d*n2d,wsave)
CALL vrfftf(nw1 ,n3,a(n1-nw1+1,i2,1),b,n1d*n2d,wsave)
ENDDO
ENDIF
c
c now we have the FFT'd array stored as described in vrfftf
c (mixed real & compex data)
c remove the norm 1/sqrt(n3) (to be compatible with cfft)
c
factor = sqrt(1.0*n3)
!CALL CPP_BLAS_sscal(n1d*n2d*n3,factor,a(1,1,1),1)
a=a*factor
c
c now, the real part of f(0) has to be moved to a(n1,n2,0) to get a purely
c complex array starting at a(1,1,0)
c
DO i1=1,n1
DO i2=1,n2
a(i1,i2,0)=a(i1,i2,1)
a(i1,i2,1)=0.0
ENDDO
ENDDO
c
ENDIF
c
c ******************** g e n e r a l p a r t *******************
c
c now perform n2*n3/2 and n1*n3/2 complex FFT's; a is assumed to be
c complex and starting at (1,1,0)
c
IF (ABS((n3/2.)-NINT(n3/2.)).gt.0.1) THEN
nup = n3
ELSE
nup = n3+1
IF (n3+1>n3d) CALL juDFT_error("n3 even & n3+1 > n3d" ,calledby
+ ="rfft")
a(:,:,n3+1)=0.0
ENDIF
IF (l_nopad) THEN
DO i3=1,nup,2
CALL cfft(a(1,1,i3-1),a(1,1,i3),n1*n2,n1,n1,isn)
CALL cfft(a(1,1,i3-1),a(1,1,i3),n1*n2,n2,n1*n2,isn)
ENDDO
ELSE
DO i3=1,nup,2
CALL cfft(a(1,1,i3-1),a(1,1,i3),(nw2+1)*n1,n1,n1,isn)
CALL cfft(a(1,n2-nw2+1,i3-1),a(1,n2-nw2+1,i3),
+ nw2*n1,n1,n1,isn)
CALL cfft(a(1,1,i3-1),a(1,1,i3),n1*n2,n2,n1*n2,isn)
ENDDO
ENDIF
c
c ******************** b a c k w a r d - t r a n s f o r m *******************
c
IF (isn.eq.-1) THEN
c
c the real part of f(0) has to be moved to a(n1,n2,1) for a correct
c setup for vrfftb (see comments therein)
c
DO i1=1,n1
DO i2=1,n2
a(i1,i2,1)=a(i1,i2,0)
a(i1,i2,0)=0.0
ENDDO
ENDDO
c
c transform n1*n2 mixed real and complex sequences of lenghth n3/2
c to n3 real values
c
CALL vrffti(n3,wsave)
IF (l_nopad) THEN
CALL vrfftb(n1*n2,n3,a(1,1,1),b,n1d*n2d,wsave)
ELSE
DO i2=1,nw2+1
CALL vrfftb(nw1+1,n3,a(1,i2,1),b,n1d*n2d,wsave)
CALL vrfftb(nw1 ,n3,a(n1-nw1+1,i2,1),b,n1d*n2d,wsave)
ENDDO
DO i2=n2-nw2+1,n2
CALL vrfftb(nw1+1,n3,a(1,i2,1),b,n1d*n2d,wsave)
CALL vrfftb(nw1 ,n3,a(n1-nw1+1,i2,1),b,n1d*n2d,wsave)
ENDDO
ENDIF
c
c remove the norm 1/sqrt(n3) (compatibility with cfft)
c
factor = sqrt(1.0*n3)
!CALL CPP_BLAS_sscal(n1d*n2d*n3,factor,a(1,1,1),1)
a=a*factor
c
ENDIF
END SUBROUTINE rfft
SUBROUTINE vrffti(n,wsave)
C***BEGIN PROLOGUE VRFFTI
C***DATE WRITTEN 860701 (YYMMDD)
C***REVISION DATE 900509 (YYMMDD)
C***CATEGORY NO. J1A1
C***KEYWORDS FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM,
C MULTIPLE SEQUENCES
C***AUTHOR SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST)
C***PURPOSE Initialization for VRFFTF and VRFFTB.
C***DESCRIPTION
C
C Subroutine VRFFTI initializes the array WSAVE which is used in
C both VRFFTF and VRFFTB. The prime factorization of N together with
C a tabulation of certain trigonometric functions are computed and
C stored in the array WSAVE.
C
C Input Parameter
C
C N the length of the sequence to be transformed. There is no
C restriction on N.
C
C Output Parameter
C
C WSAVE a work array which must be dimensioned at least N+15.
C The same work array can be used for both VRFFTF and VRFFTB
C as long as N remains unchanged. Different WSAVE arrays
C are required for different values of N. The contents of
C WSAVE must not be changed between calls of VRFFTF or VRFFTB.
C
C
C * * * * * * * * * * * * * * * * * * * * *
C * *
C * PROGRAM SPECIFICATIONS *
C * *
C * * * * * * * * * * * * * * * * * * * * *
C
C
C Dimension of R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15)
C Arguments
C
C Latest AUGUST 1, 1985
C Revision
C
C Subprograms VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3,
C Required VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2,
C VRADB3, VRADB4, VRADB5, VRADBG, PIMACH
C
C Special NONE
C Conditions
C
C Common NONE
C blocks
C
C I/O NONE
C
C Precision SINGLE
C
C Specialist ROLAND SWEET
C
C Language FORTRAN
C
C History WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE
C NATIONAL BUREAU OF STANDARDS (BOULDER).
C
C Algorithm A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION
C OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM.
C
C Portability AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77.
C THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN
C THE FUNCTION PIMACH.
C
C Required COS,SIN
C resident
C Routines
C
C
C***REFERENCES P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
C Computations, (G. Rodrigue, ed.), Academic Press, 1982,
C pp. 51-83.
C***ROUTINES CALLED VRFTI1
C***END PROLOGUE VRFFTI
C
C VRFFTPK, VERSION 1, AUGUST 1985
C
DIMENSION wsave(n+15)
C***FIRST EXECUTABLE STATEMENT VRFFTI
IF (n.EQ.1) RETURN
CALL vrfti1(n,wsave(1),wsave(n+1))
RETURN
END subroutine
SUBROUTINE vrfti1(n,wa,fac)
C
C VRFFTPK, VERSION 1, AUGUST 1985
C
USE m_constants, ONLY : pimach
DIMENSION wa(n),fac(15),ntryh(4)
DATA ntryh(1),ntryh(2),ntryh(3),ntryh(4)/4,2,3,5/
nl = n
nf = 0
j = 0
10 j = j + 1
IF ( j <=4 ) THEN
ntry = ntryh(j)
ELSE
ntry = ntry + 2
ENDIF
40 nq = nl/ntry
nr = nl - ntry*nq
IF ( nr /= 0 ) THEN
GOTO 10
ENDIF
nf = nf + 1
fac(nf+2) = ntry
nl = nq
IF (ntry.NE.2) GO TO 70
IF (nf.EQ.1) GO TO 70
DO 60 i = 2,nf
ib = nf - i + 2
fac(ib+2) = fac(ib+1)
60 CONTINUE
fac(3) = 2
70 IF (nl.NE.1) GO TO 40
fac(1) = n
fac(2) = nf
tpi = 2.*pimach()
argh = tpi/real(n)
is = 0
nfm1 = nf - 1
l1 = 1
IF (nfm1.EQ.0) RETURN
DO 100 k1 = 1,nfm1
ip = fac(k1+2)
ld = 0
l2 = l1*ip
ido = n/l2
ipm = ip - 1
DO 90 j = 1,ipm
ld = ld + l1
i = is
argld = real(ld)*argh
fi = 0.
DO 80 ii = 3,ido,2
i = i + 2
fi = fi + 1.
arg = fi*argld
wa(i-1) = cos(arg)
wa(i) = sin(arg)
80 CONTINUE
is = is + ido
90 CONTINUE
l1 = l2
100 CONTINUE
RETURN
END subroutine
SUBROUTINE vrfftf(m,n,r,rt,mdimr,wsave)
C
C***BEGIN PROLOGUE VRFFTF
C***DATE WRITTEN 850801 (YYMMDD)
C***REVISION DATE 900509 (YYMMDD)
C***CATEGORY NO. J1A1
C***KEYWORDS FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM,
C FOURIER ANALYSIS, FORWARD TRANSFORM, MULTIPLE SEQUENCES
C***AUTHOR SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST)
C***PURPOSE Forward real periodic transform, M sequences.
C***DESCRIPTION
C
C Subroutine VRFFTF computes the Fourier coefficients (forward
C transform) of a number of real periodic sequences. Specifically,
C for each sequence the subroutine claculates the independent
C Fourier coefficients described below at output parameter R.
C
C The array WSAVE which is used by subroutine VRFFTF must be
C initialized by calling subroutine VRFFTI(N,WSAVE).
C
C
C Input Parameters
C
C M the number of sequences to be transformed.
C
C N the length of the sequences to be transformed. The method
C is most efficient when N is a product of small primes,
C however n may be any positive integer.
C
C R areal two-dimensional array of size MDIMX x N containing the
C the sequences to be transformed. The sequences are stored
C in the ROWS of R. Thus, the I-th sequence to be transformed,
C X(I,J), J=0,1,...,N-1, is stored as
C
C R(I,J) = X(I,J-1) , J=1, 2, . . . , N.
C
C RT a real two-dimensional work array of size MDIMX x N.
C
C MDIMR the row (or first) dimension of the arrays R and RT exactly
C as they appear in the calling program. This parameter is
C used to specify the variable dimension of these arrays.
C
C WSAVE a real one-dimensional work array which must be dimensioned
C at least N+15. The WSAVE array must be initialized by
C calling subroutine VRFFTI. A different WSAVE array must be
C used for each different value of N. This initialization does
C not have to be repeated so long as N remains unchanged. The
C same WSAVE array may be used by VRFFTF and VRFFTB.
C
C Output Parameters
C
C R contains the Fourier coefficients F(K) for each of the M
C input sequences. Specifically, row I of R, R(I,J),
C J=1,2,..,N, contains the independent Fourier coefficients
C F(I,K), for the I-th input sequence stored as
C
C R(I,1) = REAL( F(I,0) ),
C = SQRT(1/N)*SUM(J=0,N-1)[ X(I,J) ],
C
C R(I,2*K) = REAL( F(I,K) )
C = SQRT(1/N)*SUM(J=0,N-1)[X(I,J)*COS(2J*K*PI/N)]
C
C R(I,2*K+1) = IMAG( F(I,K) )
C =-SQRT(1/N)*SUM(J=0,N-1)[X(I,J)*SIN(2J*K*PI/N)]
C
C for K = 1, 2, . . . , M-1,
C
C and, when N is even,
C
C R(I,N) = REAL( F(I,N/2) ).
C = SQRT(1/N)*SUM(J=0,N-1)[ (-1)**J*X(I,J) ].
C
C WSAVE contains results which must not be destroyed between calls
C to VRFFTF or VRFFTB.
C
C -----------------------------------------------------------------
C
C NOTE - A call of VRFFTF followed immediately by a call of
C of VRFFTB will return the original sequences R. Thus,
C VRFFTB is the correctly normalized inverse of VRFFTF.
C
C -----------------------------------------------------------------
C
C VRFFTF is a straightforward extension of the subprogram RFFTF to
C handle M simultaneous sequences. RFFTF was originally developed
C by P. N. Swarztrauber of NCAR.
C
C
C * * * * * * * * * * * * * * * * * * * * *
C * *
C * PROGRAM SPECIFICATIONS *
C * *
C * * * * * * * * * * * * * * * * * * * * *
C
C
C Dimension of R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15)
C Arguments
C
C Latest AUGUST 1, 1985
C Revision
C
C Subprograms VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3,
C Required VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2,
C VRADB3, VRADB4, VRADB5, VRADBG, PIMACH
C
C Special NONE
C Conditions
C
C Common NONE
C blocks
C
C I/O NONE
C
C Precision SINGLE
C
C Specialist ROLAND SWEET
C
C Language FORTRAN
C
C History WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE
C NATIONAL BUREAU OF STANDARDS (BOULDER).
C
C Algorithm A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION
C OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM.
C
C Portability AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77.
C THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN
C THE FUNCTION PIMACH.
C
C Required COS,SIN
C resident
C Routines
C
C
C***REFERENCES P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
C Computations, (G. Rodrigue, ed.), Academic Press, 1982,
C pp. 51-83.
C***ROUTINES CALLED VRFTF1
C***END PROLOGUE VRFFTF
C
C VRFFTPK, VERSION 1, AUGUST 1985
C
DIMENSION r(mdimr,n),rt(mdimr,n),wsave(n+15)
C***FIRST EXECUTABLE STATEMENT VRFFTF
IF (n.EQ.1) RETURN
CALL vrftf1(m,n,r,rt,mdimr,wsave(1),wsave(n+1))
RETURN
END subroutine
SUBROUTINE vradf2(mp,ido,l1,cc,ch,mdimc,wa1)
C
C VRFFTPK, VERSION 1, AUGUST 1985
C
DIMENSION ch(mdimc,ido,2,l1),cc(mdimc,ido,l1,2),wa1(ido)
DO 20 k = 1,l1
DO 10 m = 1,mp
ch(m,1,1,k) = cc(m,1,k,1) + cc(m,1,k,2)
ch(m,ido,2,k) = cc(m,1,k,1) - cc(m,1,k,2)
10 CONTINUE
20 CONTINUE
IF ( ido == 2 ) THEN
GOTO 70
ELSEIF ( ido < 2 ) THEN
GOTO 100
ENDIF
idp2 = ido + 2
DO 60 k = 1,l1
DO 50 i = 3,ido,2
ic = idp2 - i
DO 40 m = 1,mp
ch(m,i,1,k) = cc(m,i,k,1) +
+ (wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*
+ cc(m,i-1,k,2))
ch(m,ic,2,k) = (wa1(i-2)*cc(m,i,k,2)-
+ wa1(i-1)*cc(m,i-1,k,2)) - cc(m,i,k,1)
ch(m,i-1,1,k) = cc(m,i-1,k,1) +
+ (wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*
+ cc(m,i,k,2))
ch(m,ic-1,2,k) = cc(m,i-1,k,1) -
+ (wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*
+ cc(m,i,k,2))
40 CONTINUE
50 CONTINUE
60 CONTINUE
IF (mod(ido,2).EQ.1) RETURN
70 DO 90 k = 1,l1
DO 80 m = 1,mp
ch(m,1,2,k) = -cc(m,ido,k,2)
ch(m,ido,1,k) = cc(m,ido,k,1)
80 CONTINUE
90 CONTINUE
100 RETURN
END subroutine
SUBROUTINE vradf3(mp,ido,l1,cc,ch,mdimc,wa1,wa2)
C
C VRFFTPK, VERSION 1, AUGUST 1985
C
USE m_constants, ONLY : pimach
DIMENSION ch(mdimc,ido,3,l1),cc(mdimc,ido,l1,3),wa1(ido),wa2(ido)
arg = 2.*pimach()/3.
taur = cos(arg)
taui = sin(arg)
DO 20 k = 1,l1
DO 10 m = 1,mp
ch(m,1,1,k) = cc(m,1,k,1) + (cc(m,1,k,2)+cc(m,1,k,3))
ch(m,1,3,k) = taui* (cc(m,1,k,3)-cc(m,1,k,2))
ch(m,ido,2,k) = cc(m,1,k,1) +
+ taur* (cc(m,1,k,2)+cc(m,1,k,3))
10 CONTINUE
20 CONTINUE
IF (ido.EQ.1) RETURN
idp2 = ido + 2
DO 50 k = 1,l1
DO 40 i = 3,ido,2
ic = idp2 - i
DO 30 m = 1,mp
ch(m,i-1,1,k) = cc(m,i-1,k,1) +
+ ((wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*cc(m,i,
+ k,2))+ (wa2(i-2)*cc(m,i-1,k,
+ 3)+wa2(i-1)*cc(m,i,k,3)))
ch(m,i,1,k) = cc(m,i,k,1) +
+ ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k,
+ 2))+ (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m,
+ i-1,k,3)))
ch(m,i-1,3,k) = (cc(m,i-1,k,1)+
+ taur* ((wa1(i-2)*cc(m,i-1,k,
+ 2)+wa1(i-1)*cc(m,i,k,2))+ (wa2(i-2)*cc(m,
+ i-1,k,3)+wa2(i-1)*cc(m,i,k,3)))) +
+ (taui* ((wa1(i-2)*cc(m,i,k,
+ 2)-wa1(i-1)*cc(m,i-1,k,
+ 2))- (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m,
+ i-1,k,3))))
ch(m,ic-1,2,k) = (cc(m,i-1,k,1)+
+ taur* ((wa1(i-2)*cc(m,i-1,k,
+ 2)+wa1(i-1)*cc(m,i,k,
+ 2))+ (wa2(i-2)*cc(m,i-1,k,
+ 3)+wa2(i-1)*cc(m,i,k,3)))) -
+ (taui* ((wa1(i-2)*cc(m,i,k,
+ 2)-wa1(i-1)*cc(m,i-1,k,
+ 2))- (wa2(i-2)*cc(m,i,k,
+ 3)-wa2(i-1)*cc(m,i-1,k,3))))
ch(m,i,3,k) = (cc(m,i,k,1)+taur*
+ ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k,
+ 2))+ (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m,
+ i-1,k,3)))) + (taui* ((wa2(i-2)*cc(m,i-1,k,
+ 3)+wa2(i-1)*cc(m,i,k,3))- (wa1(i-2)*cc(m,
+ i-1,k,2)+wa1(i-1)*cc(m,i,k,2))))
ch(m,ic,2,k) = (taui* ((wa2(i-2)*cc(m,i-1,k,
+ 3)+wa2(i-1)*cc(m,i,k,3))- (wa1(i-2)*cc(m,
+ i-1,k,2)+wa1(i-1)*cc(m,i,k,2)))) -
+ (cc(m,i,k,1)+taur* ((wa1(i-2)*cc(m,i,k,
+ 2)-wa1(i-1)*cc(m,i-1,k,
+ 2))+ (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m,
+ i-1,k,3))))
30 CONTINUE
40 CONTINUE
50 CONTINUE
RETURN
END subroutine
SUBROUTINE vradf4(mp,ido,l1,cc,ch,mdimc,wa1,wa2,wa3)
C
C VRFFTPK, VERSION 1, AUGUST 1985
C
DIMENSION cc(mdimc,ido,l1,4),ch(mdimc,ido,4,l1),wa1(ido),wa2(ido),
+ wa3(ido)
hsqt2 = sqrt(2.)/2.
DO 20 k = 1,l1
DO 10 m = 1,mp
ch(m,1,1,k) = (cc(m,1,k,2)+cc(m,1,k,4)) +
+ (cc(m,1,k,1)+cc(m,1,k,3))
ch(m,ido,4,k) = (cc(m,1,k,1)+cc(m,1,k,3)) -
+ (cc(m,1,k,2)+cc(m,1,k,4))
ch(m,ido,2,k) = cc(m,1,k,1) - cc(m,1,k,3)
ch(m,1,3,k) = cc(m,1,k,4) - cc(m,1,k,2)
10 CONTINUE
20 CONTINUE
IF ( ido == 2 ) THEN
GOTO 70
ELSEIF ( ido < 2 ) THEN
GOTO 100
ENDIF
idp2 = ido + 2
DO 60 k = 1,l1
DO 50 i = 3,ido,2
ic = idp2 - i
DO 40 m = 1,mp
ch(m,i-1,1,k) = ((wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*cc(m,i,
+ k,2))+ (wa3(i-2)*cc(m,i-1,k,
+ 4)+wa3(i-1)*cc(m,i,k,4))) +
+ (cc(m,i-1,k,1)+ (wa2(i-2)*cc(m,i-1,k,
+ 3)+wa2(i-1)*cc(m,i,k,3)))
ch(m,ic-1,4,k) = (cc(m,i-1,k,1)+
+ (wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*cc(m,i,
+ k,3))) - ((wa1(i-2)*cc(m,i-1,k,
+ 2)+wa1(i-1)*cc(m,i,k,2))+
+ (wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*cc(m,i,
+ k,4)))
ch(m,i,1,k) = ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k,
+ 2))+ (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m,
+ i-1,k,4))) + (cc(m,i,k,1)+
+ (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m,i-1,k,
+ 3)))
ch(m,ic,4,k) = ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,
+ k,2))+ (wa3(i-2)*cc(m,i,k,
+ 4)-wa3(i-1)*cc(m,i-1,k,4))) -
+ (cc(m,i,k,1)+ (wa2(i-2)*cc(m,i,k,
+ 3)-wa2(i-1)*cc(m,i-1,k,3)))
ch(m,i-1,3,k) = ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,
+ k,2))- (wa3(i-2)*cc(m,i,k,
+ 4)-wa3(i-1)*cc(m,i-1,k,4))) +
+ (cc(m,i-1,k,1)- (wa2(i-2)*cc(m,i-1,k,
+ 3)+wa2(i-1)*cc(m,i,k,3)))
ch(m,ic-1,2,k) = (cc(m,i-1,k,1)-
+ (wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*cc(m,i,
+ k,3))) - ((wa1(i-2)*cc(m,i,k,