*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
C .. Scalar Arguments ..
DOUBLE PRECISION BNRM, ERR, PROD, R0NRM, RNRM, SNORMW, TOL
INTEGER ISYM, ITER, ITMAX, ITOL, IUNIT, JPRE, JSCAL, KMP, LGMR,
+ MAXL, MAXLP1, N, NELT, NMSL
C .. Array Arguments ..
DOUBLE PRECISION A(*), B(*), DZ(*), HES(MAXLP1, MAXL), Q(*), R(*),
+ RWORK(*), SB(*), SX(*), V(N,*), X(*), XL(*), Z(*)
INTEGER IA(*), IWORK(*), JA(*)
C .. Subroutine Arguments ..
EXTERNAL MSOLVE
C .. Arrays in Common ..
DOUBLE PRECISION SOLN(1)
C .. Local Scalars ..
DOUBLE PRECISION DXNRM, FUZZ, RAT, RATMAX, SOLNRM, TEM
INTEGER I, IELMAX
C .. External Functions ..
DOUBLE PRECISION D1MACH, DNRM2
EXTERNAL D1MACH, DNRM2
C .. External Subroutines ..
EXTERNAL DCOPY, DRLCAL, DSCAL, DXLCAL
C .. Intrinsic Functions ..
INTRINSIC ABS, MAX, SQRT
C .. Common blocks ..
COMMON /DSLBLK/ SOLN
C .. Save statement ..
SAVE SOLNRM
C***FIRST EXECUTABLE STATEMENT ISDGMR
ISDGMR = 0
IF ( ITOL.EQ.0 ) THEN
C
C Use input from DPIGMR to determine if stop conditions are met.
C
ERR = RNRM/BNRM
ENDIF
IF ( (ITOL.GT.0) .AND. (ITOL.LE.3) ) THEN
C
C Use DRLCAL to calculate the scaled residual vector.
C Store answer in R.
C
IF ( LGMR.NE.0 ) CALL DRLCAL(N, KMP, LGMR, MAXL, V, Q, R,
$ SNORMW, PROD, R0NRM)
IF ( ITOL.LE.2 ) THEN
C err = ||Residual||/||RightHandSide||(2-Norms).
ERR = DNRM2(N, R, 1)/BNRM
C
C Unscale R by R0NRM*PROD when KMP < MAXL.
C
IF ( (KMP.LT.MAXL) .AND. (LGMR.NE.0) ) THEN
TEM = 1.0D0/(R0NRM*PROD)
CALL DSCAL(N, TEM, R, 1)
ENDIF
ELSEIF ( ITOL.EQ.3 ) THEN
C err = Max |(Minv*Residual)(i)/x(i)|
C When JPRE .lt. 0, R already contains Minv*Residual.
IF ( JPRE.GT.0 ) THEN
CALL MSOLVE(N, R, DZ, NELT, IA, JA, A, ISYM, RWORK,
$ IWORK)
NMSL = NMSL + 1
ENDIF
C
C Unscale R by R0NRM*PROD when KMP < MAXL.
C
IF ( (KMP.LT.MAXL) .AND. (LGMR.NE.0) ) THEN
TEM = 1.0D0/(R0NRM*PROD)
CALL DSCAL(N, TEM, R, 1)
ENDIF
C
FUZZ = D1MACH(1)
IELMAX = 1
RATMAX = ABS(DZ(1))/MAX(ABS(X(1)),FUZZ)
DO 25 I = 2, N
RAT = ABS(DZ(I))/MAX(ABS(X(I)),FUZZ)
IF( RAT.GT.RATMAX ) THEN
IELMAX = I
RATMAX = RAT
ENDIF
25 CONTINUE
ERR = RATMAX
IF( RATMAX.LE.TOL ) ISDGMR = 1
IF( IUNIT.GT.0 ) WRITE(IUNIT,1020) ITER, IELMAX, RATMAX
RETURN
ENDIF
ENDIF
IF ( ITOL.EQ.11 ) THEN
C
C Use DXLCAL to calculate the approximate solution XL.
C
IF ( (LGMR.NE.0) .AND. (ITER.GT.0) ) THEN
CALL DXLCAL(N, LGMR, X, XL, XL, HES, MAXLP1, Q, V, R0NRM,
$ DZ, SX, JSCAL, JPRE, MSOLVE, NMSL, RWORK, IWORK,
$ NELT, IA, JA, A, ISYM)
ELSEIF ( ITER.EQ.0 ) THEN
C Copy X to XL to check if initial guess is good enough.
CALL DCOPY(N, X, 1, XL, 1)
ELSE
C Return since this is the first call to DPIGMR on a restart.
RETURN
ENDIF
C
IF ((JSCAL .EQ. 0) .OR.(JSCAL .EQ. 2)) THEN
C err = ||x-TrueSolution||/||TrueSolution||(2-Norms).
IF ( ITER.EQ.0 ) SOLNRM = DNRM2(N, SOLN, 1)
DO 30 I = 1, N
DZ(I) = XL(I) - SOLN(I)
30 CONTINUE
ERR = DNRM2(N, DZ, 1)/SOLNRM
ELSE
IF (ITER .EQ. 0) THEN
SOLNRM = 0
DO 40 I = 1,N
SOLNRM = SOLNRM + (SX(I)*SOLN(I))**2
40 CONTINUE
SOLNRM = SQRT(SOLNRM)
ENDIF
DXNRM = 0
DO 50 I = 1,N
DXNRM = DXNRM + (SX(I)*(XL(I)-SOLN(I)))**2
50 CONTINUE
DXNRM = SQRT(DXNRM)
C err = ||SX*(x-TrueSolution)||/||SX*TrueSolution|| (2-Norms).
ERR = DXNRM/SOLNRM
ENDIF
ENDIF
C
IF( IUNIT.NE.0 ) THEN
IF( ITER.EQ.0 ) THEN
WRITE(IUNIT,1000) N, ITOL, MAXL, KMP
ENDIF
WRITE(IUNIT,1010) ITER, RNRM/BNRM, ERR
ENDIF
IF ( ERR.LE.TOL ) ISDGMR = 1
C
RETURN
1000 FORMAT(' Generalized Minimum Residual(',I3,I3,') for ',
$ 'N, ITOL = ',I5, I5,
$ /' ITER',' Natural Err Est',' Error Estimate')
1010 FORMAT(1X,I4,1X,D16.7,1X,D16.7)
1020 FORMAT(1X,' ITER = ',I5, ' IELMAX = ',I5,
$ ' |R(IELMAX)/X(IELMAX)| = ',D12.5)
C------------- LAST LINE OF ISDGMR FOLLOWS ----------------------------
END