T HE UN IV ER SIT Y OF MI CHI GAN COLLEGE OF LITERATURE, SCIENCE, AND THE ARTS Department of Mathematics Technical Report No. 3 AN ALGORITHM FOR GENERALIZED MATRIX EIGENVALUE PROBLEMS C. B. Moler G.. W. Stewart ORA Project 027750 under contract with: OFFICE OF NAVAL RESEARCH DEPARTMENT OF THE NAVY CONTRACT NO. NOOO14-67-A-0181-0023 WASHINGTON, D.C. administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR February 1972 Approved for public release; distribution unlimited.

AN ALGORITHM FOR GENERALIZED MATRIX EIGENVA LU E P ROB L'EMS by C. B. Moler Computer Science Department Stanford University Department of Mathematics University of Michigan and G. W. Stewart Departments of Computer Science and Mathematics University of Texas at Austin Abstract A new method, called the QZ algorithm, is presented for the solution of the matrix eigenvalue problem Ax = XBx with general square matrices A and B. Particular attention is paid to the degeneracies which result when B is singular. No inversions of B or its submatrices are used. TChe algorithm is a generalization of the QR algorithm, and reduces to it when B = I. Problems involving higher powers of X are also mentioned. The authors were supported in part by NSF grants at Stanford University and the University of Texas at Austin, and by an ONR contract at the University of Michigan. Reproduction in whole or in part is permitted for any purpose of the United States Government. ii

AN ALGORITHM FOR GENERALIZED MATRIX EIGENVALUE PROBLEMS C. B. Moler-/ G. W. Stewart*-~ 1. Introduct ion We shall be concerned with the matrix eigenvalue problem of determining the nontrivial solutions of the equation Ax = XBx, where A and B are real matrices of order n. When B is nonsingular this problem is formally equivalent to the usual eigenvalue problem B Ax - x. When B is singular, however, such a reduction is not possible, and in fact the characteristic polynomial det(A-LB) is of degree less than n, so that there is not a complete set of eigenvalues for the problem. In some cases the missing eigenvalues may be regarded as "'infinite". In other cases the entire problem may be poorly posed. The term infinite eigenvalue is justified by the fact that if B is perturbed slightly so that it is no longer singular, there may appear a namber of large eigenvalues that grow unboundedly as the perturbation is reduced to zero. However, if det(A-XB) vanishes identically, say when A and B have a common null space, then any X may be regarded as an eigenvalue. Such problems have unusually pathological features, and we refer to them as "1ill-disposed" problems. *-/ Computer Science Department, Stanford University, and Department of Mathematics, University of Michigan. *-/ Departments of Computer Science and Mathematics, University of Texas at Austin.

In numerical work the sharp distinction between singular and nonsingular matrices is blurred, and the pathological features associated with singular B carry over to the case of nearly singular B. The object of this paper is to describe an algorithrn for computing the eigenvalues and corresponding eigenvectors that is unaffected by nearly singular B. The algorithm, the heart of which we call the QZ-algorithm, is essentially an iterative method for computing the decomposition contained in the following theorem [10]. Theorem. There are unitary matrices Q and Z so that QAZ and QBZ are both upper triangular. We say that the eigenvalue problems QAZy = \QBZy and Ax =.Bx are unitarily equivalent. The two problems obviously have the same eigenvalues, and their eigenvectors are related by the equation x = Zy The algorithm proceeds in four stages. In the first, which is a generalization of the Householder reduction of a single matrix to Hessenberg form [4, 5], A is reduced to upper Hessenberg form and at the same time B is reduced to upper triangular form. In the second step, which is a generalization of the Francis implicit double shift QR algorithm [3,8], A is reduced to quasi-triangular form while the triangular form of B is maintained. In the third stage the quasi-triangular matrix is effectively reduced to triangular form and the eigenvalues extracted. In the fourth stage the eigenvectors are obtained from the triangular matrices and then transformed back into the original coordinate system. The transformations used in reducing A and B are applied in such a way that Wilkinson's general analysis of the roundoff errors in unitary

transformations [11] shows that the computed matrices are exactly unitarily equivalent to slightly perturbed matrices A+ E and B+ F. This means that the computed eigenvalues, which are -the ratios of the diagonal elements of the final matrices, are the exact eigenvalues of the perturbed problems (A+E)x = X(B+F)x. If an eigenvalue is well conditioned in the sense that it is insensitive to small perturbations in A and B (see [10] for a detailed analysis), then it will be computed accurately. This accuracy is independent of the singularity or nonsingularity of B The use of unitary transformations in the reduction also simplifies the problem of convergence: a quantity may be set to zero if a perturbation of the same size can be tolerated in the original matrix. Our computer program does not actually produce the eigenvalues X. but instead returns Ci and., the diagonal elements of the triangular matrices QAZ and QBZ. The divisions in Xi = ci/ji become the responsibility of the program's user. We emphasize this point because the ai and Bi contain more information than the eigenvalues themselves. Since our algorithm is an extension of the QR algorithm, the well known properties of the QR algorithm apply to describe the behavior of our algorithm. In their survey article [9], Peters and Wilkinson describe another approach for the case when B is nearly singular. In their method one computes an approximate null space for B and removes it from the problem. The technique is reapplied to the deflated problem, and so on until a well conditioned problem is obtained. The method. has the crucial drawback that

one must determine the rank of B. If a wrong decision is reached, the well-conditioned eigenvalues may be seriously affected. The special case where A is symmetric and B is positive definite has been extensively treated. For the case of well-conditioned B the "Cholesky-Wilkinson" method [ 6 ] enjoys a well deserved popularity. A modification of this algorithm for band matrices is given by Crawford II A variant of the Peters-Wilkinson method for nearly semidefinite B has been given by Fix and Heiberger [ 2]. Although our method does not preserve symmnetry and is consequently more time consuming than these algorithms, its stability may make it preferable when B is nearly semidefinite. Our algorithm can also be used to solve "X-matrix" problems of the form ( frCr + kr-lCr_1 + m.. + Co ) x = o by forming the generalized block companion matrices. For example, when r = 3, C2 C1 Co CB3 0 0 A = I 0 0 B =- 0o I 0 0 I assumed to be I Note that neither Cr nor CO is assumed to be nonsingular.

2. Reduction to Hessenberg-Triangular Form In this section we shall give an algorithm whereby A is reduced to upper Hessenberg form and simultaneously B is reduced to triangular form. While a treatment of the reductions in this and the following sectncrml. can be given in terms of standard plane rotations and elementary Hermitian matrices, we find it convenient from a computational point of view to work exclusively with a modified form of the elementary Hermitians. Accordingly, we introduce the following notation. By ~fr(k) we mean the class of symmetric, orthogonal matrices of the form I + vu T where v u = -2, v is a scalar multiple of u, only components k,k+l,...,k+r-l of u are nonzero, and uk = 1. Given any vector x, it is easy to choose a member Q of /r(k) so that Qx = x + (uTx)v has its k+l,...,k+r-l components equal to zero, its k-th component changed and all other components unchanged. Since uk = 1, the computation of Qy for any y requires only 2r-1 multiplications and 2r-1 additions. (In particular, use of a matrix in V2 requires only 3 multiplications instead of the 4 required by a standard plane rotation.) For the most part, we shall use only matrices in 12 and ~3 When a matrix Q in E3 (k) premultiplies a matrix A only rows k, k+l, and k+2 in QA are changed. If the elements k, k+l, and k+2 in a column of A are zero, they remain zero in QA. Likewise, if ZcE53(k), only columns k, k+l, and k+2 are changed in AZ. If some row has elements k, k+l, and k+2 zero, then they remain zero in AZ. Similar considerations hold for the class f2 '

All our transformations will be denoted by Q's and Z's with various subscripts. The Q's will always be premultpliers, that is row operat;ioiiS. 'P'e Z' s will always Lbe postmlultipliers, or columrrm operations. The letter Q is being used in its traditional role to denote orthogonal matrices. The letter Z was chosen to denote orthogonal matrices which introduce zeros in strategic locations. The first step in the reduction is to reduce B to upper triangular form by premultiplication by Householder reflections. The details of this reduction are well known (e.g. see [4,11]) and we confine ourselves to a brief description to illustrate our notation. At the k-th stage of the reduction (illustrated below for k = 3 and n = 5 ), the elements below the first k-l diagonal elements of B are zero. x x x x x O x x x x O 0 x x x 0 O X x x o O X1 x x Each x represents an arbitrary nonzero element. Each x represents an element to be annihilated in the next step. A matrix Q C k+(k) k n-k+l is chosen to annihilate b k+, kbk+2 k...,bn,k and B is overwritten by QkB, giving a matrix of the form illustrated below. x x x x x O x x x O O x x x 0 0 0 x x 0 0 0 x x This process is repeated until k = n-l. Of course A is overwritten by Qn-iQn-2..Q1A

