function x = solve(A,x,pvt)
% usage: x = solve(A,x,pvt)
%---------------------------------------------------------
% solution of linear system A*x = b
% A is returned by *decomp* 
% do not use if decomp has detected singularity
%
% inputs
%   A:   triangularized matrix obtained from *decomp*
%   x:   right hand side vector,
%        overwritten by the solution 
%   pvt: pivot vector obtained from *decomp*
%        optional, no permutations are applied if missing
% output
%   x:   solution
%
% Note: row version
%--------------------------------------------------------

n = length(x);

% solve L*y = P*b using forward elimination

% apply permutations if pvt is available
if (nargin>2)   % apply permutation
   for k=1:n-1
      if (k~=pvt(k))
         x([k,pvt(k)]) = x([pvt(k),k]);
      end
   end
end
%
for k=2:n
   % solve for y(k), which overwrites x(k)
   x(k) = x(k) - A(k, 1:k-1)*x(1:k-1);
end

% solve U*x = y using backward substitution
x(n) = x(n)/A(n,n);
for k=n-1:-1:1
   % solve for x(k)
   x(k) = x(k) - A(k, k+1:n)*x(k+1:n);
   x(k) = x(k)/A(k,k);
end