1 SUBROUTINE BINT4(X,Y,NDATA,IBCL,IBCR,FBCL,FBCR,KNTOPT,T,BCOEF,N, 2 1 K,W) 3C***BEGIN PROLOGUE BINT4 4C***DATE WRITTEN 800901 (YYMMDD) 5C***REVISION DATE 820801 (YYMMDD) 6C***CATEGORY NO. E1A 7C***KEYWORDS B-SPLINE,DATA FITTING,INTERPOLATION,SPLINE 8C***AUTHOR AMOS, D. E., (SNLA) 9C***PURPOSE Computes the B representation of a cubic spline 10C which interpolates data (X(I),Y(I)),I=1,NDATA. 11C***DESCRIPTION 12C 13C Written by D. E. Amos, August, 1979. 14C 15C References 16C SAND-78-1968 17C 18C A Practical Guide to Splines by C. de Boor, Applied 19C Mathematics Series 27, Springer, 1979. 20C 21C SIAM J. Numerical Analysis, 14, No. 3, June 1977, pp. 441-472. 22C 23C Abstract 24C BINT4 computes the B representation (T,BCOEF,N,K) of a 25C cubic spline (K=4) which interpolates data (X(I)),Y(I))), 26C I=1,NDATA. Parameters IBCL, IBCR, FBCL, FBCR allow the 27C specification of the spline first or second derivative at 28C both X(1) and X(NDATA). When this data is not specified 29C by the problem, it is common practice to use a natural 30C spline by setting second derivatives at X(1) and X(NDATA) 31C to zero (IBCL=IBCR=2,FBCL=FBCR=0.0). The spline is defined on 32C T(4) .LE. X .LE. T(N+1) with (ordered) interior knots at X(I)) 33C values where N=NDATA+2. The knots T(1), T(2), T(3) lie to 34C the left of T(4)=X(1) and the knots T(N+2), T(N+3), T(N+4) 35C lie to the right of T(N+1)=X(NDATA) in increasing order. If 36C no extrapolation outside (X(1),X(NDATA)) is anticipated, the 37C knots T(1)=T(2)=T(3)=T(4)=X(1) and T(N+2)=T(N+3)=T(N+4)= 38C T(N+1)=X(NDATA) can be specified by KNTOPT=1. KNTOPT=2 39C selects a knot placement for T(1), T(2), T(3) to make the 40C first 7 knots symmetric about T(4)=X(1) and similarly for 41C T(N+2), T(N+3), T(N+4) about T(N+1)=X(NDATA). KNTOPT=3 42C allows the user to make his own selection, in increasing 43C order, for T(1), T(2), T(3) to the left of X(1) and T(N+2), 44C T(N+3), T(N+4) to the right of X(NDATA) in the work array 45C W(1) through W(6). In any case, the interpolation on 46C T(4) .LE. X .LE. T(N+1) by using function BVALU is unique 47C for given boundary conditions. 48C 49C BINT4 calls BSPVD, BNFAC, BNSLV, R1MACH, XERROR 50C 51C Description of Arguments 52C Input 53C X - X vector of abscissae of length NDATA, distinct 54C and in increasing order 55C Y - Y vector of ordinates of length NDATA 56C NDATA - number of data points, NDATA .GE. 2 57C IBCL - selection parameter for left boundary condition 58C IBCL = 1 constrain the first derivative at 59C X(1) to FBCL 60C = 2 constrain the second derivative at 61C X(1) to FBCL 62C IBCR - selection parameter for right boundary condition 63C IBCR = 1 constrain first derivative at 64C X(NDATA) to FBCR 65C IBCR = 2 constrain second derivative at 66C X(NDATA) to FBCR 67C FBCL - left boundary values governed by IBCL 68C FBCR - right boundary values governed by IBCR 69C KNTOPT - knot selection parameter 70C KNTOPT = 1 sets knot multiplicity at T(4) and 71C T(N+1) to 4 72C = 2 sets a symmetric placement of knots 73C about T(4) and T(N+1) 74C = 3 sets TNP)=WNP) and T(N+1+I)=w(3+I),I=1,3 75C where WNP),I=1,6 is supplied by the user 76C W - work array of dimension at least 5*(NDATA+2) 77C if KNTOPT=3, then W(1),W(2),W(3) are knot values to 78C the left of X(1) and W(4),W(5),W(6) are knot 79C values to the right of X(NDATA) in increasing 80C order to be supplied by the user 81C 82C Output 83C T - knot array of length N+4 84C BCOEF - B-spline coefficient array of length N 85C N - number of coefficients, N=NDATA+2 86C K - order of spline, K=4 87C 88C Error Conditions 89C Improper input is a fatal error 90C Singular system of equations is a fatal error 91C***REFERENCES D.E. AMOS, *COMPUTATION WITH SPLINES AND B-SPLINES*, 92C SAND78-1968,SANDIA LABORATORIES,MARCH,1979. 93C C. DE BOOR, *PACKAGE FOR CALCULATING WITH B-SPLINES*, 94C SIAM JOURNAL ON NUMERICAL ANALYSIS, VOLUME 14, NO. 3, 95C JUNE 1977, PP. 441-472. 96C C. DE BOOR, *A PRACTICAL GUIDE TO SPLINES*, APPLIED 97C MATHEMATICS SERIES 27, SPRINGER, 1979. 98C***ROUTINES CALLED BNFAC,BNSLV,BSPVD,R1MACH,XERROR 99C***END PROLOGUE BINT4 100C 101C 102 INTEGER I, IBCL, IBCR, IFLAG, ILB, ILEFT, IT, IUB, IW, IWP, J, 103 1 JW, K, KNTOPT, N, NDATA, NDM, NP, NWROW 104 REAL BCOEF,FBCL,FBCR,T, TOL,TXN,TX1,VNIKX,W,WDTOL,WORK,X, XL, 105 1 Y 106 REAL R1MACH 107 DIMENSION X(1), Y(1), T(1), BCOEF(1), W(5,1), VNIKX(4,4), WORK(15) 108C***FIRST EXECUTABLE STATEMENT BINT4 109 WDTOL = R1MACH(4) 110 TOL = SQRT(WDTOL) 111 IF (NDATA.LT.2) GO TO 200 112 NDM = NDATA - 1 113 DO 10 I=1,NDM 114 IF (X(I).GE.X(I+1)) GO TO 210 115 10 CONTINUE 116 IF (IBCL.LT.1 .OR. IBCL.GT.2) GO TO 220 117 IF (IBCR.LT.1 .OR. IBCR.GT.2) GO TO 230 118 IF (KNTOPT.LT.1 .OR. KNTOPT.GT.3) GO TO 240 119 K = 4 120 N = NDATA + 2 121 NP = N + 1 122 DO 20 I=1,NDATA 123 T(I+3) = X(I) 124 20 CONTINUE 125 GO TO (30, 50, 90), KNTOPT 126C SET UP KNOT ARRAY WITH MULTIPLICITY 4 AT X(1) AND X(NDATA) 127 30 CONTINUE 128 DO 40 I=1,3 129 T(4-I) = X(1) 130 T(NP+I) = X(NDATA) 131 40 CONTINUE 132 GO TO 110 133C SET UP KNOT ARRAY WITH SYMMETRIC PLACEMENT ABOUT END POINTS 134 50 CONTINUE 135 IF (NDATA.GT.3) GO TO 70 136 XL = (X(NDATA)-X(1))/3.0E0 137 DO 60 I=1,3 138 T(4-I) = T(5-I) - XL 139 T(NP+I) = T(NP+I-1) + XL 140 60 CONTINUE 141 GO TO 110 142 70 CONTINUE 143 TX1 = X(1) + X(1) 144 TXN = X(NDATA) + X(NDATA) 145 DO 80 I=1,3 146 T(4-I) = TX1 - X(I+1) 147 T(NP+I) = TXN - X(NDATA-I) 148 80 CONTINUE 149 GO TO 110 150C SET UP KNOT ARRAY LESS THAN X(1) AND GREATER THAN X(NDATA) TO BE 151C SUPPLIED BY USER IN WORK LOCATIONS W(1) THROUGH W(6) WHEN KNTOPT=3 152 90 CONTINUE 153 DO 100 I=1,3 154 T(4-I) = W(4-I,1) 155 JW = MAX0(1,I-1) 156 IW = MOD(I+2,5)+1 157 T(NP+I) = W(IW,JW) 158 IF (T(4-I).GT.T(5-I)) GO TO 250 159 IF (T(NP+I).LT.T(NP+I-1)) GO TO 250 160 100 CONTINUE 161 110 CONTINUE 162C 163 DO 130 I=1,5 164 DO 120 J=1,N 165 W(I,J) = 0.0E0 166 120 CONTINUE 167 130 CONTINUE 168C SET UP LEFT INTERPOLATION POINT AND LEFT BOUNDARY CONDITION FOR 169C RIGHT LIMITS 170 IT = IBCL + 1 171 CALL BSPVD(T, K, IT, X(1), K, 4, VNIKX, WORK) 172 IW = 0 173 IF (ABS(VNIKX(3,1)).LT.TOL) IW = 1 174 DO 140 J=1,3 175 W(J+1,4-J) = VNIKX(4-J,IT) 176 W(J,4-J) = VNIKX(4-J,1) 177 140 CONTINUE 178 BCOEF(1) = Y(1) 179 BCOEF(2) = FBCL 180C SET UP INTERPOLATION EQUATIONS FOR POINTS I=2 TO I=NDATA-1 181 ILEFT = 4 182 IF (NDM.LT.2) GO TO 170 183 DO 160 I=2,NDM 184 ILEFT = ILEFT + 1 185 CALL BSPVD(T, K, 1, X(I), ILEFT, 4, VNIKX, WORK) 186 DO 150 J=1,3 187 W(J+1,3+I-J) = VNIKX(4-J,1) 188 150 CONTINUE 189 BCOEF(I+1) = Y(I) 190 160 CONTINUE 191C SET UP RIGHT INTERPOLATION POINT AND RIGHT BOUNDARY CONDITION FOR 192C LEFT LIMITS(ILEFT IS ASSOCIATED WITH T(N)=X(NDATA-1)) 193 170 CONTINUE 194 IT = IBCR + 1 195 CALL BSPVD(T, K, IT, X(NDATA), ILEFT, 4, VNIKX, WORK) 196 JW = 0 197 IF (ABS(VNIKX(2,1)).LT.TOL) JW = 1 198 DO 180 J=1,3 199 W(J+1,3+NDATA-J) = VNIKX(5-J,IT) 200 W(J+2,3+NDATA-J) = VNIKX(5-J,1) 201 180 CONTINUE 202 BCOEF(N-1) = FBCR 203 BCOEF(N) = Y(NDATA) 204C SOLVE SYSTEM OF EQUATIONS 205 ILB = 2 - JW 206 IUB = 2 - IW 207 NWROW = 5 208 IWP = IW + 1 209 CALL BNFAC(W(IWP,1), NWROW, N, ILB, IUB, IFLAG) 210 IF (IFLAG.EQ.2) GO TO 190 211 CALL BNSLV(W(IWP,1), NWROW, N, ILB, IUB, BCOEF) 212 RETURN 213C 214C 215 190 CONTINUE 216 CALL XERROR( ' BINT4, THE SYSTEM OF EQUATIONS IS SINGULAR', 217 1 44, 2, 1) 218 RETURN 219 200 CONTINUE 220 CALL XERROR( ' BINT4, NDATA IS LESS THAN 2', 29, 2, 1) 221 RETURN 222 210 CONTINUE 223 CALL XERROR( ' BINT4, X VALUES ARE NOT DISTINCT OR NOT ORDERED', 224 149, 2, 1) 225 RETURN 226 220 CONTINUE 227 CALL XERROR( ' BINT4, IBCL IS NOT 1 OR 2', 27, 2, 1) 228 RETURN 229 230 CONTINUE 230 CALL XERROR( ' BINT4, IBCR IS NOT 1 OR 2', 27, 2, 1) 231 RETURN 232 240 CONTINUE 233 CALL XERROR( ' BINT4, KNTOPT IS NOT 1, 2, OR 3', 33, 2, 1) 234 RETURN 235 250 CONTINUE 236 CALL XERROR( ' BINT4, KNOT INPUT THROUGH W ARRAY IS NOT ORDERED P 237 1ROPERLY', 59, 2, 1) 238 RETURN 239 END 240