lsqr-test.c

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

/*<       DO 300 I = 1, N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          J     = (I - 1 + NDUPLC) / NDUPLC >*/
        j = (i__ - 1 + *nduplc) / *nduplc;
/*<          T     =  J * NDUPLC >*/
        t = (doublereal) (j * *nduplc);
/*<          T     =  T / N >*/
        t /= *n;
/*<          D(I)  =  T**NPOWER >*/
        d__[i__] = pow_di(&t, npower);
/*<   300 CONTINUE >*/
/* L300: */
    }
/*<       ACOND  = SQRT( (D(N)**2 + DAMPSQ) / (D(1)**2 + DAMPSQ) ) >*/
/* Computing 2nd power */
    d__1 = d__[*n];
/* Computing 2nd power */
    d__2 = d__[1];
    *acond = sqrt((d__1 * d__1 + dampsq) / (d__2 * d__2 + 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. */
/*     ------------------------------------------------------------------ */
/*<       CALL HPROD ( N, HZ, X, B ) >*/
    hprod_(n, &hz[1], &x[1], &b[1]);
/*<       DO 500 I = 1, N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          B(I)  = DAMPSQ * B(I) / D(I) >*/
        b[i__] = dampsq * b[i__] / d__[i__];
/*<   500 CONTINUE >*/
/* L500: */
    }
/*<       T      = ONE >*/
    t = 1.;
/*<       DO 600 I =   N + 1, M >*/
    i__1 = *m;
    for (i__ = *n + 1; i__ <= i__1; ++i__) {
/*<          J     =   I - N >*/
        j = i__ - *n;
/*<          B(I)  =  (T * J) / M >*/
        b[i__] = t * j / *m;
/*<          T     = - T >*/
        t = -t;
/*<   600 CONTINUE >*/
/* L600: */
    }
/*<       CALL HPROD ( M, HY, B, B ) >*/
    hprod_(m, &hy[1], &b[1], &b[1]);
/*     ------------------------------------------------------------------ */
/*     Now compute the true  B  =  RESIDUAL  +  A*X. */
/*     ------------------------------------------------------------------ */
/*<    >*/
/* Computing 2nd power */
    d__1 = dnrm2_(m, &b[1], &c__1);
/* Computing 2nd power */
    d__2 = dnrm2_(n, &x[1], &c__1);
    *rnorm = sqrt(d__1 * d__1 + dampsq * (d__2 * d__2));
/*<       CALL APROD1( M, N, X, B, D, HY, HZ, W ) >*/
    aprod1_(m, n, &x[1], &b[1], &d__[1], &hy[1], &hz[1], &w[1]);
/*     End of LSTP */
/*<       END >*/
    return 0;
} /* lstp_ */

/*<       SUBROUTINE TEST  ( M, N, NDUPLC, NPOWER, DAMP ) >*/
/* Subroutine */ int test_(integer *m, integer *n, integer *nduplc, integer *
        npower, doublereal *damp)
{
    /* System generated locals */
    integer i__1;
    doublereal d__1, d__2;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    doublereal b[200];
    integer j;
    doublereal u[200], v[100], w[100], x[100], se[100];
    integer iw[1];
    doublereal rw[600];
    integer itn, locd;
    doublereal atol, btol, etol;
    integer locw;
    extern int lstp_(integer *, integer *, integer *
            , integer *, doublereal *, doublereal *, doublereal *, doublereal 
            *, doublereal *, doublereal *, doublereal *, doublereal *, 
            doublereal *);
    integer nout;
    extern doublereal dnrm2_(integer *, doublereal *, integer *);
    doublereal acond;
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
            integer *), aprod_(integer *, integer *, integer *, doublereal *, 
            doublereal *, integer *, integer *, integer *, doublereal *,
            void*);
    doublereal anorm;
    integer lochy;
    doublereal enorm;
    integer lochz;
    extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 
            doublereal *, integer *);
    doublereal rnorm;
    integer istop;
    doublereal xnorm, xtrue[100], conlim, dampsq;
    integer itnlim;
    doublereal arnorm;
    integer ltotal;

/*<       INTEGER            M, N, NDUPLC, NPOWER >*/
/*<       DOUBLE PRECISION   DAMP >*/
/*     ------------------------------------------------------------------ */
/*     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 */
/*     ------------------------------------------------------------------ */
/*     Intrinsics and local variables */
/*<       INTRINSIC          MAX, SQRT >*/
/*<       EXTERNAL           APROD >*/
/*<       INTEGER            ISTOP, ITNLIM, J, NOUT >*/
/*<       DOUBLE PRECISION   DNRM2 >*/
/*<       PARAMETER        ( MAXM = 200,  MAXN = 100 ) >*/
/*<    >*/
/*<    >*/
/*<       PARAMETER        ( LENIW = 1,  LENRW = 600 ) >*/
/*<       INTEGER            IW(LENIW) >*/
/*<       DOUBLE PRECISION   RW(LENRW) >*/
/*<       INTEGER            LOCD, LOCHY, LOCHZ, LOCW, LTOTAL >*/
/*<       DOUBLE PRECISION   ONE >*/
/*<       PARAMETER        ( ONE = 1.0 ) >*/
/*<       CHARACTER*34       LINE >*/
/*<    >*/
/*     Set the desired solution  XTRUE. */
/*<       DO 100 J = 1, N >*/
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
/*<          XTRUE(J) = N - J >*/
        xtrue[j - 1] = (doublereal) (*n - j);
