dgpfa3f.c

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

/* temperton/dgpfa3f.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__3 = 3;

/*     fortran version of *dgpfa3* - */
/*     radix-3 section of self-sorting, in-place */
/*        generalized PFA */

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

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

    static doublereal sin60 = .866025403784437; /* constant */
    static integer lvr = 128; /* 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 *);

    /* Local variables */
    integer j, k, l, m;
    doublereal s, c1;
    integer n3;
    doublereal t1, t2, t3, u1, u2, u3;
    integer ja, jb, la, jc, jd, nb, je, jf, jg, jh, mh, kk, ji, ll, mu, nu;
    doublereal co1, co2, si1, si2, aja, ajb, ajc, bjb, bjc, bja, ajd, bjd, 
            aje, ajf, ajh, bje, bjf, bjh, aji, ajg, bji, bjg;
    integer jjj, ink, inq, 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, c1, t1, t2, t3, u1, u2, u3, co1, co2 >*/
/*<       double precision si1, si2, aja, ajb, ajc, bjb, bjc, bja, ajd, bjd >*/
/*<       double precision aje, ajf, ajh, bje, bjf, bjh, aji, ajg, bji, bjg >*/
/*<       data sin60/0.866025403784437d+0/ >*/
    /* Parameter adjustments */
    --trigs;
    --b;
    --a;

    /* Function Body */
/*<       data lvr/128/ >*/

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

/*<       n3 = 3**mm >*/
    n3 = pow_ii(&c__3, mm);
/*<       inq = n/n3 >*/
    inq = *n / n3;
/*<       jstepx = (n3-n) * inc >*/
    jstepx = (n3 - *n) * *inc;
/*<       ninc = n * inc >*/
    ninc = *n * *inc;
/*<       ink = inc * inq >*/
    ink = *inc * inq;
/*<       mu = mod(inq,3) >*/
    mu = inq % 3;
/*<       if (isign.eq.-1) mu = 3-mu >*/
    if (*isign == -1) {
        mu = 3 - mu;
    }
/*<       m = mm >*/
    m = *mm;
/*<       mh = (m+1)/2 >*/
    mh = (m + 1) / 2;
/*<       s = dfloat(isign) >*/
    s = (doublereal) (*isign);
/*<       c1 = sin60 >*/
    c1 = sin60;
/*<       if (mu.eq.2) c1 = -c1 >*/
    if (mu == 2) {
        c1 = -c1;
    }

/*<       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-3 passes */
/*  ----------------------------- */
/*<       do 160 ipass = 1 , mh >*/
        i__2 = mh;
        for (ipass = 1; ipass <= i__2; ++ipass) {
/*<       jstep = (n*inc) / (3*la) >*/
            jstep = *n * *inc / (la * 3);
/*<       jstepl = jstep - ninc >*/
            jstepl = jstep - ninc;

/*  k = 0 loop (no twiddle factors) */
/*  ------------------------------- */
/*<       do 120 jjj = 0 , (n-1)*inc , 3*jstep >*/
            i__3 = (*n - 1) * *inc;
            i__4 = jstep * 3;
            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;
                    }
/*<       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) {
/*<       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) = t2 - u3 >*/
                        a[jb + j] = t2 - u3;
/*<       b(jb+j) = u2 + t3 >*/
                        b[jb + j] = u2 + t3;
/*<       a(jc+j) = t2 + u3 >*/
                        a[jc + j] = t2 + u3;
/*<       b(jc+j) = u2 - t3 >*/
                        b[jc + 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 n3 = 3 */
/*  ------------------ */
/*<       if (n3.eq.3) go to 490 >*/
            if (n3 == 3) {
                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];

/*  loop along transform */
/*  -------------------- */
/*<       do 140 jjj = k , (n-1)*inc , 3*jstep >*/
                i__5 = (*n - 1) * *inc;
                i__6 = jstep * 3;
                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;
                        }
/*<       j = 0 >*/
                        j = 0;

/*  loop across transforms */
/*  ---------------------- */