SLATEC Routines --- ISDGMR ---
*DECK ISDGMR
INTEGER FUNCTION ISDGMR (N, B, X, XL, NELT, IA, JA, A, ISYM,
+ MSOLVE, NMSL, ITOL, TOL, ITMAX, ITER, ERR, IUNIT, R, Z, DZ,
+ RWORK, IWORK, RNRM, BNRM, SB, SX, JSCAL, KMP, LGMR, MAXL,
+ MAXLP1, V, Q, SNORMW, PROD, R0NRM, HES, JPRE)
C***BEGIN PROLOGUE ISDGMR
C***SUBSIDIARY
C***PURPOSE Generalized Minimum Residual Stop Test.
C This routine calculates the stop test for the Generalized
C Minimum RESidual (GMRES) iteration scheme. It returns a
C non-zero if the error estimate (the type of which is
C determined by ITOL) is less than the user specified
C tolerance TOL.
C***LIBRARY SLATEC (SLAP)
C***CATEGORY D2A4, D2B4
C***TYPE DOUBLE PRECISION (ISSGMR-S, ISDGMR-D)
C***KEYWORDS GMRES, LINEAR SYSTEM, SLAP, SPARSE, STOP TEST
C***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov
C Hindmarsh, Alan, (LLNL), alanh@llnl.gov
C Seager, Mark K., (LLNL), seager@llnl.gov
C Lawrence Livermore National Laboratory
C PO Box 808, L-60
C Livermore, CA 94550 (510) 423-3141
C***DESCRIPTION
C
C *Usage:
C INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, NMSL, ITOL
C INTEGER ITMAX, ITER, IUNIT, IWORK(USER DEFINED), JSCAL
C INTEGER KMP, LGMR, MAXL, MAXLP1, JPRE
C DOUBLE PRECISION B(N), X(N), XL(MAXL), A(NELT), TOL, ERR,
C $ R(N), Z(N), DZ(N), RWORK(USER DEFINED),
C $ RNRM, BNRM, SB(N), SX(N), V(N,MAXLP1),
C $ Q(2*MAXL), SNORMW, PROD, R0NRM,
C $ HES(MAXLP1,MAXL)
C EXTERNAL MSOLVE
C
C IF (ISDGMR(N, B, X, XL, NELT, IA, JA, A, ISYM, MSOLVE,
C $ NMSL, ITOL, TOL, ITMAX, ITER, ERR, IUNIT, R, Z, DZ,
C $ RWORK, IWORK, RNRM, BNRM, SB, SX, JSCAL,
C $ KMP, LGMR, MAXL, MAXLP1, V, Q, SNORMW, PROD, R0NRM,
C $ HES, JPRE) .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 Approximate solution vector as of the last restart.
C XL :OUT Double Precision XL(N)
C An array of length N used to hold the approximate
C solution as of the current iteration. Only computed by
C this routine when ITOL=11.
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", in the DGMRES,
C DSLUGM and DSDGMR routines 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 MSOLVE :EXT External.
C Name of a routine which solves a linear system Mz = r for z
C 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 NMSL :INOUT Integer.
C A counter for the number of calls to MSOLVE.
C ITOL :IN Integer.
C Flag to indicate the type of convergence criterion used.
C ITOL=0 Means the iteration stops when the test described
C below on the residual RL is satisfied. This is
C the "Natural Stopping Criteria" for this routine.
C Other values of ITOL cause extra, otherwise
C unnecessary, computation per iteration and are
C therefore much less efficient.
C ITOL=1 Means the iteration stops when the first test
C described below on the residual RL is satisfied,
C and there is either right or no preconditioning
C being used.
C ITOL=2 Implies that the user is using left
C preconditioning, and the second stopping criterion
C below is used.
C ITOL=3 Means the iteration stops when the third test
C described below on Minv*Residual is satisfied, and
C there is either left or no preconditioning begin
C used.
C ITOL=11 is often useful for checking and comparing
C different routines. For this case, the user must
C supply the "exact" solution or a very accurate
C approximation (one with an error much less than
C TOL) through a common block,
C COMMON /DSLBLK/ SOLN( )
C If ITOL=11, iteration stops when the 2-norm of the
C difference between the iterative approximation and
C the user-supplied solution divided by the 2-norm
C of the user-supplied solution is less than TOL.
C Note that this requires the user to set up the
C "COMMON /DSLBLK/ SOLN(LENGTH)" in the calling
C routine. The routine with this declaration should
C be loaded before the stop test so that the correct
C length is used by the loader. This procedure is
C not standard Fortran and may not work correctly on
C your system (although it has worked on every
C system the authors have tried). If ITOL is not 11
C then this 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 The iteration for which to check for convergence.
C ERR :OUT Double Precision.
C Error estimate of error in final approximate solution, as
C defined by ITOL. Letting norm() denote the Euclidean
C norm, ERR is defined as follows..
C
C If ITOL=0, then ERR = norm(SB*(B-A*X(L)))/norm(SB*B),
C for right or no preconditioning, and
C ERR = norm(SB*(M-inverse)*(B-A*X(L)))/
C norm(SB*(M-inverse)*B),
C for left preconditioning.
C If ITOL=1, then ERR = norm(SB*(B-A*X(L)))/norm(SB*B),
C since right or no preconditioning
C being used.
C If ITOL=2, then ERR = norm(SB*(M-inverse)*(B-A*X(L)))/
C norm(SB*(M-inverse)*B),
C since left preconditioning is being
C used.
C If ITOL=3, then ERR = Max |(Minv*(B-A*X(L)))(i)/x(i)|
C i=1,n
C If ITOL=11, then ERR = norm(SB*(X(L)-SOLN))/norm(SB*SOLN).
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 :INOUT Double Precision R(N).
C Work array used in calling routine. It contains
C information necessary to compute the residual RL = B-A*XL.
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 RNRM :IN Double Precision.
C Norm of the current residual. Type of norm depends on ITOL.
C BNRM :IN Double Precision.
C Norm of the right hand side. Type of norm depends on ITOL.
C SB :IN Double Precision SB(N).
C Scaling vector for B.
C SX :IN Double Precision SX(N).
C Scaling vector for X.
C JSCAL :IN Integer.
C Flag indicating if scaling arrays SB and SX are being
C used in the calling routine DPIGMR.
C JSCAL=0 means SB and SX are not used and the
C algorithm will perform as if all
C SB(i) = 1 and SX(i) = 1.
C JSCAL=1 means only SX is used, and the algorithm
C performs as if all SB(i) = 1.
C JSCAL=2 means only SB is used, and the algorithm
C performs as if all SX(i) = 1.
C JSCAL=3 means both SB and SX are used.
C KMP :IN Integer
C The number of previous vectors the new vector VNEW
C must be made orthogonal to. (KMP .le. MAXL)
C LGMR :IN Integer
C The number of GMRES iterations performed on the current call
C to DPIGMR (i.e., # iterations since the last restart) and
C the current order of the upper Hessenberg
C matrix HES.
C MAXL :IN Integer
C The maximum allowable order of the matrix H.
C MAXLP1 :IN Integer
C MAXPL1 = MAXL + 1, used for dynamic dimensioning of HES.
C V :IN Double Precision V(N,MAXLP1)
C The N by (LGMR+1) array containing the LGMR
C orthogonal vectors V(*,1) to V(*,LGMR).
C Q :IN Double Precision Q(2*MAXL)
C A double precision array of length 2*MAXL containing the
C components of the Givens rotations used in the QR
C decomposition of HES.
C SNORMW :IN Double Precision
C A scalar containing the scaled norm of VNEW before it
C is renormalized in DPIGMR.
C PROD :IN Double Precision
C The product s1*s2*...*sl = the product of the sines of the
C Givens rotations used in the QR factorization of the
C Hessenberg matrix HES.
C R0NRM :IN Double Precision
C The scaled norm of initial residual R0.
C HES :IN Double Precision HES(MAXLP1,MAXL)
C The upper triangular factor of the QR decomposition
C of the (LGMR+1) by LGMR upper Hessenberg matrix whose
C entries are the scaled inner-products of A*V(*,I)
C and V(*,K).
C JPRE :IN Integer
C Preconditioner type flag.
C (See description of IGWK(4) in DGMRES.)
C
C *Description
C When using the GMRES solver, the preferred value for ITOL
C is 0. This is due to the fact that when ITOL=0 the norm of
C the residual required in the stopping test is obtained for
C free, since this value is already calculated in the GMRES
C algorithm. The variable RNRM contains the appropriate
C norm, which is equal to norm(SB*(RL - A*XL)) when right or
C no preconditioning is being performed, and equal to
C norm(SB*Minv*(RL - A*XL)) when using left preconditioning.
C Here, norm() is the Euclidean norm. Nonzero values of ITOL
C require additional work to calculate the actual scaled
C residual or its scaled/preconditioned form, and/or the
C approximate solution XL. Hence, these values of ITOL will
C not be as efficient as ITOL=0.
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 This routine does not verify that ITOL has a valid value.
C The calling routine should make such a test before calling
C ISDGMR, as is done in DGMRES.
C
C***SEE ALSO DGMRES
C***ROUTINES CALLED D1MACH, DCOPY, DNRM2, DRLCAL, DSCAL, DXLCAL
C***COMMON BLOCKS DSLBLK
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 890922 Numerous changes to prologue to make closer to SLATEC
C standard. (FNF)
C 890929 Numerous changes to reduce SP/DP differences. (FNF)
C 910411 Prologue converted to Version 4.0 format. (BAB)
C 910502 Corrected conversion errors, etc. (FNF)
C 910502 Removed MSOLVE from ROUTINES CALLED list. (FNF)
C 910506 Made subsidiary to DGMRES. (FNF)
C 920407 COMMON BLOCK renamed DSLBLK. (WRB)
C 920511 Added complete declaration section. (WRB)
C 921026 Corrected D to E in output format. (FNF)
C 921113 Corrected C***CATEGORY line. (FNF)
C***END PROLOGUE ISDGMR