lsqr.c

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

/*                        ANORM should increase to roughly sqrt(n). */
/*                        A radically different value for ANORM may */
/*                        indicate an error in subroutine APROD (there */
/*                        may be an inconsistency between modes 1 and 2). */

/*     ACOND   output     An estimate of cond(Abar), the condition */
/*                        number of Abar.  A very high value of ACOND */
/*                        may again indicate an error in APROD. */

/*     RNORM   output     An estimate of the final value of norm(rbar), */
/*                        the function being minimized (see notation */
/*                        above).  This will be small if A*x = b has */
/*                        a solution. */

/*     ARNORM  output     An estimate of the final value of */
/*                        norm( Abar(transpose)*rbar ), the norm of */
/*                        the residual for the usual normal equations. */
/*                        This should be small in all cases.  (ARNORM */
/*                        will often be smaller than the true value */
/*                        computed from the output vector X.) */

/*     XNORM   output     An estimate of the norm of the final */
/*                        solution vector X. */


/*     Subroutines and functions used */
/*     ------------------------------ */

/*     USER               APROD */
/*     BLAS               DCOPY, DNRM2, DSCAL (see Lawson et al. below) */


/*     Precision */
/*     --------- */

/*     The number of iterations required by LSQR will usually decrease */
/*     if the computation is performed in higher precision.  To convert */
/*     LSQR between single and double precision, change the words */
/*                        DOUBLE PRECISION */
/*                        DCOPY, DNRM2, DSCAL */
/*     to the appropriate FORTRAN and BLAS equivalents. */
/*     Also change 'D+' or 'E+' in the PARAMETER statement. */


/*     References */
/*     ---------- */

/*     C.C. Paige and M.A. Saunders,  LSQR: An algorithm for sparse */
/*          linear equations and sparse least squares, */
/*          ACM Transactions on Mathematical Software 8, 1 (March 1982), */
/*          pp. 43-71. */

/*     C.C. Paige and M.A. Saunders,  Algorithm 583, LSQR: Sparse */
/*          linear equations and least-squares problems, */
/*          ACM Transactions on Mathematical Software 8, 2 (June 1982), */
/*          pp. 195-209. */

/*     C.L. Lawson, R.J. Hanson, D.R. Kincaid and F.T. Krogh, */
/*          Basic linear algebra subprograms for Fortran usage, */
/*          ACM Transactions on Mathematical Software 5, 3 (Sept 1979), */
/*          pp. 308-323 and 324-325. */
/* ----------------------------------------------------------------------- */


/*     LSQR development: */
/*     22 Feb 1982: LSQR sent to ACM TOMS to become Algorithm 583. */
/*     15 Sep 1985: Final F66 version.  LSQR sent to "misc" in netlib. */
/*     13 Oct 1987: Bug (Robert Davies, DSIR).  Have to delete */
/*                     IF ( (ONE + DABS(T)) .LE. ONE ) GO TO 200 */
/*                  from loop 200.  The test was an attempt to reduce */
/*                  underflows, but caused W(I) not to be updated. */
/*     17 Mar 1989: First F77 version. */
/*     04 May 1989: Bug (David Gay, AT&T).  When the second BETA is zero, */
/*                  RNORM = 0 and */
/*                  TEST2 = ARNORM / (ANORM * RNORM) overflows. */
/*                  Fixed by testing for RNORM = 0. */
/*     05 May 1989: Sent to "misc" in netlib. */

/*     Michael A. Saunders            (na.saunders @ NA-net.stanford.edu) */
/*     Department of Operations Research */
/*     Stanford University */
/*     Stanford, CA 94305-4022. */
/* ----------------------------------------------------------------------- */
/*     Intrinsics and local variables */
/*<       INTRINSIC          ABS, MOD, SQRT >*/
/*<       INTEGER            I, NCONV, NSTOP >*/
/*<       DOUBLE PRECISION   DNRM2 >*/
/*<    >*/
/*<       DOUBLE PRECISION   ZERO,           ONE >*/
/*<       PARAMETER        ( ZERO = 0.0D+0,  ONE = 1.0D+0 ) >*/
/*      CHARACTER*16       ENTER, EXIT */
/*      CHARACTER*60       MSG(0:7) */
/*      DATA               ENTER /' Enter LSQR.    '/, */
/*     $                   EXIT  /' Exit  LSQR.    '/ */
/*      DATA               MSG */
/*     $ / 'The exact solution is  X = 0', */
/*     $   'Ax - b is small enough, given ATOL, BTOL', */
/*     $   'The least-squares solution is good enough, given ATOL', */
/*     $   'The estimate of cond(Abar) has exceeded CONLIM', */
/*     $   'Ax - b is small enough for this machine', */
/*     $   'The least-squares solution is good enough for this machine', */
/*     $   'Cond(Abar) seems to be too large for this machine', */
/*     $   'The iteration limit has been reached' / */
/* ----------------------------------------------------------------------- */
/*     Initialize. */
/*      IF (NOUT .GT. 0) */
/*     $   WRITE(NOUT, 1000) ENTER, M, N, DAMP, ATOL, CONLIM, BTOL, ITNLIM */
/*<       ITN    =   0 >*/
    /* Parameter adjustments */
    --u;
    --se;
    --x;
    --w;
    --v;
    --iw;
    --rw;

    /* Function Body */
    *itn = 0;
