lsqr.c

来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 602 行 · 第 1/2 页

/*          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. */
/* ----------------------------------------------------------------------- */
/*     Initialize. */
/*      IF (NOUT .GT. 0) */
/*     $   WRITE(NOUT, 1000) ENTER, M, N, DAMP, ATOL, CONLIM, BTOL, ITNLIM */
    *itn = 0;
    *istop = 0;
    nstop = 0;
    ctol = 0.;
    if (*conlim > 0.) {
        ctol = 1. / *conlim;
    }
    *anorm = 0.;
    *acond = 0.;
    bbnorm = 0.;
    dampsq = *damp * *damp;
    ddnorm = 0.;
    res2 = 0.;
    *xnorm = 0.;
    xxnorm = 0.;
    cs2 = -1.;
    sn2 = 0.;
    z = 0.;
    for (i = 0; i < *n; ++i) {
        v[i] = 0.;
        x[i] = 0.;
        se[i] = 0.;
    }
/*     Set up the first vectors U and V for the bidiagonalization. */
/*     These satisfy  BETA*U = b,  ALFA*V = A(transpose)*U. */
    alfa = 0.;
    beta = dnrm2_(m, u, &c__1);
    if (beta > 0.) {
        d__1 = 1. / beta;
        dscal_(m, &d__1, u, &c__1);
        (*aprod)(&c__2, m, n, v, u, leniw, lenrw, iw, rw);
        alfa = dnrm2_(n, v, &c__1);
    }
    if (alfa > 0.) {
        d__1 = 1. / alfa;
        dscal_(n, &d__1, v, &c__1);
        dcopy_(n, v, &c__1, w, &c__1);
    }
    *arnorm = alfa * beta;
    if (*arnorm == 0.) {
        goto L800;
    }
    rhobar = alfa;
    phibar = beta;
    bnorm = 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. */
/*     ------------------------------------------------------------------ */
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. */
    d__1 = -alfa;
    dscal_(m, &d__1, u, &c__1);
    (*aprod)(&c__1, m, n, v, u, leniw, lenrw, iw, rw);
    beta = dnrm2_(m, u, &c__1);
    bbnorm = bbnorm + alfa * alfa + beta * beta + dampsq;
    if (beta > 0.) {
        d__1 = 1. / beta;
        dscal_(m, &d__1, u, &c__1);
        d__1 = -beta;
        dscal_(n, &d__1, v, &c__1);
        (*aprod)(&c__2, m, n, v, u, leniw, lenrw, iw, rw);
        alfa = dnrm2_(n, v, &c__1);
        if (alfa > 0.) {
            d__1 = 1. / alfa;
            dscal_(n, &d__1, v, &c__1);
        }
    }
/*     Use a plane rotation to eliminate the damping parameter. */
/*     This alters the diagonal (RHOBAR) of the lower-bidiagonal matrix.  */
    rhbar2 = rhobar * rhobar + dampsq;
    rhbar1 = sqrt(rhbar2);
    cs1 = rhobar / rhbar1;
    sn1 = *damp / rhbar1;
    psi = sn1 * phibar;
    phibar *= cs1;
/*     Use a plane rotation to eliminate the subdiagonal element (BETA) */
/*     of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.  */
    rho = sqrt(rhbar2 + beta * beta);
    cs = rhbar1 / rho;
    sn = beta / rho;
    theta = sn * alfa;
    rhobar = -cs * alfa;
    phi = cs * phibar;
    phibar *= sn;
    tau = sn * phi;
/*     Update  X, W  and the standard error estimates. */
    t1 = phi / rho;
    t2 = -theta / rho;
    t3 = 1. / rho;
    for (i = 0; i < *n; ++i) {
        t = w[i];
        x[i] += t1 * t;
        w[i] = t2 * t + v[i];
        t *= t3; t *= t;
        se[i] += t;
        ddnorm += t;
    }
/*     Use a plane rotation on the right to eliminate the */
/*     super-diagonal element (THETA) of the upper-bidiagonal matrix. */
/*     Then use the result to estimate  norm(X). */
    delta = sn2 * rho;
    gambar = -cs2 * rho;
    rhs = phi - delta * z;
    zbar = rhs / gambar;
    *xnorm = sqrt(xxnorm + zbar * zbar);
    gamma = sqrt(gambar * gambar + theta * theta);
    cs2 = gambar / gamma;
    sn2 = theta / gamma;
    z = rhs / gamma;
    xxnorm += z * z;
