function [v,u]=lufac(A)

n=size(A,1);
% Pick diag(u) or diag(v) = ones
u=eye(n);
v=u;

% let u(1,1)*v(1,1) = A(1,1)
u(1,1)=A(1,1);

for j=2:n
   u(1,j)=A(1,j)/v(1,1);
   v(j,1)=A(j,1)/u(1,1);
end

for i=2:n-1
   u(i,i)=A(i,i);
   for k=1:i-1
      u(i,i)=u(i,i)-v(i,k)*u(k,i);
   end
   u(i,i)=u(i,i)/v(i,i);
   
   for j=i+1:n
      u(i,j)=A(i,j);
      for k=1:i-1
         u(i,j)=u(i,j)-v(i,k)*u(k,j);
      end
      u(i,j)=u(i,j)/v(i,i);
      
      v(j,i)=A(j,i);
      for k=1:i-1
         v(j,i)=v(j,i)-v(j,k)*u(k,i);
      end
      v(j,i)=v(j,i)/u(i,i);
   end
end

u(n,n)=A(n,n);
for k=1:n-1
   u(n,n)=u(n,n)-v(n,k)*u(k,n);
end
 



