1 SUBROUTINE PPQAD(LDC,C,XI,LXI,K,X1,X2,PQUAD) 2C***BEGIN PROLOGUE PPQAD 3C***DATE WRITTEN 800901 (YYMMDD) 4C***REVISION DATE 820801 (YYMMDD) 5C***CATEGORY NO. H2A2A1,E3,K6 6C***KEYWORDS B-SPLINE,DATA FITTING,INTERPOLATION,QUADRATURE,SPLINE 7C***AUTHOR AMOS, D. E., (SNLA) 8C***PURPOSE Computes the integral on (X1,X2) of a K-th order B-spline 9C using the piecewise polynomial representation. 10C***DESCRIPTION 11C 12C Written by D. E. Amos, June, 1979. 13C 14C Reference SAND-79-1825 15C 16C Abstract 17C PPQAD computes the integral on (X1,X2) of a K-th order 18C B-spline using the piecewise polynomial representation 19C (C,XI,LXI,K). Here the Taylor expansion about the left 20C end point XI(J) of the J-th interval is integrated and 21C evaluated on subintervals of (X1,X2) which are formed by 22C included break points. Integration outside (XI(1),XI(LXI+1)) 23C is permitted. 24C 25C PPQAD calls INTRV 26C 27C Description of Arguments 28C Input 29C LDC - leading dimension of matrix C, LDC .GE. K 30C C(I,J) - right Taylor derivatives at XI(J), I=1,K , J=1,LXI 31C XI(*) - break point array of length LXI+1 32C LXI - number of polynomial pieces 33C K - order of B-spline, K .GE. 1 34C X1,X2 - end points of quadrature interval, normally in 35C XI(1) .LE. X .LE. XI(LXI+1) 36C 37C Output 38C PQUAD - integral of the PP representation over (X1,X2) 39C 40C Error Conditions 41C Improper input is a fatal error 42C***REFERENCES D.E. AMOS, *QUADRATURE SUBROUTINES FOR SPLINES AND 43C B-SPLINES*, SAND79-1825, SANDIA LABORATORIES, 44C DECEMBER 1979. 45C***ROUTINES CALLED INTRV,XERROR 46C***END PROLOGUE PPQAD 47C 48C 49 INTEGER I, II, IL, ILO, IL1, IL2, IM, K, LDC, LEFT, LXI, MF1, MF2 50 REAL A, AA, BB, C, DX, FLK, PQUAD, Q, S, SS, TA, TB, X, XI, X1, X2 51 DIMENSION XI(1), C(LDC,1), SS(2) 52C 53C***FIRST EXECUTABLE STATEMENT PPQAD 54 PQUAD = 0.0E0 55 IF(K.LT.1) GO TO 100 56 IF(LXI.LT.1) GO TO 105 57 IF(LDC.LT.K) GO TO 110 58 AA = AMIN1(X1,X2) 59 BB = AMAX1(X1,X2) 60 IF (AA.EQ.BB) RETURN 61 ILO = 1 62 CALL INTRV(XI, LXI, AA, ILO, IL1, MF1) 63 CALL INTRV(XI, LXI, BB, ILO, IL2, MF2) 64 Q = 0.0E0 65 DO 40 LEFT=IL1,IL2 66 TA = XI(LEFT) 67 A = AMAX1(AA,TA) 68 IF (LEFT.EQ.1) A = AA 69 TB = BB 70 IF (LEFT.LT.LXI) TB = XI(LEFT+1) 71 X = AMIN1(BB,TB) 72 DO 30 II=1,2 73 SS(II) = 0.0E0 74 DX = X - XI(LEFT) 75 IF (DX.EQ.0.0E0) GO TO 20 76 S = C(K,LEFT) 77 FLK = FLOAT(K) 78 IM = K - 1 79 IL = IM 80 DO 10 I=1,IL 81 S = S*DX/FLK + C(IM,LEFT) 82 IM = IM - 1 83 FLK = FLK - 1.0E0 84 10 CONTINUE 85 SS(II) = S*DX 86 20 CONTINUE 87 X = A 88 30 CONTINUE 89 Q = Q + (SS(1)-SS(2)) 90 40 CONTINUE 91 IF (X1.GT.X2) Q = -Q 92 PQUAD = Q 93 RETURN 94C 95C 96 100 CONTINUE 97 CALL XERROR( ' PPQAD, K DOES NOT SATISFY K.GE.1', 34, 2, 1) 98 RETURN 99 105 CONTINUE 100 CALL XERROR( ' PPQAD, LXI DOES NOT SATISFY LXI.GE.1', 38, 2, 1) 101 RETURN 102 110 CONTINUE 103 CALL XERROR( ' PPQAD, LDC DOES NOT SATISFY LDC.GE.K', 38, 2, 1) 104 RETURN 105 END 106