C*********************************************************************
C                                                                    *
C     GPFAPACK - FORTRAN IMPLEMENTATION OF THE SELF-SORTING          *
C     IN-PLACE GENERALIZED PRIME FACTOR (COMPLEX) FFT [GPFA]         *
C                                                                    *
C     WRITTEN BY CLIVE TEMPERTON                                     *
C     RECHERCHE EN PREVISION NUMERIQUE / ECMWF                       *
C                                                                    *
C     THE PACKAGE CONSISTS OF THE SETUP ROUTINE SETGPFA, TOGETHER    *
C     WITH THE ROUTINES GPFA, GPFA2F, GPFA3F, GPFA5F                 *
C                                                                    *
C*********************************************************************
C
C        SUBROUTINE 'SETGPFA'
C        SETUP ROUTINE FOR SELF-SORTING IN-PLACE
C            GENERALIZED PRIME FACTOR (COMPLEX) FFT [GPFA]
C
C        CALL SETGPFA(TRIGS,N)
C
C        INPUT :
C        -----
C        N IS THE LENGTH OF THE TRANSFORMS. N MUST BE OF THE FORM:
C          -----------------------------------
C            N = (2**IP) * (3**IQ) * (5**IR)
C          -----------------------------------
C
C        OUTPUT:
C        ------
C        TRIGS IS A TABLE OF TWIDDLE FACTORS,
C          OF LENGTH 2*IPQR (REAL) WORDS, WHERE:
C          --------------------------------------
C            IPQR = (2**IP) + (3**IQ) + (5**IR)
C          --------------------------------------
C
C        WRITTEN BY CLIVE TEMPERTON 1990
C
C----------------------------------------------------------------------
C
      SUBROUTINE SETGPFA(TRIGS,N)
C
      DIMENSION TRIGS(*)
      DIMENSION NJ(3)
C
C     DECOMPOSE N INTO FACTORS 2,3,5
C     ------------------------------
      NN = N
      IFAC = 2
C
      DO 30 LL = 1 , 3
      KK = 0
   10 CONTINUE
      IF (MOD(NN,IFAC).NE.0) GO TO 20
      KK = KK + 1
      NN = NN / IFAC
      GO TO 10
   20 CONTINUE
      NJ(LL) = KK
      IFAC = IFAC + LL
   30 CONTINUE
C
      IF (NN.NE.1) THEN
         WRITE(6,40) N
   40    FORMAT(' *** WARNING!!!',I10,' IS NOT A LEGAL VALUE OF N ***')
         RETURN
      ENDIF
C
      IP = NJ(1)
      IQ = NJ(2)
      IR = NJ(3)
C
C     COMPUTE LIST OF ROTATED TWIDDLE FACTORS
C     ---------------------------------------
      NJ(1) = 2**IP
      NJ(2) = 3**IQ
      NJ(3) = 5**IR
C
      TWOPI = 4.0 * ASIN(1.0)
      I = 1
C
      DO 60 LL = 1 , 3
      NI = NJ(LL)
      IF (NI.EQ.1) GO TO 60
C
      DEL = TWOPI / FLOAT(NI)
      IROT = N / NI
      KINK = MOD(IROT,NI)
      KK = 0
C
      DO 50 K = 1 , NI
      ANGLE = FLOAT(KK) * DEL
      TRIGS(I) = COS(ANGLE)
      TRIGS(I+1) = SIN(ANGLE)
      I = I + 2
      KK = KK + KINK
      IF (KK.GT.NI) KK = KK - NI
   50 CONTINUE
   60 CONTINUE
C
      RETURN
      END
C        SUBROUTINE 'GPFA'
C        SELF-SORTING IN-PLACE GENERALIZED PRIME FACTOR (COMPLEX) FFT
C
C        *** THIS IS THE ALL-FORTRAN VERSION ***
C            -------------------------------
C
C        CALL GPFA(A,B,TRIGS,INC,JUMP,N,LOT,ISIGN)
C
C        A IS FIRST REAL INPUT/OUTPUT VECTOR
C        B IS FIRST IMAGINARY INPUT/OUTPUT VECTOR
C        TRIGS IS A TABLE OF TWIDDLE FACTORS, PRECALCULATED
C              BY CALLING SUBROUTINE 'SETGPFA'
C        INC IS THE INCREMENT WITHIN EACH DATA VECTOR
C        JUMP IS THE INCREMENT BETWEEN DATA VECTORS
C        N IS THE LENGTH OF THE TRANSFORMS:
C          -----------------------------------
C            N = (2**IP) * (3**IQ) * (5**IR)
C          -----------------------------------
C        LOT IS THE NUMBER OF TRANSFORMS
C        ISIGN = +1 FOR FORWARD TRANSFORM
C              = -1 FOR INVERSE TRANSFORM
C
C        WRITTEN BY CLIVE TEMPERTON
C        RECHERCHE EN PREVISION NUMERIQUE
C        ATMOSPHERIC ENVIRONMENT SERVICE, CANADA
C
C----------------------------------------------------------------------
C
C        DEFINITION OF TRANSFORM
C        -----------------------
C
C        X(J) = SUM(K=0,...,N-1)(C(K)*EXP(ISIGN*2*I*J*K*PI/N))
C
C---------------------------------------------------------------------
C
C        FOR A MATHEMATICAL DEVELOPMENT OF THE ALGORITHM USED,
C        SEE:
C
C        C TEMPERTON : "A GENERALIZED PRIME FACTOR FFT ALGORITHM
C          FOR ANY N = (2**P)(3**Q)(5**R)",
C          SIAM J. SCI. STAT. COMP., MAY 1992.
C
C----------------------------------------------------------------------
C
      SUBROUTINE GPFA(A,B,TRIGS,INC,JUMP,N,LOT,ISIGN)
C
      DIMENSION A(*), B(*), TRIGS(*)
      DIMENSION NJ(3)
C
C     DECOMPOSE N INTO FACTORS 2,3,5
C     ------------------------------
      NN = N
      IFAC = 2
C
      DO 30 LL = 1 , 3
      KK = 0
   10 CONTINUE
      IF (MOD(NN,IFAC).NE.0) GO TO 20
      KK = KK + 1
      NN = NN / IFAC
      GO TO 10
   20 CONTINUE
      NJ(LL) = KK
      IFAC = IFAC + LL
   30 CONTINUE
C
      IF (NN.NE.1) THEN
         WRITE(6,40) N
   40    FORMAT(' *** WARNING!!!',I10,' IS NOT A LEGAL VALUE OF N ***')
         RETURN
      ENDIF
C
      IP = NJ(1)
      IQ = NJ(2)
      IR = NJ(3)
C
C     COMPUTE THE TRANSFORM
C     ---------------------
      I = 1
      IF (IP.GT.0) THEN
         CALL GPFA2F(A,B,TRIGS,INC,JUMP,N,IP,LOT,ISIGN)
         I = I + 2 * ( 2**IP)
      ENDIF
      IF (IQ.GT.0) THEN
         CALL GPFA3F(A,B,TRIGS(I),INC,JUMP,N,IQ,LOT,ISIGN)
         I = I + 2 * (3**IQ)
      ENDIF
      IF (IR.GT.0) THEN
         CALL GPFA5F(A,B,TRIGS(I),INC,JUMP,N,IR,LOT,ISIGN)
      ENDIF
C
      RETURN
      END
C     fortran version of *gpfa2* -
C     radix-2 section of self-sorting, in-place, generalized pfa
C     central radix-2 and radix-8 passes included
C      so that transform length can be any power of 2
C
C-------------------------------------------------------------------
C
      subroutine gpfa2f(a,b,trigs,inc,jump,n,mm,lot,isign)
      dimension a(*), b(*), trigs(*)
      data lvr/64/
C
C     ***************************************************************
C     *                                                             *
C     *  N.B. LVR = LENGTH OF VECTOR REGISTERS, SET TO 128 FOR C90. *
C     *  RESET TO 64 FOR OTHER CRAY MACHINES, OR TO ANY LARGE VALUE *
C     *  (GREATER THAN OR EQUAL TO LOT) FOR A SCALAR COMPUTER.      *
C     *                                                             *
C     ***************************************************************
C
      n2 = 2**mm
      inq = n/n2
      jstepx = (n2-n) * inc
      ninc = n * inc
      ink = inc * inq
C
      m2 = 0
      m8 = 0
      if (mod(mm,2).eq.0) then
         m = mm/2
      else if (mod(mm,4).eq.1) then
         m = (mm-1)/2
         m2 = 1
      else if (mod(mm,4).eq.3) then
         m = (mm-3)/2
         m8 = 1
      endif
      mh = (m+1)/2
C
      nblox = 1 + (lot-1)/lvr
      left = lot
      s = float(isign)
      istart = 1
C
C  loop on blocks of lvr transforms
C  --------------------------------
      do 500 nb = 1 , nblox
C
      if (left.le.lvr) then
         nvex = left
      else if (left.lt.(2*lvr)) then
         nvex = left/2
         nvex = nvex + mod(nvex,2)
      else
         nvex = lvr
      endif
      left = left - nvex
C
      la = 1
C
C  loop on type I radix-4 passes
C  -----------------------------
      mu = mod(inq,4)
      if (isign.eq.-1) mu = 4 - mu
      ss = 1.0
      if (mu.eq.3) ss = -1.0
C
      if (mh.eq.0) go to 200
C
      do 160 ipass = 1 , mh
      jstep = (n*inc) / (4*la)
      jstepl = jstep - ninc
C
C  k = 0 loop (no twiddle factors)
C  -------------------------------
      do 120 jjj = 0 , (n-1)*inc , 4*jstep
      ja = istart + jjj
C
C     "transverse" loop
C     -----------------
      do 115 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      jc = jb + jstepl
      if (jc.lt.istart) jc = jc + ninc
      jd = jc + jstepl
      if (jd.lt.istart) jd = jd + ninc
      j = 0
