gpfa2f.f

已经编译好的levoo程序

  330 continue
*
      if (n2.eq.8) go to 490
*
*  loop on nonzero k
*  -----------------
      kk = 2 * la
*
      do 350 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)
      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)
*
      do 345 jjj = k , (n-1)*inc , 8*jstep
      ja = istart + jjj
*
*     "transverse" loop
*     -----------------
      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
      aja = a(ja+j)
      aje = a(je+j)
      t0 = aja + aje
      t2 = aja - aje
      ajc = a(jc+j)
      ajg = a(jg+j)
      t1 = ajc + ajg
      t3 = c1 * ( ajc - ajg )
      bja = b(ja+j)
      bje = b(je+j)
      u0 = bja + bje
      u2 = bja - bje
      bjc = b(jc+j)
      bjg = b(jg+j)
      u1 = bjc + bjg
      u3 = c1 * ( bjc - bjg )
      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)
      ajb = a(jb+j)
      ajd = a(jd+j)
      t0 = ajb + ajd
      t2 = ajb - ajd
      ajf = a(jf+j)
      ajh = a(jh+j)
      t1 = ajf - ajh
      t3 = ajf + ajh
      bjb = b(jb+j)
      bjd = b(jd+j)
      u0 = bjb + bjd
      u2 = bjb - bjd
      bjf = b(jf+j)
      bjh = b(jh+j)
      u1 = bjf - bjh
      u3 = bjf + bjh
      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
*
      la = 8 * la
*
*  loop on type II radix-4 passes
*  ------------------------------
  400 continue
      mu = mod(inq,4)
      if (isign.eq.-1) mu = 4 - mu
      ss = 1.0
      if (mu.eq.3) ss = -1.0
*
      do 480 ipass = mh+1 , m
      jstep = (n*inc) / (4*la)
      jstepl = jstep - ninc
      laincl = la * ink - ninc
*
*  k=0 loop (no twiddle factors)
*  -----------------------------
      do 430 ll = 0 , (la-1)*ink , 4*jstep
*
      do 420 jjj = ll , (n-1)*inc , 4*la*ink
      ja = istart + jjj
*
*     "transverse" loop
*     -----------------
      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
*
*  loop across transforms
*  ----------------------
cdir$ ivdep, shortloop
      do 410 l = 1 , nvex
      aja = a(ja+j)
      ajc = a(jc+j)
      t0 = aja + ajc
      t2 = aja - ajc
      ajb = a(jb+j)
      ajd = a(jd+j)
      t1 = ajb + ajd
      t3 = ss * ( ajb - ajd )
      aji = a(ji+j)
      ajc =  aji
      bja = b(ja+j)
      bjc = b(jc+j)
      u0 = bja + bjc
      u2 = bja - bjc
      bjb = b(jb+j)
      bjd = b(jd+j)
      u1 = bjb + bjd
      u3 = ss * ( bjb - bjd )
      aje = a(je+j)
      ajb =  aje
      a(ja+j) = t0 + t1
      a(ji+j) = t0 - t1
      b(ja+j) = u0 + u1
      bjc =  u0 - u1
      bjm = b(jm+j)
      bjd =  bjm
      a(je+j) = t2 - u3
      ajd =  t2 + u3
      bjb =  u2 + t3
      b(jm+j) = u2 - t3
*----------------------
      ajg = a(jg+j)
      t0 = ajb + ajg
      t2 = ajb - ajg
      ajf = a(jf+j)
      ajh = a(jh+j)
      t1 = ajf + ajh
      t3 = ss * ( ajf - ajh )
      ajj = a(jj+j)
      ajg =  ajj
      bje = b(je+j)
      bjg = b(jg+j)
      u0 = bje + bjg
      u2 = bje - bjg
      bjf = b(jf+j)
      bjh = b(jh+j)
      u1 = bjf + bjh
      u3 = ss * ( bjf - bjh )
      b(je+j) = bjb
      a(jb+j) = t0 + t1
      a(jj+j) = t0 - t1
      bjj = b(jj+j)
      bjg =  bjj
      b(jb+j) = u0 + u1
      b(jj+j) = u0 - u1
      a(jf+j) = t2 - u3
      ajh =  t2 + u3
      b(jf+j) = u2 + t3
      bjh =  u2 - t3
*----------------------
      ajk = a(jk+j)
      t0 = ajc + ajk
      t2 = ajc - ajk
      ajl = a(jl+j)
      t1 = ajg + ajl
      t3 = ss * ( ajg - ajl )
      bji = b(ji+j)
      bjk = b(jk+j)
      u0 = bji + bjk
      u2 = bji - bjk
      ajo = a(jo+j)
      ajl =  ajo
      bjl = b(jl+j)
      u1 = bjg + bjl
      u3 = ss * ( bjg - bjl )
      b(ji+j) = bjc
      a(jc+j) = t0 + t1
      a(jk+j) = t0 - t1
      bjo = b(jo+j)
      bjl =  bjo
      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