After this reduction, A and B have the form A B x x x x x 0 x x x x x x x x x C 0 x x x x x x x x 0 0 0 x x x x x x x 0 0 0 0 x The problem now is reduce A to upper Hessenberg form while preserving the triangularity of B. This is done as follows (for k = 5 ) First QcE2(4) is determined to annihilate a1. The matrices QA and QB, which overwrite A and B then have the form x x x X X x x x x x x x x x x 0 x x x x x x x x x 0 0 x x. x x x x x x 0 0 0 x x 0 x x x x 0 0 0 x x The transformation has introduced a nonzero element on the (5,4)-position of B. However, a ZEc2(4) can be used to restore the zero without disturbing the zero introduced in A This step is typical of all the others. The elements of A are annihilated by Q's in the order illustrated below. x x x x x x x x x x x x x x x! 4 6 X X X X X x5x x

As each element of A is annihilated, it introduces a nonzero element on the subdiagonal of B, which is immediately annihilated by a suitably chosen Z. The entire algorithm, including the Householder triangularization of B may be summed up as follows. 1) For k = 1,2,...,n-l 1) Choose QC k+(k) to annihilate b +,kbk+2 bn k n-k+l k+ l,kk k+2,k' ' n,k 2) B QkB, A QkA 2) For k 1,2,...,n-2 1) For ~ = n-l,n-2,...,k+l 1) Choose Qk of2( ) to annihilate a~+l,k 2) A Qk~A, B Qk~B 3) Choose ZkQ ~CV2() to annihilate b,+l,~ 4) B ( BZkg, A - AZk The complete reduction requires about 17 n3 multiplications, 17 5n2 17 n3 additions and n square roots. If eigenvectors are also to be computed, the product of the Z's must be accumulated. This requires an additional n3 multiplications and 3 n3 additions. The product 2 t2 of the Q's is not required for the computation of eigenvectors.

In this and the next section we assume that A is upper Itessenberg and B is upper triangular. In this section we shall propose an iterative technique for reducing A to upper triangular form while maintaining the triangularity of B. The idea of our approach is to pretend that B is nonsingular and examine the standard QR algorithm for C = AB. The manipulations are then interpreted as unitary equivalences on A and B Specifically suppose that one step of the QR algorithm with shift or is applied to C. Then Q is determined as an orthogonal transformation such that the matrix (3.1) R = Q(C-olI) is upper triangular. The next iterate C' is defined as C' = RQT + I _ QCQT and is known to be upper Hessenberg. If we set A' = QAZ and B' = QBZ where Z is any unitary matrix, then A'B'-1 QAZZB Q =QAB-Q =C' The matrix Z can be chosen so that B' is upper triangular. Then, since A' = C' B' is the product of a Hessenberg and a triangular matrix, it is also upper Hessenberg. This insures that the nice distribution of zeros, introduced by the algorithm of section 2, is preserved by the QZ step. Thus a tentative form of our algorithm might read

1) Determine Q so that QC is upper triangular, 2) Determine Z so that QAZ is upper Hessenberg and QBZ is upper triangular, 3) A QAZ, B -QBZ. The problem is then to give algorithms for computing Q and Z which do not explicitly require C = AB The determination Q is relatively easy. For from (3.1) and the definition of C it follows that (3.2) Q(A-0rB) = RB - S Since R and B are upper triangular, so is S. Thus Q is the unitary matrix that reduces A-kB to upper triangular form. Since A-~B is upper Hessenberg, Q can be expressed in the form (3*3) Q Qn-lQn-2'' —1 where Qk EW2(k) To calculate Z we apply Q in its factored form (3.3) to B and determine Z in a factored form so that B stays upper triangular. Specifically Q1B has the form (n= 5) X X x X x 1 X x x x x O 0 x x x O O 0 x x 0 0 0 0 x If Q1B is postmultiplied by a suitable Z1ECV2(1) the nonzero element below the diagonal can be removed. Similarly Q2Q1BZ1 has the form 10

x x x x x o X X X X o) K X X X 0 0 0 x x 0 0 0 0 x and the offending nonzero element can be removed by a Z2cWf2(2) Proceeding in this way, we construct Z in the form Z Z1Z2... Zn_1 where ZkeCt2(k) Although QBZ is upper triangular, it is not at all clear that QAZ is upper Hessenberg. To see that it is, rewrite equation (3.2) in the form (3.4) QAZ = SZ + kQBZ From the particular form of Z and the fact that S is upper triangular, it follows that SZ is upper Hessenberg. Thus (3.4) expresses QAZ as the sum of an upper Hessenberg and an upper triangular matrix. In fact (3.4) represents a computationally convenient form for computing QAZ We summarize as follows. 1) Determine Q n-n-2= Qn_ (QkvC2(k)) so that S = Q(A-rB) is upper triangular. 2) Determine Z = Z1z2... Zn 1 (ZkeV2(k)) so that B' = QBZ is upper 'triangular 3) A' = SZ + rB' If this algorithm is applied iteratively with shifts cr!, '., there result sequences of matrices A1,A2,..., and B1,B2,.. satisfying 11

A+1 = QA, BA+1 = QBZ The matrices A, are upper Hessenberg and the B;, are upper triangular. Moreover, if B1 is nonsingular, then the matriee. -1 Cw = AVB, are the matrices which would have been obtained by applying the standard QR algorithm with shifts '1' -2'** to C1 = A1B As C tends to upper triangular form, so must Ay, since Bo1 is upper triangular. Most of the properties of the QR algorithm carry over to the QZ algorithm. The eigenvalues will tend to appear in descending order as one proceeds along the diagonal. The convergence of a(V) to zero n,n-1 may be accelerated by employing one of the conventional shifting strategies (v) Once a becomes negligible one can deflate the problem by working nn-l with the leading principal submatrices of order n-l. If some other subdiagonal element of AV, say a (,)-l becomes negligible, one can effect a further savings by working with rows and columns ~ through n Because we have used unitary transformations, an element of Av or By can be regarded as negligible if a perturbation of the same size as the element can be tolerated in A1 or B1 The algorithm given above is potentially unstable. If k is large compared with A and B, the formula (3.4) will involve subtractive cancellation and A' will be computed inaccurately. Since the shift k approximates the eigenvalue currently being found and the problem may have very large eigenvalues, there is a real possibility of encountering a large shift. Fortunately the large eigenvalues tend to be found last so that by the time a large shift emerges the small eigenvalues will have been computed stably. (The large eigenvalues are of course ill-conditioned 12

and cannot be computed accurately.) To be safe one might perform the first few iterations with a zero shift in order to give the larger eigenvalues a chance to percolate to the top. 13

4. Implicit Shifts The potential instability in the explicit algorithm results from the fact that we have used formula (3.4) rather than unitary equivalences to compute A'. One way out of this difficulty is to generalize the implicit shift method for the QR algorithm to the QZ algorithm so that both A' and B' are computed by unitary equivalences. The implicit shift technique has the additional advantage that it can be adapted to perform two shifts at a time. For real matrices this means that a double shift in which the shifts are conjugate pairs can be performed in real arithmetic. Since we are primarily interested in real matrices, we will concentrate on double shifts.' The method is based on the following observation. Suppose that A is upper Hessenberg and B is upper triangular and nonsingular. Then if Q and Z are unitary matrices such that QAZ is upper Hessenberg and QBZ is upper triangular, then Q is determined by its first row. In fact, AB-1 and QAB 1QH are both upper Hessenberg, so that, by the theorem on page 352 of [11], Q is determined by its first row. Thus we must do two things. First., find the first row of Q Second, determine Q and Z so that Q has the correct first row, QAZ is upper Hessenberg, and QBZ is upper triangular. The first part is relatively easy. The first row of Q is the first row that would be -1 obtained from a double shifted QR applied to AB. Since A is upper Hessenberg and B' upper triangular, it is easy to calculate the -1 first two columns of AB. But these, along with the shifts, completely determine the first row of Q. Only nonsingularity of the upper 2-by-2

submatrix of B is actually required here. If either bll or b22 is too small, so that this submatrix is nearly singular, a type of deflation can be carried out. We will return to this point later. The second part is a little more difficult, and is really the crux of the algorithm since it retains the Hessenberg and triangular forms. Only the first three elements of the first row of Q are nonzero. Thus, if Q1 is a matrix in ~3(1) with the same first row of Q, then Q1A and Q1B have the following form (when n = 6 ) x x x x x x x x x x x x 2 x x x x x x x x x x x x 1 1 x x x x x x x x x x x x O 0 x x x x 0 0 0 x x x O 0 0 x x x 0 0 0 0 x x O 0 0 0 x x 0 0 0 0 0 x As in the standard implicit shift QR algorithm, it is convenient to think of Q1 as the reflection which annihilates two of the three nonzero elements in a fictitious "zeroth" column of A We must reduce Q1A to upper Hessenberg and Q1B to upper triangular by unitary equivalences. However, we may not premultiply by anything which affects the first row. This is done as follows. The matrix Q1B has three nonzero elements outside the triangle. These can be annihilated by two Z's, a Zi in 135(1) which annihilates the (3,1) and (3,2) elements and then a Z'1 in X21) which annihilates the resulting (2,1) element. Let Z 1 ZZ" Then Q1BZ1 is upper triangular. Applying Z1 to Q1A gives 1AZ1 with the following form 15