C
C  loop across transforms
C  ----------------------
cdir$ ivdep, shortloop
      do 110 l = 1 , nvex
      t0 = a(ja+j) + a(jc+j)
      t2 = a(ja+j) - a(jc+j)
      t1 = a(jb+j) + a(jd+j)
      t3 = ss * ( a(jb+j) - a(jd+j) )
      u0 = b(ja+j) + b(jc+j)
      u2 = b(ja+j) - b(jc+j)
      u1 = b(jb+j) + b(jd+j)
      u3 = ss * ( b(jb+j) - b(jd+j) )
      a(ja+j) = t0 + t1
      a(jc+j) = t0 - t1
      b(ja+j) = u0 + u1
      b(jc+j) = u0 - u1
      a(jb+j) = t2 - u3
      a(jd+j) = t2 + u3
      b(jb+j) = u2 + t3
      b(jd+j) = u2 - t3
      j = j + jump
  110 continue
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  115 continue
  120 continue
C
C  finished if n2 = 4
C  ------------------
      if (n2.eq.4) go to 490
      kk = 2 * la
C
C  loop on nonzero k
C  -----------------
      do 150 k = ink , jstep-ink , ink
      co1 = trigs(kk+1)
      si1 = s*trigs(kk+2)
      co2 = trigs(2*kk+1)
      si2 = s*trigs(2*kk+2)
      co3 = trigs(3*kk+1)
      si3 = s*trigs(3*kk+2)
C
C  loop along transform
C  --------------------
      do 140 jjj = k , (n-1)*inc , 4*jstep
      ja = istart + jjj
C
C     "transverse" loop
C     -----------------
      do 135 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      jc = jb + jstepl
      if (jc.lt.istart) jc = jc + ninc
      jd = jc + jstepl
      if (jd.lt.istart) jd = jd + ninc
      j = 0
C
C  loop across transforms
C  ----------------------
cdir$ ivdep,shortloop
      do 130 l = 1 , nvex
      t0 = a(ja+j) + a(jc+j)
      t2 = a(ja+j) - a(jc+j)
      t1 = a(jb+j) + a(jd+j)
      t3 = ss * ( a(jb+j) - a(jd+j ) )
      u0 = b(ja+j) + b(jc+j)
      u2 = b(ja+j) - b(jc+j)
      u1 = b(jb+j) + b(jd+j)
      u3 = ss * ( b(jb+j) - b(jd+j) )
      a(ja+j) = t0 + t1
      b(ja+j) = u0 + u1
      a(jb+j) = co1*(t2-u3) - si1*(u2+t3)
      b(jb+j) = si1*(t2-u3) + co1*(u2+t3)
      a(jc+j) = co2*(t0-t1) - si2*(u0-u1)
      b(jc+j) = si2*(t0-t1) + co2*(u0-u1)
      a(jd+j) = co3*(t2+u3) - si3*(u2-t3)
      b(jd+j) = si3*(t2+u3) + co3*(u2-t3)
      j = j + jump
  130 continue
C-----( end of loop across transforms )
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  135 continue
  140 continue
C-----( end of loop along transforms )
      kk = kk + 2*la
  150 continue
C-----( end of loop on nonzero k )
      la = 4*la
  160 continue
C-----( end of loop on type I radix-4 passes)
C
C  central radix-2 pass
C  --------------------
  200 continue
      if (m2.eq.0) go to 300
C
      jstep = (n*inc) / (2*la)
      jstepl = jstep - ninc
C
C  k=0 loop (no twiddle factors)
C  -----------------------------
      do 220 jjj = 0 , (n-1)*inc , 2*jstep
      ja = istart + jjj
C
C     "transverse" loop
C     -----------------
      do 215 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      j = 0
C
C  loop across transforms
C  ----------------------
cdir$ ivdep, shortloop
      do 210 l = 1 , nvex
      t0 = a(ja+j) - a(jb+j)
      a(ja+j) = a(ja+j) + a(jb+j)
      a(jb+j) = t0
      u0 = b(ja+j) - b(jb+j)
      b(ja+j) = b(ja+j) + b(jb+j)
      b(jb+j) = u0
      j = j + jump
  210 continue
C-----(end of loop across transforms)
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  215 continue
  220 continue
C
C  finished if n2=2
C  ----------------
      if (n2.eq.2) go to 490
C
      kk = 2 * la
C
C  loop on nonzero k
C  -----------------
      do 260 k = ink , jstep - ink , ink
      co1 = trigs(kk+1)
      si1 = s*trigs(kk+2)
C
C  loop along transforms
C  ---------------------
      do 250 jjj = k , (n-1)*inc , 2*jstep
      ja = istart + jjj
C
C     "transverse" loop
C     -----------------
      do 245 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      j = 0
C
C  loop across transforms
C  ----------------------
      if (kk.eq.n2/2) then
cdir$ ivdep, shortloop
      do 230 l = 1 , nvex
      t0 = ss * ( a(ja+j) - a(jb+j) )
      a(ja+j) = a(ja+j) + a(jb+j)
      a(jb+j) = ss * ( b(jb+j) - b(ja+j) )
      b(ja+j) = b(ja+j) + b(jb+j)
      b(jb+j) = t0
      j = j + jump
  230 continue
C
      else
C
cdir$ ivdep, shortloop
      do 240 l = 1 , nvex
      t0 = a(ja+j) - a(jb+j)
      a(ja+j) = a(ja+j) + a(jb+j)
      u0 = b(ja+j) - b(jb+j)
      b(ja+j) = b(ja+j) + b(jb+j)
      a(jb+j) = co1*t0 - si1*u0
      b(jb+j) = si1*t0 + co1*u0
      j = j + jump
  240 continue
C
      endif
C
C-----(end of loop across transforms)
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  245 continue
  250 continue
C-----(end of loop along transforms)
      kk = kk + 2 * la
  260 continue
C-----(end of loop on nonzero k)
C-----(end of radix-2 pass)
C
      la = 2 * la
      go to 400
C
C  central radix-8 pass
C  --------------------
  300 continue
      if (m8.eq.0) go to 400
      jstep = (n*inc) / (8*la)
      jstepl = jstep - ninc
      mu = mod(inq,8)
      if (isign.eq.-1) mu = 8 - mu
      c1 = 1.0
      if (mu.eq.3.or.mu.eq.7) c1 = -1.0
      c2 = sqrt(0.5)
      if (mu.eq.3.or.mu.eq.5) c2 = -c2
      c3 = c1 * c2
C
C  stage 1
C  -------
      do 320 k = 0 , jstep - ink , ink
      do 315 jjj = k , (n-1)*inc , 8*jstep
      ja = istart + jjj
C
C     "transverse" loop
C     -----------------
      do 312 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      jc = jb + jstepl
      if (jc.lt.istart) jc = jc + ninc
      jd = jc + jstepl
      if (jd.lt.istart) jd = jd + ninc
      je = jd + jstepl
      if (je.lt.istart) je = je + ninc
      jf = je + jstepl
      if (jf.lt.istart) jf = jf + ninc
      jg = jf + jstepl
      if (jg.lt.istart) jg = jg + ninc
      jh = jg + jstepl
      if (jh.lt.istart) jh = jh + ninc
      j = 0
cdir$ ivdep, shortloop
      do 310 l = 1 , nvex
      t0 = a(ja+j) - a(je+j)
      a(ja+j) = a(ja+j) + a(je+j)
      t1 = c1 * ( a(jc+j) - a(jg+j) )
      a(je+j) = a(jc+j) + a(jg+j)
      t2 = a(jb+j) - a(jf+j)
      a(jc+j) = a(jb+j) + a(jf+j)
      t3 = a(jd+j) - a(jh+j)
      a(jg+j) = a(jd+j) + a(jh+j)
      a(jb+j) = t0
      a(jf+j) = t1
      a(jd+j) = c2 * ( t2 - t3 )
      a(jh+j) = c3 * ( t2 + t3 )
      u0 = b(ja+j) - b(je+j)
      b(ja+j) = b(ja+j) + b(je+j)
      u1 = c1 * ( b(jc+j) - b(jg+j) )
      b(je+j) = b(jc+j) + b(jg+j)
      u2 = b(jb+j) - b(jf+j)
      b(jc+j) = b(jb+j) + b(jf+j)
      u3 = b(jd+j) - b(jh+j)
      b(jg+j) = b(jd+j) + b(jh+j)
      b(jb+j) = u0
      b(jf+j) = u1
      b(jd+j) = c2 * ( u2 - u3 )
      b(jh+j) = c3 * ( u2 + u3 )
      j = j + jump
  310 continue
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  312 continue
  315 continue
  320 continue
C
C  stage 2
C  -------
C
C  k=0 (no twiddle factors)
C  ------------------------
      do 330 jjj = 0 , (n-1)*inc , 8*jstep
      ja = istart + jjj
C
C     "transverse" loop
C     -----------------
      do 328 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      jc = jb + jstepl
      if (jc.lt.istart) jc = jc + ninc
      jd = jc + jstepl
      if (jd.lt.istart) jd = jd + ninc
      je = jd + jstepl
      if (je.lt.istart) je = je + ninc
      jf = je + jstepl
      if (jf.lt.istart) jf = jf + ninc
      jg = jf + jstepl
      if (jg.lt.istart) jg = jg + ninc
      jh = jg + jstepl
      if (jh.lt.istart) jh = jh + ninc
      j = 0
