function x= lup_solv(L,U,P,b);
%%   function x= lup_solv(L,U,P,b);
%%   SOLVE THE SYSTEM(s) PLU=b  by solving Lz=Pb, Ux=z by forward and
%%   backward substition respectively. If only the coefficient matrix A
%%   is given, use [L,U,P]=lu_gescp(A) to obtain L, U and P.
%%
%%   INPUT:
%%     L - lower triangular matrix (from an LU decomposition of PA)
%%     U - upper triangular matrix (from an LU decomposition of PA)
%%     P - permutation matrix obtained for LU decomposition of PA
%%     b - right hand side of system(s), Ax=b, (PAx=Pb or PLUx=Pb)
%%         if only one system just a column vector
%%         if more than one right hand side, then a matrix with the b's
%%         as columns
%%   OUTPUT:
%%     x - solution as column vector, or, if more than one right hand
%%         side given, solutions as column vector in a matrix with same order
%%         of columns as b

b=P*b; % permute the rows of b
n=length(L(:,1));  % size of system
m=length(b(1,:));  % number of systems (right hand sides)

%% solve Lz=Pb by forward substitution;
for p=1:m;
     z(1,p)=b(1,p)/L(1,1);
     for i=2:n;
       z(i,p)=(b(i,p)-L(i,1:i-1)*z(1:i-1,p))/L(i,i);
     end;
end;

%% solve for Ux=z by backward substitution;
for p=1:m;
    x(n,p)=z(n,p)/U(n,n);
    for i=n-1:-1:1;
      x(i,p)=(z(i,p)-U(i,i+1:n)*x(i+1:n,p))/U(i,i);
    end;
end;
