1 2 SUBROUTINE BFQAD(F,T,BCOEF,N,K,ID,X1,X2,TOL,QUAD,IERR,WORK) 3C***BEGIN PROLOGUE BFQAD 4C***DATE WRITTEN 800901 (YYMMDD) 5C***REVISION DATE 820801 (YYMMDD) 6C***CATEGORY NO. H2A2A1,E3,K6 7C***KEYWORDS B-SPLINE,DATA FITTING,INTERPOLATION,QUADRATURE,SPLINE 8C***AUTHOR AMOS, D. E., (SNLA) 9C***PURPOSE Computes the integral on (X1,X2) of a product of a function 10C F and the ID-th derivative of a K-th order B-spline 11C (B-representation). 12C***DESCRIPTION 13C 14C Written by D. E. Amos, June, 1979. 15C 16C Reference SAND 79-1825 17C 18C Abstract 19C BFQAD computes the integral on (X1,X2) of a product of a 20C function F and the ID-th derivative of a K-th order B-spline, 21C using the B-representation (T,BCOEF,N,K). (X1,X2) must be 22C a subinterval of T(K) .LE. X .le. T(N+1). An integration 23C routine BSGQ8 (a modification 24C of GAUS8), integrates the product on sub- 25C intervals of (X1,X2) formed by included (distinct) knots. 26C 27C BFQAD calls INTRV,BVALU,BSGQ8,R1MACH,XERROR 28C 29C Description of Arguments 30C Input 31C F - external function of one argument for the 32C integrand BF(X)=F(X)*BVALU(T,BCOEF,N,K,ID,X,INBV, 33C WORK) 34C T - knot array of length N+K 35C BCOEF - coefficient array of length N 36C N - length of coefficient array 37C K - order of B-spline, K .GE. 1 38C ID - order of the spline derivative, 0 .LE. ID .LE. K-1 39C ID=0 gives the spline function 40C X1,X2 - end points of quadrature interval in 41C T(K) .LE. X .LE. T(N+1) 42C TOL - desired accuracy for the quadrature, suggest 43C 10.*STOL .LT. TOL .LE. 0.1 where STOL is the single 44C precision unit roundoff for the machine = R1MACH(4) 45C 46C Output 47C QUAD - integral of BF(X) on (X1,X2) 48C IERR - a status code 49C IERR=1 normal return 50C 2 some quadrature on (X1,X2) does not meet 51C the requested tolerance. 52C WORK - work vector of length 3*K 53C 54C Error Conditions 55C X1 or X2 not in T(K) .LE. X .LE. T(N+1) is a fatal error. 56C TOL not greater than the single precision unit roundoff or 57C less than 0.1 is a fatal error. 58C Some quadrature fails to meet the requested tolerance. 59C***REFERENCES D.E. AMOS, *QUADRATURE SUBROUTINES FOR SPLINES AND 60C B-SPLINES*, SAND79-1825, SANDIA LABORATORIES, 61C DECEMBER 1979. 62C***ROUTINES CALLED BSGQ8,INTRV,R1MACH,XERROR 63C***END PROLOGUE BFQAD 64C 65C 66 INTEGER ID, IERR, IFLG, ILO, IL1, IL2, K, LEFT, MFLAG, N, NPK, NP1 67 REAL A,AA,ANS,B,BB,BCOEF,Q,QUAD,T,TA,TB,TOL,WORK,WTOL, X1, 68 1 X2 69 REAL R1MACH, F 70 DIMENSION T(1), BCOEF(1), WORK(1) 71 EXTERNAL F 72C***FIRST EXECUTABLE STATEMENT BFQAD 73 IERR = 1 74 QUAD = 0.0E0 75 IF(K.LT.1) GO TO 100 76 IF(N.LT.K) GO TO 105 77 IF(ID.LT.0 .OR. ID.GE.K) GO TO 110 78 WTOL = R1MACH(4) 79 IF (TOL.LT.WTOL .OR. TOL.GT.0.1E0) GO TO 30 80 AA = AMIN1(X1,X2) 81 BB = AMAX1(X1,X2) 82 IF (AA.LT.T(K)) GO TO 20 83 NP1 = N + 1 84 IF (BB.GT.T(NP1)) GO TO 20 85 IF (AA.EQ.BB) RETURN 86 NPK = N + K 87C 88 ILO = 1 89 CALL INTRV(T, NPK, AA, ILO, IL1, MFLAG) 90 CALL INTRV(T, NPK, BB, ILO, IL2, MFLAG) 91 IF (IL2.GE.NP1) IL2 = N 92 INBV = 1 93 Q = 0.0E0 94 DO 10 LEFT=IL1,IL2 95 TA = T(LEFT) 96 TB = T(LEFT+1) 97 IF (TA.EQ.TB) GO TO 10 98 A = AMAX1(AA,TA) 99 B = AMIN1(BB,TB) 100 CALL BSGQ8(F,T,BCOEF,N,K,ID,A,B,INBV,TOL,ANS,IFLG,WORK) 101 IF (IFLG.GT.1) IERR = 2 102 Q = Q + ANS 103 10 CONTINUE 104 IF (X1.GT.X2) Q = -Q 105 QUAD = Q 106 RETURN 107C 108C 109 20 CONTINUE 110 CALL XERROR( ' BFQAD, X1 OR X2 OR BOTH DO NOT SATISFY T(K).LE.X.L 111 1E.T(N+1)', 60, 2, 1) 112 RETURN 113 30 CONTINUE 114 CALL XERROR( ' BFQAD, TOL IS LESS THAN THE SINGLE PRECISION TOLER 115 1ANCE OR GREATER THAN 0.1', 76, 2, 1) 116 RETURN 117 100 CONTINUE 118 CALL XERROR( ' BFQAD, K DOES NOT SATISFY K.GE.1', 34, 2, 1) 119 RETURN 120 105 CONTINUE 121 CALL XERROR( ' BFQAD, N DOES NOT SATISFY N.GE.K', 34, 2, 1) 122 RETURN 123 110 CONTINUE 124 CALL XERROR( ' BFQAD, ID DOES NOT SATISFY 0.LE.ID.LT.K', 125 1 41, 2, 1) 126 RETURN 127 END 128