lsqr-test.c

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

/* tests/lsqr-test.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
*/

#include "v3p_netlib.h"

#include <stdio.h>

/* Table of constant values */

static integer c__1 = 1;
static integer c__600 = 600;
static doublereal c_b53 = -1.;
static integer c__2 = 2;
static integer c__40 = 40;
static integer c__4 = 4;
static integer c__80 = 80;

/* ******************************************************** */

/*     These routines are for testing  LSQR. */

/* ******************************************************** */
/*<       SUBROUTINE APROD ( MODE, M, N, X, Y, LENIW, LENRW, IW, RW ) >*/
/* Subroutine */ int aprod_(integer *mode, integer *m, integer *n, doublereal 
        *x, doublereal *y, integer *leniw, integer *lenrw, integer *iw, 
        doublereal *rw, void* userdata)
{
    integer locd, locw, lochy, lochz;
    extern /* Subroutine */ int aprod1_(integer *, integer *, doublereal *, 
            doublereal *, doublereal *, doublereal *, doublereal *, 
            doublereal *), aprod2_(integer *, integer *, doublereal *, 
            doublereal *, doublereal *, doublereal *, doublereal *, 
            doublereal *);
    (void)leniw;
    (void)lenrw;
    (void)userdata;

/*<       INTEGER            MODE, M, N, LENIW, LENRW >*/
/*<       INTEGER            IW(LENIW) >*/
/*<       DOUBLE PRECISION   X(N), Y(M), RW(LENRW) >*/
/*     ------------------------------------------------------------------ */
/*     This is the matrix-vector product routine required by  LSQR */
/*     for a test matrix of the form  A = HY*D*HZ.  The quantities */
/*     defining D, HY, HZ are in the work array RW, followed by a */
/*     work array W.  These are passed to APROD1 and APROD2 in order to */
/*     make the code readable. */
/*     ------------------------------------------------------------------ */
/*<       INTEGER            LOCD, LOCHY, LOCHZ, LOCW >*/
/*<       LOCD   = 1 >*/
    /* Parameter adjustments */
    --y;
    --x;
    --iw;
    --rw;

    /* Function Body */
    locd = 1;
/*<       LOCHY  = LOCD  + N >*/
    lochy = locd + *n;
/*<       LOCHZ  = LOCHY + M >*/
    lochz = lochy + *m;
/*<       LOCW   = LOCHZ + N >*/
    locw = lochz + *n;
/*<    >*/
    if (*mode == 1) {
        aprod1_(m, n, &x[1], &y[1], &rw[locd], &rw[lochy], &rw[lochz], &rw[
                locw]);
    }
/*<    >*/
    if (*mode != 1) {
        aprod2_(m, n, &x[1], &y[1], &rw[locd], &rw[lochy], &rw[lochz], &rw[
                locw]);
    }
/*     End of APROD */
/*<       END >*/
    return 0;
} /* aprod_ */

/*<       SUBROUTINE APROD1( M, N, X, Y, D, HY, HZ, W ) >*/
/* Subroutine */ int aprod1_(integer *m, integer *n, doublereal *x, 
        doublereal *y, doublereal *d__, doublereal *hy, doublereal *hz, 
        doublereal *w)
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    integer i__;
    extern /* Subroutine */ int hprod_(integer *, doublereal *, doublereal *, 
            doublereal *);

/*<       INTEGER            M, N >*/
/*<       DOUBLE PRECISION   X(N), Y(M), D(N), HY(M), HZ(N), W(M) >*/
/*     ------------------------------------------------------------------ */
/*     APROD1  computes  Y = Y + A*X  for subroutine APROD, */
/*     where A is a test matrix of the form  A = HY*D*HZ, */
/*     and the latter matrices HY, D, HZ are represented by */
/*     input vectors with the same name. */
/*     ------------------------------------------------------------------ */
/*<       INTEGER            I >*/
/*<       DOUBLE PRECISION   ZERO >*/
/*<       PARAMETER        ( ZERO = 0.0 ) >*/
/*<       CALL HPROD ( N, HZ, X, W ) >*/
    /* Parameter adjustments */
    --w;
    --hy;
    --y;
    --hz;
    --d__;
    --x;

    /* Function Body */
    hprod_(n, &hz[1], &x[1], &w[1]);
/*<       DO 100 I = 1, N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          W(I)  = D(I) * W(I) >*/
        w[i__] = d__[i__] * w[i__];
