gpfa5f.c

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

/* temperton/gpfa5f.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__5 = 5;

/*     fortran version of *gpfa5* - */
/*     radix-5 section of self-sorting, in-place, */
/*        generalized pfa */

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

/*<       subroutine gpfa5f(a,b,trigs,inc,jump,n,mm,lot,isign) >*/
/* Subroutine */ int gpfa5f_(real *a, real *b, real *trigs, integer *inc, 
        integer *jump, integer *n, integer *mm, integer *lot, integer *isign)
{
    /* Initialized data */

    static real sin36 = (float).587785252292473; /* constant */
    static real sin72 = (float).951056516295154; /* constant */
    static real qrt5 = (float).559016994374947; /* 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;
    real s, c1, c2, c3;
    integer n5;
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, u1, u2, u3, u4, u5, u6, u7, u8, 
            u9;
    integer ja, jb, la, jc, jd, nb, je, jf, jg, jh;
    real t10, t11, u10, u11, ax, bx;
    integer mh, kk, ll, ji, jj, jk, mu, nu, jl, jm, jn, jo, jp, jq, jr, js, 
            jt, ju, jv, jw, jx, jy;
    real co1=0, co2=0, co3=0, co4=0, si1=0, si2=0, si3=0, si4=0,
      aja, ajb, ajc, ajd, aje, bjb,
            bje, bjc, bjd, bja, ajf, ajk, bjf, bjk, ajg, ajj, ajh, aji, ajl, 
            ajq, bjg, bjj, bjh, bji, bjl, bjq, ajo, ajm, ajn, ajr, ajw, bjo, 
            bjm, bjn, bjr, bjw, ajt, ajs, ajx, ajp, bjt, bjs, bjx, bjp, ajv, 
            ajy, aju, bjv, bjy, bju;
    integer inq, ink, jjj, ninc, left, nvex, ipass, nblox, jstep, laincl, 
            jstepl, istart, jstepx;

/*<       real a(*), b(*), trigs(*) >*/
/*<       integer inc, jump, n, mm, lot, isign >*/
/*<       real s, ax, bx, c1, c2, c3 >*/
/*<       real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11 >*/
/*<       real u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11 >*/
/*<       real co1, co2, co3, co4, si1, si2, si3, si4 >*/
/*<       real aja, ajb, ajc, ajd, aje, bjb, bje, bjc >*/
/*<       real bjd, bja, ajf, ajk, bjf, bjk, ajg, ajj >*/
/*<       real ajh, aji, ajl, ajq, bjg, bjj, bjh, bji >*/
/*<       real bjl, bjq, ajo, ajm, ajn, ajr, ajw, bjo >*/
/*<       real bjm, bjn, bjr, bjw, ajt, ajs, ajx, ajp >*/
/*<       real bjt, bjs, bjx, bjp, ajv, ajy, aju, bjv >*/
/*<       real bjy, bju >*/
/*<    >*/
    /* 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.      * */
/*     *                                                             * */
/*     *************************************************************** */

/*<       n5 = 5 ** mm >*/
    n5 = pow_ii(&c__5, mm);
/*<       inq = n / n5 >*/
    inq = *n / n5;
/*<       jstepx = (n5-n) * inc >*/
    jstepx = (n5 - *n) * *inc;
/*<       ninc = n * inc >*/
    ninc = *n * *inc;
/*<       ink = inc * inq >*/
    ink = *inc * inq;
/*<       mu = mod(inq,5) >*/
    mu = inq % 5;
/*<       if (isign.eq.-1) mu = 5 - mu >*/
    if (*isign == -1) {
        mu = 5 - mu;
    }

/*<       m = mm >*/
    m = *mm;
/*<       mh = (m+1)/2 >*/
    mh = (m + 1) / 2;
/*<       s = float(isign) >*/
    s = (real) (*isign);
/*<       c1 = qrt5 >*/
    c1 = qrt5;
/*<       c2 = sin72 >*/
    c2 = sin72;
/*<       c3 = sin36 >*/
    c3 = sin36;
/*<       if (mu.eq.2.or.mu.eq.3) then >*/
    if (mu == 2 || mu == 3) {
/*<          c1 = -c1 >*/
        c1 = -c1;
/*<          c2 = sin36 >*/
        c2 = sin36;
/*<          c3 = sin72 >*/
        c3 = sin72;
/*<       endif >*/
    }
/*<       if (mu.eq.3.or.mu.eq.4) c2 = -c2 >*/
    if (mu == 3 || mu == 4) {
        c2 = -c2;
    }
/*<       if (mu.eq.2.or.mu.eq.4) c3 = -c3 >*/
    if (mu == 2 || mu == 4) {
        c3 = -c3;
    }

/*<       nblox = 1 + (lot-1)/lvr >*/
    nblox = (*lot - 1) / lvr + 1;
/*<       left = lot >*/
    left = *lot;
/*<       s = float(isign) >*/
    s = (real) (*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-5 passes */
/*  ----------------------------- */
/*<       do 160 ipass = 1 , mh >*/
        i__2 = mh;
        for (ipass = 1; ipass <= i__2; ++ipass) {
/*<       jstep = (n*inc) / (5*la) >*/
            jstep = *n * *inc / (la * 5);
/*<       jstepl = jstep - ninc >*/
            jstepl = jstep - ninc;
/*<       kk = 0 >*/
            kk = 0;

/*  loop on k */
/*  --------- */
/*<       do 150 k = 0 , jstep-ink , ink >*/
            i__3 = jstep - ink;
            i__4 = ink;
            for (k = 0; i__4 < 0 ? k >= i__3 : k <= i__3; k += i__4) {

/*<       if (k.gt.0) then >*/
                if (k > 0) {
/*<       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];
/*<       co4 = trigs(4*kk+1) >*/
                    co4 = trigs[(kk << 2) + 1];
/*<       si4 = s*trigs(4*kk+2) >*/
                    si4 = s * trigs[(kk << 2) + 2];
/*<       endif >*/
                }

/*  loop along transform */
/*  -------------------- */
/*<       do 140 jjj = k , (n-1)*inc , 5*jstep >*/
                i__5 = (*n - 1) * *inc;
                i__6 = jstep * 5;
                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;
                        }
/*<       je = jd + jstepl >*/
                        je = jd + jstepl;
/*<       if (je.lt.istart) je = je + ninc >*/
                        if (je < istart) {
                            je += ninc;
                        }
/*<       j = 0 >*/
                        j = 0;

/*  loop across transforms */
/*  ---------------------- */
/*<       if (k.eq.0) then >*/
                        if (k == 0) {

/* dir$ ivdep, shortloop */
/*<       do 110 l = 1 , nvex >*/
                            i__8 = nvex;
                            for (l = 1; l <= i__8; ++l) {
/*<       ajb = a(jb+j) >*/
                                ajb = a[jb + j];
/*<       aje = a(je+j) >*/
                                aje = a[je + j];
/*<       t1 = ajb + aje >*/
                                t1 = ajb + aje;
/*<       ajc = a(jc+j) >*/
                                ajc = a[jc + j];
/*<       ajd = a(jd+j) >*/
                                ajd = a[jd + j];
/*<       t2 = ajc + ajd >*/
                                t2 = ajc + ajd;
/*<       t3 = ajb - aje >*/
                                t3 = ajb - aje;
/*<       t4 = ajc - ajd >*/
                                t4 = ajc - ajd;
/*<       t5 = t1 + t2 >*/
                                t5 = t1 + t2;
/*<       t6 = c1 * ( t1 - t2 ) >*/
                                t6 = c1 * (t1 - t2);
/*<       aja = a(ja+j) >*/
                                aja = a[ja + j];
/*<       t7 = aja - 0.25 * t5 >*/
                                t7 = aja - t5 * (float).25;
/*<       a(ja+j) = aja + t5 >*/
                                a[ja + j] = aja + t5;
/*<       t8 = t7 + t6 >*/
                                t8 = t7 + t6;
/*<       t9 = t7 - t6 >*/
                                t9 = t7 - t6;
/*<       t10 = c3 * t3 - c2 * t4 >*/
                                t10 = c3 * t3 - c2 * t4;
/*<       t11 = c2 * t3 + c3 * t4 >*/
                                t11 = c2 * t3 + c3 * t4;
/*<       bjb = b(jb+j) >*/
                                bjb = b[jb + j];
/*<       bje = b(je+j) >*/
                                bje = b[je + j];
/*<       u1 = bjb + bje >*/
                                u1 = bjb + bje;
/*<       bjc = b(jc+j) >*/
                                bjc = b[jc + j];
/*<       bjd = b(jd+j) >*/
                                bjd = b[jd + j];
/*<       u2 = bjc + bjd >*/
                                u2 = bjc + bjd;
/*<       u3 = bjb - bje >*/
                                u3 = bjb - bje;
/*<       u4 = bjc - bjd >*/
                                u4 = bjc - bjd;