function x = solve( LU, p, b )
  %--------------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % February 28, 1999
  %
  %
  % Usage: x = solve( LU, p, b )
  %
  % Returns:
  %	x:	solution of Ax=b
  %
  % Parameters:
  %	LU:	The permuted LU factorization of a square matrix A
  %	p:	Permutation vector
  %	b:	The RHS vector of Ax=b
  %
  % THIS IS INEFFICIENT CODE - IT WAS WRITTEN FOR DEMONSTRATION PURPOSES.
  %
  % This function performs the backward and forward solves necessary to
  % use the LU factorization of a square matrix A to solve Ax=b.  The LU
  % factorization is assumed to have been generated by the function gauss()
  % and has an associated permutation vector that indicates the true row
  % and column values in the factored matrix.
  %--------------------------------------------------------------------------

  [n, m] = size( LU );
  if ( n ~= m )
    fprintf( 'ERROR: input matrix is not square\n' );
    return;
  end

  [bn, bm] = size( b );
  if ( bn ~= n || bm ~= 1 )
    fprintf( 'ERROR: input vector not %d x 1\n', n );
    return;
  end

  % Perform forward solve

  x = zeros( n, 1 );

  for k = 1 : n - 1
    for i = k + 1 : n
      b(p(i)) = b(p(i)) - LU(p(i),k) * b(p(k));
    end
  end

  % Perform backward solve

  x(n) = b(p(n)) / LU(p(n),n);
  for i = n - 1 : -1 : 1
    sum = b(p(i));
    for j = i + 1 : n
      sum = sum - LU(p(i),j) * x(j);
    end
    x(i) = sum / LU(p(i),i);
  end

end
