Code covered by the BSD License

### Highlights fromslatec

from slatec by Ben Barrowes
The slatec library converted into matlab functions.

[x,y,t,n,k,bcoef,q,work]=dbintk(x,y,t,n,k,bcoef,q,work);
function [x,y,t,n,k,bcoef,q,work]=dbintk(x,y,t,n,k,bcoef,q,work);
%***BEGIN PROLOGUE  DBINTK
%***PURPOSE  Compute the B-representation of a spline which interpolates
%            given data.
%***LIBRARY   SLATEC
%***CATEGORY  E1A
%***TYPE      doubleprecision (BINTK-S, DBINTK-D)
%***KEYWORDS  B-SPLINE, DATA FITTING, INTERPOLATION
%***AUTHOR  Amos, D. E., (SNLA)
%***DESCRIPTION
%
%     Written by Carl de Boor and modified by D. E. Amos
%
%     Abstract    **** a doubleprecision routine ****
%
%         DBINTK is the SPLINT routine of the reference.
%
%         DBINTK produces the B-spline coefficients, BCOEF, of the
%         B-spline of order K with knots T(I), I=1,...,N+K, which
%         takes on the value Y(I) at X(I), I=1,...,N.  The spline or
%         any of its derivatives can be evaluated by calls to DBVALU.
%
%         The I-th equation of the linear system A*BCOEF = B for the
%         coefficients of the interpolant enforces interpolation at
%         X(I), I=1,...,N.  Hence, B(I) = Y(I), for all I, and A is
%         a band matrix with 2K-1 bands if A is invertible.  The matrix
%         A is generated row by row and stored, diagonal by diagonal,
%         in the rows of Q, with the main diagonal going into row K.
%         The banded system is then solved by a call to DBNFAC (which
%         constructs the triangular factorization for A and stores it
%         again in Q), followed by a call to DBNSLV (which then
%         obtains the solution BCOEF by substitution).  DBNFAC does no
%         pivoting, since the total positivity of the matrix A makes
%         this unnecessary.  The linear system to be solved is
%         (theoretically) invertible if and only if
%                 T(I) .LT. X(I) .LT. T(I+K),        for all I.
%         Equality is permitted on the left for I=1 and on the right
%         for I=N when K knots are used at X(1) or X(N).  Otherwise,
%         violation of this condition is certain to lead to an error.
%
%     Description of Arguments
%
%         Input       X,Y,T are doubleprecision
%           X       - vector of length N containing data point abscissa
%                     in strictly increasing order.
%           Y       - corresponding vector of length N containing data
%                     point ordinates.
%           T       - knot vector of length N+K
%                     Since T(1),..,T(K) .LE. X(1) and T(N+1),..,T(N+K)
%                     .GE. X(N), this leaves only N-K knots (not nec-
%                     essarily X(I) values) interior to (X(1),X(N))
%           N       - number of data points, N .GE. K
%           K       - order of the spline, K .GE. 1
%
%         Output      BCOEF,Q,WORK are doubleprecision
%           BCOEF   - a vector of length N containing the B-spline
%                     coefficients
%           Q       - a work vector of length (2*K-1)*N, containing
%                     the triangular factorization of the coefficient
%                     matrix of the linear system being solved.  The
%                     coefficients for the interpolant of an
%                     additional data set (X(I),YY(I)), I=1,...,N
%                     YY into BCOEF and then executing
%                         CALL DBNSLV (Q,2K-1,N,K-1,K-1,BCOEF)
%           WORK    - work vector of length 2*K
%
%     Error Conditions
%         Improper input is a fatal error
%         Singular system of equations is a fatal error
%
%***REFERENCES  D. E. Amos, Computation with splines and B-splines,
%                 Report SAND78-1968, Sandia Laboratories, March 1979.
%               Carl de Boor, Package for calculating with B-splines,
%                 SIAM Journal on Numerical Analysis 14, 3 (June 1977),
%                 pp. 441-472.
%               Carl de Boor, A Practical Guide to Splines, Applied
%                 Mathematics Series 27, Springer-Verlag, New York,
%                 1978.
%***ROUTINES CALLED  DBNFAC, DBNSLV, DBSPVN, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   800901  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  DBINTK
%
persistent i iflag ilp1mx iwork j jj km1 kpkm2 left lenq np1 xi ;

