function x = lusolve(U,P,b)
  %
  %  Given input upper triangular matrix U,
  %  a permutation matrix P and a vector b,
  %  returns the solution x to U*x==P*b  
  %
  Pb= permute(P,b);
  x= backsubs(U,Pb);
end

function Pb=permute(P,b)
  %
  %  Permutes the rows of b using the permutation matrix P.
  %  Warning: no checking that P is a permutation matrix is performed.  
  %
  %  Notice that no operations (add/multiply/divide) are performed here
  %
  n=length(b);
  Pb=zeros(n,1);
  for i=1:n
    for j=1:n
        if P(i,j)~=0
            break
        end
    end
    Pb(i)=b(j);
  end
end

function x=backsubs(U,y)
  %
  %  Solves the upper triangular system U*x=y and returns x. 
  %
  n=length(U);
  if length(y)~=n 
    error('The size of inputs U and y do not match')
  end
  for i = n:-1:1
    rhs=y(i);  % work with temporary variable rhs
    for j = (i+1):n
        rhs=rhs-U(i,j)*x(j);
    end
    x(i)=rhs/U(i,i);
  end
end
