function [L,U,P]=lu_gescp(A);
%%  function [L,U,P]=lu_gescp(A);
%%  
%%  LU FACTORIZATION OF THE MATRIX PA for some permutation matrix P,
%%     PA=LU where L is lower triangular, U is upper triangular and L has
%%     unit diagonal.  Accomplished using Guassian elimination and
%%     scaled column pivoting (P is obtained from the row interchanges)
%%
%%  INPUT:
%%    A   - square matrix to be factored 
%%
%%  OUTPUT:
%%    L  -  lower triangular matrix
%%    U  - upper triangular matrix
%%    P  - the permutation matrix so tha PA=LU
%%    Message that no decomposition exists

%%Determine the size of the system
n=length(A(:,1));

%%STEP 1
        NROW=1:n;
        flag=0;
        L=zeros(n,n); % start with L as zero matrix to be filled in
        U=zeros(n,n); % start with U as zero matrix to be filled in

Scale=max(abs(A'))';      %%This gives a vector of the maximums of each row
      
if min(Scale)==0 %% whether there is a row of zeros in coefficient matrix
        flag=1;
end

if flag==0;
   for i=1:n-1
         %%Determine the Max scaled quantity in the column on 
         %%and below the diagonal
     M=max(abs(A(NROW(i:n),i))./Scale(NROW(i:n)));
         %% If max is zero, then no solution
         %% otherwise, find first place where max is attained
     if M==0;
        flag=1;
     else  
        I=find(abs(A(NROW(i:n),i))./Scale(NROW(i:n))==M);
         %%(interchange rows if necessary)
            p=I(1)+i-1;
            if p>i
              NCOPY=NROW(i);
              NROW(i)=NROW(p);
              NROW(p)=NCOPY;
            end;
        L(NROW(i),i)=1; 
     end;
       if flag==0;
         for j=i+1:n
            L(NROW(j),i)=A(NROW(j),i)/A(NROW(i),i);
            A(NROW(j),i:n)=A(NROW(j),i:n)-L(NROW(j),i)*A(NROW(i),i:n);
         end;
         U(i,i:n)=A(NROW(i),i:n);
       end;
   end;
end;
if flag==0;
   if A(NROW(n),n)==0
    flag=1;
   end;
   L(NROW(n),n)=1;U(n,n)=A(NROW(n),n);
   L(1:n,:)=L(NROW(1:n),:); % rearrange rows of L
end;
if flag==0;
   Q=eye(n);
   P=Q(NROW,:);
end;
if flag==1
   U='A is singular';
end;
