gpfa2f.c

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

        c1 = 1.f;
        if (mu == 3 || mu == 7) {
            c1 = -1.f;
        }
        c2 = sqrtf(.5f);
        if (mu == 3 || mu == 5) {
            c2 = -c2;
        }
        c3 = c1 * c2;

/*  stage 1 */
        for (k = 0; ink < 0 ? k >= jstep - ink : k <= jstep - ink; k += ink) {
            i__1 = jstep << 3;
            for (jjj = k; 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 = jd + jstepl;
                    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;
                    }
                    j = 0;
/* dir$ ivdep, shortloop */
                    for (l = 0; l < nvex; ++l) {
                        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 += *jump;
                    }
                    ja += jstepx;
                    if (ja < istart) {
                        ja += ninc;
                    }
                }
            }
        }

/*  stage 2 */

/*  k=0 (no twiddle factors) */
        i__1 = jstep << 3;
        for (jjj = 0; 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 = jd + jstepl;
                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;
                }
                j = 0;
/* dir$ ivdep, shortloop */
                for (l = 0; l < nvex; ++l) {
                    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 += *jump;
                }
                ja += jstepx;
                if (ja < istart) {
                    ja += ninc;
                }
            }
        }

        if (n2 == 8) {
            goto L490;
        }

/*  loop on nonzero k */
        kk = la << 1;

        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];
            co4 = trigs[kk * 4]; si4 = s * trigs[kk * 4 + 1];
            co5 = trigs[kk * 5]; si5 = s * trigs[kk * 5 + 1];
            co6 = trigs[kk * 6]; si6 = s * trigs[kk * 6 + 1];
            co7 = trigs[kk * 7]; si7 = s * trigs[kk * 7 + 1];

            i__1 = jstep << 3;
            for (jjj = k; 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 = jd + jstepl;
                    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;
                    }
                    j = 0;
/* dir$ ivdep, shortloop */
                    for (l = 0; l < nvex; ++l) {
                        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 += *jump;
                    }
                    ja += jstepx;
                    if (ja < istart) {
                        ja += ninc;
                    }
                }
            }
            kk += la << 1;
        }

        la <<= 3;

/*  loop on type II radix-4 passes */
L400:
        mu = inq % 4;
        if (*isign == -1) {
            mu = 4 - mu;
        }
        ss = 1.f;
        if (mu == 3) {
            ss = -1.f;
        }

        for (ipass = mh; ipass < m; ++ipass) {
            jstep = *n * *inc / (la << 2);
            jstepl = jstep - ninc;
            laincl = la * ink - ninc;

/*  k=0 loop (no twiddle factors) */
            i__1 = jstep << 2;
            for (ll = 0; 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;