lsqr-test.c

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

/*     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                                     */
/*     ------------------------------------------------------------------ */
/*     Make two vectors of norm 1.0 for the Householder transformations.  */
/*     FOURPI  need not be exact.                                         */
/*     ------------------------------------------------------------------ */

    dampsq = *damp * *damp;
    fourpi = 12.566368000000001f;
    alfa = fourpi / *m;
    beta = fourpi / *n;
    for (i = 0; i < *m; ++i)
        hy[i] = sin((i+1) * alfa);
    for (i = 0; i < *n; ++i)
        hz[i] = cos((i+1) * beta);
    alfa = dnrm2_(m, hy, &c__1);
    beta = dnrm2_(n, hz, &c__1);
    d__1 = -1./alfa; dscal_(m, &d__1, hy, &c__1);
    d__1 = -1./beta; dscal_(n, &d__1, hz, &c__1);

/*     ------------------------------------------------------------------ */
/*     Set the diagonal matrix  D.  These are the singular values of  A.  */
/*     ------------------------------------------------------------------ */
    for (i = 0; i < *n; ++i) {
        j = (i + *nduplc) / *nduplc;
        t = (doublereal) (j * *nduplc);
        t /= *n;
        d[i] = pow_di(&t, npower);
    }
    *acond = sqrt((d[*n-1] * d[*n-1] + dampsq) / (d[0] * d[0] + dampsq));
/*     ------------------------------------------------------------------ */
/*     Compute the residual vector, storing it in  B.                     */
/*     It takes the form  HY*( s )                                        */
/*                           ( t )                                        */
/*     where  s  is obtained from  D*s = DAMP**2 * HZ * X                 */
/*     and    t  can be anything.                                         */
/*     ------------------------------------------------------------------ */
    hprod_(n, hz, x, b);
    for (i = 0; i < *n; ++i)
        b[i] = dampsq * b[i] / d[i];
    t = 1.;
    for (i = *n; i < *m; ++i) {
        j = i + 1 - *n;
        b[i] = t * j / *m;
        t = -t;
    }
    hprod_(m, hy, b, b);
/*     ------------------------------------------------------------------ */
/*     Now compute the true  B  =  RESIDUAL  +  A*X.                      */
/*     ------------------------------------------------------------------ */
    d__1 = dnrm2_(m, b, &c__1);
    d__2 = dnrm2_(n, x, &c__1);
    *rnorm = sqrt(d__1 * d__1 + dampsq * (d__2 * d__2));
    aprod1_(m, n, x, b, d, hy, hz, w);
} /* lstp_ */

