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