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