/* Subroutine */ void test_(m, n, nduplc, npower, damp)
integer *m, *n, *nduplc, *npower;
doublereal *damp;
{
    /* System generated locals */
    doublereal d__1;

    /* Local variables */
    static doublereal atol, btol;
    static integer nout;
    static doublereal b[200];
    static integer j;
    static doublereal acond, u[200], v[100], w[100], x[100];
    static doublereal anorm;
    static doublereal enorm;
    static doublereal rnorm;
    static integer istop;
    static doublereal xnorm, xtrue[100], se[100];
    static integer iw[1];
    static doublereal rw[600], conlim, dampsq;
    static integer itnlim;
    static doublereal arnorm;
    static integer ltotal, itn;

/*     ------------------------------------------------------------------ */
/*     This is an example driver routine for running LSQR.                */
/*     It generates a test problem, solves it, and examines the results.  */
/*     Note that subroutine APROD must be declared EXTERNAL               */
/*     if it is used only in the call to LSQR.                            */

/*     Functions and subroutines                                          */
/*                                                                        */
/*     TESTPROB           APROD                                           */
/*     BLAS               DCOPY, DNRM2, DSCAL                             */
/*     ------------------------------------------------------------------ */

/*     Set the desired solution  XTRUE. */
    for (j = 0; j < *n; ++j) {
        xtrue[j] = (doublereal) (*n - j - 1);
    }
/*     Generate the specified test problem.                               */
/*     The workspace array  IW  is not needed in this application.        */
/*     The workspace array  RW  is used for the following vectors:        */
/*     D(N), HY(M), HZ(N), W(MAX(M,N)).                                   */
/*     The vectors  D, HY, HZ  will define the test matrix A.             */
/*     W is needed for workspace in APROD1 and APROD2.                    */
    ltotal = *n + *m + *n + max(*m,*n);
    if (ltotal > 600) {
        /* Not enough workspace. */
        printf("\n XXX  Insufficient workspace.  The length of  RW  should be at least %d\n", ltotal);
        return;
    }
    lstp_(m, n, nduplc, npower, damp, xtrue, b, rw, &rw[*n], &rw[*n + *m], &rw[*n + *m + *n], &acond, &rnorm);
/*     Solve the problem defined by APROD, DAMP and B.                    */
/*     Copy the rhs vector B into U  (LSQR will overwrite U)              */
/*     and set the other input parameters for LSQR.                       */
    dcopy_(m, b, &c__1, u, &c__1);
    atol = 1e-10f;
    btol = atol;
    conlim = acond * 10.f;
    itnlim = *m + *n + 50;
    printf("\n\n --------------------------------------------------------------------\n");
    printf(" Least-Squares Test Problem      P( %d %d %d %d %12.2e )\n\n", *m,*n,*nduplc,*npower,*damp);
    printf(" Condition no. =%12.4e     Residual function =%17.9e\n", acond, rnorm);
    printf(" --------------------------------------------------------------------\n");
    lsqr_(m, n, aprod_, damp, &c__1, &c__600, iw, rw, u, v, w, x, se, &atol, &btol, &conlim,
          &itnlim, &nout, &istop, &itn, &anorm, &acond, &rnorm, &arnorm, &xnorm);
/*     Examine the results.                                               */
/*     We print the residual norms  RNORM  and  ARNORM  given by LSQR,    */
/*     and then compute their true values in terms of the solution  X     */
/*     obtained by  LSQR.  At least one of them should be small.          */
    dampsq = *damp * *damp;
    printf("\n\n                       Residual norm    Residual norm    Solution norm\n");
    printf("                      (Abar X - bbar)   (Normal eqns)         (X)\n");
    printf(" Estimated by LSQR%17.5e%17.5e%17.5e\n", rnorm, arnorm, xnorm);
/*     Compute  U = A*X - B.                                              */
/*     This is the negative of the usual residual vector.                 */
/*     It will be close to zero only if  B  is a compatible rhs           */
/*     and  X  is close to a solution.                                    */
    dcopy_(m, b, &c__1, u, &c__1);
    dscal_(m, &c_m1, u, &c__1);
    aprod_(&c__1, m, n, x, u, &c__1, &c__600, iw, rw);
/*     Compute  V = A(transpose)*U  +  DAMP**2 * X.                       */
/*     This will be close to zero in all cases                            */
/*     if  X  is close to a solution.                                     */
    dcopy_(n, x, &c__1, v, &c__1);
    dscal_(n, &dampsq, v, &c__1);
    aprod_(&c__2, m, n, v, u, &c__1, &c__600, iw, rw);
/*     Compute the norms associated with  X, U, V. */
    xnorm = dnrm2_(n, x, &c__1);
    d__1 = dnrm2_(m, u, &c__1);
    rnorm = sqrt(d__1 * d__1 + dampsq * xnorm * xnorm);
    arnorm = dnrm2_(n, v, &c__1);
    printf(" Computed from  X %17.5e%17.5e%17.5e\n", rnorm, arnorm, xnorm);
/*     Print the solution and standard error estimates from  LSQR. */
    printf("\n\n Solution  X\n");
    for (j = 0; j < *n; ++j)
        printf("%6d%14.6g\n", j+1, x[j]);
    printf("\n\n Standard errors  SE\n");
    for (j = 0; j < *n; ++j)
        printf("%6d%14.6g\n", j+1, se[j]);
    printf("\n");
/*     Print a clue about whether the solution looks OK. */
    for (j = 0; j < *n; ++j)
        w[j] = x[j] - xtrue[j];
    enorm = dnrm2_(n, w, &c__1) / (dnrm2_(n, xtrue, &c__1) + 1.);
    if (enorm <= 1e-5)
        printf("\n LSQR  appears to be successful.     Relative error in  X  =%10.2e\n", enorm);
    else
        printf("\n LSQR  appears to have failed.       Relative error in  X  =%10.2e\n", enorm);
} /* test_ */

/*     ------------- */
/*     Main program. */
/*     ------------- */
int main(void)
{
    static doublereal zero = 0.f;
/*  static doublereal damp1 = .1f; */
    static doublereal damp2 = .01f;
/*  static doublereal damp3 = .001f; */
/*  static doublereal damp4 = 1e-4f; */
    test_(&c__1, &c__1, &c__1, &c__1, &zero);
    test_(&c__2, &c__1, &c__1, &c__1, &zero);
    test_(&c__40, &c__40, &c__4, &c__4, &zero);
    test_(&c__40, &c__40, &c__4, &c__4, &damp2);
    test_(&c__80, &c__40, &c__4, &c__4, &damp2);
    return 0;
} /* End of main program for testing LSQR */