function x=slvblk(blokmat,b)
% SLVBLK	Solve an almost block-diagonal linear system.
%
%        x=slvblk(blokmat,b)
%
% returns the solution (if any) of the linear system  Ax=b, with the matrix 
%  A  stored in  blokmat  in the spline almost block diagonal form.
%
%  If the system is overdetermined (i.e., has more equations than unknowns),
% the least-squares solution is returned.

% C. de Boor / latest change: Jan.14, 1990
% C. de Boor / latest change: Nov.29, 1990 (correct signum problem)
% C. de Boor / latest change: Feb. 1, 1991 (correct signum problem further)
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.


% 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;

f=1;l=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;
   [-x(f-1+[elim+1:ncols],:);eye(bcol)];
   x(f:l,:)=blocks(f:l,1:elim)\(blocks(f:l,elim+1:ccols)*ans);
end
x=x(1:nu,:);
