function sp=spapi(knots,x,y)
% SPAPI	Spline interpolation.
%
%         spapi(knots,x,y)
%
%  returns in  sp  the spline  f  (if any) of order
%               k := length(knots) - length(x)
%  with knot sequence  knots  for which
%                 y = f(x) .
% This is taken in the osculatory sense in case some  x  are repeated, i.e.,
% in the sense that  D^m(i) f(x(i))=y(i)  in case the  x  are in
% nondecreasing order, with 
%   m(i) := #{j<i: x(j) = x(i)}.

% C. de Boor / latest change: December 24, 1989
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.

if (~isempty(find(diff(knots)<0))),
   error('The knot sequence should be nondecreasing!')
end

n=length(x);npk=length(knots);k=npk-n;
if (k<1),
   fprintf('number of data points, %.0f, should be less than ',n)
   fprintf('the number, %.0f, of knots!\n',npk)
   error('')
end

% It is assumed that  y  has the same sense as  x  which is important in case
% the entries of  y  are vectors rather than scalars
[rx,cx]=size(x); [ry,cy]=size(y);
if (rx>1),
   if (ry~=rx), error(' y  should contain exactly as many entries as  x ')
   else, y=y';end
else,
   if (cy~=cx), error(' y  should contain exactly as many entries as  x ')
   end
end

%  sort the given abscissae.
[x,index]=sort(x);y=y(:,index);

%  Generate the collocation matrix and divide it into the possibly reordered
%  sequence of given ordinates to generate the B-spline coefficients of the 
%  interpolant, then put it all together into  sp .

sp=spmak(knots,slvblk(spcol(knots,k,x,1),y')');
