1 subroutine DSBASD(KORDER, LEFT, T, X, IDERIV, BDERIV) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5c>> 1994-10-20 DSBASD Krogh Changes to use M77CON 6c>> 1993-01-12 DSBASD CLL Bringing all computation inline. 7c>> 1992-11-02 C. L. Lawson, JPL 8c>> 1988-03-22 C. L. Lawson, JPL 9c 10c Given the knot-set T(i), i = 1, ..., NT, this subr computes 11c values at X of the derivative of order IDERIV of each of the 12c the KORDER B-spline basis functions of order KORDER 13c that are nonzero on the open interval (T(LEFT), T(LEFT+1). 14c Require T(LEFT) < T(LEFT+1) and KORDER .le. LEFT .le. 15c NT-KORDER. 16c For proper evaluation, X must lie in the closed interval 17c [T(LEFT), T(LEFT+1], otherwise the computation constitutes 18c extrapolation. 19c The basis functions whose derivatives are returned will be those 20c indexed from LEFT+1-KORDER through LEFT, and their values will be 21c stored in BDERIV(I), I = 1, ..., KORDER. 22c ------------------------------------------------------------------ 23c Method: 24c 25c In general there are two stages: First compute values of basis 26c functions of order KORDER-IDERIV. Then, if IDERIV > 0, transform 27c these values to produce values of higher order derivatives of 28c higher order basis functions, ultimately obtaining values of 29c the IDERIV derivative of basis functions of order KORDER. 30c 31c The first stage uses the recursion: 32C 33C X - T(I) T(I+J+1) - X 34C B(I,J+1)(X) = -----------B(I,J)(X) + ---------------B(I+1,J)(X) 35C T(I+J)-T(I) T(I+J+1)-T(I+1) 36C 37c where B(I,J) denotes the basis function of order J and index I. 38c For order J, the only basis functions that can be nonzero on 39c [T(LEFT), T(LEFT+1] are those indexed from LEFT+1-J to LEFT. 40c For order 1, the only basis function nonzero on this interval is 41c is B(LEFT,1), whose value is 1. The organization of the calculation 42c using the above recursion follows Algorithm (8) in Chapter X of the 43c book by DeBoor. 44c 45c For the second stage, let B(ID, K, J) denote the value at X of the 46c IDth derivative of the basis function of order K, with index J. 47c From the first stage we will have values of B(0, KORDER-IDERIV, j) 48c for j = LEFT+1-KORDER+IDERIV), ..., LEFT, stored in BDERIV(i), for 49c i = 1+IDERIV, ..., KORDER. 50c 51c Loop for ID = 1, ..., IDERIV 52c Replace the contents of BDERIV() by values of 53c B(ID, KORDER-IDERIV+ID, j) for 54c j = LEFT+1-KORDER+IDERIV-ID), ..., LEFT, storing these in 55c BDERIV(i), i = 1+IDERIV-ID, ..., KORDER. 56c End loop 57c 58c The above loop uses formula (10), p. 138, of 59c A PRACTICAL GUIDE TO SPLINES by Carl DeBoor, Springer-Verlag, 60c 1978, when ID = 1, and successive derivatives of that formula 61c when ID > 1. Note that we are using Formula (10) in the special 62c case in which (Alpha sub i) = 1, and all other Alpha's are zero. 63c 64c This approach and implementation by C. L. Lawson, JPL, March 1988. 65C ------------------------------------------------------------------ 66C KORDER [in] Gives both the order of the B-spline basis functions 67c whose derivatives are to be evaluated, 68c and the number of such functions to be evaluated. 69c Note that the polynomial degree of these functions is one less 70c than the order. Example: For cubic splines the order is 4. 71c Require 1 .le. KORDER .le. KMAX, where KMAX is an internal 72c parameter. 73c LEFT [in] Index identifying the left end of the interval 74c [T(LEFT), T(LEFT+1)] relative to which the B-spline basis 75c functions will be evaluated. 76c Require KORDER .le. LEFT .le. NT-KORDER 77c and T(LEFT) < T(LEFT+1) 78c The evaluation is proper if T(LEFT) .le. X .le. T(LEFT+1), and 79c otherwise constitutes extrapolation. 80C DIVISION BY ZERO will result if T(LEFT) = T(LEFT+1) 81C T() [in] Knot sequence, indexed from 1 to NT, with 82c NT .ge. LEFT+KORDER. Knot values must be nonincreasing. 83c Repetition of values is permitted and has the effect of 84c decreasing the order of contimuity at the repeated knot. 85c Proper function representation by splines of order K is 86c supported on the interval from T(K) to T(NT+1-K). 87c Extrapolation can be done outside this interval. 88C X [in] The abcissa at which the B-spline basis functions are to 89c be evaluated. The evaluation is proper if T(KORDER) .le. X .le 90c T(NT+1-KORDER), and constitutes extrapolation otherwise. 91c IDERIV [in] Order of derivative requested. Require IDERIV .ge. 0. 92c Derivatives of order .ge. KORDER are zero. 93C BDERIV() [out] On normal return, this array will contain in 94c BDERIV(i), i = 1, ...,KORDER, the values computed by this subr. 95c These are values at X of the derivative of order IDERIV of each 96c of the basis functions indexed LEFT+1-KORDER through LEFT. 97C ------------------------------------------------------------------ 98c--D replaces "?": ?SBASD 99c Both versions use IERM1, IERV1 100C ------------------------------------------------------------------ 101 integer KMAX 102 parameter(KMAX = 20) 103 integer I, I1, ID, IDERIV, ISTART, IT1, IT2 104 integer J, JP1, K1, KORDER, LEFT 105 double precision BDERIV(KORDER), DELTAL(KMAX), DELTAR(KMAX), FLK1 106 double precision SAVED, T(LEFT+KORDER), TERM, X 107C ------------------------------------------------------------------ 108 ID = IDERIV 109 if(ID .lt. 0 ) then 110 call IERM1('DSBASD',2,2,'Require IDERIV .ge. 0', 111 * 'IDERIV',ID,'.') 112 return 113 endif 114 if(ID .ge. KORDER) then 115 do 5 I = 1, KORDER 116 BDERIV(I) = 0.0d0 117 5 continue 118 return 119 endif 120 if(KORDER .gt. KMAX) then 121 call IERM1('DSBASD',1,2,'Require KORDER .le. KMAX.', 122 * 'KORDER',KORDER,',') 123 call IERV1('KMAX',KMAX,'.') 124 return 125 endif 126 ISTART = 1+ID 127 K1 = KORDER - ID 128c 129c Evaluate K1 basis functions of order K1. Store values in 130c BDERIV(ISTART:ISTART+K1-1) 131c 132 BDERIV(ISTART) = 1.0d0 133 do 30 J = 1, K1-1 134 JP1 = J + 1 135 DELTAR(J) = T(LEFT+J) - X 136 DELTAL(J) = X - T(LEFT+1-J) 137 SAVED = 0.0d0 138 do 26 I=1,J 139 TERM = BDERIV(ID+I)/(DELTAR(I) + DELTAL(JP1-I)) 140 BDERIV(ID+I) = SAVED + DELTAR(I)*TERM 141 SAVED = DELTAL(JP1-I)*TERM 142 26 continue 143 BDERIV(ID+JP1) = SAVED 144 30 continue 145c 146c Loop IDERIV times, each time advancing the order of the 147c derivative by one, so final results are values of derivatives 148c of order IDERIV. 149c 150 do 70 I1 = ISTART, 2, -1 151 FLK1 = dble(K1) 152 IT1 = LEFT+1-I1 153 IT2 = IT1 - K1 154 do 50 I = I1, KORDER 155 BDERIV(I) = BDERIV(I) * FLK1 /(T(IT1+I) - T(IT2+I)) 156 50 continue 157c 158 BDERIV(I1-1) = -BDERIV(I1) 159 do 60 I = I1, KORDER-1 160 BDERIV(I) = BDERIV(I) - BDERIV(I+1) 161 60 continue 162 K1 = K1 + 1 163 70 continue 164 return 165 end 166