gpfa5f.c
来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 911 行 · 第 1/3 页
C
911 行
#include "f2c.h"
#include "netlib.h"
/* Table of constant values */
static integer c__5 = 5;
/* fortran version of *gpfa5* - */
/* radix-5 section of self-sorting, in-place, */
/* generalized pfa */
/* ------------------------------------------------------------------- */
/* Subroutine */ void gpfa5f_(real *a, real *b, const real *trigs, const integer *inc,
const integer *jump, const integer *n, const integer *mm, const integer *lot, const integer *isign)
{
/* Initialized data */
static real sin36 = .587785252292473f;
static real sin72 = .951056516295154f;
static real qrt5 = .559016994374947f;
static integer lvr = 128;
/* System generated locals */
integer i__3, i__4, i__5, i__6, i__7, i__8;
/* Local variables */
static integer ninc, left, nvex, j, k, l, m;
static real s;
static integer ipass, nblox;
static real c1, c2, c3;
static integer jstep, n5;
static real t1, t2, t3, t4, t5, t6, t7, t8, t9, u1, u2, u3, u4, u5, u6, u7, u8, u9;
static integer ja, jb, jc, jd, je, jf, jg, jh, ji, jj, jk, jl, jm, jn, jo, jp, jq, jr, js, jt, ju, jv, jw, jx, jy;
static real t10, t11, u10, u11;
static integer ll, mu, nu, laincl, la, nb, mh, kk;
static real ax, bx;
static integer jstepl;
static real co1, co2, co3;
static integer istart;
static real co4;
static integer jstepx;
static real si1, si2, si3, si4, aja, ajb, ajc, ajd, aje, bjb, bje, bjc, bjd, bja, ajf, ajk, bjf, bjk, ajg, ajj, ajh;
static integer jjj;
static real aji, ajl;
static integer ink;
static real ajq, bjg, bjj, bjh, bji;
static integer inq;
static real bjl, bjq, ajo, ajm, ajn, ajr, ajw, bjo, bjm, bjn, bjr, bjw,
ajt, ajs, ajx, ajp, bjt, bjs, bjx, bjp, ajv, ajy, aju, bjv, bjy, bju;
/* *************************************************************** */
/* * * */
/* * 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. * */
/* * * */
/* *************************************************************** */
n5 = pow_ii(&c__5, mm);
inq = *n / n5;
jstepx = (n5 - *n) * *inc;
ninc = *n * *inc;
ink = *inc * inq;
mu = inq % 5;
if (*isign == -1) {
mu = 5 - mu;
}
m = *mm;
mh = (m + 1) / 2;
s = (real) (*isign);
c1 = qrt5;
c2 = sin72;
c3 = sin36;
if (mu == 2 || mu == 3) {
c1 = -c1;
c2 = sin36;
c3 = sin72;
}
if (mu == 3 || mu == 4) {
c2 = -c2;
}
if (mu == 2 || mu == 4) {
c3 = -c3;
}
nblox = (*lot - 1) / lvr + 1;
left = *lot;
s = (real) (*isign);
istart = 0;
/* loop on blocks of lvr transforms */
/* -------------------------------- */
for (nb = 1; nb <= nblox; ++nb) {
if (left <= lvr) {
nvex = left;
} else if (left < lvr << 1) {
nvex = left / 2;
nvex += nvex % 2;
} else {
nvex = lvr;
}
left -= nvex;
la = 1;
/* loop on type I radix-5 passes */
/* ----------------------------- */
for (ipass = 0; ipass < mh; ++ipass) {
jstep = *n * *inc / (la * 5);
jstepl = jstep - ninc;
kk = 0;
/* loop on k */
/* --------- */
i__3 = jstep - ink;
for (k = 0; ink < 0 ? k >= i__3 : k <= i__3; k += ink) {
if (k > 0) {
co1 = trigs[kk]; si1 = s * trigs[kk + 1];
co2 = trigs[kk * 2]; si2 = s * trigs[kk * 2 + 1];
co3 = trigs[kk * 3]; si3 = s * trigs[kk * 3 + 1];
co4 = trigs[kk * 4]; si4 = s * trigs[kk * 4 + 1];
}
/* loop along transform */
/* -------------------- */
i__5 = (*n - 1) * *inc;
i__6 = jstep * 5;
for (jjj = k; i__6 < 0 ? jjj >= i__5 : jjj <= i__5; jjj += i__6) {
ja = istart + jjj;
/* "transverse" loop */
/* ----------------- */
for (nu = 1; nu <= inq; ++nu) {
jb = ja + jstepl;
if (jb < istart) {
jb += ninc;
}
jc = jb + jstepl;
if (jc < istart) {
jc += ninc;
}
jd = jc + jstepl;
if (jd < istart) {
jd += ninc;
}
je = jd + jstepl;
if (je < istart) {
je += ninc;
}
j = 0;
/* loop across transforms */
/* ---------------------- */
if (k == 0) {
/* dir$ ivdep, shortloop */
for (l = 1; l <= nvex; ++l) {
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 - t5 * .25f;
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 - u5 * .25f;
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 += *jump;
}
} else {
/* dir$ ivdep,shortloop */
for (l = 1; l <= nvex; ++l) {
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 - t5 * .25f;
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 - u5 * .25f;
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 += *jump;
}
}
ja += jstepx;
if (ja < istart) {
ja += ninc;
}
}
}
kk += la << 1;
}
la *= 5;
}
if (*n == 5) {
goto L490;
}
/* loop on type II radix-5 passes */
/* ------------------------------ */
for (ipass = mh; ipass < m; ++ipass) {
jstep = *n * *inc / (la * 5);
jstepl = jstep - ninc;
laincl = la * ink - ninc;
kk = 0;
/* loop on k */
/* --------- */
i__4 = jstep - ink;
for (k = 0; ink < 0 ? k >= i__4 : k <= i__4; k += ink) {
if (k > 0) {
co1 = trigs[kk]; si1 = s * trigs[kk + 1];
co2 = trigs[kk * 2]; si2 = s * trigs[kk * 2 + 1];
co3 = trigs[kk * 3]; si3 = s * trigs[kk * 3 + 1];
co4 = trigs[kk * 4]; si4 = s * trigs[kk * 4 + 1];
}
/* double loop along first transform in block */
/* ------------------------------------------ */
i__6 = (la - 1) * ink;
i__5 = jstep * 5;
for (ll = k; i__5 < 0 ? ll >= i__6 : ll <= i__6; ll += i__5) {
i__7 = (*n - 1) * *inc;
i__8 = la * 5 * ink;
for (jjj = ll; i__8 < 0 ? jjj >= i__7 : jjj <= i__7; jjj += i__8) {
ja = istart + jjj;
/* "transverse" loop */
/* ----------------- */
for (nu = 1; nu <= inq; ++nu) {
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?