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