1 FUNCTION BVALU(T,A,N,K,IDERIV,X,INBV,WORK) 2C***BEGIN PROLOGUE BVALU 3C***DATE WRITTEN 800901 (YYMMDD) 4C***REVISION DATE 820801 (YYMMDD) 5C***CATEGORY NO. E3,K6 6C***KEYWORDS B-SPLINE,DATA FITTING,INTERPOLATION,SPLINE 7C***AUTHOR AMOS, D. E., (SNLA) 8C***PURPOSE Evaluates the B-representation of a B-spline at X for the 9C function value or any of its derivatives. 10C***DESCRIPTION 11C 12C Written by Carl de Boor and modified by D. E. Amos 13C 14C Reference 15C SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. 16C 17C Abstract 18C BVALU is the BVALUE function of the reference. 19C 20C BVALU evaluates the B-representation (T,A,N,K) of a B-spline 21C at X for the function value on IDERIV = 0 or any of its 22C derivatives on IDERIV = 1,2,...,K-1. Right limiting values 23C (right derivatives) are returned except at the right end 24C point X=T(N+1) where left limiting values are computed. The 25C spline is defined on T(K) .LE. X .LE. T(N+1). BVALU returns 26C a fatal error message when X is outside of this interval. 27C 28C To compute left derivatives or left limiting values at a 29C knot T(I), replace N by I-1 and set X=T(I), I=K+1,N+1. 30C 31C BVALU calls INTRV 32C 33C Description of Arguments 34C Input 35C T - knot vector of length N+K 36C A - B-spline coefficient vector of length N 37C N - number of B-spline coefficients 38C N = sum of knot multiplicities-K 39C K - order of the B-spline, K .GE. 1 40C IDERIV - order of the derivative, 0 .LE. IDERIV .LE. K-1 41C IDERIV=0 returns the B-spline value 42C X - argument, T(K) .LE. X .LE. T(N+1) 43C INBV - an initialization parameter which must be set 44C to 1 the first time BVALU is called. 45C 46C Output 47C INBV - INBV contains information for efficient process- 48C ing after the initial call and INBV must not 49C be changed by the user. Distinct splines require 50C distinct INBV parameters. 51C WORK - work vector of length 3*K. 52C BVALU - value of the IDERIV-th derivative at X 53C 54C Error Conditions 55C An improper input is a fatal error 56C***REFERENCES C. DE BOOR, *PACKAGE FOR CALCULATING WITH B-SPLINES*, 57C SIAM JOURNAL ON NUMERICAL ANALYSIS, VOLUME 14, NO. 3, 58C JUNE 1977, PP. 441-472. 59C***ROUTINES CALLED INTRV,XERROR 60C***END PROLOGUE BVALU 61C 62C 63 INTEGER I,IDERIV,IDERP1,IHI,IHMKMJ,ILO,IMK,IMKPJ, INBV, IPJ, 64 1 IP1, IP1MJ, J, JJ, J1, J2, K, KMIDER, KMJ, KM1, KPK, MFLAG, N 65 REAL A, FKMJ, T, WORK, X 66C DIMENSION T(N+K), WORK(3*K) 67 DIMENSION T(1), A(N), WORK(1) 68C***FIRST EXECUTABLE STATEMENT BVALU 69 BVALU = 0.0E0 70 IF(K.LT.1) GO TO 102 71 IF(N.LT.K) GO TO 101 72 IF(IDERIV.LT.0 .OR. IDERIV.GE.K) GO TO 110 73 KMIDER = K - IDERIV 74C 75C *** FIND *I* IN (K,N) SUCH THAT T(I) .LE. X .LT. T(I+1) 76C (OR, .LE. T(I+1) IF T(I) .LT. T(I+1) = T(N+1)). 77 KM1 = K - 1 78 CALL INTRV(T, N+1, X, INBV, I, MFLAG) 79 IF (X.LT.T(K)) GO TO 120 80 IF (MFLAG.EQ.0) GO TO 20 81 IF (X.GT.T(I)) GO TO 130 82 10 IF (I.EQ.K) GO TO 140 83 I = I - 1 84 IF (X.EQ.T(I)) GO TO 10 85C 86C *** DIFFERENCE THE COEFFICIENTS *IDERIV* TIMES 87C WORK(I) = AJ(I), WORK(K+I) = DP(I), WORK(K+K+I) = DM(I), I=1.K 88C 89 20 IMK = I - K 90 DO 30 J=1,K 91 IMKPJ = IMK + J 92 WORK(J) = A(IMKPJ) 93 30 CONTINUE 94 IF (IDERIV.EQ.0) GO TO 60 95 DO 50 J=1,IDERIV 96 KMJ = K - J 97 FKMJ = FLOAT(KMJ) 98 DO 40 JJ=1,KMJ 99 IHI = I + JJ 100 IHMKMJ = IHI - KMJ 101 WORK(JJ) = (WORK(JJ+1)-WORK(JJ))/(T(IHI)-T(IHMKMJ))*FKMJ 102 40 CONTINUE 103 50 CONTINUE 104C 105C *** COMPUTE VALUE AT *X* IN (T(I),(T(I+1)) OF IDERIV-TH DERIVATIVE, 106C GIVEN ITS RELEVANT B-SPLINE COEFF. IN AJ(1),...,AJ(K-IDERIV). 107 60 IF (IDERIV.EQ.KM1) GO TO 100 108 IP1 = I + 1 109 KPK = K + K 110 J1 = K + 1 111 J2 = KPK + 1 112 DO 70 J=1,KMIDER 113 IPJ = I + J 114 WORK(J1) = T(IPJ) - X 115 IP1MJ = IP1 - J 116 WORK(J2) = X - T(IP1MJ) 117 J1 = J1 + 1 118 J2 = J2 + 1 119 70 CONTINUE 120 IDERP1 = IDERIV + 1 121 DO 90 J=IDERP1,KM1 122 KMJ = K - J 123 ILO = KMJ 124 DO 80 JJ=1,KMJ 125 WORK(JJ) = (WORK(JJ+1)*WORK(KPK+ILO)+WORK(JJ) 126 1 *WORK(K+JJ))/(WORK(KPK+ILO)+WORK(K+JJ)) 127 ILO = ILO - 1 128 80 CONTINUE 129 90 CONTINUE 130 100 BVALU = WORK(1) 131 RETURN 132C 133C 134 101 CONTINUE 135 CALL XERROR( ' BVALU, N DOES NOT SATISFY N.GE.K',34,2,1) 136 RETURN 137 102 CONTINUE 138 CALL XERROR( ' BVALU, K DOES NOT SATISFY K.GE.1',34,2,1) 139 RETURN 140 110 CONTINUE 141 CALL XERROR( ' BVALU, IDERIV DOES NOT SATISFY 0.LE.IDERIV.LT.K', 142 1 49, 2, 1) 143 RETURN 144 120 CONTINUE 145 CALL XERROR( ' BVALU, X IS N0T GREATER THAN OR EQUAL TO T(K)', 146 1 47, 2, 1) 147 RETURN 148 130 CONTINUE 149 CALL XERROR( ' BVALU, X IS NOT LESS THAN OR EQUAL TO T(N+1)', 150 1 46, 2, 1) 151 RETURN 152 140 CONTINUE 153 CALL XERROR( ' BVALU, A LEFT LIMITING VALUE CANN0T BE OBTAINED AT 154 1 T(K)', 57, 2, 1) 155 RETURN 156 END 157