dir.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	dir (N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MSOLVE, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, Z, DZ, RWORK, IWORK)
integer function	isdir (N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, Z, DZ, RWORK, IWORK, BNRM, SOLNRM)

Function/Subroutine Documentation

dir()

subroutine dir	(	integer	N,
		double precision, dimension(n)	B,
		double precision, dimension(n)	X,
		integer	NELT,
		integer, dimension(nelt)	IA,
		integer, dimension(nelt)	JA,
		double precision, dimension(nelt)	A,
		integer	ISYM,
		external	MATVEC,
		external	MSOLVE,
		integer	ITOL,
		double precision	TOL,
		integer	ITMAX,
		integer	ITER,
		double precision	ERR,
		integer	IERR,
		integer	IUNIT,
		double precision, dimension(n)	R,
		double precision, dimension(n)	Z,
		double precision, dimension(n)	DZ,
		double precision, dimension(*)	RWORK,
		integer, dimension(*)	IWORK
	)

C***BEGIN PROLOGUE  DIR
C***PURPOSE  Preconditioned Iterative Refinement Sparse Ax = b Solver.
C            Routine to solve a general linear system  Ax = b  using
C            iterative refinement with a matrix splitting.
C***LIBRARY   SLATEC (SLAP)
C***CATEGORY  D2A4, D2B4
C***TYPE      DOUBLE PRECISION (SIR-S, DIR-D)
C***KEYWORDS  ITERATIVE PRECONDITION, LINEAR SYSTEM, SLAP, SPARSE
C***AUTHOR  Greenbaum, Anne, (Courant Institute)
C           Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-60
C             Livermore, CA 94550 (510) 423-3141
C             seager@llnl.gov
C***DESCRIPTION
C
C *Usage:
C     INTEGER  N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX
C     INTEGER  ITER, IERR, IUNIT, IWORK(USER DEFINED)
C     DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, R(N), Z(N), DZ(N)
C     DOUBLE PRECISION RWORK(USER DEFINED)
C     EXTERNAL MATVEC, MSOLVE
C
C     CALL DIR(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MSOLVE, ITOL,
C    $     TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, Z, DZ, RWORK, IWORK)
C
C *Arguments:
C N      :IN       Integer.
C         Order of the Matrix.
C B      :IN       Double Precision B(N).
C         Right-hand side vector.
C X      :INOUT    Double Precision X(N).
C         On input X is your initial guess for solution vector.
C         On output X is the final approximate solution.
C NELT   :IN       Integer.
C         Number of Non-Zeros stored in A.
C IA     :IN       Integer IA(NELT).
C JA     :IN       Integer JA(NELT).
C A      :IN       Double Precision A(NELT).
C         These arrays contain the matrix data structure for A.
C         It could take any form.  See "Description", below,
C         for more details.
C ISYM   :IN       Integer.
C         Flag to indicate symmetric storage format.
C         If ISYM=0, all non-zero entries of the matrix are stored.
C         If ISYM=1, the matrix is symmetric, and only the upper
C         or lower triangle of the matrix is stored.
C MATVEC :EXT      External.
C         Name of a routine which performs the matrix vector multiply
C         Y = A*X given A and X.  The name of the MATVEC routine must
C         be declared external in the calling program.  The calling
C         sequence to MATVEC is:
C             CALL MATVEC( N, X, Y, NELT, IA, JA, A, ISYM )
C         Where N is the number of unknowns, Y is the product A*X
C         upon return, X is an input vector, NELT is the number of
C         non-zeros in the SLAP IA, JA, A storage for the matrix A.
C         ISYM is a flag which, if non-zero, denotes that A is
C         symmetric and only the lower or upper triangle is stored.
C MSOLVE :EXT      External.
C         Name of a routine which solves a linear system MZ = R for
C         Z given R with the preconditioning matrix M (M is supplied via
C         RWORK and IWORK arrays).  The name of the MSOLVE routine must
C         be declared external in the calling program.  The calling
C         sequence to MSOLVE is:
C             CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK)
C         Where N is the number of unknowns, R is the right-hand side
C         vector and Z is the solution upon return.  NELT, IA, JA, A and
C         ISYM are defined as above.  RWORK is a double precision array
C         that can be used to pass necessary preconditioning information
C         and/or workspace to MSOLVE.  IWORK is an integer work array
C         for the same purpose as RWORK.
C ITOL   :IN       Integer.
C         Flag to indicate type of convergence criterion.
C         If ITOL=1, iteration stops when the 2-norm of the residual
C         divided by the 2-norm of the right-hand side is less than TOL.
C         If ITOL=2, iteration stops when the 2-norm of M-inv times the
C         residual divided by the 2-norm of M-inv times the right hand
C         side is less than TOL, where M-inv is the inverse of the
C         diagonal of A.
C         ITOL=11 is often useful for checking and comparing different
C         routines.  For this case, the user must supply the "exact"
C         solution or a very accurate approximation (one with an error
C         much less than TOL) through a common block,
C             COMMON /DSLBLK/ SOLN( )
C         If ITOL=11, iteration stops when the 2-norm of the difference
C         between the iterative approximation and the user-supplied
C         solution divided by the 2-norm of the user-supplied solution
C         is less than TOL.  Note that this requires the user to set up
C         the "COMMON /DSLBLK/ SOLN(LENGTH)" in the calling routine.
C         The routine with this declaration should be loaded before the
C         stop test so that the correct length is used by the loader.
C         This procedure is not standard Fortran and may not work
C         correctly on your system (although it has worked on every
C         system the authors have tried).  If ITOL is not 11 then this
C         common block is indeed standard Fortran.
C TOL    :INOUT    Double Precision.
C         Convergence criterion, as described above.  (Reset if IERR=4.)
C ITMAX  :IN       Integer.
C         Maximum number of iterations.
C ITER   :OUT      Integer.
C         Number of iterations required to reach convergence, or
C         ITMAX+1 if convergence criterion could not be achieved in
C         ITMAX iterations.
C ERR    :OUT      Double Precision.
C         Error estimate of error in final approximate solution, as
C         defined by ITOL.
C IERR   :OUT      Integer.
C         Return error flag.
C           IERR = 0 => All went well.
C           IERR = 1 => Insufficient space allocated for WORK or IWORK.
C           IERR = 2 => Method failed to converge in ITMAX steps.
C           IERR = 3 => Error in user input.
C                       Check input values of N, ITOL.
C           IERR = 4 => User error tolerance set too tight.
C                       Reset to 500*D1MACH(3).  Iteration proceeded.
C           IERR = 5 => Preconditioning matrix, M, is not positive
C                       definite.  (r,z) < 0.
C           IERR = 6 => Matrix A is not positive definite.  (p,Ap) < 0.
C IUNIT  :IN       Integer.
C         Unit number on which to write the error at each iteration,
C         if this is desired for monitoring convergence.  If unit
C         number is 0, no writing will occur.
C R      :WORK     Double Precision R(N).
C Z      :WORK     Double Precision Z(N).
C DZ     :WORK     Double Precision DZ(N).
C         Double Precision arrays used for workspace.
C RWORK  :WORK     Double Precision RWORK(USER DEFINED).
C         Double Precision array that can be used by  MSOLVE.
C IWORK  :WORK     Integer IWORK(USER DEFINED).
C         Integer array that can be used by  MSOLVE.
C
C *Description:
C       The basic algorithm for iterative refinement (also known as
C       iterative improvement) is:
C
C            n+1    n    -1       n
C           X    = X  + M  (B - AX  ).
C
C           -1   -1
C       If M =  A then this  is the  standard  iterative  refinement
C       algorithm and the "subtraction" in the  residual calculation
C       should be done in double precision (which it is  not in this
C       routine).
C       If M = DIAG(A), the diagonal of A, then iterative refinement
C       is  known  as  Jacobi's  method.   The  SLAP  routine  DSJAC
C       implements this iterative strategy.
C       If M = L, the lower triangle of A, then iterative refinement
C       is known as Gauss-Seidel.   The SLAP routine DSGS implements
C       this iterative strategy.
C
C       This routine does  not care  what matrix data   structure is
C       used for  A and M.  It simply   calls  the MATVEC and MSOLVE
C       routines, with  the arguments as  described above.  The user
C       could write any type of structure and the appropriate MATVEC
C       and MSOLVE routines.  It is assumed  that A is stored in the
C       IA, JA, A  arrays in some fashion and  that M (or INV(M)) is
C       stored  in  IWORK  and  RWORK)  in  some fashion.   The SLAP
C       routines DSJAC and DSGS are examples of this procedure.
C
C       Two  examples  of  matrix  data structures  are the: 1) SLAP
C       Triad  format and 2) SLAP Column format.
C
C       =================== S L A P Triad format ===================
C
C       In  this   format only the  non-zeros are  stored.  They may
C       appear  in *ANY* order.   The user  supplies three arrays of
C       length NELT, where  NELT  is the number  of non-zeros in the
C       matrix:  (IA(NELT), JA(NELT),  A(NELT)).  For each  non-zero
C       the  user puts   the row  and  column index   of that matrix
C       element in the IA and JA arrays.  The  value of the non-zero
C       matrix  element is  placed in  the corresponding location of
C       the A  array.  This is  an extremely easy data  structure to
C       generate.  On  the other hand it  is  not too  efficient  on
C       vector  computers   for the  iterative  solution  of  linear
C       systems.  Hence, SLAP  changes this input  data structure to
C       the SLAP   Column  format for the  iteration (but   does not
C       change it back).
C
C       Here is an example of the  SLAP Triad   storage format for a
C       5x5 Matrix.  Recall that the entries may appear in any order.
C
C           5x5 Matrix      SLAP Triad format for 5x5 matrix on left.
C                              1  2  3  4  5  6  7  8  9 10 11
C       |11 12  0  0 15|   A: 51 12 11 33 15 53 55 22 35 44 21
C       |21 22  0  0  0|  IA:  5  1  1  3  1  5  5  2  3  4  2
C       | 0  0 33  0 35|  JA:  1  2  1  3  5  3  5  2  5  4  1
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C
C       =================== S L A P Column format ==================
C
C       In  this format   the non-zeros are    stored counting  down
C       columns (except  for the diagonal  entry, which must  appear
C       first  in each "column") and are  stored in the  double pre-
C       cision array  A. In  other  words,  for each  column  in the
C       matrix  first put  the diagonal entry in A.  Then put in the
C       other non-zero  elements going  down the column  (except the
C       diagonal)  in order.  The IA array  holds the  row index for
C       each non-zero.  The JA array  holds the offsets into the IA,
C       A  arrays  for  the  beginning  of  each  column.  That  is,
C       IA(JA(ICOL)),A(JA(ICOL)) are the first elements of the ICOL-
C       th column in IA and A, and IA(JA(ICOL+1)-1), A(JA(ICOL+1)-1)
C       are  the last elements of the ICOL-th column.   Note that we
C       always have JA(N+1)=NELT+1, where N is the number of columns
C       in the matrix  and NELT  is the number  of non-zeros  in the
C       matrix.
C
C       Here is an example of the  SLAP Column  storage format for a
C       5x5 Matrix (in the A and IA arrays '|'  denotes the end of a
C       column):
C
C           5x5 Matrix      SLAP Column format for 5x5 matrix on left.
C                              1  2  3    4  5    6  7    8    9 10 11
C       |11 12  0  0 15|   A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35
C       |21 22  0  0  0|  IA:  1  2  5 |  2  1 |  3  5 |  4 |  5  1  3
C       | 0  0 33  0 35|  JA:  1  4  6    8  9   12
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C
C *Examples:
C       See the SLAP routines DSJAC, DSGS
C
C *Cautions:
C     This routine will attempt to write to the Fortran logical output
C     unit IUNIT, if IUNIT .ne. 0.  Thus, the user must make sure that
C     this logical unit is attached to a file or terminal before calling
C     this routine with a non-zero value for IUNIT.  This routine does
C     not check for the validity of a non-zero IUNIT unit number.
C
C***SEE ALSO  DSJAC, DSGS
C***REFERENCES  1. Gene Golub and Charles Van Loan, Matrix Computations,
C                  Johns Hopkins University Press, Baltimore, Maryland,
C                  1983.
C               2. Mark K. Seager, A SLAP for the Masses, in
C                  G. F. Carey, Ed., Parallel Supercomputing: Methods,
C                  Algorithms and Applications, Wiley, 1989, pp.135-155.
C***ROUTINES CALLED  D1MACH, ISDIR
C***REVISION HISTORY  (YYMMDD)
C   890404  DATE WRITTEN
C   890404  Previous REVISION DATE
C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
C   890921  Removed TeX from comments.  (FNF)
C   890922  Numerous changes to prologue to make closer to SLATEC
C           standard.  (FNF)
C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
C   891004  Added new reference.
C   910411  Prologue converted to Version 4.0 format.  (BAB)
C   910502  Removed MATVEC and MSOLVE from ROUTINES CALLED list.  (FNF)
C   920407  COMMON BLOCK renamed DSLBLK.  (WRB)
C   920511  Added complete declaration section.  (WRB)
C   920929  Corrected format of references.  (FNF)
C   921019  Changed 500.0 to 500 to reduce SP/DP differences.  (FNF)
C***END PROLOGUE  DIR
C     .. Scalar Arguments ..
      DOUBLE PRECISION err, tol
      INTEGER ierr, isym, iter, itmax, itol, iunit, n, nelt