cdir$ ivdep, shortloop
      do 325 l = 1 , nvex
      t0 = a(ja+j) + a(je+j)
      t2 = a(ja+j) - a(je+j)
      t1 = a(jc+j) + a(jg+j)
      t3 = c1 * ( a(jc+j) - a(jg+j) )
      u0 = b(ja+j) + b(je+j)
      u2 = b(ja+j) - b(je+j)
      u1 = b(jc+j) + b(jg+j)
      u3 = c1 * ( b(jc+j) - b(jg+j ) )
      a(ja+j) = t0 + t1
      a(je+j) = t0 - t1
      b(ja+j) = u0 + u1
      b(je+j) = u0 - u1
      a(jc+j) = t2 - u3
      a(jg+j) = t2 + u3
      b(jc+j) = u2 + t3
      b(jg+j) = u2 - t3
      t0 = a(jb+j) + a(jd+j)
      t2 = a(jb+j) - a(jd+j)
      t1 = a(jf+j) - a(jh+j)
      t3 = a(jf+j) + a(jh+j)
      u0 = b(jb+j) + b(jd+j)
      u2 = b(jb+j) - b(jd+j)
      u1 = b(jf+j) - b(jh+j)
      u3 = b(jf+j) + b(jh+j)
      a(jb+j) = t0 - u3
      a(jh+j) = t0 + u3
      b(jb+j) = u0 + t3
      b(jh+j) = u0 - t3
      a(jd+j) = t2 + u1
      a(jf+j) = t2 - u1
      b(jd+j) = u2 - t1
      b(jf+j) = u2 + t1
      j = j + jump
  325 continue
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  328 continue
  330 continue
C
      if (n2.eq.8) go to 490
C
C  loop on nonzero k
C  -----------------
      kk = 2 * la
C
      do 350 k = ink , jstep - ink , ink
C
      co1 = trigs(kk+1)
      si1 = s * trigs(kk+2)
      co2 = trigs(2*kk+1)
      si2 = s * trigs(2*kk+2)
      co3 = trigs(3*kk+1)
      si3 = s * trigs(3*kk+2)
      co4 = trigs(4*kk+1)
      si4 = s * trigs(4*kk+2)
      co5 = trigs(5*kk+1)
      si5 = s * trigs(5*kk+2)
      co6 = trigs(6*kk+1)
      si6 = s * trigs(6*kk+2)
      co7 = trigs(7*kk+1)
      si7 = s * trigs(7*kk+2)
C
      do 345 jjj = k , (n-1)*inc , 8*jstep
      ja = istart + jjj
C
C     "transverse" loop
C     -----------------
      do 342 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      jc = jb + jstepl
      if (jc.lt.istart) jc = jc + ninc
      jd = jc + jstepl
      if (jd.lt.istart) jd = jd + ninc
      je = jd + jstepl
      if (je.lt.istart) je = je + ninc
      jf = je + jstepl
      if (jf.lt.istart) jf = jf + ninc
      jg = jf + jstepl
      if (jg.lt.istart) jg = jg + ninc
      jh = jg + jstepl
      if (jh.lt.istart) jh = jh + ninc
      j = 0
cdir$ ivdep, shortloop
      do 340 l = 1 , nvex
      t0 = a(ja+j) + a(je+j)
      t2 = a(ja+j) - a(je+j)
      t1 = a(jc+j) + a(jg+j)
      t3 = c1 * ( a(jc+j) - a(jg+j) )
      u0 = b(ja+j) + b(je+j)
      u2 = b(ja+j) - b(je+j)
      u1 = b(jc+j) + b(jg+j)
      u3 = c1 * ( b(jc+j) - b(jg+j ) )
      a(ja+j) = t0 + t1
      b(ja+j) = u0 + u1
      a(je+j) = co4*(t0-t1) - si4*(u0-u1)
      b(je+j) = si4*(t0-t1) + co4*(u0-u1)
      a(jc+j) = co2*(t2-u3) - si2*(u2+t3)
      b(jc+j) = si2*(t2-u3) + co2*(u2+t3)
      a(jg+j) = co6*(t2+u3) - si6*(u2-t3)
      b(jg+j) = si6*(t2+u3) + co6*(u2-t3)
      t0 = a(jb+j) + a(jd+j)
      t2 = a(jb+j) - a(jd+j)
      t1 = a(jf+j) - a(jh+j)
      t3 = a(jf+j) + a(jh+j)
      u0 = b(jb+j) + b(jd+j)
      u2 = b(jb+j) - b(jd+j)
      u1 = b(jf+j) - b(jh+j)
      u3 = b(jf+j) + b(jh+j)
      a(jb+j) = co1*(t0-u3) - si1*(u0+t3)
      b(jb+j) = si1*(t0-u3) + co1*(u0+t3)
      a(jh+j) = co7*(t0+u3) - si7*(u0-t3)
      b(jh+j) = si7*(t0+u3) + co7*(u0-t3)
      a(jd+j) = co3*(t2+u1) - si3*(u2-t1)
      b(jd+j) = si3*(t2+u1) + co3*(u2-t1)
      a(jf+j) = co5*(t2-u1) - si5*(u2+t1)
      b(jf+j) = si5*(t2-u1) + co5*(u2+t1)
      j = j + jump
  340 continue
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  342 continue
  345 continue
      kk = kk + 2 * la
  350 continue
C
      la = 8 * la
C
C  loop on type II radix-4 passes
C  ------------------------------
  400 continue
      mu = mod(inq,4)
      if (isign.eq.-1) mu = 4 - mu
      ss = 1.0
      if (mu.eq.3) ss = -1.0
C
      do 480 ipass = mh+1 , m
      jstep = (n*inc) / (4*la)
      jstepl = jstep - ninc
      laincl = la * ink - ninc
C
C  k=0 loop (no twiddle factors)
C  -----------------------------
      do 430 ll = 0 , (la-1)*ink , 4*jstep
C
      do 420 jjj = ll , (n-1)*inc , 4*la*ink
      ja = istart + jjj
C
C     "transverse" loop
C     -----------------
      do 415 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      jc = jb + jstepl
      if (jc.lt.istart) jc = jc + ninc
      jd = jc + jstepl
      if (jd.lt.istart) jd = jd + ninc
      je = ja + laincl
      if (je.lt.istart) je = je + ninc
      jf = je + jstepl
      if (jf.lt.istart) jf = jf + ninc
      jg = jf + jstepl
      if (jg.lt.istart) jg = jg + ninc
      jh = jg + jstepl
      if (jh.lt.istart) jh = jh + ninc
      ji = je + laincl
      if (ji.lt.istart) ji = ji + ninc
      jj = ji + jstepl
      if (jj.lt.istart) jj = jj + ninc
      jk = jj + jstepl
      if (jk.lt.istart) jk = jk + ninc
      jl = jk + jstepl
      if (jl.lt.istart) jl = jl + ninc
      jm = ji + laincl
      if (jm.lt.istart) jm = jm + ninc
      jn = jm + jstepl
      if (jn.lt.istart) jn = jn + ninc
      jo = jn + jstepl
      if (jo.lt.istart) jo = jo + ninc
      jp = jo + jstepl
      if (jp.lt.istart) jp = jp + ninc
      j = 0
C
C  loop across transforms
C  ----------------------
cdir$ ivdep, shortloop
      do 410 l = 1 , nvex
      t0 = a(ja+j) + a(jc+j)
      t2 = a(ja+j) - a(jc+j)
      t1 = a(jb+j) + a(jd+j)
      t3 = ss * ( a(jb+j) - a(jd+j) )
      a(jc+j) = a(ji+j)
      u0 = b(ja+j) + b(jc+j)
      u2 = b(ja+j) - b(jc+j)
      u1 = b(jb+j) + b(jd+j)
      u3 = ss * ( b(jb+j) - b(jd+j) )
      a(jb+j) = a(je+j)
      a(ja+j) = t0 + t1
      a(ji+j) = t0 - t1
      b(ja+j) = u0 + u1
      b(jc+j) = u0 - u1
      b(jd+j) = b(jm+j)
      a(je+j) = t2 - u3
      a(jd+j) = t2 + u3
      b(jb+j) = u2 + t3
      b(jm+j) = u2 - t3
C----------------------
      t0 = a(jb+j) + a(jg+j)
      t2 = a(jb+j) - a(jg+j)
      t1 = a(jf+j) + a(jh+j)
      t3 = ss * ( a(jf+j) - a(jh+j) )
      a(jg+j) = a(jj+j)
      u0 = b(je+j) + b(jg+j)
      u2 = b(je+j) - b(jg+j)
      u1 = b(jf+j) + b(jh+j)
      u3 = ss * ( b(jf+j) - b(jh+j) )
      b(je+j) = b(jb+j)
      a(jb+j) = t0 + t1
      a(jj+j) = t0 - t1
      b(jg+j) = b(jj+j)
      b(jb+j) = u0 + u1
      b(jj+j) = u0 - u1
      a(jf+j) = t2 - u3
      a(jh+j) = t2 + u3
      b(jf+j) = u2 + t3
      b(jh+j) = u2 - t3
C----------------------
      t0 = a(jc+j) + a(jk+j)
      t2 = a(jc+j) - a(jk+j)
      t1 = a(jg+j) + a(jl+j)
      t3 = ss * ( a(jg+j) - a(jl+j) )
      u0 = b(ji+j) + b(jk+j)
      u2 = b(ji+j) - b(jk+j)
      a(jl+j) = a(jo+j)
      u1 = b(jg+j) + b(jl+j)
      u3 = ss * ( b(jg+j) - b(jl+j) )
      b(ji+j) = b(jc+j)
      a(jc+j) = t0 + t1
      a(jk+j) = t0 - t1
      b(jl+j) = b(jo+j)
      b(jc+j) = u0 + u1
      b(jk+j) = u0 - u1
      a(jg+j) = t2 - u3
      a(jo+j) = t2 + u3
      b(jg+j) = u2 + t3
      b(jo+j) = u2 - t3