/*<       ISTOP  =   0 >*/
    *istop = 0;
/*<       NSTOP  =   0 >*/
    nstop = 0;
/*<       CTOL   =   ZERO >*/
    ctol = 0.;
/*<       IF (CONLIM .GT. ZERO) CTOL = ONE / CONLIM >*/
    if (*conlim > 0.) {
        ctol = 1. / *conlim;
    }
/*<       ANORM  =   ZERO >*/
    *anorm = 0.;
/*<       ACOND  =   ZERO >*/
    *acond = 0.;
/*<       BBNORM =   ZERO >*/
    bbnorm = 0.;
/*<       DAMPSQ =   DAMP**2 >*/
/* Computing 2nd power */
    d__1 = *damp;
    dampsq = d__1 * d__1;
/*<       DDNORM =   ZERO >*/
    ddnorm = 0.;
/*<       RES2   =   ZERO >*/
    res2 = 0.;
/*<       XNORM  =   ZERO >*/
    *xnorm = 0.;
/*<       XXNORM =   ZERO >*/
    xxnorm = 0.;
/*<       CS2    = - ONE >*/
    cs2 = -1.;
/*<       SN2    =   ZERO >*/
    sn2 = 0.;
/*<       Z      =   ZERO >*/
    z__ = 0.;
/*<       DO 10  I = 1, N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          V(I)  =  ZERO >*/
        v[i__] = 0.;
/*<          X(I)  =  ZERO >*/
        x[i__] = 0.;
/*<         SE(I)  =  ZERO >*/
        se[i__] = 0.;
/*<    10 CONTINUE >*/
/* L10: */
    }
/*     Set up the first vectors U and V for the bidiagonalization. */
/*     These satisfy  BETA*U = b,  ALFA*V = A(transpose)*U. */
/*<       ALFA   =   ZERO >*/
    alfa = 0.;
/*<       BETA   =   DNRM2 ( M, U, 1 ) >*/
    beta = dnrm2_(m, &u[1], &c__1);
/*<       IF (BETA .GT. ZERO) THEN >*/
    if (beta > 0.) {
/*<          CALL DSCAL ( M, (ONE / BETA), U, 1 ) >*/
        d__1 = 1. / beta;
        dscal_(m, &d__1, &u[1], &c__1);
/*<          CALL APROD ( 2, M, N, V, U, LENIW, LENRW, IW, RW ) >*/
        (*aprod)(&c__2, m, n, &v[1], &u[1], leniw, lenrw, &iw[1], &rw[1],
                 userdata);
/*<          ALFA   =   DNRM2 ( N, V, 1 ) >*/
        alfa = dnrm2_(n, &v[1], &c__1);
/*<       END IF >*/
    }
/*<       IF (ALFA .GT. ZERO) THEN >*/
    if (alfa > 0.) {
/*<          CALL DSCAL ( N, (ONE / ALFA), V, 1 ) >*/
        d__1 = 1. / alfa;
        dscal_(n, &d__1, &v[1], &c__1);
/*<          CALL DCOPY ( N, V, 1, W, 1 ) >*/
        dcopy_(n, &v[1], &c__1, &w[1], &c__1);
/*<       END IF >*/
    }
/*<       ARNORM =   ALFA * BETA >*/
    *arnorm = alfa * beta;
/*<       IF (ARNORM .EQ. ZERO) GO TO 800 >*/
    if (*arnorm == 0.) {
        goto L800;
    }
/*<       RHOBAR =   ALFA >*/
    rhobar = alfa;
/*<       PHIBAR =   BETA >*/
    phibar = beta;
/*<       BNORM  =   BETA >*/
    bnorm = beta;
