function [L,D,U]=chap_2_sect_6_prob_25(A)
% slu Square ldu factorization with now row exchanges
%
% Written by:
% -- 
% John L. Weatherwax                2006-03-28
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

[n,n]=size(A);
tol=1.e-6; 
for k=1:n,
  if( abs(A(k,k))< tol )
    disp(['Small pivot in column', int2str(k)]); 
  end % cannot proceed without a row exchange
  L(k,k)=1;
  for i=k+1:n
    L(i,k) = A(i,k)/A(k,k); % the kth pivot is now A(k,k)
    for j=k+1:n
      A(i,j) = A(i,j) - L(i,k)*A(k,j); 
    end
  end
  for j=k:n
    U(k,j) = A(k,j); % rename eliminated A to call it U
  end
end
% factor out D from the matrix U:
for i=1:n,
  D(i,i)=U(i,i);
end
for i=1:n,
  for j=i:n,
    U(i,j)=U(i,j)/D(i,i);
  end
end