Code covered by the BSD License

# LU factorization with complete pivoting.

### Nick Henderson (view profile)

14 Apr 2010 (Updated )

An implementation of algorithm 3.4.2 from Matrix Computations.

lucp(A,tol,pm_opt)
```function [L U p q] = lucp(A,tol,pm_opt)
%LUCP     LU factorization with complete pivoting.
%
% To compute the LU factorization under default settings:
%
%  [L U p q] = lucp(A)
%
% This produces a factorization such that L*U = A(p,q).  Vectors p and q
% permute the rows and columns, respectively.
%
% The pivot tolerance can be controlled:
%
%  [L U p q] = lucp(A,tol)
%
% The algorithm will terminate if the absolute value of the pivot is less
% than tol.
%
% Permutation matrices can be generated:
%
%  [L U P Q] = lucp(A,tol,'matrix')
%  [L U P Q] = lucp(A,tol,'sparse')
%
% The first will generate full permutation matrices P and Q such that
% L*U = P*A*Q.  The second generates sparse P and Q.
%
% If A is sparse, L and U will be sparse.
%
% This function works on non-square matrices.
%
% Input:
%  A = matrix
%  tol = pivot tolerance (defualt is 1e-10)
%  pm_opt = permutation output options
%         = 'vector' for permutation vectors, L*U = A(p,q), defualt
%         = 'matrix' for full permutation matrices, L*U = P*A*Q
%         = 'sparse' for sparse permutation matrices, L*U = P*A*Q
%
% Output:
%  L = lower triangular factor
%  U = upper triangular factor
%  p = row permutation
%  q = column permutation
%
% Reference:
%  Algorithm 3.4.2, Matrix Computations, Golub and Van Loan. (3rd ed)
%
% Other Implementations:
%  Gaussian Elimination using Complete Pivoting.  Alvaro Moraes.
%  http://www.mathworks.com/matlabcentral/fileexchange/25758
%
%  Gauss elimination with complete pivoting.  Nickolas Cheilakos.
%  http://www.mathworks.com/matlabcentral/fileexchange/13451
%  (Does not work with rectangular A)
%
%  Rank Revealing Code.  Leslie Foster.
%  http://www.math.sjsu.edu/~foster/rankrevealingcode.html
%  (Uses mex libraries for computation)
%

%
% 2010-03-28 (nwh) first version.
% 2010-04-14 (nwh) added more references.
% 2010-04-24 (nwh) perform final column swap so U is well conditioned
%

%
%

% handle optional inputs
if nargin < 2 || isempty(tol)
tol = 1e-10;
end

if nargin < 3 || isempty(pm_opt)
pm_opt = 'vector';
end

if strcmp(pm_opt,'vector')
pm_flag = false;
sp_flag = false;
elseif strcmp(pm_opt,'matrix')
pm_flag = true;
sp_flag = false;
elseif strcmp(pm_opt,'sparse')
pm_flag = true;
sp_flag = true;
else
error('lucp:invalidOption','''%s'' is an un recognized option.',pm_opt)
end

[n m] = size(A);

% pivot vectors
p = (1:n)';
q = (1:m)';

% temp storage
rt = zeros(m,1); % row temp
ct = zeros(n,1); % col temp
t = 0; % scalar temp

for k = 1:min(n-1,m)
% determine pivot
[cv ri] = max(abs(A(k:n,k:m)));
[rv ci] = max(cv);
rp = ri(ci)+k-1;
cp = ci+k-1;

% swap row
t = p(k);
p(k) = p(rp);
p(rp) = t;
rt = A(k,:);
A(k,:) = A(rp,:);
A(rp,:) = rt;

% swap col
t = q(k);
q(k) = q(cp);
q(cp) = t;
ct = A(:,k);
A(:,k) = A(:,cp);
A(:,cp) = ct;

if abs(A(k,k)) >= tol
rows = (k+1):n;
cols = (k+1):m;
A(rows,k) = A(rows,k)/A(k,k);
A(rows,cols) = A(rows,cols)-A(rows,k)*A(k,cols);
else
% stop factorization if pivot is too small
break
end
end

% final column swap if m > n
if m > n
% determine col pivot
[cv ci] = max(abs(A(n,n:m)));
cp = ci+n-1;

% swap col
t = q(n);
q(n) = q(cp);
q(cp) = t;
ct = A(:,n);
A(:,n) = A(:,cp);
A(:,cp) = ct;
end

% produce L and U matrices
% these are sparse if L and U are sparse
l = min(n,m);
L = tril(A(1:n,1:l),-1) + speye(n,l);
U = triu(A(1:l,1:m));

% produce sparse permutation matrices if desired
if pm_flag
p = sparse(1:n,p,1);
q = sparse(q,1:m,1);
end

% produce full permutation matrices if desired
if ~sp_flag
p = full(p);
q = full(q);
end
```