function [L,U]=LUbyRank1(A);
% LUBYRANK1 computes the LU factorization
%   [L,U]=LUbyRank1(A); computes the LU-factorization of A with
%   diagonal pivoting as n rank-one modifications. The implementation
%   here is for didactic purposes only.

n=max(size(A));
L=[]; U=[];
for j=1:n
  u=[zeros(j-1,1);1; A(j+1:n,j)/A(j,j)];
  v=[zeros(j-1,1); A(j,j:n)'];
  A=A-u*v';
  L=[L u]; U=[U;v'];
end