SUBROUTINE SYMMLQ( N, B, R1, R2, V, W, X, Y,
1 APROD, MSOLVE, PRECON, SHIFT,
2 NOUT, ITNLIM, EPS, RTOL,
3 ISTOP, ITN, ANORM, ACOND, RNORM, XNORM )
C
EXTERNAL APROD, MSOLVE
INTEGER N, NOUT, ITNLIM, ISTOP, ITN
LOGICAL PRECON
DOUBLE PRECISION SHIFT, EPS, RTOL, ANORM, ACOND, RNORM, XNORM
DOUBLE PRECISION B(N), R1(N), R2(N), V(N), W(N), X(N), Y(N)
C ------------------------------------------------------------------
C
C SYMMLQ is designed to solve the system of linear equations
C
C A*x = b
C
C where A is an n*n symmetric matrix and b is a given vector.
C The matrix A is not required to be positive definite.
C (If A is known to be definite, the method of conjugate gradients
C may be used -- it will require about the same number of iterations
C as SYMMLQ but slightly less work per iteration.)
C
C
C The matrix A is intended to be large and sparse. It is accessed
C by means of a subroutine call of the form
C
C CALL APROD( n,x,y )
C
C which must return the product y = A*x for any given vector x.
C
C
C More generally, SYMMLQ is designed to solve the system
C
C (A - shift*I)*x = b
C
C where shift is a specified scalar value. If shift and b
C are suitably chosen, the computed vector x may approximate an
C (unnormalized) eigenvector of A, as in the methods of
C inverse iteration and/or Rayleigh-quotient iteration.
C Again, the matrix (A - shift*I) need not be positive definite.
C The work per iteration is very slightly less if shift = 0.
C
C
C A further option is that of preconditioning, which may reduce
C the number of iterations required. If M is a positive
C definite matrix that is known to approximate (A - shift*I)
C in some sense, and if systems of the form M*y = x can be
C solved efficiently, the parameters PRECON and MSOLVE may be
C used (see below). When PRECON = .true. , SYMMLQ will
C implicitly solve the system of equations
C
C P*( A - shift*I )*P*y = P*b, where P = M**(-1/2),
C
C and return the solution x = P*y.
C
C
C Parameters
C ----------
C
C N input n, the dimension of the matrix A.
C
C B(N) input The rhs vector b.
C
C R1(N) output Returns the final residual vector,
C r = b - (A - shift*I)*x.
C
C R2(N) workspace
C V(N) workspace
C W(N) workspace
C
C X(N) output Returns the computed solution x.
C
C Y(N) workspace
C
C APROD external A subroutine defining the matrix A.
C For a given vector x, the statement
C
C CALL APROD( n,x,y )
C
C must return the product y = A*x
C without altering the vector x.
C
C MSOLVE external An optional subroutine defining a
C preconditioning matrix M, which should
C approximate (A - shift*I) in some sense ---
C typically having many similar eigenvalues.
C M must be positive definite.
C For a given vector x, the statement
C
C CALL MSOLVE( n,x,y )
C
C must solve the linear system M*y = x
C without altering the vector x.
C
C NOTE. The program calling SYMMLQ must
C declare APROD and MSOLVE to be external.
C
C PRECON input If PRECON = .true. , preconditioning will
C be invoked. Otherwise, subroutine MSOLVE
C will not be referenced; in this case the
C actual parameter corresponding to MSOLVE may
C be the same as that corresponding to APROD.
C
C SHIFT input Should be zero if the system A*x = b is to
C be solved. Otherwise, it could be an
C approximation to an eigenvalue of A, such as
C the Rayleigh quotient b(t)*A*b/(b(t)*b)
C corresponding to the vector b.
C If b is sufficiently like an eigenvector
C corresponding to an eigenvalue near SHIFT,
C then the computed x may have very large
C components. When normalized, x may be
C closer to an eigenvector than b.
C
C NOUT input The output file number. If positive,
C a summary will be printed on unit NOUT.
C
C ITNLIM input An upper limit on the number of iterations.
C
C EPS input The relative machine precision.
C For example,
C Burroughs B6700 EPS = 2.0**(-37)
C CDC 6600, 7600 EPS = 2.0**(-47)
C IBM 370 (single) EPS = 16.0**(-5)
C IBM 370 (double) EPS = 16.0**(-13)
C
C RTOL input A user-specified tolerance. SYMMLQ terminates
C if it appears that norm(Rbar) is smaller than
C RTOL * norm(Abar) * norm(y), where
C Abar is the transformed matrix operator
C
C Abar = P * (A - shift*I) * P
C
C and Rbar is the transformed residual vector
C
C Rbar = P * ( b - (A - shift*I)*x ).
C
C If shift = 0 and PRECON = .false., SYMMLQ
C terminates if norm(b - A*x) is smaller than
C RTOL * norm(A) * norm(x).
C
C ISTOP output An integer giving the reason for termination...
C
C 0 b = 0, so the exact solution is x = 0.
C No iterations were performed.
C
C 1 Norm(Rbar) appears to be less than
C the value RTOL * norm(Abar) * norm(y).
C The solution in X should be acceptable.
C
C 2 Norm(Rbar) appears to be less than
C the value EPS * norm(Abar) * norm(y).
C This means that the residual is as small as
C seems reasonable on this machine.
C
C 3 Norm(Abar)*norm(y) exceeds norm(b)/EPS,
C which should indicate that x has essentially
C converged to an eigenvector of A
C corresponding to the eigenvalue shift.
C
C 4 ACOND (see below) has exceeded 0.1/EPS, so
C the matrix Abar must be very ill-conditioned.
C X may not contain an acceptable solution.
C
C 5 The iteration limit was reached before any of
C the previous criteria were satisfied.
C
C 6 An inner product of the form x(t)*M**(-1)*x
C was not positive, so the preconditioning matrix
C M does not appear to be positive definite.
C X will not contain an acceptable solution.
C
C ITN output The number of iterations performed.
C
C ANORM output An estimate of the norm of the matrix operator
C Abar = P*(A - shift*I)*P, where P = M**(-1/2).
C
C ACOND output An estimate of the condition of Abar above.
C This will usually be a substantial
C under-estimate of the true condition.
C
C RNORM output The norm of the final residual vector,
C b - (A - shift*I)*x.
C
C XNORM output The norm of the final solution vector x.
C
C
C
C To change precision
C -------------------
C
C Alter the words
C DABS, DMAX1, DMIN1, DOUBLE PRECISION,
C DAXPY, DCOPY, DDOT, DNRM2, DSQRT
C throughout.
C
C Also make sure EPS is set correctly in the calling program.
C ------------------------------------------------------------------
C
C
C This routine is an implementation of the algorithm described in
C the following references:
C
C C.C. Paige and M.A. Saunders, Solution of Sparse Indefinite
C Systems of Linear Equations,
C SIAM J. Numer. Anal. 12, 4, September 1975, pp. 617-629.
C
C J.G. Lewis, Algorithms for Sparse Matrix Eigenvalue Problems,
C Report STAN-CS-77-595, Computer Science Department,
C Stanford University, Stanford, California, March 1977.
C
C Applications of SYMMLQ and the theory of preconditioning
C are described in the following references:
C
C D.B. Szyld and O.B. Widlund, Applications of Conjugate Gradient
C Type Methods to Eigenvalue Calculations,
C in R. Vichnevetsky and R.S. Steplman (editors),
C Advances in Computer Methods for Partial Differential
C Equations -- III, IMACS, 1979, 167-173.
C
C D.B. Szyld, A Two-level Iterative Method for Large Sparse
C Generalized Eigenvalue Calculations,
C Ph. D. dissertation, Department of Mathematics,
C New York University, New York, October 1983.
C ------------------------------------------------------------------
C
C
C SYMMLQ. This version dated 15 Sept 1985.
C
C Michael A. Saunders
C Department of Operations Research
C Stanford University
C ------------------------------------------------------------------
C
C
C Subroutines and functions
C
C USER APROD, MSOLVE
C BLAS DAXPY, DCOPY, DDOT , DNRM2
C FORTRAN DABS , DMAX1, DMIN1, DSQRT, MOD
C
C
C Functions and local variables
C
DOUBLE PRECISION DDOT, DNRM2
DOUBLE PRECISION DABS, DMAX1, DMIN1, DSQRT
INTEGER MOD
DOUBLE PRECISION S,T,Z,B1,CS,SN,ONE,ALFA,BETA,DBAR,DIAG,
1 EPSA,EPSR,EPSX,GBAR,GMAX,GMIN,
2 OLDB,RHS1,RHS2,X1CG,X1LQ,ZBAR,ZERO,
3 BETA1,BSTEP,DELTA,DENOM,EPSLN,GAMMA,GPERT,
4 TNORM,YNORM,
5 CGNORM,ELQNRM,QRNORM,SNPROD,YNORM2
INTEGER I
C ------------------------------------------------------------------
C
C
C Print heading and initialize.
C
IF (NOUT .GT. 0)
1 WRITE(NOUT, 1000) N, PRECON, SHIFT, NOUT, ITNLIM, EPS, RTOL
ZERO = 0.0
ONE = 1.0
ISTOP = 0
ITN = 0
ANORM = ZERO
ACOND = ZERO
RNORM = ZERO
XNORM = ZERO
YNORM = ZERO
C
DO 50 I = 1, N
X(I) = ZERO
W(I) = ZERO
50 CONTINUE
C
C Set up V, the first vector in the Lanczos sequence.
C
CALL DCOPY ( N,B,1,Y,1 )
CALL DCOPY ( N,B,1,R1,1 )
IF (PRECON) CALL MSOLVE( N,R1,Y )
B1 = Y(1)
BETA1 = DDOT( N,R1,1,Y,1 )
IF (BETA1 .LT. ZERO) ISTOP = 6
IF (BETA1 .LE. ZERO) GO TO 900
BETA1 = DSQRT( BETA1 )
S = ONE/BETA1
C
DO 100 I = 1, N
V(I) = S*Y(I)
100 CONTINUE
C
C Set up Y for the second Lanczos vector.
C
CALL APROD ( N,V,Y )
CALL DAXPY ( N,(-SHIFT),V,1,Y,1 )
ALFA = DDOT( N,V,1,Y,1 )
CALL DAXPY ( N,(-ALFA/BETA1),R1,1,Y,1 )
C
C Make sure R2 will be orthogonal to the first V.
C
Z = DDOT( N,V,1,Y,1 )
S = DDOT( N,V,1,V,1 )
CALL DAXPY ( N,(-Z/S),V,1,Y,1 )
CALL DCOPY ( N,Y,1,R2,1 )
IF (PRECON) CALL MSOLVE( N,R2,Y )
OLDB = BETA1
BETA = DDOT( N,R2,1,Y,1 )
IF (BETA .LE. ZERO) ISTOP = 6
IF (BETA .LE. ZERO) GO TO 900
BETA = DSQRT( BETA )
C
C See if the local reorthogonalization achieved anything.
C
DENOM = DSQRT(S) * DNRM2( N,R2,1 )
S = Z/DENOM
T = DDOT( N,V,1,R2,1 )/DENOM
IF (NOUT .GT. 0) WRITE(NOUT, 1100) BETA1,ALFA
IF (NOUT .GT. 0) WRITE(NOUT, 1120) S,T
C
C Initialize other quantities.
C
CGNORM = BETA1
GBAR = ALFA
DBAR = BETA
RHS1 = BETA1
RHS2 = ZERO
BSTEP = ZERO
SNPROD = ONE
TNORM = ALFA**2
YNORM2 = ZERO
GMAX = ZERO
GMIN = ONE/EPS
C
C ------------------------------------------------------------------
C Main iteration loop.
C ------------------------------------------------------------------
C
C Estimate various norms and test for convergence.
C
300 ANORM = DSQRT(TNORM)
YNORM = DSQRT(YNORM2)
EPSA = ANORM*EPS
EPSX = ANORM*YNORM*EPS
EPSR = ANORM*YNORM*RTOL
DIAG = GBAR
IF (DIAG .EQ. ZERO) DIAG = EPSA
C
ELQNRM = DSQRT(RHS1**2 + RHS2**2)
QRNORM = SNPROD*BETA1
CGNORM = QRNORM*BETA / DABS(DIAG)
C
C Estimate COND(A).
C In this version we look at the diagonals of L in the
C factorization of the tridiagonal matrix, T = L*Q.
C
DENOM = DMIN1( GMIN, DABS(DIAG) )
ACOND = GMAX/DENOM
C
C See if any of the stopping criteria are satisfied.
C
IF (ITN .GE. ITNLIM ) ISTOP = 5
IF (ACOND .GE. 0.1/EPS) ISTOP = 4
IF (EPSX .GE. BETA1 ) ISTOP = 3
IF (CGNORM .LE. EPSX ) ISTOP = 2
IF (CGNORM .LE. EPSR ) ISTOP = 1
C ==================================================================
C
C See if it is time to print something.
C
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 (CGNORM .LE. 10.0*EPSX) GO TO 400
IF (CGNORM .LE. 10.0*EPSR) GO TO 400
IF (ACOND .GE. 0.01/EPS ) GO TO 400
IF (ISTOP .NE. 0) GO TO 400
GO TO 600
C
C Print a line for this iteration.
C
400 ZBAR = RHS1/DIAG
Z = (SNPROD*ZBAR + BSTEP)/BETA1
X1LQ = X(1) + B1*BSTEP/BETA1
X1CG = X(1) + W(1)*ZBAR + B1*Z
C
IF ( ITN .EQ. 0) WRITE(NOUT, 1200)
WRITE(NOUT, 1300) ITN,X1CG,CGNORM,BSTEP,ANORM,ACOND
IF (MOD(ITN,10) .EQ. 0) WRITE(NOUT, 1500)
C ==================================================================
C
C
C Obtain the current Lanczos vector V = (1/BETA)*Y
C and set up Y for the next iteration.
C
600 IF (ISTOP .NE. 0) GO TO 800
S = ONE/BETA
C
DO 620 I = 1, N
V(I) = S*Y(I)
620 CONTINUE
C
CALL APROD ( N,V,Y )
CALL DAXPY ( N,(-SHIFT),V,1,Y,1 )
CALL DAXPY ( N,(-BETA/OLDB),R1,1,Y,1 )
ALFA = DDOT( N,V,1,Y,1 )
TNORM = TNORM + (ALFA**2) + 2.0*(BETA**2)
CALL DAXPY ( N,(-ALFA/BETA),R2,1,Y,1 )
CALL DCOPY ( N,R2,1,R1,1 )
CALL DCOPY ( N,Y,1,R2,1 )
IF (PRECON) CALL MSOLVE( N,R2,Y )
OLDB = BETA
BETA = DDOT( N,R2,1,Y,1 )
IF (BETA .LE. ZERO) ISTOP = 6
IF (BETA .LE. ZERO) GO TO 800
BETA = DSQRT( BETA )
C
C Compute the next plane rotation for Q.
C
GAMMA = DSQRT(GBAR**2 + OLDB**2)
CS = GBAR/GAMMA
SN = OLDB/GAMMA
DELTA = CS*DBAR + SN*ALFA
GBAR = SN*DBAR - CS*ALFA
EPSLN = SN*BETA
DBAR = - CS*BETA
C
C Update X.
C
Z = RHS1/GAMMA
S = Z*CS
T = Z*SN
C
DO 700 I = 1, N
X(I) = (W(I)*S + V(I)*T) + X(I)
W(I) = W(I)*SN - V(I)*CS
700 CONTINUE
C
C Accumulate the step along the direction B,
C and go round again.
C
BSTEP = SNPROD*CS*Z + BSTEP
SNPROD = SNPROD*SN
GMAX = DMAX1( GMAX, GAMMA )
GMIN = DMIN1( GMIN, GAMMA )
YNORM2 = Z**2 + YNORM2
RHS1 = RHS2 - DELTA*Z
RHS2 = - EPSLN*Z
ITN = ITN + 1
GO TO 300
C
C ------------------------------------------------------------------
C End of main iteration loop.
C ------------------------------------------------------------------
C
C Move to the CG point.
C (In this version of SYMMLQ, we never stop at an LQ point.)
C
800 ZBAR = RHS1/DIAG
BSTEP = SNPROD*ZBAR + BSTEP
YNORM = DSQRT(YNORM2 + ZBAR**2)
CALL DAXPY ( N,ZBAR,W,1,X,1 )
C
C Add the step along B.
C
BSTEP = BSTEP/BETA1
CALL DCOPY ( N,B,1,Y,1 )
IF (PRECON) CALL MSOLVE( N,B,Y )
CALL DAXPY ( N,BSTEP,Y,1,X,1 )
C
C Compute the final residual, R1 = B - (A - SHIFT*I)*X.
C
CALL APROD ( N,X,Y )
CALL DAXPY ( N,(-SHIFT),X,1,Y,1 )
C
DO 850 I = 1, N
R1(I) = B(I) - Y(I)
850 CONTINUE
C
RNORM = DNRM2( N,R1,1 )
XNORM = DNRM2( N,X ,1 )
C
C ==================================================================
C Display final status.
C ==================================================================
900 IF (NOUT .LE. 0) GO TO 950
WRITE(NOUT, 2000) ISTOP, ITN, ANORM, ACOND, RNORM, XNORM
IF (ISTOP .EQ. 0) WRITE(NOUT, 3000)
IF (ISTOP .EQ. 1) WRITE(NOUT, 3100)
IF (ISTOP .EQ. 2) WRITE(NOUT, 3200)
IF (ISTOP .EQ. 3) WRITE(NOUT, 3300)
IF (ISTOP .EQ. 4) WRITE(NOUT, 3400)
IF (ISTOP .EQ. 5) WRITE(NOUT, 3500)
IF (ISTOP .EQ. 6) WRITE(NOUT, 3600)
950 RETURN
C
C ------------------------------------------------------------------
1000 FORMAT(// ' SYMMLQ. Solution of symmetric Ax = b'
1 / ' N =', I7, 6X, 'PRECON =', L6, 6X, 'SHIFT =', 1PE23.14
2 / ' NOUT =', I7, 6X, 'ITNLIM =', I6, 6X, 'EPS =', 1PE11.2, 5X,
3 ' RTOL =', 1PE11.2)
1100 FORMAT(/ ' beta1 =', 1PE12.2 / ' alpha1 =', 1PE12.2)
1120 FORMAT(/ ' (V1,V2) before and after ', 1PE15.2
1 / ' local reorthogonalization', 1PE15.2 )
1200 FORMAT(// 5X, 'ITN', 7X, 'X1(CG)',
1 10X, 'Norm(R)', 5X, 'BSTEP', 7X, 'Norm(A)', 3X, 'Cond(A)')
1300 FORMAT(I8, 1PE19.10, 1PE11.2, 1PE14.5, 1P2E10.2)
1500 FORMAT(1X)
2000 FORMAT(/ ' EXIT SYMMLQ.', 14X, 'ISTOP =', I3, 18X, 'ITN =', I8
1 / ' ANORM =', 1PE13.5, 6X, 'ACOND =', 1PE13.5, 5X,
2 ' RNORM =', 1PE13.5, 6X, 'XNORM =', 1PE13.5)
3000 FORMAT(/ ' The exact solution is X = 0.')
3100 FORMAT(/ ' Requested accuracy achieved, as determined by RTOL.')
3200 FORMAT(/ ' Reasonable accuracy achieved, given EPS.')
3300 FORMAT(/ ' X has converged to an eigenvector.')
3400 FORMAT(/ ' ACOND has exceeded 0.1/EPS.')
3500 FORMAT(/ ' The iteration limit was reached.')
3600 FORMAT(/ ' MSOLVE does not define a positive definite',
1 ' preconditioner.')
C ------------------------------------------------------------------
C END OF SYMMLQ
END
C ******************************************************
C
C WARNING. Delete the following imitation BLAS routines
C if a genuine BLAS library is available.
C
C ******************************************************
C
SUBROUTINE DAXPY ( N,A,X,INCX,Y,INCY )
DOUBLE PRECISION A,X(N),Y(N)
C
C This may be replaced by the corresponding BLAS routine.
C The following is a simple version for use with SYMMLQ.
C
IF (A .EQ. 0.0) RETURN
DO 10 I = 1, N
Y(I) = A*X(I) + Y(I)
10 CONTINUE
RETURN
END
SUBROUTINE DCOPY ( N,X,INCX,Y,INCY )
DOUBLE PRECISION X(N),Y(N)
C
C This may be replaced by the corresponding BLAS routine.
C The following is a simple version for use with SYMMLQ.
C
DO 10 I = 1, N
Y(I) = X(I)
10 CONTINUE
RETURN
END
FUNCTION DDOT ( N,X,INCX,Y,INCY )
DOUBLE PRECISION DDOT,X(N),Y(N),SUM
C
C This may be replaced by the corresponding BLAS routine.
C The following is a simple version for use with SYMMLQ.
C
SUM = 0.0
DO 10 I = 1, N
SUM = X(I)*Y(I) + SUM
10 CONTINUE
DDOT = SUM
RETURN
END
FUNCTION DNRM2 ( N,X,INCX )
DOUBLE PRECISION DNRM2,X(N),DSQRT,SUM
C
C This may be replaced by the corresponding BLAS routine.
C The following is a simple version for use with SYMMLQ.
C
SUM = 0.0
DO 10 I = 1, N
SUM = X(I)**2 + SUM
10 CONTINUE
DNRM2 = DSQRT(SUM)
RETURN
END
C ***************************************************
C
C The remaining routines are for testing SYMMLQ.
C
C ***************************************************
C
SUBROUTINE APROD ( N,X,Y )
IMPLICIT REAL*8(A-H,O-Z)
DOUBLE PRECISION X(N),Y(N)
C
C APROD computes Y = A*X for some matrix A.
C This is a simple example for testing SYMMLQ.
C
DO 10 I = 1, N
D = I * 1.01
D = D / N
Y(I) = D * X(I)
10 CONTINUE
RETURN
C
C END OF APROD
END
SUBROUTINE MSOLVE( N,X,Y )
IMPLICIT REAL*8(A-H,O-Z)
DOUBLE PRECISION X(N),Y(N)
COMMON /MSHIFT/ SHIFTM
C
C MSOLVE solves M*Y = X for some symmetric pos-def matrix M.
C This is a simple example for testing SYMMLQ.
C
DO 10 I = 1, N
D = I * 1.1
D = D / N
D = DABS(D - SHIFTM)
Y(I) = X(I) / D
10 CONTINUE
RETURN
C
C END OF MSOLVE
END
SUBROUTINE TEST ( N, PRECON, SHIFT )
IMPLICIT REAL*8(A-H,O-Z)
LOGICAL PRECON
COMMON /MSHIFT/ SHIFTM
C
EXTERNAL APROD, MSOLVE
DOUBLE PRECISION B(100), R1(100), R2(100), V(100), W(100)
DOUBLE PRECISION X(100), Y(100), XTRUE(100)
C
C
SHIFTM = SHIFT
NOUT = 6
WRITE(NOUT, 1000)
C
C Compute machine precision. The call to DAXPY is
C intended to fool compilers that use extra-length registers.
C
ONE = 1
TWO = 2
EPS = ONE
C
10 EPS = EPS/TWO
X(1) = EPS
Y(1) = ONE
CALL DAXPY ( 1, ONE, X, 1, Y, 1 )
IF (Y(1) .GT. ONE) GO TO 10
EPS = EPS*TWO
C
C Set the true solution and the rhs
C so that (A - shift*I)*XTRUE = B.
C
DO 100 I = 1, N
XTRUE(I) = N + 1 - I
100 CONTINUE
CALL APROD ( N,XTRUE,B )
CALL DAXPY ( N,(-SHIFT),XTRUE,1,B,1 )
C
C Set other parameters and solve.
C
ITNLIM = N*2
RTOL = EPS*10.0
C
CALL SYMMLQ( N,B,R1,R2,V,W,X,Y,
1 APROD,MSOLVE,PRECON,SHIFT,
2 NOUT,ITNLIM,EPS,RTOL,
3 ISTOP,ITN,ANORM,ACOND,RNORM,XNORM )
C
C Print the solution and some clue about whether it is OK.
C
WRITE(NOUT, 2000) (J,X(J), J=1,N)
C
DO 500 J = 1, N
W(J) = X(J) - XTRUE(J)
500 CONTINUE
ENORM = DNRM2 ( N, W, 1 ) / 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
C
1000 FORMAT('1' / ' Test of SYMMLQ.')
2000 FORMAT(/ ' Solution X' / 5(I6, 1PE18.10))
3000 FORMAT(/ ' The test of SYMMLQ appears to be successful.',
1 ' Relative error in X =', 1PE10.2)
3100 FORMAT(/ ' The test of SYMMLQ appears to have failed. ',
1 ' Relative error in X =', 1PE10.2)
C End of TEST
END
C -------------
C Main program.
C -------------
IMPLICIT REAL*8(A-H,O-Z)
LOGICAL NORMAL, PRECON
C
N = 50
NORMAL = .FALSE.
PRECON = .TRUE.
ZERO = 0.0
ONE = 1.0
SHIFT = ONE / 9.0
C
CALL TEST( N, NORMAL, ZERO )
CALL TEST( N, PRECON, ZERO )
CALL TEST( N, NORMAL, SHIFT )
CALL TEST( N, PRECON, SHIFT )
STOP
END