C----------------------
      t0 = a(jm+j) + a(jl+j)
      t2 = a(jm+j) - a(jl+j)
      t1 = a(jn+j) + a(jp+j)
      t3 = ss * ( a(jn+j) - a(jp+j) )
      a(jm+j) = a(jd+j)
      u0 = b(jd+j) + b(jl+j)
      u2 = b(jd+j) - b(jl+j)
      u1 = b(jn+j) + b(jp+j)
      u3 = ss * ( b(jn+j) - b(jp+j) )
      a(jn+j) = a(jh+j)
      a(jd+j) = t0 + t1
      a(jl+j) = t0 - t1
      b(jd+j) = u0 + u1
      b(jl+j) = u0 - u1
      b(jn+j) = b(jh+j)
      a(jh+j) = t2 - u3
      a(jp+j) = t2 + u3
      b(jh+j) = u2 + t3
      b(jp+j) = u2 - t3
      j = j + jump
  410 continue
C-----( end of loop across transforms )
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  415 continue
  420 continue
  430 continue
C-----( end of double loop for k=0 )
C
C  finished if last pass
C  ---------------------
      if (ipass.eq.m) go to 490
C
      kk = 2*la
C
C     loop on nonzero k
C     -----------------
      do 470 k = ink , jstep-ink , ink
      co1 = trigs(kk+1)
      si1 = s*trigs(kk+2)
      co2 = trigs(2*kk+1)
      si2 = s*trigs(2*kk+2)
      co3 = trigs(3*kk+1)
      si3 = s*trigs(3*kk+2)
C
C  double loop along first transform in block
C  ------------------------------------------
      do 460 ll = k , (la-1)*ink , 4*jstep
C
      do 450 jjj = ll , (n-1)*inc , 4*la*ink
      ja = istart + jjj
C
C     "transverse" loop
C     -----------------
      do 445 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      jc = jb + jstepl
      if (jc.lt.istart) jc = jc + ninc
      jd = jc + jstepl
      if (jd.lt.istart) jd = jd + ninc
      je = ja + laincl
      if (je.lt.istart) je = je + ninc
      jf = je + jstepl
      if (jf.lt.istart) jf = jf + ninc
      jg = jf + jstepl
      if (jg.lt.istart) jg = jg + ninc
      jh = jg + jstepl
      if (jh.lt.istart) jh = jh + ninc
      ji = je + laincl
      if (ji.lt.istart) ji = ji + ninc
      jj = ji + jstepl
      if (jj.lt.istart) jj = jj + ninc
      jk = jj + jstepl
      if (jk.lt.istart) jk = jk + ninc
      jl = jk + jstepl
      if (jl.lt.istart) jl = jl + ninc
      jm = ji + laincl
      if (jm.lt.istart) jm = jm + ninc
      jn = jm + jstepl
      if (jn.lt.istart) jn = jn + ninc
      jo = jn + jstepl
      if (jo.lt.istart) jo = jo + ninc
      jp = jo + jstepl
      if (jp.lt.istart) jp = jp + ninc
      j = 0
C
C  loop across transforms
C  ----------------------
cdir$ ivdep, shortloop
      do 440 l = 1 , nvex
      t0 = a(ja+j) + a(jc+j)
      t2 = a(ja+j) - a(jc+j)
      t1 = a(jb+j) + a(jd+j)
      t3 = ss * ( a(jb+j) - a(jd+j) )
      a(jc+j) = a(ji+j)
      u0 = b(ja+j) + b(jc+j)
      u2 = b(ja+j) - b(jc+j)
      u1 = b(jb+j) + b(jd+j)
      u3 = ss * ( b(jb+j) - b(jd+j) )
      a(jb+j) = a(je+j)
      a(ja+j) = t0 + t1
      b(ja+j) = u0 + u1
      a(je+j) = co1*(t2-u3) - si1*(u2+t3)
      b(jb+j) = si1*(t2-u3) + co1*(u2+t3)
      b(jd+j) = b(jm+j)
      a(ji+j) = co2*(t0-t1) - si2*(u0-u1)
      b(jc+j) = si2*(t0-t1) + co2*(u0-u1)
      a(jd+j) = co3*(t2+u3) - si3*(u2-t3)
      b(jm+j) = si3*(t2+u3) + co3*(u2-t3)
C----------------------------------------
      t0 = a(jb+j) + a(jg+j)
      t2 = a(jb+j) - a(jg+j)
      t1 = a(jf+j) + a(jh+j)
      t3 = ss * ( a(jf+j) - a(jh+j) )
      a(jg+j) = a(jj+j)
      u0 = b(je+j) + b(jg+j)
      u2 = b(je+j) - b(jg+j)
      u1 = b(jf+j) + b(jh+j)
      u3 = ss * ( b(jf+j) - b(jh+j) )
      b(je+j) = b(jb+j)
      a(jb+j) = t0 + t1
      b(jb+j) = u0 + u1
      b(jg+j) = b(jj+j)
      a(jf+j) = co1*(t2-u3) - si1*(u2+t3)
      b(jf+j) = si1*(t2-u3) + co1*(u2+t3)
      a(jj+j) = co2*(t0-t1) - si2*(u0-u1)
      b(jj+j) = si2*(t0-t1) + co2*(u0-u1)
      a(jh+j) = co3*(t2+u3) - si3*(u2-t3)
      b(jh+j) = si3*(t2+u3) + co3*(u2-t3)
C----------------------------------------
      t0 = a(jc+j) + a(jk+j)
      t2 = a(jc+j) - a(jk+j)
      t1 = a(jg+j) + a(jl+j)
      t3 = ss * ( a(jg+j) - a(jl+j) )
      u0 = b(ji+j) + b(jk+j)
      u2 = b(ji+j) - b(jk+j)
      a(jl+j) = a(jo+j)
      u1 = b(jg+j) + b(jl+j)
      u3 = ss * ( b(jg+j) - b(jl+j) )
      b(ji+j) = b(jc+j)
      a(jc+j) = t0 + t1
      b(jc+j) = u0 + u1
      b(jl+j) = b(jo+j)
      a(jg+j) = co1*(t2-u3) - si1*(u2+t3)
      b(jg+j) = si1*(t2-u3) + co1*(u2+t3)
      a(jk+j) = co2*(t0-t1) - si2*(u0-u1)
      b(jk+j) = si2*(t0-t1) + co2*(u0-u1)
      a(jo+j) = co3*(t2+u3) - si3*(u2-t3)
      b(jo+j) = si3*(t2+u3) + co3*(u2-t3)
C----------------------------------------
      t0 = a(jm+j) + a(jl+j)
      t2 = a(jm+j) - a(jl+j)
      t1 = a(jn+j) + a(jp+j)
      t3 = ss * ( a(jn+j) - a(jp+j) )
      a(jm+j) = a(jd+j)
      u0 = b(jd+j) + b(jl+j)
      u2 = b(jd+j) - b(jl+j)
      a(jn+j) = a(jh+j)
      u1 = b(jn+j) + b(jp+j)
      u3 = ss * ( b(jn+j) - b(jp+j) )
      b(jn+j) = b(jh+j)
      a(jd+j) = t0 + t1
      b(jd+j) = u0 + u1
      a(jh+j) = co1*(t2-u3) - si1*(u2+t3)
      b(jh+j) = si1*(t2-u3) + co1*(u2+t3)
      a(jl+j) = co2*(t0-t1) - si2*(u0-u1)
      b(jl+j) = si2*(t0-t1) + co2*(u0-u1)
      a(jp+j) = co3*(t2+u3) - si3*(u2-t3)
      b(jp+j) = si3*(t2+u3) + co3*(u2-t3)
      j = j + jump
  440 continue
C-----(end of loop across transforms)
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  445 continue
  450 continue
  460 continue
C-----( end of double loop for this k )
      kk = kk + 2*la
  470 continue
C-----( end of loop over values of k )
      la = 4*la
  480 continue
C-----( end of loop on type II radix-4 passes )
C-----( nvex transforms completed)
  490 continue
      istart = istart + nvex * jump
  500 continue
C-----( end of loop on blocks of transforms )
C
      return
      end
C     fortran version of *gpfa3* -
C     radix-3 section of self-sorting, in-place
C        generalized PFA
C
C-------------------------------------------------------------------
C
      subroutine gpfa3f(a,b,trigs,inc,jump,n,mm,lot,isign)
      dimension a(*), b(*), trigs(*)
      data sin60/0.866025403784437/
      data lvr/64/
C
C     ***************************************************************
C     *                                                             *
C     *  N.B. LVR = LENGTH OF VECTOR REGISTERS, SET TO 128 FOR C90. *
C     *  RESET TO 64 FOR OTHER CRAY MACHINES, OR TO ANY LARGE VALUE *
C     *  (GREATER THAN OR EQUAL TO LOT) FOR A SCALAR COMPUTER.      *
C     *                                                             *
C     ***************************************************************
C
      n3 = 3**mm
      inq = n/n3
      jstepx = (n3-n) * inc
      ninc = n * inc
      ink = inc * inq
      mu = mod(inq,3)
      if (isign.eq.-1) mu = 3-mu
      m = mm
      mh = (m+1)/2
      s = float(isign)
      c1 = sin60
      if (mu.eq.2) c1 = -c1
C
      nblox = 1 + (lot-1)/lvr
      left = lot
      s = float(isign)
      istart = 1
C
C  loop on blocks of lvr transforms
C  --------------------------------
      do 500 nb = 1 , nblox
C
      if (left.le.lvr) then
         nvex = left
      else if (left.lt.(2*lvr)) then
         nvex = left/2
         nvex = nvex + mod(nvex,2)
      else
         nvex = lvr
      endif
      left = left - nvex
C
      la = 1
C
C  loop on type I radix-3 passes
C  -----------------------------
      do 160 ipass = 1 , mh
      jstep = (n*inc) / (3*la)
      jstepl = jstep - ninc
C
C  k = 0 loop (no twiddle factors)
C  -------------------------------
      do 120 jjj = 0 , (n-1)*inc , 3*jstep
      ja = istart + jjj
C
C  "transverse" loop
C  -----------------
      do 115 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      jc = jb + jstepl
      if (jc.lt.istart) jc = jc + ninc
      j = 0
C
C  loop across transforms
C  ----------------------
cdir$ ivdep, shortloop
      do 110 l = 1 , nvex
      t1 = a(jb+j) + a(jc+j)
      t2 = a(ja+j) - 0.5 * t1
      t3 = c1 * ( a(jb+j) - a(jc+j) )
      u1 = b(jb+j) + b(jc+j)
      u2 = b(ja+j) - 0.5 * u1
      u3 = c1 * ( b(jb+j) - b(jc+j) )
      a(ja+j) = a(ja+j) + t1
      b(ja+j) = b(ja+j) + u1
      a(jb+j) = t2 - u3
      b(jb+j) = u2 + t3
      a(jc+j) = t2 + u3
      b(jc+j) = u2 - t3
      j = j + jump
  110 continue
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  115 continue
  120 continue
C
C  finished if n3 = 3
C  ------------------
      if (n3.eq.3) go to 490
      kk = 2 * la
C
C  loop on nonzero k
C  -----------------
      do 150 k = ink , jstep-ink , ink
      co1 = trigs(kk+1)
      si1 = s*trigs(kk+2)
      co2 = trigs(2*kk+1)
      si2 = s*trigs(2*kk+2)
C
C  loop along transform
C  --------------------
      do 140 jjj = k , (n-1)*inc , 3*jstep
      ja = istart + jjj
C
C  "transverse" loop
C  -----------------
      do 135 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      jc = jb + jstepl
      if (jc.lt.istart) jc = jc + ninc
      j = 0
C
C  loop across transforms
C  ----------------------
cdir$ ivdep,shortloop
      do 130 l = 1 , nvex
      t1 = a(jb+j) + a(jc+j)
      t2 = a(ja+j) - 0.5 * t1
      t3 = c1 * ( a(jb+j) - a(jc+j) )
      u1 = b(jb+j) + b(jc+j)
      u2 = b(ja+j) - 0.5 * u1
      u3 = c1 * ( b(jb+j) - b(jc+j) )
      a(ja+j) = a(ja+j) + t1
      b(ja+j) = b(ja+j) + u1
      a(jb+j) = co1*(t2-u3) - si1*(u2+t3)
      b(jb+j) = si1*(t2-u3) + co1*(u2+t3)
      a(jc+j) = co2*(t2+u3) - si2*(u2-t3)
      b(jc+j) = si2*(t2+u3) + co2*(u2-t3)
      j = j + jump
  130 continue
C-----( end of loop across transforms )
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  135 continue
  140 continue
C-----( end of loop along transforms )
      kk = kk + 2*la
  150 continue
C-----( end of loop on nonzero k )
      la = 3*la
  160 continue
C-----( end of loop on type I radix-3 passes)
C
C  loop on type II radix-3 passes
C  ------------------------------
  400 continue
C
      do 480 ipass = mh+1 , m
      jstep = (n*inc) / (3*la)
      jstepl = jstep - ninc
      laincl = la*ink - ninc
C
C  k=0 loop (no twiddle factors)
C  -----------------------------
      do 430 ll = 0 , (la-1)*ink , 3*jstep
C
      do 420 jjj = ll , (n-1)*inc , 3*la*ink
      ja = istart + jjj
C
C  "transverse" loop
C  -----------------
      do 415 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      jc = jb + jstepl
      if (jc.lt.istart) jc = jc + ninc
      jd = ja + laincl
      if (jd.lt.istart) jd = jd + ninc
      je = jd + jstepl
      if (je.lt.istart) je = je + ninc
      jf = je + jstepl
      if (jf.lt.istart) jf = jf + ninc
      jg = jd + laincl
      if (jg.lt.istart) jg = jg + ninc
      jh = jg + jstepl
      if (jh.lt.istart) jh = jh + ninc
      ji = jh + jstepl
      if (ji.lt.istart) ji = ji + ninc
      j = 0
C
C  loop across transforms
C  ----------------------
cdir$ ivdep, shortloop
      do 410 l = 1 , nvex
      t1 = a(jb+j) + a(jc+j)
      t2 = a(ja+j) - 0.5 * t1
      t3 = c1 * ( a(jb+j) - a(jc+j) )
      a(jb+j) = a(jd+j)
      u1 = b(jb+j) + b(jc+j)
      u2 = b(ja+j) - 0.5 * u1
      u3 = c1 * ( b(jb+j) - b(jc+j) )
      b(jb+j) = b(jd+j)
      a(ja+j) = a(ja+j) + t1
      b(ja+j) = b(ja+j) + u1
      a(jd+j) = t2 - u3
      b(jd+j) = u2 + t3
      a(jc+j) = t2 + u3
      b(jc+j) = u2 - t3
C----------------------
      t1 = a(je+j) + a(jf+j)
      t2 = a(jb+j) - 0.5 * t1
      t3 = c1 * ( a(je+j) - a(jf+j) )
      a(jf+j) = a(jh+j)
      u1 = b(je+j) + b(jf+j)
      u2 = b(jb+j) - 0.5 * u1
      u3 = c1 * ( b(je+j) - b(jf+j) )
      b(jf+j) = b(jh+j)
      a(jb+j) = a(jb+j) + t1
      b(jb+j) = b(jb+j) + u1
      a(je+j) = t2 - u3
      b(je+j) = u2 + t3
      a(jh+j) = t2 + u3
      b(jh+j) = u2 - t3
C----------------------
      t1 = a(jf+j) + a(ji+j)
      t2 = a(jg+j) - 0.5 * t1
      t3 = c1 * ( a(jf+j) - a(ji+j) )
      t1 = a(jg+j) + t1
      a(jg+j) = a(jc+j)
      u1 = b(jf+j) + b(ji+j)
      u2 = b(jg+j) - 0.5 * u1
      u3 = c1 * ( b(jf+j) - b(ji+j) )
      u1 = b(jg+j) + u1
      b(jg+j) = b(jc+j)
      a(jc+j) = t1
      b(jc+j) = u1
      a(jf+j) = t2 - u3
      b(jf+j) = u2 + t3
      a(ji+j) = t2 + u3
      b(ji+j) = u2 - t3
      j = j + jump
  410 continue
C-----( end of loop across transforms )
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  415 continue
  420 continue
  430 continue
C-----( end of double loop for k=0 )
C
C  finished if last pass
C  ---------------------
      if (ipass.eq.m) go to 490
C
      kk = 2*la
C
C     loop on nonzero k
C     -----------------
      do 470 k = ink , jstep-ink , ink
      co1 = trigs(kk+1)
      si1 = s*trigs(kk+2)
      co2 = trigs(2*kk+1)
      si2 = s*trigs(2*kk+2)
C
C  double loop along first transform in block
C  ------------------------------------------
      do 460 ll = k , (la-1)*ink , 3*jstep
C
      do 450 jjj = ll , (n-1)*inc , 3*la*ink
      ja = istart + jjj
C
C  "transverse" loop
C  -----------------
      do 445 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      jc = jb + jstepl
      if (jc.lt.istart) jc = jc + ninc
      jd = ja + laincl
      if (jd.lt.istart) jd = jd + ninc
      je = jd + jstepl
      if (je.lt.istart) je = je + ninc
      jf = je + jstepl
      if (jf.lt.istart) jf = jf + ninc
      jg = jd + laincl
      if (jg.lt.istart) jg = jg + ninc
      jh = jg + jstepl
      if (jh.lt.istart) jh = jh + ninc
      ji = jh + jstepl
      if (ji.lt.istart) ji = ji + ninc
      j = 0
C
C  loop across transforms
C  ----------------------
cdir$ ivdep, shortloop
      do 440 l = 1 , nvex
      t1 = a(jb+j) + a(jc+j)
      t2 = a(ja+j) - 0.5 * t1
      t3 = c1 * ( a(jb+j) - a(jc+j) )
      a(jb+j) = a(jd+j)
      u1 = b(jb+j) + b(jc+j)
      u2 = b(ja+j) - 0.5 * u1
      u3 = c1 * ( b(jb+j) - b(jc+j) )
      b(jb+j) = b(jd+j)
      a(ja+j) = a(ja+j) + t1
      b(ja+j) = b(ja+j) + u1
      a(jd+j) = co1*(t2-u3) - si1*(u2+t3)
      b(jd+j) = si1*(t2-u3) + co1*(u2+t3)
      a(jc+j) = co2*(t2+u3) - si2*(u2-t3)
      b(jc+j) = si2*(t2+u3) + co2*(u2-t3)
C----------------------
      t1 = a(je+j) + a(jf+j)
      t2 = a(jb+j) - 0.5 * t1
      t3 = c1 * ( a(je+j) - a(jf+j) )
      a(jf+j) = a(jh+j)
      u1 = b(je+j) + b(jf+j)
      u2 = b(jb+j) - 0.5 * u1
      u3 = c1 * ( b(je+j) - b(jf+j) )
      b(jf+j) = b(jh+j)
      a(jb+j) = a(jb+j) + t1
      b(jb+j) = b(jb+j) + u1
      a(je+j) = co1*(t2-u3) - si1*(u2+t3)
      b(je+j) = si1*(t2-u3) + co1*(u2+t3)
      a(jh+j) = co2*(t2+u3) - si2*(u2-t3)
      b(jh+j) = si2*(t2+u3) + co2*(u2-t3)
