function [LU,row,det1] = lufact(A)
% [LU,row,det1] = lufact(A)
% A is an nxn matrix, input.
% LU is an nxn matrix, output.
% row is the row permutation information, output.
% det1 is the determinant of A, output.
% Partial pivoting is used.
[n n] = size(A);
det1 = 1;
row = 1:n; % Initialize pointer vector.
for p = 1:n-1,
[max1 prow] = max(abs(A(p:n,p))); % Find the pivot row.
if p < prow+p-1,
temp = row(p);
row(p) = row(prow+p-1);
row(prow+p-1) = temp;
det1 = - det1;
end
det1 = det1*A(row(p),p);
if det1 == 0, break, end
for k = p+1:n,
mult = A(row(k),p)/A(row(p),p);
A(row(k),p) = mult;
A(row(k),p+1:n) = A(row(k),p+1:n) - mult*A(row(p),p+1:n);
end
end
det1 = det1*A(row(n),n);
LU = A;