function knots = optknt(tau,k,maxiter)
%OPTKNT Optimal knot distribution.
%
%   OPTKNT(TAU,K)  returns the optimal knot sequence (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.
%   OPTKNT(TAU,K,MAXITER) permits specification of the number MAXITER of
%   iterations to be tried, the default being 10.
%
%   The interior knots are the  n-K  sign-changes in the absolutely constant
%   function  h  that satisfies
%
%          integral{ f(x)h(x) : TAU(1) < x < TAU(n) } = 0
%
%   for all splines  f  of order K with knot sequence TAU.
%
%   For example, Micchelli/Rivlin/Winograd and Gaffney/Powell would 
%   approve of the following way of constructing a spline interpolant
%   of order k to given data  x, y :
%
%      sp = spapi(optknt(x,k),x,y);
%
%   See also APTKNT, NEWKNT, AVEKNT.

%   Carl de Boor 23 mar 91
%   cb :  9 may '95 (use .' instead of ')
%   cb : 19 nov '95 (use fprintf for detailed error message)
%   matlab : 19mar96 (use sprintf and %d in error message)
%   cb : 11sep96 (compute (-1)^i via rem(i,2) and use sparse matrix)
%   cb : 5 oct 97 (replace use of ANS, to help compilation)
%   cb : 6dec97 (replicate without multiply; standardize the help)
%   cb : 13feb99 (use *damped* Newton to maintain SW conditions; make MAXITER
%        an optional third argument)
%   cb : 20feb00 (also provide appropriate end knots, obviating use of augknt)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.11 $

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

% Newton iteration, a la `Computational aspects of optimal recovery'
for iter = 1:maxiter
   ais = [signs*spcol(tauext,k+1,xi,'sp') -1/2];     % a_1,...,a_n  from (18)
   temp = cumsum(ais(n:-1:1).').'; cumais = temp(n:-1:1);% sum{a_j: j = i,...,n}
   % now solve for the changes in  xi  (see eqns (20) and (21))
   mmdflag = spparms('autommd'); spparms('autommd',0);    % suppress pivoting
   dx = signs.* ...
      ((-cumais(1:nmk).*(tau(k+[1:nmk])-tau([1:nmk]))/k)/spcol(tau,k,xi,'sp'));
   spparms('autommd',mmdflag);
   maxdx = max(abs(dx)); xi = xi+dx;
   % make sure the change does not destroy order nor the Schoenberg-Whitney
   % conditions:
   for checks=1:20
      if all(diff(xi)>0)&all(xi>tau(1:nmk))&all(xi<tau(k+1:n))
         checks=0; break, end
      dx = dx/2; xi = xi-dx;
   end
   if checks>0, error('Damped Newton iteration failed to converge.'), end
   if maxdx<tol, knots = augknt([tau(1) xi tau(end)],k); return, end
end
error(sprintf('Newton iteration failed to converge in %d steps.',maxiter));