C     .. Array Arguments ..
      DOUBLE PRECISION a(nelt), b(n), dz(n), r(n), rwork(*), x(n), z(n)
      INTEGER ia(nelt), iwork(*), ja(nelt)
C     .. Subroutine Arguments ..
      EXTERNAL matvec, msolve
C     .. Local Scalars ..
      DOUBLE PRECISION bnrm, solnrm, tolmin
      INTEGER i, k
C     .. External Functions ..
      DOUBLE PRECISION d1mach
      INTEGER isdir
      EXTERNAL d1mach, isdir
C***FIRST EXECUTABLE STATEMENT  DIR
C
C         Check some of the input data.
C
      iter = 0
      ierr = 0
      IF( n.LT.1 ) THEN
         ierr = 3
         RETURN
      ENDIF
      tolmin = 500*d1mach(3)
      IF( tol.LT.tolmin ) THEN
         tol = tolmin
         ierr = 4
      ENDIF
C
C         Calculate initial residual and pseudo-residual, and check
C         stopping criterion.
      CALL matvec(n, x, r, nelt, ia, ja, a, isym)
      DO 10 i = 1, n
         r(i) = b(i) - r(i)
 10   CONTINUE
      CALL msolve(n, r, z, nelt, ia, ja, a, isym, rwork, iwork)