x x X x X x X X X X X x x x x x x x x 0 x x x 0 0 0 0 x x This is multiplied by Q2 in 3(2) that annihilates the (3,1) and (4,1) elements. Then Q2Q1AZ1 and Q2Q1BZ1 have the forms X x x x x X X x x X X X x x x x x x 0 x x x x x 0 x x x x x 0 x x x x x O x x x x x o 1 1 x O O X x x 0 0 0 0 X X O O 0 0 x x 0 0 0 0 0 X The first columns are now in the desired form. The nonzero elements outside the desired structure have been "chased" into the lower 5-by-5 submatric es. Now, postmultiply by Z2, a product of a matrix in /3(2) and a matrix in M2(2) that reduces the current B to triangular form. Then premultiply by Q3 in?/3(3) to annihilate two elements outside the Hessenberg structure of the resulting A The process continues in a similar way, chasing the unwanted nonzero elements towards the lower, right-hand corners. It ends with a slightly simpler step which uses Q in 2(n-1) to annihilate the (n,n-2) element of the current A, thereby producing a Hessenberg matrix, and Zn-2 in V2(n-l1) which annihilates the (n,n-1) element of the current B, producing a triangular B but not destroying the Hessenberg A 16

The fictitious zeroth column of A is determined in part by the shifts. In analogy with the implicit double shift algorithm, we take the shifts o1 and a2 to be the two zeros of the two-by-two problem det(A-~$) = 0 where an-l n-l an-ln bn-l, n-1 n-l, n n,( I Bn-n an, n-l an,,n 0 bn, n It is not desirable to compute F1 and 02 explicitly, or even to find the coefficients in the quadratic polynomial det(A —B). Instead, following the techniques used in " hqr2 " [8 ], we obtain ratios of the three nonzero elements of the first column of (AB -F1I) (AB -P 2I) directly from formulas which involve only the differences of diagonal elements. This insures that small, but non-negligible, offdiagonal elements are not lost in the shift calculation. The formulas are (m = n-l) a a a a a a b b alo = [(Tmm 11 nn -11 - mmn nnm11 1 b -b 10 __b b b b () ( ( )( ) ( ) ] ( MnuM nn LI nn mm mm nn b11 21 b22 11 22 (I41) a aaa 21121 a a b 22 12 12 lnn 11 a2bll 22 b b b 20 2 11 (11 (22 bumm bl nn bi mm nn a32 a30 b22 17

We are now in a position to summarize the double implicit shift method. It is understood that A and B are to be overwritten by the transformed matrices as they are generated. 1) Compute alo, a20, and a30 by (4.1). 2) For k = l,2,...,n-2 a) Determine QkcV35(k) to annihilate ak+lk-l and ak+2,kb) Determine Zecl,3'(k) to annihilate bk+2 k+1 and b c) Determine Z"VC2(k) to annihilate bk+l k 3) Determine QnlE 2(n-1) to annihilate a 2 n, n-2 )4) Determine Znlaj2 (n-l) to annihilate b For each k, determination of Qk requires a few multiplications and one square root. Application of Qk to both A and B requires about lO0(n-k) multiplications. The work involved with each Z' is the same. Application of Z" requires only about 6(n-k) multiplications. The number of additions is about the same. Summing these for k from 1 to n-1 gives a total of about 131n multiplications, 115n additions and 3n square roots per double iteration. By way of comparison, for the double shift QR algorithm as implemented in " hqr" Z becomes simply k and Z" is not used. Furthermore, the k blk k transformations are carried out on only one matrix. Consequently, each double iteration requires about 5n2 multiplications, 5n2 additions and n square roots. Thus the QZ algorithm applied on two matrices can be expected to require roughly 2.6 times as much work per iteration as the QR algorithm on a single matrix. 18

In order to obtain eigenvectors, the Q's are ignored and the Z's 2 accumulated. This requires about 8n more multiplications and 8n more additions per double iteration. There is one difficulty. The formulas for a10, a20 and a< are not defined when bll and b2 are zero. Moreover, as bl and 11 22 1 b2 approach zero the terms that determine the shift (terms involving 22 ann n b,, etc.) become negligible compared to the other terms, so that the effect of the shift is felt only weakly. Part of the solution to this difficulty is to deflate from the top. If bl is negligible it may be set to zero to give the forms for A and B (n = 4) x x x x 0 x x x x x x x 0 x x x 0 x x x 0 0 x x 0 0 x x 0 0 0 x A Q in 9/2(1) can then be used to annihilate the (2,1) element of A which deflates the problem. The rest of the solution lies in recognizing that there is not much of a problem. If bll and b22 are small then the problem has large eigenvalues. We have already observed that the larger eigenvalues tend to emerge at the upper left, and the larger the eigenvalue, the swifter its emergence. Moreover the speed will not be affected by a small shift. This means thaet whenever the implicit shift is diluted by a small b or b22, the algorithm is none the less profitably employed in finding a large eigenvalue. 19

5. Further Reduction of the Quasi-Triangular Form The result of the algorithm described so far is in an upper triangular matrix B and a quasi-upper triangular matrix A in which no two consecutive subdiagonal elements are nonzero. This means that the original problem decomposes into one by one and two by two subproblerns. The eigenvalues of the one by one problems are the ratios of the corresponding diagonal elements of A and B. The eigenvalues of the two by two problems might be calculated as the roots of a quadratic equation, and may be complex even for real A and B There are two good reasons for not using the quadratic directly, but instead reducing the two by two problems. First, when A and B are real, the calculation of eigenvectors is greatly facilitated if all the real eigenvalues are contained in one by one problems. A more important second reason is that the one by one problems contain more information then the eigenvalues alone. For example, if all and bll are small then the eigenvalue k1 = all/bll is ill conditioned, however reasonable it may appear. This reason obviously applies to complex eigenvalues as well as real ones. Accordingly, we recommend that the two by two problems be reduced to one by one problems and that the diagonal elements, rather than the eigenvalues, be reported. Without loss of generality we may consider the problem of reducing two by two matrices A and B simultaneously to upper triangular form by unitary equivalences. For our purposes we may assume that B is upper triangular. Two special cases may be disposed of immediately. If bll is zero, then a Qc~2(1) may be chosen to reduce a21 to zero. The zero elements 20

of Q13 are notI; di.isturbed. S:milarl.y, if ,; zero, a Z,(. may be chosen to reduce a21 to zero without disturbing b21 In the general two-by-two case, it is not difficult to write down formulas for the elements of A' = QAZ and B' - QBZ for any Q and Z Moreover, these formulas can be arranged so that numerically one of a2 21 or b' is effectively zero. It is not obvious, however, that the other 21 element is numerically zero, and the effect of assuming that it is by setting it to zero could be disastrous. Consequently, we must consider a somewhat more complicated procedure. The theoretical procedure for reducing A to triangular form may be described as follows. Let X be an eigenvalue of the problem and form the matrix E = A-RB. Choose a ZcV2(1) to annihilate either ell or e21. Since the rows of E are parallel, it follows that whichever of ell or e21 is annihilated the other must also be annihilated. Now choose Qc&f2(1) so that either QAZ or QBZ is upper triangular. Since the first column of QEZ is zero and QEZ = QAZ - -QBZ, it follows that, however Q is chosen, both QAZ and QBZ must be upper triangular. In the presence of rounding error the method of computing. and the choice of Z and Q are critical to the stability of the process. A rigorous rounding error analysis will show that, under a reasonable assumption concerning the computed h, the process described below is stable. However, to avoid excessive detail, we only outline the analysis. We assume that all computations are done in floating point arithmetic with t base B digits and that the problem has been so scaled that underflows and overflows do not occur. We further assume that a21 is not negligible in the sense that a211 < -t/ AlA, where I1| denotes, say, the row sum norm. 21

The algorithm for computing k amounts to making an appropriate origin shift and computing an eigenvalue from the characteristic equation. It goes as follows. all = bl a12 = a12 - pb12 a22 = a22 - pb22 1 a22 b12a21 n - _ b 2 ~ 22 11 22 a21a12 q = b2 11 22 2 r =p + q (5-1) k = p + p + sign(p).17 (complex if r < ) We must now assume that the computed X satisfies the equation det(A' -kB') = 0, where IIA-A'/ _ < %A/A and!\B-B'l _<BI~ 1 with and B mall constants of order f. Define E' = A' - XB' and let E denote the computed value E = fl(A -XB) Then E' = E+H with IIH1! < c maxQJlAll, Ihl I|BI|} with U of order -t We claim that, approximately, (5.2) IEii > -t maxtlIAI1, Ij\ 1/BIl} 22

First we note t;lat (5. 9) >l till,e la > Pt by the assumption that a21 is significant. Now assume that llEf\ < B-% I I /B//. Then subtractive cancellation must occur in the computation of ell, el2, and e22. Thus all %,bll ~ a12 " %b12 and a22 \.b22. Hence we have I/A/ ~ Hl IIBII, and, from (5.3), iEl ~ H-tl liiBl a, a contradiction. Now 0 = det(E') = det(E)+ (el1+hll)h (e+h)h22 (+h2)h21+ hlle22- h12e21 Hence!det(E) I S 'pl1lEll max{(|A|, 1% IIBIII }+ P2[maxllAtl, X i IfBI }]2 where P1 and p2 are.of order -. From (5.2) it then follows that |det(E) I < p|lE|| max[/IA|I, |k I|BIt } -t where p is of order f. Now consider the determination of Z. Assume that the second row of E is larger than the first. Then ZcV2(1) is chosen to annihilate e2 Let F = EZ. Then f21 is essentially zero. Furthermore, since Z is unitary I f f22 -- det(E) I ~< p|l E'| maxAillA/, IxlflB/jV But If221 lie21l and, since e2 was assumed to be the larger row, lie211 = ||EII. Hence we have approximately 1fl1 < p max~lJA/l, Ill/\B1 I To choose Q, let C = AZ, D = BZ, 23

and let f.; l, c1,and dl be tlc fir,^st columns o i, C, and I Let Iq, denrlote vthe second row o-l' Q.:If lIA/ > |X| 11 1,We C W oosa (A to annihilate d21. Numerically this means that -t T where C is a constant of the order of f. We must show that q2c! is negligible. But IsT2 J- 1dT T |q ~l = l%2 fl + X qo dl| < jfl1il + 1 Xi l/qc dll < p maxQ[fAlA, III11BIJ + crlli BIIB < (P + o) IAli If, on the other hand, iXl\|B|\ > |A\\, we choose Q so that Iq TCl < cllAll~ It then follows that q2 dil - qT f-2 Cl / \xI < Plj- max(AA1/, IX!I\B/1} + Olx -l1iAi| < (p+ o)\B\ In summary, X is computed using (5.1), Z is chosen to annihilate the first element of the larger of the two rows of A-xB and Q is chosen to annihilate the (2,1) element of the smaller of the two matrices AZ and %BZ. In this way, we can be sure that the computed (2,1) elements of both QAZ and QBZ are negligible. In practice with matrices of any order, if the transformations are real, they are applied to the entire matrices. If the transformations are 21e4

complex, they are used to compute the diagonal elements that would result;, but are not actually applied. We thus obtain a quasi-triangular problem in which each two-by-two block is known to correspond to a pair of complex eigenvalues. The generalized eigenvectors of this reduced problem can be found by a back-substitution process which is a straightforward extension of the method used in " hqr2 " [ 8 ]. The vectors of the original problem are then found by applying the accumulated Z's 25

6. Some Numerical Results The entire process described above has been implemented in a Fortran program [ 7 ]. There are four main subroutines: the initial reduction to Hessenberg-triangular form, the iteration itself, the computation of the final diagonal elements, and the computation of the eigenvectors. The complete program contains about 600 Fortran statements, although this could be reduced somewhat at the expense of some clarity. The numerical properties observed experimentally are consistent with the use of unitary transformations. The eigenvalues are always found to whatever accuracy is justified by their condition. If an eigenvalue and eigenvector are not too "ill-disposed", then they produce a small relative residual. Similar numerical properties can not generally be expected from any algorithm which inverts B or any submatrix of B. This is even true of 2-by-2 submatrices, as illustrated by the following example due to Wilkinson..1.2 1 A B =.3.4 0 Here A is about the square root of the machine precision, that is, I is not negligible compared to 1, but 2 is. There is one eigenvalue near -2. Small relative changes in the elements of the matrices cause only small relative changes in this eigenvalue. The other eigenvalue becomes infinite as A approaches zero. Great care must be taken in solving this problem so that the mild instability of the one eigenvalue does not cause an inaccurate result for the other, stable eigenvalue.

Of course, the use of unitary transformations makes our tectlmique somewhat slower than others which might be considered. But the added cost is not very great. In testing our program, we solve problems of order 50 regularly. A few problems of orders greater than 100 have been run, but these become somewhat expensive when they are merely tests. One typical example of order 50 requires 45 seconds on Stanford's IBM 360 model 67. Of this, 13 seconds are spent in the initial reduction 29 seconds are used for the 61 double iterations required, and 3 seconds are needed for the diagonal elements and eigenvectors. If the eigenvectors are not needed and so the transformations not saved, the total time is reduced to 27 seconds. By way of comparison, formation of B A a la Peters and Wilkinson [9] and use of Fortran versions [12] of "orthes" [5] and " hqr2 " [8] requires a total of 27 seconds for this example. (All of these times are for code generated by the IBM Fortran IV compiler, H level, with the optimization parameter set to 2.) In the examples we have seen so far, the total number of double iterations required is usually about 1.2 or 1.3 times the order of the matrices. This figure is fairly constant, although it is not difficult to find examples which require many fewer or many more iterations. As a rule of thumb, for a matrix of order n the time required on the model 67 is about.36 n3 milliseconds if vectors are computed,.22 n3 milliseconds if they are not. 27

The example in Table 1 is not typical, but it does illustrate several interesting points. It was generated by applying non-orthogonal rank one modifications of the identity to direct sums of companion matrices. The companion matrices were chosen so that the resulting problem has three double roots, kl = 2 co 1 2 D5 5 2+ 2 i 4 = 6 - 26 2 The double root at o results from the fact that B has a double zero eigenvalue. All three roots are associated with quadratic elementary divisors; i.e., each root has only one corresponding eigenvector. The computed diagonals of the triangularized matrices are given in the table. Note that the four finite eigenvalues are obtained with a relative accuracy of about 10. This is about the square root of the machine precision and is the expected behavior for eigenvalues with quadratic elementary divisors. The singularity of B does not cause any further deterioration in their accuracy. Furthermore, the infinite eigenvalues are obtained from the reciprocals of quantities which are roughly the square root of the machine precision times the norm of B. Consequently we are somewhat justified if we claim to have computed the square root of infinity./ This prompts us to recall the limerick which introduces George Gamow's One, Two, Three, Infinity: There was a young fellow from Trinity Who tried vC But the number of digits Gave him such fidgets That he gave up Math for Divinity. 28

'?0 -o 50( -2( 6 6 16 5 5 5 -6 Do) -<'3 L -14 D 5 5 LO 5 ) -6 5 2'7 -I't 2'7 -17 5 5 5 16 -6 = B= 27 -28 38 -17 5 5 5 5 5 16 -6 5 27 -28 27 -17 16 5 5 5 5 5 -6 16 27 -28 27 -17 5 16 6 6 6 6 -5 6 25. 768670843143.2637605112.10-6 -12.821841071323.1312405807. 0-6 5.814535434181 + 10.071071345641 i 11.62907102873 0 5.800765071150 - 10.047220375909 i 11.601530302268 5.736511506410 + 9.935928843473 i 11.473022854605 5.510819468089 - 9.545122710676 i 11.021758784186 0.976972281.108 -0.976972290.108 0.49999999310489 + 0.86602543924271 i 0.49999999310489 - 0.86602543924271 i 0.50000000689511 + 0.86602536832617 i 0.50000000689511 - 0.86602536832617 i Table 1 29

