function colloc = spcol(knots,k,tau,arg1,arg2)
%SPCOL B-spline collocation matrix.
%
%   SPCOL(KNOTS,K,TAU)  returns the matrix   (D^m(i)B_j(TAU_i)) , with B_j
%   the j-th B-spline of order K for the knot sequence KNOTS, TAU a
%   sequence of points, both assumed to be  n o n d e c r e a s i n g ,
%   and   m(i) := #{ j<i : TAU(j) = TAU(i) }.
%
%   If one of the optional arguments is a string with the same first two
%   characters as 'slvblk', the matrix is returned in the almost block-
%   diagonal format (specialized for splines) required  by SLVBLK.m .
%   (For backward compatibility, this also happens if the fourth argument is
%   present and is not a string).
%
%   If one of the optional arguments is a string with the same first two
%   characters as 'sparse', the matrix is returned as a sparse matrix.
%
%   If one of the optional arguments is a string with the same first character
%   as 'noderiv', then multiplicities are ignored, i.e., m(i) := 1 for all i.
%   (For backward compatibility, this also happens if the fifth argument is
%   present and is not a string).
%
%   The recurrence relations are used to generate the entries of the
%   matrix.
%
%   For example, the statement
%
%      x = linspace(t(1),t(end),101); c = spcol(t,3,x);
%
%   provides, in c(:,j), a fine sequence of values of the j-th quadratic 
%   B-spline for the given knot sequence t .
%
%   See also SLVBLK, SPARSE, SPAPI, SPAP2, BSPLINE.

%   Carl de Boor 13 may 89
%   cb : 19 apr 1993; correct output for case  n < 1
%   cb :  7 jun 1994; add option to ignore multiplicities
%   cb :  9 may '95 (use .' instead of ')
%   cb :  8 dec '95 (update help, add check on TAU and on M)
%   cb : 22 may '96 (make optional arguments more flexible)
%   cb : 11sep96 (add the option of returning a sparse matrix)
%   cb : 5 oct 97 (replace use of ANS, to help compilation)
%   cb : 19may98 (standardize the help)
%   cb : 28feb99 (use repmat; accept TAU as row or column)
%   cb : 24aug99 (correct mispelling in help)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc. 
%   $Revision: 1.12 $

if ~isempty(find(diff(knots)<0))
   error('The knot sequence KNOTS should be nondecreasing.')
end
if ~isempty(find(diff(tau)<0))
   error('The point sequence TAU should be nondecreasing.')
end

%  Compute the number  n  of B-splines of order K supported by the given
%  knot sequence and return empty matrix in case there aren't any.

npk=length(knots); n=npk-k;
if n<1 fprintf('There are no B-splines for the given input.\n')
   colloc = []; return
end

