lsqr-test.c

InsightToolkit-1.4.0(有大量的优化算法程序)

/* lsqr-test.f -- translated by f2c (version of 23 April 1993  18:34:30).
   You must link the resulting object file with the libraries:
        -lf2c -lm   (in that order)
*/
/* Peter Vanroose - May 2002 - re-inserted FORTRAN output format, now using printf() */

#include "../f2c.h"
#include "../netlib.h"
#include <stdio.h>
extern double sin(double), cos(double), sqrt(double);

/* Subroutine */ void lstp_(int*m, int*n, int*nduplc, int*npower, double*damp, double*x, double*b, double*d,
                            double*hy, double*hz, double*w, double*acond, double*rnorm);

/* Subroutine */ void aprod_(int*mode, int*m, int*n, double*x, double*y, int*leniw, int*lenrw, int*iw, double*rw);
/* Subroutine */ void hprod_(int*n, double*hz, double*x, double*y);
/* Subroutine */ void aprod1_(int*m, int*n, double*x, double*y, double*d, double*hy, double*hz, double*w);
/* Subroutine */ void aprod2_(int*m, int*n, double*x, double*y, double*d, double*hy, double*hz, double*w);

/* Subroutine */ void test_(int*m, int*n, int*nduplc, int*npower, double*damp);

/* Table of constant values */
static integer c__1 = 1;
static integer c__600 = 600;
static doublereal c_m1 = -1.;
static integer c__2 = 2;
static integer c__40 = 40;
static integer c__4 = 4;
static integer c__80 = 80;

/*     ****************************************************** */
/*                                                            */
/*     WARNING.  Delete the following imitation BLAS routines */
/*               if a genuine BLAS library is available.      */
/*                                                            */
/*     ****************************************************** */
/*                                                            */
/*      SUBROUTINE DCOPY ( N,X,INCX,Y,INCY )                  */
/*      INTEGER            N,INCX,INCY                        */
/*      DOUBLE PRECISION   X(N),Y(N)                          */
/*                                                            */
/*     This may be replaced by the corresponding BLAS routine.*/
/*     The following is a simple version for use with  LSQR.  */
/*                                                            */
/*      DO 10 I = 1, N                                        */
/*         Y(I) = X(I)                                        */
/*   10 CONTINUE                                              */
/*      RETURN                                                */
/*                                                            */
/*     END OF DCOPY                                           */
/*      END                                                   */
/*      DOUBLE PRECISION   FUNCTION DNRM2 ( N,X,INCX )        */
/*      INTEGER            N,INCX                             */
/*      DOUBLE PRECISION   X(N)                               */
/*                                                            */
/*     This may be replaced by the corresponding BLAS routine.*/
/*     The following is a simple version for use with  LSQR.  */
/*                                                            */
/*      INTEGER            I                                  */
/*      DOUBLE PRECISION   D, DSQRT                           */
/*                                                            */
/*      D     = 0.0                                           */
/*      DO 10 I = 1, N                                        */
/*         D    = D + X(I)**2                                 */
/*   10 CONTINUE                                              */
/*      DNRM2 = DSQRT(D)                                      */
/*      RETURN                                                */
/*                                                            */
/*     END OF DNRM2                                           */
/*      END                                                   */
/*      SUBROUTINE DSCAL ( N,A,X,INCX )                       */
/*      INTEGER            N,INCX                             */
/*      DOUBLE PRECISION   A,X(N)                             */
/*                                                            */
/*     This may be replaced by the corresponding BLAS routine.*/
/*     The following is a simple version for use with  LSQR.  */
/*                                                            */
/*      DO 10 I = 1, N                                        */
/*         X(I) = A*X(I)                                      */
/*   10 CONTINUE                                              */
/*      RETURN                                                */
/*                                                            */
/*     END OF DSCAL                                           */
/*      END                                                   */

/* ******************************************************** */
/*                                                          */
/*     These routines are for testing  LSQR.                */
/*                                                          */
/* ******************************************************** */

/* Subroutine */ void aprod_(mode, m, n, x, y, leniw, lenrw, iw, rw)
integer *mode, *m, *n;
doublereal *x, *y;
integer *leniw, *lenrw, *iw;
doublereal *rw;
{
/*     ------------------------------------------------------------------ */
/*     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.                                            */
/*     ------------------------------------------------------------------ */
  (void)leniw;
  (void)lenrw;
  (void)iw;
    if (*mode == 1)
        aprod1_(m, n, x, y, rw, &rw[*n], &rw[*n + *m], &rw[*n + *m + *n]);
    else
        aprod2_(m, n, x, y, rw, &rw[*n], &rw[*n + *m], &rw[*n + *m + *n]);
} /* aprod_ */

/* Subroutine */ void aprod1_(m, n, x, y, d, hy, hz, w)
integer *m, *n;
doublereal *x, *y, *d, *hy, *hz, *w;
{
    /* Local variables */
    static integer i;

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

    hprod_(n, hz, x, w);
    for (i =  0; i < *n; ++i) w[i] = d[i] * w[i];
    for (i = *n; i < *m; ++i) w[i] = 0.;
    hprod_(m, hy, w, w);
    for (i =  0; i < *m; ++i) y[i] += w[i];
} /* aprod1_ */

/* Subroutine */ void aprod2_(m, n, x, y, d, hy, hz, w)
integer *m, *n;
doublereal *x, *y, *d, *hy, *hz, *w;
{
    /* Local variables */
    static integer i;

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

    hprod_(m, hy, y, w);
    for (i = 0; i < *n; ++i) w[i] = d[i] * w[i];
    hprod_(n, hz, w, w);
    for (i = 0; i < *n; ++i) x[i] += w[i];
} /* aprod2_ */

/* Subroutine */ void hprod_(n, hz, x, y)
integer *n;
doublereal *hz, *x, *y;
{
    /* Local variables */
    static integer i;
    static doublereal s;

/*     ------------------------------------------------------------------ */
/*     HPROD  applies a Householder transformation stored in  HZ          */
/*     to get  Y = ( I - 2*HZ*HZ(transpose) ) * X.                        */
/*     ------------------------------------------------------------------ */

    s = 0.0;
    for (i = 0; i < *n; ++i)
        s += hz[i] * x[i];
    s += s;
    for (i = 0; i < *n; ++i)
        y[i] = x[i] - s * hz[i];
} /* hprod_ */

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

    /* Local variables */
    static doublereal alfa, beta;
    static integer i, j;
    static doublereal t;
    static doublereal dampsq, fourpi;

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