# The Matrix Computation Toolbox

### Nick Higham (view profile)

11 Sep 2002 (Updated )

A collection of M-files for carrying out various numerical linear algebra tasks.

gep(A, piv)
```function [L, U, P, Q, rho, ncomp] = gep(A, piv)
%GEP    Gaussian elimination with pivoting: none, complete, partial or rook.
%       [L, U, P, Q, RHO] = GEP(A, piv) computes the factorization P*A*Q = L*U
%       of the m-by-n matrix A, where m >= n,
%       where L is m-by-n unit lower triangular, U is n-by-n upper triangular,
%       and P and Q are permutation matrices.  RHO is the growth factor.
%       PIV controls the pivoting strategy:
%          PIV = 'c': complete pivoting,
%          PIV = 'p': partial pivoting,
%          PIV = 'r': rook pivoting.
%       The default is no pivoting (PIV = '').
%       For PIV = 'r' only, NCOMP is the total number of comparisons.
%
%       By itself, GEP(A) returns the final reduced matrix from the
%       elimination containing both L and U.

%       Reference:
%       N. J. Higham, Accuracy and Stability of Numerical Algorithms,
%       Second edition, Society for Industrial and Applied Mathematics,
%       Philadelphia, PA, 2002; chap. 9.

[m, n] = size(A);
if m < n, error('Matrix must be m-by-n with m >= n.'), end
if nargin < 2, piv = ''; end
pp = 1:m; qq = 1:n;
if nargout >= 5
maxA = norm(A(:), inf);
rho = maxA;
end
ncomp = 0;

for k = 1:min(m-1,n)

if findstr(piv, 'cpr')
if strcmp(piv, 'c')

% Find largest element in remaining square submatrix.
% Note: when tie for max, no guarantee which element is chosen.
[colmaxima, rowindices] = max( abs(A(k:m, k:n)) );
[biggest, colindex] = max(colmaxima);
row = rowindices(colindex)+k-1; col = colindex+k-1;

elseif strcmp(piv, 'p')

% Find largest element in k'th column.
[colmaxima, rowindices] = max( abs(A(k:m, k)) );
row = rowindices(1)+k-1; col = k;

elseif strcmp(piv, 'r')

% Find element that is largest in its row and its column.
col_last = k;
for it = 1:inf
[colmaxima, rowindices] = max( abs(A(k:m, col_last)) );
ncomp = ncomp + m-k;
row = rowindices(1)+k-1;
new_abs = abs(A(row,col_last));
if it > 1
if new_abs == last_abs
row = row_last;
break
end
end
last_abs = new_abs;
row_last = row;
[rowmaxima, colindices] = max( abs(A(row, k:n)) );
ncomp = ncomp + n-k;
col = colindices(1)+k-1;
new_abs = abs(A(row,col));
if new_abs == last_abs
col = col_last;
break
end
last_abs = new_abs;
col_last = col;
end

end

% Permute largest element into pivot position.
A( [k, row], : ) = A( [row, k], : );
A( :, [k, col] ) = A( :, [col, k] );
pp( [k, row] ) = pp( [row, k] ); qq( [k, col] ) = qq( [col, k] );
end

if A(k,k) == 0
if findstr(piv, 'c')
break
elseif strcmp(piv, '') % Zero pivot is problem only for no pivoting.
error('Elimination breaks down with zero pivot.  Quitting...')
end
end

i = k+1:m;
if A(k,k) ~= 0  % Simplest way to handle zero pivot for partial and rook.
A(i,k) = A(i,k)/A(k,k);          % Multipliers.
end

if k+1 <= n
% Elimination
j = k+1:n;
A(i,j) = A(i,j) - A(i,k) * A(k,j);
if nargout >= 5, rho = max( rho, max(max(abs(A(i,j)))) ); end
end

end

if nargout <= 1
L = A;
return
end

L = tril(A,-1) + eye(m,n);
U = triu(A);
U = U(1:n,:);

if nargout >= 3, P = eye(m); P = P(pp,:); end
if nargout >= 4, Q = eye(n); Q = Q(:,qq); end
if nargout >= 5, rho = rho/maxA; end
```