% LU Factorization
% Jason Cox, CS371, Spring 1996, Lab 7
%
% I selected this lab because Jason's was the most well-documented
% of the examples where someone 'vectorized' the inner loop of
% the factorization (see the section called 'forward elimination'.) --JMH

function [L,U,l] = lufact(A);
U = A;
L = zeros(size(A));

[m n] = size(A);
if n~=m
	error('Bad matrix dimensions');
end

% initialize pivot vector
l = [1:n];

% initialize scale vector
s = [];
for i=l,
	s = [s max(abs(U(i,:)))];
end

% loop through columns/rows
for k=1:n-1,
	% determine pivot row for this column
	rmax = 0;
	for i=k:n,
		r = abs(U(l(i), k)) / s(l(i));
		if r > rmax
			rmax = r;
			j = i;
		end
	end

	% place pivot row in proper order
	x = l(j);
	l(j) = l(k);
	l(k) = x;

	% forward elimination
	for i=k+1:n,  % go through the other remaining rows
		xmult = U(l(i), k) ./ U(x,k);
		U(l(i),:) = U(l(i),:) - xmult .* U(x,:);
		L(l(i),k) = xmult;
	end
end


% Build permutation matrix
i = 0;
p = zeros(n,n);
for j=l,
	i = i + 1;
	p(i,j) = 1;
end

% Switch around rows!
U = p * U;
L = p * L;
L = L + diag(ones(1,n),0);