lsqr.c

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

    cs = rhbar1 / rho;
/*<       SN     =   BETA   / RHO >*/
    sn = beta / rho;
/*<       THETA  =   SN * ALFA >*/
    theta = sn * alfa;
/*<       RHOBAR = - CS * ALFA >*/
    rhobar = -cs * alfa;
/*<       PHI    =   CS * PHIBAR >*/
    phi = cs * phibar;
/*<       PHIBAR =   SN * PHIBAR >*/
    phibar = sn * phibar;
/*<       TAU    =   SN * PHI >*/
    tau = sn * phi;
/*     Update  X, W  and the standard error estimates. */
/*<       T1     =   PHI   / RHO >*/
    t1 = phi / rho;
/*<       T2     = - THETA / RHO >*/
    t2 = -theta / rho;
/*<       T3     =   ONE   / RHO >*/
    t3 = 1. / rho;
/*<       DO 200  I =  1, N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          T      =  W(I) >*/
        t = w[i__];
/*<          X(I)   =  T1*T  +  X(I) >*/
        x[i__] = t1 * t + x[i__];
/*<          W(I)   =  T2*T  +  V(I) >*/
        w[i__] = t2 * t + v[i__];
/*<          T      = (T3*T)**2 >*/
/* Computing 2nd power */
        d__1 = t3 * t;
        t = d__1 * d__1;
/*<          SE(I)  =  T     +  SE(I) >*/
        se[i__] = t + se[i__];
/*<          DDNORM =  T     +  DDNORM >*/
        ddnorm = t + ddnorm;
/*<   200 CONTINUE >*/
/* L200: */
    }
/*     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 >*/
    delta = sn2 * rho;
/*<       GAMBAR = - CS2 * RHO >*/
    gambar = -cs2 * rho;
/*<       RHS    =   PHI    - DELTA * Z >*/
    rhs = phi - delta * z__;
/*<       ZBAR   =   RHS    / GAMBAR >*/
    zbar = rhs / gambar;
/*<       XNORM  =   SQRT( XXNORM    + ZBAR **2 ) >*/
/* Computing 2nd power */
    d__1 = zbar;
    *xnorm = sqrt(xxnorm + d__1 * d__1);
/*<       GAMMA  =   SQRT( GAMBAR**2 + THETA**2 ) >*/
/* Computing 2nd power */
    d__1 = gambar;
/* Computing 2nd power */
    d__2 = theta;
    gamma = sqrt(d__1 * d__1 + d__2 * d__2);
/*<       CS2    =   GAMBAR / GAMMA >*/
    cs2 = gambar / gamma;
/*<       SN2    =   THETA  / GAMMA >*/
    sn2 = theta / gamma;
/*<       Z      =   RHS    / GAMMA >*/
    z__ = rhs / gamma;
/*<       XXNORM =   XXNORM + Z**2 >*/
/* Computing 2nd power */
    d__1 = z__;
    xxnorm += d__1 * d__1;
/*     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 ) >*/
    *anorm = sqrt(bbnorm);
/*<       ACOND  =   ANORM * SQRT( DDNORM ) >*/
    *acond = *anorm * sqrt(ddnorm);
/*<       RES1   =   PHIBAR**2 >*/
/* Computing 2nd power */
    d__1 = phibar;
    res1 = d__1 * d__1;
/*<       RES2   =   RES2  +  PSI**2 >*/
/* Computing 2nd power */
    d__1 = psi;
    res2 += d__1 * d__1;
/*<       RNORM  =   SQRT( RES1 + RES2 ) >*/
    *rnorm = sqrt(res1 + res2);
/*<       ARNORM =   ALFA  * ABS( TAU ) >*/
    *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 >*/
    test1 = *rnorm / bnorm;
/*<       TEST2  =   ZERO >*/
    test2 = 0.;
/*<       IF (RNORM .GT. ZERO) TEST2 = ARNORM / (ANORM * RNORM) >*/
    if (*rnorm > 0.) {
        test2 = *arnorm / (*anorm * *rnorm);
    }
/*<       TEST3  =   ONE   /  ACOND >*/
    test3 = 1. / *acond;
/*<       T1     =   TEST1 / (ONE  +  ANORM * XNORM / BNORM) >*/
    t1 = test1 / (*anorm * *xnorm / bnorm + 1.);
/*<       RTOL   =   BTOL  +  ATOL *  ANORM * XNORM / BNORM >*/
    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     =   ONE + TEST3 >*/
    t3 = test3 + 1.;
/*<       T2     =   ONE + TEST2 >*/
    t2 = test2 + 1.;
/*<       T1     =   ONE + T1 >*/
    t1 += 1.;
/*<       IF (ITN .GE. ITNLIM) ISTOP = 7 >*/
    if (*itn >= *itnlim) {
        *istop = 7;
    }
/*<       IF (T3  .LE. ONE   ) ISTOP = 6 >*/
    if (t3 <= 1.) {
        *istop = 6;
    }