C
      IF( isdir(n, b, x, nelt, ia, ja, a, isym, msolve, itol, tol,
     $     itmax, iter, err, ierr, iunit, r, z, dz, rwork,
     $     iwork, bnrm, solnrm) .NE. 0 ) GO TO 200
      IF( ierr.NE.0 ) RETURN
C
C         ***** iteration loop *****
C
      DO 100 k=1,itmax
         iter = k
C
C         Calculate new iterate x, new residual r, and new
C         pseudo-residual z.
         DO 20 i = 1, n
            x(i) = x(i) + z(i)
 20      CONTINUE
         CALL matvec(n, x, r, nelt, ia, ja, a, isym)
         DO 30 i = 1, n
            r(i) = b(i) - r(i)
 30      CONTINUE
         CALL msolve(n, r, z, nelt, ia, ja, a, isym, rwork, iwork)
C
C         check stopping criterion.
         IF( isdir(n, b, x, nelt, ia, ja, a, isym, msolve, itol, tol,
     $        itmax, iter, err, ierr, iunit, r, z, dz, rwork,
     $        iwork, bnrm, solnrm) .NE. 0 ) GO TO 200
C
 100  CONTINUE
C
C         *****   end of loop  *****
C         Stopping criterion not satisfied.
      iter = itmax + 1
      ierr = 2
