Subroutine FFT661 & ! in ( ISIGN, & ! in INC, & ! in JUMP, & ! in LOT, & ! in N, & ! in IFAX, & ! in TRIGS, & ! in A, & ! inout DIM ) ! in ! ! ! Description: ! MULTIPLE FAST SINE TRANSFORM ! (Originally called FFT661, then Var_SinTrans) ! author: clive temperton, may 1988 ! (slightly modified for all-fortran version) ! ! Sine 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. ! ! Method: ! ! 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=1,...,n-1)(z(j)*sin(i*j*pi/n)) ! ! isign=-1: z(j)=(2/n)*sum(i=1,...,n-1)(x(i)*sin(i*j*pi/n)) ! ! Current Code Owner: Andrew Lorenc ! ! History: ! Version Date Comment ! ------- ---- ------- ! 0.1 14/12/93 Original code. Phil Andrews ! 0.2 16/09/94 Small Modifications for the ! incorporation in the VAR project. HB ! 1.1 21/04/95 placed under control. JB ! 1.2 01/06/95 Tracing added. JB ! ! Code Description: ! NB BECAUSE OF THE TRICKY NESTED LOOPS ! ORIGINAL CODE F77 IS HARDLY TOUCHED !!! Implicit none ! Subroutine arguments 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 ! vectors (which include zeros ! at the endpoints) Integer , intent (in) :: DIM ! dimension workspace 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 ! No descriptions given Integer :: NFAX,NX,NH Integer :: NBLOX,NVEX,NB Integer :: K,JA,JB,IA,IB,IGO,LA,J Integer :: IFAC,IERR,ISTART Real :: SI,T1,T2,SCALE, vectorlength Real :: WORK(DIM) ! size (n+1)*min(lot,VectorLength) 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 ! ------------- DO 120 K=1,NH-1 JA=K+1 JB=N+1-K IA=ISTART+K*INC IB=ISTART+(JB-1)*INC SI=TRIGS(2*N+K) IF (MOD(NFAX,2).EQ.0) THEN !DIR$ IVDEP DO 110 J=1,NVEX WORK(JA) = SI*(A(IA)+A(IB)) + 0.5*(A(IA)-A(IB)) WORK(JB) = SI*(A(IA)+A(IB)) - 0.5*(A(IA)-A(IB)) IA=IA+JUMP IB=IB+JUMP JA=JA+NX JB=JB+NX 110 CONTINUE ELSE !DIR$ IVDEP DO 115 J=1,NVEX T1 = SI*(A(IA)+A(IB)) + 0.5*(A(IA)-A(IB)) T2 = SI*(A(IA)+A(IB)) - 0.5*(A(IA)-A(IB)) A(IA) = T1 A(IB) = T2 IA=IA+JUMP IB=IB+JUMP 115 CONTINUE ENDIF 120 CONTINUE JA=1 JB=NH+1 IA=ISTART IB=ISTART+NH*INC IF (MOD(NFAX,2).EQ.0) THEN !DIR$ IVDEP DO 130 J=1,NVEX WORK(JA)=0.0 WORK(JB)=2.0*A(IB) IB=IB+JUMP JA=JA+NX JB=JB+NX 130 CONTINUE IGO = +1 ELSE !DIR$ IVDEP DO 135 J=1,NVEX A(IA)=0.0 A(IB)=2.0*A(IB) IA=IA+JUMP IB=IB+JUMP 135 CONTINUE IGO = -1 ENDIF ! ! PERIODIC FOURIER ANALYSIS ! ------------------------- IA=ISTART LA=N ! 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(LA*INC+IA), & TRIGS,1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR) ELSE IF (IGO.EQ.-1) THEN CALL qpassm(A(IA),A(IFAC*LA*INC+IA),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) JA=ISTART JB=JA+N*INC IA=1 !DIR$ IVDEP DO 150 J=1,NVEX A(JA)=0.0 A(JA+INC)=0.5*SCALE*WORK(IA) A(JB)=0.0 IA=IA+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+1) A(JA+INC)=SCALE*WORK(IA)+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,I5,A)') 'NVEX=', NVEX ,' GREATER THAN VectorLength' GO TO 570 530 CONTINUE WRITE(UNIT=0,FMT='(A,I5,A)') 'IFAC=', IFAC, 'NOT CATERED FOR' GO TO 570 550 CONTINUE WRITE(UNIT=0,FMT='(A,I5,A)') 'IFAC=', IFAC, ' ONLY CATERED FOR IF LA*IFAC=N' 570 CONTINUE RETURN End subroutine FFT661