function [X,Y] = lusolv(A,B,row)
% [X,Y] = lufact(A,B,row)
% A is the nxn matrix obtained from lufact, input.
% The matrix A must be in factored LU form,
% and row is the row permutation information, input.
% B is any nx1 vector, input.
% The solution to AX=B is X, output.
% The solution to LY=PB is Y, output.
% Forward substitution followed
% by back substitution is used.
[n n] = size(A);
Y = zeros(n,1);
Y(1) = B(row(1));
for k = 2:n,
Y(k) = B(row(k)) - A(row(k),1:k-1)*Y(1:k-1);
end
X = zeros(n,1);
X(n) = Y(n)/A(row(n),n);
for k = n-1:-1:1,
X(k) = (Y(k) - A(row(k),k+1:n)*X(k+1:n))/A(row(k),k);
end