C
 200  RETURN
C------------- LAST LINE OF DIR FOLLOWS -------------------------------

isdir()

integer function isdir	(	integer	N,
		double precision, dimension(n)	B,
		double precision, dimension(n)	X,
		integer	NELT,
		integer, dimension(nelt)	IA,
		integer, dimension(nelt)	JA,
		double precision, dimension(nelt)	A,
		integer	ISYM,
		external	MSOLVE,
		integer	ITOL,
		double precision	TOL,
		integer	ITMAX,
		integer	ITER,
		double precision	ERR,
		integer	IERR,
		integer	IUNIT,
		double precision, dimension(n)	R,
		double precision, dimension(n)	Z,
		double precision, dimension(n)	DZ,
		double precision, dimension(*)	RWORK,
		integer, dimension(*)	IWORK,
		double precision	BNRM,
		double precision	SOLNRM
	)

C***BEGIN PROLOGUE  ISDIR
C***SUBSIDIARY
C***PURPOSE  Preconditioned Iterative Refinement Stop Test.
C            This routine calculates the stop test for the iterative
C            refinement iteration scheme.  It returns a non-zero if the
C            error estimate (the type of which is determined by ITOL)
C            is less than the user specified tolerance TOL.
C***LIBRARY   SLATEC (SLAP)
C***CATEGORY  D2A4, D2B4
C***TYPE      DOUBLE PRECISION (ISSIR-S, ISDIR-D)
C***KEYWORDS  LINEAR SYSTEM, SLAP, SPARSE, STOP TEST
C***AUTHOR  Greenbaum, Anne, (Courant Institute)
C           Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-60
C             Livermore, CA 94550 (510) 423-3141
C             seager@llnl.gov
C***DESCRIPTION
C
C *Usage:
C     INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX, ITER
C     INTEGER IERR, IUNIT, IWORK(USER DEFINED)
C     DOUBLE PRECISION B(N), X(N), A(N), TOL, ERR, R(N), Z(N), DZ(N)
C     DOUBLE PRECISION RWORK(USER DEFINED), BNRM, SOLNRM
C     EXTERNAL MSOLVE
C
C     IF( ISDIR(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL,
C    $     ITMAX, ITER, ERR, IERR, IUNIT, R, Z, DZ, RWORK, IWORK,
C    $     BNRM, SOLNRM) .NE. 0 ) THEN ITERATION DONE
C
C *Arguments:
C N      :IN       Integer.
C         Order of the Matrix.
C B      :IN       Double Precision B(N).
C         Right-hand side vector.
C X      :IN       Double Precision X(N).
C         The current approximate solution vector.
C NELT   :IN       Integer.
C         Number of Non-Zeros stored in A.
C IA     :IN       Integer IA(NELT).
C JA     :IN       Integer JA(NELT).
C A      :IN       Double Precision A(NELT).
C         These arrays contain the matrix data structure for A.
C         It could take any form.  See "C *Description" in the
C         DIR routine.
C ISYM   :IN       Integer.
C         Flag to indicate symmetric storage format.
C         If ISYM=0, all non-zero entries of the matrix are stored.
C         If ISYM=1, the matrix is symmetric, and only the upper
C         or lower triangle of the matrix is stored.
C MSOLVE :EXT      External.
C         Name of a routine which solves a linear system Mz = r for
C         z given r with the preconditioning matrix M (M is supplied via
C         RWORK and IWORK arrays.  The name of the MSOLVE routine must
C         be declared external in the calling program.  The calling
C         sequence to MSOLVE is:
C             CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK)
C         Where N is the number of unknowns, R is the right-hand side
C         vector and Z is the solution upon return.  NELT, IA, JA, A and
C         ISYM are defined as above.  RWORK is a double precision array
C         that can be used to pass necessary preconditioning information
C         and/or workspace to MSOLVE.  IWORK is an integer work array
C         for the same purpose as RWORK.
C ITOL   :IN       Integer.
C         Flag to indicate type of convergence criterion.
C         If ITOL=1, iteration stops when the 2-norm of the residual
C         divided by the 2-norm of the right-hand side is less than TOL.
C         If ITOL=2, iteration stops when the 2-norm of M-inv times the
C         residual divided by the 2-norm of M-inv times the right hand
C         side is less than TOL, where M-inv is the inverse of the
C         diagonal of A.
C         ITOL=11 is often useful for checking and comparing different
C         routines.  For this case, the user must supply the "exact"
C         solution or a very accurate approximation (one with an error
C         much less than TOL) through a common block,
C             COMMON /DSLBLK/ SOLN( )
C         If ITOL=11, iteration stops when the 2-norm of the difference
C         between the iterative approximation and the user-supplied
C         solution divided by the 2-norm of the user-supplied solution
C         is less than TOL.  Note that this requires the user to set up
C         the "COMMON /DSLBLK/ SOLN(LENGTH)" in the calling routine.
C         The routine with this declaration should be loaded before the
C         stop test so that the correct length is used by the loader.
C         This procedure is not standard Fortran and may not work
C         correctly on your system (although it has worked on every
C         system the authors have tried).  If ITOL is not 11 then this
C         common block is indeed standard Fortran.
C TOL    :IN       Double Precision.
C         Convergence criterion, as described above.
C ITMAX  :IN       Integer.
C         Maximum number of iterations.
C ITER   :IN       Integer.
C         Current iteration count.  (Must be zero on first call.)
C ERR    :OUT      Double Precision.
C         Error estimate of error in the X(N) approximate solution, as
C         defined by ITOL.
C IERR   :OUT      Integer.
C         Error flag.  IERR is set to 3 if ITOL is not one of the
C         acceptable values, see above.
C IUNIT  :IN       Integer.
C         Unit number on which to write the error at each iteration,
C         if this is desired for monitoring convergence.  If unit
C         number is 0, no writing will occur.
C R      :IN       Double Precision R(N).
C         The residual R = B-AX.
C Z      :WORK     Double Precision Z(N).
C         Workspace used to hold the pseudo-residual M z = r.
C DZ     :WORK     Double Precision DZ(N).
C         Workspace used to hold temporary vector(s).
C RWORK  :WORK     Double Precision RWORK(USER DEFINED).
C         Double Precision array that can be used by  MSOLVE.
C IWORK  :WORK     Integer IWORK(USER DEFINED).
C         Integer array that can be used by MSOLVE.
C BNRM   :INOUT    Double Precision.
C         Norm of the right hand side.  Type of norm depends on ITOL.
C         Calculated only on the first call.
C SOLNRM :INOUT    Double Precision.
C         2-Norm of the true solution, SOLN.  Only computed and used
C         if ITOL = 11.
C
C *Function Return Values:
C       0 : Error estimate (determined by ITOL) is *NOT* less than the
C           specified tolerance, TOL.  The iteration must continue.
C       1 : Error estimate (determined by ITOL) is less than the
C           specified tolerance, TOL.  The iteration can be considered
C           complete.
C
C *Cautions:
C     This routine will attempt to write to the Fortran logical output
C     unit IUNIT, if IUNIT .ne. 0.  Thus, the user must make sure that
C     this logical unit is attached to a file or terminal before calling
C     this routine with a non-zero value for IUNIT.  This routine does
C     not check for the validity of a non-zero IUNIT unit number.
C
C***SEE ALSO  DIR, DSJAC, DSGS
C***ROUTINES CALLED  D1MACH, DNRM2
C***COMMON BLOCKS    DSLBLK
C***REVISION HISTORY  (YYMMDD)
C   871119  DATE WRITTEN
C   880320  Previous REVISION DATE
C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
C   890922  Numerous changes to prologue to make closer to SLATEC
C           standard.  (FNF)
C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
C   891003  Removed C***REFER TO line, per MKS.
C   910411  Prologue converted to Version 4.0 format.  (BAB)
C   910502  Removed MSOLVE from ROUTINES CALLED list.  (FNF)
C   910506  Made subsidiary to DIR.  (FNF)
C   920407  COMMON BLOCK renamed DSLBLK.  (WRB)
C   920511  Added complete declaration section.  (WRB)
C   921026  Changed 1.0E10 to D1MACH(2) and corrected E to D in
C           output format.  (FNF)
C***END PROLOGUE  ISDIR
C     .. Scalar Arguments ..
      DOUBLE PRECISION bnrm, err, solnrm, tol
      INTEGER ierr, isym, iter, itmax, itol, iunit, n, nelt
