function lugausst(A, b) 
%LUGAUSS LU factorization using Gaussian elimination
%it is assumed that a zero never occurs on the diagonal during the
%process or row elimination
%[L,U] = lugauss(A) produces the LU factorization of a matrix A using 
%Gaussian elimination without pivoting : A = LU. L is a lower triangular
%and U is an upper triangular matrix. The solution to Ax = b is then
%obtained using forward followed by back substitution.
%Input  : Matrix A
%Output : Matrices L and U 

[m,n] = size(A);
if m ~= n
	disp('The matrix is not square.');
	return;
end

L = zeros(n,n);
for i = 1:n
	L(i,i) = 1;
end;

disp('Initial value of A');
A
disp('Initial value of L');
L
pause;

for i = 1:n-1
	if (A(i,i) == 0)
		disp('The algorithm has encountered a zero pivot');
		% make L and U empty matrices
		L=[];
		U=[];
		return;
	end;

	%replace the elements in column i, rows i+1 to n by
	%the multipliers that we will store in L
	A
	fprintf('The algorithm is at the diagonal element A(%d,%d)\n', i, i);
	fprintf('Zero-out entries A(%d,%d) to A(%d,%d)\n', i+1,i,n,i);
	disp('The multipliers that will be used are:');
	A(i+1:n,i)/ A(i,i)
	pause;
	fprintf('Store these multipliers in L at the same indices.\n');
	L(i+1:n,i) =  A(i+1:n,i)/ A(i,i)
	pause;
	fprintf('Modify the elements in rows %d to %d, columns %d to %d by subtracting multiples of row %d\n',...
			  i+1,n,i+1,n,i);
	A(i+1:n,i:n)  = A(i+1:n,i:n)  - L(i+1:n,i) * A(i,i:n);
	A
	pause;
end;

disp('L and U are');
L
U = triu(A)
pause;
disp('Forward substitution:');
y = zeros(n,1);
y(1) = b(1)/L(1,1);
fprintf('y(%d) = %f\n', 1, y(1));
pause;
for i=2:n
	y(i) = (b(i) - L(i, 1:i-1)*y(1:i-1))/L(i, i);
	fprintf('y(%d) = %f\n', i, y(i));
	pause;
end

disp('Back substitution');
x = zeros(n,1);
% Solve Ux = b by back substitution
x(n) = b(n)/U(n,n);
fprintf('x(%d) = %f\n', n, x(n));
pause;

for i = n-1:-1:1
	x(i) = (b(i) - U(i, i+1:n)*x(i+1:n))/U(i,i);
	fprintf('x(%d) = %f\n', i, x(i));
	pause;
end
