function [L, U, P] = ludecompt(A)
% Factor square matrix into lower triangular matrix L and
% upper triangular matrix U using partial pivoting to reduce
% round-off error. P is the permutation matrix that reflects
% the row exchanges required by pivoting

% make sure A is a square matrix
[m, n] = size(A);
if m ~= n
	error('Matrix A must be square');
end
disp('Initial value of A');
A
% L must have 1's on its diagonal
disp('L at the start');
L = eye(n)
pause;
% Starting with I, we swap rows whenever pivoting swaps rows of A

disp('P at the start');
P = eye(n)
pause;

% Perform forward elimination to compute L, U, and P
for i = 1:n-1
	% Perform partial pivoting
	% Find the largest entry in absolute value over the index
	% range (i,i) to (n, i)
	[maxabs, k] = max(abs(A(i:n, i)));
	  if (maxabs == 0.0)
			continue;
	  end

	pivotindex = i + k -1;
	fprintf('Row %d is the pivot row\n', pivotindex);
	pause;
	if pivotindex ~= i
		% Swap rows so the largest element in magnitude is at A(i,i).
		% We don't have to exchange the zeros in columns
		% 1 through i-1 in each row
		tmp = A(i, i:n);	% Save A(i,i), A(i,i+1), ..., A(i,n)
		A(i, i:n) = A(pivotindex, i:n);
		A(pivotindex, i:n) = tmp;
		disp('Matrix after exchanging rows:');
		A
		pause
		% Swap the same rows in P.
		% We must swap whole rows
		tmp = P(i, 1:n);
		P(i, 1:n) = P(pivotindex, 1:n);
		P(pivotindex, 1:n) = tmp;
		disp('Permutation matrix after exchanging rows:');
		P
		pause
		% swap rows of L also, but only in columns 1 through i-1
		tmp = L(i, 1:i-1);
		L(i, 1:i-1) = L(pivotindex, 1:i-1);
		L(pivotindex, 1:i-1) = tmp;
		disp('Lower triangular matrix after exchanging rows:');
		L
		pause
	else
		disp('There is no need to exchange rows');
		pause;
	end

	% Now make all entries
	%       A(i+1,i), A(i+2,i), ..., A(n, i) 
	% zero. At the same time, store the factors used in L
	factors = A(i+1:n,i)/A(i,i);
	  % Use submatrix calculations instead of a loop to perform
	  % the row operations on the entris of rows i+1 to n to the
	  % right of the diagonal
	A(i+1:n, i+1:n) = A(i+1:n, i+1:n) - repmat(A(i, i+1:n), n-i, 1).*repmat(factors,1,n-i);
	  % The entries of column i, rows i+1 to n are zero 
	A(i+1:n,i) = zeros(n-i,1);
	disp('Matrix after applying the pivot element');
	A
	pause;
	  % Assign the factors to L
	L(i+1:n,i) = factors;
	disp('Lower triangular matrix after row operations:');
	L
	pause;
end

% A is the upper-triangular matrix we require
 magdet = abs(prod(diag(A)));
 if magdet == 0
	  fprintf('The matrix is singular to working precision');
 elseif magdet < 1.0e-6
	  fprintf('The matrix may be ill-conditioned\n');
 end

U = A;
disp('Final results:');
L
U
P
end
