📄 lsqr.f
字号:
* 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 ALFA, BBNORM, BETA, BNORM,
$ CS, CS1, CS2, CTOL, DAMPSQ, DDNORM, DELTA,
$ GAMMA, GAMBAR, PHI, PHIBAR, PSI,
$ RES1, RES2, RHO, RHOBAR, RHBAR1, RHBAR2,
$ RHS, RTOL, SN, SN1, SN2,
$ T, TAU, TEST1, TEST2, TEST3,
$ THETA, T1, T2, T3, XXNORM, Z, ZBAR
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
C CHARACTER*16 ENTER, EXIT
C CHARACTER*60 MSG(0:7)
C DATA ENTER /' Enter LSQR. '/,
C $ EXIT /' Exit LSQR. '/
C DATA MSG
C $ / 'The exact solution is X = 0',
C $ 'Ax - b is small enough, given ATOL, BTOL',
C $ 'The least-squares solution is good enough, given ATOL',
C $ 'The estimate of cond(Abar) has exceeded CONLIM',
C $ 'Ax - b is small enough for this machine',
C $ 'The least-squares solution is good enough for this machine',
C $ 'Cond(Abar) seems to be too large for this machine',
C $ 'The iteration limit has been reached' /
*-----------------------------------------------------------------------
* Initialize.
C IF (NOUT .GT. 0)
C $ WRITE(NOUT, 1000) ENTER, M, N, DAMP, ATOL, CONLIM, BTOL, ITNLIM
ITN = 0
ISTOP = 0
NSTOP = 0
CTOL = ZERO
IF (CONLIM .GT. ZERO) CTOL = ONE / CONLIM
ANORM = ZERO
ACOND = ZERO
BBNORM = ZERO
DAMPSQ = DAMP**2
DDNORM = ZERO
RES2 = ZERO
XNORM = ZERO
XXNORM = ZERO
CS2 = - ONE
SN2 = ZERO
Z = ZERO
DO 10 I = 1, N
V(I) = ZERO
X(I) = ZERO
SE(I) = ZERO
10 CONTINUE
* Set up the first vectors U and V for the bidiagonalization.
* These satisfy BETA*U = b, ALFA*V = A(transpose)*U.
ALFA = ZERO
BETA = DNRM2 ( M, U, 1 )
IF (BETA .GT. ZERO) THEN
CALL DSCAL ( M, (ONE / BETA), U, 1 )
CALL APROD ( 2, M, N, V, U, LENIW, LENRW, IW, RW )
ALFA = DNRM2 ( N, V, 1 )
END IF
IF (ALFA .GT. ZERO) THEN
CALL DSCAL ( N, (ONE / ALFA), V, 1 )
CALL DCOPY ( N, V, 1, W, 1 )
END IF
ARNORM = ALFA * BETA
IF (ARNORM .EQ. ZERO) GO TO 800
RHOBAR = ALFA
PHIBAR = BETA
BNORM = BETA
RNORM = BETA
C IF (NOUT .GT. 0 ) THEN
C IF (DAMPSQ .EQ. ZERO) THEN
C WRITE(NOUT, 1200)
C ELSE
C WRITE(NOUT, 1300)
C END IF
C TEST1 = ONE
C TEST2 = ALFA / BETA
C WRITE(NOUT, 1500) ITN, X(1), RNORM, TEST1, TEST2
C WRITE(NOUT, 1600)
C END IF
* ------------------------------------------------------------------
* Main iteration loop.
* ------------------------------------------------------------------
100 ITN = ITN + 1
* 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 )
CALL APROD ( 1, M, N, V, U, LENIW, LENRW, IW, RW )
BETA = DNRM2 ( M, U, 1 )
BBNORM = BBNORM + ALFA**2 + BETA**2 + DAMPSQ
IF (BETA .GT. ZERO) THEN
CALL DSCAL ( M, (ONE / BETA), U, 1 )
CALL DSCAL ( N, (- BETA), V, 1 )
CALL APROD ( 2, M, N, V, U, LENIW, LENRW, IW, RW )
ALFA = DNRM2 ( N, V, 1 )
IF (ALFA .GT. ZERO) THEN
CALL DSCAL ( N, (ONE / ALFA), V, 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
RHBAR1 = SQRT( RHBAR2 )
CS1 = RHOBAR / RHBAR1
SN1 = DAMP / RHBAR1
PSI = SN1 * 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 )
CS = RHBAR1 / RHO
SN = BETA / RHO
THETA = SN * ALFA
RHOBAR = - CS * ALFA
PHI = CS * PHIBAR
PHIBAR = SN * PHIBAR
TAU = SN * PHI
* Update X, W and the standard error estimates.
T1 = PHI / RHO
T2 = - THETA / RHO
T3 = ONE / RHO
DO 200 I = 1, N
T = W(I)
X(I) = T1*T + X(I)
W(I) = T2*T + V(I)
T = (T3*T)**2
SE(I) = T + SE(I)
DDNORM = T + DDNORM
200 CONTINUE
* 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 **2 )
GAMMA = SQRT( GAMBAR**2 + THETA**2 )
CS2 = GAMBAR / GAMMA
SN2 = THETA / GAMMA
Z = RHS / GAMMA
XXNORM = XXNORM + Z**2
* 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**2
RES2 = RES2 + PSI**2
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 = ZERO
IF (RNORM .GT. ZERO) TEST2 = ARNORM / (ANORM * RNORM)
TEST3 = ONE / ACOND
T1 = TEST1 / (ONE + 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
T2 = ONE + TEST2
T1 = ONE + T1
IF (ITN .GE. ITNLIM) ISTOP = 7
IF (T3 .LE. ONE ) ISTOP = 6
IF (T2 .LE. ONE ) ISTOP = 5
IF (T1 .LE. ONE ) ISTOP = 4
* Allow for tolerances set by the user.
IF (TEST3 .LE. CTOL) ISTOP = 3
IF (TEST2 .LE. ATOL) ISTOP = 2
IF (TEST1 .LE. RTOL) ISTOP = 1
* ==================================================================
* See if it is time to print something.
IF (NOUT .LE. 0 ) GO TO 600
IF (N .LE. 40 ) GO TO 400
IF (ITN .LE. 10 ) GO TO 400
IF (ITN .GE. ITNLIM-10) GO TO 400
IF (MOD(ITN,10) .EQ. 0 ) GO TO 400
IF (TEST3 .LE. 2.0*CTOL) GO TO 400
IF (TEST2 .LE. 10.0*ATOL) GO TO 400
IF (TEST1 .LE. 10.0*RTOL) GO TO 400
IF (ISTOP .NE. 0 ) GO TO 400
GO TO 600
* Print a line for this iteration.
400 IF (1 .EQ. 1) GO TO 600
C 400 WRITE(NOUT, 1500) ITN, X(1), RNORM, TEST1, TEST2, ANORM, ACOND
C 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
IF (ISTOP .EQ. 0) GO TO 100
NCONV = 1
NSTOP = NSTOP + 1
IF (NSTOP .LT. NCONV .AND. ITN .LT. ITNLIM) ISTOP = 0
IF (ISTOP .EQ. 0) GO TO 100
* ------------------------------------------------------------------
* End of iteration loop.
* ------------------------------------------------------------------
* Finish off the standard error estimates.
T = ONE
IF (M .GT. N ) T = M - N
IF (DAMPSQ .GT. ZERO) T = M
T = RNORM / SQRT( T )
DO 700 I = 1, N
SE(I) = T * SQRT( SE(I) )
700 CONTINUE
* Print the stopping condition.
800 IF (1 .EQ. 1) GO TO 900
C 800 IF (NOUT .GT. 0) THEN
C WRITE(NOUT, 2000) EXIT, ISTOP, ITN,
C $ EXIT, ANORM, ACOND,
C $ EXIT, RNORM, ARNORM,
C $ EXIT, BNORM, XNORM
C WRITE(NOUT, 3000) EXIT, MSG(ISTOP)
C END IF
900 RETURN
* ------------------------------------------------------------------
1000 FORMAT(// 1P, A, ' Least-squares solution of A*x = b'
$ / ' The matrix A has', I7, ' rows and', I7, ' columns'
$ / ' The damping parameter is DAMP =', E10.2
$ / ' ATOL =', E10.2, 15X, 'CONLIM =', E10.2
$ / ' BTOL =', E10.2, 15X, 'ITNLIM =', I10)
1200 FORMAT(// ' Itn x(1) Function',
$ ' Compatible LS Norm A Cond A' /)
1300 FORMAT(// ' Itn x(1) Function',
$ ' Compatible LS Norm Abar Cond Abar' /)
1500 FORMAT(1P, I6, 2E17.9, 4E10.2)
1600 FORMAT(1X)
2000 FORMAT(/ 1P, A, 6X, 'ISTOP =', I3, 16X, 'ITN =', I9
$ / A, 6X, 'ANORM =', E13.5, 6X, 'ACOND =', E13.5
$ / A, 6X, 'RNORM =', E13.5, 6X, 'ARNORM =', E13.5,
$ / A, 6X, 'BNORM =', E13.5, 6X, 'XNORM =', E13.5)
3000 FORMAT( A, 6X, A )
* ------------------------------------------------------------------
* End of LSQR
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -