% Solves the PLU decomposed system PLUx = b. 
% If b contains several right hand sides, then each
% column of x represents the corresponding solution.
function x = plusolve(P,L,U,b)

[n,m] = size(b);

% Permute the rows of b appropriately.  This could be done
% by multiplying b by P transpose, but this involves
% unnecessary multiplies by zero values.
for j = 1:n
   bt(j,:) = b(find(P(:,j)==1),:);
end
b = bt;

% Forward substitute
y = zeros(n,m);
for k=1:n
   y(k,:) = (b(k,:) - L(k,1:k-1)*y(1:k-1,:))/L(k,k);
end

% Back substitute
x = zeros(n,m);
for k=n:-1:1
   x(k,:) = y(k,:)-U(k,k+1:n)*x(k+1:n,:);
end