function [L,U] = luscalar(A)
% function [L,U] = luscalar(A)
%
% lu factorization of A with fully scalar implementation

[n,n] = size(A);
for k=1:n-1
  for j=k+1:n
      A(j,k) = A(j,k)/A(k,k);
      for m = k+1:n
         A(j,m) = A(j,m) - A(j,k)*A(k,m);
      end
  end
end
U = triu(A);
L = eye(n)+tril(A,-1);