5 * Revision 1.1.1.1 1996/04/01 15:02:26 mclareni
10 FUNCTION DSPPS1(K,M,NDER,X,T,C,NERR)
12 #include "gen/imp64.inc"
13 DIMENSION T(*),C(*),B(27)
16 PARAMETER (NAME = 'DSPPS1')
18 ************************************************************************
19 * NORBAS, VERSION: 15.03.1993
20 ************************************************************************
22 * DSPPS1 COMPUTES FUNCTION VALUES, VALUES OF DERIVATIVES, AND THE
23 * VALUE OF INTEGRAL, RESPECTIVELY, OF A POLYNOMIAL SPLINE S(X) IN
24 * B-SPLINE REPRESENTATION
26 * S(X) = SUMME(I=1,...,M-K-1) C(I) * B(I,K)(X) .
28 * THE FUNCTIONS B(I,K)(X) ARE NORMALIZED B-SPLINES OF DEGREE K
29 * (0<= K <= 25) WITH INDEX I (1 <= I <= M-K-1) OVER A SET OF
31 * T(1),T(2), ... ,T(M) ( M >= 2*K+2 )
32 * (KNOTS IN ASCENDING ORDER, WITH MULTIPLICITIES NOT GREATER
35 * THE FUNCTION VALUE OF THE NORMALIZED B-SPLINE B(I,K)(X) IS
36 * IDENTICALLY ZERO OUTSIDE THE INTERVAL T(I) <= X < T(I+K+1).
38 * THE NORMALIZATION OF N(X) = B(I,K)(X) IS SUCH THAT THE INTGRAL OF
39 * N(X) OVER THE WHOLE X-RANGE EQUALS
40 * ( T(I+K+1) - T(I) ) / (K+1) .
42 * C(1),...,C(M-K-NDER-1) MUST CONTAIN THE COEFFICIENTS OF THE
43 * POLYNOMIAL SPLINE S(X) OR ITS DERIVATIVE, REPRESENTED BY
44 * NORMALIZED B-SPLINES OF DEGREE K-NDER .
46 * FOR TRANSFORMATION THE COEFFICIENTS OF THE POLYNOMIAL SPLINE S(X)
47 * TO THE CORRESPONDING COEFFICIENTS OF THE NDER-TH DERIVATIVE OF
48 * S(X) THE ROUTINE DSPCD1 MAY BE USED.
50 * ESPECIALLY FOR COMPUTING THE COEFFICIENTS C(1),...,C(M-K-1) OF
51 * THE POLYNOMIAL VARIATION DIMINISHING SPLINE APPROXIMATION OF A USER
52 * SUPPLIED FUNCTION F(X) THE ROUTINE DSPVD1 MAY BE USED.
56 * K (INTEGER) DEGREE (= ORDER - 1) OF B-SPLINES.
57 * M (INTEGER) NUMBER OF KNOTS FOR THE B-SPLINES.
58 * NDER (INTEGER) ON ENTRY, NDER MUST CONTAIN AN INTEGER VALUE .GE. -1
59 * = -1: DSPPS1 COMPUTES THE INTEGRAL OF S(TAU) OVER THE
61 * = 0: DSPPS1 COMPUTES THE FUNCTION VALUE OF THE POLYNOMIAL
63 * >= 1: DSPPS1 COMPUTES THE VALUE OF THE NDER-TH DERIVATIVE OF
64 * S(X) AT X (IF NDER > K ZERO RETURNS).
65 * X (DOUBLE PRECISION) ON ENTRY, X MUST CONTAIN THE VALUE OF THE
66 * INDEPENDENT VARIABLE X OF S(X).
67 * T (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER M CONTAINING THE
69 * C (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER M-K-NDER-1,
70 * ON ENTRY C(J) CONTAINS THE J-TH COEFFICIENT OF THE B-SPLINE
71 * REPRESENTATION OF S(X) OR OF ITS DERIVATIVE.
72 * NERR (INTEGER) ERROR INDICATOR. ON EXIT:
73 * = 0: NO ERROR DETECTED
74 * = 1: AT LEAST ONE OF THE CONSTANTS K , M , NDER IS ILLEGAL
79 * Y = DSPPS1(K,M,NDER,X,T,C,NERR)
81 * - THE VALUE OF THE INTEGRAL (NDER = -1) OR
82 * - THE FUNCTION VALUE (NDER = 0 ) OR
83 * - THE VALUE OF THE NDER-TH DERIVATIVE (NDER > 0 )
84 * OF THE POLYNOMIAL SPLINE (IN B-SPLINE REPRESENTATION)
85 * S(X) = SUMME(I=1,...,M-K-1) C(I)*B(I,K)(X)
86 * OF DEGREE K WITH THE SET OF KNOTS T(1),...,T(M) AT X.
90 * IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT-
91 * PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED:
96 ************************************************************************
98 PARAMETER (Z0 = 0 , Z1 = 1)
101 IF(K .LT. 0 .OR. K .GT. 25) THEN
102 WRITE(ERRTXT,101) 'K',K
103 CALL MTLPRT(NAME,'E210.1',ERRTXT)
104 ELSEIF(M .LT. 2*K+2) THEN
105 WRITE(ERRTXT,101) 'M',M
106 CALL MTLPRT(NAME,'E210.2',ERRTXT)
107 ELSEIF(NDER .LT. -1) THEN
108 WRITE(ERRTXT,101) 'NDER',NDER
109 CALL MTLPRT(NAME,'E210.5',ERRTXT)
115 + .OR. X .GT. T(M-K) .AND. NDER .GE. 0
116 + .OR. K .LT. NDER ) RETURN
118 IF(NDER .EQ. -1) THEN
119 IF(X .GE. T(M-K)) THEN
122 10 R=R+C(I)*(T(I+K+1)-T(I))
125 KK=LKKSPL(X,T(K+1),M-2*K-1)+K
128 20 R=R+C(I)*(T(I+K+1)-T(I))
134 DO 50 I=MAX(1,KK-K-1),KK-1
135 CALL DVSET(K+1,Z0,B(1),B(2))
136 B(KK-I)=1/(T(KK)-T(KK-1))
139 DO 30 J=MAX(1,KK-I-L),MIN(K+1-L,KK-I)
140 DIF=T(I+J+L)-T(I+J-1)
142 IF(DIF .NE. 0) B0=((X-T(I+J-1))*B(J)+(T(I+J+L)-X)*B(J+1))/DIF
146 40 S=S+(X-T(I+L-1))*B(L)
147 S=S*(T(I+K+1)-T(I))/(K+1)
153 KK=LKKSPL(X,T(K+1),M-2*K-1)+K
157 60 B(J)=E1*B(J-1)/(T(KK-2+J)-T(KK-1))
158 IF(KK .NE. 0 .OR. NDER .NE. 0) THEN
162 B(1)=E2*B(1)/(T(KK)-T(KK-J))
164 70 B(L)=E3*B(L-1)/(T(KK-2+L)-T(KK-1-J))+
165 + (T(KK-1+L)-X)*B(L)/(T(KK-1+L)-T(KK-J))
167 R=DVMPY(K-NDER+1,C(KK-1-K),C(KK-K),B(1),B(2))
173 101 FORMAT(1X,A5,' =',I6,' NOT IN RANGE')