/*<   100 CONTINUE >*/
/* L100: */
    }
/*     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. */
/*<       LOCD   = 1 >*/
    locd = 1;
/*<       LOCHY  = LOCD  + N >*/
    lochy = locd + *n;
/*<       LOCHZ  = LOCHY + M >*/
    lochz = lochy + *m;
/*<       LOCW   = LOCHZ + N >*/
    locw = lochz + *n;
/*<       LTOTAL = LOCW  + MAX(M,N) - 1 >*/
    ltotal = locw + max(*m,*n) - 1;
/*<       IF (LTOTAL .GT. LENRW) GO TO 900 >*/
    if (ltotal > 600) {
        goto L900;
    }
/*<    >*/
    lstp_(m, n, nduplc, npower, damp, xtrue, b, &rw[locd - 1], &rw[lochy - 1],
             &rw[lochz - 1], &rw[locw - 1], &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. */
/*<       CALL DCOPY ( M, B, 1, U, 1 ) >*/
    dcopy_(m, b, &c__1, u, &c__1);
/*<       ATOL   = 1.0E-10 >*/
    atol = (float)1e-10;
/*<       BTOL   = ATOL >*/
    btol = atol;
/*<       CONLIM = 10.0 * ACOND >*/
    conlim = acond * (float)10.;
/*<       ITNLIM = M + N + 50 >*/
    itnlim = *m + *n + 50;
/*<       NOUT   = 6 >*/
    nout = 6;
/*<    >*/
    printf("\n\n --------------------------------------------------------------------\n");
    printf(" Least-Squares Test Problem      P( %ld %ld %ld %ld %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, 0);
/*     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**2 >*/
/* Computing 2nd power */
    d__1 = *damp;
    dampsq = d__1 * d__1;
/*<       WRITE(NOUT, 2000) >*/
/*
 2000 FORMAT(
     $ // 22X, ' Residual norm    Residual norm    Solution norm'
     $  / 22X, '(Abar X - bbar)   (Normal eqns)         (X)' /)
*/
    printf("\n\n                       Residual norm    Residual norm    Solution norm\n");
    printf("                      (Abar X - bbar)   (Normal eqns)         (X)\n");
/*<       WRITE(NOUT, 2100) RNORM, ARNORM, XNORM >*/
/*
 2100 FORMAT(1P, ' Estimated by LSQR', 3E17.5)
*/
    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. */
/*<       CALL DCOPY ( M, B, 1, U, 1 ) >*/
    dcopy_(m, b, &c__1, u, &c__1);
/*<       CALL DSCAL ( M, (-ONE), U, 1 ) >*/
    dscal_(m, &c_b53, u, &c__1);
/*<       CALL APROD ( 1, M, N, X, U, LENIW, LENRW, IW, RW ) >*/
    aprod_(&c__1, m, n, x, u, &c__1, &c__600, iw, rw, 0);
/*     Compute  V = A(transpose)*U  +  DAMP**2 * X. */
/*     This will be close to zero in all cases */
/*     if  X  is close to a solution. */
/*<       CALL DCOPY ( N, X, 1, V, 1 ) >*/
    dcopy_(n, x, &c__1, v, &c__1);
/*<       CALL DSCAL ( N, DAMPSQ, V, 1 ) >*/
    dscal_(n, &dampsq, v, &c__1);
/*<       CALL APROD ( 2, M, N, V, U, LENIW, LENRW, IW, RW ) >*/
    aprod_(&c__2, m, n, v, u, &c__1, &c__600, iw, rw, 0);
/*     Compute the norms associated with  X, U, V. */
/*<       XNORM  = DNRM2 ( N, X, 1 ) >*/
    xnorm = dnrm2_(n, x, &c__1);
/*<       RNORM  = SQRT( DNRM2 ( M, U, 1 )**2  +  DAMPSQ * XNORM**2 ) >*/
/* Computing 2nd power */
    d__1 = dnrm2_(m, u, &c__1);
/* Computing 2nd power */
    d__2 = xnorm;
    rnorm = sqrt(d__1 * d__1 + dampsq * (d__2 * d__2));
/*<       ARNORM = DNRM2 ( N, V, 1 ) >*/
    arnorm = dnrm2_(n, v, &c__1);
/*<       WRITE(NOUT, 2200) RNORM, ARNORM, XNORM >*/
/*
 2200 FORMAT(1P, ' Computed from  X ', 3E17.5) 
*/
    printf(" Computed from  X %17.5e%17.5e%17.5e\n", rnorm, arnorm, xnorm);
/*     Print the solution and standard error estimates from  LSQR. */
/*<       WRITE(NOUT, 2500) (J, X(J),  J = 1, N) >*/
/*
 2500 FORMAT(//' Solution  X' / 4(I6, G14.6))
*/
    printf("\n\n Solution  X\n");
    for (j = 0; j < *n; ++j)
        printf("%6ld%14.6g\n", j+1, x[j]);
/*<       WRITE(NOUT, 2600) (J, SE(J), J = 1, N) >*/
/*
 2600 FORMAT(/ ' Standard errors  SE' / 4(I6, G14.6))
*/
    printf("\n\n Standard errors  SE\n");
    for (j = 0; j < *n; ++j)
        printf("%6ld%14.6g\n", j+1, se[j]);
    printf("\n");
/*     Print a clue about whether the solution looks OK. */
/*<       DO 500 J = 1, N >*/
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
/*<          W(J)  = X(J) - XTRUE(J) >*/
        w[j - 1] = x[j - 1] - xtrue[j - 1];
/*<   500 CONTINUE >*/
/* L500: */
    }
/*<       ENORM    = DNRM2 ( N, W, 1 ) / (ONE  +  DNRM2 ( N, XTRUE, 1 )) >*/
    enorm = dnrm2_(n, w, &c__1) / (dnrm2_(n, xtrue, &c__1) + 1.);
/*<       ETOL     = 1.0D-5 >*/
    etol = 1e-5;
/*<       IF (ENORM .LE. ETOL) WRITE(NOUT, 3000) ENORM >*/
/*
 3000 FORMAT(1P / ' LSQR  appears to be successful.',
     $        '     Relative error in  X  =', E10.2)
*/
    if (enorm <= etol) {
        printf("\n LSQR  appears to be successful.     Relative error in  X  =%10.2e\n", enorm);
    }
/*<       IF (ENORM .GT. ETOL) WRITE(NOUT, 3100) ENORM >*/
/*
 3100 FORMAT(1P / ' LSQR  appears to have failed.  ',
     $        '     Relative error in  X  =', E10.2)
*/
    if (enorm > etol) {
        printf("\n LSQR  appears to have failed.       Relative error in  X  =%10.2e\n", enorm);
    }
/*<       RETURN >*/
    return 0;
/*     Not enough workspace. */
/*<   900 WRITE(NOUT, 9000) LTOTAL >*/
/*
 9000 FORMAT(/ ' XXX  Insufficient workspace.',
     $        '  The length of  RW  should be at least', I6)
*/
L900:
    printf("\n XXX  Insufficient workspace."
           "  The length of  RW  should be at least %ld\n", ltotal);
/*<       RETURN >*/
    return 0;
/*<  1 >*/
/*<  2 >*/
/*<  2100 FORMAT(1P, ' Estimated by LSQR', 3E17.5) >*/
/*<  2200 FORMAT(1P, ' Computed from  X ', 3E17.5)  >*/
/*<  2500 FORMAT(//' Solution  X' / 4(I6, G14.6)) >*/
/*<  2600 FORMAT(/ ' Standard errors  SE' / 4(I6, G14.6)) >*/
/*<  3 >*/
/*<  3 >*/
/*<  9 >*/
/*     End of TEST */
/*<       END >*/
} /* test_ */

