function x = forwardbueler(A,b)
% FORWARDBUELER  Solve lower triangular system by forward substitution.
% Check the size of the inputs, and does checks before division by zero.
% Also checks that the input matrix A is lower triangular.
%
% Example:
%   >> A = tril(randn(3,3))    % lower triangular 3x3 matrix
%   >> b = randn(3,1)
%   >> x = forwardbueler(A,b)
%   >> A * x - b               % should be nearly zero

[n,m] = size(A);       % for any size of matrix, but must be square (mxm)
if m ~= n,  error('A must be square (n x n)'), end
if max(max(abs(triu(A,1))))>0, error('A must be lower triangular'), end

[p,q] = size(b);
if q ~= 1, error('b must be a column vector'), end
if p ~= n, error('b must be have same number of rows as A'), end

x = zeros(size(b));    % create  x  as a column vector like  b

if A(1,1) == 0.0, error('zero in A(1,1) position'), end
x(1) = b(1) / A(1,1);

for i = 2:n
  if A(i,i) == 0.0, error('zero in A(%d,%d) position',i,i), end
  % next line does dot product of i-1 values:
  x(i) = (b(i) - A(i,1:i-1) * x(1:i-1)) / A(i,i);
end