/*<       IF (T2  .LE. ONE   ) ISTOP = 5 >*/
    if (t2 <= 1.) {
        *istop = 5;
    }
/*<       IF (T1  .LE. ONE   ) ISTOP = 4 >*/
    if (t1 <= 1.) {
        *istop = 4;
    }
/*     Allow for tolerances set by the user. */
/*<       IF (TEST3 .LE. CTOL) ISTOP = 3 >*/
    if (test3 <= ctol) {
        *istop = 3;
    }
/*<       IF (TEST2 .LE. ATOL) ISTOP = 2 >*/
    if (test2 <= *atol) {
        *istop = 2;
    }
/*<       IF (TEST1 .LE. RTOL) ISTOP = 1 >*/
    if (test1 <= rtol) {
        *istop = 1;
    }
/*     ================================================================== */
/*     See if it is time to print something. */
/*<       IF (NOUT  .LE.  0       ) GO TO 600 >*/
    if (*nout <= 0) {
        goto L600;
    }
/*<       IF (N     .LE. 40       ) GO TO 400 >*/
    if (*n <= 40) {
        goto L400;
    }
/*<       IF (ITN   .LE. 10       ) GO TO 400 >*/
    if (*itn <= 10) {
        goto L400;
    }
/*<       IF (ITN   .GE. ITNLIM-10) GO TO 400 >*/
    if (*itn >= *itnlim - 10) {
        goto L400;
    }
/*<       IF (MOD(ITN,10) .EQ. 0  ) GO TO 400 >*/
    if (*itn % 10 == 0) {
        goto L400;
    }
/*<       IF (TEST3 .LE.  2.0*CTOL) GO TO 400 >*/
    if (test3 <= ctol * (float)2.) {
        goto L400;
    }
/*<       IF (TEST2 .LE. 10.0*ATOL) GO TO 400 >*/
    if (test2 <= *atol * (float)10.) {
        goto L400;
    }
/*<       IF (TEST1 .LE. 10.0*RTOL) GO TO 400 >*/
    if (test1 <= rtol * (float)10.) {
        goto L400;
    }
/*<       IF (ISTOP .NE.  0       ) GO TO 400 >*/
    if (*istop != 0) {
        goto L400;
    }
/*<       GO TO 600 >*/
    goto L600;
/*     Print a line for this iteration. */
/*<   400 IF (1 .EQ. 1) GO TO 600 >*/
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. */
/*<   600 IF (ISTOP .EQ. 0) NSTOP = 0 >*/
L600:
    if (*istop == 0) {
        nstop = 0;
    }
/*<       IF (ISTOP .EQ. 0) GO TO 100 >*/
    if (*istop == 0) {
        goto L100;
    }
/*<       NCONV  =   1 >*/
    nconv = 1;
/*<       NSTOP  =   NSTOP + 1 >*/
    ++nstop;
/*<       IF (NSTOP .LT. NCONV  .AND.  ITN .LT. ITNLIM) ISTOP = 0 >*/
    if (nstop < nconv && *itn < *itnlim) {
        *istop = 0;
    }
/*<       IF (ISTOP .EQ. 0) GO TO 100 >*/
    if (*istop == 0) {
        goto L100;
    }
/*     ------------------------------------------------------------------ */
/*     End of iteration loop. */
/*     ------------------------------------------------------------------ */
/*     Finish off the standard error estimates. */
/*<       T    =   ONE >*/
    t = 1.;
/*<       IF (M      .GT.   N )  T = M - N >*/
    if (*m > *n) {
        t = (doublereal) (*m - *n);
    }
/*<       IF (DAMPSQ .GT. ZERO)  T = M >*/
    if (dampsq > 0.) {
        t = (doublereal) (*m);
    }
/*<       T    =   RNORM / SQRT( T ) >*/
    t = *rnorm / sqrt(t);
/*<       DO 700  I = 1, N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          SE(I)  = T * SQRT( SE(I) ) >*/
        se[i__] = t * sqrt(se[i__]);
/*<   700 CONTINUE >*/
/* L700: */
    }
/*     Print the stopping condition. */
/*<   800 IF (1 .EQ. 1) GO TO 900 >*/
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 */
/*<   900 RETURN >*/
L900:
    return 0;
/*     ------------------------------------------------------------------ */
/*<  1 >*/
/* L1000: */
/*<  1 >*/
/* L1200: */
/*<  1 >*/
/* L1300: */
/*<  1500 FORMAT(1P, I6, 2E17.9, 4E10.2) >*/
/* L1500: */
/*<  1600 FORMAT(1X) >*/
/* L1600: */
/*<  2 >*/
/* L2000: */
/*<  3000 FORMAT( A, 6X, A ) >*/
/* L3000: */
/*     ------------------------------------------------------------------ */
/*     End of LSQR */
/*<       END >*/
} /* lsqr_ */

#ifdef __cplusplus
        }
#endif