function x = backsub(LR, b)
%BACKSUB Backsubstitution to solve triangular systems of equations.
%	X = BACKSUB(L,B) performs backsubstitution to solve
%	the equation L*X = B, where L is lower triangular.
%	BACKSUB is provided for pedagogical purposes; in most cases
%	it is faster to use X = L\B.

%	Copyright (C) 1988 the MathWorks, Inc.

% Check if any diagonal elements are essentially 0.
if any(abs(diag(LR)) < eps*norm(LR))
	error('Matrix is singular.')
end
[m,n] = size(LR);
if m ~= n
	error('Matrix must be square.')
end

% Find out if LR is L or R
lower = 0;
if norm(LR - tril(LR),1) == 0
	lower = 1;
elseif norm(LR - triu(LR),1) == 0
	lower = 0;
else
	error('Matrix must be triangular')
end

% Reorder rows of b and LR if system is upper triangular
if ~lower
	b = b(n:-1:1);
	LR = LR(n:-1:1,n:-1:1);
end

% initialize x to force it to be a column vector.
x = zeros(n,1);

% perform backsubstitution
x(1) = b(1) / LR(1,1);
for i = 2:n
	x(i) = (b(i) - LR(i,1:i-1) * x(1:i-1)) / LR(i,i);
end

% reorder answer if system is upper triangular
if ~lower
	x = x(n:-1:1);
end
