1*DECK BSPVN 2 SUBROUTINE BSPVN (T, JHIGH, K, INDEX, X, ILEFT, VNIKX, WORK, 3 + IWORK) 4C***BEGIN PROLOGUE BSPVN 5C***PURPOSE Calculate the value of all (possibly) nonzero basis 6C functions at X. 7C***LIBRARY SLATEC 8C***CATEGORY E3, K6 9C***TYPE SINGLE PRECISION (BSPVN-S, DBSPVN-D) 10C***KEYWORDS EVALUATION OF B-SPLINE 11C***AUTHOR Amos, D. E., (SNLA) 12C***DESCRIPTION 13C 14C Written by Carl de Boor and modified by D. E. Amos 15C 16C Abstract 17C BSPVN is the BSPLVN routine of the reference. 18C 19C BSPVN calculates the value of all (possibly) nonzero basis 20C functions at X of order MAX(JHIGH,(J+1)*(INDEX-1)), where 21C T(K) .LE. X .LE. T(N+1) and J=IWORK is set inside the routine 22C on the first call when INDEX=1. ILEFT is such that T(ILEFT) 23C .LE. X .LT. T(ILEFT+1). A call to INTRV(T,N+1,X,ILO,ILEFT, 24C MFLAG) produces the proper ILEFT. BSPVN calculates using the 25C basic algorithm needed in BSPVD. If only basis functions are 26C desired, setting JHIGH=K and INDEX=1 can be faster than 27C calling BSPVD, but extra coding is required for derivatives 28C (INDEX=2) and BSPVD is set up for this purpose. 29C 30C Left limiting values are set up as described in BSPVD. 31C 32C Description of Arguments 33C Input 34C T - knot vector of length N+K, where 35C N = number of B-spline basis functions 36C N = sum of knot multiplicities-K 37C JHIGH - order of B-spline, 1 .LE. JHIGH .LE. K 38C K - highest possible order 39C INDEX - INDEX = 1 gives basis functions of order JHIGH 40C = 2 denotes previous entry with WORK, IWORK 41C values saved for subsequent calls to 42C BSPVN. 43C X - argument of basis functions, 44C T(K) .LE. X .LE. T(N+1) 45C ILEFT - largest integer such that 46C T(ILEFT) .LE. X .LT. T(ILEFT+1) 47C 48C Output 49C VNIKX - vector of length K for spline values. 50C WORK - a work vector of length 2*K 51C IWORK - a work parameter. Both WORK and IWORK contain 52C information necessary to continue for INDEX = 2. 53C When INDEX = 1 exclusively, these are scratch 54C variables and can be used for other purposes. 55C 56C Error Conditions 57C Improper input is a fatal error. 58C 59C***REFERENCES Carl de Boor, Package for calculating with B-splines, 60C SIAM Journal on Numerical Analysis 14, 3 (June 1977), 61C pp. 441-472. 62C***ROUTINES CALLED XERMSG 63C***REVISION HISTORY (YYMMDD) 64C 800901 DATE WRITTEN 65C 890831 Modified array declarations. (WRB) 66C 890831 REVISION DATE from Version 3.2 67C 891214 Prologue converted to Version 4.0 format. (BAB) 68C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 69C 900326 Removed duplicate information from DESCRIPTION section. 70C (WRB) 71C 920501 Reformatted the REFERENCES section. (WRB) 72C***END PROLOGUE BSPVN 73C 74 INTEGER ILEFT, IMJP1, INDEX, IPJ, IWORK, JHIGH, JP1, JP1ML, K, L 75 REAL T, VM, VMPREV, VNIKX, WORK, X 76C DIMENSION T(ILEFT+JHIGH) 77 DIMENSION T(*), VNIKX(*), WORK(*) 78C CONTENT OF J, DELTAM, DELTAP IS EXPECTED UNCHANGED BETWEEN CALLS. 79C WORK(I) = DELTAP(I), WORK(K+I) = DELTAM(I), I = 1,K 80C***FIRST EXECUTABLE STATEMENT BSPVN 81 IF(K.LT.1) GO TO 90 82 IF(JHIGH.GT.K .OR. JHIGH.LT.1) GO TO 100 83 IF(INDEX.LT.1 .OR. INDEX.GT.2) GO TO 105 84 IF(X.LT.T(ILEFT) .OR. X.GT.T(ILEFT+1)) GO TO 110 85 GO TO (10, 20), INDEX 86 10 IWORK = 1 87 VNIKX(1) = 1.0E0 88 IF (IWORK.GE.JHIGH) GO TO 40 89C 90 20 IPJ = ILEFT + IWORK 91 WORK(IWORK) = T(IPJ) - X 92 IMJP1 = ILEFT - IWORK + 1 93 WORK(K+IWORK) = X - T(IMJP1) 94 VMPREV = 0.0E0 95 JP1 = IWORK + 1 96 DO 30 L=1,IWORK 97 JP1ML = JP1 - L 98 VM = VNIKX(L)/(WORK(L)+WORK(K+JP1ML)) 99 VNIKX(L) = VM*WORK(L) + VMPREV 100 VMPREV = VM*WORK(K+JP1ML) 101 30 CONTINUE 102 VNIKX(JP1) = VMPREV 103 IWORK = JP1 104 IF (IWORK.LT.JHIGH) GO TO 20 105C 106 40 RETURN 107C 108C 109 90 CONTINUE 110 CALL XERMSG ('SLATEC', 'BSPVN', 'K DOES NOT SATISFY K.GE.1', 2, 111 + 1) 112 RETURN 113 100 CONTINUE 114 CALL XERMSG ('SLATEC', 'BSPVN', 115 + 'JHIGH DOES NOT SATISFY 1.LE.JHIGH.LE.K', 2, 1) 116 RETURN 117 105 CONTINUE 118 CALL XERMSG ('SLATEC', 'BSPVN', 'INDEX IS NOT 1 OR 2', 2, 1) 119 RETURN 120 110 CONTINUE 121 CALL XERMSG ('SLATEC', 'BSPVN', 122 + 'X DOES NOT SATISFY T(ILEFT).LE.X.LE.T(ILEFT+1)', 2, 1) 123 RETURN 124 END 125