function x = bs0(U,c)
% BS0  backsubstitution with no fancy stuff
% usage:   >> x = bs0(U,c)     % should be equivalent to  x = U\c

% FIXME:  unprotected divisions!
% FIXME:  no check on sizes
% FIXME:  check U is upper triangular

n = size(U,1);

x = zeros(n,1);

x(n) = c(n) / U(n,n);

for j = n-1:-1:1                % j is the row number
  sum = 0;
  for k = j+1:n
    sum = sum + U(j,k) * x(k);
  end
  x(j) = (c(j) - sum) / U(j,j); % U(j,j) is diagonal element in row j
end
