function [L, U, P, interchanges] = ludecomp(A)
%ludecomp LU factorization using Gaussian elimination with partial
%pivoting.
%   [P U P interchanges] = ludecomp(A) factors a square matrix so that
%   P*A = L*U. U is an upper triangular matrix, L is a lower
%   triangular matrix, and P is a permutation matrix that reflects
%   the row exchanges required by the partial pivoting used to reduce
%   round-off error. In the event it is useful, interchanges is the number
%   of interchanges required.

% Make sure A is a square matrix
[m, n] = size(A);
if m ~= n
	error('Matrix A must be square');
end

% L must have 1's on its diagonal
L = eye(n);
% Starting with I, we swap rows whenever pivoting swaps rows of A
P = eye(n);
interchanges = 0;
	
% Perform Gaussian 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 <= eps)
		continue;
	end

	pivotindex = i + k - 1;
	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;
		% 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;
		% 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;
		interchanges = interchanges + 1;
	end
	
   % compute the multipliers
	multipliers = A(i+1:n,i)/A(i,i);
	% Use submatrix calculations instead of a loop to perform
	% the row operations on the entries 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) - multipliers*A(i,i+1:n);
	% The entries of column i, rows i+1 to n are zero 
	A(i+1:n,i) = zeros(n-i,1);
	% Assign the multipliers to L
	L(i+1:n,i) = multipliers;
end

% A is the upper-triangular matrix we require
U = A;
