function [x,L,U,NROW]=gescplu(A,b);
%%% Gaussian Elimination with Scaled Column Pivoting
%%% and saving the LU decomposition and row permutation
%%Determine the size of the system
n=length(b);
x=zeros(n,1);
flag=0;

Scale=max(abs(A'))';%%This gives a vector of the maximums of each row
%%determine whether there is a row of zeros in coefficient matrix
if min(Scale)==0
        message='no unique solution'
        flag=1;
        else
        NROW=1:n;
end
A=[A,b];i=1;L=zeros(n);
while ( flag==0 & i<n);
   if i<n
         %%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
        message='no unique solution'
        flag=1;
     else  
        I=find(abs(A(NROW(i:n),i))./Scale(NROW(i:n))==M);
         %%(interchange rows if necessary)
            if I(1)~=1
              message='row interchange necessary'
              p=I(1)+i-1;
              NCOPY=NROW(i);
              NROW(i)=NROW(p);
              NROW(p)=NCOPY;
            end;
     end;
         for j=i+1:n
            L(NROW(j),i)=A(NROW(j),i)/A(NROW(i),i);
            A(NROW(j),i:n+1)=A(NROW(j),i:n+1)-L(NROW(j),i)*A(NROW(i),i:n+1);
         end;
     i=i+1;
   end;
   if A(NROW(n),n)==0
    message='no unique solution'
    flag=1;
   end;
end;
%%begin backward substitution
if flag==0;
  x(n,1)=A(NROW(n),n+1)/A(NROW(n),n);
   for i=n-1:-1:1;
     x(i,1)=(A(NROW(i),n+1)-sum(x(i+1:n,1)'.*A(NROW(i),i+1:n)))/A(NROW(i),i); 
   end;
end;
L=L(NROW(1:n),:)+eye(n);U=triu(A(NROW(1:n),1:n));
