MODULE FFT_235 ! Code converted using TO_F90 by Alan Miller ! Date: 2003-08-12 Time: 22:29:05 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) CONTAINS ! SUBROUTINE 'SETGPFA' ! SETUP ROUTINE FOR SELF-SORTING IN-PLACE ! GENERALIZED PRIME FACTOR (COMPLEX) FFT [GPFA] ! CALL SETGPFA(TRIGS,N) ! INPUT : ! ----- ! N IS THE LENGTH OF THE TRANSFORMS. N MUST BE OF THE FORM: ! ----------------------------------- ! N = (2**IP) * (3**IQ) * (5**IR) ! ----------------------------------- ! OUTPUT: ! ------ ! TRIGS IS A TABLE OF TWIDDLE FACTORS, ! OF LENGTH 2*IPQR (REAL) WORDS, WHERE: ! -------------------------------------- ! IPQR = (2**IP) + (3**IQ) + (5**IR) ! -------------------------------------- ! WRITTEN BY CLIVE TEMPERTON 1990 !---------------------------------------------------------------------- SUBROUTINE setgpfa(trigs, n) IMPLICIT NONE REAL (dp), INTENT(OUT) :: trigs(:) INTEGER, INTENT(IN) :: n INTEGER :: i, ifac, ip, iq, ir, irot, k, kink, kk, ll, ni, nj(3), nn REAL (dp) :: angle, del, twopi ! DECOMPOSE N INTO FACTORS 2,3,5 ! ------------------------------ nn = n ifac = 2 DO ll = 1, 3 kk = 0 10 IF (MOD(nn,ifac) /= 0) GO TO 20 kk = kk + 1 nn = nn / ifac GO TO 10 20 nj(ll) = kk ifac = ifac + ll END DO IF (nn /= 1) THEN WRITE(6,40) n 40 FORMAT(' *** WARNING!!!', i10, ' IS NOT A LEGAL VALUE OF N ***') RETURN END IF ip = nj(1) iq = nj(2) ir = nj(3) ! COMPUTE LIST OF ROTATED TWIDDLE FACTORS ! --------------------------------------- nj(1) = 2**ip nj(2) = 3**iq nj(3) = 5**ir twopi = 4.0_dp * ASIN(1.0_dp) i = 1 DO ll = 1, 3 ni = nj(ll) IF (ni == 1) CYCLE del = twopi / ni irot = n / ni kink = MOD(irot,ni) kk = 0 DO k = 1, ni angle = kk * del trigs(i) = COS(angle) trigs(i+1) = SIN(angle) i = i + 2 kk = kk + kink IF (kk > ni) kk = kk - ni END DO END DO RETURN END SUBROUTINE setgpfa ! SUBROUTINE 'GPFA' ! SELF-SORTING IN-PLACE GENERALIZED PRIME FACTOR (COMPLEX) FFT ! *** THIS IS THE ALL-FORTRAN VERSION *** ! ------------------------------- ! CALL GPFA(A, B, TRIGS, INC, JUMP, N, LOT, ISIGN) ! A IS FIRST REAL INPUT/OUTPUT VECTOR ! B IS FIRST IMAGINARY INPUT/OUTPUT VECTOR ! TRIGS IS A TABLE OF TWIDDLE FACTORS, PRECALCULATED ! BY CALLING SUBROUTINE 'SETGPFA' ! INC IS THE INCREMENT WITHIN EACH DATA VECTOR ! JUMP IS THE INCREMENT BETWEEN DATA VECTORS ! N IS THE LENGTH OF THE TRANSFORMS: ! ----------------------------------- ! N = (2**IP) * (3**IQ) * (5**IR) ! ----------------------------------- ! LOT IS THE NUMBER OF TRANSFORMS ! ISIGN = +1 FOR FORWARD TRANSFORM ! = -1 FOR INVERSE TRANSFORM ! WRITTEN BY CLIVE TEMPERTON ! RECHERCHE EN PREVISION NUMERIQUE ! ATMOSPHERIC ENVIRONMENT SERVICE, CANADA !---------------------------------------------------------------------- ! DEFINITION OF TRANSFORM ! ----------------------- ! X(J) = SUM(K=0,...,N-1)(C(K)*EXP(ISIGN*2*I*J*K*PI/N)) !--------------------------------------------------------------------- ! FOR A MATHEMATICAL DEVELOPMENT OF THE ALGORITHM USED, ! SEE: ! C TEMPERTON : "A GENERALIZED PRIME FACTOR FFT ALGORITHM ! FOR ANY N = (2**P)(3**Q)(5**R)", ! SIAM J. SCI. STAT. COMP., MAY 1992. !---------------------------------------------------------------------- SUBROUTINE gpfa(a, b, trigs, inc, jump, n, lot, ISIGN) IMPLICIT NONE REAL (dp), INTENT(IN OUT) :: a(:) REAL (dp), INTENT(IN OUT) :: b(:) REAL (dp), INTENT(IN) :: trigs(:) INTEGER, INTENT(IN) :: inc INTEGER, INTENT(IN) :: jump INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN OUT) :: lot INTEGER, INTENT(IN) :: ISIGN INTEGER :: i, ifac, ip, iq, ir, kk, ll, nj(3), nn ! DECOMPOSE N INTO FACTORS 2,3,5 ! ------------------------------ nn = n ifac = 2 DO ll = 1, 3 kk = 0 10 IF (MOD(nn,ifac) /= 0) GO TO 20 kk = kk + 1 nn = nn / ifac GO TO 10 20 nj(ll) = kk ifac = ifac + ll END DO IF (nn /= 1) THEN WRITE(6,40) n 40 FORMAT(' *** WARNING!!!', i10, ' IS NOT A LEGAL VALUE OF N ***') RETURN END IF ip = nj(1) iq = nj(2) ir = nj(3) ! COMPUTE THE TRANSFORM ! --------------------- i = 1 IF (ip > 0) THEN CALL gpfa2f(a, b, trigs, inc, jump, n, ip, lot, ISIGN) i = i + 2 * ( 2**ip) END IF IF (iq > 0) THEN CALL gpfa3f(a, b, trigs(i), inc, jump, n, iq, lot, ISIGN) i = i + 2 * (3**iq) END IF IF (ir > 0) THEN CALL gpfa5f(a, b, trigs(i), inc, jump, n, ir, lot, ISIGN) END IF RETURN END SUBROUTINE gpfa ! fortran version of *gpfa2* - ! radix-2 section of self-sorting, in-place, generalized pfa ! central radix-2 and radix-8 passes included ! so that transform length can be any power of 2 !------------------------------------------------------------------- SUBROUTINE gpfa2f(a, b, trigs, inc, jump, n, mm, lot, ISIGN) IMPLICIT NONE REAL (dp), INTENT(IN OUT) :: a(*) REAL (dp), INTENT(IN OUT) :: b(*) REAL (dp), INTENT(IN) :: trigs(*) INTEGER, INTENT(IN) :: inc INTEGER, INTENT(IN) :: jump INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: mm INTEGER, INTENT(IN) :: lot INTEGER, INTENT(IN) :: ISIGN INTEGER :: lvr = 1024 INTEGER :: ink, inq, ipass, istart, & j, ja, jb, jc, jd, je, jf, jg, jh, ji, jj, jjj, & jk, jl, jm, jn, jo, jp, jstep, jstepl, jstepx, & k, kk, l, la, laincl, left, ll, & m, m2, m8, mh, mu, n2, nb, nblox, ninc, nu, nvex REAL (dp) :: aja, ajb, ajc, ajd, aje, ajf, ajg, ajh, aji, ajj, & ajk, ajl, ajm, ajn, ajo, ajp, & bja, bjb, bjc, bjd, bjf, bjg, bje, bjh, bji, bjj, & bjk, bjl, bjm, bjn, bjo, bjp, & c1, c2, c3, co1, co2, co3, co4, co5, co6, co7, & s, si1, si2, si3, si4, si5, si6, si7, ss, t0, t1, t2, t3, & u0, u1, u2, u3 ! *************************************************************** ! * * ! * N.B. LVR = LENGTH OF VECTOR REGISTERS, SET TO 128 FOR C90. * ! * RESET TO 64 FOR OTHER CRAY MACHINES, OR TO ANY LARGE VALUE * ! * (GREATER THAN OR EQUAL TO LOT) FOR A SCALAR COMPUTER. * ! * * ! *************************************************************** n2 = 2**mm inq = n/n2 jstepx = (n2-n) * inc ninc = n * inc ink = inc * inq m2 = 0 m8 = 0 IF (MOD(mm,2) == 0) THEN m = mm/2 ELSE IF (MOD(mm,4) == 1) THEN m = (mm-1)/2 m2 = 1 ELSE IF (MOD(mm,4) == 3) THEN m = (mm-3)/2 m8 = 1 END IF mh = (m+1)/2 nblox = 1 + (lot-1)/lvr left = lot s = ISIGN istart = 1 ! loop on blocks of lvr transforms ! -------------------------------- DO nb = 1, nblox IF (left <= lvr) THEN nvex = left ELSE IF (left < (2*lvr)) THEN nvex = left/2 nvex = nvex + MOD(nvex,2) ELSE nvex = lvr END IF left = left - nvex la = 1 ! loop on type I radix-4 passes ! ----------------------------- mu = MOD(inq,4) IF (ISIGN == -1) mu = 4 - mu ss = 1.0_dp IF (mu == 3) ss = -1.0_dp IF (mh == 0) GO TO 200 DO ipass = 1, mh jstep = (n*inc) / (4*la) jstepl = jstep - ninc ! k = 0 loop (no twiddle factors) ! ------------------------------- DO jjj = 0, (n-1)*inc, 4*jstep ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc jc = jb + jstepl IF (jc < istart) jc = jc + ninc jd = jc + jstepl IF (jd < istart) jd = jd + ninc j = 0 ! loop across transforms ! ---------------------- ! DO 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 ) 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 ) 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 END DO ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO ! finished if n2 = 4 ! ------------------ IF (n2 == 4) GO TO 490 kk = 2 * la ! loop on nonzero k ! ----------------- DO 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) ! loop along transform ! -------------------- DO jjj = k, (n-1)*inc, 4*jstep ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc jc = jb + jstepl IF (jc < istart) jc = jc + ninc jd = jc + jstepl IF (jd < istart) jd = jd + ninc j = 0 ! loop across transforms ! ---------------------- ! DO 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 ) 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 ) 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 END DO !-----( end of loop across transforms ) ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO !-----( end of loop along transforms ) kk = kk + 2*la END DO !-----( end of loop on nonzero k ) la = 4*la END DO !-----( end of loop on type I radix-4 passes) ! central radix-2 pass ! -------------------- 200 IF (m2 == 0) GO TO 300 jstep = (n*inc) / (2*la) jstepl = jstep - ninc ! k=0 loop (no twiddle factors) ! ----------------------------- DO jjj = 0, (n-1)*inc, 2*jstep ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc j = 0 ! loop across transforms ! ---------------------- ! DO l = 1, nvex aja = a(ja+j) ajb = a(jb+j) t0 = aja - ajb a(ja+j) = aja + ajb a(jb+j) = t0 bja = b(ja+j) bjb = b(jb+j) u0 = bja - bjb b(ja+j) = bja + bjb b(jb+j) = u0 j = j + jump END DO !-----(end of loop across transforms) ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO ! finished if n2=2 ! ---------------- IF (n2 == 2) GO TO 490 kk = 2 * la ! loop on nonzero k ! ----------------- DO k = ink, jstep - ink, ink co1 = trigs(kk+1) si1 = s*trigs(kk+2) ! loop along transforms ! --------------------- DO jjj = k, (n-1)*inc, 2*jstep ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc j = 0 ! loop across transforms ! ---------------------- IF (kk == n2/2) THEN ! DO l = 1, nvex aja = a(ja+j) ajb = a(jb+j) t0 = ss * ( aja - ajb ) a(ja+j) = aja + ajb bjb = b(jb+j) bja = b(ja+j) a(jb+j) = ss * ( bjb - bja ) b(ja+j) = bja + bjb b(jb+j) = t0 j = j + jump END DO ELSE DO l = 1, nvex aja = a(ja+j) ajb = a(jb+j) t0 = aja - ajb a(ja+j) = aja + ajb bja = b(ja+j) bjb = b(jb+j) u0 = bja - bjb b(ja+j) = bja + bjb a(jb+j) = co1*t0 - si1*u0 b(jb+j) = si1*t0 + co1*u0 j = j + jump END DO END IF !-----(end of loop across transforms) ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO !-----(end of loop along transforms) kk = kk + 2 * la END DO !-----(end of loop on nonzero k) !-----(end of radix-2 pass) la = 2 * la GO TO 400 ! central radix-8 pass ! -------------------- 300 IF (m8 == 0) GO TO 400 jstep = (n*inc) / (8*la) jstepl = jstep - ninc mu = MOD(inq,8) IF (ISIGN == -1) mu = 8 - mu c1 = 1.0_dp IF (mu == 3 .OR. mu == 7) c1 = -1.0_dp c2 = SQRT(0.5_dp) IF (mu == 3 .OR. mu == 5) c2 = -c2 c3 = c1 * c2 ! stage 1 ! ------- DO k = 0, jstep - ink, ink DO jjj = k, (n-1)*inc, 8*jstep ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc jc = jb + jstepl IF (jc < istart) jc = jc + ninc jd = jc + jstepl IF (jd < istart) jd = jd + ninc je = jd + jstepl IF (je < istart) je = je + ninc jf = je + jstepl IF (jf < istart) jf = jf + ninc jg = jf + jstepl IF (jg < istart) jg = jg + ninc jh = jg + jstepl IF (jh < istart) jh = jh + ninc j = 0 ! DO l = 1, nvex aja = a(ja+j) aje = a(je+j) t0 = aja - aje a(ja+j) = aja + aje ajc = a(jc+j) ajg = a(jg+j) t1 = c1 * ( ajc - ajg ) a(je+j) = ajc + ajg ajb = a(jb+j) ajf = a(jf+j) t2 = ajb - ajf a(jc+j) = ajb + ajf ajd = a(jd+j) ajh = a(jh+j) t3 = ajd - ajh a(jg+j) = ajd + ajh a(jb+j) = t0 a(jf+j) = t1 a(jd+j) = c2 * ( t2 - t3 ) a(jh+j) = c3 * ( t2 + t3 ) bja = b(ja+j) bje = b(je+j) u0 = bja - bje b(ja+j) = bja + bje bjc = b(jc+j) bjg = b(jg+j) u1 = c1 * ( bjc - bjg ) b(je+j) = bjc + bjg bjb = b(jb+j) bjf = b(jf+j) u2 = bjb - bjf b(jc+j) = bjb + bjf bjd = b(jd+j) bjh = b(jh+j) u3 = bjd - bjh b(jg+j) = bjd + bjh b(jb+j) = u0 b(jf+j) = u1 b(jd+j) = c2 * ( u2 - u3 ) b(jh+j) = c3 * ( u2 + u3 ) j = j + jump END DO ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO END DO ! stage 2 ! ------- ! k=0 (no twiddle factors) ! ------------------------ DO jjj = 0, (n-1)*inc, 8*jstep ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc jc = jb + jstepl IF (jc < istart) jc = jc + ninc jd = jc + jstepl IF (jd < istart) jd = jd + ninc je = jd + jstepl IF (je < istart) je = je + ninc jf = je + jstepl IF (jf < istart) jf = jf + ninc jg = jf + jstepl IF (jg < istart) jg = jg + ninc jh = jg + jstepl IF (jh < istart) jh = jh + ninc j = 0 ! DO 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 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 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) = 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 END DO ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO IF (n2 == 8) GO TO 490 ! loop on nonzero k ! ----------------- kk = 2 * la DO 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 jjj = k, (n-1)*inc, 8*jstep ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc jc = jb + jstepl IF (jc < istart) jc = jc + ninc jd = jc + jstepl IF (jd < istart) jd = jd + ninc je = jd + jstepl IF (je < istart) je = je + ninc jf = je + jstepl IF (jf < istart) jf = jf + ninc jg = jf + jstepl IF (jg < istart) jg = jg + ninc jh = jg + jstepl IF (jh < istart) jh = jh + ninc j = 0 ! DO 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 END DO ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO kk = kk + 2 * la END DO la = 8 * la ! loop on type II radix-4 passes ! ------------------------------ 400 mu = MOD(inq,4) IF (ISIGN == -1) mu = 4 - mu ss = 1.0_dp IF (mu == 3) ss = -1.0_dp DO ipass = mh+1, m jstep = (n*inc) / (4*la) jstepl = jstep - ninc laincl = la * ink - ninc ! k=0 loop (no twiddle factors) ! ----------------------------- DO ll = 0, (la-1)*ink, 4*jstep DO jjj = ll, (n-1)*inc, 4*la*ink ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc jc = jb + jstepl IF (jc < istart) jc = jc + ninc jd = jc + jstepl IF (jd < istart) jd = jd + ninc je = ja + laincl IF (je < istart) je = je + ninc jf = je + jstepl IF (jf < istart) jf = jf + ninc jg = jf + jstepl IF (jg < istart) jg = jg + ninc jh = jg + jstepl IF (jh < istart) jh = jh + ninc ji = je + laincl IF (ji < istart) ji = ji + ninc jj = ji + jstepl IF (jj < istart) jj = jj + ninc jk = jj + jstepl IF (jk < istart) jk = jk + ninc jl = jk + jstepl IF (jl < istart) jl = jl + ninc jm = ji + laincl IF (jm < istart) jm = jm + ninc jn = jm + jstepl IF (jn < istart) jn = jn + ninc jo = jn + jstepl IF (jo < istart) jo = jo + ninc jp = jo + jstepl IF (jp < istart) jp = jp + ninc j = 0 ! loop across transforms ! ---------------------- ! DO 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 END DO !-----( end of loop across transforms ) ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO END DO !-----( end of double loop for k=0 ) ! finished if last pass ! --------------------- IF (ipass == m) EXIT kk = 2*la ! loop on nonzero k ! ----------------- DO 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 ll = k, (la-1)*ink, 4*jstep DO jjj = ll, (n-1)*inc, 4*la*ink ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc jc = jb + jstepl IF (jc < istart) jc = jc + ninc jd = jc + jstepl IF (jd < istart) jd = jd + ninc je = ja + laincl IF (je < istart) je = je + ninc jf = je + jstepl IF (jf < istart) jf = jf + ninc jg = jf + jstepl IF (jg < istart) jg = jg + ninc jh = jg + jstepl IF (jh < istart) jh = jh + ninc ji = je + laincl IF (ji < istart) ji = ji + ninc jj = ji + jstepl IF (jj < istart) jj = jj + ninc jk = jj + jstepl IF (jk < istart) jk = jk + ninc jl = jk + jstepl IF (jl < istart) jl = jl + ninc jm = ji + laincl IF (jm < istart) jm = jm + ninc jn = jm + jstepl IF (jn < istart) jn = jn + ninc jo = jn + jstepl IF (jo < istart) jo = jo + ninc jp = jo + jstepl IF (jp < istart) jp = jp + ninc j = 0 ! loop across transforms ! ---------------------- ! DO 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 END DO !-----(end of loop across transforms) ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO END DO !-----( end of double loop for this k ) kk = kk + 2*la END DO !-----( end of loop over values of k ) la = 4*la END DO !-----( end of loop on type II radix-4 passes ) !-----( nvex transforms completed) 490 istart = istart + nvex * jump END DO !-----( end of loop on blocks of transforms ) RETURN END SUBROUTINE gpfa2f ! fortran version of *gpfa3* - ! radix-3 section of self-sorting, in-place ! generalized PFA !------------------------------------------------------------------- SUBROUTINE gpfa3f(a, b, trigs, inc, jump, n, mm, lot, ISIGN) REAL (dp), INTENT(IN OUT) :: a(*) REAL (dp), INTENT(IN OUT) :: b(*) REAL (dp), INTENT(IN) :: trigs(*) INTEGER, INTENT(IN) :: inc INTEGER, INTENT(IN) :: jump INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: mm INTEGER, INTENT(IN) :: lot INTEGER, INTENT(IN) :: ISIGN REAL (dp), PARAMETER :: sin60 = 0.866025403784437_dp INTEGER, PARAMETER :: lvr = 128 ! *************************************************************** ! * * ! * N.B. LVR = LENGTH OF VECTOR REGISTERS, SET TO 128 FOR C90. * ! * RESET TO 64 FOR OTHER CRAY MACHINES, OR TO ANY LARGE VALUE * ! * (GREATER THAN OR EQUAL TO LOT) FOR A SCALAR COMPUTER. * ! * * ! *************************************************************** INTEGER :: ink, inq, ipass, istart, & j, ja, jb, jc, jd, je, jf, jg, jh, ji, jjj, jstep, jstepl, jstepx, & k, kk, l, la, laincl, left, ll, & m, mh, mu, n3, nb, nblox, ninc, nu, nvex REAL (dp) :: aja, ajb, ajc, ajd, aje, ajf, ajg, ajh, aji, & bja, bjb, bjc, bjd, bje, bjf, bjg, bjh, bji, c1, co1, co2, & s, si1, si2, t1, t2, t3, u1, u2, u3 n3 = 3**mm inq = n/n3 jstepx = (n3-n) * inc ninc = n * inc ink = inc * inq mu = MOD(inq,3) IF (ISIGN == -1) mu = 3-mu m = mm mh = (m+1)/2 s = ISIGN c1 = sin60 IF (mu == 2) c1 = -c1 nblox = 1 + (lot-1)/lvr left = lot s = ISIGN istart = 1 ! loop on blocks of lvr transforms ! -------------------------------- DO nb = 1, nblox IF (left <= lvr) THEN nvex = left ELSE IF (left < (2*lvr)) THEN nvex = left/2 nvex = nvex + MOD(nvex,2) ELSE nvex = lvr END IF left = left - nvex la = 1 ! loop on type I radix-3 passes ! ----------------------------- DO ipass = 1, mh jstep = (n*inc) / (3*la) jstepl = jstep - ninc ! k = 0 loop (no twiddle factors) ! ------------------------------- DO jjj = 0, (n-1)*inc, 3*jstep ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc jc = jb + jstepl IF (jc < istart) jc = jc + ninc j = 0 ! loop across transforms ! ---------------------- ! DO l = 1, nvex ajb = a(jb+j) ajc = a(jc+j) t1 = ajb + ajc aja = a(ja+j) t2 = aja - 0.5_dp * t1 t3 = c1 * ( ajb - ajc ) bjb = b(jb+j) bjc = b(jc+j) u1 = bjb + bjc bja = b(ja+j) u2 = bja - 0.5_dp * u1 u3 = c1 * ( bjb - bjc ) a(ja+j) = aja + t1 b(ja+j) = bja + u1 a(jb+j) = t2 - u3 b(jb+j) = u2 + t3 a(jc+j) = t2 + u3 b(jc+j) = u2 - t3 j = j + jump END DO ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO ! finished if n3 = 3 ! ------------------ IF (n3 == 3) GO TO 490 kk = 2 * la ! loop on nonzero k ! ----------------- DO 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) ! loop along transform ! -------------------- DO jjj = k, (n-1)*inc, 3*jstep ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc jc = jb + jstepl IF (jc < istart) jc = jc + ninc j = 0 ! loop across transforms ! ---------------------- ! DO l = 1, nvex ajb = a(jb+j) ajc = a(jc+j) t1 = ajb + ajc aja = a(ja+j) t2 = aja - 0.5_dp * t1 t3 = c1 * ( ajb - ajc ) bjb = b(jb+j) bjc = b(jc+j) u1 = bjb + bjc bja = b(ja+j) u2 = bja - 0.5_dp * u1 u3 = c1 * ( bjb - bjc ) a(ja+j) = aja + t1 b(ja+j) = bja + 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 END DO !-----( end of loop across transforms ) ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO !-----( end of loop along transforms ) kk = kk + 2*la END DO !-----( end of loop on nonzero k ) la = 3*la END DO !-----( end of loop on type I radix-3 passes) ! loop on type II radix-3 passes ! ------------------------------ DO ipass = mh+1, m jstep = (n*inc) / (3*la) jstepl = jstep - ninc laincl = la*ink - ninc ! k=0 loop (no twiddle factors) ! ----------------------------- DO ll = 0, (la-1)*ink, 3*jstep DO jjj = ll, (n-1)*inc, 3*la*ink ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc jc = jb + jstepl IF (jc < istart) jc = jc + ninc jd = ja + laincl IF (jd < istart) jd = jd + ninc je = jd + jstepl IF (je < istart) je = je + ninc jf = je + jstepl IF (jf < istart) jf = jf + ninc jg = jd + laincl IF (jg < istart) jg = jg + ninc jh = jg + jstepl IF (jh < istart) jh = jh + ninc ji = jh + jstepl IF (ji < istart) ji = ji + ninc j = 0 ! loop across transforms ! ---------------------- ! DO l = 1, nvex ajb = a(jb+j) ajc = a(jc+j) t1 = ajb + ajc aja = a(ja+j) t2 = aja - 0.5_dp * t1 t3 = c1 * ( ajb - ajc ) ajd = a(jd+j) ajb = ajd bjb = b(jb+j) bjc = b(jc+j) u1 = bjb + bjc bja = b(ja+j) u2 = bja - 0.5_dp * u1 u3 = c1 * ( bjb - bjc ) bjd = b(jd+j) bjb = bjd a(ja+j) = aja + t1 b(ja+j) = bja + u1 a(jd+j) = t2 - u3 b(jd+j) = u2 + t3 ajc = t2 + u3 bjc = u2 - t3 !---------------------- aje = a(je+j) ajf = a(jf+j) t1 = aje + ajf t2 = ajb - 0.5_dp * t1 t3 = c1 * ( aje - ajf ) ajh = a(jh+j) ajf = ajh bje = b(je+j) bjf = b(jf+j) u1 = bje + bjf u2 = bjb - 0.5_dp * u1 u3 = c1 * ( bje - bjf ) bjh = b(jh+j) bjf = bjh a(jb+j) = ajb + t1 b(jb+j) = bjb + u1 a(je+j) = t2 - u3 b(je+j) = u2 + t3 a(jh+j) = t2 + u3 b(jh+j) = u2 - t3 !---------------------- aji = a(ji+j) t1 = ajf + aji ajg = a(jg+j) t2 = ajg - 0.5_dp * t1 t3 = c1 * ( ajf - aji ) t1 = ajg + t1 a(jg+j) = ajc bji = b(ji+j) u1 = bjf + bji bjg = b(jg+j) u2 = bjg - 0.5_dp * u1 u3 = c1 * ( bjf - bji ) u1 = bjg + u1 b(jg+j) = bjc 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 END DO !-----( end of loop across transforms ) ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO END DO !-----( end of double loop for k=0 ) ! finished if last pass ! --------------------- IF (ipass == m) EXIT kk = 2*la ! loop on nonzero k ! ----------------- DO 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) ! double loop along first transform in block ! ------------------------------------------ DO ll = k, (la-1)*ink, 3*jstep DO jjj = ll, (n-1)*inc, 3*la*ink ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc jc = jb + jstepl IF (jc < istart) jc = jc + ninc jd = ja + laincl IF (jd < istart) jd = jd + ninc je = jd + jstepl IF (je < istart) je = je + ninc jf = je + jstepl IF (jf < istart) jf = jf + ninc jg = jd + laincl IF (jg < istart) jg = jg + ninc jh = jg + jstepl IF (jh < istart) jh = jh + ninc ji = jh + jstepl IF (ji < istart) ji = ji + ninc j = 0 ! loop across transforms ! ---------------------- ! DO l = 1, nvex ajb = a(jb+j) ajc = a(jc+j) t1 = ajb + ajc aja = a(ja+j) t2 = aja - 0.5_dp * t1 t3 = c1 * ( ajb - ajc ) ajd = a(jd+j) ajb = ajd bjb = b(jb+j) bjc = b(jc+j) u1 = bjb + bjc bja = b(ja+j) u2 = bja - 0.5_dp * u1 u3 = c1 * ( bjb - bjc ) bjd = b(jd+j) bjb = bjd a(ja+j) = aja + t1 b(ja+j) = bja + u1 a(jd+j) = co1*(t2-u3) - si1*(u2+t3) b(jd+j) = si1*(t2-u3) + co1*(u2+t3) ajc = co2*(t2+u3) - si2*(u2-t3) bjc = si2*(t2+u3) + co2*(u2-t3) !---------------------- aje = a(je+j) ajf = a(jf+j) t1 = aje + ajf t2 = ajb - 0.5_dp * t1 t3 = c1 * ( aje - ajf ) ajh = a(jh+j) ajf = ajh bje = b(je+j) bjf = b(jf+j) u1 = bje + bjf u2 = bjb - 0.5_dp * u1 u3 = c1 * ( bje - bjf ) bjh = b(jh+j) bjf = bjh a(jb+j) = ajb + t1 b(jb+j) = bjb + 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) !---------------------- aji = a(ji+j) t1 = ajf + aji ajg = a(jg+j) t2 = ajg - 0.5_dp * t1 t3 = c1 * ( ajf - aji ) t1 = ajg + t1 a(jg+j) = ajc bji = b(ji+j) u1 = bjf + bji bjg = b(jg+j) u2 = bjg - 0.5_dp * u1 u3 = c1 * ( bjf - bji ) u1 = bjg + u1 b(jg+j) = bjc 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 END DO !-----(end of loop across transforms) ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO END DO !-----( end of double loop for this k ) kk = kk + 2*la END DO !-----( end of loop over values of k ) la = 3*la END DO !-----( end of loop on type II radix-3 passes ) !-----( nvex transforms completed) 490 istart = istart + nvex * jump END DO !-----( end of loop on blocks of transforms ) RETURN END SUBROUTINE gpfa3f ! fortran version of *gpfa5* - ! radix-5 section of self-sorting, in-place, ! generalized pfa !------------------------------------------------------------------- SUBROUTINE gpfa5f(a, b, trigs, inc, jump, n, mm, lot, ISIGN) IMPLICIT NONE REAL (dp), INTENT(IN OUT) :: a(*) REAL (dp), INTENT(IN OUT) :: b(*) REAL (dp), INTENT(IN) :: trigs(*) INTEGER, INTENT(IN) :: inc INTEGER, INTENT(IN) :: jump INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: mm INTEGER, INTENT(IN) :: lot INTEGER, INTENT(IN) :: ISIGN REAL (dp), PARAMETER :: sin36 = 0.587785252292473_dp, & sin72 = 0.951056516295154_dp, & qrt5 = 0.559016994374947_dp INTEGER, PARAMETER :: lvr = 128 ! *************************************************************** ! * * ! * N.B. LVR = LENGTH OF VECTOR REGISTERS, SET TO 128 FOR C90. * ! * RESET TO 64 FOR OTHER CRAY MACHINES, OR TO ANY LARGE VALUE * ! * (GREATER THAN OR EQUAL TO LOT) FOR A SCALAR COMPUTER. * ! * * ! *************************************************************** INTEGER :: ink, inq, ipass, istart, & j, ja, jb, jc, jd, je, jf, jg, jh, ji, jj, jjj, jk, jl, jm, jn, & jo, jp, jq, jr, js, jstep, jstepl, jstepx, jt, ju, jv, jw, jx, & jy, k, kk, l, la, laincl, left, ll, m, mh, mu, & n5, nb, nblox, ninc, nu, nvex REAL (dp) :: aja, ajb, ajc, ajd, aje, ajf, ajg, ajh, aji, ajj, ajk, ajl, & ajm, ajn, & ajo, ajp, ajq, ajr, ajs, ajt, aju, ajv, ajw, ajx, ajy, ax, & bja, bjb, bjc, bjd, bje, bjg, bjf, bjh, bji, bjj, bjk, bjl, & bjm, bjn, bjo, bjp, bjq, bjr, bjs, bjt, bju, bjv, bjw, bjx, & bjy, bx, & c1, c2, c3, co1, co2, co3, co4, s, si1, si2, si3, si4, & t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, & u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11 n5 = 5 ** mm inq = n / n5 jstepx = (n5-n) * inc ninc = n * inc ink = inc * inq mu = MOD(inq,5) IF (ISIGN == -1) mu = 5 - mu m = mm mh = (m+1)/2 s = ISIGN c1 = qrt5 c2 = sin72 c3 = sin36 IF (mu == 2.OR.mu == 3) THEN c1 = -c1 c2 = sin36 c3 = sin72 END IF IF (mu == 3.OR.mu == 4) c2 = -c2 IF (mu == 2.OR.mu == 4) c3 = -c3 nblox = 1 + (lot-1)/lvr left = lot s = ISIGN istart = 1 ! loop on blocks of lvr transforms ! -------------------------------- DO nb = 1, nblox IF (left <= lvr) THEN nvex = left ELSE IF (left < (2*lvr)) THEN nvex = left/2 nvex = nvex + MOD(nvex,2) ELSE nvex = lvr END IF left = left - nvex la = 1 ! loop on type I radix-5 passes ! ----------------------------- DO ipass = 1, mh jstep = (n*inc) / (5*la) jstepl = jstep - ninc kk = 0 ! loop on k ! --------- DO k = 0, jstep-ink, ink IF (k > 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) END IF ! loop along transform ! -------------------- DO jjj = k, (n-1)*inc, 5*jstep ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc jc = jb + jstepl IF (jc < istart) jc = jc + ninc jd = jc + jstepl IF (jd < istart) jd = jd + ninc je = jd + jstepl IF (je < istart) je = je + ninc j = 0 ! loop across transforms ! ---------------------- IF (k == 0) THEN ! DO l = 1, nvex ajb = a(jb+j) aje = a(je+j) t1 = ajb + aje ajc = a(jc+j) ajd = a(jd+j) t2 = ajc + ajd t3 = ajb - aje t4 = ajc - ajd t5 = t1 + t2 t6 = c1 * ( t1 - t2 ) aja = a(ja+j) t7 = aja - 0.25_dp * t5 a(ja+j) = aja + t5 t8 = t7 + t6 t9 = t7 - t6 t10 = c3 * t3 - c2 * t4 t11 = c2 * t3 + c3 * t4 bjb = b(jb+j) bje = b(je+j) u1 = bjb + bje bjc = b(jc+j) bjd = b(jd+j) u2 = bjc + bjd u3 = bjb - bje u4 = bjc - bjd u5 = u1 + u2 u6 = c1 * ( u1 - u2 ) bja = b(ja+j) u7 = bja - 0.25_dp * u5 b(ja+j) = bja + 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 END DO ELSE DO l = 1, nvex ajb = a(jb+j) aje = a(je+j) t1 = ajb + aje ajc = a(jc+j) ajd = a(jd+j) t2 = ajc + ajd t3 = ajb - aje t4 = ajc - ajd t5 = t1 + t2 t6 = c1 * ( t1 - t2 ) aja = a(ja+j) t7 = aja - 0.25_dp * t5 a(ja+j) = aja + t5 t8 = t7 + t6 t9 = t7 - t6 t10 = c3 * t3 - c2 * t4 t11 = c2 * t3 + c3 * t4 bjb = b(jb+j) bje = b(je+j) u1 = bjb + bje bjc = b(jc+j) bjd = b(jd+j) u2 = bjc + bjd u3 = bjb - bje u4 = bjc - bjd u5 = u1 + u2 u6 = c1 * ( u1 - u2 ) bja = b(ja+j) u7 = bja - 0.25_dp * u5 b(ja+j) = bja + 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 END DO END IF !-----( end of loop across transforms ) ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO !-----( end of loop along transforms ) kk = kk + 2*la END DO !-----( end of loop on nonzero k ) la = 5*la END DO !-----( end of loop on type I radix-5 passes) IF (n == 5) GO TO 490 ! loop on type II radix-5 passes ! ------------------------------ DO ipass = mh+1, m jstep = (n*inc) / (5*la) jstepl = jstep - ninc laincl = la * ink - ninc kk = 0 ! loop on k ! --------- DO k = 0, jstep-ink, ink IF (k > 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) END IF ! double loop along first transform in block ! ------------------------------------------ DO ll = k, (la-1)*ink, 5*jstep DO jjj = ll, (n-1)*inc, 5*la*ink ja = istart + jjj ! "transverse" loop ! ----------------- DO nu = 1, inq jb = ja + jstepl IF (jb < istart) jb = jb + ninc jc = jb + jstepl IF (jc < istart) jc = jc + ninc jd = jc + jstepl IF (jd < istart) jd = jd + ninc je = jd + jstepl IF (je < istart) je = je + ninc jf = ja + laincl IF (jf < istart) jf = jf + ninc jg = jf + jstepl IF (jg < istart) jg = jg + ninc jh = jg + jstepl IF (jh < istart) jh = jh + ninc ji = jh + jstepl IF (ji < istart) ji = ji + ninc jj = ji + jstepl IF (jj < istart) jj = jj + ninc jk = jf + laincl IF (jk < istart) jk = jk + ninc jl = jk + jstepl IF (jl < istart) jl = jl + ninc jm = jl + jstepl IF (jm < istart) jm = jm + ninc jn = jm + jstepl IF (jn < istart) jn = jn + ninc jo = jn + jstepl IF (jo < istart) jo = jo + ninc jp = jk + laincl IF (jp < istart) jp = jp + ninc jq = jp + jstepl IF (jq < istart) jq = jq + ninc jr = jq + jstepl IF (jr < istart) jr = jr + ninc js = jr + jstepl IF (js < istart) js = js + ninc jt = js + jstepl IF (jt < istart) jt = jt + ninc ju = jp + laincl IF (ju < istart) ju = ju + ninc jv = ju + jstepl IF (jv < istart) jv = jv + ninc jw = jv + jstepl IF (jw < istart) jw = jw + ninc jx = jw + jstepl IF (jx < istart) jx = jx + ninc jy = jx + jstepl IF (jy < istart) jy = jy + ninc j = 0 ! loop across transforms ! ---------------------- IF (k == 0) THEN DO l = 1, nvex ajb = a(jb+j) aje = a(je+j) t1 = ajb + aje ajc = a(jc+j) ajd = a(jd+j) t2 = ajc + ajd t3 = ajb - aje t4 = ajc - ajd ajf = a(jf+j) ajb = ajf t5 = t1 + t2 t6 = c1 * ( t1 - t2 ) aja = a(ja+j) t7 = aja - 0.25_dp * t5 a(ja+j) = aja + t5 t8 = t7 + t6 t9 = t7 - t6 ajk = a(jk+j) ajc = ajk t10 = c3 * t3 - c2 * t4 t11 = c2 * t3 + c3 * t4 bjb = b(jb+j) bje = b(je+j) u1 = bjb + bje bjc = b(jc+j) bjd = b(jd+j) u2 = bjc + bjd u3 = bjb - bje u4 = bjc - bjd bjf = b(jf+j) bjb = bjf u5 = u1 + u2 u6 = c1 * ( u1 - u2 ) bja = b(ja+j) u7 = bja - 0.25_dp * u5 b(ja+j) = bja + u5 u8 = u7 + u6 u9 = u7 - u6 bjk = b(jk+j) bjc = bjk u10 = c3 * u3 - c2 * u4 u11 = c2 * u3 + c3 * u4 a(jf+j) = t8 - u11 b(jf+j) = u8 + t11 aje = t8 + u11 bje = u8 - t11 a(jk+j) = t9 - u10 b(jk+j) = u9 + t10 ajd = t9 + u10 bjd = u9 - t10 !---------------------- ajg = a(jg+j) ajj = a(jj+j) t1 = ajg + ajj ajh = a(jh+j) aji = a(ji+j) t2 = ajh + aji t3 = ajg - ajj t4 = ajh - aji ajl = a(jl+j) ajh = ajl t5 = t1 + t2 t6 = c1 * ( t1 - t2 ) t7 = ajb - 0.25_dp * t5 a(jb+j) = ajb + t5 t8 = t7 + t6 t9 = t7 - t6 ajq = a(jq+j) aji = ajq t10 = c3 * t3 - c2 * t4 t11 = c2 * t3 + c3 * t4 bjg = b(jg+j) bjj = b(jj+j) u1 = bjg + bjj bjh = b(jh+j) bji = b(ji+j) u2 = bjh + bji u3 = bjg - bjj u4 = bjh - bji bjl = b(jl+j) bjh = bjl u5 = u1 + u2 u6 = c1 * ( u1 - u2 ) u7 = bjb - 0.25_dp * u5 b(jb+j) = bjb + u5 u8 = u7 + u6 u9 = u7 - u6 bjq = b(jq+j) bji = bjq u10 = c3 * u3 - c2 * u4 u11 = c2 * u3 + c3 * u4 a(jg+j) = t8 - u11 b(jg+j) = u8 + t11 ajj = t8 + u11 bjj = u8 - t11 a(jl+j) = t9 - u10 b(jl+j) = u9 + t10 a(jq+j) = t9 + u10 b(jq+j) = u9 - t10 !---------------------- ajo = a(jo+j) t1 = ajh + ajo ajm = a(jm+j) ajn = a(jn+j) t2 = ajm + ajn t3 = ajh - ajo t4 = ajm - ajn ajr = a(jr+j) ajn = ajr t5 = t1 + t2 t6 = c1 * ( t1 - t2 ) t7 = ajc - 0.25_dp * t5 a(jc+j) = ajc + t5 t8 = t7 + t6 t9 = t7 - t6 ajw = a(jw+j) ajo = ajw t10 = c3 * t3 - c2 * t4 t11 = c2 * t3 + c3 * t4 bjo = b(jo+j) u1 = bjh + bjo bjm = b(jm+j) bjn = b(jn+j) u2 = bjm + bjn u3 = bjh - bjo u4 = bjm - bjn bjr = b(jr+j) bjn = bjr u5 = u1 + u2 u6 = c1 * ( u1 - u2 ) u7 = bjc - 0.25_dp * u5 b(jc+j) = bjc + u5 u8 = u7 + u6 u9 = u7 - u6 bjw = b(jw+j) bjo = bjw 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 !---------------------- ajt = a(jt+j) t1 = aji + ajt ajs = a(js+j) t2 = ajn + ajs t3 = aji - ajt t4 = ajn - ajs ajx = a(jx+j) ajt = ajx t5 = t1 + t2 t6 = c1 * ( t1 - t2 ) ajp = a(jp+j) t7 = ajp - 0.25_dp * t5 ax = ajp + t5 t8 = t7 + t6 t9 = t7 - t6 a(jp+j) = ajd t10 = c3 * t3 - c2 * t4 t11 = c2 * t3 + c3 * t4 a(jd+j) = ax bjt = b(jt+j) u1 = bji + bjt bjs = b(js+j) u2 = bjn + bjs u3 = bji - bjt u4 = bjn - bjs bjx = b(jx+j) bjt = bjx u5 = u1 + u2 u6 = c1 * ( u1 - u2 ) bjp = b(jp+j) u7 = bjp - 0.25_dp * u5 bx = bjp + u5 u8 = u7 + u6 u9 = u7 - u6 b(jp+j) = bjd 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 !---------------------- ajv = a(jv+j) ajy = a(jy+j) t1 = ajv + ajy t2 = ajo + ajt t3 = ajv - ajy t4 = ajo - ajt a(jv+j) = ajj t5 = t1 + t2 t6 = c1 * ( t1 - t2 ) aju = a(ju+j) t7 = aju - 0.25_dp * t5 ax = aju + t5 t8 = t7 + t6 t9 = t7 - t6 a(ju+j) = aje t10 = c3 * t3 - c2 * t4 t11 = c2 * t3 + c3 * t4 a(je+j) = ax bjv = b(jv+j) bjy = b(jy+j) u1 = bjv + bjy u2 = bjo + bjt u3 = bjv - bjy u4 = bjo - bjt b(jv+j) = bjj u5 = u1 + u2 u6 = c1 * ( u1 - u2 ) bju = b(ju+j) u7 = bju - 0.25_dp * u5 bx = bju + u5 u8 = u7 + u6 u9 = u7 - u6 b(ju+j) = bje 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 END DO ELSE DO l = 1, nvex ajb = a(jb+j) aje = a(je+j) t1 = ajb + aje ajc = a(jc+j) ajd = a(jd+j) t2 = ajc + ajd t3 = ajb - aje t4 = ajc - ajd ajf = a(jf+j) ajb = ajf t5 = t1 + t2 t6 = c1 * ( t1 - t2 ) aja = a(ja+j) t7 = aja - 0.25_dp * t5 a(ja+j) = aja + t5 t8 = t7 + t6 t9 = t7 - t6 ajk = a(jk+j) ajc = ajk t10 = c3 * t3 - c2 * t4 t11 = c2 * t3 + c3 * t4 bjb = b(jb+j) bje = b(je+j) u1 = bjb + bje bjc = b(jc+j) bjd = b(jd+j) u2 = bjc + bjd u3 = bjb - bje u4 = bjc - bjd bjf = b(jf+j) bjb = bjf u5 = u1 + u2 u6 = c1 * ( u1 - u2 ) bja = b(ja+j) u7 = bja - 0.25_dp * u5 b(ja+j) = bja + u5 u8 = u7 + u6 u9 = u7 - u6 bjk = b(jk+j) bjc = bjk 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) aje = co4*(t8+u11) - si4*(u8-t11) bje = si4*(t8+u11) + co4*(u8-t11) a(jk+j) = co2*(t9-u10) - si2*(u9+t10) b(jk+j) = si2*(t9-u10) + co2*(u9+t10) ajd = co3*(t9+u10) - si3*(u9-t10) bjd = si3*(t9+u10) + co3*(u9-t10) !---------------------- ajg = a(jg+j) ajj = a(jj+j) t1 = ajg + ajj ajh = a(jh+j) aji = a(ji+j) t2 = ajh + aji t3 = ajg - ajj t4 = ajh - aji ajl = a(jl+j) ajh = ajl t5 = t1 + t2 t6 = c1 * ( t1 - t2 ) t7 = ajb - 0.25_dp * t5 a(jb+j) = ajb + t5 t8 = t7 + t6 t9 = t7 - t6 ajq = a(jq+j) aji = ajq t10 = c3 * t3 - c2 * t4 t11 = c2 * t3 + c3 * t4 bjg = b(jg+j) bjj = b(jj+j) u1 = bjg + bjj bjh = b(jh+j) bji = b(ji+j) u2 = bjh + bji u3 = bjg - bjj u4 = bjh - bji bjl = b(jl+j) bjh = bjl u5 = u1 + u2 u6 = c1 * ( u1 - u2 ) u7 = bjb - 0.25_dp * u5 b(jb+j) = bjb + u5 u8 = u7 + u6 u9 = u7 - u6 bjq = b(jq+j) bji = bjq 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) ajj = co4*(t8+u11) - si4*(u8-t11) bjj = 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) !---------------------- ajo = a(jo+j) t1 = ajh + ajo ajm = a(jm+j) ajn = a(jn+j) t2 = ajm + ajn t3 = ajh - ajo t4 = ajm - ajn ajr = a(jr+j) ajn = ajr t5 = t1 + t2 t6 = c1 * ( t1 - t2 ) t7 = ajc - 0.25_dp * t5 a(jc+j) = ajc + t5 t8 = t7 + t6 t9 = t7 - t6 ajw = a(jw+j) ajo = ajw t10 = c3 * t3 - c2 * t4 t11 = c2 * t3 + c3 * t4 bjo = b(jo+j) u1 = bjh + bjo bjm = b(jm+j) bjn = b(jn+j) u2 = bjm + bjn u3 = bjh - bjo u4 = bjm - bjn bjr = b(jr+j) bjn = bjr u5 = u1 + u2 u6 = c1 * ( u1 - u2 ) u7 = bjc - 0.25_dp * u5 b(jc+j) = bjc + u5 u8 = u7 + u6 u9 = u7 - u6 bjw = b(jw+j) bjo = bjw 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) !---------------------- ajt = a(jt+j) t1 = aji + ajt ajs = a(js+j) t2 = ajn + ajs t3 = aji - ajt t4 = ajn - ajs ajx = a(jx+j) ajt = ajx t5 = t1 + t2 t6 = c1 * ( t1 - t2 ) ajp = a(jp+j) t7 = ajp - 0.25_dp * t5 ax = ajp + t5 t8 = t7 + t6 t9 = t7 - t6 a(jp+j) = ajd t10 = c3 * t3 - c2 * t4 t11 = c2 * t3 + c3 * t4 a(jd+j) = ax bjt = b(jt+j) u1 = bji + bjt bjs = b(js+j) u2 = bjn + bjs u3 = bji - bjt u4 = bjn - bjs bjx = b(jx+j) bjt = bjx u5 = u1 + u2 u6 = c1 * ( u1 - u2 ) bjp = b(jp+j) u7 = bjp - 0.25_dp * u5 bx = bjp + u5 u8 = u7 + u6 u9 = u7 - u6 b(jp+j) = bjd 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) !---------------------- ajv = a(jv+j) ajy = a(jy+j) t1 = ajv + ajy t2 = ajo + ajt t3 = ajv - ajy t4 = ajo - ajt a(jv+j) = ajj t5 = t1 + t2 t6 = c1 * ( t1 - t2 ) aju = a(ju+j) t7 = aju - 0.25_dp * t5 ax = aju + t5 t8 = t7 + t6 t9 = t7 - t6 a(ju+j) = aje t10 = c3 * t3 - c2 * t4 t11 = c2 * t3 + c3 * t4 a(je+j) = ax bjv = b(jv+j) bjy = b(jy+j) u1 = bjv + bjy u2 = bjo + bjt u3 = bjv - bjy u4 = bjo - bjt b(jv+j) = bjj u5 = u1 + u2 u6 = c1 * ( u1 - u2 ) bju = b(ju+j) u7 = bju - 0.25_dp * u5 bx = bju + u5 u8 = u7 + u6 u9 = u7 - u6 b(ju+j) = bje 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 END DO END IF !-----(end of loop across transforms) ja = ja + jstepx IF (ja < istart) ja = ja + ninc END DO END DO END DO !-----( end of double loop for this k ) kk = kk + 2*la END DO !-----( end of loop over values of k ) la = 5*la END DO !-----( end of loop on type II radix-5 passes ) !-----( nvex transforms completed) 490 istart = istart + nvex * jump END DO !-----( end of loop on blocks of transforms ) RETURN END SUBROUTINE gpfa5f END MODULE FFT_235