dgpfa3f.c

InsightToolkit-1.4.0(有大量的优化算法程序)

                        if (jf < istart) {
                            jf += ninc;
                        }
                        jg = jd + laincl;
                        if (jg < istart) {
                            jg += ninc;
                        }
                        jh = jg + jstepl;
                        if (jh < istart) {
                            jh += ninc;
                        }
                        ji = jh + jstepl;
                        if (ji < istart) {
                            ji += ninc;
                        }
                        j = 0;

/*  loop across transforms */
/*  ---------------------- */
/* dir$ ivdep, shortloop */
                        for (l = 1; l <= nvex; ++l) {
                            ajb = a[jb + j];
                            ajc = a[jc + j];
                            t1 = ajb + ajc;
                            aja = a[ja + j];
                            t2 = aja - t1 * .5;
                            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 - u1 * .5;
                            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 - t1 * .5;
                            t3 = c1 * (aje - ajf);
                            ajh = a[jh + j];
                            ajf = ajh;
                            bje = b[je + j];
                            bjf = b[jf + j];
                            u1 = bje + bjf;
                            u2 = bjb - u1 * .5;
                            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 - t1 * .5;
                            t3 = c1 * (ajf - aji);
                            t1 += ajg;
                            a[jg + j] = ajc;
                            bji = b[ji + j];
                            u1 = bjf + bji;
                            bjg = b[jg + j];
                            u2 = bjg - u1 * .5;
                            u3 = c1 * (bjf - bji);
                            u1 += bjg;
                            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 += *jump;
                        }
/* -----( end of loop across transforms ) */
                        ja += jstepx;
                        if (ja < istart) {
                            ja += ninc;
                        }
                    }
                }
            }
/* -----( end of double loop for k=0 ) */

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

            kk = la << 1;

/*     loop on nonzero k */
/*     ----------------- */
            i__4 = jstep - ink;
            for (k = ink; ink < 0 ? k >= i__4 : k <= i__4; k += ink) {
                co1 = trigs[kk];     si1 = s * trigs[kk + 1];
                co2 = trigs[kk * 2]; si2 = s * trigs[kk * 2 + 1];

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

                    i__7 = (*n - 1) * *inc;
                    i__8 = la * 3 * 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) {
                            jb = ja + jstepl;
                            if (jb < istart) {
                                jb += ninc;
                            }
                            jc = jb + jstepl;
                            if (jc < istart) {
                                jc += ninc;
                            }
                            jd = ja + laincl;
                            if (jd < istart) {
                                jd += ninc;
                            }
                            je = jd + jstepl;
                            if (je < istart) {
                                je += ninc;
                            }
                            jf = je + jstepl;
                            if (jf < istart) {
                                jf += ninc;
                            }
                            jg = jd + laincl;
                            if (jg < istart) {
                                jg += ninc;
                            }
                            jh = jg + jstepl;
                            if (jh < istart) {
                                jh += ninc;
                            }
                            ji = jh + jstepl;
                            if (ji < istart) {
                                ji += ninc;
                            }
                            j = 0;

/*  loop across transforms */
/*  ---------------------- */
/* dir$ ivdep, shortloop */
                            for (l = 1; l <= nvex; ++l) {
                                ajb = a[jb + j];
                                ajc = a[jc + j];
                                t1 = ajb + ajc;
                                aja = a[ja + j];
                                t2 = aja - t1 * .5;
                                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 - u1 * .5;
                                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 - t1 * .5;
                                t3 = c1 * (aje - ajf);
                                ajh = a[jh + j];
                                ajf = ajh;
                                bje = b[je + j];
                                bjf = b[jf + j];
                                u1 = bje + bjf;
                                u2 = bjb - u1 * .5;
                                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 - t1 * .5;
                                t3 = c1 * (ajf - aji);
                                t1 += ajg;
                                a[jg + j] = ajc;
                                bji = b[ji + j];
                                u1 = bjf + bji;
                                bjg = b[jg + j];
                                u2 = bjg - u1 * .5;
                                u3 = c1 * (bjf - bji);
                                u1 += bjg;
                                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 += *jump;
                            }
/* -----(end of loop across transforms) */
                            ja += jstepx;
                            if (ja < istart) {
                                ja += ninc;
                            }
                        }
                    }
                }
/* -----( end of double loop for this k ) */
                kk += la << 1;
            }
/* -----( end of loop over values of k ) */
            la *= 3;
        }
/* -----( end of loop on type II radix-3 passes ) */
/* -----( nvex transforms completed) */
L490:
        istart += nvex * *jump;
    }
/* -----( end of loop on blocks of transforms ) */

} /* dgpfa3f_ */