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 permutations
   for k=1:n-1
      if (k~=pvt(k))
         x([k,pvt(k)]) = x([pvt(k),k]);
      endif
   endfor
endif
#
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);
endfor

# 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);
endfor

endfunction