1 SUBROUTINE DBSPVD(T,K,NDERIV,X,ILEFT,LDVNIK,VNIKX,WORK) 2C***BEGIN PROLOGUE DBSPVD 3C***DATE WRITTEN 800901 (YYMMDD) 4C***REVISION DATE 820801 (YYMMDD) 5C***REVISION HISTORY (YYMMDD) 6C 000330 Modified array declarations. (JEC) 7C 8C***CATEGORY NO. E3,K6 9C***KEYWORDS B-SPLINE,DATA FITTING,DOUBLE PRECISION,INTERPOLATION, 10C SPLINE 11C***AUTHOR AMOS, D. E., (SNLA) 12C***PURPOSE Calculates the value and all derivatives of order less than 13C NDERIV of all basis functions which do not vanish at X. 14C***DESCRIPTION 15C 16C Written by Carl de Boor and modified by D. E. Amos 17C 18C Reference 19C SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. 20C 21C Abstract **** a double precision routine **** 22C 23C DBSPVD is the BSPLVD routine of the reference. 24C 25C DBSPVD calculates the value and all derivatives of order 26C less than NDERIV of all basis functions which do not 27C (possibly) vanish at X. ILEFT is input such that 28C T(ILEFT) .LE. X .LT. T(ILEFT+1). A call to INTRV(T,N+1,X, 29C ILO,ILEFT,MFLAG) will produce the proper ILEFT. The output of 30C DBSPVD is a matrix VNIKX(I,J) of dimension at least (K,NDERIV) 31C whose columns contain the K nonzero basis functions and 32C their NDERIV-1 right derivatives at X, I=1,K, J=1,NDERIV. 33C These basis functions have indices ILEFT-K+I, I=1,K, 34C K .LE. ILEFT .LE. N. The nonzero part of the I-th basis 35C function lies in (T(I),T(I+K)), I=1,N). 36C 37C If X=T(ILEFT+1) then VNIKX contains left limiting values 38C (left derivatives) at T(ILEFT+1). In particular, ILEFT = N 39C produces left limiting values at the right end point 40C X=T(N+1). To obtain left limiting values at T(I), I=K+1,N+1, 41C set X= next lower distinct knot, call INTRV to get ILEFT, 42C set X=T(I), and then call DBSPVD. 43C 44C DBSPVD calls DBSPVN 45C 46C Description of Arguments 47C Input T,X are double precision 48C T - knot vector of length N+K, where 49C N = number of B-spline basis functions 50C N = sum of knot multiplicities-K 51C K - order of the B-spline, K .GE. 1 52C NDERIV - number of derivatives = NDERIV-1, 53C 1 .LE. NDERIV .LE. K 54C X - argument of basis functions, 55C T(K) .LE. X .LE. T(N+1) 56C ILEFT - largest integer such that 57C T(ILEFT) .LE. X .LT. T(ILEFT+1) 58C LDVNIK - leading dimension of matrix VNIKX 59C 60C Output VNIKX,WORK are double precision 61C VNIKX - matrix of dimension at least (K,NDERIV) contain- 62C ing the nonzero basis functions at X and their 63C derivatives columnwise. 64C WORK - a work vector of length (K+1)*(K+2)/2 65C 66C Error Conditions 67C Improper input is a fatal error 68C***REFERENCES C. DE BOOR, *PACKAGE FOR CALCULATING WITH B-SPLINES*, 69C SIAM JOURNAL ON NUMERICAL ANALYSIS, VOLUME 14, NO. 3, 70C JUNE 1977, PP. 441-472. 71C***ROUTINES CALLED DBSPVN,XERROR 72C***END PROLOGUE DBSPVD 73C 74C 75 INTEGER I,IDERIV,ILEFT,IPKMD,J,JJ,JLOW,JM,JP1MID,K,KMD, KP1, L, 76 1 LDUMMY, M, MHIGH, NDERIV 77 DOUBLE PRECISION FACTOR, FKMD, T, V, VNIKX, WORK, X 78C DIMENSION T(ILEFT+K), WORK((K+1)*(K+2)/2) 79C A(I,J) = WORK(I+J*(J+1)/2), I=1,J+1 J=1,K-1 80C A(I,K) = W0RK(I+K*(K-1)/2) I=1.K 81C WORK(1) AND WORK((K+1)*(K+2)/2) ARE NOT USED. 82 DIMENSION T(*), VNIKX(LDVNIK,NDERIV), WORK(*) 83C***FIRST EXECUTABLE STATEMENT DBSPVD 84 IF(K.LT.1) GO TO 200 85 IF(NDERIV.LT.1 .OR. NDERIV.GT.K) GO TO 205 86 IF(LDVNIK.LT.K) GO TO 210 87 IDERIV = NDERIV 88 KP1 = K + 1 89 JJ = KP1 - IDERIV 90 CALL DBSPVN(T, JJ, K, 1, X, ILEFT, VNIKX, WORK, IWORK) 91 IF (IDERIV.EQ.1) GO TO 100 92 MHIGH = IDERIV 93 DO 20 M=2,MHIGH 94 JP1MID = 1 95 DO 10 J=IDERIV,K 96 VNIKX(J,IDERIV) = VNIKX(JP1MID,1) 97 JP1MID = JP1MID + 1 98 10 CONTINUE 99 IDERIV = IDERIV - 1 100 JJ = KP1 - IDERIV 101 CALL DBSPVN(T, JJ, K, 2, X, ILEFT, VNIKX, WORK, IWORK) 102 20 CONTINUE 103C 104 JM = KP1*(KP1+1)/2 105 DO 30 L = 1,JM 106 WORK(L) = 0.0D0 107 30 CONTINUE 108C A(I,I) = WORK(I*(I+3)/2) = 1.0 I = 1,K 109 L = 2 110 J = 0 111 DO 40 I = 1,K 112 J = J + L 113 WORK(J) = 1.0D0 114 L = L + 1 115 40 CONTINUE 116 KMD = K 117 DO 90 M=2,MHIGH 118 KMD = KMD - 1 119 FKMD = FLOAT(KMD) 120 I = ILEFT 121 J = K 122 JJ = J*(J+1)/2 123 JM = JJ - J 124 DO 60 LDUMMY=1,KMD 125 IPKMD = I + KMD 126 FACTOR = FKMD/(T(IPKMD)-T(I)) 127 DO 50 L=1,J 128 WORK(L+JJ) = (WORK(L+JJ)-WORK(L+JM))*FACTOR 129 50 CONTINUE 130 I = I - 1 131 J = J - 1 132 JJ = JM 133 JM = JM - J 134 60 CONTINUE 135C 136 DO 80 I=1,K 137 V = 0.0D0 138 JLOW = MAX0(I,M) 139 JJ = JLOW*(JLOW+1)/2 140 DO 70 J=JLOW,K 141 V = WORK(I+JJ)*VNIKX(J,M) + V 142 JJ = JJ + J + 1 143 70 CONTINUE 144 VNIKX(I,M) = V 145 80 CONTINUE 146 90 CONTINUE 147 100 RETURN 148C 149C 150 200 CONTINUE 151 CALL XERROR( ' DBSPVD, K DOES NOT SATISFY K.GE.1', 35, 2, 1) 152 RETURN 153 205 CONTINUE 154 CALL XERROR( ' DBSPVD, NDERIV DOES NOT SATISFY 1.LE.NDERIV.LE.K', 155 1 50, 2, 1) 156 RETURN 157 210 CONTINUE 158 CALL XERROR( ' DBSPVD, LDVNIK DOES NOT SATISFY LDVNIK.GE.K',45, 159 1 2, 1) 160 RETURN 161 END 162