On Tridiagonalizing and Diagonalizing Symmetric Matrices with Repeated Eigenvalues.
We describe a divide-and-conquer tridiagonalization approach for  matrices with repeated eigenvalues. Our algorithm hinges on the  fact that, under easily constructively verifiable conditions,   a symmetric matrix with band width  $b$ and $k$ distinct eigenvalues must be block diagonal with  diagonal blocks of size at most $b k$. A slight modification of  the usual orthogonal band-reduction algorithm allows us to reveal  this structure, which then leads to potential parallelism in the form of  independent diagonal blocks.  Compared to the usual Householder  reduction algorithm, the new approach exhibits improved data  locality, significantly more scope for parallelism, and the potential to  reduce arithmetic complexity by close to 50% for matrices that have only  two numerically distinct eigenvalues. The actual improvement   depends to a large extent on the number of distinct eigenvalues and  a good estimate thereof. However, at worst the algorithms behave like  a successive band-reduction approach to tridiagonalization.  Moreover, we provide   a numerically reliable and effective algorithm for computing  the eigenvalue decomposition of a   symmetric matrix with two numerically distinct eigenvalues.  Such matrices arise, for example, in invariant subspace decomposition  approaches to the symmetric eigenvalue problem.
Introduction
Let A be an n \Theta n symmetric matrix. Our goal is to compute an orthogonal-tridiagonal decomposition
of A, AQ=QT , where Q is orthogonal and T is tridiagonal. Reduction to tridiagonal form
is a standard preprocessing step in dense eigensolvers based on QR iteration, bisection, or Cuppen's
method [16]. The conventional tridiagonalization procedure [16, p. 419] reduces A one column at a
time through a Householder transformation at a cost of O(4n 3 =3) flops for the reduction of A, and an
This work was supported by the Advanced Research Projects Agency, under contract DM28E04120 and P-95006,
and the Mathematical, Information, and Computational Sciences Division subprogram of the Office of Computational
and Technology Research, U.S. Department of Energy, under Contract W-31-109-Eng-38.
y The work of this author was partially performed while she was a postdoctoral associate at Argonne National
Laboratory.
z All PRISM Working Notes can be retrieved via anonymous ftp from the pub/prism directory at ftp.super.org.
additional O(4n 3 =3) flops if the orthogonal matrix is accumulated at the same time. This algorithm
employs mainly matrix-vector multiplications and symmetric rank-one updates, which require more
memory references than the matrix-matrix operations [9,8,14].
The block tridiagonalization algorithm in [5, 15] combines sets of p successive symmetric rank-
updates into one symmetric rank-p update, at the cost of O(2pn 2 ) extra flops. As a result,
this algorithm exhibits improved data locality and hence is likely to be preferable on cache-based
architectures. This block algorithm has been incorporated into the LAPACK library of portable
linear algebra codes for high-performance architectures [1, 2]. Parallel versions for distributed-memory
machines of the standard algorithm and of the block algorithm are described in [12] and
in [13], respectively. A different approach to tridiagonalization is the so-called successive band
reduction (SBR) method, which completes the tridiagonal reduction through a sequence of band
reductions [10,7]. This approach leads to algorithms that exhibit an even greater degree of memory
locality, among other desirable features.
In this paper we show that if the number k (say) of distinct eigenvalues of a symmetric matrix
A is small, considerable scope exists for further savings in tridiagonalization algorithms. As will be
demonstrated, A can be cheaply reduced to a block diagonal banded form through a slightly modified
SBR approach. The final tridiagonal form is then achieved by applying the algorithm recursively
on the subblocks on the diagonal. Compared with the conventional approach, this approach has the
following advantages.
Improved data locality: The tridiagonalization process can employ mainly matrix-matrix oper-
ations, both in the reduction of A and in the update of the transformation matrix Q (see
also [10,7]).
Enhanced scope for parallelism: In the traditional algorithm, the scope for the exploitation
of parallelism in the reduction of A is limited to the application of the rank-1 update (for
the unblocked algorithm) or the rank-p update (for the blocked algorithm), and the scope
for parallelism decreases as subproblems become smaller. In contrast, our algorithm generates
independent subproblems during the reduction of A, which can be worked on independently and
whose number increases as the iteration proceeds. Thus, a shift occurs from data parallelism
(updates of large matrices) to functional parallelism (several independent subproblems), but
at any stage there is sufficient parallelism to exploit.
Reduced complexity: Depending on the number of distinct eigenvalues, we may almost halve the
number of floating-point operations. In addition, the need for data movement is reduced.
One particular situation where repeated eigenvalues arise is in the context of invariant-subspace
methods for eigenvalue problems [3,19,6,4], where a matrix with only two distinct, predetermined,
eigenvalues is generated either by repeated application of incomplete beta functions [19] or the
matrix sign function [4]. In exact arithmetic, our tridiagonalization procedure would result in a block
diagonal matrix with diagonal blocks of order no larger than 2. Hence the eigenvalue decomposition
could be computed easily by independently diagonalizing the 2 \Theta 2 blocks on the diagonal. In the
presence of roundoff errors, the computed tridiagonal matrix may not have this desirable structure.
However, we can prove that such a tridiagonal matrix can be diagonalized as reliably as with any
other method by two "cleanup sweeps," where each sweep solves at most n=2 independent 2 \Theta 2
eigenvalue problems.
The paper is organized as follows. We show in Section 2 that, under certain conditions that
can be verified easiliy, a banded symmetric matrix with bandwidth b and k distinct eigenvalues
is block diagonal with diagonal blocks of order at most bk. In Section 3, we present a reduction
algorithm to achieve the desired banded block-diagonal structure, through a slight modification of
the conventional band reduction procedure. This approach is then employed to develop a divide-and-
conquer tridiagonalization algorithm. An inexpensive algorithm for decoupling invariant subspaces of
matrices with eigenvalue clusters at 0 and 1 is given and verified in Section 4. Numerical experiments
with a Matlab implementation are reported in Section 5. Lastly, we summarize our results.
2 The Structure of Band Matrices with Repeated Eigenvalues
A tridiagonal matrix whose off-diagonal entries are all nonzero is called unreduced. It is well
known [18, p. 66] that an unreduced tridiagonal matrix does not have multiple eigenvalues. Conse-
quently, if an n\Thetan tridiagonal matrix has only k-n distinct eigenvalues, it must be block diagonal,
and the largest block cannot be larger than k\Thetak. The generalization of this fact to banded matrices
underpins the algorithm we propose, yet it is not as straightforward as it might seem.
Assuming that A is an n \Theta n symmetric matrix, we define the ith row bandwidth of A, denoted
by band row(i), as
That is, band row(i) is the distance of the first nonzero element in row i from the ith diagonal
element. Further, we say that A is nonincreasing in row bandwidth from b if
In particular, a banded matrix that is all zero below the bth subdiagonal, and all nonzero on the bth
subdiagonal is nonincreasing in row bandwidth from b.
With these definitions, we can now prove the following theorem.
Theorem 1 Let T be symmetric matrix with k distinct eigenvalues. If T is block diagonal, with
each diagonal block nonincreasing in bandwidth from at most b, the size of the largest block cannot
exceed kb.
Proof. Assume T has a diagonal block D of size p ? kb. By assumption, D is nonincreasing in
bandwidth from b; that is, D has p\Gammab rows with their first nonzero elements in different columns to
the left of the diagonal. Thus, for any -, rank(D \Gamma -I) is not less than
On the other hand, since p?kb and D has at most k distinct eigenvalues, D has an eigenvalue
- with multiplicity greater than b. Hence, rank(D \Gamma -I) is less than b. The contradiction verifies
the result of the theorem.
The following example shows the necessity of the "nonincreasing bandwidth" restriction in Theorem
1. Let
columns and
A=QQ T is symmetric with only 0 and 1 as eigenvalues. In fact,
\Theta \Theta \Theta \Theta 0 0
\Theta \Theta \Theta \Theta 0 0
\Theta \Theta \Theta 0 0 \Theta 0 0
\Theta \Theta 0 \Theta 0 \Theta 0 0
\Theta \Theta
\Theta \Theta
We see that A is banded with semi-bandwidth b= 3, but it is not block diagonal with blocks of size at
most since the "nonincreasing bandwidth condition" is violated by a(5; 2)=a(7; 4)= 0.
3 A Divide-and-Conquer Tridiagonalization Approach
The example in the preceding section showed that the standard Householder band reduction
algorithm will not necessarily reveal the block-diagonal structure. For example, if we had applied
the standard algorithm for reduction to bandwidth 3 to the matrix of example (3), the matrix would
have remained unchanged. Fortunately, a minor modification of the standard algorithm enforces
nonincreasing row-bandwidth, and hence the prerequisites of Theorem 1.
Let us consider the conventional reduction approach, where the matrix is reduced one column at
a time to semibandwidth b. In each reduction, the pivot row is always b rows below the diagonal,
whether the reduction of the previous column was skipped (i.e., the transformation was an identity)
or not. For example, if we reduce the matrix A in (3) to semibandwidth 3, row number 4 is the pivot
row for the reduction of the second column, and, since a(4 : 8; reduction is skipped. We
then proceed to column 3, using row 5 as pivot row, and the row-bandwidth increases. If, instead,
we employ a Householder transformation acting on a(4 : 8; 3) to eliminate a(5 : 8; 3), keeping row 4
as pivot row, we obtain
~
\Theta \Theta \Theta \Theta 0 0
\Theta \Theta \Theta \Theta 0 0
\Theta \Theta \Theta 0 \Theta 0 0 0
\Theta \Theta 0 \Theta \Theta 0 0 0
\Theta \Theta
\Theta \Theta
A is decoupled into two diagonal blocks of size at most 6 \Theta 6.
This example shows that nonincreasing bandwidth can be obtained easily if we do not increase
the pivot row when the previous reduction was skipped. For computational purposes, we define the
row bandwidth with respect to a threshold - :
band
That is, given a tolerance threshold - , a column a(i : n) is considered numerically zero if its 2-
norm is at most - . The Matlab function bred in Figure 3 shows the conventional bandreduction
algorithm augmented with (1) a threshold criterion for the generation of a Householder vector and
(2) a modified pivot row selection strategy, which does not change the pivot row if a transformation
is skipped.
The subroutines gen hh, pre hh, post hh, and sym hh generate a Householder vector and apply it
from the left, right, and symmetrically, respectively. Note that for simplicity the algorithm presented
here does not exploit the symmetry of A. However, if we wish to do so, we have sym hh work only
with a triangular part of A and omit the post hh (pre hh) call when working only with the lower
(upper) triangle. We also note that all the algorithms presented in this paper are available via
anonymous ftp from the pub/prism directory at ftp.super.org.
If no transformations are skipped, the procedure is identical to the conventional band reduction
procedure; otherwise, it may terminate earlier when the reduction reaches the last column of the
first diagonal block, and the problem is decoupled. Since we drop pivot columns whose norm is
O(- ), the decomposition will be accurate up to a residual of order - .
For simplicity we omitted an optimization in Figure 3. If the reduction of the first column of A
results in a bandwidth ~ b, say, where ~ b ! b, due to the small size of entries a( ~ b 1), we can
directly pursue a reduction of the trailing block to nonincreasing bandwidth ~ b, in the same fashion
as shown above.
If the parameter b is chosen such that kb ! n, where k is the number of distinct eigenvalues
of A, Theorem 1 predicts a decoupling of the problem, with the leading block being of size no
larger than kb. In particular, if b is chosen such that n=2, we can expect bred to generate
two decoupled subproblems of about the same size. We can then recursively divide the problem
until the transformed matrix becomes tridiagonal (i.e., 1). Figure 3 is a serial implementation
of tridiagonalization based on this approach. Note that the various subproblems can be dealt with
independently and simultaneously. The subroutine blk diag, which is called in tri sbr, is shown
in

Figure

3 and reduces a matrix to block diagonal form with a given bandwidth.
For example, if we reduce a 12 \Theta 12 matrix A with only two eigenvalues to bandwidth 3, no
diagonal block can be larger than 6\Theta6. Thus, if a(4; 1), a(5; 2), and a(6; are all nonzero after the
reductions in the first three columns have been completed, the next three columns must already be
reduced, and the (partially reduced) matrix A is of the
\Theta \Theta \Theta \Theta 0 0
\Theta \Theta \Theta \Theta \Theta 0 0
\Theta \Theta \Theta \Theta \Theta \Theta 0 0
\Theta \Theta \Theta \Theta \Theta \Theta 0 0
\Theta \Theta \Theta \Theta \Theta 0 0
\Theta \Theta 0 0
\Theta \Theta \Theta \Theta
\Theta \Theta \Theta \Theta
\Theta \Theta \Theta \Theta
\Theta \Theta \Theta \Theta
\Theta \Theta \Theta \Theta
\Theta \Theta \Theta \Theta
As a result, we do not need to perform the reductions that would otherwise have occurred in
columns 4 through 6. The complexity of the algorithm for the case for the
reduction of A, and O(1:25 n 3 ) for the update of Q, as compared with O(4n 3 =3) for both these
operations in the usual approach. The savings for Q are minor, since updates at later stages still
involve vectors of length n, whereas only diagonal subblocks are affected in A. In addition, we can
work in parallel on independent problems. If the estimate k of the number of distinct eigenvalues
is inaccurate, the algorithm becomes either the standard eigenvalue algorithm (for k ? n=2) or the
SBR tridiagonalization procedure suggested in [10], but in either case, it will return numerically
accurate results.
4 Invariant Subspace Splitting
The computational cost and the degree of parallelism in the algorithm depend on k, the number
of distinct eigenvalues. One particularly intriguing case is matrices that have only two eigenvalues,
which arise in eigensolvers based on variant subspace decompositions [3, 19, 4]. We may assume
without loss of generality that the eigenvalues are at 1 and 0 (any other two eigenvalues can be
mapped to 0 and 1 by shifting and scaling). The following corollary is a special case of Theorem 1.
Corollary 2 Let A be a matrix with two distinct eigenvalues, and let A=Q T TQ be a tridiagonal-
ization of A. Then T is block diagonal with diagonal blocks of size at most 2 \Theta 2.
Corollary 2 implies that one can determine the range space, R(A), and the null space, N (A), in
essence by tridiagonalizing A. Let AQ=QT be the orthogonal-tridiagonal decomposition of A. For
a 1 \Theta 1 diagonal block T (j; j),
Since the eigenvalues of A and T are the same, a 2 \Theta 2 diagonal block T (j:j +1; j:j +1) must have
eigenvalues 0 and 1. Since the trace is the sum of the eigenvalues, and the offdiagonal entry is
nonzero, we have
we conclude
\Gamma-
One can see that the separation of the range and null subspaces of A, and in fact its eigenvalue
decomposition, can be affected by diagonalizing (potentially in parallel) the 2 \Theta 2 subproblems still
occurring in the block tridiagonal decomposition.
In the presence of rounding errors, a computed tridiagonal matrix may, however, not exhibit the
block structure we could expect from Corollary 2, because of perturbations in the eigenvalues. That
-]g, and a repeated eigenvalue numerically manifests itself as an
eigenvalue cluster.
function [A, block1,
Given a symmetric matrix A, a bandwidth b, and a threshold tau, bred
% computes an orthogonal-banded matrix decomposition,
5 % where O(tau) denotes a matrix with a two-norm of order tau, and
% W is an orthogonal matrix.
% The output matrix A-output will be a 2x2 block diagonal matrix,
% where the first diagonal block A-output(1:block1,1:block1)
% is banded with bandwidth nonincreasing from b, and the second block
may be empty.
current pivot row
if (piv-row == n)
if (piv-row == j), break, end
row and column sets involved in current transformation
generate HH matrix to annihilate A(piv-row+1:n,j)
update jth row and column of A
A( rows,
% if the reduction is not "skipped", perform symmetric
update of A, update Q if required, and shift the pivot row
A(rows,
A(cols,
A( rows, rows
increase pivot row if A(piv-row,j) is not negligible
if (j ==
if (piv-row == j+1),
else
return; end

Figure

1: Nonincreasing Row-Bandwidth-Preserving Bandreduction Algorithm
function [A,
% produces an orthogonal-tridiagonal decomposition of
% a symmetric matrix A such that
% where A-new is tridiagonal and Q is orthogonal.
% The number k is a guess at the number of numerically distinct
% eigenvalues of A.
% Matrices are successively reduced to smaller bandwidth in an
% attempt to exploit the divide-and-conquer nature becoming
apparent in the successive bandreduction algorithm when the number
chosen is a good guess at the actual number of numerically distinct
eigenvalues.
if (block1 == n) % If problem didn't decouple, just reduce to
tridiagonal form
else
first subproblem is not tridiagonal yet
second subproblem is nontrivial
return;

Figure

2: Divide-and-Conquer Tridiagonalization
function
Given a symmetric matrix A, a desired bandwidth b, and a threshold tau,
produces an orthogonal-block-diagonal decomposition
% where O(tau) denotes a matrix whose norm is of order tay, and
% W is an orthogonal matrix.
will be a block diagonal matrix with each block banded with
nonincreasing bandwidth b. The i-th diagonal block starts
% at (blkvec(i), blkvec(i)).
If Q is not the empty matrix on input, Q is postmultiplied by W,
return; end

Figure

3: Reduction to Block Diagonal Form
Example 3 The matrix
e
ffl), has eigenvalues -(T
Hence, it seems that, for numerically relevant computations, we now would be faced with computing
the eigenvalue decomposition of a tridiagonal matrix. This is not the case, however. By
exploiting the special structure of the tridiagonal matrix, we can diagonalize it in two "sweeps" that
compute the eigendecomposition of all "even" or "odd" 2\Theta2 blocks on the diagonal (simultaneously),
respectively. As we show below, the fill-ins generated by these sweeps are of the same order as the
perturbation in the eigenvalues and hence can be considered negligible.
Lemma 4 Let T be a symmetric tridiagonal matrix with
Proof. Let Q be orthogonal and respectively, such that
I
Thus,
-.
The next lemma gives bounds on the elements of the Givens rotation we will choose to diagonalize
a 2 \Theta 2 block and minimize the size of fill-ins.
c s
\Gammas c
be a Givens rotation that diagonalizes a 2 \Theta 2 symmetric matrix
. Assume without loss of generality that fi ? 0, and define oe - 0 by
Then, s and c can be chosen such that
and 1
Proof. Let we wish to eliminate the off-diagonal elements in
G
G T , we obtain
If we choose
with oe as defined in (6), then
and
2oe
as claimed.
In the following theorem we now show that, employing these Givens rotations, we can limit
the size of the fill-in entries generated when applying these rotations to a tridiagonal matrix with
eigenvalue clusters around 0 and 1.
Theorem 6 Let T and - be as in Lemma 4. Let
c s
\Gammas c
I) be the Givens rotation
that diagonalizes one 2 \Theta 2 diagonal block of T , namely,
. fi
.   fl
* ~
where we assume fi ? 0 without loss of generality. If fi ?
7 - and c and s are chosen as suggested
in Lemma 5, then
Proof. Comparing corresponding entries in T 2 and T and invoking Lemma 4, we know that there
exist ffl, ffl, and ffl
-, such that
Using these identities, we have
and hence we can express oe 2 defined as in (6) as
Thus,
Now let - 1 be chosen such that fi ? -. Then
Equations (11) together with s - fi
2oe
imply
2oe
and
2oe
Using (12), we can easily show that -
2oe
- and hence the result of the theorem.
As a consequence of Theorem 6, we are able to compute the eigenvalue decomposition of a 2 \Theta 2
diagonal block in a tridiagonal matrix T with eigenvalue clusters at 0 and 1 such that the generated
fill-in is negligible compared with the eigenvalue perturbation. Thus, the diagonalization of T can
be done by two sweeps of (potentially concurrent) 2 \Theta 2 eigenvalue problems as shown in Figure 4.
In the first sweep, we diagonalize an "odd-even" 2 \Theta 2 problem if the off-diagonal entry is not too
small, and set the fill-in entries to zero, or otherwise just set the off-diagonal entry to zero to zero.
In the second sweep, we diagonalize the "even-odd" blocks. Since no more rotations follow, there is
no need to set the fill-in entries to zero.
Theorem 6 shows that the Frobenius norm of the fill-in matrix introduced by the algorithm
rr diagshown in Figure 4 is bounded by 3
n-, which is of the same order as the perturbation in
eigenvalues. The subroutine diag2, which is not shown here, computes the diagonalizing rotations
as outlined in Lemma 5. Hence, Algorithm rr diag is as numerically reliable as any other approach
for diagonalizing T , albeit much cheaper because it exploits the special structure of T .
5 Numerical Experiments
In this section we report on some numerical experiments with the algorithms presented in this
paper. All experiments were performed with Matlab Version 4.2a on a Sun Sparcstation iPX. For
function [Q,
Given a tridiagonal matrix A with eigenvalues 1 and 0, with
contained in [1-tau,1+tau] or [-tau,tau]
rr-diag computes an approximate eigendecomposition
diagonal 2x2 matrices
negligible fill-ins
diagonal 2x2 matrices
no more need to zero fill-ins
return; end

Figure

4: Diagonalization of a Tridiagonal Matrix with Eigenvalue Clusters at 0 and 1
readers wishing to experiment on their own, the Matlab files employed to generate these results can
be retrieved via anonymous ftp from the pub/prism directory at ftp.super.org.
First, we apply the bandreduction algorithm bred of Figure 3 recursively to the trailing subblock
of a 200 \Theta 200 matrix with two eigenvalue clusters of size 50 each at 1g. The radius
of each cluster is ffl 1:0e 3 , where ffl is the machine precision. The drop threshold tau in bred is set to
, and at each step the bandwidth is chosen so as to decouple the problem in the middle.
The succession of matrices generated is shown in Figure 5. The title of each picture shows the
current matrix size being worked on and the bandwidth to which it is to be reduced. At each
step, we compute the residual
current
We observe ffi - 7:2e \Gamma13 , which, given a machine precision consistent with our theory.
The same experiment, employing a matrix with 100 eigenvalues at 0 and 1 each, with the same
eigenvalue perturbation and drop threshold, is shown in Figure 6. Note that it is sufficient to reduce
the matrix to half the bandwidth chosen in Figure 5 to achieve decoupling. We observe ffi - 2:7e \Gamma13 .
We also note that in both cases, the first, third, and fourth splits occur at row (and column) 100,
176, and 188, respectively. The second split occurs at row 152 for Figure 5 and at row 150 for

Figure

6.
To test the behavior of our rank-revealing tridiagonalization (RRDG), we compare it with the
standard eigenvalue decomposition (EIG) and the QR factorization with column pivoting (QR). Our
test matrices are
1. tridiagonal matrices with eigenvalue clusters of radius p ffl generated by inserting random off-diagonal
perturbations of the order p p ffl in the matrix shown in Example 3, and
2. matrices generated by symmetrically multiplying the matrices from Example 3 with orthogonal
matrices generated via the QR factorization of a random matrix.
In the first case, we call rr diag shown in Figure 4, in the second case, we precede the call to
rr diag by a call to tri sbr as shown in Figure 3. The drop threshold for the divide-and-conquer
tridiagonalization is set to
7p ffl, which is the same threshold employed in the two final diagonalization
sweeps. For each of 100, we run 50 test cases each with matrix sizes 125, 250, and 375.
RRDG and EIG both compute an eigenvalue decomposition Q D, with D diagonal. We compute
~
round(D), that is, round each diagonal entry to the nearest integer, and we report both
the relative eigenvalue residual kQ T A \Gamma ~
n=2 as well as the relative orthogonality residual
n. Note that
n=2 is an estimate of kAkF . In the case of the QR factorization with
pivoting, which computes AP = QR for a permutation matrix P and an upper triangular matrix R,
we compute the rank
and ~
We then report
First Divide:
50 100 150 2004080120160200
Second Divide:
50 100 150 2004080120160200
Third Divide:
50 100 150 2004080120160200
Fourth Divide:
50 100 150 2004080120160200

Figure

5: Bandreduction Applied to Trailing Subblock of a 200 \Theta 200 Matrix with Four
First Divide:
50 100 150 2004080120160200
Second Divide:
50 100 150 2004080120160200
Third Divide:
50 100 150 2004080120160200
Fourth Divide:
50 100 150 2004080120160200

Figure

Bandreduction Applied to Trailing Subblock of a 200 \Theta 200 Matrix with Two

Table

1: Relative Residual in Subspace Splitting

Table

2: Relative Residual in Orthogonality
which should be small, since Q(1 : is a basis for the range space of A. For each case, we report
the worst residual.
We see that the divide-and-conquer tridiagonalization, followed by the two cleanup sweeps over
the resulting tridiagonal matrix, performs just as well as a full-fledged eigenvalue decomposition.
In both cases, the residual in the subspace splitting is of O(p ffl), as expected. The residual for QR
factorization does not include the perturbation at the eigenvalue 1, as do the other two approaches,
and therefore is smaller in all cases. In any case, the computed orthogonal matrices are orthogonal
up to machine precision. The Q computed by the eig function in Matlab is slightly less orthogonal,
since eig involves more transformations and as a result accumulates more rounding errors. Note
that all three approaches are worse for a full matrix in the case because the roundoff errors
in the orthogonal reductions are of the same order of machine precision. When p is bigger, the
roundoff errors are dominated by the perturbation in the eigenvalues, and hence RRDG and EIG
behave about the same for tridiagonal and full matrices.
6 Conclusions
This paper introduced an algorithm for reducing a symmetric matrix with repeated eigenvalues to
tridiagonal form. The algorithm progresses through a series of band reductions, each band reduction
stage forcing a decoupling of the band matrix into independent subblocks. Compared with the usual
Householder tridiagonalization procedure, this approach can save up to 50% of the floating-point
operations. We also developed a robust and inexpensive numerical procedure for diagonalizing the
resulting tridiagonal matrix in the case where the matrix has only two eigenvalue clusters around
zero and one. This case arises in eigenvalue decomposition algorithms based on invariant subspace
approaches. Taken together, these two algorithms allow for an efficient diagonalization of such
matrices.
The algorithm can be generalized immediately to the reduction of unsymmetric matrices to Hessenberg
form. The same irreducibility argument underlying Theorem 1 goes through for Hessenberg
matrices. We also note that in exact arithmetic, conjugate transposed eigenvalue pairs would end
up in the same block. However, since one triangle of a Hessenberg matrix is still full, the potential
for computational savings is greatly reduced.
Apart from its divide-and-conquer nature and the resulting potential for parallelism, as well as
its reduced operation count, our divide-and-conquer algorithm has another attractive feature. Since
our algorithm (at least in the early stages) reduces matrices to banded form with a relatively wide
band, it is easy to block the Householder transformations by using the WY representation [11] or
the compact WY representation [20], as has, for example, been described in [17]. In this fashion,
one can easily capitalize on the favorable memory transfer characteristics of block algorithms.



--R



A divide-and-conquer algorithm for the eigenproblem via complementary invariant subspace decomposition
Design of a parallel nonsymmetric toolbox
Structured second- and higher-order derivatives through univariate Taylor series
The PRISM project: Infrastructure and algorithms for parallel eigensolvers.
Parallel tridiagonalization through two-step band reduction
Fundamental Linear Algebra Computations on High-Performance Com- puters
A project for developing a linear algebra library for high-performance computers
A framework for band reduction and tridiagonalization of symmetric matrices.
The WY representation for products of Householder matrices.
A parallel Householder tridiagonalization stratagem using scattered square decomposition.
Reduction to condensed form for the eigenvalue problem on distributed-memory architectures
Evolution of Numerical Software for Dense Linear Algebra
Block reduction of matrices to condensed form for eigenvalue computations.
Matrix Computations.
Solution of large
The Theory of Matrices in Numerical Analysis.
A parallelizable eigensolver for real diagonalizable matrices with real eigenvalues.
A storage efficient WY representation for products of Householder transformations.
--TR
repeated eigenvalues;eigenvalue decomposition;tridiagonalization