function [L,R] = lr(A)

% desc:	[L,R] = lr(A) stores an upper triangular matrix in R and a
%	lower triangular matrix in L, so
%	that A = L*R.  A must be square.
%	NO pivoting, error if any pivot is 0

n=size(A); n=n(1);	% A is n \times n
L=eye(n);  R=zeros(n); 
for k=1:n-1 		% loop over rows (without last)
  if A(k,k)==0 error(['0 pivot in row ',int2str(k)]); end
  for i=k+1:n 		% loop over columns
    L(i,k) = A(i,k)/A(k,k);
    A(i,k:n) = A(i,k:n) - A(k,k:n)*L(i,k);	% get new row
  end
end
R=A;
return;
