var/da/da_ffts/fft661.inc

   1 Subroutine FFT661  & ! in
   2  ( ISIGN,              & ! in
   3    INC,                & ! in
   4    JUMP,               & ! in
   5    LOT,                & ! in
   6    N,                  & ! in
   7    IFAX,               & ! in
   8    TRIGS,              & ! in
   9    A,                  & ! inout
  10    DIM )                 ! in
  11 !
  12 !
  13 ! Description:
  14 !     MULTIPLE FAST SINE TRANSFORM
  15 !     (Originally called FFT661, then Var_SinTrans)
  16 !      author: clive temperton, may 1988
  17 !       (slightly modified for all-fortran version)
  18 !
  19 !     Sine transform of length n is converted to
  20 !     Real periodic transform of length n by pre- and post-
  21 !     processing. Real periodic transform is performed by
  22 !     pruning redundant operations from complex transform.
  23 !
  24 !     see for example paul swarztrauber, "symmetric fft's",
  25 !     math. comp. 47 (1986), 323-346.
  26 !
  27 ! Method:
  28 !
  29 !     ordering of coefficients:   z(0) , z(1) , z(2) , ... , z(n)
  30 !     ordering of data:           x(0) , x(1) , x(2) , ... , x(n)
  31 !
  32 !    vectorization is achieved on cray by doing the transforms
  33 !    in parallel
  34 !
  35 !    N must be composed of factors 2,3 & 5 and must be even
  36 !
  37 !    definition of transforms:
  38 !     -------------------------
  39 !
  40 !     isign=+1: x(i)=sum(j=1,...,n-1)(z(j)*sin(i*j*pi/n))
  41 !
  42 !     isign=-1: z(j)=(2/n)*sum(i=1,...,n-1)(x(i)*sin(i*j*pi/n))
  43 !
  44 ! Current Code Owner: Andrew Lorenc
  45 !
  46 !   History:
  47 ! Version   Date     Comment
  48 ! -------   ----     -------
  49 ! 0.1       14/12/93 Original code. Phil Andrews
  50 ! 0.2       16/09/94 Small Modifications for the
  51 !                    incorporation in the VAR project. HB
  52 ! 1.1       21/04/95 placed under control. JB
  53 ! 1.2       01/06/95 Tracing added. JB
  54 !
  55 ! Code Description:
  56 !    NB   BECAUSE OF THE TRICKY NESTED LOOPS
  57 !         ORIGINAL CODE F77 IS HARDLY TOUCHED !!!
  58
  59 Implicit none
  60
  61 ! Subroutine arguments
  62  Integer , intent (in)    :: ISIGN         ! Switch forward (-1) or inverse (+1)
  63  Integer , intent (in)    :: INC           ! increment within each data
  64                                            ! vector  (e.g. INC=1 for
  65                                            ! consecutively stored data)
  66  Integer , intent (in)    :: Jump          ! increment between start of
  67                                            ! data vectors
  68  Integer , intent (in)    :: LOT           ! Number of data vectors
  69  Integer , intent (in)    :: N             ! N+1 is the length of the data
  70                                            !  vectors  (which include zeros
  71                                            ! at the endpoints)
  72  Integer , intent (in)    :: DIM           ! dimension workspace
  73  Integer , intent (in)    :: IFAX(10)      ! previously prepared list of
  74                                            ! factors of N
  75
  76  Real    , intent (in)    :: TRIGS(3*N)    ! previously prepared list of
  77                                            ! trigonometric function values
  78  Real    , intent (inout) :: A( INC*(N+1) + JUMP*(LOT-1) ) ! data array
  79
  80                                                     ! No descriptions given
  81  Integer                  :: NFAX,NX,NH
  82  Integer                  :: NBLOX,NVEX,NB
  83  Integer                  :: K,JA,JB,IA,IB,IGO,LA,J
  84  Integer                  :: IFAC,IERR,ISTART
  85
  86  Real                     :: SI,T1,T2,SCALE, vectorlength
  87  Real                     :: WORK(DIM)     ! size (n+1)*min(lot,VectorLength)
  88
  89       VectorLength = LOT
  90       NFAX=IFAX(1)
  91       NX=N+1
  92       NH=N/2
  93       NBLOX=1+(LOT-1)/VectorLength
  94       NVEX=LOT-(NBLOX-1)*VectorLength
  95       ISTART=1
  96 !
  97       DO 200 NB=1,NBLOX
  98 !
  99 !     PREPROCESSING
 100 !     -------------
 101       DO 120 K=1,NH-1
 102       JA=K+1
 103       JB=N+1-K
 104       IA=ISTART+K*INC
 105       IB=ISTART+(JB-1)*INC
 106       SI=TRIGS(2*N+K)
 107       IF (MOD(NFAX,2).EQ.0) THEN
 108 !DIR$ IVDEP
 109          DO 110 J=1,NVEX
 110          WORK(JA) = SI*(A(IA)+A(IB)) + 0.5*(A(IA)-A(IB))
 111          WORK(JB) = SI*(A(IA)+A(IB)) - 0.5*(A(IA)-A(IB))
 112          IA=IA+JUMP
 113          IB=IB+JUMP
 114          JA=JA+NX
 115          JB=JB+NX
 116   110    CONTINUE
 117       ELSE
 118 !DIR$ IVDEP
 119          DO 115 J=1,NVEX
 120          T1 = SI*(A(IA)+A(IB)) + 0.5*(A(IA)-A(IB))
 121          T2 = SI*(A(IA)+A(IB)) - 0.5*(A(IA)-A(IB))
 122          A(IA) = T1
 123          A(IB) = T2
 124          IA=IA+JUMP
 125          IB=IB+JUMP
 126   115    CONTINUE
 127       ENDIF
 128   120 CONTINUE
 129
 130       JA=1
 131       JB=NH+1
 132       IA=ISTART
 133       IB=ISTART+NH*INC
 134       IF (MOD(NFAX,2).EQ.0) THEN
 135 !DIR$ IVDEP
 136          DO 130 J=1,NVEX
 137          WORK(JA)=0.0
 138          WORK(JB)=2.0*A(IB)
 139          IB=IB+JUMP
 140          JA=JA+NX
 141          JB=JB+NX
 142   130    CONTINUE
 143          IGO = +1
 144       ELSE
 145 !DIR$ IVDEP
 146          DO 135 J=1,NVEX
 147          A(IA)=0.0
 148          A(IB)=2.0*A(IB)
 149          IA=IA+JUMP
 150          IB=IB+JUMP
 151   135    CONTINUE
 152          IGO = -1
 153       ENDIF
 154 !
 155 !     PERIODIC FOURIER ANALYSIS
 156 !     -------------------------
 157       IA=ISTART
 158       LA=N
 159 !
 160       DO 140 K=1,NFAX
 161       IFAC=IFAX(NFAX+2-K)
 162       LA=LA/IFAC
 163       IERR=-1
 164       IF (IGO.EQ.+1) THEN
 165         CALL qpassm(WORK,WORK(IFAC*LA+1),A(IA),A(LA*INC+IA), &
 166                     TRIGS,1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR)
 167       ELSE IF (IGO.EQ.-1) THEN
 168         CALL qpassm(A(IA),A(IFAC*LA*INC+IA),WORK,WORK(LA+1), &
 169                     TRIGS,INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR)
 170       ENDIF
 171       IF (IERR.NE.0) GO TO 500
 172       IGO=-IGO
 173   140 CONTINUE
 174 !
 175 !     POSTPROCESSING
 176 !     --------------
 177       SCALE=2.0
 178       IF (ISIGN.EQ.+1) SCALE = FLOAT(N)
 179       JA=ISTART
 180       JB=JA+N*INC
 181       IA=1
 182 !DIR$ IVDEP
 183       DO 150 J=1,NVEX
 184       A(JA)=0.0
 185       A(JA+INC)=0.5*SCALE*WORK(IA)
 186       A(JB)=0.0
 187       IA=IA+NX
 188       JA=JA+JUMP
 189       JB=JB+JUMP
 190   150 CONTINUE
 191 !
 192       DO 170 K=2,N-2,2
 193       JA=ISTART+K*INC
 194       IA=K
 195 !DIR$ IVDEP
 196       DO 160 J=1,NVEX
 197       A(JA)=-SCALE*WORK(IA+1)
 198       A(JA+INC)=SCALE*WORK(IA)+A(JA-INC)
 199       IA=IA+NX
 200       JA=JA+JUMP
 201   160 CONTINUE
 202   170 CONTINUE
 203 !
 204       ISTART=ISTART+NVEX*JUMP
 205       NVEX=VectorLength
 206   200 CONTINUE
 207
 208       Go To 570
 209 !
 210 !     ERROR MESSAGES
 211 !     --------------
 212   500 CONTINUE
 213       GO TO (510,530,550) IERR
 214   510 CONTINUE
 215       WRITE(UNIT=0,FMT='(A,I5,A)') 'NVEX=', NVEX ,' GREATER THAN VectorLength'
 216       GO TO 570
 217   530 CONTINUE
 218       WRITE(UNIT=0,FMT='(A,I5,A)') 'IFAC=', IFAC, 'NOT CATERED FOR'
 219       GO TO 570
 220   550 CONTINUE
 221       WRITE(UNIT=0,FMT='(A,I5,A)') 'IFAC=', IFAC, ' ONLY CATERED FOR IF LA*IFAC=N'
 222   570 CONTINUE
 223
 224
 225       RETURN
 226
 227
 228 End subroutine FFT661
 229
 230
 231