C----------------------
      t1 = a(jf+j) + a(ji+j)
      t2 = a(jg+j) - 0.5 * t1
      t3 = c1 * ( a(jf+j) - a(ji+j) )
      t1 = a(jg+j) + t1
      a(jg+j) = a(jc+j)
      u1 = b(jf+j) + b(ji+j)
      u2 = b(jg+j) - 0.5 * u1
      u3 = c1 * ( b(jf+j) - b(ji+j) )
      u1 = b(jg+j) + u1
      b(jg+j) = b(jc+j)
      a(jc+j) = t1
      b(jc+j) = u1
      a(jf+j) = co1*(t2-u3) - si1*(u2+t3)
      b(jf+j) = si1*(t2-u3) + co1*(u2+t3)
      a(ji+j) = co2*(t2+u3) - si2*(u2-t3)
      b(ji+j) = si2*(t2+u3) + co2*(u2-t3)
      j = j + jump
  440 continue
C-----(end of loop across transforms)
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  445 continue
  450 continue
  460 continue
C-----( end of double loop for this k )
      kk = kk + 2*la
  470 continue
C-----( end of loop over values of k )
      la = 3*la
  480 continue
C-----( end of loop on type II radix-3 passes )
C-----( nvex transforms completed)
  490 continue
      istart = istart + nvex * jump
  500 continue
C-----( end of loop on blocks of transforms )
C
      return
      end
C     fortran version of *gpfa5* -
C     radix-5 section of self-sorting, in-place,
C        generalized pfa
C
C-------------------------------------------------------------------
C
      subroutine gpfa5f(a,b,trigs,inc,jump,n,mm,lot,isign)
      dimension a(*), b(*), trigs(*)
      data sin36/0.587785252292473/, sin72/0.951056516295154/,
     *      qrt5/0.559016994374947/
      data lvr/64/
C
C     ***************************************************************
C     *                                                             *
C     *  N.B. LVR = LENGTH OF VECTOR REGISTERS, SET TO 128 FOR C90. *
C     *  RESET TO 64 FOR OTHER CRAY MACHINES, OR TO ANY LARGE VALUE *
C     *  (GREATER THAN OR EQUAL TO LOT) FOR A SCALAR COMPUTER.      *
C     *                                                             *
C     ***************************************************************
C
      n5 = 5 ** mm
      inq = n / n5
      jstepx = (n5-n) * inc
      ninc = n * inc
      ink = inc * inq
      mu = mod(inq,5)
      if (isign.eq.-1) mu = 5 - mu
C
      m = mm
      mh = (m+1)/2
      s = float(isign)
      c1 = qrt5
      c2 = sin72
      c3 = sin36
      if (mu.eq.2.or.mu.eq.3) then
         c1 = -c1
         c2 = sin36
         c3 = sin72
      endif
      if (mu.eq.3.or.mu.eq.4) c2 = -c2
      if (mu.eq.2.or.mu.eq.4) c3 = -c3
C
      nblox = 1 + (lot-1)/lvr
      left = lot
      s = float(isign)
      istart = 1
C
C  loop on blocks of lvr transforms
C  --------------------------------
      do 500 nb = 1 , nblox
C
      if (left.le.lvr) then
         nvex = left
      else if (left.lt.(2*lvr)) then
         nvex = left/2
         nvex = nvex + mod(nvex,2)
      else
         nvex = lvr
      endif
      left = left - nvex
C
      la = 1
C
C  loop on type I radix-5 passes
C  -----------------------------
      do 160 ipass = 1 , mh
      jstep = (n*inc) / (5*la)
      jstepl = jstep - ninc
      kk = 0
C
C  loop on k
C  ---------
      do 150 k = 0 , jstep-ink , ink
C
      if (k.gt.0) then
      co1 = trigs(kk+1)
      si1 = s*trigs(kk+2)
      co2 = trigs(2*kk+1)
      si2 = s*trigs(2*kk+2)
      co3 = trigs(3*kk+1)
      si3 = s*trigs(3*kk+2)
      co4 = trigs(4*kk+1)
      si4 = s*trigs(4*kk+2)
      endif
C
C  loop along transform
C  --------------------
      do 140 jjj = k , (n-1)*inc , 5*jstep
      ja = istart + jjj
C
C     "transverse" loop
C     -----------------
      do 135 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      jc = jb + jstepl
      if (jc.lt.istart) jc = jc + ninc
      jd = jc + jstepl
      if (jd.lt.istart) jd = jd + ninc
      je = jd + jstepl
      if (je.lt.istart) je = je + ninc
      j = 0
C
C  loop across transforms
C  ----------------------
      if (k.eq.0) then
C
cdir$ ivdep, shortloop
      do 110 l = 1 , nvex
      t1 = a(jb+j) + a(je+j)
      t2 = a(jc+j) + a(jd+j)
      t3 = a(jb+j) - a(je+j)
      t4 = a(jc+j) - a(jd+j)
      t5 = t1 + t2
      t6 = c1 * ( t1 - t2 )
      t7 = a(ja+j) - 0.25 * t5
      a(ja+j) = a(ja+j) + t5
      t8 = t7 + t6
      t9 = t7 - t6
      t10 = c3 * t3 - c2 * t4
      t11 = c2 * t3 + c3 * t4
      u1 = b(jb+j) + b(je+j)
      u2 = b(jc+j) + b(jd+j)
      u3 = b(jb+j) - b(je+j)
      u4 = b(jc+j) - b(jd+j)
      u5 = u1 + u2
      u6 = c1 * ( u1 - u2 )
      u7 = b(ja+j) - 0.25 * u5
      b(ja+j) = b(ja+j) + u5
      u8 = u7 + u6
      u9 = u7 - u6
      u10 = c3 * u3 - c2 * u4
      u11 = c2 * u3 + c3 * u4
      a(jb+j) = t8 - u11
      b(jb+j) = u8 + t11
      a(je+j) = t8 + u11
      b(je+j) = u8 - t11
      a(jc+j) = t9 - u10
      b(jc+j) = u9 + t10
      a(jd+j) = t9 + u10
      b(jd+j) = u9 - t10
      j = j + jump
  110 continue
C
      else
C
cdir$ ivdep,shortloop
      do 130 l = 1 , nvex
      t1 = a(jb+j) + a(je+j)
      t2 = a(jc+j) + a(jd+j)
      t3 = a(jb+j) - a(je+j)
      t4 = a(jc+j) - a(jd+j)
      t5 = t1 + t2
      t6 = c1 * ( t1 - t2 )
      t7 = a(ja+j) - 0.25 * t5
      a(ja+j) = a(ja+j) + t5
      t8 = t7 + t6
      t9 = t7 - t6
      t10 = c3 * t3 - c2 * t4
      t11 = c2 * t3 + c3 * t4
      u1 = b(jb+j) + b(je+j)
      u2 = b(jc+j) + b(jd+j)
      u3 = b(jb+j) - b(je+j)
      u4 = b(jc+j) - b(jd+j)
      u5 = u1 + u2
      u6 = c1 * ( u1 - u2 )
      u7 = b(ja+j) - 0.25 * u5
      b(ja+j) = b(ja+j) + u5
      u8 = u7 + u6
      u9 = u7 - u6
      u10 = c3 * u3 - c2 * u4
      u11 = c2 * u3 + c3 * u4
      a(jb+j) = co1*(t8-u11) - si1*(u8+t11)
      b(jb+j) = si1*(t8-u11) + co1*(u8+t11)
      a(je+j) = co4*(t8+u11) - si4*(u8-t11)
      b(je+j) = si4*(t8+u11) + co4*(u8-t11)
      a(jc+j) = co2*(t9-u10) - si2*(u9+t10)
      b(jc+j) = si2*(t9-u10) + co2*(u9+t10)
      a(jd+j) = co3*(t9+u10) - si3*(u9-t10)
      b(jd+j) = si3*(t9+u10) + co3*(u9-t10)
      j = j + jump
  130 continue
C
      endif
C
C-----( end of loop across transforms )
C
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  135 continue
  140 continue
C-----( end of loop along transforms )
      kk = kk + 2*la
  150 continue
C-----( end of loop on nonzero k )
      la = 5*la
  160 continue
C-----( end of loop on type I radix-5 passes)
C
      if (n.eq.5) go to 490
C
C  loop on type II radix-5 passes
C  ------------------------------
  400 continue
C
      do 480 ipass = mh+1 , m
      jstep = (n*inc) / (5*la)
      jstepl = jstep - ninc
      laincl = la * ink - ninc
      kk = 0
C
C     loop on k
C     ---------
      do 470 k = 0 , jstep-ink , ink
C
      if (k.gt.0) then
      co1 = trigs(kk+1)
      si1 = s*trigs(kk+2)
      co2 = trigs(2*kk+1)
      si2 = s*trigs(2*kk+2)
      co3 = trigs(3*kk+1)
      si3 = s*trigs(3*kk+2)
      co4 = trigs(4*kk+1)
      si4 = s*trigs(4*kk+2)
      endif
C
C  double loop along first transform in block
C  ------------------------------------------
      do 460 ll = k , (la-1)*ink , 5*jstep
C
      do 450 jjj = ll , (n-1)*inc , 5*la*ink
      ja = istart + jjj
