function x = plu_solve ( P, L, U, b )
%
%  function x = plu_solve ( P, L, U, b )
%
%  Purpose:
%
%    PLU_SOLVE solves a linear system A*x=b, given the PLU factors of A.
%
%  Discussion:
%
%    A is assumed to be a square nonsingular matrix.
%
%    The PLU factors of A might come from the routine GE_PP, for instance.
%
%    The order of the system, N, is determined by the rows of B.
%    The number of right hand sides, NRHS, is determined by the columns of B.
%
%    It's best to think of the code in terms of the steps:
%
%      1: Solve P * z = b
%      2: Solve L * y = z.
%      3: Solve U * x = y.
%
%    where separate variables X, Y and Z are used.  However, it turns
%    out that we can actually use the same vector X for all three
%    variables X, Y and Z, so that's what we do here.
%
%  Modified:
%
%    22 February 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P(N,N), the permutation factor.
%
%    Input, real L(N,N), the unit lower triangular factor.
%
%    Input, real U(N,N), the upper triangular factor.
%
%    Input, real B(N,NRHS), the right hand side.  
%
%    Output, real X(N), the solution.
%
[ n, nrhs ] = size ( b );

x = b;
%
%  Step 1: Solve P * z = b
%
x = P' * x;
%
%  Step 2: Solve L * y = z.
%
for k = 1 : nrhs
  for j = 1 : n
    for i = j+1 : n
      x(i,k) = x(i,k) - L(i,j) * x(j,k);
    end
  end
end
%
%  Step 3: Solve U * x = y.
%
for k = 1 : nrhs
  for j = n : -1 : 1
    x(j,k) = x(j,k) / U(j,j);
    for i = 1 : j-1
      x(i,k) = x(i,k) - U(i,j) * x(j,k);
    end
  end
end
