*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