/*     Test for convergence. */
/*     First, estimate the norm and condition of the matrix  Abar, */
/*     and the norms of  rbar  and  Abar(transpose)*rbar. */
    *anorm = sqrt(bbnorm);
    *acond = *anorm * sqrt(ddnorm);
    res1 = phibar * phibar;
    res2 += psi * psi;
    *rnorm = sqrt(res1 + res2);
    *arnorm = alfa * abs(tau);
/*     Now use these norms to estimate certain other quantities, */
/*     some of which will be small near a solution. */
    test1 = *rnorm / bnorm;
    test2 = 0.;
    if (*rnorm > 0.) {
        test2 = *arnorm / (*anorm * *rnorm);
    }
    test3 = 1. / *acond;
    t1 = test1 / (*anorm * *xnorm / bnorm + 1.);
    rtol = *btol + *atol * *anorm * *xnorm / bnorm;
/*     The following tests guard against extremely small values of */
/*     ATOL, BTOL  or  CTOL.  (The user may have set any or all of */
/*     the parameters  ATOL, BTOL, CONLIM  to zero.) */
/*     The effect is equivalent to the normal tests using */
/*     ATOL = RELPR,  BTOL = RELPR,  CONLIM = 1/RELPR. */
    t3 = test3 + 1.;
    t2 = test2 + 1.;
    t1 += 1.;
    if (*itn >= *itnlim) {
        *istop = 7;
    }
    if (t3 <= 1.) {
        *istop = 6;
    }
    if (t2 <= 1.) {
        *istop = 5;
    }
    if (t1 <= 1.) {
        *istop = 4;
    }
/*     Allow for tolerances set by the user. */
    if (test3 <= ctol) {
        *istop = 3;
    }
    if (test2 <= *atol) {
        *istop = 2;
    }
    if (test1 <= rtol) {
        *istop = 1;
    }
/*     ================================================================== */
/*     See if it is time to print something. */
    if (*nout <= 0) {
        goto L600;
    }
    if (*n <= 40) {
        goto L400;
    }
    if (*itn <= 10) {
        goto L400;
    }
    if (*itn >= *itnlim - 10) {
        goto L400;
    }
    if (*itn % 10 == 0) {
        goto L400;
    }
    if (test3 <= ctol * 2.f) {
        goto L400;
    }
    if (test2 <= *atol * 10.f) {
        goto L400;
    }
    if (test1 <= rtol * 10.f) {
        goto L400;
    }
    if (*istop != 0) {
        goto L400;
    }
    goto L600;
/*     Print a line for this iteration. */
L400:
    if (TRUE_) {
        goto L600;
    }
/*  400 WRITE(NOUT, 1500) ITN, X(1), RNORM, TEST1, TEST2, ANORM, ACOND */
/*  IF (MOD(ITN,10) .EQ. 0) WRITE(NOUT, 1600) */
/*     ================================================================== */
/*     Stop if appropriate. */
/*     The convergence criteria are required to be met on  NCONV */
/*     consecutive iterations, where  NCONV  is set below. */
/*     Suggested value:  NCONV = 1, 2  or  3. */
L600:
    if (*istop == 0) {
        nstop = 0;
    }
    if (*istop == 0) {
        goto L100;
    }
    nconv = 1;
    ++nstop;
    if (nstop < nconv && *itn < *itnlim) {
        *istop = 0;
    }
    if (*istop == 0) {
        goto L100;
    }
/*     ------------------------------------------------------------------ */
/*     End of iteration loop. */
/*     ------------------------------------------------------------------ */
/*     Finish off the standard error estimates. */
    t = 1.;
    if (*m > *n) {
        t = (doublereal) (*m - *n);
    }
    if (dampsq > 0.) {
        t = (doublereal) (*m);
    }
    t = *rnorm / sqrt(t);
    for (i = 0; i < *n; ++i) {
        se[i] = t * sqrt(se[i]);
    }
/*     Print the stopping condition. */
L800:
    if (TRUE_) {
        goto L900;
    }
/*  800 IF (NOUT .GT. 0) THEN */
/*         WRITE(NOUT, 2000) EXIT, ISTOP, ITN, */
/*     $                     EXIT, ANORM, ACOND, */
/*     $                     EXIT, RNORM, ARNORM, */
/*     $                     EXIT, BNORM, XNORM */
/*         WRITE(NOUT, 3000) EXIT, MSG(ISTOP) */
/*      END IF */
L900:
    return;
/*     ------------------------------------------------------------------ */
/*     ------------------------------------------------------------------ */
} /* lsqr_ */