gpfa3f.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/*<       do 130 l = 1 , nvex >*/
                        i__8 = nvex;
                        for (l = 1; l <= i__8; ++l) {
/*<       ajb = a(jb+j) >*/
                            ajb = a[jb + j];
/*<       ajc = a(jc+j) >*/
                            ajc = a[jc + j];
/*<       t1 = ajb + ajc >*/
                            t1 = ajb + ajc;
/*<       aja = a(ja+j) >*/
                            aja = a[ja + j];
/*<       t2 = aja - 0.5 * t1 >*/
                            t2 = aja - t1 * (float).5;
/*<       t3 = c1 * ( ajb - ajc ) >*/
                            t3 = c1 * (ajb - ajc);
/*<       bjb = b(jb+j) >*/
                            bjb = b[jb + j];
/*<       bjc = b(jc+j) >*/
                            bjc = b[jc + j];
/*<       u1 = bjb + bjc >*/
                            u1 = bjb + bjc;
/*<       bja = b(ja+j) >*/
                            bja = b[ja + j];
/*<       u2 = bja - 0.5 * u1 >*/
                            u2 = bja - u1 * (float).5;
/*<       u3 = c1 * ( bjb - bjc ) >*/
                            u3 = c1 * (bjb - bjc);
/*<       a(ja+j) = aja + t1 >*/
                            a[ja + j] = aja + t1;
/*<       b(ja+j) = bja + u1 >*/
                            b[ja + j] = bja + 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*(t2+u3) - si2*(u2-t3) >*/
                            a[jc + j] = co2 * (t2 + u3) - si2 * (u2 - t3);
/*<       b(jc+j) = si2*(t2+u3) + co2*(u2-t3) >*/
                            b[jc + j] = si2 * (t2 + u3) + co2 * (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 = 3*la >*/
            la *= 3;
/*<   160 continue >*/
/* L160: */
        }
/* -----( end of loop on type I radix-3 passes) */

/*  loop on type II radix-3 passes */
/*  ------------------------------ */
/*<   400 continue >*/
/* L400: */

/*<       do 480 ipass = mh+1 , m >*/
        i__2 = m;
        for (ipass = mh + 1; ipass <= i__2; ++ipass) {
/*<       jstep = (n*inc) / (3*la) >*/
            jstep = *n * *inc / (la * 3);
/*<       jstepl = jstep - ninc >*/
            jstepl = jstep - ninc;
/*<       laincl = la*ink - ninc >*/
            laincl = la * ink - ninc;

/*  k=0 loop (no twiddle factors) */
/*  ----------------------------- */
/*<       do 430 ll = 0 , (la-1)*ink , 3*jstep >*/
            i__3 = (la - 1) * ink;
            i__4 = jstep * 3;
            for (ll = 0; i__4 < 0 ? ll >= i__3 : ll <= i__3; ll += i__4) {

/*<       do 420 jjj = ll , (n-1)*inc , 3*la*ink >*/
                i__6 = (*n - 1) * *inc;
                i__5 = la * 3 * ink;
                for (jjj = ll; i__5 < 0 ? jjj >= i__6 : jjj <= i__6; jjj += 
                        i__5) {
/*<       ja = istart + jjj >*/
                    ja = istart + jjj;

/*  "transverse" loop */
/*  ----------------- */
/*<       do 415 nu = 1 , inq >*/
                    i__7 = inq;
                    for (nu = 1; nu <= i__7; ++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 = ja + laincl >*/
                        jd = ja + laincl;
/*<       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;
                        }
/*<       jg = jd + laincl >*/
                        jg = jd + laincl;
/*<       if (jg.lt.istart) jg = jg + ninc >*/
                        if (jg < istart) {
                            jg += ninc;
                        }
/*<       jh = jg + jstepl >*/
                        jh = jg + jstepl;
/*<       if (jh.lt.istart) jh = jh + ninc >*/
                        if (jh < istart) {
                            jh += ninc;
                        }
/*<       ji = jh + jstepl >*/
                        ji = jh + jstepl;
/*<       if (ji.lt.istart) ji = ji + ninc >*/
                        if (ji < istart) {
                            ji += ninc;
                        }
/*<       j = 0 >*/
                        j = 0;

/*  loop across transforms */
/*  ---------------------- */
/* dir$ ivdep, shortloop */
/*<       do 410 l = 1 , nvex >*/
                        i__8 = nvex;
                        for (l = 1; l <= i__8; ++l) {
/*<       ajb = a(jb+j) >*/
                            ajb = a[jb + j];
/*<       ajc = a(jc+j) >*/
                            ajc = a[jc + j];
/*<       t1 = ajb + ajc >*/
                            t1 = ajb + ajc;
/*<       aja = a(ja+j) >*/
                            aja = a[ja + j];
/*<       t2 = aja - 0.5 * t1 >*/
                            t2 = aja - t1 * (float).5;
/*<       t3 = c1 * ( ajb - ajc ) >*/
                            t3 = c1 * (ajb - ajc);
/*<       ajd = a(jd+j) >*/
                            ajd = a[jd + j];
/*<       ajb =  ajd >*/
                            ajb = ajd;
/*<       bjb = b(jb+j) >*/
                            bjb = b[jb + j];
/*<       bjc = b(jc+j) >*/
                            bjc = b[jc + j];
/*<       u1 = bjb + bjc >*/
                            u1 = bjb + bjc;
/*<       bja = b(ja+j) >*/
                            bja = b[ja + j];
/*<       u2 = bja - 0.5 * u1 >*/
                            u2 = bja - u1 * (float).5;
/*<       u3 = c1 * ( bjb - bjc ) >*/
                            u3 = c1 * (bjb - bjc);
/*<       bjd = b(jd+j) >*/
                            bjd = b[jd + j];
/*<       bjb =  bjd >*/
                            bjb = bjd;
/*<       a(ja+j) = aja + t1 >*/
                            a[ja + j] = aja + t1;
/*<       b(ja+j) = bja + u1 >*/
                            b[ja + j] = bja + u1;
/*<       a(jd+j) = t2 - u3 >*/
                            a[jd + j] = t2 - u3;
/*<       b(jd+j) = u2 + t3 >*/
                            b[jd + j] = u2 + t3;
/*<       ajc =  t2 + u3 >*/
                            ajc = t2 + u3;
/*<       bjc =  u2 - t3 >*/
                            bjc = u2 - t3;
/* ---------------------- */
/*<       aje = a(je+j) >*/
                            aje = a[je + j];
/*<       ajf = a(jf+j) >*/
                            ajf = a[jf + j];
/*<       t1 = aje + ajf >*/
                            t1 = aje + ajf;
/*<       t2 = ajb - 0.5 * t1 >*/
                            t2 = ajb - t1 * (float).5;
/*<       t3 = c1 * ( aje - ajf ) >*/
                            t3 = c1 * (aje - ajf);
/*<       ajh = a(jh+j) >*/
                            ajh = a[jh + j];
/*<       ajf =  ajh >*/
                            ajf = ajh;
/*<       bje = b(je+j) >*/
                            bje = b[je + j];
/*<       bjf = b(jf+j) >*/
                            bjf = b[jf + j];
/*<       u1 = bje + bjf >*/
                            u1 = bje + bjf;
/*<       u2 = bjb - 0.5 * u1 >*/
                            u2 = bjb - u1 * (float).5;
/*<       u3 = c1 * ( bje - bjf ) >*/
                            u3 = c1 * (bje - bjf);
/*<       bjh = b(jh+j) >*/
                            bjh = b[jh + j];
/*<       bjf =  bjh >*/
                            bjf = bjh;
/*<       a(jb+j) = ajb + t1 >*/
                            a[jb + j] = ajb + t1;
/*<       b(jb+j) = bjb + u1 >*/
                            b[jb + j] = bjb + u1;
/*<       a(je+j) = t2 - u3 >*/
                            a[je + j] = t2 - u3;
/*<       b(je+j) = u2 + t3 >*/
                            b[je + j] = u2 + t3;
/*<       a(jh+j) = t2 + u3 >*/
                            a[jh + j] = t2 + u3;
/*<       b(jh+j) = u2 - t3 >*/
                            b[jh + j] = u2 - t3;
/* ---------------------- */
/*<       aji = a(ji+j) >*/
                            aji = a[ji + j];
/*<       t1 = ajf + aji >*/
                            t1 = ajf + aji;
/*<       ajg = a(jg+j) >*/
                            ajg = a[jg + j];
/*<       t2 = ajg - 0.5 * t1 >*/
                            t2 = ajg - t1 * (float).5;
/*<       t3 = c1 * ( ajf - aji ) >*/
                            t3 = c1 * (ajf - aji);
/*<       t1 = ajg + t1 >*/
                            t1 = ajg + t1;
/*<       a(jg+j) = ajc >*/
                            a[jg + j] = ajc;
/*<       bji = b(ji+j) >*/
                            bji = b[ji + j];
/*<       u1 = bjf + bji >*/
                            u1 = bjf + bji;
/*<       bjg = b(jg+j) >*/
                            bjg = b[jg + j];
/*<       u2 = bjg - 0.5 * u1 >*/
                            u2 = bjg - u1 * (float).5;
/*<       u3 = c1 * ( bjf - bji ) >*/
                            u3 = c1 * (bjf - bji);
/*<       u1 = bjg + u1 >*/
                            u1 = bjg + u1;
/*<       b(jg+j) = bjc >*/
                            b[jg + j] = bjc;
/*<       a(jc+j) = t1 >*/
                            a[jc + j] = t1;
/*<       b(jc+j) = u1 >*/
                            b[jc + j] = u1;
/*<       a(jf+j) = t2 - u3 >*/
                            a[jf + j] = t2 - u3;
/*<       b(jf+j) = u2 + t3 >*/
                            b[jf + j] = u2 + t3;
/*<       a(ji+j) = t2 + u3 >*/
                            a[ji + j] = t2 + u3;
/*<       b(ji+j) = u2 - t3 >*/
                            b[ji + j] = u2 - t3;
/*<       j = j + jump >*/
                            j += *jump;
/*<   410 continue >*/
/* L410: */
                        }
/* -----( end of loop across transforms ) */
/*<       ja = ja + jstepx >*/
                        ja += jstepx;