A c lit iowl. ect';mnl t s Nlost. o:' Mof:ll:e't s work was clone dlnrinr5 a visit lo Stanfordl nivcr- LtyI whlere lhe received support from the Computer Science Department, the Stanford Linear Accelerator Center and National Science Foundation grant GJ-1158. Partial support was also obtained at the University of Michigan from the Office of Naval Research, contract NR-044-377. Stewart received support at the University of Texas from NSF grant GP-23655. W. Kahan and J. H. Wilkinson made several helpful comments. Linda Kaufman of Stanford has recently shown how to carry out the corresponding generalization of the LR algorithm and has written a program that accepts general complex matrices. 30

References [1] Charles Crawford, "The numerical solution of the generalized eigenvalue problem," Ph.D. Thesis, University of Michigan) 1970. Also to appear in Comm. ACM. [2] G. Fix and R. Heiberger, "An algorithnm for the ill-conditioned generalized eigenvalue problem," to appear in Numer. Math. [3] J. G. F. Francis, "The QR transformation -- a unitary analogu.e to the LR transforrration," Computer Journal 4, 265-271, 332-345 (1961-62). [4] A. S. Householder, "Unitary triangularization of a nonsymmetric matrix.' J. Assoc. Comput. Mach. 5, 339-34'2 (1958). [5] R. S. Martin and J. H. Wilkinson, "Similarity reduction of a general matrix to Hessenberg form," Numer. Math. 12, 349-368 (1968). [6] R. S. Martin and J. H. Wilkinson, "Reduction of the symmetric eigenproblem Ax = kBx and related problems to standard form," Numer. Math. 11, 99-110 (1968). [7] C. B. Moler and G. W. Stewart, "The QZ algorithm for Ax = XBx," submitted for publication. [8] G. Peters and J. H. Wilkinson, "Eigenvectors of real and complex matrices by LR and QR triangularizations," Numer. Math. 16, 181-204, (1970). [9] G. Peters and J. H. Wilkinson, " Ax = kBx and the generalized eigenproblem," SIAM J. Numer. Anal. 7, 479-492 (1970). [10] G. W. Stewart, "On the sensitivity of the eigenvalue problem Ax = XBx," submitted for publication. [11] J. H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford UniversityPress, 1965. [12] Fortran translations of the.Algol procedures in [5] ancd. r8] may be obtained from the NATS Project, Applied Mathematics Division, Argonne National Laboratories, Argonne, Illinois. 31

