1*DECK ISSBCG 2 INTEGER FUNCTION ISSBCG (N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, 3 + ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, RR, ZZ, PP, 4 + DZ, RWORK, IWORK, AK, BK, BNRM, SOLNRM) 5C***BEGIN PROLOGUE ISSBCG 6C***SUBSIDIARY 7C***PURPOSE Preconditioned BiConjugate Gradient Stop Test. 8C This routine calculates the stop test for the BiConjugate 9C Gradient iteration scheme. It returns a non-zero if the 10C error estimate (the type of which is determined by ITOL) 11C is less than the user specified tolerance TOL. 12C***LIBRARY SLATEC (SLAP) 13C***CATEGORY D2A4, D2B4 14C***TYPE SINGLE PRECISION (ISSBCG-S, ISDBCG-D) 15C***KEYWORDS ITERATIVE PRECONDITION, NON-SYMMETRIC LINEAR SYSTEM, SLAP, 16C SPARSE, STOP TEST 17C***AUTHOR Greenbaum, Anne, (Courant Institute) 18C Seager, Mark K., (LLNL) 19C Lawrence Livermore National Laboratory 20C PO BOX 808, L-60 21C Livermore, CA 94550 (510) 423-3141 22C seager@llnl.gov 23C***DESCRIPTION 24C 25C *Usage: 26C INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX, ITER 27C INTEGER IERR, IUNIT, IWORK(USER DEFINED) 28C REAL B(N), X(N), A(N), TOL, ERR, R(N), Z(N), P(N) 29C REAL RR(N), ZZ(N), PP(N), DZ(N) 30C REAL RWORK(USER DEFINED), AK, BK, BNRM, SOLNRM 31C EXTERNAL MSOLVE 32C 33C IF( ISSBCG(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL, 34C $ ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, RR, ZZ, PP, DZ, 35C $ RWORK, IWORK, AK, BK, BNRM, SOLNRM) .NE. 0 ) 36C $ THEN ITERATION DONE 37C 38C *Arguments: 39C N :IN Integer 40C Order of the Matrix. 41C B :IN Real B(N). 42C Right-hand side vector. 43C X :INOUT Real X(N). 44C On input X is your initial guess for solution vector. 45C On output X is the final approximate solution. 46C NELT :IN Integer. 47C Number of Non-Zeros stored in A. 48C IA :IN Integer IA(NELT). 49C JA :IN Integer JA(NELT). 50C A :IN Real A(NELT). 51C These arrays contain the matrix data structure for A. 52C It could take any form. See "Description", in the SLAP 53C routine SBCG for more details. 54C ISYM :IN Integer. 55C Flag to indicate symmetric storage format. 56C If ISYM=0, all non-zero entries of the matrix are stored. 57C If ISYM=1, the matrix is symmetric, and only the upper 58C or lower triangle of the matrix is stored. 59C MSOLVE :EXT External. 60C Name of a routine which solves a linear system MZ = R for Z 61C given R with the preconditioning matrix M (M is supplied via 62C RWORK and IWORK arrays). The name of the MSOLVE routine 63C must be declared external in the calling program. The 64C calling sequence of MSOLVE is: 65C CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) 66C Where N is the number of unknowns, R is the right-hand side 67C vector, and Z is the solution upon return. NELT, IA, JA, A, 68C and ISYM define the SLAP matrix data structure. 69C RWORK is a real array that can be used to pass necessary 70C preconditioning information and/or workspace to MSOLVE. 71C IWORK is an integer work array for the same purpose as RWORK. 72C ITOL :IN Integer. 73C Flag to indicate type of convergence criterion. 74C If ITOL=1, iteration stops when the 2-norm of the residual 75C divided by the 2-norm of the right-hand side is less than TOL. 76C If ITOL=2, iteration stops when the 2-norm of M-inv times the 77C residual divided by the 2-norm of M-inv times the right hand 78C side is less than TOL, where M-inv is the inverse of the 79C diagonal of A. 80C ITOL=11 is often useful for checking and comparing different 81C routines. For this case, the user must supply the "exact" 82C solution or a very accurate approximation (one with an error 83C much less than TOL) through a common block, 84C COMMON /SSLBLK/ SOLN( ) 85C If ITOL=11, iteration stops when the 2-norm of the difference 86C between the iterative approximation and the user-supplied 87C solution divided by the 2-norm of the user-supplied solution 88C is less than TOL. Note that this requires the user to set up 89C the "COMMON /SSLBLK/ SOLN(LENGTH)" in the calling routine. 90C The routine with this declaration should be loaded before the 91C stop test so that the correct length is used by the loader. 92C This procedure is not standard Fortran and may not work 93C correctly on your system (although it has worked on every 94C system the authors have tried). If ITOL is not 11 then this 95C common block is indeed standard Fortran. 96C TOL :IN Real. 97C Convergence criterion, as described above. 98C ITMAX :IN Integer. 99C Maximum number of iterations. 100C ITER :IN Integer. 101C Current iteration count. (Must be zero on first call.) 102C ERR :OUT Real. 103C Error estimate of error in final approximate solution, as 104C defined by ITOL. 105C IERR :OUT Integer. 106C Error flag. IERR is set to 3 if ITOL is not one of the 107C acceptable values, see above. 108C IUNIT :IN Integer. 109C Unit number on which to write the error at each iteration, 110C if this is desired for monitoring convergence. If unit 111C number is 0, no writing will occur. 112C R :IN Real R(N). 113C The residual r = b - Ax. 114C Z :WORK Real Z(N). 115C P :DUMMY Real P(N). 116C RR :DUMMY Real RR(N). 117C ZZ :DUMMY Real ZZ(N). 118C PP :DUMMY Real PP(N). 119C Real arrays used for workspace. 120C DZ :WORK Real DZ(N). 121C If ITOL.eq.0 then DZ is used to hold M-inv * B on the first 122C call. If ITOL.eq.11 then DZ is used to hold X-SOLN. 123C RWORK :WORK Real RWORK(USER DEFINED). 124C Real array that can be used for workspace in MSOLVE 125C and MTSOLV. 126C IWORK :WORK Integer IWORK(USER DEFINED). 127C Integer array that can be used for workspace in MSOLVE 128C and MTSOLV. 129C AK :IN Real. 130C Current iterate BiConjugate Gradient iteration parameter. 131C BK :IN Real. 132C Current iterate BiConjugate Gradient iteration parameter. 133C BNRM :INOUT Real. 134C Norm of the right hand side. Type of norm depends on ITOL. 135C Calculated only on the first call. 136C SOLNRM :INOUT Real. 137C 2-Norm of the true solution, SOLN. Only computed and used 138C if ITOL = 11. 139C 140C *Function Return Values: 141C 0 : Error estimate (determined by ITOL) is *NOT* less than the 142C specified tolerance, TOL. The iteration must continue. 143C 1 : Error estimate (determined by ITOL) is less than the 144C specified tolerance, TOL. The iteration can be considered 145C complete. 146C 147C *Cautions: 148C This routine will attempt to write to the Fortran logical output 149C unit IUNIT, if IUNIT .ne. 0. Thus, the user must make sure that 150C this logical unit is attached to a file or terminal before calling 151C this routine with a non-zero value for IUNIT. This routine does 152C not check for the validity of a non-zero IUNIT unit number. 153C 154C***SEE ALSO SBCG 155C***ROUTINES CALLED R1MACH, SNRM2 156C***COMMON BLOCKS SSLBLK 157C***REVISION HISTORY (YYMMDD) 158C 871119 DATE WRITTEN 159C 881213 Previous REVISION DATE 160C 890915 Made changes requested at July 1989 CML Meeting. (MKS) 161C 890922 Numerous changes to prologue to make closer to SLATEC 162C standard. (FNF) 163C 890929 Numerous changes to reduce SP/DP differences. (FNF) 164C 891003 Removed C***REFER TO line, per MKS. 165C 910411 Prologue converted to Version 4.0 format. (BAB) 166C 910502 Removed MSOLVE from ROUTINES CALLED list. (FNF) 167C 910506 Made subsidiary to SBCG. (FNF) 168C 920407 COMMON BLOCK renamed SSLBLK. (WRB) 169C 920511 Added complete declaration section. (WRB) 170C 920930 Corrected to not print AK,BK when ITER=0. (FNF) 171C 921026 Changed 1.0E10 to R1MACH(2). (FNF) 172C 921113 Corrected C***CATEGORY line. (FNF) 173C***END PROLOGUE ISSBCG 174C .. Scalar Arguments .. 175 REAL AK, BK, BNRM, ERR, SOLNRM, TOL 176 INTEGER IERR, ISYM, ITER, ITMAX, ITOL, IUNIT, N, NELT 177C .. Array Arguments .. 178 REAL A(NELT), B(N), DZ(N), P(N), PP(N), R(N), RR(N), RWORK(*), 179 + X(N), Z(N), ZZ(N) 180 INTEGER IA(NELT), IWORK(*), JA(NELT) 181C .. Subroutine Arguments .. 182 EXTERNAL MSOLVE 183C .. Arrays in Common .. 184 REAL SOLN(1) 185C .. Local Scalars .. 186 INTEGER I 187C .. External Functions .. 188 REAL R1MACH, SNRM2 189 EXTERNAL R1MACH, SNRM2 190C .. Common blocks .. 191 COMMON /SSLBLK/ SOLN 192C***FIRST EXECUTABLE STATEMENT ISSBCG 193 ISSBCG = 0 194C 195 IF( ITOL.EQ.1 ) THEN 196C err = ||Residual||/||RightHandSide|| (2-Norms). 197 IF(ITER .EQ. 0) BNRM = SNRM2(N, B, 1) 198 ERR = SNRM2(N, R, 1)/BNRM 199 ELSE IF( ITOL.EQ.2 ) THEN 200C -1 -1 201C err = ||M Residual||/||M RightHandSide|| (2-Norms). 202 IF(ITER .EQ. 0) THEN 203 CALL MSOLVE(N, B, DZ, NELT, IA, JA, A, ISYM, RWORK, IWORK) 204 BNRM = SNRM2(N, DZ, 1) 205 ENDIF 206 ERR = SNRM2(N, Z, 1)/BNRM 207 ELSE IF( ITOL.EQ.11 ) THEN 208C err = ||x-TrueSolution||/||TrueSolution|| (2-Norms). 209 IF(ITER .EQ. 0) SOLNRM = SNRM2(N, SOLN, 1) 210 DO 10 I = 1, N 211 DZ(I) = X(I) - SOLN(I) 212 10 CONTINUE 213 ERR = SNRM2(N, DZ, 1)/SOLNRM 214 ELSE 215C 216C If we get here ITOL is not one of the acceptable values. 217 ERR = R1MACH(2) 218 IERR = 3 219 ENDIF 220C 221 IF(IUNIT .NE. 0) THEN 222 IF( ITER.EQ.0 ) THEN 223 WRITE(IUNIT,1000) N, ITOL 224 WRITE(IUNIT,1010) ITER, ERR 225 ELSE 226 WRITE(IUNIT,1010) ITER, ERR, AK, BK 227 ENDIF 228 ENDIF 229 IF(ERR .LE. TOL) ISSBCG = 1 230C 231 RETURN 232 1000 FORMAT(' Preconditioned BiConjugate Gradient for N, ITOL = ', 233 $ I5,I5,/' ITER',' Error Estimate',' Alpha', 234 $ ' Beta') 235 1010 FORMAT(1X,I4,1X,E16.7,1X,E16.7,1X,E16.7) 236C------------- LAST LINE OF ISSBCG FOLLOWS ---------------------------- 237 END 238