C     .. Array Arguments ..
      DOUBLE PRECISION a(nelt), b(n), dz(n), r(n), rwork(*), x(n), z(n)
      INTEGER ia(nelt), iwork(*), ja(nelt)
C     .. Subroutine Arguments ..
      EXTERNAL msolve
C     .. Arrays in Common ..
      DOUBLE PRECISION soln(1)
C     .. Local Scalars ..
      INTEGER i
C     .. External Functions ..
      DOUBLE PRECISION d1mach, dnrm2
      EXTERNAL d1mach, dnrm2
C     .. Common blocks ..
      COMMON /dslblk/ soln
C***FIRST EXECUTABLE STATEMENT  ISDIR
      isdir = 0
      IF( itol.EQ.1 ) THEN
C         err = ||Residual||/||RightHandSide|| (2-Norms).
         IF(iter .EQ. 0) bnrm = dnrm2(n, b, 1)
         err = dnrm2(n, r, 1)/bnrm
      ELSE IF( itol.EQ.2 ) THEN
C                  -1              -1
C         err = ||M  Residual||/||M  RightHandSide|| (2-Norms).
         IF(iter .EQ. 0) THEN
            CALL msolve(n, b, dz, nelt, ia, ja, a, isym, rwork, iwork)
            bnrm = dnrm2(n, dz, 1)
         ENDIF
         err = dnrm2(n, z, 1)/bnrm
      ELSE IF( itol.EQ.11 ) THEN
C         err = ||x-TrueSolution||/||TrueSolution|| (2-Norms).
         IF( iter.EQ.0 ) solnrm = dnrm2(n, soln, 1)
         DO 10 i = 1, n
            dz(i) = x(i) - soln(i)
 10      CONTINUE
         err = dnrm2(n, dz, 1)/solnrm
      ELSE
C
C         If we get here ITOL is not one of the acceptable values.
         err = d1mach(2)
         ierr = 3
      ENDIF
C
      IF( iunit.NE.0 ) THEN
         WRITE(iunit,1000) iter,err
      ENDIF
C
      IF( err.LE.tol ) isdir = 1
C
      RETURN
 1000 FORMAT(5x,'ITER = ',i4,' Error Estimate = ',d16.7)
C------------- LAST LINE OF ISDIR FOLLOWS -----------------------------