*----------------------
      ajm = a(jm+j)
      t0 = ajm + ajl
      t2 = ajm - ajl
      ajn = a(jn+j)
      ajp = a(jp+j)
      t1 = ajn + ajp
      t3 = ss * ( ajn - ajp )
      a(jm+j) = ajd
      u0 = bjd + bjl
      u2 = bjd - bjl
      bjn = b(jn+j)
      bjp = b(jp+j)
      u1 = bjn + bjp
      u3 = ss * ( bjn - bjp )
      a(jn+j) = ajh
      a(jd+j) = t0 + t1
      a(jl+j) = t0 - t1
      b(jd+j) = u0 + u1
      b(jl+j) = u0 - u1
      b(jn+j) = bjh
      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
*-----( end of loop across transforms )
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  415 continue
  420 continue
  430 continue
*-----( end of double loop for k=0 )
*
*  finished if last pass
*  ---------------------
      if (ipass.eq.m) go to 490
*
      kk = 2*la
*
*     loop on nonzero k
*     -----------------
      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)
*
*  double loop along first transform in block
*  ------------------------------------------
      do 460 ll = k , (la-1)*ink , 4*jstep
*
      do 450 jjj = ll , (n-1)*inc , 4*la*ink
      ja = istart + jjj
*
*     "transverse" loop
*     -----------------
      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
*
*  loop across transforms
*  ----------------------
cdir$ ivdep, shortloop
      do 440 l = 1 , nvex
      aja = a(ja+j)
      ajc = a(jc+j)
      t0 = aja + ajc
      t2 = aja - ajc
      ajb = a(jb+j)
      ajd = a(jd+j)
      t1 = ajb + ajd
      t3 = ss * ( ajb - ajd )
      aji = a(ji+j)
      ajc =  aji
      bja = b(ja+j)
      bjc = b(jc+j)
      u0 = bja + bjc
      u2 = bja - bjc
      bjb = b(jb+j)
      bjd = b(jd+j)
      u1 = bjb + bjd
      u3 = ss * ( bjb - bjd )
      aje = a(je+j)
      ajb =  aje
      a(ja+j) = t0 + t1
      b(ja+j) = u0 + u1
      a(je+j) = co1*(t2-u3) - si1*(u2+t3)
      bjb =  si1*(t2-u3) + co1*(u2+t3)
      bjm = b(jm+j)
      bjd =  bjm
      a(ji+j) = co2*(t0-t1) - si2*(u0-u1)
      bjc =  si2*(t0-t1) + co2*(u0-u1)
      ajd =  co3*(t2+u3) - si3*(u2-t3)
      b(jm+j) = si3*(t2+u3) + co3*(u2-t3)
*----------------------------------------
      ajg = a(jg+j)
      t0 = ajb + ajg
      t2 = ajb - ajg
      ajf = a(jf+j)
      ajh = a(jh+j)
      t1 = ajf + ajh
      t3 = ss * ( ajf - ajh )
      ajj = a(jj+j)
      ajg =  ajj
      bje = b(je+j)
      bjg = b(jg+j)
      u0 = bje + bjg
      u2 = bje - bjg
      bjf = b(jf+j)
      bjh = b(jh+j)
      u1 = bjf + bjh
      u3 = ss * ( bjf - bjh )
      b(je+j) = bjb
      a(jb+j) = t0 + t1
      b(jb+j) = u0 + u1
      bjj = b(jj+j)
      bjg =  bjj
      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)
      ajh =  co3*(t2+u3) - si3*(u2-t3)
      bjh =  si3*(t2+u3) + co3*(u2-t3)
*----------------------------------------
      ajk = a(jk+j)
      t0 = ajc + ajk
      t2 = ajc - ajk
      ajl = a(jl+j)
      t1 = ajg + ajl
      t3 = ss * ( ajg - ajl )
      bji = b(ji+j)
      bjk = b(jk+j)
      u0 = bji + bjk
      u2 = bji - bjk
      ajo = a(jo+j)
      ajl =  ajo
      bjl = b(jl+j)
      u1 = bjg + bjl
      u3 = ss * ( bjg - bjl )
      b(ji+j) = bjc
      a(jc+j) = t0 + t1
      b(jc+j) = u0 + u1
      bjo = b(jo+j)
      bjl =  bjo
      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)
*----------------------------------------
      ajm = a(jm+j)
      t0 = ajm + ajl
      t2 = ajm - ajl
      ajn = a(jn+j)
      ajp = a(jp+j)
      t1 = ajn + ajp
      t3 = ss * ( ajn - ajp )
      a(jm+j) = ajd
      u0 = bjd + bjl
      u2 = bjd - bjl
      a(jn+j) = ajh
      bjn = b(jn+j)
      bjp = b(jp+j)
      u1 = bjn + bjp
      u3 = ss * ( bjn - bjp )
      b(jn+j) = bjh
      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
*-----(end of loop across transforms)
      ja = ja + jstepx
      if (ja.lt.istart) ja = ja + ninc
  445 continue
  450 continue
  460 continue
*-----( end of double loop for this k )
      kk = kk + 2*la
  470 continue
*-----( end of loop over values of k )
      la = 4*la
  480 continue
*-----( end of loop on type II radix-4 passes )
*-----( nvex transforms completed)
  490 continue
      istart = istart + nvex * jump
  500 continue
*-----( end of loop on blocks of transforms )
*
      return
      end