/*<       RNORM  =   BETA >*/
    *rnorm = beta;
/*      IF (NOUT   .GT.  0  ) THEN */
/*         IF (DAMPSQ .EQ. ZERO) THEN */
/*             WRITE(NOUT, 1200) */
/*         ELSE */
/*             WRITE(NOUT, 1300) */
/*         END IF */
/*         TEST1  = ONE */
/*         TEST2  = ALFA / BETA */
/*         WRITE(NOUT, 1500) ITN, X(1), RNORM, TEST1, TEST2 */
/*         WRITE(NOUT, 1600) */
/*      END IF */
/*     ------------------------------------------------------------------ */
/*     Main iteration loop. */
/*     ------------------------------------------------------------------ */
/*<   100 ITN    = ITN + 1 >*/
L100:
    ++(*itn);
/*     Perform the next step of the bidiagonalization to obtain the */
/*     next  BETA, U, ALFA, V.  These satisfy the relations */
/*                BETA*U  =  A*V  -  ALFA*U, */
/*                ALFA*V  =  A(transpose)*U  -  BETA*V. */
/*<       CALL DSCAL ( M, (- ALFA), U, 1 ) >*/
    d__1 = -alfa;
    dscal_(m, &d__1, &u[1], &c__1);
/*<       CALL APROD ( 1, M, N, V, U, LENIW, LENRW, IW, RW ) >*/
    (*aprod)(&c__1, m, n, &v[1], &u[1], leniw, lenrw, &iw[1], &rw[1],
             userdata);
/*<       BETA   =   DNRM2 ( M, U, 1 ) >*/
    beta = dnrm2_(m, &u[1], &c__1);
/*<       BBNORM =   BBNORM  +  ALFA**2  +  BETA**2  +  DAMPSQ >*/
/* Computing 2nd power */
    d__1 = alfa;
/* Computing 2nd power */
    d__2 = beta;
    bbnorm = bbnorm + d__1 * d__1 + d__2 * d__2 + dampsq;
/*<       IF (BETA .GT. ZERO) THEN >*/
    if (beta > 0.) {
/*<          CALL DSCAL ( M, (ONE / BETA), U, 1 ) >*/
        d__1 = 1. / beta;
        dscal_(m, &d__1, &u[1], &c__1);
/*<          CALL DSCAL ( N, (- BETA), V, 1 ) >*/
        d__1 = -beta;
        dscal_(n, &d__1, &v[1], &c__1);
/*<          CALL APROD ( 2, M, N, V, U, LENIW, LENRW, IW, RW ) >*/
        (*aprod)(&c__2, m, n, &v[1], &u[1], leniw, lenrw, &iw[1], &rw[1],
                 userdata);
/*<          ALFA   =   DNRM2 ( N, V, 1 ) >*/
        alfa = dnrm2_(n, &v[1], &c__1);
/*<          IF (ALFA .GT. ZERO) THEN >*/
        if (alfa > 0.) {
/*<             CALL DSCAL ( N, (ONE / ALFA), V, 1 ) >*/
            d__1 = 1. / alfa;
            dscal_(n, &d__1, &v[1], &c__1);
/*<          END IF >*/
        }
/*<       END IF >*/
    }
/*     Use a plane rotation to eliminate the damping parameter. */
/*     This alters the diagonal (RHOBAR) of the lower-bidiagonal matrix. */
/*<       RHBAR2 = RHOBAR**2  +  DAMPSQ >*/
/* Computing 2nd power */
    d__1 = rhobar;
    rhbar2 = d__1 * d__1 + dampsq;
/*<       RHBAR1 = SQRT( RHBAR2 ) >*/
    rhbar1 = sqrt(rhbar2);
/*<       CS1    = RHOBAR / RHBAR1 >*/
    cs1 = rhobar / rhbar1;
/*<       SN1    = DAMP   / RHBAR1 >*/
    sn1 = *damp / rhbar1;
/*<       PSI    = SN1 * PHIBAR >*/
    psi = sn1 * phibar;
/*<       PHIBAR = CS1 * PHIBAR >*/
    phibar = cs1 * phibar;
/*     Use a plane rotation to eliminate the subdiagonal element (BETA) */
/*     of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix. */
/*<       RHO    =   SQRT( RHBAR2  +  BETA**2 ) >*/
/* Computing 2nd power */
    d__1 = beta;
    rho = sqrt(rhbar2 + d__1 * d__1);
/*<       CS     =   RHBAR1 / RHO >*/