# 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.

ldlt_symm(A, piv)
```function [L, D, P, rho, ncomp] = ldlt_symm(A, piv)
%LDLT_SYMM  Block LDL^T factorization for a symmetric indefinite matrix.
%           Given a Hermitian matrix A,
%           [L, D, P, RHO, NCOMP] = LDLT_SYMM(A, PIV) computes a permutation P,
%           a unit lower triangular L, and a real block diagonal D
%           with 1x1 and 2x2 diagonal blocks, such that  P*A*P' = L*D*L'.
%           PIV controls the pivoting strategy:
%             PIV = 'p': partial pivoting (Bunch and Kaufman),
%             PIV = 'r': rook pivoting (Ashcraft, Grimes and Lewis).
%           The default is partial pivoting.
%           RHO is the growth factor.
%           For PIV = 'r' only, NCOMP is the total number of comparisons.

%           References:
%           J. R. Bunch and L. Kaufman, Some stable methods for calculating
%              inertia and solving symmetric linear systems, Math. Comp.,
%              31(137):163-179, 1977.
%           C. Ashcraft, R. G. Grimes and J. G. Lewis, Accurate symmetric
%              indefinite linear equation solvers. SIAM J. Matrix Anal. Appl.,
%              20(2):513-561, 1998.
%           N. J. Higham, Accuracy and Stability of Numerical Algorithms,
%              Second edition, Society for Industrial and Applied Mathematics,
%              Philadelphia, PA, 2002; chap. 11.

%           This routine does not exploit symmetry and is not designed to be
%           efficient.

if ~isequal(triu(A,1)',tril(A,-1)), error('Must supply Hermitian matrix.'), end
if nargin < 2, piv = 'p'; end

n = length(A);
k = 1;
D = eye(n); L = eye(n);  if n == 1, D = A; end
pp = 1:n;
if nargout >= 4
maxA = norm(A(:), inf);
rho = maxA;
end
ncomp = 0;

alpha = (1 + sqrt(17))/8;

while k < n
[lambda, r] = max( abs(A(k+1:n,k)) );
r = r(1) + k;

if lambda > 0
swap = 0;
if abs(A(k,k)) >= alpha*lambda
s = 1;
else
if piv == 'p'
temp = A(k:n,r); temp(r-k+1) = 0;
sigma = norm(temp, inf);
if alpha*lambda^2 <= abs(A(k,k))*sigma
s = 1;
elseif abs(A(r,r)) >= alpha*sigma
swap = 1;
m1 = k; m2 = r;
s = 1;
else
swap = 1;
m1 = k+1; m2 = r;
s = 2;
end
if swap
A( [m1, m2], : ) = A( [m2, m1], : );
L( [m1, m2], : ) = L( [m2, m1], : );
A( :, [m1, m2] ) = A( :, [m2, m1] );
L( :, [m1, m2] ) = L( :, [m2, m1] );
pp( [m1, m2] ) = pp( [m2, m1] );
end
elseif piv == 'r'
j = k;
pivot = 0;
lambda_j = lambda;
while ~pivot
[temp, r] = max( abs(A(k:n,j)) );
ncomp = ncomp + n-k;
r = r(1) + k - 1;
temp = A(k:n,r); temp(r-k+1) = 0;
lambda_r = max( abs(temp) );
ncomp = ncomp + n-k;
if alpha*lambda_r <= abs(A(r,r))
pivot = 1;
s = 1;
A( [k, r], : ) = A( [r, k], : );
L( [k, r], : ) = L( [r, k], : );
A( :, [k, r] ) = A( :, [r, k] );
L( :, [k, r] ) = L( :, [r, k] );
pp( [k, r] ) = pp( [r, k] );
elseif lambda_j == lambda_r
pivot = 1;
s = 2;
A( [k, j], : ) = A( [j, k], : );
L( [k, j], : ) = L( [j, k], : );
A( :, [k, j] ) = A( :, [j, k] );
L( :, [k, j] ) = L( :, [j, k] );
pp( [k, j] ) = pp( [j, k] );
k1 = k+1;
A( [k1, r], : ) = A( [r, k1], : );
L( [k1, r], : ) = L( [r, k1], : );
A( :, [k1, r] ) = A( :, [r, k1] );
L( :, [k1, r] ) = L( :, [r, k1] );
pp( [k1, r] ) = pp( [r, k1] );
else
j = r;
lambda_j = lambda_r;
end
end
end
end

if s == 1

D(k,k) = A(k,k);
A(k+1:n,k) = A(k+1:n,k)/A(k,k);
L(k+1:n,k) = A(k+1:n,k);
i = k+1:n;
A(i,i) = A(i,i) - A(i,k) * A(k,i);
A(i,i) = 0.5 * (A(i,i) + A(i,i)');

elseif s == 2

E = A(k:k+1,k:k+1);
D(k:k+1,k:k+1) = E;
C = A(k+2:n,k:k+1);
temp = C/E;
L(k+2:n,k:k+1) = temp;
A(k+2:n,k+2:n) = A(k+2:n,k+2:n) - temp*C';
A(k+2:n,k+2:n) = 0.5 * (A(k+2:n,k+2:n) + A(k+2:n,k+2:n)');

end

% Ensure diagonal real (see LINPACK User's Guide, p. 5.17).
for i=k+s:n
A(i,i) = real(A(i,i));
end

if nargout >= 4 & k+s <= n
rho = max(rho, max(max(abs(A(k+s:n,k+s:n)))) );
end

else  % Nothing to do.

s = 1;
D(k,k) = A(k,k);

end

k = k + s;

if k == n
D(n,n) = A(n,n);
break
end

end

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