C
C     "transverse" loop
C     -----------------
      do 445 nu = 1 , inq
      jb = ja + jstepl
      if (jb.lt.istart) jb = jb + ninc
      jc = jb + jstepl
      if (jc.lt.istart) jc = jc + ninc
      jd = jc + jstepl
      if (jd.lt.istart) jd = jd + ninc
      je = jd + jstepl
      if (je.lt.istart) je = je + ninc
      jf = ja + laincl
      if (jf.lt.istart) jf = jf + ninc
      jg = jf + jstepl
      if (jg.lt.istart) jg = jg + ninc
      jh = jg + jstepl
      if (jh.lt.istart) jh = jh + ninc
      ji = jh + jstepl
      if (ji.lt.istart) ji = ji + ninc
      jj = ji + jstepl
      if (jj.lt.istart) jj = jj + ninc
      jk = jf + laincl
      if (jk.lt.istart) jk = jk + ninc
      jl = jk + jstepl
      if (jl.lt.istart) jl = jl + ninc
      jm = jl + jstepl
      if (jm.lt.istart) jm = jm + ninc
      jn = jm + jstepl
      if (jn.lt.istart) jn = jn + ninc
      jo = jn + jstepl
      if (jo.lt.istart) jo = jo + ninc
      jp = jk + laincl
      if (jp.lt.istart) jp = jp + ninc
      jq = jp + jstepl
      if (jq.lt.istart) jq = jq + ninc
      jr = jq + jstepl
      if (jr.lt.istart) jr = jr + ninc
      js = jr + jstepl
      if (js.lt.istart) js = js + ninc
      jt = js + jstepl
      if (jt.lt.istart) jt = jt + ninc
      ju = jp + laincl
      if (ju.lt.istart) ju = ju + ninc
      jv = ju + jstepl
      if (jv.lt.istart) jv = jv + ninc
      jw = jv + jstepl
      if (jw.lt.istart) jw = jw + ninc
      jx = jw + jstepl
      if (jx.lt.istart) jx = jx + ninc
      jy = jx + jstepl
      if (jy.lt.istart) jy = jy + ninc
      j = 0
C
C  loop across transforms
C  ----------------------
      if (k.eq.0) then
C
cdir$ ivdep, shortloop
      do 410 l = 1 , nvex
      t1 = a(jb+j) + a(je+j)
      t2 = a(jc+j) + a(jd+j)
      t3 = a(jb+j) - a(je+j)
      t4 = a(jc+j) - a(jd+j)
      a(jb+j) = a(jf+j)
      t5 = t1 + t2
      t6 = c1 * ( t1 - t2 )
      t7 = a(ja+j) - 0.25 * t5
      a(ja+j) = a(ja+j) + t5
      t8 = t7 + t6
      t9 = t7 - t6
      a(jc+j) = a(jk+j)
      t10 = c3 * t3 - c2 * t4
      t11 = c2 * t3 + c3 * t4
      u1 = b(jb+j) + b(je+j)
      u2 = b(jc+j) + b(jd+j)
      u3 = b(jb+j) - b(je+j)
      u4 = b(jc+j) - b(jd+j)
      b(jb+j) = b(jf+j)
      u5 = u1 + u2
      u6 = c1 * ( u1 - u2 )
      u7 = b(ja+j) - 0.25 * u5
      b(ja+j) = b(ja+j) + u5
      u8 = u7 + u6
      u9 = u7 - u6
      b(jc+j) = b(jk+j)
      u10 = c3 * u3 - c2 * u4
      u11 = c2 * u3 + c3 * u4
      a(jf+j) = t8 - u11
      b(jf+j) = u8 + t11
      a(je+j) = t8 + u11
      b(je+j) = u8 - t11
      a(jk+j) = t9 - u10
      b(jk+j) = u9 + t10
      a(jd+j) = t9 + u10
      b(jd+j) = u9 - t10
C----------------------
      t1 = a(jg+j) + a(jj+j)
      t2 = a(jh+j) + a(ji+j)
      t3 = a(jg+j) - a(jj+j)
      t4 = a(jh+j) - a(ji+j)
      a(jh+j) = a(jl+j)
      t5 = t1 + t2
      t6 = c1 * ( t1 - t2 )
      t7 = a(jb+j) - 0.25 * t5
      a(jb+j) = a(jb+j) + t5
      t8 = t7 + t6
      t9 = t7 - t6
      a(ji+j) = a(jq+j)
      t10 = c3 * t3 - c2 * t4
      t11 = c2 * t3 + c3 * t4
      u1 = b(jg+j) + b(jj+j)
      u2 = b(jh+j) + b(ji+j)
      u3 = b(jg+j) - b(jj+j)
      u4 = b(jh+j) - b(ji+j)
      b(jh+j) = b(jl+j)
      u5 = u1 + u2
      u6 = c1 * ( u1 - u2 )
      u7 = b(jb+j) - 0.25 * u5
      b(jb+j) = b(jb+j) + u5
      u8 = u7 + u6
      u9 = u7 - u6
      b(ji+j) = b(jq+j)
      u10 = c3 * u3 - c2 * u4
      u11 = c2 * u3 + c3 * u4
      a(jg+j) = t8 - u11
      b(jg+j) = u8 + t11
      a(jj+j) = t8 + u11
      b(jj+j) = u8 - t11
      a(jl+j) = t9 - u10
      b(jl+j) = u9 + t10
      a(jq+j) = t9 + u10
      b(jq+j) = u9 - t10
C----------------------
      t1 = a(jh+j) + a(jo+j)
      t2 = a(jm+j) + a(jn+j)
      t3 = a(jh+j) - a(jo+j)
      t4 = a(jm+j) - a(jn+j)
      a(jn+j) = a(jr+j)
      t5 = t1 + t2
      t6 = c1 * ( t1 - t2 )
      t7 = a(jc+j) - 0.25 * t5
      a(jc+j) = a(jc+j) + t5
      t8 = t7 + t6
      t9 = t7 - t6
      a(jo+j) = a(jw+j)
      t10 = c3 * t3 - c2 * t4
      t11 = c2 * t3 + c3 * t4
      u1 = b(jh+j) + b(jo+j)
      u2 = b(jm+j) + b(jn+j)
      u3 = b(jh+j) - b(jo+j)
      u4 = b(jm+j) - b(jn+j)
      b(jn+j) = b(jr+j)
      u5 = u1 + u2
      u6 = c1 * ( u1 - u2 )
      u7 = b(jc+j) - 0.25 * u5
      b(jc+j) = b(jc+j) + u5
      u8 = u7 + u6
      u9 = u7 - u6
      b(jo+j) = b(jw+j)
      u10 = c3 * u3 - c2 * u4
      u11 = c2 * u3 + c3 * u4
      a(jh+j) = t8 - u11
      b(jh+j) = u8 + t11
      a(jw+j) = t8 + u11
      b(jw+j) = u8 - t11
      a(jm+j) = t9 - u10
      b(jm+j) = u9 + t10
      a(jr+j) = t9 + u10
      b(jr+j) = u9 - t10
C----------------------
      t1 = a(ji+j) + a(jt+j)
      t2 = a(jn+j) + a(js+j)
      t3 = a(ji+j) - a(jt+j)
      t4 = a(jn+j) - a(js+j)
      a(jt+j) = a(jx+j)
      t5 = t1 + t2
      t6 = c1 * ( t1 - t2 )
      t7 = a(jp+j) - 0.25 * t5
      ax = a(jp+j) + t5
      t8 = t7 + t6
      t9 = t7 - t6
      a(jp+j) = a(jd+j)
      t10 = c3 * t3 - c2 * t4
      t11 = c2 * t3 + c3 * t4
      a(jd+j) = ax
      u1 = b(ji+j) + b(jt+j)
      u2 = b(jn+j) + b(js+j)
      u3 = b(ji+j) - b(jt+j)
      u4 = b(jn+j) - b(js+j)
      b(jt+j) = b(jx+j)
      u5 = u1 + u2
      u6 = c1 * ( u1 - u2 )
      u7 = b(jp+j) - 0.25 * u5
      bx = b(jp+j) + u5
      u8 = u7 + u6
      u9 = u7 - u6
      b(jp+j) = b(jd+j)
      u10 = c3 * u3 - c2 * u4
      u11 = c2 * u3 + c3 * u4
      b(jd+j) = bx
      a(ji+j) = t8 - u11
      b(ji+j) = u8 + t11
      a(jx+j) = t8 + u11
      b(jx+j) = u8 - t11
      a(jn+j) = t9 - u10
      b(jn+j) = u9 + t10
      a(js+j) = t9 + u10
      b(js+j) = u9 - t10
C----------------------
      t1 = a(jv+j) + a(jy+j)
      t2 = a(jo+j) + a(jt+j)
      t3 = a(jv+j) - a(jy+j)
      t4 = a(jo+j) - a(jt+j)
      a(jv+j) = a(jj+j)
      t5 = t1 + t2
      t6 = c1 * ( t1 - t2 )
      t7 = a(ju+j) - 0.25 * t5
      ax = a(ju+j) + t5
      t8 = t7 + t6
      t9 = t7 - t6
      a(ju+j) = a(je+j)
      t10 = c3 * t3 - c2 * t4
      t11 = c2 * t3 + c3 * t4
      a(je+j) = ax
      u1 = b(jv+j) + b(jy+j)
      u2 = b(jo+j) + b(jt+j)
      u3 = b(jv+j) - b(jy+j)
      u4 = b(jo+j) - b(jt+j)
      b(jv+j) = b(jj+j)
      u5 = u1 + u2
      u6 = c1 * ( u1 - u2 )
      u7 = b(ju+j) - 0.25 * u5
      bx = b(ju+j) + u5
      u8 = u7 + u6
      u9 = u7 - u6
      b(ju+j) = b(je+j)
      u10 = c3 * u3 - c2 * u4
      u11 = c2 * u3 + c3 * u4
      b(je+j) = bx
      a(jj+j) = t8 - u11
      b(jj+j) = u8 + t11
      a(jy+j) = t8 + u11
      b(jy+j) = u8 - t11
      a(jo+j) = t9 - u10
      b(jo+j) = u9 + t10
      a(jt+j) = t9 + u10
      b(jt+j) = u9 - t10
      j = j + jump
  410 continue
C
      else
