function x = slvblk(blokmat,b,w)
%SLVBLK Solve almost block-diagonal linear system.
%
%   SLVBLK(BLOKMAT,B)  returns the solution (if any) of the linear system  
%   A*X=B, with the matrix A stored in BLOKMAT in the spline almost block
%   diagonal form (as generated, e.g., in SPCOL).
%
%   If the system is overdetermined (i.e., has more equations than
%   unknowns), the least-squares solution is returned.  In that case, it
%   may be useful to use the optional third argument, i.e., to use the
%   following:
%
%   SLVBLK(BLOKMAT,B,W)  returns the vector X that minimizes the 
%   w e i g h t e d  l_2 sum
%
%      sum_j W(j)*( (A*X-B)(j) )^2 .
%
%   The default for W is the sequence [1,1,1,...].
%
%   See also SPCOL, SPAPS, SPAPI, SPAP2.

%   Carl de Boor 14 jan 90
%   cb : Nov.29, 1990 (correct signum problem)
%   cb : Feb. 1, 1991 (correct signum problem further)
%   cb : 21 Jun, 1992 (change .. to full ellipsis ...)
%   cb : 15 sep, 1994 (add optional weight)
%   cb : 11 nov, 1994 (remove wasteful multiplication by identity)
%   cb : 13sep96 (add sparse matrix option)
%   cb : 10may98 (standardize the help)
%   cb : 13jul98 (initialize ELIM to avoid confusion with any ELIM.M)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc. 
%   $Revision: 1.9 $

% If BLOKMAT is sparse, handle the problem sparsely:
if issparse(blokmat)
   if nargin>2
      n = length(w); spw = sparse(1:n,1:n,sqrt(w));
      x = (spw*blokmat)\(spw*b);
   else
      x = blokmat\b;
   end
   return
end

% get the basic information
[nb,rows,ncols,last,blocks] = bkbrk(blokmat);

ne = sum(rows);nu = sum(last);
if any(cumsum(rows)<cumsum(last))|any(last>ncols)
   error('The coefficient matrix has a nontrivial nullspace.')
end

[brow,bcol] = size(b);
if(ne~=brow)
   error('Matrix and right side are incompatible.')
end

blocks = [blocks b];
ccols = ncols+bcol;
if nargin>2, w = sqrt(w); blocks = (w(:)*ones(1,ccols)).*blocks; end

f = 1; l = 0; elim = 0;
for j=1:nb
   if (f<=l) % shift the rows still remaining from previous block
      blocks(f:l,:) = ...
         [blocks(f:l,elim+1:ncols) zeros(l+1-f,elim),blocks(f:l,ncols+1:ccols)];
   end
   l = l+rows(j);

   elim = last(j);
   % ideally, one would now use
   %   [q,r] = qr(blocks(f:l,1:elim));
   % followed up by
   %   blocks(f:l,:) = q'*blocks(f:l,:);
   %   f = f+elim;
   % but, unfortunately, this generates the possibly very large square matrix q
   % The unhappy alternative is to do the elimination explicitly here, using
   % Householder reflections (and an additional inner loop):
   for k=1:elim
      a = norm(blocks(f:l,k));
      vv = abs(blocks(f,k))+a;
      c = vv*a;
      if blocks(f,k)<0, vv = -vv; end
      q = [vv;blocks(f+1:l,k)];
      blocks(f:l,:) = ...
       blocks(f:l,:)-((q/c)*ones(1,ccols)).*(ones(l+1-f,1)*(q'*blocks(f:l,:)));
      f = f+1;
   end
end

% now we are ready for back-substitution
x = zeros(f-elim-1+ncols,bcol);

for j=nb:-1:1
   elim = last(j); l = f-1; f = f-elim;
   % here is another occasion where empty matrices of various sizes would help;
   % instead, use an if statement:
   if elim<ncols, blocks(f:l,ncols+1:ccols) = blocks(f:l,ncols+1:ccols) ...
                    - blocks(f:l,elim+1:ncols)*x(f-1+[elim+1:ncols],:); end
   x(f:l,:) = blocks(f:l,1:elim) \ blocks(f:l,ncols+1:ccols);
end
x = x(1:nu,:);