/*     ------------- */
/*     Main program. */
/*     ------------- */
/*<       DOUBLE PRECISION   DAMP1, DAMP2, DAMP3, DAMP4, ZERO >*/
/* Main program */ int main()
{
    /* Local variables */
    doublereal zero;
    extern /* Subroutine */ int test_(integer *, integer *, integer *, 
            integer *, doublereal *);
    doublereal damp1, damp2, damp3, damp4;


/*<       ZERO   = 0.0 >*/
    zero = (float)0.;
/*<       DAMP1  = 0.1 >*/
    damp1 = (float).1;
/*<       DAMP2  = 0.01 >*/
    damp2 = (float).01;
/*<       DAMP3  = 0.001 >*/
    damp3 = (float).001;
/*<       DAMP4  = 0.0001 >*/
    damp4 = (float)1e-4;
/*<       CALL TEST  (  1,  1, 1, 1, ZERO  ) >*/
    test_(&c__1, &c__1, &c__1, &c__1, &zero);
/*<       CALL TEST  (  2,  1, 1, 1, ZERO  ) >*/
    test_(&c__2, &c__1, &c__1, &c__1, &zero);
/*<       CALL TEST  ( 40, 40, 4, 4, ZERO  ) >*/
    test_(&c__40, &c__40, &c__4, &c__4, &zero);
/*<       CALL TEST  ( 40, 40, 4, 4, DAMP2 ) >*/
    test_(&c__40, &c__40, &c__4, &c__4, &damp2);
/*<       CALL TEST  ( 80, 40, 4, 4, DAMP2 ) >*/
    test_(&c__80, &c__40, &c__4, &c__4, &damp2);
/*<       STOP >*/
/*     End of main program for testing LSQR */
/*<       END >*/
    return 0;
} /* MAIN__ */