/*<   100 CONTINUE >*/
/* L100: */
    }
/*<       DO 200 I = N + 1, M >*/
    i__1 = *m;
    for (i__ = *n + 1; i__ <= i__1; ++i__) {
/*<          W(I)  = ZERO >*/
        w[i__] = 0.;
/*<   200 CONTINUE >*/
/* L200: */
    }
/*<       CALL HPROD ( M, HY, W, W ) >*/
    hprod_(m, &hy[1], &w[1], &w[1]);
/*<       DO 600 I = 1, M >*/
    i__1 = *m;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          Y(I)  = Y(I) + W(I) >*/
        y[i__] += w[i__];
/*<   600 CONTINUE >*/
/* L600: */
    }
/*     End of APROD1 */
/*<       END >*/
    return 0;
} /* aprod1_ */

/*<       SUBROUTINE APROD2( M, N, X, Y, D, HY, HZ, W ) >*/
/* Subroutine */ int aprod2_(integer *m, integer *n, doublereal *x, 
        doublereal *y, doublereal *d__, doublereal *hy, doublereal *hz, 
        doublereal *w)
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    integer i__;
    extern /* Subroutine */ int hprod_(integer *, doublereal *, doublereal *, 
            doublereal *);

/*<       INTEGER            M, N >*/
/*<       DOUBLE PRECISION   X(N), Y(M), D(N), HY(M), HZ(N), W(M) >*/
/*     ------------------------------------------------------------------ */
/*     APROD2  computes  X = X + A(T)*Y  for subroutine APROD, */
/*     where  A  is a test matrix of the form  A = HY*D*HZ, */
/*     and the latter matrices  HY, D, HZ  are represented by */
/*     input vectors with the same name. */
/*     ------------------------------------------------------------------ */
/*<       INTEGER            I >*/
/*<       CALL HPROD ( M, HY, Y, W ) >*/
    /* Parameter adjustments */
    --w;
    --hy;
    --y;
    --hz;
    --d__;
    --x;

    /* Function Body */
    hprod_(m, &hy[1], &y[1], &w[1]);
/*<       DO 100 I = 1, N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          W(I)  = D(I)*W(I) >*/
        w[i__] = d__[i__] * w[i__];
/*<   100 CONTINUE >*/
/* L100: */
    }
/*<       CALL HPROD ( N, HZ, W, W ) >*/
    hprod_(n, &hz[1], &w[1], &w[1]);
/*<       DO 600 I = 1, N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          X(I)  = X(I) + W(I) >*/
        x[i__] += w[i__];
/*<   600 CONTINUE >*/
/* L600: */
    }
/*     End of APROD2 */
/*<       END >*/
    return 0;
} /* aprod2_ */

/*<       SUBROUTINE HPROD ( N, HZ, X, Y ) >*/
/* Subroutine */ int hprod_(integer *n, doublereal *hz, doublereal *x, 
        doublereal *y)
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    integer i__;
    doublereal s;

/*<       INTEGER            N >*/
/*<       DOUBLE PRECISION   HZ(N), X(N), Y(N) >*/
/*     ------------------------------------------------------------------ */
/*     HPROD  applies a Householder transformation stored in  HZ */
/*     to get  Y = ( I - 2*HZ*HZ(transpose) ) * X. */
/*     ------------------------------------------------------------------ */
/*<       INTEGER            I >*/
/*<       DOUBLE PRECISION   S >*/
/*<       S      = 0.0 >*/
    /* Parameter adjustments */
    --y;
    --x;
    --hz;

    /* Function Body */
    s = (float)0.;
/*<       DO 100 I = 1, N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          S     = HZ(I) * X(I)  +  S >*/
        s = hz[i__] * x[i__] + s;
/*<   100 CONTINUE >*/
/* L100: */
    }
/*<       S      = S + S >*/
    s += s;
/*<       DO 200 I = 1, N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          Y(I)  = X(I)  -  S * HZ(I) >*/
        y[i__] = x[i__] - s * hz[i__];
/*<   200 CONTINUE >*/
/* L200: */
    }
/*     End of HPROD */
/*<       END >*/
    return 0;
} /* hprod_ */