if isempty(iflag), iflag=0; end;
if isempty(iwork), iwork=0; end;
if isempty(i), i=0; end;
if isempty(ilp1mx), ilp1mx=0; end;
if isempty(j), j=0; end;
if isempty(jj), jj=0; end;
if isempty(km1), km1=0; end;
if isempty(kpkm2), kpkm2=0; end;
if isempty(left), left=0; end;
if isempty(lenq), lenq=0; end;
if isempty(np1), np1=0; end;
bcoef_shape=size(bcoef);bcoef=reshape(bcoef,1,[]);
y_shape=size(y);y=reshape(y,1,[]);
q_shape=size(q);q=reshape(q,1,[]);
t_shape=size(t);t=reshape(t,1,[]);
x_shape=size(x);x=reshape(x,1,[]);
if isempty(xi), xi=0; end;
work_shape=size(work);work=reshape(work,1,[]);
%     DIMENSION Q(2*K-1,N), T(N+K)
%***FIRST EXECUTABLE STATEMENT  DBINTK
if( k<1 )
xermsg('SLATEC','DBINTK','K DOES NOT SATISFY K.GE.1',2,1);
elseif( n<k ) ;
xermsg('SLATEC','DBINTK','N DOES NOT SATISFY N.GE.K',2,1);
else;
jj = fix(n - 1);
if( jj~=0 )
for i = 1 : jj;
if( x(i)>=x(i+1) )
xermsg('SLATEC','DBINTK','X(I) DOES NOT SATISFY X(I).LT.X(I+1) FOR SOME I',2,1);
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
end; i = fix(jj+1);
end;
np1 = fix(n + 1);
km1 = fix(k - 1);
kpkm2 = fix(2.*km1);
left = fix(k);
%                ZERO OUT ALL ENTRIES OF Q
lenq = fix(n.*(k+km1));
for i = 1 : lenq;
q(i) = 0.0d0;
end; i = fix(lenq+1);
%
%  ***   LOOP OVER I TO CONSTRUCT THE  N  INTERPOLATION EQUATIONS
for i = 1 : n;
xi = x(i);
ilp1mx = fix(min(i+k,np1));
%        *** FIND  LEFT  IN THE CLOSED INTERVAL (I,I+K-1) SUCH THAT
%                T(LEFT) .LE. X(I) .LT. T(LEFT+1)
%        MATRIX IS SINGULAR IF THIS IS NOT POSSIBLE
left = fix(max(left,i));
if( xi<t(left) )
%
%
xermsg('SLATEC','DBINTK',['SOME ABSCISSA WAS NOT IN THE SUPPORT OF THE CORRESPONDING ','BASIS FUNCTION AND THE SYSTEM IS SINGULAR.'],2,1);
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
else;
while( xi>=t(left+1) );
left = fix(left + 1);
if( left>=ilp1mx )
left = fix(left - 1);
if( xi<=t(left+1) )
break;
end;
xermsg('SLATEC','DBINTK',['SOME ABSCISSA WAS NOT IN THE SUPPORT OF THE CORRESPONDING ','BASIS FUNCTION AND THE SYSTEM IS SINGULAR.'],2,1);
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
end;
%        *** THE I-TH EQUATION ENFORCES INTERPOLATION AT XI, HENCE
%        A(I,J) = B(J,K,T,XI), ALL J. ONLY THE  K  ENTRIES WITH  J =
%        LEFT-K+1,...,LEFT ACTUALLY MIGHT BE NONZERO. THESE  K  NUMBERS
%        ARE RETURNED, IN  BCOEF (USED FOR TEMP. STORAGE HERE), BY THE
%        FOLLOWING
k_orig=k;    [t,k,dumvar3,dumvar4,xi,left,bcoef,work,iwork]=dbspvn(t,k,k,1,xi,left,bcoef,work,iwork);    k(dumvar3~=k_orig)=dumvar3(dumvar3~=k_orig);
%        WE THEREFORE WANT  BCOEF(J) = B(LEFT-K+J,XI) TO GO INTO
%        A(I,LEFT-K+J), I.E., INTO  Q(I-(LEFT+J)+2*K,(LEFT+J)-K) SINCE
%        A(I+J,J)  IS TO GO INTO  Q(I+K,J), ALL I,J,  IF WE CONSIDER  Q
%        AS A TWO-DIM. ARRAY , WITH  2*K-1  ROWS (SEE COMMENTS IN
%        DBNFAC). IN THE PRESENT PROGRAM, WE TREAT  Q  AS AN EQUIVALENT
%        ONE-DIMENSIONAL ARRAY (BECAUSE OF FORTRAN RESTRICTIONS ON
%        DIMENSION STATEMENTS) . WE THEREFORE WANT  BCOEF(J) TO GO INTO
%        ENTRY
%            I -(LEFT+J) + 2*K + ((LEFT+J) - K-1)*(2*K-1)
%                   =  I-LEFT+1 + (LEFT -K)*(2*K-1) + (2*K-2)*J
%        OF  Q .
jj = fix(i - left + 1 +(left-k).*(k+km1));
for j = 1 : k;
jj = fix(jj + kpkm2);
q(jj) = bcoef(j);
end; j = fix(k+1);
end;
end;
%
%     ***OBTAIN FACTORIZATION OF  A  , STORED AGAIN IN  Q.
km1_orig=km1;    [q,dumvar2,n,km1,dumvar5,iflag]=dbnfac(q,k+km1,n,km1,km1,iflag);    km1(dumvar5~=km1_orig)=dumvar5(dumvar5~=km1_orig);
if( iflag==2 )
xermsg('SLATEC','DBINTK',['THE SYSTEM OF SOLVER DETECTS A SINGULAR SYSTEM ALTHOUGH ','THE THEORETICAL CONDITIONS FOR A SOLUTION WERE SATISFIED.'],8,1);
else;
%     *** SOLVE  A*BCOEF = Y  BY BACKSUBSTITUTION
for i = 1 : n;
bcoef(i) = y(i);
end; i = fix(n+1);
km1_orig=km1;    [q,dumvar2,n,km1,dumvar5,bcoef]=dbnslv(q,k+km1,n,km1,km1,bcoef);    km1(dumvar5~=km1_orig)=dumvar5(dumvar5~=km1_orig);
end;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end %subroutine dbintk
%DECK DBKIAS