📄 lsqr-test.f
字号:
*********************************************************
*
* These routines are for testing LSQR.
*
*********************************************************
SUBROUTINE APROD ( MODE, M, N, X, Y, LENIW, LENRW, IW, RW )
INTEGER MODE, M, N, LENIW, LENRW
INTEGER IW(LENIW)
DOUBLE PRECISION X(N), Y(M), RW(LENRW)
* ------------------------------------------------------------------
* This is the matrix-vector product routine required by LSQR
* for a test matrix of the form A = HY*D*HZ. The quantities
* defining D, HY, HZ are in the work array RW, followed by a
* work array W. These are passed to APROD1 and APROD2 in order to
* make the code readable.
* ------------------------------------------------------------------
INTEGER LOCD, LOCHY, LOCHZ, LOCW
LOCD = 1
LOCHY = LOCD + N
LOCHZ = LOCHY + M
LOCW = LOCHZ + N
IF (MODE .EQ. 1) CALL APROD1( M, N, X, Y,
$ RW(LOCD), RW(LOCHY), RW(LOCHZ), RW(LOCW) )
IF (MODE .NE. 1) CALL APROD2( M, N, X, Y,
$ RW(LOCD), RW(LOCHY), RW(LOCHZ), RW(LOCW) )
* End of APROD
END
SUBROUTINE APROD1( M, N, X, Y, D, HY, HZ, W )
INTEGER M, N
DOUBLE PRECISION X(N), Y(M), D(N), HY(M), HZ(N), W(M)
* ------------------------------------------------------------------
* APROD1 computes Y = Y + A*X for subroutine APROD,
* where A is a test matrix of the form A = HY*D*HZ,
* and the latter matrices HY, D, HZ are represented by
* input vectors with the same name.
* ------------------------------------------------------------------
INTEGER I
DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0 )
CALL HPROD ( N, HZ, X, W )
DO 100 I = 1, N
W(I) = D(I) * W(I)
100 CONTINUE
DO 200 I = N + 1, M
W(I) = ZERO
200 CONTINUE
CALL HPROD ( M, HY, W, W )
DO 600 I = 1, M
Y(I) = Y(I) + W(I)
600 CONTINUE
* End of APROD1
END
SUBROUTINE APROD2( M, N, X, Y, D, HY, HZ, W )
INTEGER M, N
DOUBLE PRECISION X(N), Y(M), D(N), HY(M), HZ(N), W(M)
* ------------------------------------------------------------------
* APROD2 computes X = X + A(T)*Y for subroutine APROD,
* where A is a test matrix of the form A = HY*D*HZ,
* and the latter matrices HY, D, HZ are represented by
* input vectors with the same name.
* ------------------------------------------------------------------
INTEGER I
CALL HPROD ( M, HY, Y, W )
DO 100 I = 1, N
W(I) = D(I)*W(I)
100 CONTINUE
CALL HPROD ( N, HZ, W, W )
DO 600 I = 1, N
X(I) = X(I) + W(I)
600 CONTINUE
* End of APROD2
END
SUBROUTINE HPROD ( N, HZ, X, Y )
INTEGER N
DOUBLE PRECISION HZ(N), X(N), Y(N)
* ------------------------------------------------------------------
* HPROD applies a Householder transformation stored in HZ
* to get Y = ( I - 2*HZ*HZ(transpose) ) * X.
* ------------------------------------------------------------------
INTEGER I
DOUBLE PRECISION S
S = 0.0
DO 100 I = 1, N
S = HZ(I) * X(I) + S
100 CONTINUE
S = S + S
DO 200 I = 1, N
Y(I) = X(I) - S * HZ(I)
200 CONTINUE
* End of HPROD
END
SUBROUTINE LSTP ( M, N, NDUPLC, NPOWER, DAMP, X,
$ B, D, HY, HZ, W, ACOND, RNORM )
INTEGER M, N, MAXMN, NDUPLC, NPOWER
DOUBLE PRECISION DAMP, ACOND, RNORM
DOUBLE PRECISION B(M), X(N), D(N), HY(M), HZ(N), W(M)
* ------------------------------------------------------------------
* LSTP generates a sparse least-squares test problem of the form
*
* ( A )*X = ( B )
* ( DAMP*I ) ( 0 )
*
* having a specified solution X. The matrix A is constructed
* in the form A = HY*D*HZ, where D is an M by N diagonal matrix,
* and HY and HZ are Householder transformations.
*
* The first 6 parameters are input to LSTP. The remainder are
* output. LSTP is intended for use with M .GE. N.
*
*
* Functions and subroutines
*
* TESTPROB APROD1, HPROD
* BLAS DNRM2
* ------------------------------------------------------------------
* Intrinsics and local variables
INTRINSIC COS, SIN, SQRT
INTEGER I, J
DOUBLE PRECISION DNRM2
DOUBLE PRECISION ALFA, BETA, DAMPSQ, FOURPI, T
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0, ONE = 1.0 )
* ------------------------------------------------------------------
* Make two vectors of norm 1.0 for the Householder transformations.
* FOURPI need not be exact.
* ------------------------------------------------------------------
DAMPSQ = DAMP**2
FOURPI = 4.0 * 3.141592
ALFA = FOURPI / M
BETA = FOURPI / N
DO 100 I = 1, M
HY(I) = SIN( I * ALFA )
100 CONTINUE
DO 200 I = 1, N
HZ(I) = COS( I * BETA )
200 CONTINUE
ALFA = DNRM2 ( M, HY, 1 )
BETA = DNRM2 ( N, HZ, 1 )
CALL DSCAL ( M, (- ONE / ALFA), HY, 1 )
CALL DSCAL ( N, (- ONE / BETA), HZ, 1 )
*
* ------------------------------------------------------------------
* Set the diagonal matrix D. These are the singular values of A.
* ------------------------------------------------------------------
DO 300 I = 1, N
J = (I - 1 + NDUPLC) / NDUPLC
T = J * NDUPLC
T = T / N
D(I) = T**NPOWER
300 CONTINUE
ACOND = SQRT( (D(N)**2 + DAMPSQ) / (D(1)**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 )
DO 500 I = 1, N
B(I) = DAMPSQ * B(I) / D(I)
500 CONTINUE
T = ONE
DO 600 I = N + 1, M
J = I - N
B(I) = (T * J) / M
T = - T
600 CONTINUE
CALL HPROD ( M, HY, B, B )
* ------------------------------------------------------------------
* Now compute the true B = RESIDUAL + A*X.
* ------------------------------------------------------------------
RNORM = SQRT( DNRM2 ( M, B, 1 )**2
$ + DAMPSQ * DNRM2 ( N, X, 1 )**2 )
CALL APROD1( M, N, X, B, D, HY, HZ, W )
* End of LSTP
END
SUBROUTINE TEST ( M, N, NDUPLC, NPOWER, DAMP )
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 )
DOUBLE PRECISION B(MAXM), U(MAXM),
$ V(MAXN), W(MAXN), X(MAXN),
$ SE(MAXN), XTRUE(MAXN)
DOUBLE PRECISION ATOL, BTOL, CONLIM,
$ ANORM, ACOND, RNORM, ARNORM,
$ DAMPSQ, ENORM, ETOL, XNORM
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
DATA LINE
$ /'----------------------------------'/
* Set the desired solution XTRUE.
DO 100 J = 1, N
XTRUE(J) = N - J
100 CONTINUE
* 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
LOCHY = LOCD + N
LOCHZ = LOCHY + M
LOCW = LOCHZ + N
LTOTAL = LOCW + MAX(M,N) - 1
IF (LTOTAL .GT. LENRW) GO TO 900
CALL LSTP ( M, N, NDUPLC, NPOWER, DAMP, XTRUE,
$ B, RW(LOCD), RW(LOCHY), RW(LOCHZ), RW(LOCW),
$ 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 )
ATOL = 1.0E-10
BTOL = ATOL
CONLIM = 10.0 * ACOND
ITNLIM = M + N + 50
NOUT = 6
WRITE(NOUT, 1000) LINE, LINE,
$ M, N, NDUPLC, NPOWER, DAMP, ACOND, RNORM,
$ LINE, LINE
CALL LSQR ( M, N, APROD, DAMP,
$ LENIW, LENRW, IW, RW,
$ U, V, W, X, SE,
$ ATOL, BTOL, CONLIM, ITNLIM, NOUT,
$ ISTOP, ITN, ANORM, ACOND, RNORM, ARNORM, XNORM )
* 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
WRITE(NOUT, 2000)
WRITE(NOUT, 2100) 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 )
CALL DSCAL ( M, (-ONE), U, 1 )
CALL APROD ( 1, M, N, X, U, LENIW, LENRW, IW, RW )
* 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 )
CALL DSCAL ( N, DAMPSQ, V, 1 )
CALL APROD ( 2, M, N, V, U, LENIW, LENRW, IW, RW )
* Compute the norms associated with X, U, V.
XNORM = DNRM2 ( N, X, 1 )
RNORM = SQRT( DNRM2 ( M, U, 1 )**2 + DAMPSQ * XNORM**2 )
ARNORM = DNRM2 ( N, V, 1 )
WRITE(NOUT, 2200) RNORM, ARNORM, XNORM
* Print the solution and standard error estimates from LSQR.
WRITE(NOUT, 2500) (J, X(J), J = 1, N)
WRITE(NOUT, 2600) (J, SE(J), J = 1, N)
* Print a clue about whether the solution looks OK.
DO 500 J = 1, N
W(J) = X(J) - XTRUE(J)
500 CONTINUE
ENORM = DNRM2 ( N, W, 1 ) / (ONE + DNRM2 ( N, XTRUE, 1 ))
ETOL = 1.0D-5
IF (ENORM .LE. ETOL) WRITE(NOUT, 3000) ENORM
IF (ENORM .GT. ETOL) WRITE(NOUT, 3100) ENORM
RETURN
* Not enough workspace.
900 WRITE(NOUT, 9000) LTOTAL
RETURN
1000 FORMAT(1P
$ // 1X, 2A
$ / ' Least-Squares Test Problem P(', 4I5, E12.2, ' )'
$ // ' Condition no. =', E12.4, ' Residual function =', E17.9
$ / 1X, 2A)
2000 FORMAT(
$ // 22X, ' Residual norm Residual norm Solution norm'
$ / 22X, '(Abar X - bbar) (Normal eqns) (X)' /)
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))
3000 FORMAT(1P / ' LSQR appears to be successful.',
$ ' Relative error in X =', E10.2)
3100 FORMAT(1P / ' LSQR appears to have failed. ',
$ ' Relative error in X =', E10.2)
9000 FORMAT(/ ' XXX Insufficient workspace.',
$ ' The length of RW should be at least', I6)
* End of TEST
END
* -------------
* Main program.
* -------------
DOUBLE PRECISION DAMP1, DAMP2, DAMP3, DAMP4, ZERO
*
ZERO = 0.0
DAMP1 = 0.1
DAMP2 = 0.01
DAMP3 = 0.001
DAMP4 = 0.0001
CALL TEST ( 1, 1, 1, 1, ZERO )
CALL TEST ( 2, 1, 1, 1, ZERO )
CALL TEST ( 40, 40, 4, 4, ZERO )
CALL TEST ( 40, 40, 4, 4, DAMP2 )
CALL TEST ( 80, 40, 4, 4, DAMP2 )
STOP
* End of main program for testing LSQR
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -