function xi = optknt(tau,k)
%OPTKNT	Optimal knots.
%
%        xi = optknt(tau,k)
%
% returns the (interior) optimal knots (according to the optimal recovery 
% theory of Micchelli/Rivlin/Winograd and Gaffney/Powell) for interpolation 
% at data points  tau(1), ..., tau(n)  by splines of order  k .
%
% These are the  n-k  sign-changes in the absolutely constant function  h  
% which satisfies  integral{ f(x)h(x) : tau(1) < x < tau(n) } = 0  for all
% splines  f  of order  k  with knot sequence  tau .

% C. de Boor: March 23, 1991
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.

n = length(tau);nmk=n-k;
signs=(-ones(1,nmk)).^[nmk-1:-1:0];
tauext=[tau, tau(n)*ones(1,k)];
% initial guess:
xi = aveknt(tau,k);  
tol = 1.e-7*(tau(n)-tau(1))/(nmk);maxiter=10;

% Newton iteration, a la `Computational aspects of optimal recovery'
for iter=1:maxiter,
   ais = [signs*spcol(tauext,k+1,xi) -1/2];     % a_1,...,a_n  from (18)
   cumsum(ais(n:-1:1)')'; cumais=ans(n:-1:1);      % sum{a_j: j=i,...,n}
   % now solve for the changes in  xi  (see eqns (20) and (21))
   dx = ...
  signs.*((-cumais(1:nmk).*(tau(k+[1:nmk])-tau([1:nmk]))/k)/spcol(tau,k,xi));
   % make sure that change does not destroy order
   xi = sort(xi + dx);
   if (max(abs(dx))<tol), return, end
end
error('Newton iteration failed to converge in %f steps\n',maxiter)