% Settle the options:
slvblk=0; noderiv=0;
if nargin>3
   if isstr(arg1)
      if     arg1(1)=='s', slvblk =1;
         if length(arg1)>1, if arg1(2)=='p', slvblk = 2; end, end
      elseif arg1(1)=='n', noderiv=1;
      else error(['The second optional argument should be ''sl'', ''sp'', ',...
      '''n'', or a number.'])
      end
   else slvblk=1;
   end
end
if nargin>4
   if isstr(arg2)
      if     arg2(1)=='s', slvblk =1;
         if length(arg2)>2, if arg2(2)=='p', slvblk = 2; end, end
      elseif arg2(1)=='n', noderiv=1;
      else error(['The second optional argument should be ''sl'', ''sp'', ',...
      '''n'', or a number.'])
      end
   else noderiv=1;
   end
end

% If  NODERIV==0, remove all multiplicities from TAU and generate repetitions
% of rows instead of rows containing values of successive derivatives.
nrows = length(tau); tau = reshape(tau,1,nrows);
if noderiv
   index = 1:nrows; m = ones(1,nrows); nd = 1; pts = tau;
else
   index = [1 find(diff(tau)>0)+1];
   m = diff([index nrows+1]); nd = max(m);
   if nd>k
      error(sprintf('Point multiplicity should not exceed the given order %g.',k));
   end
   pts = tau(index);
end

%set some abbreviations
km1 = k-1;

%  augment knot sequence to provide a K-fold knot at each end, in order to avoid
% struggles near ends of basic interval,  [KNOTS(1) .. KNOTS(npk)] .
% The resulting additional B-splines, if any, will NOT appear in the output.

[augknot,addl] = augknt(knots,k); naug = length(augknot)-k;
pts = pts(:); augknot = augknot(:);

%  For each  i , determine  savl(i)  so that  K <= savl(i) < naug+1 , and,
% within that restriction,
%        augknot(savl(i)) <= pts(i) < augknot(savl(i)+1) .

savl = max([sorted(augknot(1:naug),pts); k*ones(1,length(pts))]);

b = zeros(nrows,k);

% first do those without derivatives
index1 = find(m==1);
if ~isempty(index1)
   pt1s = pts(index1); savls = savl(index1); lpt1 = length(index1);
   % initialize the  b  array.
   lpt1s = index(index1); b(lpt1s,1) = ones(lpt1,1);

   % run the recurrence simultaneously for all  pt1(i) .
   for j=1:km1
      saved = zeros(lpt1,1);
      for r=1:j
         tr = augknot(savls+r)-pt1s;
         tl = pt1s-augknot(savls+r-j);
         term = b(lpt1s,r)./(tr+tl);
         b(lpt1s,r) = saved+tr.*term;
         saved = tl.*term;
      end
      b(lpt1s,j+1) = saved;
   end
end

% then do those with derivatives, if any:
if nd>1
   indexm=find(m>1);ptss=pts(indexm);savls=savl(indexm);lpts=length(indexm);
   % initialize the  bb  array.
   temp = [1 zeros(1,km1)]; bb = temp(ones(nd*lpts,1),:);
   lptss = nd*[1:lpts];

   % run the recurrence simultaneously for all  pts(i) .
   % First, bring it up to the intended level:
   for j=1:k-nd
      saved = zeros(lpts,1);
      for r=1:j
         tr = augknot(savls+r)-ptss;
         tl = ptss-augknot(savls+r-j);
         term = bb(lptss,r)./(tr+tl);
         bb(lptss,r) = saved+tr.*term;
         saved = tl.*term;
      end
      bb(lptss,j+1) = saved;
   end

   % save the B-spline values in successive blocks in  bb .

   for jj=1:nd-1
      j = k-nd+jj; saved = zeros(lpts,1); lptsn = lptss-1;
      for r=1:j
         tr = augknot(savls+r)-ptss;
         tl = ptss-augknot(savls+r-j);
         term = bb(lptss,r)./(tr+tl);
         bb(lptsn,r) = saved+tr.*term;
         saved = tl.*term;
      end
      bb(lptsn,j+1) = saved; lptss = lptsn;
   end

   % now use the fact that derivative values can be obtained by differencing:

   for jj=nd-1:-1:1
      j = k-jj;
      temp = [jj:nd-1].'*ones(1,lpts)+ones(nd-jj,1)*lptsn; lptss=temp(:);
      for r=j:-1:1
         temp = ones(nd-jj,1)*(augknot(savls+r)-augknot(savls+r-j)).'/j;
         bb(lptss,r) = -bb(lptss,r)./temp(:);
         bb(lptss,r+1) = bb(lptss,r+1) - bb(lptss,r);
      end
   end

   % finally, combine appropriately with  b  by interspersing the multiple
   % point conditions appropriately:
   dtau = diff([tau(1)-1 tau(:).' tau(nrows)+1]);
   index=find(min(dtau(2:nrows+1),dtau(1:nrows))==0); % Determines all rows
                                                    % involving multiple tau.
   dtau=diff(tau(index));index2=find(dtau>0)+1;     % We need to make sure to
   index3=[1 (dtau==0)];                            % skip unwanted derivs:
   if ~isempty(index2)
             index3(index2)=1+nd-m(indexm(1:length(indexm)-1));end
   b(index,:)=bb(cumsum(index3),:);

   % ... and appropriately enlarge  savl
   index = cumsum([1 (diff(tau)>0)]);
   savl = savl(index);
end

% Finally, zero out all rows of  b  corresponding to TAU outside the basic
% interval,  [knots(1) .. knots(npk)] .

index = find(tau<knots(1)|tau>knots(npk));
if ~isempty(index)
   temp = repmat([1-nd:0].',1,length(index))+repmat(nd*index(:).',nd,1);
   b(temp(:),:) = zeros(nd*length(index),k);
end

% The first B-spline of interest begins at KNOTS(1), i.e., at  augknot(1+addl)
% (since  augknot's  first knot has exact multiplicity K). If  addl<0 ,
% this refers to a nonexistent index and means that the first  -addl  columns
% of the collocation matrix should be trivial.  This we manage by setting
savl = savl+max(0,-addl);

if slvblk     % return the collocation matrix in almost block diagonal form.
              % For this, make the blocks out of the entries with the same
              %  SAVL(i) , with  LAST  computed from the differences.
   % There are two issues, the change in the B-splines considered because of
   % the use of  AUGKNOT  instead of  KNOTS , and the possible drop of B-splines
   % because the extreme  TAU  fail to involve the extreme knot intervals.

   % SAVL(j) is the index in  AUGKNOT  of the left knot for  TAU(j) , hence the
   % corresponding row involves  B-splines to index  savl(j) wrto augknot, i.e.,
   % B-splines to index  savl(j)-addl  wrto  KNOTS.
   % Those with negative index are removed by cutting out their columns (i.e.,
   % shifting the blocks in which they lie appropriately). Those with index
   % greater than  n  will be ignored because of  last .

   if addl>0   % if B-splines were added on the left, remove them now:
      width = km1+k;cc = zeros(nrows*width,1);
      index = min(k,savl-addl); 
      temp = +repmat(nrows*[0:km1],nrows,1);
    cc(repmat(([1-nrows:0]+nrows*index).',1,k)+repmat(nrows*[0:km1],nrows,1))=b;
      b(:)=cc(repmat([1-nrows:0].',1,k)+repmat(nrows*(k+[0:km1]),nrows,1));
      savl=savl+k-index;
   end
   ds=(diff(savl));
   index=[0 find(ds>0) nrows];
   rows=diff(index);
   nb=length(index)-1;
   last=ds(index(2:nb));
   if addl<0  nb=nb+1; rows=[0 rows]; last=[-addl last]; end
   addr=naug-n-addl;
   if addr<0  nb=nb+1; rows=[rows 0]; last=[last -addr]; end
   if slvblk==1
      colloc=[41 nb rows k last n-sum(last) b(:).'];
   else   % return the equivalent sparse matrix (cf BKBRK)
      nr = (1:nrows).'; nc = 1:k; nrnc = nrows*k;
      ncc = zeros(1,nrows); ncc(1+cumsum(rows(1:(nb-1)))) = last;
      ncc = reshape(cumsum(ncc),nrows,1);
      ijs = [reshape(nr(:,ones(1,k)),nrnc,1), ...
           reshape( ncc(:,ones(1,k))+nc(ones(nrows,1),:) , nrnc,1), ...
           reshape(b,nrnc,1)];
      index = find(ijs(:,2)>n);
      if ~isempty(index), ijs(index,:) = []; end
      colloc = sparse(ijs(:,1),ijs(:,2),ijs(:,3),nrows,n);
   end
else          % return the collocation matrix in standard matrix form
   width = max([n,naug])+km1+km1;
   cc = zeros(nrows*width,1);
   cc(repmat([1-nrows:0].',1,k)+ ...
              repmat(nrows*savl.',1,k)+repmat(nrows*[-km1:0],nrows,1))=b;
   % (This uses the fact that, for a column vector  v  and a matrix  A ,
   %  v(A)(i,j)=v(A(i,j)), all i,j.)
   colloc = reshape(cc(repmat([1-nrows:0].',1,n) + ...
                    repmat(nrows*(max(0,addl)+[1:n]),nrows,1)), nrows,n);
end