/*<    >*/
/* Subroutine */ int lstp_(integer *m, integer *n, integer *nduplc, integer *
        npower, doublereal *damp, doublereal *x, doublereal *b, doublereal *
        d__, doublereal *hy, doublereal *hz, doublereal *w, doublereal *acond,
         doublereal *rnorm)
{
    /* System generated locals */
    integer i__1;
    doublereal d__1, d__2;

    /* Builtin functions */
    double sin(doublereal), cos(doublereal), pow_di(doublereal *, integer *), 
            sqrt(doublereal);

    /* Local variables */
    integer i__, j;
    doublereal t, alfa, beta;
    extern doublereal dnrm2_(integer *, doublereal *, integer *);
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
            integer *), hprod_(integer *, doublereal *, doublereal *, 
            doublereal *), aprod1_(integer *, integer *, doublereal *, 
            doublereal *, doublereal *, doublereal *, doublereal *, 
            doublereal *);
    doublereal dampsq, fourpi;

/*<       INTEGER            M, N, MAXMN, NDUPLC, NPOWER >*/
/*<       DOUBLE PRECISION   DAMP, ACOND, RNORM >*/
/*<       DOUBLE PRECISION   B(M), X(N), D(N), HY(M), HZ(N), W(M) >*/
/*     ------------------------------------------------------------------ */
/*     LSTP  generates a sparse least-squares test problem of the form */

/*                (   A    )*X = ( B ) */
/*                ( DAMP*I )     ( 0 ) */

/*     having a specified solution X.  The matrix A is constructed */
/*     in the form  A = HY*D*HZ,  where D is an M by N diagonal matrix, */
/*     and HY and HZ are Householder transformations. */

/*     The first 6 parameters are input to LSTP.  The remainder are */
/*     output.  LSTP is intended for use with M .GE. N. */


/*     Functions and subroutines */

/*     TESTPROB           APROD1, HPROD */
/*     BLAS               DNRM2 */
/*     ------------------------------------------------------------------ */
/*     Intrinsics and local variables */
/*<       INTRINSIC          COS,  SIN, SQRT >*/
/*<       INTEGER            I, J >*/
/*<       DOUBLE PRECISION   DNRM2 >*/
/*<       DOUBLE PRECISION   ALFA, BETA, DAMPSQ, FOURPI, T >*/
/*<       DOUBLE PRECISION   ZERO,        ONE >*/
/*<       PARAMETER        ( ZERO = 0.0,  ONE = 1.0 ) >*/
/*     ------------------------------------------------------------------ */
/*     Make two vectors of norm 1.0 for the Householder transformations. */
/*     FOURPI  need not be exact. */
/*     ------------------------------------------------------------------ */
/*<       DAMPSQ = DAMP**2 >*/
    /* Parameter adjustments */
    --w;
    --hy;
    --b;
    --hz;
    --d__;
    --x;

    /* Function Body */
/* Computing 2nd power */
    d__1 = *damp;
    dampsq = d__1 * d__1;
/*<       FOURPI = 4.0 * 3.141592 >*/
    fourpi = (float)12.566368000000001;
/*<       ALFA   = FOURPI / M >*/
    alfa = fourpi / *m;
/*<       BETA   = FOURPI / N >*/
    beta = fourpi / *n;
/*<       DO 100 I = 1, M >*/
    i__1 = *m;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          HY(I) = SIN( I * ALFA ) >*/
        hy[i__] = sin(i__ * alfa);
/*<   100 CONTINUE >*/
/* L100: */
    }
/*<       DO 200 I = 1, N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          HZ(I) = COS( I * BETA ) >*/
        hz[i__] = cos(i__ * beta);
/*<   200 CONTINUE                 >*/
/* L200: */
    }
/*<       ALFA   = DNRM2 ( M, HY, 1 ) >*/
    alfa = dnrm2_(m, &hy[1], &c__1);
/*<       BETA   = DNRM2 ( N, HZ, 1 ) >*/
    beta = dnrm2_(n, &hz[1], &c__1);
/*<       CALL DSCAL ( M, (- ONE / ALFA), HY, 1 ) >*/
    d__1 = -1. / alfa;
    dscal_(m, &d__1, &hy[1], &c__1);
/*<       CALL DSCAL ( N, (- ONE / BETA), HZ, 1 ) >*/
    d__1 = -1. / beta;
    dscal_(n, &d__1, &hz[1], &c__1);

/*     ------------------------------------------------------------------ */
/*     Set the diagonal matrix  D.  These are the singular values of  A. */
/*     ------------------------------------------------------------------ */