function Xvec = forwardSubstitutionC(LOWER,Bvec,n)
% Solving a lower triangular system LOWER Xvec = Bvec by forward-substitution.
% Input matrix LOWER is an n by n lower triangular matrix
% Input vector Bvec is n by 1
% Input scalar n specifies the dimensions of the arrays
% Output vector Xvec is the solution to the linear system
% This version uses column processing (good for Fortran).
% K. Ming Leung, 01/25/03, 02/16/05 (better variabe names)

Xvec = zeros(n,1);  % initialization for efficiency purpose only
for k = 1:n
    if (LOWER(k,k)==0) error('Matrix is singular!'); end;
    Xvec(k) = Bvec(k)/LOWER(k,k);
    % Column processing: involves the jth column of L 
    % below the main diagonal
    j = k+1:n;
    Bvec(j) = Bvec(j)-LOWER(j,k)*Xvec(k); % involves a scalar multiplying a vector
end       
