*DECK DPPQAD SUBROUTINE DPPQAD (LDC, C, XI, LXI, K, X1, X2, PQUAD) C***BEGIN PROLOGUE DPPQAD C***PURPOSE Compute the integral on (X1,X2) of a K-th order B-spline C using the piecewise polynomial (PP) representation. C***LIBRARY SLATEC C***CATEGORY H2A2A1, E3, K6 C***TYPE DOUBLE PRECISION (PPQAD-S, DPPQAD-D) C***KEYWORDS B-SPLINE, DATA FITTING, INTERPOLATION, QUADRATURE, SPLINES C***AUTHOR Amos, D. E., (SNLA) C***DESCRIPTION C C Abstract **** a double precision routine **** C DPPQAD computes the integral on (X1,X2) of a K-th order C B-spline using the piecewise polynomial representation C (C,XI,LXI,K). Here the Taylor expansion about the left C end point XI(J) of the J-th interval is integrated and C evaluated on subintervals of (X1,X2) which are formed by C included break points. Integration outside (XI(1),XI(LXI+1)) C is permitted. C C Description of Arguments C Input C,XI,X1,X2 are double precision C LDC - leading dimension of matrix C, LDC .GE. K C C(I,J) - right Taylor derivatives at XI(J), I=1,K , J=1,LXI C XI(*) - break point array of length LXI+1 C LXI - number of polynomial pieces C K - order of B-spline, K .GE. 1 C X1,X2 - end points of quadrature interval, normally in C XI(1) .LE. X .LE. XI(LXI+1) C C Output PQUAD is double precision C PQUAD - integral of the PP representation over (X1,X2) C C Error Conditions C Improper input is a fatal error C C***REFERENCES D. E. Amos, Quadrature subroutines for splines and C B-splines, Report SAND79-1825, Sandia Laboratories, C December 1979. C***ROUTINES CALLED DINTRV, XERMSG C***REVISION HISTORY (YYMMDD) C 800901 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890911 Removed unnecessary intrinsics. (WRB) C 890911 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DPPQAD C INTEGER I, II, IL, ILO, IL1, IL2, IM, K, LDC, LEFT, LXI, MF1, MF2 DOUBLE PRECISION A,AA,BB,C,DX,FLK,PQUAD,Q,S,SS,TA,TB,X,XI,X1,X2 DIMENSION XI(*), C(LDC,*), SS(2) C C***FIRST EXECUTABLE STATEMENT DPPQAD PQUAD = 0.0D0 IF(K.LT.1) GO TO 100 IF(LXI.LT.1) GO TO 105 IF(LDC.LT.K) GO TO 110 AA = MIN(X1,X2) BB = MAX(X1,X2) IF (AA.EQ.BB) RETURN ILO = 1 CALL DINTRV(XI, LXI, AA, ILO, IL1, MF1) CALL DINTRV(XI, LXI, BB, ILO, IL2, MF2) Q = 0.0D0 DO 40 LEFT=IL1,IL2 TA = XI(LEFT) A = MAX(AA,TA) IF (LEFT.EQ.1) A = AA TB = BB IF (LEFT.LT.LXI) TB = XI(LEFT+1) X = MIN(BB,TB) DO 30 II=1,2 SS(II) = 0.0D0 DX = X - XI(LEFT) IF (DX.EQ.0.0D0) GO TO 20 S = C(K,LEFT) FLK = K IM = K - 1 IL = IM DO 10 I=1,IL S = S*DX/FLK + C(IM,LEFT) IM = IM - 1 FLK = FLK - 1.0D0 10 CONTINUE SS(II) = S*DX 20 CONTINUE X = A 30 CONTINUE Q = Q + (SS(1)-SS(2)) 40 CONTINUE IF (X1.GT.X2) Q = -Q PQUAD = Q RETURN C C 100 CONTINUE CALL XERMSG ('SLATEC', 'DPPQAD', 'K DOES NOT SATISFY K.GE.1', 2, + 1) RETURN 105 CONTINUE CALL XERMSG ('SLATEC', 'DPPQAD', 'LXI DOES NOT SATISFY LXI.GE.1', + 2, 1) RETURN 110 CONTINUE CALL XERMSG ('SLATEC', 'DPPQAD', 'LDC DOES NOT SATISFY LDC.GE.K', + 2, 1) RETURN END