Lt&PLICI~ a' u.*,J (A-ai0- 6),li'al,l',i N A (NL.),iDj), t: (NU,.I)J,ALE' (N), AL L i(Nj, '.iA (i,),X (i, 0,,'b) LUG ICAL WA N1' C A 1) A i, it - ti- I- AL 4ATLiCES, SP dkieD IN AIuAY$ 4dli'lf NI LCUiwj. Ci PS IS TIE t LATivd O'LuCI'LfiJ OF' IL Li1ENi3 Ox' A ANJ o. C 'liiULS i PAI~aL OF.JALiahS, (ALL'A (M),lTA i(il)).oi 2'tih C uS'i'A (i) * A * - iA,'A (thi) * ~:, SiWGlULA. C '~t& JI GN;ALU"oS - Or' A4* - LA.ti)DA*134* CAL jg DJol Ai.iEJ b C iL.VIDING ALFi (At) Di L4i A (li) 1, Ei 21T 3E I,( ) L.i tf iu i;a hO. C IF (WAAILX) ALS) r'LNDS COi;RE ';dNEJu ~NiC U.. c: L S. OiLI UNj.L'AL(i i{ANOa[i/%A'LuS LU INVa... I &, C S3 ETiTdEr ii Uit J (Uii diUThl) ioAi Lid J ia,;ULa'i. C oE TA( ) IJ, i.LAL. C.rA (i1) z15 CJiM19LzA, kir AiL At4i~.LAAILNvut id Li-A'J 1i d an () Ii!i )..a1"i (lI) C LC-A,i( Phi5 UkCCJii wiTti ALi.A (MM)/ kti i A) V A ( ) /aLLAtl/A tfi+ 1 ) C CXO tM A.I;C CU N J GAL'S L V L Tl JUG i A L d4'"(i1) A N AL iLA( +) At 1; L L C i6 Ct: jAAitL O LJJ( ii i. U' ajj e, ijAL N1 / Ail' 'li.1C I t' A AD b lta.L ai.duC'D 'O 1i'iANGOLa. OtOLIa UitiliAh/ I iVA~L~JiN, C: I~ALkA A/J) B3iTA 'OJUL 'i: Iilt UAGCNit. jL:. it AiaD Ai!o ACIUAILLl Yih/UCLD ONLY Ti,) UA;I-T''iaAadLJAi( i L-, a Wi Thi C 1-d3Y- I AND 2-di-2 DLUC;KS U)N DiaGONiAL Ct A. C lr' A LFA (ei) IS. 4UT i/AL, Ti,-itiN 6EBTA() I. NUI tO. 1T 1l 'ii 1 litvUUbLI; riA1ZA1Ua AtL) Clndhe' Cbu N rK C Lt ({'r.TR (ij.i:..) L lVi(f'IHlNG IS Ctk. i'ii- (ci) ~i Ndni Ok' i~-AEhA'~LN3. iiD'Di FO( i L-l'ia:iCL;NVALU.* _ IL' (' {l i(1) 1Td ITn1:'i(l'i). r. -1) 1~a 1 1 i'E; lA' LtAiJ i3Ua f1-:itii '.'.L, -l,lV.LUa J.Lv. iOT,C;NVh iiGE AND ALEA (1) '~it U AF Ae) AN 1 C u1. (l) -lATd LU EIA ( ii) AR: PaOIABL~ l[.acddli',. 'i i' (,AiaL&) A i.,id) iS TIm, li-rTA dEAL:iOLNW'kCi'ud. C( A. (-,!) AID j,1+1) Li Ttki ilALt AND l' I iNt/ii Ft'itL C )' ~,ia iLi LlL PkLA:iGEVECTOFi. AC.,.1) ANo -A (., a i ) aNi) IlI i RL:AIL AND) I AliJ aAL.i tiaIl'S C OJF TI'i (i + 1)-3S LEU', LA El V lCiU. VLCJL i4J)f.1A~AAlii.LL 50 SU 11A~,.Ai(gS r C:CitCN.''NT 1. 1r 1.+).1L Ci USL f4Jia PlldALiY bUbSuu'INEIS, QZtiES, Z'T, IZVAL AND,iVLt. C U j ktuUd ALA iALi.LIAX SUBOU'I'NS, l1tn3, hSt'ii, CaiSt-2 aNL) C).IV. C ~j, ro5 Sr NvDAtUIu ktULLONS, UShT' AilD ) bAu5. C A Uliui'(S: C. A UiLAu, STANF'OaD, ANJ) i. U. S 1'tWaIth, U. ut' TEAAS C- 'it'iL V'.JiON L4TrD) 1/ t/72. C 'ALL, lS (, L, h, A,b, AN'L X,X) CLL QZIT (N.D,d,A,B,L~$, PSA,EPS~, iL ER, iAJTk,A) Ca2,..3L V AL (ND,;J,A,,EPS;A, ~2$SBA LF Aii{,Ai' lEiA, NANrX,A) iL; iANTX) CALL,~V (mD,&,A,U,rPSA.,/Af,,AL.E.Ls]AX).\ iU tL ~NI)

N PNT o. F:) o 6 (( - (f T) P N'lI=T 06 O0 )' c- = (7'!") ni!nT, l I IT,. (. L '" Tr) q*j. + (f' T) ' - (Pr") v '1=T OL 0C or~.1/,T- = T,q fin T RTt: 03 Oq (r'T)Y, "'TT)cT + 3, = 3.'! 7-T (C0 0 * (r'T);1 1t:,'7=T Ot 00?fNT;Trnn C 'i * + ( 1 = J = 7 =T (I )) On T.' =., '. t= r T- Oh ('lI) [ + (7l! " = (T.) -' - 'l (r 0 dT *C = r ((7'1)FT) ~Cl"f!' =P n r 001 oI 0o (OPn (-r"s) AT 7*' TqO T) cq + = c;I"'T=T C; 7O l',q + - = OOIt ( " TT ( ' 1f - N T'[! P.'N 7 HPiT' T=T. p rr2 c = (r'0 = t't= 'r 7' 0( NLT r 0(1 t -J F Y:aN,. yl ).ST VT)',T )' T l V T t P.' 1JT T ('TTI X X Pt PtM M v iJ T e l (n s ne1X',O ij n. ( TNT ' I '') n!yT e; j, ~tV~ ~ ~ ~ (-')]'IJI]~I:J''J',.V. OJ lq' ~y9r~ 1.,:t?]lTIJN

100 CONTINUE C aMDUCE A TO UPPER fE SSENBEB G, KEEP b TBIANGULAR C r (Pi.El-) -;O lu 170 DO 16'0 'K1#,I2 [.:1:&.+ t.X1 = N-K-1 DO 150 LD=1,K[1 L = a-LB LI 5 L+1 CALL USH2(A(L.,K),A(L1, K)U1,U2,V1,V2) IF (U1.E.l1.) GO TO 125 DO 110 JK,em T 3 AIL,J) + 02*A(1,.J) A(L,J) = A(If,J) + T*V1 A(L1,J) = A(L1,J) + T*V2 10 CONTIUUE A(L1,K) = 0. D00 120 J2L,* T = B(L,J) + Ue*8(LlJ) B(LJ) = B (L,J) + T*V1 BL1,iJ) = iL1,J) + T*V, 120 Ct04l'I.NUE i2'..q. CA1L, i &- (B (lJlrl. 11) #i; (L;1, L) Odf IJ J v' ) i i (.UI. 1. 1.Jt; Tv 15i0 DO 1 1= Iil T: (iL1) + U1:2 (.,L) (i,,1) = U(s (I,L1) + T*V1 B(II,L): U(IL) + T*V2 130 COSt i'N U Li (iL, = 0. DO 143 I1l,. T = A(t,L1) + U2*A(,L,*) A(I,L'ItL: AtI,L1) + T*V I A.('I,) - A (],L) 'T "eV2 140 CONTINUE,fk' I* O'. WANTX) GO TU 150 DO 145 I= 1, T i -X(,,IL1 + U2*k (i,L) -KlX(I.J1) A [IK, L1) + T*V1, I(.L1):= ''X (I,L) + T*V2 1:4S -COlTIlUE 1.50 coolixut 150 CO#L;MUa 170 RaetSiN MND

S UBROdT1 J QZOi.IT ( D,,A, Bv EPS, ESAEPS bIThi WAN'Xs, A) I&IPLICIT RHEAL*d (A- t,O-Z) DIMENSIN I1 A(IDN). LOGICAL WANIX,tID C IN1TIALIZE IZIhR, CClh13UT'E E'SA,pEP$S C ANoal 0. DO 185 1=1,N ITXR (I) - 0 ANi = 0. IF (I.NE. ) abi = BaidSitAi,1-1)) BNI = 0. DO 180 J=i,C ANI - aNI + DAUS (A (,J)) -bNI - Bi + BS (B (i, )) 10u CUNIliNU '- (AN '. AG.ANR) aiNOhd = AdL IF (iqI.Gti.dNu/LM) tibJLd - = J 1d5 CONTINUE IF (ALUOlid. EQ.O.) ACfhiM = EPS 1F? ( iLRM. LQ.O.) JNORd = EVS i&PSA = EPS*ANdilM1t iPS = EL)PS*biORti C dE&JUCIH A TO QUASI-TRI1ANGULAt:, KEEP 8 ThEii NGUL&R C OO F Li. L &. ) ~;'o () j30 C C CHECK fOd CGNVEHGICE OUR IEDUCIBILITY DU 220 LB=1,1 L = + 1-LE iF (L.EQ.1) GO TC 260 IF (ABS (A(L, L-1). LB. EPSA) GO I 2 J0 220 CONTINUE 230 A(L,L- 1) = 0. IF (L.LT.ii'-1);O 10 260 a = L-1 60 TO 200 C,;tr1CK 'OR SlA LL 1G C Gt d C 2b0 tI (tJA8S(d(L,L)).GT.ECS2B) ' TO 300 d (L,L) = 0. L = L+1 CALL ilSI12 (A (L,L),(L 1,L), J, U21, V,2) It' (F1l..NE.1.) GO TO 280 DO 270 J=L,N T = A(L,J) + U2*A(L1,J) A(L,J) = A(L,J) + T+V1 A (i1,J) = A(1,J) + T*V2 ',~ = t(L,J) + U2*13(L1,J) i tL,J) = B(L,J) + T*V1 d(LIJ) = tliL1,J) + T*V2 270 CCNTINU1E

230 L = L1 (;o rJ 2.30 C JE1GIN ONI O Z N 'rE, II ElATrioN Si1tATEGY 300 A1 = d - 1 L1 = L + 1 C ST'~ = 0.75 L~raR (M) = ITEH (d) + 1 iF (I:EaE (d).:2. 1) GO TO 335 iF (DA13StA(It,LA-1). L'L.CONJT*OLD1) GiO TC 305 IF (DAUi (A (d-1,t-2)).LT. CCL'ST*OLD') GC TO 305 iF (LTEL (Ml).i',.l3) GO T') 310 iF (LTEJd(l). U j30) o ro JdO C C, EkIfl COLULdl oU" A C 305 b I11 = (L, L) d22 = B (L1,L1) le. (BAdS (dJ2)..~.~~b) 1.12: EPSb F (IA5 (23).Lr.SLSb) B3J = EPSd dJJ4 = J(a,1i) 3 44 4 (Ml, i) ik (AS AS(d'4j. L.t.S J) L8 44 EPSb All = A (L,L)/B'i11 ai12 = A(L,I), l/iZ2,'L1 = A (L1,L)/dl I tk,2 = A (L 1 L )/i 2,t 33 = A(dl,t 1 ) /J3 A 34 -= (Li 1 l) /,)4 Li Ai 4 3 = A (, l) /dJJj A44 = A (t, M)/ b44 d I= 13 (L,LI /i Z t10 = a3J33-Il l) A4 4- A 1 1) - 34* A43 + 43*3d34*a11 )/Ail 1 + i'1l2 - All*b12 Ag0 = (A22-A11-Az1* 1i) - (A33-A11) - IA44-A11) + A43*B34 A30J = A (L+2,A1)/J22 GO TO 315 iUAD ilOC SlIFI C 310 AlO 0. AgO = 1. 1AJO = 1. l05 315 uLVI = DABS(Atdrf -11) CLI2 = DtaBSA(M-l,M-2) Li' (. NOT.WANTX) LL0 1 = L ii (WANTI) LOfll = 1 lk' (.NOT'~.iANT) MUOiRN = M ~LF (WANTX) 11tOa = N C C bEGIkN liAIN LO0P C DO 360 (=L,M1 MID = K.NE.al K1 =.K+1 K2 = K+2 K3 = K+3

if (K3.GT. 1) Ki = i KMll = K-1 IT tKl-1.LL.Lt),.11 = L C c ZSZROU A(K+I,K-1) AND A(t(+2,K-1) IF (K.EQ. L) CALL liSii3 (A 10,A 0A30,U1, tu2, 3,V1,V2, V3) IF (K.GT.L.AtiDL.LI[.i1l CALL 1S1i31A(iA,KAn1),A( K1t,Ki)A(K2,KLl1),A U1,U2,U3,V1,V2,V3) IF (K.EQ.M.1) CAiL hi12JA (KKM1),A(K1, Kl),UlU2, V1,V2) I ( (U1. E. I.) O I'0.325 [)O 320 J=K:1, 6LcdN ' = A(K,J) + U2*A (K1,J) i {MIID) I' = I + U3*A(K2,J) A (K,J) - A (,J) + T*V1 A (k,J) = A(K1,J) + f*;V2 1k (Alv) al t2,J) = A(K2,J) + T*V3 l' = i3u(tJ) + U'2*3U(n1,J) IE (Mli) L = 2 + UJ3*j(i(2,J) di(Ks,J) -= ( (,J) + TV1 B(Ki,J) = Ui(K11) + I*VJ2 I1F (MlI) L (K2,J) - d(K2,J) + T*VJ:3_20 CON rINUE IF (K.' G.oL),O L'u 325 A (KI,K-1) 0. IF (MiD) A (i42,K-1) = 0. C C:Zid{O AB(K+,K+I) Ai)D b(A 2,(K) 325 If (K.1U.M 1) jO 0 340 CaLL dSH3 (I(K, IJ),8(K2,K1, (,, U,) 1,UZ, 3,V I,VV3) If (U1. N:. 1.; Gu 10 34q DO 330 1.LOit'1 K3 T = A(i,r 2) + U4A (,Ki1) + U3*A{i,K) A (I, K2) = A t1,n2) + I*V1 A (1, 1=) A (,1) + 'L'+ vV 2 A(Z, K) = a (1,K) + Tl*3 r = Ou,'K) + d.U*B(L,K1) + GJ3 (iK) 3(i,1K2) - li,K2) + '*V 1 k (I, K1): o i,K 1) + T'*V2 8 (I,K) = ( (I,K) + T*V3 330 CONTINUE B (K2, K) 0. {(K2,K1) = 0. IF (.NUT. A d i) GO TO 340 DO 335 1=1,SJ ' = -X(.L,K2) + * + U3(X(i,K) (I, K2) = (LL,;2z) + ~*V1 X (I,K1) = X (i,K1) + t*V2 K(I,K) = lt(I,K) + T*V3 335 CONTItNU i C hAltO i3(K+l,t) C 34i CALL tlJli2t(3(iK1,t1),3(K1,ki,U1,[(1,1,V2) IF (U 1.. 1.) 3 'iti 3o9 DO 350 I-.LtOd1,K3 X = A(I,K1) + U2*A(I,K) A(.LK1) = A(l,K1) + 't*V1l A (I,K) = A(I,K) + T*~z

T = B(l,K1) + +U*i (,,K) B (I,K1) = 3J(~,K 1) + ~LV1 U (I, K) + T'IV2 350 tONTINUE 3(K1,K) = O. if(. N.T.IAdTIX) GU LO juO Do 355 I=1,,j I = X (~,:1 + U2+X(1,K) X(I, K1) = X(i,K1; + L*V1 X (Ik) =A (sAy) + TWi2 355 CONTINUE 360 CONTINUE C c END iA dAIN LUOt' C C Es.d ONE i 2 bri c 380 JO 385 1=l,f iTLiE (1) = -1 385 CCNTINUE 390 LATUtdil ciD)

JjLMPIC1T ']AL*0 (A-l,-gZ) i)LzIiLSUON A (LJ, Ni),0 (, ) ~i)), ALL"i tJ), ikLrl {i),L~kTA (;4),& (NV,NL)) LtGICAL wAiNTX, Fa. 1 C C kfiN EIGENifALUEi iJr UAJl-T1IAANGULAd LCIAisICE.j C i iO 400 THiRU 49SJ l:OUt, = ij S''L (-1 Gi — ) UN'I'L 1 C 400 COINTIhIUOi IE ( M.. U1) J; TO q 10 Iie (A(d,i Il- 1). Nt.) GO '! TO 420O C OL-BY-CNE SUbLihATILs CA IHEAL HC1T 410 ALFd ('.) = A Ld,i) UETA(i) = bili,M) A Ltl ( d) = O. A = =-1 UO:0 490 C TTWO-BY-T 3G SjUdibA lX C 420 L = 8I-1 I- (LAS(.8(L(t,tL)). T.k.ESi) Go IoU 425 3 (L,L) = 0 CALL tlS1tI2A(tL,L),A i1,L),U1,12,VIV2) 0O TO 460 425 IF (LABS(iJ(Lil,)) UTi'.' ) GO T.i 430 H (MP1, lS) J. CALL ti, i{tAi aI),A (ta, L),U1I,2,VV'I p Jj) UN = O. 430 AN D= uAbS (A (LL, ) ) +LAs (L,) ) )+ 3 (A (L,) j +A (A (, ) All A(L,L)/AN A12 = A (L, M) /i A21 = A (iH,L)/AN A2i = A( i,il)/A 1 811 B3(L,L)/blN b12 = B (L,e)/b N tt22 = i3(d,i1)/B1N E - Al l/B11 C = ((A2;-*d22) /U22 - tA2 1* 12) / (ll *i) ) /2. 1 = C*( + (A21*iAC1 '-t* L) )/(t11*E52) Ii (. L''.O.) GO TO 4 0O C C TrO RHEAL ROUlS 2: LHERO i)OT1i A (,L) AbLD {it,L) iF (c.GE.0.) h = E + (C + DS'nT(D)) If (C.LT.O.) ' = E + (c- - s t;) A11 = All - E*B11 A12 = A12 - ~*t 12 A22 = A22 - F*lw2 ffLIP = (LAd (A1 1) +uDAbS (A12) ).r. (;AUS (A 1) + AZ SiAA22) ) IF (IfLIP) CALL S 1 A 1 1 dSL 12(A12,A,I, U 2, V 1,V2) L' (.NOL. i'LIP) CALL iS t 2 (Ad222, A2 1,U1,U2, V 1,V2)

435 IE' (dIl.NL. 1.) GO i') 145 Do 440 1-l,rI T = A(l,dl) + *i*t, l,, A (I, L) - alil.e) + V1#'iT A(I,i.) h i1,L) + V2*1' I.= B( 1 i *4 U Z *t Li i i,) 8 (I,) o (i, + V *T B (i,L) = d (1.,-) + V4*T 440 CUs'~ItUr.f (.NOT. w ~AaT1) oTO '4 50 DO i44 iL=I,L T = A(l,) + U X IL) X (I, t) A (1, i) + V 1 'I X LI, L). (L) t 'V"' 4uS ' ida.~,J.J.) (uG '( 417. FLIP -= AN. lJ. blS I-) { b l) lF (FLIkP) Ca.L itiz L dJ (i t,L),{ (; L,), i1, L i 1, V2) iL- (.NT. I LiL') L ALL drii d (A (i., L) A (M,, L, 1 U Z, V 1,; 460 iF (U1. NE. 1.) D O I~ 4 7T DO 470 J=L, N T = A L,J) + L,*1I (M1,J) A (L, J) = i (L,J) + V1*'i A (nJ) = A (h,J) + V2*'.' T' = i(L,J) + u'Z*b(ii,J) 13(.1,,J) = JL,J) + V1+1' E3(id,J) = l (l,dJ) + V4'r 470 CN~'INULJ 4175 A (i,L) = 0. ti,(M,.) O. ALFR (L) = A L,L) ALFia(i1) = L(,, h) iBE'TA tL) = BtLL) Bi'LA (;:1 = (' 0, M) A tF1 (,i) = 0. A 11"F 1 ( L): O. M = b-z, Go 1,i 493 C iWO COitLLLX fliCtiLS 4bO J.. + C; = i),S,2, i (- D) A111 = I* '11 Alg2 = Alz - - i*B12 A 1l = EI*BlZ A211( = A21 A21. = O. A22R = A22 - zii*fJ2 A22.1 = ii*B22 FLIP = (iAB (A 1 1 R) + ti (A11I) +DABS (A 1 ) +DAGBS {(l/) ). k. (DMA JS (i 1 t) +J)A3as ([i2z2L) +3ABS (A 221) iF (FLI ) CA LL C1 d2(A1, AA 1 I,-A11,-i 'I I,CS,C iS L, S 1i) IF (.NOI'.fLI2) CAI.L C;iSHL2(A22R, A22I,2-Ai,2 -A21i,C4iL,,a(,SZ1i) FLIP = AN *.dz. (D (iBS (I) +DAJS (EI) ) *BN IF (FLIPk) CALL CdStitZ! (CZ4lJ11+SZi*1B11z, SZi12,312 ~1 +SM[*]322, SZI*2S2, CIz, SQi, S.I) If (.NOT.EFLLP) CALL C tiSdL2(C*A1l1+SMa*a i1, SZi*A1I, ~1 C3$~CA21+S+iSt*AZ, SZ;I*A22 C, 3!i;, SQl)

S::;d = Syaojfi + SiS$ShGi Si. = -iltabs - S j~iSote f,~ = tQ*Cg 11 1+ 4S$IC*CA1L + SR(*C(*A1 + SzEH*Ah2Z i'. = C S *i,'Ji*, -- Swii*CZ*A21 + SS1*AZ2 ria = c C, 4 u, 1 + C,:~i*algi S + SH+i3Z2 dAJl C= e*i 12 * lb + SSi 22 d = 1) S tii' t jkW t * + i) i4%JL) ltiA (L) = Li N4iLt kai ( ) -!\A A 4 ('T;a LJ k + I1*bjT)/ ALFL L) = AN* ('1i.*bbi - " iL*bU')/h '~R -.j*Tal l -.iC;*A12 1- C * $Zti*Az + C *CLZ*A2 TI = - iS'.*A I I - $ij*C *Al + C.*Szi*zl R 5,4%II _ $.S'11-:;&aC,,.:1 + *C; * 00.1 = - 'i.itA ol1 - ij.,':.'Ci;*oLz 1 A.L'tI li) (it + J A LCF (PI) bl L 430 1F (A.GiT.)o) GiO u 400 ( Lf O N _N iD

.L,1eLL.` dREALr\d (hi-i,6 DiMEii43IOU A (Ai, D b),tj (NI),),Ai) ALeh (N), ALiI (N),ft'rA (N),X (Mh, iN&) LUG iAL AJ 'A4TX,r'L. i C C I1NL EIGENVEi'CT'dLS 01 UhAU i-'i'lANLGU,.Alt iAhlktC-S tY 1itkCSUbS'i1Idi'Lj C US B3 D ' u IV N fEh,)llA1ti S.eORAt(;, 8-TI VkCiU)a Nz u (. L) C DOU 500 ~ItRU b)9 iUn,i = i STEP (-1 b -) JiNTiL 1 d = N 500 CONTiNUt if (ALF' (~i).Li.(0.) GO TO 55U C L;:AL VIC' Iii ALJ'd = ALYR(wl) B L(, L) = 1. C U1 51l0 T'l'U 54C O, i iL: l- ''~i:' (-1 i (O -) iNT1L 1 L = L-1 iF (L.,u.O)' O iC' 4)9 5 10 CONTINUE =l1 - L+1 SL = 0. DO:15 J=.1,. SL - 3L + (olMI'[t.(L,J) -a.i.. '*t (L,tJ *btlJ,) 515 CUNTINUi I' iL...1) GO 10 520 L1" (Bi.*ai`LL-1; ~ N.O. ) L'O 30. C dAdL 1-if-1 bLULK 520 0 = L4:lI.{. (i,L) -ALFN*b (L,L) j' (o. L. ) J = (thSA+E~,d) 12. d (L,) -O/l L = L-1 G$O TU 5 4J C; d;AIL 2-i-z3 BJUCK 530 K = L-1 S = 0. Jo 535 J=L1,i 535 C(SCON~it.lo TKK = dEPM*A(K,K) A.k'M*d(k(,K) TKL = Bh r2I*A(K,L) - A.Lit'b(K,L) T~LK = bET1*A (L,,K) TLL = b '~*i*A L,L) - ALFM*r (L, L) 0 = IKK*i'LL - 'KL*TI 1I? (C.tu.0.) i = (EL'S +;SB) /2. o (L, i1) L'J-;- KK*SL) /D El.IP = DADS (TKK).G. oAB~ {'TLK) Ii" (FL1P) d(K(,L1) = -(3Ki + 'T'L*B (L,ti) /Tai iF (.NCI '.Fl ) B(Kit) = -(S + LL (L,)/L L =L-.

54U Ie (L.G.1. U),JU I') 51 C C " CON-.PLEA. V ECu~E C 550 AiM.R: ALt(11-1) BETM = BE ETA(d - ) dl = M fi LF[ d-T1 CUildUi ENNi' = (j.1.** $0 "~dil B iS 'LLtAN&JUiJAi C i1 - 1) 5/ = -ibk rn (';'i C'l — AL n*'di i fl,) ) * (i- tH) / ( lDTe~Ll*H (ia, a1' B(a-l ai) = il,~ rl i*A ([, a) - A LMR*i (a, t') ) / (bE ] a'*A (a, rl)- 1) ) L: H- i" (L..OI) GO '.I.'C 505 560 CONTINJ4 U L1 = L+i kJt1,JTi) = -1. 1DO 565 J.=L1,,C TIi = BL',1*A (1, ) - aLl'l*O (Lo) rI - -Aadli$fl(L~,) SLii =.iLa + rhl*b(id#1) -: L1'~ (Ot.i) SLI = SLi ~ l'B>+ fIJ ) + T'i*u (Je(,f) 56u CONTINdJE IOI' (5L6 5 1 GG J 5170 I; (JiTl*' A IL, L- 1). iE'. i.) tcc o1 b75 C c CLUeLkX 1-Li'- I 0i, O. CALL C.; IVi-SLIi, - LL i)t Ui1 {(i lair), b(L~L'I) I L = L-1 Go '~ I5f5 C COdtPLiA 2-}Y-Zi ULOCK C 575 K = L-1 SKit 0. SNI = 0. DO 580 J=1I,m Te = -i:;i~- (K, J) ALMIi*I(K,J) TI = -ALl*Mib ]K, J) SKR = bB + 'IAJ*Jid aR) - f1*bj(J,!'-i) SKI = SKi + ~d* J, 1i) + Ti*E (, tiidj 580 CONTINU E TKKl = -iALiJ-* (&,K) TKLR = Brt;Id*h(B, L) - ALftII*EB (K,L) TIKLL - -Alr lII* K,L)

Cig,,*r (F/ (0 'r 4 1 +Z* (T ) -V' T) X) ) IFP * -" EG ~~9 01 09 (0'0TI.']") AT S?, OT OD ('O':'~) AT ( ( Y IT)) SI V!I + ((-,' T) Y r SI = = U c/~vT~r) v0: (W 'T) Y N'L =T 0t9 0C 0 -, T I n V.e ' 9 (U'T)X:.qf Of 09. (F'T, I'P,JT lW = C,'I=i 0 P9 0nO Q r)I.(t,~ J O Ft3TT) ~~~\00T.: nsIn h c;7 %A 01 J 0 9 )r ) T = 5f 9,J ol9f (0 ',9 r:',,) AT t-~?: *,, _ T;.]'t -rr oL O T N 0,' = rT 61 L~:~C~fl1 J'1..?7T 07V9 07 P I- N = W 00 ' 9(C',T, 9'Pt) AT C0(; o'' %: (', 0 (C *.T,9) P1. ci r~- =: (fTW' )? ' ('(T,,?'),~ '! TW ' '.T 7 () 7)q *T1.' T'TT,- (TWT)R.r. — T ' - T '(TWPr)P*T rT'TJ (Ts p. r) q*7T"T TIT,-..-)A.,-) AI 11 1r7 ('T'.; 9OTO') JT ' (p W ' I q T ) tT'; - ( T W T 7 I J T* )N J -T - T(T'? * T) rt t71 l + (P J *T) U H T T, JT-'P - ) ll T.T) I 73 I rGIVA -.I T F T I) sgPr gO a ( frTYp) vSnVrr V I T) sn r = d ra T((V wr IT+1''T) 'T 'T T I '"S T 7 'T 7S * I'M V I - T rT S;Jrl)-I Y S taNt T 5, ~7 (F.a,3 + vr F = 1, (T o''' M v 'C, ' (I,I q -.T T. 's 51'-T*T"[r TJ - )11"]J*rtYYJ + Tr 17l*FTYYL = T-.T":Tt* 'T' T, - /T7[ TJTV, - PT T*t, ' I Yv' = F(T (?" 1) n, t IV - (t'") 1 '. W,;;q = it VrX - U '= J -? P'I

Li X A( 1, - 1) DI - X(L,M) 655 CiN/ LNUI/ CALL CIV(A (1M-'I) # (M),L~1 bbu C. Li;J T 1' I N % b90 It' (I.Gi. 0) GO TO 0o30 C 70U iAETURN iEN)

S L U3 C 1A'l K Li A Li * ( ii tI 1AJ U I j 2, 3 Vl V v ) ) ia.oPLICiI' LtA-L*6 (i-,JC C t'iNSz it hOdiiitLOEi) iiANSI'JR&"'AT.i4N ~'AT llL /d1, i AZ AN L A3 C P = ~ + (V1, V2,fJJ) *(II1,Uz, UJ) *+4' C IF (A2.&.j..O. f A3 9.L3. ) GO '~) 10 s = i)AbS(A1) + AibSiA2) + ' A3dSi(AJJ d1i = A 1/S U2 = A2/S UJ = Aj/.5 V.3z - 3/ i{ = DSl' ({ll1*Ul I+ 0Z*i d+dJ*UJ) Vl = -(U1 + 1.)/a V2 = -U /a V3 = -O.3/.~ UI= 1, U3 = VJ/ V1 Li tT ild dd 1i ul = 0. C iit'T UI, 'I~LLC~i: Bk A.*d ( - h,- Z) i:C, P I + (V1,V 2) * J1,U2)*T -i;' i(A2.dO.) GO '1U 1T s DA S (A1) + DAbS 6('2) ul = A1/S U = A2/S it = SIafR (U 1 8 1.. d *J 2) lE (1i.iJ..)J d - -i Vi = -UZ/i U1 = 1. UZ = -2/V1.ia ii U NN 1 ~ l = OL [? I. N l kESV

'JSUBC1tUTiN C$i 1Zi (I i,fAIIvA2s A 2 '- I,,sS s Si) "WL j Aj;jAL4d (i Lfl-Oi,0-Z) C C COaiM1LEX tIOUSEtiOLibd Tiail WILL ZERU AZ (C s*) C, = kS -C ), C AiLAL. C0MIPI.EA C if (A2d*SQ..,D. A 2i. '.0.) uG ~0 1(0 if (A1R. zR 0.. A L41. ti 1 -0 ) ~ J0 zO e = oi)Str (A lAl I. tA 1 1i A 1 ~ ) Sd - (A 1dt*A L(+AI A Zi) /d Lt = DSL1T (C'T *C+Z*S*ul i *I L) C = C/8 10 C = 1. 6X = 0. Si = 0. Et rT r1 td N 20 c = 0.;Lt = 1. Si = 0. d TU LdN I 'lPI.L Ci d.iA.dJ (A-~,t 0-.) C Cu 'LA DJX TVl.E.. - / Ir (DA1$ (YL).L'. bAbS ({L)) ){o TL 10 ll -= YI/tk ) _ id + tid*lj, = (A4R + i*Xl)/L (-. - it*XR) /D L iET U I N L - + i Y ZR = {hs&ALt + XiL) /I) ZL (-L*XL - XA)I/L IND

Unclass ifiedc A _Y _,_S.,I f. r1i I.-. D OCCU i, NT CO.TR:OL DATh - R & D c, tjrity c.si.i l icati on ol title, h Li, e l' af cht irai 't t dad ind:xintl an Inotltficon ntj:;.t be c tct rtrd when tIhe overall rin, r? i:r I cl., 'sifiu, i G i GINA TI, N; ACTIVITY (Corporate aut),hor) '2-a. F:PORT SECURITY C._Ab.SIFICAl11,O Department of Mathematics Unclasssi fied University of Michigan h. GROUP i'.Ann Arbor. Michigan 48103 [,3. REPORT TITLE An Algorithm for Generalized Matrix Eigenvalue Problems 4 DfEscIriPTIVE NOTES (Type of report atld itlclisiveC a tes) Technical Report, January 1972 ii 5. AUTHOR(S) (F'irst name, middle initial, last name) C. B. Moler G. W. Stewart 6. REPORT DATE 7a. TOTAL NO. Of PAGES t7b. NO. o REFS January 1972 51 8a. CON TRAC T OR GRANT NO. Sa. ORIGINATOR'S REPORT NUME.R(S) N00014-67-A- 0181-0023 b. P.ROJECT NO. Technical Report No. 3 NR 044-377 gi C. | 9b. OTHER REPORT NO(S) (Any other numbers that may be assigned C.|~~~~~~~~ [ ~~~~~~~~~this report) d. 10. DISTRIBUTION STATEMENT b Approved for public release; distribution unlimited. II' 1. SUPPLEMENTARY NOTES 12. SPONSORING MILITARY ACTIVITY Mathematics Program Office of Naval Research ______________________ Arlington, Virginia 22217 tf 13. ABSTRACT A new method, called the QZ algorithm, is presented for the solution of the matrix eigenvalue problem Ax = XBx with the general square matrices A and B. Particular attention is paid to the degeneracies which result when B is singular. No inversions of B or its submatrices are used. The algorithm is a generalization of the QR algorithm, and reduces to it when B = I. Problems involving powers of X are also mentioned. I "'- _* FOR.. ' L' I (PAGE 1) sified S/,N 1.NOvr., 0 Sec Unclas sified S/N 010 10 807- 680 1 Security Classification

Unclas sified 14~~~~~~~~~~~~~~ 1 1 (1i\l'tV '1t NIf~>iCtIC~triu. ~4. A.... K E Y W O R D S ~~~~~~~~~L I N K A L. I ~. K ~5 ~,, ',- CKEY WORDS ROLE WT RO L WT F:O L -: V"_ S.'~~~ I s., matrix eigenvalues QR algorithm i, I! i! rl!i~ ~ ~ ~ ~ ~ ~ ~ fi I I r,. 'i:' I I ii~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ r i.. ~ ii i~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ I i~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -.O.-', _, o /,~ (BACK) Unclas sified (PAGE' 2) Security Classification