gpfa2f.c

来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 1,107 行 · 第 1/3 页

                        if (jj < istart) {
                            jj += ninc;
                        }
                        jk = jj + jstepl;
                        if (jk < istart) {
                            jk += ninc;
                        }
                        jl = jk + jstepl;
                        if (jl < istart) {
                            jl += ninc;
                        }
                        jm = ji + laincl;
                        if (jm < istart) {
                            jm += ninc;
                        }
                        jn = jm + jstepl;
                        if (jn < istart) {
                            jn += ninc;
                        }
                        jo = jn + jstepl;
                        if (jo < istart) {
                            jo += ninc;
                        }
                        jp = jo + jstepl;
                        if (jp < istart) {
                            jp += ninc;
                        }
                        j = 0;

/*  loop across transforms */
                        for (l = 0; l < nvex; ++l) {
                            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 += *jump;
                        }
                        ja += jstepx;
                        if (ja < istart) {
                            ja += ninc;
                        }
                    }
                }
            }

/*  finished if last pass */
            if (ipass == m-1) {
                goto L490;
            }

            kk = la << 1;

/*     loop on nonzero k */
            for (k = ink; ink < 0 ? k >= jstep - ink : k <= jstep - ink; k += ink) {
                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];

/*  double loop along first transform in block */
                i__1 = jstep << 2;
                for (ll = k; i__1 < 0 ? ll >= (la - 1) * ink : ll <= (la - 1) * ink; ll += i__1) {

                    i__1 = (la << 2) * ink;
                    for (jjj = ll; i__1 < 0 ? jjj >= (*n - 1) * *inc : jjj <= (*n - 1) * *inc; jjj += i__1) {
                        ja = istart + jjj;

/*     "transverse" loop */
                        for (nu = 0; 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 = ja + laincl;
                            if (je < istart) {
                                je += ninc;
                            }
                            jf = je + jstepl;
                            if (jf < istart) {
                                jf += ninc;
                            }
                            jg = jf + jstepl;
                            if (jg < istart) {
                                jg += ninc;
                            }
                            jh = jg + jstepl;
                            if (jh < istart) {
                                jh += ninc;
                            }
                            ji = je + laincl;
                            if (ji < istart) {
                                ji += ninc;
                            }
                            jj = ji + jstepl;
                            if (jj < istart) {
                                jj += ninc;
                            }
                            jk = jj + jstepl;
                            if (jk < istart) {
                                jk += ninc;
                            }
                            jl = jk + jstepl;
                            if (jl < istart) {
                                jl += ninc;
                            }
                            jm = ji + laincl;
                            if (jm < istart) {
                                jm += ninc;
                            }
                            jn = jm + jstepl;
                            if (jn < istart) {
                                jn += ninc;
                            }
                            jo = jn + jstepl;
                            if (jo < istart) {
                                jo += ninc;
                            }
                            jp = jo + jstepl;
                            if (jp < istart) {
                                jp += ninc;
                            }
                            j = 0;

/*  loop across transforms */
/* dir$ ivdep, shortloop */
                            for (l = 0; l < nvex; ++l) {
                                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 += *jump;
                            }
                            ja += jstepx;
                            if (ja < istart) {
                                ja += ninc;
                            }
                        }
                    }
                }
                kk += la << 1;
            }
            la <<= 2;
        }
L490:
        istart += nvex * *jump;
    }
} /* gpfa2f_ */