function [L,U] = LUNaive(A)
% LUNaive: Create an LU decomposition via naive Gauss elimination

% Although it is conceptually straightforward it takes 
% about 5 times more flops than strictly necessary

% This will differ from the MATLAB lu function in that it doesn't use pivoting
% Because of this the script will break if it encounters 
% a zero in a pivot position.
% input:
% A = coefficient matrix
% output:
% L = Lower Triangular left factor
% U = Upper Triangular right factor

% The first time through the loop we calculate U
% and the second time through we calculate L
% Essentially, the second loop calculates the inverse of the 
% elementary row operations used to find U.

[m,n] = size(A);
if m~=n, error('Matrix A must be square'); end

Aug = [A eye(size(A))];
% forward elimination to find U
for k = 1:n-1
for i = k+1:n

factor = Aug(i,k)/Aug(k,k);

Aug(i,k:2*n) = Aug(i,k:2*n)-factor*Aug(k,k:2*n);

end
end

U=Aug(:,1:n);

% forward elimination to find L

A=Aug(:,n+1:2*n);

Aug = [A eye(size(A))];
% forward elimination

for k = 1:n-1
for i = k+1:n

factor = Aug(i,k)/Aug(k,k);

Aug(i,k:2*n) = Aug(i,k:2*n)-factor*Aug(k,k:2*n);

end
end

L=Aug(:,n+1:2*n);

end