📄 gpfa2f.c
字号:
ajd = a[jd + j];
/*< t1 = ajb + ajd >*/
t1 = ajb + ajd;
/*< t3 = ss * ( ajb - ajd ) >*/
t3 = ss * (ajb - ajd);
/*< bja = b(ja+j) >*/
bja = b[ja + j];
/*< bjc = b(jc+j) >*/
bjc = b[jc + j];
/*< u0 = bja + bjc >*/
u0 = bja + bjc;
/*< u2 = bja - bjc >*/
u2 = bja - bjc;
/*< bjb = b(jb+j) >*/
bjb = b[jb + j];
/*< bjd = b(jd+j) >*/
bjd = b[jd + j];
/*< u1 = bjb + bjd >*/
u1 = bjb + bjd;
/*< u3 = ss * ( bjb - bjd ) >*/
u3 = ss * (bjb - bjd);
/*< a(ja+j) = t0 + t1 >*/
a[ja + j] = t0 + t1;
/*< b(ja+j) = u0 + u1 >*/
b[ja + j] = u0 + u1;
/*< a(jb+j) = co1*(t2-u3) - si1*(u2+t3) >*/
a[jb + j] = co1 * (t2 - u3) - si1 * (u2 + t3);
/*< b(jb+j) = si1*(t2-u3) + co1*(u2+t3) >*/
b[jb + j] = si1 * (t2 - u3) + co1 * (u2 + t3);
/*< a(jc+j) = co2*(t0-t1) - si2*(u0-u1) >*/
a[jc + j] = co2 * (t0 - t1) - si2 * (u0 - u1);
/*< b(jc+j) = si2*(t0-t1) + co2*(u0-u1) >*/
b[jc + j] = si2 * (t0 - t1) + co2 * (u0 - u1);
/*< a(jd+j) = co3*(t2+u3) - si3*(u2-t3) >*/
a[jd + j] = co3 * (t2 + u3) - si3 * (u2 - t3);
/*< b(jd+j) = si3*(t2+u3) + co3*(u2-t3) >*/
b[jd + j] = si3 * (t2 + u3) + co3 * (u2 - t3);
/*< j = j + jump >*/
j += *jump;
/*< 130 continue >*/
/* L130: */
}
/* -----( end of loop across transforms ) */
/*< ja = ja + jstepx >*/
ja += jstepx;
/*< if (ja.lt.istart) ja = ja + ninc >*/
if (ja < istart) {
ja += ninc;
}
/*< 135 continue >*/
/* L135: */
}
/*< 140 continue >*/
/* L140: */
}
/* -----( end of loop along transforms ) */
/*< kk = kk + 2*la >*/
kk += la << 1;
/*< 150 continue >*/
/* L150: */
}
/* -----( end of loop on nonzero k ) */
/*< la = 4*la >*/
la <<= 2;
/*< 160 continue >*/
/* L160: */
}
/* -----( end of loop on type I radix-4 passes) */
/* central radix-2 pass */
/* -------------------- */
/*< 200 continue >*/
L200:
/*< if (m2.eq.0) go to 300 >*/
if (m2 == 0) {
goto L300;
}
/*< jstep = (n*inc) / (2*la) >*/
jstep = *n * *inc / (la << 1);
/*< jstepl = jstep - ninc >*/
jstepl = jstep - ninc;
/* k=0 loop (no twiddle factors) */
/* ----------------------------- */
/*< do 220 jjj = 0 , (n-1)*inc , 2*jstep >*/
i__2 = (*n - 1) * *inc;
i__3 = jstep << 1;
for (jjj = 0; i__3 < 0 ? jjj >= i__2 : jjj <= i__2; jjj += i__3) {
/*< ja = istart + jjj >*/
ja = istart + jjj;
/* "transverse" loop */
/* ----------------- */
/*< do 215 nu = 1 , inq >*/
i__4 = inq;
for (nu = 1; nu <= i__4; ++nu) {
/*< jb = ja + jstepl >*/
jb = ja + jstepl;
/*< if (jb.lt.istart) jb = jb + ninc >*/
if (jb < istart) {
jb += ninc;
}
/*< j = 0 >*/
j = 0;
/* loop across transforms */
/* ---------------------- */
/* dir$ ivdep, shortloop */
/*< do 210 l = 1 , nvex >*/
i__6 = nvex;
for (l = 1; l <= i__6; ++l) {
/*< aja = a(ja+j) >*/
aja = a[ja + j];
/*< ajb = a(jb+j) >*/
ajb = a[jb + j];
/*< t0 = aja - ajb >*/
t0 = aja - ajb;
/*< a(ja+j) = aja + ajb >*/
a[ja + j] = aja + ajb;
/*< a(jb+j) = t0 >*/
a[jb + j] = t0;
/*< bja = b(ja+j) >*/
bja = b[ja + j];
/*< bjb = b(jb+j) >*/
bjb = b[jb + j];
/*< u0 = bja - bjb >*/
u0 = bja - bjb;
/*< b(ja+j) = bja + bjb >*/
b[ja + j] = bja + bjb;
/*< b(jb+j) = u0 >*/
b[jb + j] = u0;
/*< j = j + jump >*/
j += *jump;
/*< 210 continue >*/
/* L210: */
}
/* -----(end of loop across transforms) */
/*< ja = ja + jstepx >*/
ja += jstepx;
/*< if (ja.lt.istart) ja = ja + ninc >*/
if (ja < istart) {
ja += ninc;
}
/*< 215 continue >*/
/* L215: */
}
/*< 220 continue >*/
/* L220: */
}
/* finished if n2=2 */
/* ---------------- */
/*< if (n2.eq.2) go to 490 >*/
if (n2 == 2) {
goto L490;
}
/*< kk = 2 * la >*/
kk = la << 1;
/* loop on nonzero k */
/* ----------------- */
/*< do 260 k = ink , jstep - ink , ink >*/
i__3 = jstep - ink;
i__2 = ink;
for (k = ink; i__2 < 0 ? k >= i__3 : k <= i__3; k += i__2) {
/*< co1 = trigs(kk+1) >*/
co1 = trigs[kk + 1];
/*< si1 = s*trigs(kk+2) >*/
si1 = s * trigs[kk + 2];
/* loop along transforms */
/* --------------------- */
/*< do 250 jjj = k , (n-1)*inc , 2*jstep >*/
i__4 = (*n - 1) * *inc;
i__6 = jstep << 1;
for (jjj = k; i__6 < 0 ? jjj >= i__4 : jjj <= i__4; jjj += i__6) {
/*< ja = istart + jjj >*/
ja = istart + jjj;
/* "transverse" loop */
/* ----------------- */
/*< do 245 nu = 1 , inq >*/
i__5 = inq;
for (nu = 1; nu <= i__5; ++nu) {
/*< jb = ja + jstepl >*/
jb = ja + jstepl;
/*< if (jb.lt.istart) jb = jb + ninc >*/
if (jb < istart) {
jb += ninc;
}
/*< j = 0 >*/
j = 0;
/* loop across transforms */
/* ---------------------- */
/*< if (kk.eq.n2/2) then >*/
if (kk == n2 / 2) {
/* dir$ ivdep, shortloop */
/*< do 230 l = 1 , nvex >*/
i__7 = nvex;
for (l = 1; l <= i__7; ++l) {
/*< aja = a(ja+j) >*/
aja = a[ja + j];
/*< ajb = a(jb+j) >*/
ajb = a[jb + j];
/*< t0 = ss * ( aja - ajb ) >*/
t0 = ss * (aja - ajb);
/*< a(ja+j) = aja + ajb >*/
a[ja + j] = aja + ajb;
/*< bjb = b(jb+j) >*/
bjb = b[jb + j];
/*< bja = b(ja+j) >*/
bja = b[ja + j];
/*< a(jb+j) = ss * ( bjb - bja ) >*/
a[jb + j] = ss * (bjb - bja);
/*< b(ja+j) = bja + bjb >*/
b[ja + j] = bja + bjb;
/*< b(jb+j) = t0 >*/
b[jb + j] = t0;
/*< j = j + jump >*/
j += *jump;
/*< 230 continue >*/
/* L230: */
}
/*< else >*/
} else {
/* dir$ ivdep, shortloop */
/*< do 240 l = 1 , nvex >*/
i__7 = nvex;
for (l = 1; l <= i__7; ++l) {
/*< aja = a(ja+j) >*/
aja = a[ja + j];
/*< ajb = a(jb+j) >*/
ajb = a[jb + j];
/*< t0 = aja - ajb >*/
t0 = aja - ajb;
/*< a(ja+j) = aja + ajb >*/
a[ja + j] = aja + ajb;
/*< bja = b(ja+j) >*/
bja = b[ja + j];
/*< bjb = b(jb+j) >*/
bjb = b[jb + j];
/*< u0 = bja - bjb >*/
u0 = bja - bjb;
/*< b(ja+j) = bja + bjb >*/
b[ja + j] = bja + bjb;
/*< a(jb+j) = co1*t0 - si1*u0 >*/
a[jb + j] = co1 * t0 - si1 * u0;
/*< b(jb+j) = si1*t0 + co1*u0 >*/
b[jb + j] = si1 * t0 + co1 * u0;
/*< j = j + jump >*/
j += *jump;
/*< 240 continue >*/
/* L240: */
}
/*< endif >*/
}
/* -----(end of loop across transforms) */
/*< ja = ja + jstepx >*/
ja += jstepx;
/*< if (ja.lt.istart) ja = ja + ninc >*/
if (ja < istart) {
ja += ninc;
}
/*< 245 continue >*/
/* L245: */
}
/*< 250 continue >*/
/* L250: */
}
/* -----(end of loop along transforms) */
/*< kk = kk + 2 * la >*/
kk += la << 1;
/*< 260 continue >*/
/* L260: */
}
/* -----(end of loop on nonzero k) */
/* -----(end of radix-2 pass) */
/*< la = 2 * la >*/
la <<= 1;
/*< go to 400 >*/
goto L400;
/* central radix-8 pass */
/* -------------------- */
/*< 300 continue >*/
L300:
/*< if (m8.eq.0) go to 400 >*/
if (m8 == 0) {
goto L400;
}
/*< jstep = (n*inc) / (8*la) >*/
jstep = *n * *inc / (la << 3);
/*< jstepl = jstep - ninc >*/
jstepl = jstep - ninc;
/*< mu = mod(inq,8) >*/
mu = inq % 8;
/*< if (isign.eq.-1) mu = 8 - mu >*/
if (*isign == -1) {
mu = 8 - mu;
}
/*< c1 = 1.0 >*/
c1 = (float)1.;
/*< if (mu.eq.3.or.mu.eq.7) c1 = -1.0 >*/
if (mu == 3 || mu == 7) {
c1 = (float)-1.;
}
/*< c2 = sqrt(0.5) >*/
c2 = sqrt((float).5);
/*< if (mu.eq.3.or.mu.eq.5) c2 = -c2 >*/
if (mu == 3 || mu == 5) {
c2 = -c2;
}
/*< c3 = c1 * c2 >*/
c3 = c1 * c2;
/* stage 1 */
/* ------- */
/*< do 320 k = 0 , jstep - ink , ink >*/
i__2 = jstep - ink;
i__3 = ink;
for (k = 0; i__3 < 0 ? k >= i__2 : k <= i__2; k += i__3) {
/*< do 315 jjj = k , (n-1)*inc , 8*jstep >*/
i__6 = (*n - 1) * *inc;
i__4 = jstep << 3;
for (jjj = k; i__4 < 0 ? jjj >= i__6 : jjj <= i__6; jjj += i__4) {
/*< ja = istart + jjj >*/
ja = istart + jjj;
/* "transverse" loop */
/* ----------------- */
/*< do 312 nu = 1 , inq >*/
i__5 = inq;
for (nu = 1; nu <= i__5; ++nu) {
/*< jb = ja + jstepl >*/
jb = ja + jstepl;
/*< if (jb.lt.istart) jb = jb + ninc >*/
if (jb < istart) {
jb += ninc;
}
/*< jc = jb + jstepl >*/
jc = jb + jstepl;
/*< if (jc.lt.istart) jc = jc + ninc >*/
if (jc < istart) {
jc += ninc;
}
/*< jd = jc + jstepl >*/
jd = jc + jstepl;
/*< if (jd.lt.istart) jd = jd + ninc >*/
if (jd < istart) {
jd += ninc;
}
/*< je = jd + jstepl >*/
je = jd + jstepl;
/*< if (je.lt.istart) je = je + ninc >*/
if (je < istart) {
je += ninc;
}
/*< jf = je + jstepl >*/
jf = je + jstepl;
/*< if (jf.lt.istart) jf = jf + ninc >*/
if (jf < istart) {
jf += ninc;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -