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