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