C
cdir$ ivdep, shortloop
      do 440 l = 1 , nvex
      t1 = a(jb+j) + a(je+j)
      t2 = a(jc+j) + a(jd+j)
      t3 = a(jb+j) - a(je+j)
      t4 = a(jc+j) - a(jd+j)
      a(jb+j) = a(jf+j)
      t5 = t1 + t2
      t6 = c1 * ( t1 - t2 )
      t7 = a(ja+j) - 0.25 * t5
      a(ja+j) = a(ja+j) + t5
      t8 = t7 + t6
      t9 = t7 - t6
      a(jc+j) = a(jk+j)
      t10 = c3 * t3 - c2 * t4
      t11 = c2 * t3 + c3 * t4
      u1 = b(jb+j) + b(je+j)
      u2 = b(jc+j) + b(jd+j)
      u3 = b(jb+j) - b(je+j)
      u4 = b(jc+j) - b(jd+j)
      b(jb+j) = b(jf+j)
      u5 = u1 + u2
      u6 = c1 * ( u1 - u2 )
      u7 = b(ja+j) - 0.25 * u5
      b(ja+j) = b(ja+j) + u5
      u8 = u7 + u6
      u9 = u7 - u6
      b(jc+j) = b(jk+j)
      u10 = c3 * u3 - c2 * u4
      u11 = c2 * u3 + c3 * u4
      a(jf+j) = co1*(t8-u11) - si1*(u8+t11)
      b(jf+j) = si1*(t8-u11) + co1*(u8+t11)
      a(je+j) = co4*(t8+u11) - si4*(u8-t11)
      b(je+j) = si4*(t8+u11) + co4*(u8-t11)
      a(jk+j) = co2*(t9-u10) - si2*(u9+t10)
      b(jk+j) = si2*(t9-u10) + co2*(u9+t10)
      a(jd+j) = co3*(t9+u10) - si3*(u9-t10)
      b(jd+j) = si3*(t9+u10) + co3*(u9-t10)
C----------------------
      t1 = a(jg+j) + a(jj+j)
      t2 = a(jh+j) + a(ji+j)
      t3 = a(jg+j) - a(jj+j)
      t4 = a(jh+j) - a(ji+j)
      a(jh+j) = a(jl+j)
      t5 = t1 + t2
      t6 = c1 * ( t1 - t2 )
      t7 = a(jb+j) - 0.25 * t5
      a(jb+j) = a(jb+j) + t5
      t8 = t7 + t6
      t9 = t7 - t6
      a(ji+j) = a(jq+j)
      t10 = c3 * t3 - c2 * t4
      t11 = c2 * t3 + c3 * t4
      u1 = b(jg+j) + b(jj+j)
      u2 = b(jh+j) + b(ji+j)
      u3 = b(jg+j) - b(jj+j)
      u4 = b(jh+j) - b(ji+j)
      b(jh+j) = b(jl+j)
      u5 = u1 + u2
      u6 = c1 * ( u1 - u2 )
      u7 = b(jb+j) - 0.25 * u5
      b(jb+j) = b(jb+j) + u5
      u8 = u7 + u6
      u9 = u7 - u6
      b(ji+j) = b(jq+j)
      u10 = c3 * u3 - c2 * u4
      u11 = c2 * u3 + c3 * u4
      a(jg+j) = co1*(t8-u11) - si1*(u8+t11)
      b(jg+j) = si1*(t8-u11) + co1*(u8+t11)
      a(jj+j) = co4*(t8+u11) - si4*(u8-t11)
      b(jj+j) = si4*(t8+u11) + co4*(u8-t11)
      a(jl+j) = co2*(t9-u10) - si2*(u9+t10)
      b(jl+j) = si2*(t9-u10) + co2*(u9+t10)
      a(jq+j) = co3*(t9+u10) - si3*(u9-t10)
      b(jq+j) = si3*(t9+u10) + co3*(u9-t10)
C----------------------
      t1 = a(jh+j) + a(jo+j)
      t2 = a(jm+j) + a(jn+j)
      t3 = a(jh+j) - a(jo+j)
      t4 = a(jm+j) - a(jn+j)
      a(jn+j) = a(jr+j)
      t5 = t1 + t2
      t6 = c1 * ( t1 - t2 )
      t7 = a(jc+j) - 0.25 * t5
      a(jc+j) = a(jc+j) + t5
      t8 = t7 + t6
      t9 = t7 - t6
      a(jo+j) = a(jw+j)
      t10 = c3 * t3 - c2 * t4
      t11 = c2 * t3 + c3 * t4
      u1 = b(jh+j) + b(jo+j)
      u2 = b(jm+j) + b(jn+j)
      u3 = b(jh+j) - b(jo+j)
      u4 = b(jm+j) - b(jn+j)
      b(jn+j) = b(jr+j)
      u5 = u1 + u2
      u6 = c1 * ( u1 - u2 )
      u7 = b(jc+j) - 0.25 * u5
      b(jc+j) = b(jc+j) + u5
      u8 = u7 + u6
      u9 = u7 - u6
      b(jo+j) = b(jw+j)
      u10 = c3 * u3 - c2 * u4
      u11 = c2 * u3 + c3 * u4
      a(jh+j) = co1*(t8-u11) - si1*(u8+t11)
      b(jh+j) = si1*(t8-u11) + co1*(u8+t11)
      a(jw+j) = co4*(t8+u11) - si4*(u8-t11)
      b(jw+j) = si4*(t8+u11) + co4*(u8-t11)
      a(jm+j) = co2*(t9-u10) - si2*(u9+t10)
      b(jm+j) = si2*(t9-u10) + co2*(u9+t10)
      a(jr+j) = co3*(t9+u10) - si3*(u9-t10)
      b(jr+j) = si3*(t9+u10) + co3*(u9-t10)
C----------------------
      t1 = a(ji+j) + a(jt+j)
      t2 = a(jn+j) + a(js+j)
      t3 = a(ji+j) - a(jt+j)
      t4 = a(jn+j) - a(js+j)
      a(jt+j) = a(jx+j)
      t5 = t1 + t2
      t6 = c1 * ( t1 - t2 )
      t7 = a(jp+j) - 0.25 * t5
      ax = a(jp+j) + t5
      t8 = t7 + t6
      t9 = t7 - t6
      a(jp+j) = a(jd+j)
      t10 = c3 * t3 - c2 * t4
      t11 = c2 * t3 + c3 * t4
      a(jd+j) = ax
      u1 = b(ji+j) + b(jt+j)
      u2 = b(jn+j) + b(js+j)
      u3 = b(ji+j) - b(jt+j)
      u4 = b(jn+j) - b(js+j)
      b(jt+j) = b(jx+j)
      u5 = u1 + u2
      u6 = c1 * ( u1 - u2 )
      u7 = b(jp+j) - 0.25 * u5
      bx = b(jp+j) + u5
      u8 = u7 + u6
      u9 = u7 - u6
      b(jp+j) = b(jd+j)
      u10 = c3 * u3 - c2 * u4
      u11 = c2 * u3 + c3 * u4
      b(jd+j) = bx
      a(ji+j) = co1*(t8-u11) - si1*(u8+t11)
      b(ji+j) = si1*(t8-u11) + co1*(u8+t11)
      a(jx+j) = co4*(t8+u11) - si4*(u8-t11)
      b(jx+j) = si4*(t8+u11) + co4*(u8-t11)
      a(jn+j) = co2*(t9-u10) - si2*(u9+t10)
      b(jn+j) = si2*(t9-u10) + co2*(u9+t10)
      a(js+j) = co3*(t9+u10) - si3*(u9-t10)
      b(js+j) = si3*(t9+u10) + co3*(u9-t10)
C----------------------
      t1 = a(jv+j) + a(jy+j)
      t2 = a(jo+j) + a(jt+j)
      t3 = a(jv+j) - a(jy+j)
      t4 = a(jo+j) - a(jt+j)
      a(jv+j) = a(jj+j)
      t5 = t1 + t2
      t6 = c1 * ( t1 - t2 )
      t7 = a(ju+j) - 0.25 * t5
      ax = a(ju+j) + t5
      t8 = t7 + t6
      t9 = t7 - t6
      a(ju+j) = a(je+j)
      t10 = c3 * t3 - c2 * t4
      t11 = c2 * t3 + c3 * t4
      a(je+j) = ax
      u1 = b(jv+j) + b(jy+j)
      u2 = b(jo+j) + b(jt+j)
      u3 = b(jv+j) - b(jy+j)
      u4 = b(jo+j) - b(jt+j)
      b(jv+j) = b(jj+j)
      u5 = u1 + u2
      u6 = c1 * ( u1 - u2 )
      u7 = b(ju+j) - 0.25 * u5
      bx = b(ju+j) + u5
      u8 = u7 + u6
      u9 = u7 - u6
      b(ju+j) = b(je+j)
      u10 = c3 * u3 - c2 * u4
      u11 = c2 * u3 + c3 * u4
      b(je+j) = bx
      a(jj+j) = co1*(t8-u11) - si1*(u8+t11)
      b(jj+j) = si1*(t8-u11) + co1*(u8+t11)
      a(jy+j) = co4*(t8+u11) - si4*(u8-t11)
      b(jy+j) = si4*(t8+u11) + co4*(u8-t11)
      a(jo+j) = co2*(t9-u10) - si2*(u9+t10)
      b(jo+j) = si2*(t9-u10) + co2*(u9+t10)
      a(jt+j) = co3*(t9+u10) - si3*(u9-t10)
      b(jt+j) = si3*(t9+u10) + co3*(u9-t10)
      j = j + jump
  440 continue
C
      endif
C
C-----(end of loop across transforms)
C
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  445 continue
  450 continue
  460 continue
C-----( end of double loop for this k )
      kk = kk + 2*la
  470 continue
C-----( end of loop over values of k )
      la = 5*la
  480 continue
C-----( end of loop on type II radix-5 passes )
C-----( nvex transforms completed)
  490 continue
      istart = istart + nvex * jump
  500 continue
C-----( end of loop on blocks of transforms )
C
      return
      end


