dgpfa2f.c

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

/* temperton/dgpfa2f.f -- translated by f2c (version 20050501).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
        or, if you install libf2c.a in a standard place, with -lf2c -lm
        -- in that order, at the end of the command line, as in
                cc *.o -lf2c -lm
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

                http://www.netlib.org/f2c/libf2c.zip
*/

#ifdef __cplusplus
extern "C" {
#endif
#include "v3p_netlib.h"

/* Table of constant values */

static integer c__2 = 2;

/*     fortran version of *dgpfa2* - */
/*     radix-2 section of self-sorting, in-place, generalized pfa */
/*     central radix-2 and radix-8 passes included */
/*      so that transform length can be any power of 2 */

/* ------------------------------------------------------------------- */

/*<       subroutine dgpfa2f(a,b,trigs,inc,jump,n,mm,lot,isign) >*/
/* Subroutine */ int dgpfa2f_(doublereal *a, doublereal *b, doublereal *trigs,
         integer *inc, integer *jump, integer *n, integer *mm, integer *lot, 
        integer *isign)
{
    /* Initialized data */

    static integer lvr = 1024; /* constant */

    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5, i__6, i__7, i__8, i__9, i__10;

    /* Builtin functions */
    integer pow_ii(integer *, integer *);
    double sqrt(doublereal);

    /* Local variables */
    integer j, k, l, m=0;
    doublereal s, c1, c2, c3;
    integer m2, n2;
    doublereal t0, t1, t2, t3, u0, u1, u2, u3;
    integer m8, ja, jb, la, jc, jd, nb, je, jf, jg, jh, mh, kk, ji, ll, jj, 
            jk, jl, jm, jn, jo, jp, mu, nu;
    doublereal ss, co1, co2, co3, co4, co5, co6, co7, si1, si2, si3, si4, si5,
             si6, si7, aja, ajb, ajc, ajd, bja, bjc, bjb, bjd, aje, ajg, ajf, 
            ajh, bje, bjg, bjf, bjh, aji, bjm, ajj, bjj, ajk, ajl, bji, bjk, 
            ajo, bjl, bjo, ajm, ajn, ajp, bjn, bjp;
    integer inq, ink, jjj, ninc, left, nvex, ipass, nblox, jstep, laincl, 
            jstepl, istart, jstepx;

/*<       double precision a(*), b(*), trigs(*) >*/
/*<       integer inc, jump, n, mm, lot, isign >*/
/*<       double precision s, ss, c1, c2, c3, t0, t1, t2, t3, u0, u1, u2, u3 >*/
/*<       double precision co1, co2, co3, co4, co5, co6, co7 >*/
/*<       double precision si1, si2, si3, si4, si5, si6, si7 >*/
/*<       double precision aja, ajb, ajc, ajd, bja, bjc, bjb, bjd >*/
/*<       double precision aje, ajg, ajf, ajh, bje, bjg, bjf, bjh, aji >*/
/*<       double precision bjm, ajj, bjj, ajk, ajl, bji, bjk >*/
/*<       double precision ajo, bjl, bjo, ajm, ajn, ajp, bjn, bjp >*/
/*<       data lvr/1024/ >*/
    /* Parameter adjustments */
    --trigs;
    --b;
    --a;

    /* Function Body */

/*     *************************************************************** */
/*     *                                                             * */
/*     *  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.      * */
/*     *                                                             * */
/*     *************************************************************** */

/*<       n2 = 2**mm >*/
    n2 = pow_ii(&c__2, mm);
/*<       inq = n/n2 >*/
    inq = *n / n2;
/*<       jstepx = (n2-n) * inc >*/
    jstepx = (n2 - *n) * *inc;
/*<       ninc = n * inc >*/
    ninc = *n * *inc;
/*<       ink = inc * inq >*/
    ink = *inc * inq;

/*<       m2 = 0 >*/
    m2 = 0;
/*<       m8 = 0 >*/
    m8 = 0;
/*<       if (mod(mm,2).eq.0) then >*/
    if (*mm % 2 == 0) {
/*<          m = mm/2 >*/
        m = *mm / 2;
/*<       else if (mod(mm,4).eq.1) then >*/
    } else if (*mm % 4 == 1) {
/*<          m = (mm-1)/2 >*/
        m = (*mm - 1) / 2;
/*<          m2 = 1 >*/
        m2 = 1;
/*<       else if (mod(mm,4).eq.3) then >*/
    } else if (*mm % 4 == 3) {
/*<          m = (mm-3)/2 >*/
        m = (*mm - 3) / 2;
/*<          m8 = 1 >*/
        m8 = 1;
/*<       endif >*/
    }
/*<       mh = (m+1)/2 >*/
    mh = (m + 1) / 2;

/*<       nblox = 1 + (lot-1)/lvr >*/
    nblox = (*lot - 1) / lvr + 1;
/*<       left = lot >*/
    left = *lot;
/*<       s = dfloat(isign) >*/
    s = (doublereal) (*isign);
/*<       istart = 1 >*/
    istart = 1;

/*  loop on blocks of lvr transforms */
/*  -------------------------------- */
/*<       do 500 nb = 1 , nblox >*/
    i__1 = nblox;
    for (nb = 1; nb <= i__1; ++nb) {

/*<       if (left.le.lvr) then >*/
        if (left <= lvr) {
/*<          nvex = left >*/
            nvex = left;
/*<       else if (left.lt.(2*lvr)) then >*/
        } else if (left < lvr << 1) {
/*<          nvex = left/2 >*/
            nvex = left / 2;
/*<          nvex = nvex + mod(nvex,2) >*/
            nvex += nvex % 2;
/*<       else >*/
        } else {
/*<          nvex = lvr >*/
            nvex = lvr;
/*<       endif >*/
        }
/*<       left = left - nvex >*/
        left -= nvex;

/*<       la = 1 >*/
        la = 1;

/*  loop on type I radix-4 passes */
/*  ----------------------------- */
/*<       mu = mod(inq,4) >*/
        mu = inq % 4;
/*<       if (isign.eq.-1) mu = 4 - mu >*/
        if (*isign == -1) {
            mu = 4 - mu;
        }
/*<       ss = 1.0 >*/
        ss = (float)1.;
/*<       if (mu.eq.3) ss = -1.0 >*/
        if (mu == 3) {
            ss = (float)-1.;
        }

/*<       if (mh.eq.0) go to 200 >*/
        if (mh == 0) {
            goto L200;
        }

/*<       do 160 ipass = 1 , mh >*/
        i__2 = mh;
        for (ipass = 1; ipass <= i__2; ++ipass) {
/*<       jstep = (n*inc) / (4*la) >*/
            jstep = *n * *inc / (la << 2);
/*<       jstepl = jstep - ninc >*/
            jstepl = jstep - ninc;

/*  k = 0 loop (no twiddle factors) */
/*  ------------------------------- */
/*<       do 120 jjj = 0 , (n-1)*inc , 4*jstep >*/
            i__3 = (*n - 1) * *inc;
            i__4 = jstep << 2;
            for (jjj = 0; i__4 < 0 ? jjj >= i__3 : jjj <= i__3; jjj += i__4) {
/*<       ja = istart + jjj >*/
                ja = istart + jjj;

/*     "transverse" loop */
/*     ----------------- */
/*<       do 115 nu = 1 , inq >*/
                i__5 = inq;
                for (nu = 1; nu <= i__5; ++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 = jc + jstepl >*/
                    jd = jc + jstepl;
/*<       if (jd.lt.istart) jd = jd + ninc >*/
                    if (jd < istart) {
                        jd += ninc;
                    }
/*<       j = 0 >*/
                    j = 0;

/*  loop across transforms */
/*  ---------------------- */
/* dir$ ivdep, shortloop */
/*<       do 110 l = 1 , nvex >*/
                    i__6 = nvex;
                    for (l = 1; l <= i__6; ++l) {
/*<       aja = a(ja+j) >*/
                        aja = a[ja + j];
/*<       ajc = a(jc+j) >*/
                        ajc = a[jc + j];
/*<       t0 = aja + ajc >*/
                        t0 = aja + ajc;
/*<       t2 = aja - ajc >*/
                        t2 = aja - ajc;
/*<       ajb = a(jb+j) >*/
                        ajb = a[jb + j];
/*<       ajd = a(jd+j) >*/
                        ajd = a[jd + j];
/*<       t1 = ajb + ajd >*/
                        t1 = ajb + ajd;
/*<       t3 = ss * ( ajb - ajd ) >*/
                        t3 = ss * (ajb - ajd);
/*<       bja = b(ja+j) >*/
                        bja = b[ja + j];
/*<       bjc = b(jc+j) >*/
                        bjc = b[jc + j];
/*<       u0 = bja + bjc >*/
                        u0 = bja + bjc;
/*<       u2 = bja - bjc >*/
                        u2 = bja - bjc;
/*<       bjb = b(jb+j) >*/
                        bjb = b[jb + j];
/*<       bjd = b(jd+j) >*/
                        bjd = b[jd + j];
/*<       u1 = bjb + bjd >*/
                        u1 = bjb + bjd;
/*<       u3 = ss * ( bjb - bjd ) >*/
                        u3 = ss * (bjb - bjd);
/*<       a(ja+j) = t0 + t1 >*/
                        a[ja + j] = t0 + t1;
/*<       a(jc+j) = t0 - t1 >*/
                        a[jc + j] = t0 - t1;
/*<       b(ja+j) = u0 + u1 >*/
                        b[ja + j] = u0 + u1;
/*<       b(jc+j) = u0 - u1 >*/
                        b[jc + j] = u0 - u1;
/*<       a(jb+j) = t2 - u3 >*/
                        a[jb + j] = t2 - u3;
/*<       a(jd+j) = t2 + u3 >*/
                        a[jd + j] = t2 + u3;
/*<       b(jb+j) = u2 + t3 >*/
                        b[jb + j] = u2 + t3;
/*<       b(jd+j) = u2 - t3 >*/
                        b[jd + j] = u2 - t3;
/*<       j = j + jump >*/
                        j += *jump;
/*<   110 continue >*/
/* L110: */
                    }
/*<       ja = ja + jstepx >*/
                    ja += jstepx;
/*<       if (ja.lt.istart) ja = ja + ninc >*/
                    if (ja < istart) {
                        ja += ninc;
                    }
/*<   115 continue >*/
/* L115: */
                }
/*<   120 continue >*/
/* L120: */
            }

/*  finished if n2 = 4 */
/*  ------------------ */
/*<       if (n2.eq.4) go to 490 >*/
            if (n2 == 4) {
                goto L490;
            }
/*<       kk = 2 * la >*/
            kk = la << 1;

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

/*  loop along transform */
/*  -------------------- */
/*<       do 140 jjj = k , (n-1)*inc , 4*jstep >*/
                i__5 = (*n - 1) * *inc;
                i__6 = jstep << 2;
                for (jjj = k; i__6 < 0 ? jjj >= i__5 : jjj <= i__5; jjj += 
                        i__6) {
/*<       ja = istart + jjj >*/
                    ja = istart + jjj;

/*     "transverse" loop */
/*     ----------------- */
/*<       do 135 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 = jc + jstepl >*/
                        jd = jc + jstepl;
/*<       if (jd.lt.istart) jd = jd + ninc >*/
                        if (jd < istart) {
                            jd += ninc;
                        }
/*<       j = 0 >*/
                        j = 0;

/*  loop across transforms */
/*  ---------------------- */
/* dir$ ivdep,shortloop */
/*<       do 130 l = 1 , nvex >*/
                        i__8 = nvex;
                        for (l = 1; l <= i__8; ++l) {
/*<       aja = a(ja+j) >*/
                            aja = a[ja + j];
/*<       ajc = a(jc+j) >*/
                            ajc = a[jc + j];
/*<       t0 = aja + ajc >*/
                            t0 = aja + ajc;
/*<       t2 = aja - ajc >*/
                            t2 = aja - ajc;
/*<       ajb = a(jb+j) >*/