<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
!DECK FFT551
! SUBROUTINE 'FFT551' - MULTIPLE FAST COSINE TRANSFORM
!
! AUTHOR: CLIVE TEMPERTON, MAY 1988
! [ALL-FORTRAN VERSION: C.T., OCTOBER 1995]
!
! COSINE TRANSFORM OF LENGTH N IS CONVERTED TO
! REAL PERIODIC TRANSFORM OF LENGTH N BY PRE- AND POST-
! PROCESSING. REAL PERIODIC TRANSFORM IS PERFORMED BY
! PRUNING REDUNDANT OPERATIONS FROM COMPLEX TRANSFORM.
!
! SEE FOR EXAMPLE PAUL SWARZTRAUBER, "SYMMETRIC FFT'S",
! MATH. COMP. 47 (1986), 323-346.
!
! A IS THE ARRAY CONTAINING INPUT & OUTPUT DATA
! WORK IS AN AREA OF SIZE (N+1)*MIN(LOT,64)
! TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
! IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N
! INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR'
! (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)
! JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
! N+1 IS THE LENGTH OF THE DATA VECTORS
! (WHICH INCLUDE NONZERO VALUES AT BOTH ENDPOINTS)
! LOT IS THE NUMBER OF DATA VECTORS
! ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
! = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
!
! ORDERING OF COEFFICIENTS: Z(0) , Z(1) , Z(2) , ... , Z(N)
!
! ORDERING OF DATA: X(0) , X(1) , X(2) , ... , X(N)
!
! VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS
! IN PARALLEL
!
! N MUST BE COMPOSED OF FACTORS 2,3 & 5 AND MUST BE EVEN
!
! DEFINITION OF TRANSFORMS:
! -------------------------
!
! ISIGN=+1: X(I)=SUM(J=0,...,N)(E(J)*Z(J)*COS(I*J*PI/N))
! WHERE E(J)=0.5 FOR J=0,N --- ELSE E(J)=1
!
! ISIGN=-1: Z(J)=(2/N)*SUM(I=0,...,N)(E(I)*X(I)*COS(I*J*PI/N))
!
! N.B. FFT551 has an unusual definition of the FFTs,
! such that the the coeff of wave0 is NOT the mean.
!
!---------------------------------------------------------------------
<A NAME='FFT551'><A href='../../html_code/ffts/fft551.inc.html#FFT551' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
Subroutine FFT551 & ! in 8,2
( ISIGN, & ! in
INC, & ! in
JUMP, & ! in
LOT, & ! in
N, & ! in
IFAX, & ! in
TRIGS, & ! in
A, & ! inout
IDIM ) ! in
! Code Description: ORIGINAL CODE F77 IS HARDLY TOUCHED !!!
Integer , intent (in) :: ISIGN ! Switch forward (-1) or inverse (+1)
Integer , intent (in) :: INC ! increment within each data
! vector (e.g. INC=1 for
! consecutively stored data)
Integer , intent (in) :: Jump ! increment between start of
! data vectors
Integer , intent (in) :: LOT ! Number of data vectors
Integer , intent (in) :: N ! N+1 is the length of the data
Integer , intent (in) :: IFAX(10) ! previously prepared list of
! factors of N
Real , intent (in) :: TRIGS(3*N) ! previously prepared list of
! trigonometric function values
Real , intent (inout) :: A( INC*(N+1) + JUMP*(LOT-1) ) ! data array ! vectors (which include zeros
! at the endpoints)
Integer , intent (in) :: IDIM ! dimension workspace
Real :: WORK(IDIM) ! size (n+1)*min(lot,VectorLength)
Integer :: NFAX,NX,NH
Integer :: NBLOX,NVEX,NB
Integer :: K, IC, J, LA, IGO, JA,JB,IA,IB
Integer :: IFAC,IERR,ISTART
Real :: CO,S, t1,t2,si,scale, vectorlength
CHARACTER (LEN=*), PARAMETER :: RoutineName = "Var_FFT551"
VectorLength = LOT
NFAX=IFAX(1)
NX=N+1
NH=N/2
NBLOX=1+(LOT-1)/VectorLength
NVEX=LOT-(NBLOX-1)*VectorLength
ISTART=1
!
DO 200 NB=1,NBLOX
!
! PREPROCESSING
! -------------
IA=ISTART
IB=IA+NH*INC
IC=IA+N*INC
JA=1
JB=NH+1
IF (MOD(NFAX,2).EQ.1) THEN
!DIR$ IVDEP
DO 105 J=1,NVEX
T1=0.5*(A(IA)+A(IC))
T2=0.5*(A(IA)-A(IC))
A(IA)=T1
A(IC)=T2
IA=IA+JUMP
IC=IC+JUMP
105 CONTINUE
ELSE
!DIR$ IVDEP
DO 110 J=1,NVEX
WORK(JA)=0.5*(A(IA)+A(IC))
WORK(JB)=A(IB)
A(IC)=0.5*(A(IA)-A(IC))
IA=IA+JUMP
IB=IB+JUMP
IC=IC+JUMP
JA=JA+NX
JB=JB+NX
110 CONTINUE
ENDIF
!
DO 130 K=1,NH-1
JA=K+1
JB=N+1-K
IA=ISTART+K*INC
IB=ISTART+(JB-1)*INC
IC=ISTART+N*INC
SI=TRIGS(2*N+K)
CO=TRIGS(2*N+NH-K)
IF (MOD(NFAX,2).EQ.1) THEN
!DIR$ IVDEP
DO 115 J=1,NVEX
T1 = 0.5*(A(IA)+A(IB)) - SI*(A(IA)-A(IB))
T2 = 0.5*(A(IA)+A(IB)) + SI*(A(IA)-A(IB))
A(IC) = A(IC) + CO*(A(IA)-A(IB))
A(IA) = T1
A(IB) = T2
IA=IA+JUMP
IB=IB+JUMP
IC=IC+JUMP
115 CONTINUE
ELSE
!DIR$ IVDEP
DO 120 J=1,NVEX
WORK(JA) = 0.5*(A(IA)+A(IB)) - SI*(A(IA)-A(IB))
WORK(JB) = 0.5*(A(IA)+A(IB)) + SI*(A(IA)-A(IB))
A(IC) = A(IC) + CO*(A(IA)-A(IB))
IA=IA+JUMP
IB=IB+JUMP
IC=IC+JUMP
JA=JA+NX
JB=JB+NX
120 CONTINUE
ENDIF
130 CONTINUE
!
! PERIODIC FOURIER ANALYSIS
! -------------------------
IA=ISTART
LA=N
IGO=1-2*MOD(NFAX,2)
!
DO 140 K=1,NFAX
IFAC=IFAX(NFAX+2-K)
LA=LA/IFAC
IERR=-1
IF (IGO.EQ.+1) THEN
CALL qpassm
(WORK,WORK(IFAC*LA+1),A(IA),A(IA+LA*INC), &
TRIGS,1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR)
ELSE IF (IGO.EQ.-1) THEN
CALL qpassm(A(IA),A(IA+IFAC*LA*INC),WORK,WORK(LA+1), &
TRIGS,INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR)
ENDIF
IF (IERR.NE.0) GO TO 500
IGO=-IGO
140 CONTINUE
!
! POSTPROCESSING
! --------------
SCALE=2.0
IF (ISIGN.EQ.+1) SCALE = FLOAT(N)
S=1.0
IF (ISIGN.EQ.-1) S = 2.0/FLOAT(N)
JA=ISTART
JB=JA+N*INC
IA=1
IB=N
!DIR$ IVDEP
DO 150 J=1,NVEX
A(JA)=SCALE*WORK(IA)
A(JA+INC)=S*A(JB)
A(JB)=SCALE*WORK(IB)
IA=IA+NX
IB=IB+NX
JA=JA+JUMP
JB=JB+JUMP
150 CONTINUE
!
DO 170 K=2,N-2,2
JA=ISTART+K*INC
IA=K
!DIR$ IVDEP
DO 160 J=1,NVEX
A(JA)=SCALE*WORK(IA)
A(JA+INC)=-SCALE*WORK(IA+1)+A(JA-INC)
IA=IA+NX
JA=JA+JUMP
160 CONTINUE
170 CONTINUE
!
ISTART=ISTART+NVEX*JUMP
NVEX=VectorLength
200 CONTINUE
GO TO 570
!
! ERROR MESSAGES
! --------------
500 CONTINUE
GO TO (510,530,550) IERR
510 CONTINUE
WRITE(UNIT=0,FMT='(A,I4,A)') &
'VECTOR LENGTH =',NVEX,', GREATER THAN VectorLength'
GO TO 570
530 CONTINUE
WRITE(UNIT=0,FMT='(A,I3,A)') &
'FACTOR =',IFAC,', NOT CATERED FOR'
GO TO 570
550 CONTINUE
WRITE(UNIT=0,FMT='(A,I3,A)') &
'FACTOR =',IFAC,', ONLY CATERED FOR IF LA*IFAC=N'
570 CONTINUE
RETURN
END SUBROUTINE FFT551