1      SUBROUTINE BSGQ8(FUN,XT,BC,N,KK,ID,A,B,INBV,ERR,ANS,IERR,WORK)
2C***BEGIN PROLOGUE  BSGQ8
3C***REFER TO  BFQAD
4C
5C     Written by R.E. Jones and modified by D.E. Amos
6C
7C     Abstract
8C        BSGQ8, a modification of GAUS8, integrates the
9C        product of FUN(X) by the ID-th derivative of a spline
10C        BVALU(XT,BC,N,KK,ID,X,INBV,WORK)  between limits A and B.
11C
12C        BSGQ8 calls BVALU, INTRV, I1MACH, R1MACH, XERROR
13C
14C     Description of Arguments
15C
16C        INPUT--
17C        FUN - Name of external function of one argument which
18C              multiplies BVALU.
19C        XT  - Knot array for BVALU
20C        BC  - B-coefficient array for BVALU
21C        N   - Number of B-coefficients for BVALU
22C        KK  - Order of the spline, KK.GE.1
23C        ID  - Order of the spline derivative, 0.LE.ID.LE.KK-1
24C        A   - Lower limit of integral
25C        B   - Upper limit of integral (may be less than A)
26C        INBV- Initialization parameter for BVALU
27C        ERR - Is a requested pseudorelative error tolerance.  Normally
28C              pick a value of ABS(ERR).LT.1E-3.  ANS will normally
29C              have no more error than ABS(ERR) times the integral of
30C              the absolute value of FUN(X)*BVALU(XT,BC,N,KK,X,ID,
31C              INBV,WORK).
32C
33C
34C        OUTPUT--
35C        ERR - Will be an estimate of the absolute error in ans if the
36C              input value of ERR was negative.  (ERR is unchanged if
37C              the input value of ERR was nonnegative.)  The estimated
38C              error is solely for information to the user and should
39C              not be used as a correction to the computed integral.
40C        ANS - Computed value of integral
41C        IERR- A status code
42C            --Normal Codes
43C               1 ANS most likely meets requested error tolerance,
44C                 or A=B.
45C              -1 A and B are too nearly equal to allow normal
46C                 integration.  ANS is set to zero.
47C            --Abnormal Code
48C               2 ANS probably does not meet requested error tolerance.
49C        WORK- Work vector of length 3*K for BVALU
50C***ROUTINES CALLED  BVALU,I1MACH,R1MACH,XERROR
51C***END PROLOGUE  BSGQ8
52C
53      INTEGER ICALL,ID,IERR,INBV,K, KK, KML, KMX, L, LMN, LMX, LR, MXL,
54     1 N, NBITS, NIB, NLMN, NLMX
55      INTEGER I1MACH
56      REAL A, AA, AE, ANIB, ANS, AREA, B, BC, C, CE, EE, EF, EPS, ERR,
57     1 EST,GL,GLR,GR,HH,SQ2,TOL,VL,VR,WORK,W1, W2, W3, W4, XT, X1,
58     2 X2, X3, X4, X, H
59      REAL R1MACH, BVALU, G8, FUN
60      DIMENSION XT(1), BC(1)
61      DIMENSION AA(30), HH(30), LR(30), VL(30), GR(30)
62      DATA X1, X2, X3, X4/
63     1     1.83434642495649805E-01,     5.25532409916328986E-01,
64     2     7.96666477413626740E-01,     9.60289856497536232E-01/
65      DATA W1, W2, W3, W4/
66     1     3.62683783378361983E-01,     3.13706645877887287E-01,
67     2     2.22381034453374471E-01,     1.01228536290376259E-01/
68      DATA ICALL  /  0  /
69      DATA SQ2/1.41421356E0/
70      DATA NLMN/1/,KMX/5000/,KML/6/
71      G8(X,H)=H*((W1*(FUN(X-X1*H)*BVALU(XT,BC,N,KK,ID,X-X1*H,INBV,WORK)+
72     1                FUN(X+X1*H)*BVALU(XT,BC,N,KK,ID,X+X1*H,INBV,WORK))
73     2           +W2*(FUN(X-X2*H)*BVALU(XT,BC,N,KK,ID,X-X2*H,INBV,WORK)+
74     3              FUN(X+X2*H)*BVALU(XT,BC,N,KK,ID,X+X2*H,INBV,WORK)))
75     4          +(W3*(FUN(X-X3*H)*BVALU(XT,BC,N,KK,ID,X-X3*H,INBV,WORK)+
76     5                FUN(X+X3*H)*BVALU(XT,BC,N,KK,ID,X+X3*H,INBV,WORK))
77     6           +W4*(FUN(X-X4*H)*BVALU(XT,BC,N,KK,ID,X-X4*H,INBV,WORK)+
78     7             FUN(X+X4*H)*BVALU(XT,BC,N,KK,ID,X+X4*H,INBV,WORK))))
79C
80C     INITIALIZE
81C
82C***FIRST EXECUTABLE STATEMENT  BSGQ8
83      IF (ICALL.NE.0) CALL XERROR( 'BSGQ8- BSGQ8 CALLED RECURSIVELY.  RE
84     1CURSIVE CALLS ARE ILLEGAL IN FORTRAN.', 73, 7, 2)
85      ICALL = 1
86      K = I1MACH(11)
87      ANIB = R1MACH(5)*FLOAT(K)/0.30102000E0
88      NBITS = INT(ANIB)
89      NLMX = (NBITS*5)/8
90      ANS = 0.0E0
91      IERR = 1
92      CE = 0.0E0
93      IF (A.EQ.B) GO TO 140
94      LMX = NLMX
95      LMN = NLMN
96      IF (B.EQ.0.0E0) GO TO 10
97      IF (SIGN(1.0E0,B)*A.LE.0.0E0) GO TO 10
98      C = ABS(1.0E0-A/B)
99      IF (C.GT.0.1E0) GO TO 10
100      IF (C.LE.0.0E0) GO TO 140
101      ANIB = 0.5E0 - ALOG(C)/0.69314718E0
102      NIB = INT(ANIB)
103      LMX = MIN0(NLMX,NBITS-NIB-7)
104      IF (LMX.LT.1) GO TO 130
105      LMN = MIN0(LMN,LMX)
106   10 TOL = AMAX1(ABS(ERR),2.0E0**(5-NBITS))/2.0E0
107      IF (ERR.EQ.0.0E0) TOL = SQRT(R1MACH(4))
108      EPS = TOL
109      HH(1) = (B-A)/4.0E0
110      AA(1) = A
111      LR(1) = 1
112      L = 1
113      EST = G8(AA(L)+2.0E0*HH(L),2.0E0*HH(L))
114      K = 8
115      AREA = ABS(EST)
116      EF = 0.5E0
117      MXL = 0
118C
119C     COMPUTE REFINED ESTIMATES, ESTIMATE THE ERROR, ETC.
120C
121   20 GL = G8(AA(L)+HH(L),HH(L))
122      GR(L) = G8(AA(L)+3.0E0*HH(L),HH(L))
123      K = K + 16
124      AREA = AREA + (ABS(GL)+ABS(GR(L))-ABS(EST))
125      GLR = GL + GR(L)
126      EE = ABS(EST-GLR)*EF
127      AE = AMAX1(EPS*AREA,TOL*ABS(GLR))
128      IF (EE-AE) 40, 40, 50
129   30 MXL = 1
130   40 CE = CE + (EST-GLR)
131      IF (LR(L)) 60, 60, 80
132C
133C     CONSIDER THE LEFT HALF OF THIS LEVEL
134C
135   50 IF (K.GT.KMX) LMX = KML
136      IF (L.GE.LMX) GO TO 30
137      L = L + 1
138      EPS = EPS*0.5E0
139      EF = EF/SQ2
140      HH(L) = HH(L-1)*0.5E0
141      LR(L) = -1
142      AA(L) = AA(L-1)
143      EST = GL
144      GO TO 20
145C
146C     PROCEED TO RIGHT HALF AT THIS LEVEL
147C
148   60 VL(L) = GLR
149   70 EST = GR(L-1)
150      LR(L) = 1
151      AA(L) = AA(L) + 4.0E0*HH(L)
152      GO TO 20
153C
154C     RETURN ONE LEVEL
155C
156   80 VR = GLR
157   90 IF (L.LE.1) GO TO 120
158      L = L - 1
159      EPS = EPS*2.0E0
160      EF = EF*SQ2
161      IF (LR(L)) 100, 100, 110
162  100 VL(L) = VL(L+1) + VR
163      GO TO 70
164  110 VR = VL(L+1) + VR
165      GO TO 90
166C
167C      EXIT
168C
169  120 ANS = VR
170      IF ((MXL.EQ.0) .OR. (ABS(CE).LE.2.0E0*TOL*AREA)) GO TO 140
171      IERR = 2
172      CALL XERROR( 'BSGQ8- ANS IS PROBABLY INSUFFICIENTLY ACCURATE.',
173     1 47, 3, 1)
174      GO TO 140
175  130 IERR = -1
176      CALL XERROR( 'BSGQ8- THE FOLLOWING TEMPORARY DIAGNOSTIC WILL APPEA
177     1R ONLY ONCE.  A AND B ARE TOO NEARLY EQUAL TO ALLOW NORMAL INTEGRA
178     2TION.  ANS IS SET TO ZERO, AND IERR=-1.', 157, 1, -1)
179  140 ICALL = 0
180      IF (ERR.LT.0.0E0) ERR = CE
181      RETURN
182      END
183