1*DECK GAUS8 2 SUBROUTINE GAUS8 (FUN, A, B, ERR, ANS, IERR) 3C***BEGIN PROLOGUE GAUS8 4C***PURPOSE Integrate a real function of one variable over a finite 5C interval using an adaptive 8-point Legendre-Gauss 6C algorithm. Intended primarily for high accuracy 7C integration or integration of smooth functions. 8C***LIBRARY SLATEC 9C***CATEGORY H2A1A1 10C***TYPE SINGLE PRECISION (GAUS8-S, DGAUS8-D) 11C***KEYWORDS ADAPTIVE QUADRATURE, AUTOMATIC INTEGRATOR, 12C GAUSS QUADRATURE, NUMERICAL INTEGRATION 13C***AUTHOR Jones, R. E., (SNLA) 14C***DESCRIPTION 15C 16C Abstract 17C GAUS8 integrates real functions of one variable over finite 18C intervals using an adaptive 8-point Legendre-Gauss algorithm. 19C GAUS8 is intended primarily for high accuracy integration 20C or integration of smooth functions. 21C 22C Description of Arguments 23C 24C Input-- 25C FUN - name of external function to be integrated. This name 26C must be in an EXTERNAL statement in the calling program. 27C FUN must be a REAL function of one REAL argument. The 28C value of the argument to FUN is the variable of 29C integration which ranges from A to B. 30C A - lower limit of integration 31C B - upper limit of integration (may be less than A) 32C ERR - is a requested pseudorelative error tolerance. Normally 33C pick a value of ABS(ERR) so that STOL .LT. ABS(ERR) .LE. 34C 1.0E-3 where STOL is the single precision unit roundoff 35C R1MACH(4). ANS will normally have no more error than 36C ABS(ERR) times the integral of the absolute value of 37C FUN(X). Usually, smaller values for ERR yield more 38C accuracy and require more function evaluations. 39C 40C A negative value for ERR causes an estimate of the 41C absolute error in ANS to be returned in ERR. Note that 42C ERR must be a variable (not a constant) in this case. 43C Note also that the user must reset the value of ERR 44C before making any more calls that use the variable ERR. 45C 46C Output-- 47C ERR - will be an estimate of the absolute error in ANS if the 48C input value of ERR was negative. (ERR is unchanged if 49C the input value of ERR was non-negative.) The estimated 50C error is solely for information to the user and should 51C not be used as a correction to the computed integral. 52C ANS - computed value of integral 53C IERR- a status code 54C --Normal codes 55C 1 ANS most likely meets requested error tolerance, 56C or A=B. 57C -1 A and B are too nearly equal to allow normal 58C integration. ANS is set to zero. 59C --Abnormal code 60C 2 ANS probably does not meet requested error tolerance. 61C 62C***REFERENCES (NONE) 63C***ROUTINES CALLED I1MACH, R1MACH, XERMSG 64C***REVISION HISTORY (YYMMDD) 65C 810223 DATE WRITTEN 66C 890531 Changed all specific intrinsics to generic. (WRB) 67C 890531 REVISION DATE from Version 3.2 68C 891214 Prologue converted to Version 4.0 format. (BAB) 69C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 70C 900326 Removed duplicate information from DESCRIPTION section. 71C (WRB) 72C***END PROLOGUE GAUS8 73 INTEGER IERR, K, KML, KMX, L, LMN, LMX, LR, MXL, NBITS, 74 1 NIB, NLMN, NLMX 75 INTEGER I1MACH 76 REAL A, AA, AE, ANIB, ANS, AREA, B, C, CE, EE, EF, EPS, ERR, EST, 77 1 GL, GLR, GR, HH, SQ2, TOL, VL, VR, W1, W2, W3, W4, X1, X2, X3, 78 2 X4, X, H 79 REAL R1MACH, G8, FUN 80 DIMENSION AA(30), HH(30), LR(30), VL(30), GR(30) 81 SAVE X1, X2, X3, X4, W1, W2, W3, W4, SQ2, 82 1 NLMN, KMX, KML 83 DATA X1, X2, X3, X4/ 84 1 1.83434642495649805E-01, 5.25532409916328986E-01, 85 2 7.96666477413626740E-01, 9.60289856497536232E-01/ 86 DATA W1, W2, W3, W4/ 87 1 3.62683783378361983E-01, 3.13706645877887287E-01, 88 2 2.22381034453374471E-01, 1.01228536290376259E-01/ 89 DATA SQ2/1.41421356E0/ 90 DATA NLMN/1/,KMX/5000/,KML/6/ 91 G8(X,H)=H*((W1*(FUN(X-X1*H) + FUN(X+X1*H)) 92 1 +W2*(FUN(X-X2*H) + FUN(X+X2*H))) 93 2 +(W3*(FUN(X-X3*H) + FUN(X+X3*H)) 94 3 +W4*(FUN(X-X4*H) + FUN(X+X4*H)))) 95C***FIRST EXECUTABLE STATEMENT GAUS8 96C 97C Initialize 98C 99 K = I1MACH(11) 100 ANIB = R1MACH(5)*K/0.30102000E0 101 NBITS = ANIB 102 NLMX = MIN(30,(NBITS*5)/8) 103 ANS = 0.0E0 104 IERR = 1 105 CE = 0.0E0 106 IF (A .EQ. B) GO TO 140 107 LMX = NLMX 108 LMN = NLMN 109 IF (B .EQ. 0.0E0) GO TO 10 110 IF (SIGN(1.0E0,B)*A .LE. 0.0E0) GO TO 10 111 C = ABS(1.0E0-A/B) 112 IF (C .GT. 0.1E0) GO TO 10 113 IF (C .LE. 0.0E0) GO TO 140 114 ANIB = 0.5E0 - LOG(C)/0.69314718E0 115 NIB = ANIB 116 LMX = MIN(NLMX,NBITS-NIB-7) 117 IF (LMX .LT. 1) GO TO 130 118 LMN = MIN(LMN,LMX) 119 10 TOL = MAX(ABS(ERR),2.0E0**(5-NBITS))/2.0E0 120 IF (ERR .EQ. 0.0E0) TOL = SQRT(R1MACH(4)) 121 EPS = TOL 122 HH(1) = (B-A)/4.0E0 123 AA(1) = A 124 LR(1) = 1 125 L = 1 126 EST = G8(AA(L)+2.0E0*HH(L),2.0E0*HH(L)) 127 K = 8 128 AREA = ABS(EST) 129 EF = 0.5E0 130 MXL = 0 131C 132C Compute refined estimates, estimate the error, etc. 133C 134 20 GL = G8(AA(L)+HH(L),HH(L)) 135 GR(L) = G8(AA(L)+3.0E0*HH(L),HH(L)) 136 K = K + 16 137 AREA = AREA + (ABS(GL)+ABS(GR(L))-ABS(EST)) 138C IF (L .LT. LMN) GO TO 11 139 GLR = GL + GR(L) 140 EE = ABS(EST-GLR)*EF 141 AE = MAX(EPS*AREA,TOL*ABS(GLR)) 142 IF (EE-AE) 40, 40, 50 143 30 MXL = 1 144 40 CE = CE + (EST-GLR) 145 IF (LR(L)) 60, 60, 80 146C 147C Consider the left half of this level 148C 149 50 IF (K .GT. KMX) LMX = KML 150 IF (L .GE. LMX) GO TO 30 151 L = L + 1 152 EPS = EPS*0.5E0 153 EF = EF/SQ2 154 HH(L) = HH(L-1)*0.5E0 155 LR(L) = -1 156 AA(L) = AA(L-1) 157 EST = GL 158 GO TO 20 159C 160C Proceed to right half at this level 161C 162 60 VL(L) = GLR 163 70 EST = GR(L-1) 164 LR(L) = 1 165 AA(L) = AA(L) + 4.0E0*HH(L) 166 GO TO 20 167C 168C Return one level 169C 170 80 VR = GLR 171 90 IF (L .LE. 1) GO TO 120 172 L = L - 1 173 EPS = EPS*2.0E0 174 EF = EF*SQ2 175 IF (LR(L)) 100, 100, 110 176 100 VL(L) = VL(L+1) + VR 177 GO TO 70 178 110 VR = VL(L+1) + VR 179 GO TO 90 180C 181C Exit 182C 183 120 ANS = VR 184 IF ((MXL.EQ.0) .OR. (ABS(CE).LE.2.0E0*TOL*AREA)) GO TO 140 185 IERR = 2 186 CALL XERMSG ('SLATEC', 'GAUS8', 187 + 'ANS is probably insufficiently accurate.', 3, 1) 188 GO TO 140 189 130 IERR = -1 190 CALL XERMSG ('SLATEC', 'GAUS8', 191 + 'A and B are too nearly equal to allow normal integration. $$' 192 + // 'ANS is set to zero and IERR to -1.', 1, -1) 193 140 IF (ERR .LT. 0.0E0) ERR = CE 194 RETURN 195 END 196