a I E> hahy 
 
 OF THE 
 U N IVLRSITY 
 Of ILLINOIS 
 
 510.84 
 IZQv 
 
 no. 257-264 
 cop. 2 
 
CENTRAL CIRCULATION AND BOOKSTACKS 
 
 The person borrowing this material is re- 
 sponsible for its renewal or return before 
 the Latest Date stamped below. You may 
 be charged a minimum fee of $75.00 for 
 each non-returned or lost item. 
 
 Theft, mutilation, or defacement of library materials can be 
 causes for student disciplinary action. All materials owned by 
 the University of Illinois Library are the property of the State 
 of Illinois and are protected by Article 16B of Illinois Criminal 
 Law and Procedure. 
 
 TO RENEW, CALL (217) 333-8400. 
 University of Illinois Library at Urbana-Champaign 
 
 JUN 2 8 1999 
 
 When renewing by phone, write new due date 
 
 Kal/Mit e\*>air1/Mlr> Hka ylnlA 11 f\/ 
 
 below previous due date 
 
 L162 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/eigenvalueproble258same 
 
6/J 
 
 Report No. 258 
 
 April k, 1968 
 
 ILLIAC IV Doc. No, 182 
 EIGEN -VALUE PROBLEMS 
 
 by we mm 0F THE 
 
 Ahmed Sameh AUG 15 !°C I 
 
 6NIVEISHY Of ILUNOIS 
 
 Luke Han 
 
Report No. 258 
 ILLIAC IV Doc. No. 182 
 
 EIGEN- VALUE PROBLEMS 
 
 by 
 
 Ahmed Sameh 
 Luke Han 
 
 April h, 1968 
 
 Department of Computer Science 
 University of Illinois 
 Urbana, Illinois 618OI 
 
 (This work was supported in part by the Department of Computer Science, 
 University of Illinois, Urbana, Illinois, and in part by the Advanced 
 Research Projects Agency as administered by the Rome Air Development 
 Center under Contract No. US AF 30(602)^UM. ) 
 
ABSTRACT 
 
 The most used methods for solving algebraic eigen-value problems 
 are discussed in detail in this report. These methods are Jacobi's, 
 Householder's and QR-algorithms. Attempts are made to adapt or even 
 modify the algorithm to take advantage of parallel computations. 
 
 A. Jacobi's Method 
 
 It proved to be one of the most effective methods on a parallel 
 computer. However, it is limited to symmetric matrices with dominant 
 principal diagonals. Evaluation of the eigen -values and eigen-vector 
 is rather straightforward once the method converges. 
 
 B. Householder's Method 
 
 This method reduces the symmetric matrix to the tridiagonal 
 form by means of elementary Hermitian orthogonal matrices. Once the 
 matrix is in the tridiagonal form, evaluation of the eigen-values and 
 eigen-vectors is performed by computational techniques that make almost 
 full use of parallelism. 
 
 C. QR -Algorithm 
 
 The QR-algorithm is used for finding the eigen-values of 
 unsymmetric matrices. For best results, the matrix is first reduced 
 to an upper Hessenberg form, using Householder's method. Convergence 
 is accelerated by performing a single origin shift; however, when the 
 matrix has some complex conjugate eigen-values, a double origin shift 
 is used in order to avoid complex conjugate shifts. 
 
 Flow charts and Tranquil language programs for all algorithms 
 are included in the report. 
 
TABLE OF CONTENTS 
 
 Page 
 
 1. Jacobi's Method 1 
 
 1.1 Introduction 1 
 
 1.2 Modification to the Classical Jacobi's Method 1 
 
 1.3 High-Level Language Code and Flow Chart 6 
 
 2. Householder's Method 6 
 
 2.1 Introduction 6 
 
 2.2 Theoretical Background 7 
 
 2.3 Computing Eigen-values and Eigen-vectors 10 
 
 2.3.1 Eigen-values 10 
 
 2.3.2 Eigen-vectors 12 
 
 2.k High-Level Language Program li+ 
 
 2.5 Flow Chart. . . 19 
 
 3. An Algorithm for Finding Eigen-values and Eigen-vectors 
 
 for Non-Symmetric Matrices 23 
 
 3.1 Introduction 23 
 
 3.2 Theoretical Background (QP -Algorithm) 26 
 
 3.3 Numerical Solution for QR Algorithm 30 
 
 3.^ High-Level Language Program 3^- 
 
 3.5 Flow Chart 37 
 
 k. References Ul 
 
 APPENDIX A (High-Level Language Program for Finding EMAX). . k2 
 
 APPENDIX B (QP -Algorithm With Double Origin Shift) kk 
 
1. Jacobi's Method 
 
 1.1 Introduction 
 
 There are many methods for finding eigenvalues and eigenvectors 
 of a symmetric matrix. Two of the most popular methods are Jacobi's 
 method and Householder -Givens method. In order to apply them on a 
 parallel computer effectively, we have made several modifications to 
 both of the methods. In this section, we will discuss Jacobi's method. 
 
 1.2 Modification to the Classical Jacobi's Method 
 
 Instead of using the Classical Jacobi's Method, i.e., going 
 through all the trouble to search for a pair (since it is symmetric) 
 of largest off diagonal element and to eliminate them, we modified 
 the method by eliminating all the off-diagonal elements of each 2x2 
 submatrices along the diagonal through orthogonal transformation. 
 At each iteration, this modification will remove n elements of a 
 n x n matrix„ We eliminate the elements by determining an orthogonal 
 matrix cp such that cpA 1 cp = A where A ~ having zeros in the approp- 
 riate locations. To eliminate those elements, we choose 
 
 cp = Diag. (T x , T 2 , ..., T^ 2 ) 
 
 [cos a sin a ~| 
 K = 1, 2, , n/2 
 sin a, -cos a ' ' ' ' 
 
 t t t 
 
 So T = T , Cp = cp and TT = I. 
 
 If matrix A is partitioned into 2x2 blocks as follows: 
 
 - 1 - 
 
A 1 = 
 
 A ll A 12 A 13 * ' ' A l n/2 
 , j , r _ 
 
 A 12 A 22 A 23 ' * * ! A 2 n /2 
 r- r p ----- r _ - 
 
 A t ! A* ! A A 
 
 H 13 | 23 ; 33 j • * * i 3 n/2 
 
 1- 1 1 1 
 
 • I • I • I la 
 
 • I • I • I I • 
 
 • • I •!•••! • 
 
 • I • I • I I • 
 
 f. , , , __ 
 
 A l n/2 ! A 2 n/2 j A 3 n/2 j ' ' ■ j A n /2 n/2 
 
 Then A 1+ = CpA 1 cp will be 
 
 i+1 
 
 T 1 A 11 T 1 
 
 TAT 
 1 12 2 
 
 
 T l A ln/2 T n/2 
 
 T A t T* 
 2 12 1 
 
 T A T 
 2 22 2 
 
 
 T 2 A 2n/2 T n/2 
 
 T n/2 A ln/2 T l 
 
 T A* T* 
 n/2 A 2n/2 1 2 
 
 • 
 • 
 
 T n/2 A n/2n/2 T n/2 
 
 Considering the diagonal siibmat rices A, , = T. A, . T, 
 to to kk k kk k 
 
 cos OL sin a, 
 sin a, -cos OL 
 
 a 2k-l,2k-l a 2k-l,2k 
 J*2k-l,2k a 2k,2k J 
 
 cos Oi sin cc 
 sin <X -cos OL 
 
 b 2k-l,2k-l b 2k-l,2k 
 _ b 2k-l,2k b 2k,2k J 
 
 where ^_ 1>2 ^_ ± u *2k-l,2k-l cos ~< 
 
 + a 2k,2k Sin2 \> 
 
 k "2k-l,2k Sin °k 
 
 - 2 - 
 
b 2k,2k = a 2k-l,2k-l Sin2 a k " a 2k-l,2k, ^ \ 
 + a 2k,2k C ° s2 \ 
 
 and b 2k-l,2k " 2 ( a 2k-l,2k-l ' a 2k,2k] 
 
 sin 2 a k - a^ cos 2 ^ 
 
 Now to eliminate the off -diagonal terms, i.e., b = 
 tan a k = ( a 2k _ 1;2 k ) /2 (a 2k-l,2k-l " a 2k,2k ) 
 
 1 x "I 
 a - — tan 
 
 2 a 
 
 2k- l,2k 
 
 ^ a 2k-l,2k-l " a 2k,2k> 
 
 If we define that 
 
 z _ I 
 
 Z k ~ 2 \ 
 
 1 / ([Ua2 2k-l,2k/ (a 2k,2k " a 2k-l,2k-l )2] +1 ^ 
 
 Then the rotation angles will be 
 
 sin <X = 
 k 
 
 — + Z. ; and cos a, 
 
 2k 7 1 
 
 \2 
 
 i-z. 
 
 and 
 
 ^k-l^k-l " 2 ( a 2k-l,2k-l + a 2k,2k) - 
 2k ,2k 
 
 \ 
 
 2 / a 2k-k,2k-l a 2k,2k x2 
 
 a 2k-l,2k + [ 2 ; 
 
 we also notice that 
 
 i+l> 
 
 tr (A 1 ) = tr (A 1+± ) 
 
 and since the norm is invariant under an orthogonal transformation of 
 the matrix. Therefore, 
 
 N 2 (A 1+1 ) = H 2 (A 1 ) 
 
 - 3 - 
 
Since 
 
 N 2 (A 1 ) = tr (A 1 A 1 ) 
 
 where A is the absolute of A. Also 
 
 2 2 2 2 2 
 
 therefore, b 2k-l,2k-l + b 2k,2k = a 2k-l,2k-l + a 2k,2k + 2 a 2k-l,2k 
 
 Hence by one orthogonal transformation, the sum of squares of 
 the diagonal coefficients increases by the sum of squares of the 
 (n) vanished terms. 
 
 Now, after one orthogonal transformation, the matrix will 
 look as follows: 
 
 A 
 
 i+1 
 
 \i° ! b i 3 \h\- ■ — |vv 
 
 b 22 " b 2 3 b 2U -! b 2,n-l b 2,n 
 
 J- 1 1 
 
 h 13 h 23 j h 33 o j ; 
 
 \k h 2k 1 ° \k I" " " -' " 
 
 .L 1 i. 
 
 J_ 1 J_ 
 
 _l_ _ _ _l_ - - -lh 
 
 1 j n-l,n-l 
 
 -■- - - -i- - - -0 h 
 
 1 ' ' n,n 
 
 Therefore, "by shifting the second row to the "bottom and 
 the second column to the far right, we bring to the diagonal new 
 submatrices with non-zero off-diagonal elements, and hence the 
 matrix is ready for the second transformation . This shifting is 
 repeated for n - 1 transformations. Then, we will shift 
 
 - h - 
 
the first row to the bottom and first column to the far right 
 repeat the above mentioned shuffling process for m transformation 
 (the diagonal terms are sufficiently dominating). Therefore, the 
 diagonal elements will be equal to the eigen-values, since the eigen- 
 values are not changed by orthorgonal transformations. For comparing 
 number of transformations required by this revised method with the 
 classical Jacobi's method, see Appendix F, first Q.P.R. The eigen- 
 vectors will be the column vectors of B, where 
 
 B = <P l t *2 ••• C 
 
 because if we let the final matrix be F, then 
 
 t t t „ 
 
 cp m . . . cp 3 cp 2 cp 1 A cp x cp 2 . . . 9 m = F 
 
 but AX = XX 
 
 Since F is almost diagonal with X's being its diagonal elements, 
 
 . * . AX = FX where X is the matrix with eigen-vectors as 
 
 its columns 
 
 t t t _ t t t 
 
 Hence A cp x cp 2 • • • <f> m = F cp 1 cp 2 ... <P m 
 
 and B = cp., cp^ ... CD is the matrix with the corresponding 
 
 T l ^2 T m 
 
 eigen-vectors. 
 
 Thus, in practice one can follow the following four steps: 
 
 (i) Calculate orthogonal transformation matrix cp by 
 finding n rotation angles and construct 2X2 submatrices on the 
 diagonal and zeroes everywhere else. 
 
 (ii) Pre and post-multiply matrix A at i-th step by 
 . transformation n 
 ;et B by B cp where B , initial B, is identity. 
 
 orthogonal transformation matrix cp to get A . At the same time, 
 
 - 5 - 
 
(iii) Testing whether the ratio between the sum of the squares 
 of the diagonal elements and the sum of the squares of the off -diagonal 
 elements is greater than (l/£), where £ is a given small number (.001 - 
 .01). If greater, print out the eigen-values and the eigen-vectors from 
 the present matrices A and B respectively, otherwise go to next step. 
 
 (iv) Shuffle the second row and second column to start another 
 transformation. After (n-l) such shuffling, shuffle the first row and 
 first column and start all over again. If after (n-l) shuffling of the 
 first row and column, the diagonal elements are still not sufficiently 
 dominant, the method fails. 
 
 Jacobi ' s method has proved to be suitable for diagonally- 
 dominant matrices. Therefore, although Householder's method works for 
 any symmetric matrix, due to the simplicity of Jacobi' s method, it will 
 be used for finding eigen-values and eigen-vectors for diagonally-dominant 
 matrices. 
 
 1.3 High-Level Language Code and Flow Char t 
 
 For flow chart and high-level language code see ILLIAC IV 
 Document Number 165. 
 
 2. Householder's Method 
 2.1 Introduction 
 
 This method is one of the most effective methods that exist 
 for reducing a symmetric matrix to tri-diagonal form using elementary 
 Hermitian orthogonal matrices. For a matrix of order n, there are 
 
 (n-2) steps in this reduction, in the r of which [n-(r+l)] zeroes 
 
 th th 
 
 are introduced in the r row and r column without destroying the 
 
 zeroes introduced in the previous steps. 
 
 - 6 - 
 
2.2 Theoretical Background 
 
 The present problem is that, for a given vector x of order 
 n, it is required to determine an elementary Hermitian matrix P such 
 that pre-multiplication by P leaves the first component of x and 
 eliminates the rest, i.e., 
 
 Px - ke (2.1) 
 
 where 
 
 M 
 
 constant 
 
 e = ( 1, 0, 0, ...., 0} (column vector) 
 
 Let P = I - 2ww t (2.2) 
 
 where w is a column vector of order n and w w = 1.0. P is not 
 only elementary Hermitian but also orthogonal for 
 
 P* P = (I - 2wwV (I - 2ww t ) 
 = I - Ww J + Uw(w w)w 
 
 = I 
 
 from Equ. (2.1), it can be proved that the Euclidean norm of x is 
 
 invariant, i.e. 
 
 2 2 2 2 2 
 S = x 1 + x 2 + x n = k , 
 
 f t 
 
 2-Norm of Px = (Ex) (Ex) = (Ex) (Ex) 
 
 = x t P* F = xVpx 
 x 
 
 = x x which is the 2-norm of x, 
 
 Therefore, k = + S (2.3) 
 
 hence, equation (2.1) yields, 
 
 x - 2w K = + S 
 
 x. - 2w. K = (i = 2,3, ...,n) 
 
 (2.U) 
 
 - 7 - 
 
where K = w x which is a scalar. Squaring and adding equations (2.4), 
 the value of the scalar K can be obtained 
 
 K 
 
 \I77 
 
 x. S 
 
 (2.5) 
 
 If u = (x +- S, x , x„, ...., x ), equations (2.k) yield, 
 
 w = 
 
 2K 
 
 and 
 
 P = I 
 
 uu 
 
 2K 
 
 (2.6) 
 (2.7) 
 
 For a matrix A of order n, suppose that the configuration of the matrix 
 just before the rth step is as follows: 
 
 A 
 
 r-1 
 
 
 
 
 1 
 
 
 
 
 c 
 
 r-1 
 
 1 
 1 
 
 1 
 
 ^Vx 
 
 
 
 T 
 1 
 
 b 
 r- 
 
 _ T 
 
 ■l 1 
 
 r-1 
 
 1 — 
 
 1 
 
 
 1 
 
 
 (2.8) 
 
 where C n is a tri-diagonal matrix of size r, b , is a column vector 
 r-1 r-1 
 
 of order (n-r) and B n is a matrix of size (n-r). The matrix P of 
 the rth transformation may be expressed in the form, 
 
 " I » " 
 1 
 
 
 ' I ' "1 
 
 I 
 
 . o ' Q r _ 
 
 
 • I-2v v t 
 l. r r -* 
 
 (2.9) 
 
 where v is a unit vector having (n-r) components, w = (0 | v ). 
 Hence, we have 
 
 A =P A _ P 
 rs r r-1 r 
 
 'r-1 
 
 _ _ r _ 
 
 i C 
 
 r r-1 r 
 
 (2.10) 
 
 - 8 - 
 
where c = Q b , 
 r r r-1 
 
 If we chosse v such that all components of c vanish except for 11 
 first component, then A is a tri -diagonal matrix as far as its first 
 (r-t-l) rows and columns are concerned. 
 
 Making use of equations (2.2) through (2.7), the method 
 can be organized in a fast systematic manner, 
 
 P = I - (u u t /2K 2 ) 
 
 where 
 
 where 
 
 u. r = (i = 1, 2, ...., r) 
 
 u ^ = a J _ 1 ~ S (2.11) 
 
 r,r+l r,r+l r 
 
 u. = a . (i = r + 2, . . . . , n) 
 lr ri 
 
 n 1 
 
 S ( E a^ .)2 
 r . . ri 
 i=r+l 
 
 1 
 
 K = |"(S 2 - a j _. S )/2] 2 
 r r -t- r,r+l r ' 
 
 In equations (2.11) the sign associated with S is to be 
 
 taken as that of a , . 
 r,r+l 
 
 The method discussed above involves quite a large number of 
 multiplications which may be necessary in case of reducing a non- 
 symmetric matrix to an upper Hessenberg form. However, for symmetric 
 matrices, the above method can be modified slightly in order to take 
 advantage of the present symmetry. 
 
 A=PA,P (P t = Pj 
 
 r r r-1 r r r 
 
 using equation (2.7), and denoting that, 
 
 - 9 - 
 
(A . u )/2K 2 = p and (u t A J/2K 2 = p t (2.12) 
 
 r-1 r ' r r r r-1'' r r v ; 
 
 .'. A = A - - p u t -u pSu (u t p ) u t /2K 2 (2.13) 
 
 r r-1 r r r r r r r r ' r ' 
 
 again defining, 
 
 then, 
 
 4 . =P. "|u (u* p /2K 2 ) (2.14) 
 
 r ^r 2 r r r r 
 
 A r - \_i - % \ - (q r \^ (2-15) 
 
 The matrix of interest is (Q B Q ), Equ. (2.10), which 
 
 is of order n-r, since the elements in the first (r-l) rows and 
 
 columns are the same as the corresponding elements of A n . Therefore, 
 
 r-1 
 
 the transformation matrix P is designed as to eliminate a . (i = r-f-2, 
 
 . ...,n) and to convert a _ to + S , where S is defined before. 
 
 r,r+l — r' r 
 
 Hence the vector p has its first (r-l) elements equal to zero, Equ. 
 (2.12), and the rth element is not needed. Therefore there are 
 (n-r) multiplications required for the evaluation of p . Similarly, 
 it can be proved that there are (n-r) multiplications involved in 
 
 the computation of q . Finally the total number of multiplications 
 
 r 2 
 
 involved in the computations of the elements of A is [2(n-r) + 
 
 (n-r)], but since the matrix is symmetric, the number of multiplica- 
 
 tions concerning the upper triangle only is (n-r) + (n-r)/2, the 
 
 second term may be neglected compared to the first, .'. No. of 
 
 2 
 multiplication for one transformation ~ (n-r) . Hence, the total 
 
 number of multiplications to reduce the matrix to a tri-diagonal one 
 
 2 3 
 is essentially rr n which is half of that of Givens' method. 
 
 2.3 Computing Eigen-Values and Eigen-Vectors 
 
 2.3.1 Eigen-Values 
 
 Once the matrix has been reduced to the tri-diagonal form, 
 there are so many techniques for getting the eigen-values, for example 
 (Bisecting method, Sturm sequence or Gerschgorin discs). 
 
 10 - 
 
An effective method for computing eigen-values that are most 
 
 suitable for parallel computers can be summarized as follows. The 
 
 largest and smallest eigen-values can be obtained by the power method. 
 
 Dividing the interval between \ . = \. and \ = \ into (n-l) 
 
 nun 1 max n ' 
 
 divisions and get \. (j=l,2,....,n). For each \., the determinate 
 
 J J 
 
 of the tri -diagonal matrix can be obtained, using the method of prin- 
 cipal minors 
 
 f (\.) = 1 
 o y 
 
 f l ^j) - Ca n - V,) 
 
 (2.16) 
 
 f l < X J> = ^Sk- V f k-l ( ^ } " Vl,k WV 
 
 supposing that the value of each determinate corresponding to X . is 
 
 J 
 R. ( j=l,2, . . . . ,n) a plot can be made between \. and R., Figure a. 
 
 R 
 
 T 
 
 mm 
 
 Every R=0 indicates the 
 position of an Eigen-value 
 
 max 
 
 Figure (a) 
 
 11 - 
 
Each zero crossing indicates a correct eigen-value, supposing that 
 
 X. gave use to a positive residual R.,_ while X. gave a negative 
 J+l J+l J 
 
 R., then, again the interval between \. and \. can be divided to 
 
 J J J -^- 
 
 get the exact eigen- value with enough accuracy. 
 
 2.3-2 Eigen-Vectors 
 
 If x is the eigen-vector of the tridiagonal matrix (A ) 
 corresponding to an eigen-value \, then the corresponding eigen-vector, 
 V , of the original matrix A is given by, 
 
 V = P 1 P 2 ••■ P „-2 X 
 
 To avoid storing the matrices P., the vector v is calculated 
 from the relations, 
 
 P x = (x) 
 n-2 v 7 n-2 
 
 P n _3 (x) n-2 = W n -3 
 v = p 1 ( x ) 2 = ( x \ 
 
 (2.17) 
 
 So, it remains to find the eigen-vectors of the tri-diagonal 
 
 SS' 
 
 in the form, 
 
 matrix. Assume that the tri-diagonal matrix (A „) can be expressed 
 
 n-2 
 
 n-2 
 
 c l 1 
 
 b l i 
 
 
 
 
 
 
 
 
 
 
 
 
 
 L i 
 
 \ 
 
 
 
 
 
 
 
 
 
 
 
 c s 
 
 • 
 
 
 
 
 
 
 
 
 
 
 . 
 
 . 
 
 
 
 
 
 
 
 
 
 • 
 
 c. 
 
 1 
 
 b. 
 
 l 
 
 
 
 
 
 
 
 
 
 b. 
 
 1 
 
 c i + i 
 
 b. , 
 l+l 
 
 
 
 
 
 
 
 
 
 b i + i 
 
 • 
 
 • 
 
 
 
 
 
 
 
 
 
 • 
 
 • 
 
 
 
 
 
 
 
 
 
 
 
 # 
 
 - 12 - 
 
for an eigen-value X and an eigen-vector x, where x = (x , x , ..., 
 x ) , of (A ,,) the following relation can be obtained, 
 
 (C ± - X) x ± + b 1 x 2 = 
 
 Vl x i-l + (C i " X ) x i + b i x i + i = ° U=2,3,...,n-1) 
 Vl x n-l + K " M x n = (2.18) 
 
 If x = 0, the rest of the components of the vector x are zero, so 
 x will be assumed as 1. From the equations (2.18) and (2.l6) an 
 expression for x (r=2,...,n) may be obtained in form, 
 
 x r = (-l) r_1 f (X)/(b 1 b 2 ... b r ) (r=2,...,n) (2.19) 
 
 Similarly all the eigen-vectors of the original matrix A may be 
 
 o 
 
 obtained, each corresponding to an eigen-value. 
 
 A flow chart and a higher level language program for this 
 algorithm to find all eigen-values are in sections 2.k. 
 
 - 13 
 
.4 High Level Language Program 
 
 #BEGIN 
 
 #LABEL 
 
 #ARRAY 
 
 START 
 
 L00P1 : 
 
 #REAL 
 
 #REAL 
 
 #SET 
 #READ 
 
 START, L00P1, LOOPE, LP, LLP, LI, LL1, L2, SEUP, REPE, 
 
 FINAL, FN>L, EEND; 
 
 #REAL #FLOAT #SHORT #STRAIGHT 
 
 A[ 614,64], MQ[64,64], TMQ[ 64,64], FNT[64,64], FCN[64,64], 
 
 UT[64,l], U[l,64], P[l,64], Q[l,64], EVALUE[64], 
 
 FNTO[64], FNTl[64], FTO[64], FTl[64], NEVAL[64]; 
 
 #FLOAT #SHORT 
 
 S2, S, T, TEM, GRAG, INTV, STRT, FINI, RANG, STP, 
 
 SSTP, INIA, FINB; 
 
 #FIXED 
 
 J, CNT, NCNT, EV, TK, JK, NUM, NUM1, H, L, K, I, N, M; 
 
 HH, LL, KK, II, JJ, IJ, JI, 
 
 A; 
 
 I, NN; 
 
 J *- 2; 
 
 HH : : [J,J+1, , 64]; 
 
 LL :: [1, 2, ...., J-l] ; 
 
 KK : : [J+1, J+2, , 64]; 
 
 S2 <- #FOR (H) #SLM (HH)#DO #SUM(A[ J-l,H]t 2) ; 
 
 S <- SQRT (S2); 
 
 T ^ S2 - A[J-1,J] x S; 
 
 #FOR (L) #SIM (LL)#DO UT[L,l] - 0.; 
 
 UT[J,1] - A[J-1,J] - S; 
 
 #FOR (K)#SIM (KK)#DO UT[K,l] - A[j-1,K]; 
 
 U «- TRANSPOSE [UT]; 
 
 P ^ A x U; 
 
 P - P/T; 
 
 TEM <- UT x P; 
 
 TEM *- TEM/T; 
 
 Q *- .5 x TEM x U; 
 
 Q - P-Qi 
 
 - 14 - 
 
MQ «- Q x UT; 
 MQT «- TRANSPOSE[MQ]; 
 A - A - MQ - MQT; 
 J «- J+lj 
 #IF J < 63 #THEN #GOTO L00P1; 
 NCNT «- 0; 
 EV «- 1; 
 
 GPAG «- EMAX[A] - EMIN[A]; (SEE APPENDIX A) 
 INTV <- GRAG/8. ', 
 STRT <- EMAX[A] ; 
 II :: [1,2,... .,64]; 
 JJ :: [3,h,....,6k]; 
 LOOPE: FINI «- STRT - INTV; 
 
 RANG <- STRT - FINI; 
 STP <- RANG/ Gk.; 
 
 #FOR (I) #SIM (II) #D0 
 
 #BEGIN EVALUE[I] «- FUNI + STR x (l-l); 
 
 FNT0[I] <- 1.; 
 
 FNT1[I] <- A[l,l] - EVALUE[I]; 
 #END; 
 
 #FOR (K) #SEQ (JJ) #D0 
 
 #FOR (I) #SIM (II)#D0 
 
 #BEGIN 
 
 FNT[K,I] «- (A[K,K]-EVALUE[I]) x FNTl[l] 
 
 A[K-1,K]T2-FNT0[I]; 
 
 FN0[I] - FNT1[I]; 
 
 FNT1[I] *• FNT[K,I]; 
 #END 
 
 NCNT - NCNT + 1; 
 TK «- 1; 
 CNT «- 0; 
 
 15 - 
 
LP : #IF TK > Gk #THEN #GOTO FINAL #ELSE . 
 
 #IF FNT[6U,TK] < 0. #THEN #G0T0 LI #ELSE 
 
 #IF FNT[6U,TK] > 0. #THEN #G0T0 L2 #ELSE #D0 
 
 #BEGIN 
 
 #PRINT EVALUE [TK]; 
 EV «- EV + 1; 
 
 #IF EV > 63 #THEN #G0T0 EEND #ELSE #D0 
 #BEGIN 
 
 TK <- TK + 1; 
 #G0T0 LP; 
 #END; 
 #END; 
 
 LI : #IF TK > Gh #THEN #G0T0 FINAL #ELSE 
 
 #IF FNT[64,TK+l] < 0. #THEN #D0 
 
 #BEGIN 
 
 TK <- TK+1; #G0T0 LI; 
 #END #ELSE 
 #IF FNT[6U,TK+l] > 0. #THEN #D0 
 
 #BEGIN 
 
 TK - TK+1; #G0T0 SEUP; 
 #END #ELSE #D0 
 #BEGIN 
 
 #PRINT EVALUE [TK+1]; 
 EV - EV + 1; 
 
 #IF EV > Gh #THEN #G0T0 EEND #ELSE #D0 
 #BEGIN 
 
 TK <- TK+1; #G0T0 LI; 
 #END; 
 #END; 
 
 L2: #IF TK > 6^ #THEN #G0T0 FINAL #ELSE 
 
 #IF FNT[64,TK+l] < 0. #THEN #D0 
 
 #BEGIN 
 
 TK *- TK+1; #G0T0 SEUP; 
 #END #ELSE 
 
 - 16 - 
 
#IF FNT[6>+,TK+l] > 0. #THEN #D0 
 
 #BEGIN 
 
 TK ♦- TK+1; #G0T0 L2; 
 #END #ELSE #D0 
 #BEGIN 
 
 #PRINT E VALUE [TK+l]; 
 EV *- EV + 1; 
 
 #IF EV > 6k #THEN #G0T0 EEND #ELSE 
 #BEGIN 
 
 TK <- TK+1; #G0T0 L2; 
 #END; 
 #END; 
 
 SEUP MM «- CNT x 16 + 1; 
 
 NUM1 - HUM + 15; 
 
 NN : : [1,2, ,16]; 
 
 MM :: [NUM, WUM+1, ,NUM1] ; 
 
 INIA «- E VALUE [TK-l]; 
 FHB <- EVALUE [TK]; 
 SSTP - (FINE - INIA)/l6. ; 
 #F0R (N) #SIM (NN) #D0 
 #F0R (M) #SIM (MM) #D0 
 #BEGIN 
 
 NEVAL [M] - INIA + SSTP x (N-l); 
 FT0[M] <- 1. ; 
 
 FT1[M] *- A[l,l] - NEVAL[M]; 
 #END; 
 CNT «- CNT + 1; 
 
 #IF CNT < k #THEN #G0T0 LP #ELSE 
 
 REPE TEM <- CNT x 16 
 
 IJ : : [1,2, ,TEM]; 
 
 JI :: [3,k, ,TEM]; 
 
 #F0R (K) #SEQ (JI) #D0 
 #F0R (I) #SEQ (IJ) #D0 
 #BEGIN 
 
 FCN[K,I]«- (A[K,K]-NEVAL[I]) x FTl[ I] - 
 A[K-l,K]t2 x FT0[I]; 
 - 17 - 
 
LLP 
 
 LL1 
 
 FINAL 
 FNAL: 
 
 
 FT0[I] «- FTl[l]; 
 
 
 FT1[I] - FCN[K,I]; 
 
 
 #END; 
 
 
 JK *- 1 
 
 #IF 
 
 JK > Gh #THEN #GOTO FINAL #ELSE 
 
 #IF 
 
 FCN[6^,JK] = 0. #THEN #D0 
 
 
 #BEGIN 
 
 
 #PRINT NEVAL [JK]; 
 
 
 EV <- EV + 1; 
 
 
 #IF EV > 6k #THEN #G0T0 EENL #ELSE #D0 
 
 
 #BEGIN 
 
 
 JK *- JK.+ 1; #G0T0 LL1; 
 
 
 #END; 
 
 
 #ENL #ELSE #D0 
 
 
 #BEGIN 
 
 
 JK - JK + 1; #G0T0 LLP; 
 
 
 #END; 
 
 #IF 
 
 FCN[64,JK] = #THEN #G0T0 LL1 #ELSE #D0 
 
 
 #BEGIN 
 
 
 JK <- JK+1; #G0T0 LLP; 
 
 
 #END; 
 
 #IF 
 
 CNT = #THEN #G0T0 FNAL #ELSE #G0T0 REPE; 
 
 #IF 
 
 NCNT < 8 #THEN #D0 
 
 
 #BEGIN 
 
 
 STRT «- FUJI; #G0T0 L00PE; 
 
 
 #END #ELSE #D0 
 
 
 #BEGIN 
 
 
 #PRINT "SOME ROOTS ARE EITHER TOO CLOSE 
 
 
 TOGETHER OR THERE EXIST MULTIPLE ROOTS"; 
 
 
 #STOP; 
 
 
 #END; 
 
 EEND : #END 
 
 18 
 
2.5 Flow Chart 
 
 HOUSEHOLDER TRIDIAGONALIZATION 
 
 S^ = S2 
 
 2 
 2K = T 
 
 H «- 2 
 
 2 64 
 S^-Z [A((H-1),J)] 2 
 J=H 
 
 2K 2 - S 2 -A(H-l,H)x S 
 
 L«1,2,...,H-1 > 
 
 UT(L) «- 0. 
 
 UT(H) <- A(H-1,H)-S 
 
 K=H+1,H+2,...,N » 
 
 UT(K) - A(H-1,K) 
 
 U *- TRANSPOSE (UT) 
 
 [II 
 
 <D 
 
 - 19 
 
p *- a x u/ac' 
 
 n td ! tt / UTxP N 
 Q «- P-^- x U x ( — s-) 
 
 2K 
 
 A - A - Q x UT - TRANSPOSE^ x UT) 
 
 H «- H + 1 
 
 H < 63 
 
 yes 
 
 <D 
 
 FINDING EIGEN -VALUES 
 
 no 
 
 FIND EMAX AND EMIN 
 
 DIVIDE INTO 8 SEGMENTS 
 BETWEEN EMAX AND EMIN 
 
 NCNT «- 
 
 EV <- 1 
 
 20 
 
3=1,2,... ,6k 
 
 K=2,3,...,61+ 
 
 NCNT *- NCNT + 1 
 
 th 
 DIVIDE THE NCNT 
 
 SEGMENT INTO 6k INTERVALS 
 
 f (E.)=l 
 
 f (E )=(A(1,1)-E.) 
 
 f k (E.) = (A(K,K)-E.)xf k _ 1 (E.) 
 -A 2 (K-l,K)xf k _ 2 (E j .) 
 
 PRINT THE EIGENVALUES FOR WHICH 
 f 6[f (E)=0 AND INCREMENT EV 
 
 WHEN EACH EIGENVALUE IS BEING FOUND 
 
 yes 
 
 STOP 
 
 No 
 
 DIVIDE EACH PAIR OF E . and E ' s 
 
 FOR WHICH ^(E.) HAVING DIFFERENT 
 
 SIGH FROM fV^E.+l) INTO l6 INTERVALE 
 
 AND DO FOUR SUCH PAIR AT A TIME 
 
 © 
 
 21 - 
 
0* 
 
 No 
 
 
 <& 
 
 No 
 
 F (E.) = l 
 o v j y 
 
 F 1 (B j )=(A(l,l)-E j ) 
 
 F K (E J )=(A(K,K)-E J )xF K _ 1 (E J ) 
 -A 2 (K-l,K)xF K _ 2 (E J ) 
 
 PRINT THE EIGENVALUES FOR WHICH 
 f 6 ^(e)=o AND INCREMENT EV AT 
 
 EACH EIGENVALUES 
 
 SJ?ECIAL._CASE 
 yes 
 
 NCNT > 8 
 
 yes 
 
 PRINT 
 "SOME ROOTS ARE 
 MULTIPLE ROOTS" 
 
 AFTER 
 
 "YsTOPj 
 
 H 
 H 
 
 O 
 
 o 
 
 H 
 1-3 
 
 I 
 
 0^ 
 
 -p- 
 
 9 
 
 22 
 
3 • An A l gorithm For Finding Eigen-values and Eigen-vect or: i 
 For N on - Symmetric Matrices 
 
 3 • 1 Int roduction 
 
 There are two algorithms for dealing with unsymmetric 
 matrices. The first one is (LR) which was developed by Rutishauser 
 (1953). The LR-algorithm has proved to be a powerful method for 
 finding the eigen-values of symmetric band matrices. However, it 
 proved to be inferior for handling the ei gen- value problem of large 
 non- symmetric matrices, the most important difficulties are, 
 
 (i) Triangular decomposition may not always be 
 possible. 
 
 (ii) The amount of computation required by the method 
 is very great. The LR-algorithm is essentially 
 an iterative procedure, where the method consists 
 of a series of similarity transformations on the 
 original matrix, 
 
 A ± = L x R x (3-1) 
 
 where, 
 
 A is the matrix under consideration (size = n) . 
 
 L is a lower triangular matrix. 
 R is an upper triangular matrix 
 
 i.e. , 
 
 L = 
 
 1,1 
 
 'n,l 
 
 where I. . 
 11 
 
 1, i = 1, n 
 
 -i 
 
 n,n 
 
 r r 
 
 11 In 
 
 nn 
 
 - 23 - 
 
by similarity transformation, 
 
 hence 
 similarly 
 
 hence 
 
 let, 
 
 A 2 = Ll A x L x 
 
 A 2 = R l L l 
 
 Vl BS ^k\\ ••' L- 1 )(A 1 )(L 1 L 2 ... L k ) 
 
 (L X L 2 ...L k ) A k+1= A X ( Ll L 2 ... L k ) 
 
 \ = L l L 2 
 
 k 
 
 \ = W-i ••• R i 
 
 therefore, T fc U R = L Lg ••• 1^ (i^ \) \_ ± - • • \ ^ 
 
 = A l L l L 2 • ' * L k-1 \-l ' ' * R 2 R l 
 
 and so forth, 
 
 finally 
 
 T . U k = < A 1> 
 
 (3-2) 
 
 (3-3) 
 
 (3-4) 
 
 and it was shown that, if certain conditions are fulfilled then, as 
 k -* 00, then (A ) tends to an upper- triangular matrix, where the 
 diagonal elements are the eigenvalues in order of their modulus, 
 (the first being the largest), i.e., 
 
 (A x ) k = 
 
 X 
 
 o \ 
 
 ■\ 
 
 n 
 
 (3.5) 
 
 follows: 
 
 The restrictions on the LR-algorithm nay be summarized as 
 
 1. Eigenvalues must be distinct since convergence depends 
 upon the ratio (\ /\ ) and will be very slow if 
 
 - 2k - 
 
separation of the eigenvalues is poor. 
 
 2. Triangular decomposition of A must exist, (it fails 
 if any of the principal minors of A vanishes). 
 
 3« The leading principal minors of the modal matrix, 
 (matrix whose columns are the eigenvectors), and 
 its inverse should not vanish. 
 
 h. Even if the triangular decomposition exists, there is 
 a large class of matrices where triangular decomposi- 
 tion is numerically unstable which may lead to gross 
 errors in the values of the computed eigenvalues. 
 
 2 3 
 5. The volume of computation is very large, — n multi- 
 plications for each step of iteration, half of them 
 in triangulation and the other half in post-multipli- 
 cation. 
 
 Due to the above mentioned restriction on LR, QR proved to be 
 a much more efficient algorithm, and it will be discussed here in more 
 detail. 
 
 A more stable method of triangulation may be achieved by usin^ 
 elementary unitary transformations (QR-algorithm). The QR-algorithm is 
 defined by the relations, 
 
 A k = \ R k 
 
 , (3-6) 
 
 A k+1 = \ A k \ = R k \ 
 
 where R is an upper-triangular matrix, 
 Q is a unitary matrix, i.e., 
 Q = Q _1 
 
 The unitary-triangular decomposition of any square matrix exists. 
 
 - 25 - 
 
3.2 Theoretical Background (QR -Algorithm) 
 
 (l) For any matrix A (order n), there exists a unitary matrix 
 Q such that 
 
 \~ w 
 
 where, R, is an upper triangular matrix which has real positive diagonal 
 
 elements, and Q, = IL N n . . • W (3-7) 
 
 ' k 1 2 n w ' ' 
 
 where 
 
 N. = 
 
 l 
 
 J i-i i ° 
 
 • M. 
 
 (3-8) 
 
 M. is of order [n-(i-l)], the unitary matrix M can operate on any 
 vector b such that, 
 
 M b = | |b | | e 
 
 (3-9) 
 
 where 
 
 | To | | is the Eulerian norm of the vector b, e, = {1,0,0,...0} 
 
 M is an elementary unitary matrix, that is a matrix which 
 
 differs from an identity matrix at most in one principal 2x2 submatrix. 
 
 Therefore, if the vector b is of order m, M = T T -, T„ T n (3-10) 
 
 7 ' m m-1 d 1 
 
 where, 
 
 t- n t. 
 
 11 lr 
 
 T = 
 r 
 
 1 
 
 i N i 
 
 i \ I i 
 
 -1 
 
 t n t 
 
 rl rr 
 
 - 26 - 
 
The elements of T can be easily expressed in terms of the components 
 
 r 
 
 of vector b, 
 
 <o r -(t„ ) r = [i . JMi.2] I 
 
 ll 7 rr' L |b_ 
 
 p<r » 
 
 (3.11) 
 
 (t n r = -(t j 1 = b / zib 
 
 \ ]_ r / v rl' r' ' p 
 
 p<r 
 
 Hence, if we consider that the matrix A, consists of the 
 
 K. 
 
 column vectors C., j =1, • • • , n, therefore, 
 J 
 
 N 1 !^ N* N^ A. = R. (3-12) 
 
 n n-1 2 1 Tc Is. \~> / 
 
 where N. is defined by Eq. (3*8), and r.. = |b.||, where 
 
 b. = {a.., a. . . , ..... a .}, i = 1. .... n 
 
 i l n' i+l,i' ' n,i J ' ' 
 
 (2) Before starting this step, it is of importance to prove 
 that Q is unique if A is non-singular. 
 
 Proof: Suppose that 
 
 A = Q l R l = Q 2 R 2 
 
 therefore, if A is non-singular, R and R p are non-singular 
 
 R l R ~2 = Q l" L Q 2 = Q l Q 2 
 
 therefore, (R, R ? ) is unitary also, hence, (R, Rp ) '" = (R-, Rp ) (3 '13) 
 
 since R R„ is an upper-triangular matrix, Eq. (3'13) is satisfied if 
 
 -1 -1 . 
 
 and only if R R is a diagonal matrix, therefore, R R ? = I, i.e. 
 
 R 1 = R 2 , hence Q = Q . 
 
 - 27 - 
 
Step (l) is completed by post multiplying (R ) from Eq. (3-12) 
 
 k 
 
 by Q = N x N 2 ... N n to get A = R k Q, and the process is repeated all 
 over again until all the elements a. . , i > j of the matrix A vanish if 
 the eigenvalues of A are distinct. Similar to the LR-algorithm the 
 process may be looked upon as follows, 
 
 if P k = Ql Q 2 ...Q k 
 
 therefore, 
 
 S k " R k R k-1 • • • R 2 R l 
 
 P k S k = Ql Q 2 ... Q^ (Q^) R^ ... R 2Rl 
 
 = A l Q l 9 2 ••■ \J\-1 \-l> \-2 ••• R 2 R l 
 
 and so forth, finally 
 
 T k ~k " XJi l 
 
 P. S v = (Aj k ( 3 .1U) 
 
 Since P, is unitary and S is an upper -triangular, therefore the unitary- 
 triangular decomposition of a non-singular matrix is unique, according 
 to the lemma proved at the beginning of step (2). Francis (1961) proved 
 the convergence of P S = (A) and showed that, "If any non-singular 
 matrix A has eigenvalues of distinct modulus, then under the QR trans- 
 formation the elements below the principal diagonal tend to zero, the 
 moduli of those above it tend to fixed values, and the elements on the 
 principal diagonal itself tend to the eigenvalues." 
 
 However, for eigenvalues of the equal modulus, it can be shown 
 that the matrix A, becomes split into independent principal submatrices 
 coupled only in so far as the eigenvectors are concerned. Usually, each 
 of these submatrices is associated with groups of eigenvalues of the same 
 modulus. 
 
 28 - 
 
(3) The amount of computations that the method involves in 
 
 3 
 one iteration is proportional to n . However if the matrix is first 
 
 reduced to the upper Hessenberg form by any stable unitary transformation 
 
 technique, i.e., Householder's method, the amount of computations become 
 
 2 
 proportional to n . Reducing the matrix to the upper Hessenberg form 
 
 has some advantages, 
 
 (i) If any element a. . of the matrix A is zero, the 
 matrix can be partitioned at this element and the 
 rest of the submatrices can be operated on separ- 
 ately so far as the eigenvalues are concerned. 
 
 (ii) The upper-Hessenberg form is preserved under the 
 QR-transformation. 
 
 (iii) If an almost triangular matrix A is non-singular 
 
 and has non-zero elements, then under the QR-trans- 
 formation the principal diagonal elements of A 
 converge to the eigenvalues in order of size, if 
 they are not equal. 
 
 (iv) For distinct eigenvalues, the elements (a ) i > J, 
 
 ij k 
 
 K k 
 below the principal diagonal converge to zero as , i,. . This is the 
 
 ^ \. } 
 3 
 
 same rate of convergence as the LR-transformation which is rather slow. 
 
 Taking a closer look at the process where, for example, (a ) approaches 
 
 the exact value of the eigenvalue \ , it is found that 
 
 t> n' 
 
 e k+l = € k * ^ 
 
 Vi 
 
 where e k = (a^ - ^ and \ < 1 
 
 Vl 
 
 To speed up convergence, the origin of all eigenvalues is 
 
 r the value t before 
 3 k 
 
 after the iteration, therefore 
 
 shifted by the value £ before an iteration and then shifted back again 
 
 - 29 - 
 
n - b k 
 
 € k+l € k * X - t 
 
 n-1 s k 
 
 which is quite an improved rate of convergence for t, , = \ . 
 
 k n 
 
 Therefore, the generalized QR-transformation becomes, 
 
 Vi = \ \\ 
 
 where q* ( A)£ - E l)_ = Ek 
 
 In Section 3 methods of choosing the origin shifts (\ are 
 
 k 
 
 discussed for "both distinct and equal modulus eigenvalues especially 
 conjugate complex pairs. Once K is found by sufficient accuracy, the 
 order of the matrix may be reduced by omitting the last row and column, 
 
 3.3 Numerical Solution for QR Algorithm 
 
 The numerical solution for Q-R algorithm can be divided into 
 two major steps, viz., 
 
 (a) Applying Householder 's technique to reduce a matrix 
 to an upper Hessenberg form. 
 
 (b) Using unitary transformation with the help of origin- 
 shift to iterate, i.e., 
 
 A« . A 
 
 A (K + 1) . q (K)t A (K) Q (K) 
 
 ~(K)t . (K) produces an upper-diagonal matrix, 
 where Q, A 
 
 Since the Householder's technique has been discussed 
 in detail already, therefore we will only concentrate on part (b). 
 
 - 30 - 
 
The very first step in part (b) is to determine the value 
 of the origin shift. To find this value, we first find the roots 
 of the last principal 2x2 submatrix and then distinguish the follow- 
 ing two cases. 
 
 (l) Roots are real: choose one that differs least from 
 
 the last diagonal element a v ' (before kth iteration) 
 
 / , \ nn i \r\ 
 
 and call this root X. , then compare a. with the 
 previous one, viz. a/ '(assuming aA '=0) if 
 
 1 
 
 „(io . ^(k-D 
 
 IW 
 
 <2 
 
 set E^ = K^\ otherwise set E^ = where E^ 
 is the value of origin shift for kth iteration. 
 
 (2) Roots are complex: we do one of two things, 
 
 (i) Set A. = |r|, where |r| is the modulus 
 
 of the complex roots and then compare with 
 
 (ii) Do double origin-shift by setting A. = R and 
 
 * ( k ) ^ -+v, 4.1, • \ ( k - X ) 
 
 A. = R, compare with the previous A. 
 
 and \ p , (assuming a. = 0, a. ' = 0; . 
 
 After the value of origin shift has been found, we subtract 
 all diagonal elements from this value and then start iteration; add 
 the value back after each iteration, until the last subdiagonal 
 element -* 0, then deflate the matrix by taking out the last row 
 and column, the last diagonal element is one of the eigen-values. 
 
 In doing QR iteration, we simply produce a series of 2 x 2 
 submatrices to operate on two rows at a time on matrix A and eliminate one 
 subdiagonal element at each time during the pre-multiplication. 
 operate on two columns at a time during the post-multiplication accord- 
 ing to the following algorithm: 
 
 31 - 
 
(i) PRE-multiply: 
 
 fi 
 
 M- V 
 
 N l 
 
 1 
 J 
 
 
 r r r 
 
 r 
 
 r- - 
 
 - 
 
 - r 
 
 
 r r 
 
 r 
 
 r- - 
 
 - 
 
 - r 
 
 
 P 
 
 y 
 
 y- - 
 
 - 
 
 - y 
 
 
 a 
 
 X 
 
 x- - 
 
 - 
 
 - X 
 
 
 
 a 
 
 a- - 
 
 - 
 
 - a 
 
 1 
 
 
 
 a- - 
 
 V 
 
 a 
 
 - a 
 
 a a 
 a a 
 
 r r 
 
 r- - 
 
 - 
 
 - 
 
 r r 
 
 r r 
 
 r- - 
 
 - 
 
 - 
 
 r r 
 
 k 
 
 y*y' 
 
 - 
 
 - 
 
 yy 
 
 
 x'x' 
 
 - 
 
 - 
 
 x'x' 
 
 
 a a- 
 
 - 
 
 - 
 
 a a 
 
 
 a- 
 
 V 
 
 - 
 
 a a 
 
 
 
 \ 
 
 a 
 
 a a 
 a a 
 
 where the elements r, r' and k are those of the final triangular matrix 
 and the elements a, x, and 0! are in their original state. The elements 
 in the two rows that are changed are given by: 
 
 k = (|p| 2 + \cc\ 2 )2 
 x ' = ux - vy 
 y'= uy - vx 
 u = p/k 
 v = a/k 
 
 (ii) POST-multiply: 
 
 a x y r r- 
 
 a x y r r- 
 
 a x y r r- 
 
 k r r- 
 
 r r- 
 
 N _ 
 
 x r 
 
 f\ 
 
 \1 -V 
 V [1 
 
 \1 
 
 a a x'y'r r - 
 
 a a x'y'r r - 
 
 a x'y'r r - 
 
 a p r r - 
 
 r r - 
 
 r - 
 
 - r 
 
 - r 
 
 - r 
 
 - r 
 
 - r 
 
 - r 
 
 Where the elements r, y and k are so far unaltered from the triangular 
 state, and the elements a, x' and a are in their final state. The two 
 columns that are changed are given by: 
 
 - 32 - 
 
a = vk 
 P = jlk 
 
 x' = \ix + vy 
 
 y T = uy - vx 
 
 The higher level language program with a. detailed flow 
 chart for single origin-shift is in section 3«^+ and 3«5. 
 
 - 33 - 
 
3«^ High-Level Language Program 
 
 START: 
 
 L00P1 : 
 
 #BEGIN 
 
 #LABEL 
 
 #ARRAY 
 
 #REAL 
 
 #REAL 
 
 #SET 
 #READ 
 
 #IF 
 
 START, LOOP1, L00P2, LOOP3, LIP, L2P, L3P, L^P, 
 LOPM, LOPE; 
 
 #REAL #FLOAT #SHORT #SKWD 
 
 A[64,64], P[64,64], IM[6^,64], UT[6U,l], TEM[6U,l], 
 
 V[61+], \j[6k]; 
 
 #FLOAT #SHORT 
 
 S2, S, T, El, TEMPI, XI, X2, LI, L2, E2, E, K; 
 
 #FIXED 
 
 J, H, L, Kl, N, M, I; 
 
 HE, LL, KK, II, IJ, JI; 
 
 A; 
 
 J - 2; 
 
 HH: : [J, J+l, , Gh] ; 
 
 LL:: [l, 2, ...., J-l] ; 
 
 KK: : [J+l, J+2, , Gk] ; 
 
 S2 «- #FOR (H) #SIM (HH) #DO #SUM(A[H, J-l] 1 ' 2) ; 
 
 S «- SQRT(S2); 
 
 T «- S2 - A[J,J-1] x S; 
 
 #FOR (L) #SIM (LL) #DO UT[L,1>0.; 
 
 UT[J,1] - A[J,J-1] - S; 
 
 #FOR (Kl) #SIM (KK) #DO UT[K,l] «-A[K,J-l]; 
 
 P «- IM - TRANSPOSE (UT) x UT/T; 
 
 A <- P x A; 
 
 A - A x TRANSPOSE(P); 
 
 J «- J+l; 
 
 J < 63 #THEN #GOTO LOOPl; 
 
 N *- 6h; 
 
 - 3h - 
 
L00P2 
 
 L00P3 
 
 LIP 
 
 L2P 
 
 LOPM: 
 
 II: : [1, ?-, , N] 
 
 El «- 0.; 
 
 TEMPI - SQRT((A[N,N] + A[N-l,N-l] )t 2-k x(A[N-l,N-l] x 
 
 A[N,N] - A[N,N-1] x A[N-1,N])); 
 XI «- .5 x (A[N,N] + A[N-1,N-1] + TEMPI); 
 X2 - .5 x (A[N,N] + A[N-1,N-1] -TEMPI); 
 LI «- ABS(A[N,N] - XI); 
 L2 «- ABS(A[N,N] - X2); 
 #IF LI < L2 #THEN #D0 #BEGIN 
 
 E2 <- LI; 
 #GOTO LIP; 
 #END; 
 E2 - L2; 
 
 #IF ABS((E2-E1)/E2) < .5 #THEN #D0 #BEGIN 
 
 E <- E2; 
 El - E; 
 #GOTO L2P; 
 #END; 
 E *- 0.; 
 
 #FOR (I) #SLM (II) #D0 
 A[I,I] - A[I,I] - E; 
 M <- 1} 
 
 JI: : [M+l,M+2, ,M] ; 
 
 K ♦- SQRT(A[M,M] t2 + A[M+1,M] t2); 
 
 U[M] «- A[M,M]/K; 
 
 V[M] <- A[M+1,M]/K; 
 
 #FOR (I) #SIM (JI) #D0 #BEGIN 
 
 TEM[M,I] <- U[M] x A[M,I] - V[M] x A[M+1,I]; 
 
 A[M+1,I] - U[M] x A[M+1,I] + V[M] x A[M,I]; 
 
 A[M,I] *- TEM[M,I]; 
 
 #MD; 
 
 A[M,M] «- K; 
 
 A[M+1,M] *- 0.; 
 #IF M > N-l #THEN #GOTO L3P; 
 
 - 35 - 
 
M *~ M+l; 
 #GOTO LOPM; 
 
 L3P: M «- 1; 
 
 LOPN: IJ:: [1,2, ,M] ; 
 
 A[M+1,M] *- V[M] x A[M+l,M+l]; 
 A[ M+l, M+l] <- U[M] x A[M+l,M+l]; 
 #FOR (I) #SIM (IJ) #DO #BEGIN 
 
 TEM[I,M] *- U[M] x A[I,M] + V[M] x A[l,M+l]j 
 
 A[I,M+1] - U[M] x A[I,M+1] + V[M] x A[I,M]; 
 
 A[I,M] <- TEM[l,M]j 
 
 #EM); 
 
 M <- M+l; 
 #IF M > N #THEN #GOTO L^P; 
 
 #GOTO LOPN; 
 
 L4P: #FOR (I) #SIM (II) #DO 
 
 A[I,I] - A[I,I] + E; 
 #IF ABS(A[N,N-1]) > .OCOOOl #THEN #GOTO L00P3; 
 #PRINT A[N,N]j 
 N <- N-l; 
 #IF N > 1 #THEN #GOTO L00P2; 
 #PRINT A[l,l]; 
 #END 
 
 - 36 - 
 
3-5 Flow Chart 
 
 TO GET UPPER HESSENBERG FORM 
 
 H «- 2 
 
 2 
 S = S2 
 
 ? 
 
 2K~= T 
 
 m — 
 
 S 2 -Z [A(j,(H-l))] 2 
 J=H 
 
 2K _ s 2 - A(H,H-l)xS 
 
 L— x,2.f • • • jH-1 
 
 UT(L) «- 
 UT(H) «= A(H,H-l)-S 
 
 K=H+1,H+2,...,N 
 
 UT(K) - A(K,H-1) 
 
 I is 64x6*4 
 
 identity matrix 
 
 P^I -TRANS POSE ( UT ) xUT/ 2K' 
 
 A <- PxA 
 A *- AxTRANSPOSE(P) 
 
 H «- H+l 
 
 yes 
 
 no 
 
 © 
 
 - 37 - 
 
FIND ORIGIN-SHIFT VALUE 
 
 III) • 
 
 © 
 
 N «- 6k 
 
 E1--0 
 
 Xl4(a +a . ,+SQRT((a +a , _ Y 
 2 X n,n n-l,n-l ^ vv n,n n-l,n-l' 
 
 -i+(a . *a -a *a . ))) 
 
 n-l,n-l n,n n,n-l n-ljrr'' 
 
 X2«-i(a +a . ,-SQRT((a +a _ . )' 
 2 V n,n n-l,n-l vv n,n n-1,11-1' 
 
 -Ma i n* a -a n *a . ))) 
 n-l,n-l n,n n,n-l n-l,n 
 
 E - E2 
 
 El «- E 
 
 LI 
 L2 
 
 a -XI 
 
 a -X2| 
 n,n ' 
 
 yes 
 
 E2 «- L2 
 
 yes 
 fc 
 
 E - 
 
 A(l,l)-E-A(l,l) 
 
 ) E2 <- LI 
 
 a =A(N,N) 
 n,n v ' ' 
 
 1=1,2,..., N 
 
 - 38 - 
 
PREMULTIPLY: GET UPPER DIAGONAL MATRIX 
 
 M *- 1 
 
 K-SQRT(A 2 (m,m)+A 2 (m+l ; m)) 
 
 U(M) 
 V(M) 
 
 A(M;M)/K 
 
 a(m+i,m)/k 
 
 J=M+l,M+2, ,N 
 
 a(m,j)«-u(m)*a(m,j)-v(m)*a(m+i,j) 
 
 A(M+1,J) <-U(M)*A(M+l,J)+V(M)*A(M,j) 
 
 A(M,M)*^ 
 
 a(m+i,m)*o 
 
 Yes 
 
 No 
 
 * M^M+1 
 
 - 39 
 
POST -MULTI PLY AND TESTING 
 
 © 
 
 
 v 
 
 
 M - 1 
 
 
 
 
 
 V 
 
 
 a(m+i,mK-v(m)*a(m+i,m+i) 
 
 a(m+i,m+i)«-u(m)*a(m+i,m+i) 
 
 
 \ 
 
 1 
 
 
 a(j,mKu(m)*a(j,m)+v(m)*a(j,m+i) 
 a(j,m+i)^u(m)*a(j,m+i)-v(m)*a(j,m) 
 
 
 \ 
 
 ' 
 
 
 M «- M+l 
 
 i 
 
 
 
 No 
 
 M > N 
 
 Yes 
 
 A(l,l)-A(l,l)+E 
 
 Yes 
 
 N*-N-l 
 PRT.A(N,N) 
 
 J=l,2, ...,M 
 
 1=1,2,..., N 
 
 n < r 
 
 No 
 
 Yes 
 
 PRT.A(1,1) 
 
 kO - 
 
k. References 
 
 Bodewig, E., 1959j "Matrix Calculus". Interscience Publishers, 
 New York; North-Holland Publishing Company, Amsterdam. 
 
 Francis, J. G. F., I96.I, .1962, "The QR Transformation, Part I 
 and II." Computer J. 1+, 265-271, 332-3^5 • 
 
 Wilkinson, J. H-, 1960a, "Householder's Method for the Solution 
 of the Algebraic Eigenproblem. " Computer J. 3> 23-27- 
 
 Wilkinson, J. H. , 19&5, "The Algebraic Eigenvalue Problem." 
 Clarendon Press, Oxford. 
 
 Young, David M. Jr., i960, "Numerical Methods for Solving 
 Problems in Linear Algebra." The University of Texas 
 Computer Center, Austin. 
 
 - kl - 
 
APPENDIX A 
 
 High-Level Language Program For Finding MAX 
 
 #BEGIN 
 
 #LABEL START, LOOP, LEND; 
 
 #ARRAY #REAL #FLOAT #SHORT #SKWD 
 
 A[6k,6k], B[6k, 1], Cl6k, 1]; 
 #REAL #FLOAT #SHORT 
 
 LIM, CMAX; 
 #REAL #FIXED 
 
 i; 
 
 #SET II, JJ; 
 
 START: #READ A 
 
 II:: [1,2,...., 6k]; 
 
 #FOR (I) #SIM (II) #D0 
 B[I,1] - 1. ; 
 LIM ^ 1. ; 
 
 LOOP: C - A x B; 
 
 CMAX «- C[l,l]; 
 
 JJ:: [1,2,. ...,63]; 
 
 #FOR (I) #SEQ (JJ) #DO #BEGIN 
 
 #LT ABS(CMAX) - ABS(C[I+1, 1] ) > #THEN #GOTO LEND 
 
 #ELSE CMAX *- C [1+1,1]; 
 
 LEND: #END; 
 
 #IF ABS(ABS(CMAX) - ABS(LIM)) > .00001 #THEN #D0 #BEGIN 
 LIM «- CMAX; 
 
 #F0R (I) #SIM (II) #D0 
 B[I,1] - C[I,1]/LIM; 
 #G0T0 LOOP; 
 #END #ELSE 
 #PRINT CMAX; 
 
 #END 
 
 - k2 
 
READ IN A 
 
 B[I,1] - I- 
 
 < A is N x I (6k here) 
 
 < 1=1,2,... ,M (64 here) 
 
 LIM «- 1. 
 
 C <- A x 
 
 FIND MAX(C) -> CMAX 
 
 yes 
 
 B[I,1] = C[I,1]/LIM 
 
 PRINT CMAX 
 
 STOP 
 
 ^3 
 
APPENDIX B 
 
 QP-Algorithm With Double Origin Shift 
 
 When a matrix has a complex-pairs of eigen-values, the double 
 origin shift is needed. In this scheme, we combine two iterations into 
 one with a pair of origin shifts since single iteration may not be 
 
 sufficient enough due to the existence of a complex -pair eigen-values. 
 
 th 
 So at k stage, we will have 
 
 A ( k+2 ) = Q (k+i)* Q (k)» A (k) Q (k) Q ( k+ i) 
 
 Let W = Q (k) Q (k+1) and P = (A (k) - E (k) l) (A (k) - E (k+l) l) where 
 
 E ' and E are the pair of origin shifts, such that W*T gives an 
 upper triangular matrix. W* here is unitary and W = N N ... N 
 
 N. = 
 
 i 
 
 form: 
 
 •i-1 
 
 M. 
 
 l- 
 
 in the previous discussion. T has the following 
 
 x x x x x x 
 
 x x x x x x 
 
 x x x x x x 
 
 x x x x x 
 
 x x x x 
 
 X X X 
 
 N 
 
 X 
 
 x XXX 
 
 Let the first column be r then N should be defined so that 
 
 N r = 
 1 1 
 
 r ll l- e l 
 
 - kk - 
 
Hence, N A v J N 
 
 takes the form: 
 
 X X X X X X X 
 
 X X X X X X X 
 
 X X X X X X X 
 
 X X X X X X X 
 
 X X X X 
 
 X X X 
 
 V 
 
 \ 
 
 \ 
 
 XXX 
 
 X X 
 
 Where the new matrix has few more elements then before, and this is 
 not what we want to do, therefore the general technique will be as 
 follows: 
 
 (l) Perform an initial transformation (N a 
 
 ( v\ 
 A with double origin shifts. 
 
 (k) 
 
 N ) on 
 
 (2) Reduce the resulting matrix to almost triangular form 
 by a method, such as Householder's to obtain A 
 
 A typical stage in the iteration will take the forms: 
 Before Householder's transformation: (ii) After Householder's transformation: 
 
 XXX X'Y'Y'Y Y Y 
 
 XXX X'Y'Y'Y Y Y 
 
 X X X'Y'Y'Y Y Y 
 
 X X'Y'Y'Y Y Y 
 
 Z'Z'Z'Y Y Y 
 
 Z'Z'Z'Y Y Y 
 
 " Z'Z'Z'T T T 
 
 T T T 
 
 T T 
 
 X 
 
 X 
 
 X 
 
 Y 
 
 Y 
 
 Y 
 
 Y 
 
 Y 
 
 Y 
 
 X 
 
 X 
 
 X 
 
 Y 
 
 Y 
 
 Y 
 
 Y 
 
 Y 
 
 Y 
 
 
 X 
 
 X 
 
 Y 
 
 Y 
 
 Y 
 
 Y 
 
 Y 
 
 Y 
 
 
 
 z 
 
 Z 
 
 Z 
 
 Y 
 
 Y 
 
 Y 
 
 Y 
 
 
 
 z 
 
 Z 
 
 Z 
 
 Y 
 
 Y 
 
 Y 
 
 Y 
 
 
 
 z 
 
 z 
 
 Z 
 
 T 
 T 
 
 T 
 T 
 
 T 
 
 T 
 T 
 T 
 T 
 
 T 
 T 
 T 
 T 
 
 where 3 rows and 3 columns being affected. X and T are those of the 
 
 • • + ■ 1 * «- 1 i. i «(k) A (k+2) 
 
 initial and final matrices A and A res 
 
 Y's and Z's are changed by the transformation. 
 
 - ^5 
 
(k) 
 Since the first transformation of A "by N is similar to 
 
 the subsequent transformations, we can set up the initial condition 
 
 for the iteration by finding r , r_ , and r given by 
 
 r n = a n (a n " a) + a i2 ' a 2i + p 
 
 r 21 = a 21 (a ll + a 22 " a) 
 
 r 31 == a 21 * a 32 
 
 where 
 
 a = E 1 +■ E 2 
 
 E 1 • E 2 
 
 the transformation matrices of Householder's method are of the form: 
 
 N. - N.* = I - 2t.t.* 
 
 11 11 
 
 t. is such that ti = 1 
 l i i i i 
 
 t. . = for j < i and j > i+2 
 
 1 Z 10 i 
 
 t. . 
 
 J 20 
 
 i,i+l 2kt 
 
 li 
 
 t. . 
 
 J 30 
 
 i,i+2 2kt, 
 
 li 
 
 2*i 
 
 where k =(sign Z 1Q ] • (Z 1Q + Z 2Q + Z^ )2 
 
 where Z 
 
 rs 
 
 Z 10 Z ll Z 12 
 
 7 7 7 
 
 20 21 22 
 
 7 7 7 
 
 30 31 32 
 
 corresponds to the Z's of the figure above, 
 
 - k6 
 
In order to reduce the number of multiplication in the Householder's 
 
 transformation let H. = [0, . . . ,0,1.x. ,X o ,0, . . . ,0] = t — • t. 
 
 l 1 2 t.. l 
 
 11 
 
 obtaining N. = I - (l + cp) H. H. 
 
 i 11 
 
 Z 10 Z 20 Z 30 
 
 * = T" ' X l = k?Z^ ' X 2 = k^ and a11 < 1 
 
 Hence at i stage, the operation on columns a., a. _ and a. „ is [a., a. ...a. _] 
 
 r -, "I i i+l i+2 i' l+l' H-2 J 
 
 [a i , a 1+1 , a i+2 ].[I - (l + cp) 
 
 {a. + X a + X 2 a. +2 ] 
 
 therefore, a. = a. - n 
 li 
 
 X l 
 X 2 
 
 [1, X 1 , X 2 ]] Let r) = (1 +■ cp) . 
 
 Vl = Vl - V 11 
 
 a i + 2 = a i + 2 " X 2^ 
 
 Similarly for the row operation. Z_ and Z~ n are eliminated and replace 
 Z n by - k in each transformation. 
 
 For origin shift, at k stage, we need to consider the last 
 2x2 sub-diagonal matrix. First, find the two roots of this 2x2 matrix, 
 
 X , X which differ least from the last two diagonal elements 
 
 (k-2) 
 
 respectively. Then compare with the previous two roots X and 
 
 X , and calculate 
 
 (k) (k-2) (k+1) (k-1) 
 
 I- ~r-^ 1 and I- / '■ v 1 
 
 I JJk) ' ' jjk+l) I 
 
 if both > | then set E^ = e^ 1 ' = 0; if both < | , then set E^ k ' = 
 
 . (k) „(k+l) . (k+l) ,, . . . , , _(k) , „(k+l) , , .. -, 
 \ v , E v ' = \ s ' ; otherwise set both E v andE v ' to be the real 
 
 part of either X^ ' or X whichever corresponds to the quantity 
 
 < p . Then iterate again until either or both of the last two sub- 
 diagonal elements are near zero, then deflate the matrix accordingly; 
 record the eigen-value or eigen-values; find new origin shifts; iterate 
 again in the same way until all eigen-values are found. 
 
 - ^7 - 
 
''TRANQUIL FOR Q-R ALGORITHM WITH DOUBLE ORIGIN SHIFT 
 
 #BEGIN 
 
 #LABEL START, LOOPl, L00P2, LOOP3, ROWOP, CLMOP; 
 
 #ARRAY #REAL #FLOAT #SHORT #SKWD 
 
 k[6k,6k], v[6k,6k], m[6k,6k], ur[6i+,i]; 
 
 START: 
 
 LOOP1: 
 
 L00P2: 
 
 #REAL #FLOAT #SHORT 
 
 S2, S, T, El, E2, TEMPI, XI, X2, LI, L2, EE1, EE2, 
 ESUM, EPDT, Rl, R2, R3, K, ALP, PI, P2, AN, RT1, RT2; 
 #REAL #FIXED 
 
 J, H, Kl, L, N, I; 
 #SET HH, LL, KK, II, JJ; 
 #READ A; 
 
 J- 2; 
 HH::[J,J+l,...,6l4]; 
 LL::[l,2,..,,J-l]; 
 KK:: [J+l, J+2, . . . ,6k] ; 
 
 S2 «- #OR (H) #SIM (HH) #DO #SUM(A[H, J-l]t 2) ; 
 S - SQRT(S2); 
 T <- S2 - A[J,J-l]* S; 
 #FOR (L) #SIM (LL) #DO UT[L,1] «- 0.; 
 UT[J,1] - A[J,J-1]- S; 
 
 #FOR (Kl) #SIM (KK) #D0 UT[Kl,l] *- A[Kl,J-l]; 
 P «- IM-TRANSPOSE (UT) x UT/T; 
 A <- PxA; 
 
 A ^ A x TRANSPOSE (P); 
 J <- J+l; 
 #IF J < 63 #THEN #GOTO LOOPl; 
 N *- 6k; 
 
 II::[1,2,...,N] 
 El <- E2 <- 0.; 
 
 k8 
 
L00P3: TEMPI «- SQRT((A[N,N] + A[N-l,N-l])t 2-l+*(A[N-l,N-l]* 
 
 A[N,N] - A[N,N-l]* A[N-1,N])); 
 XI *■ .5* (A[N,N] + A[N-l,N-l] + TEMPI); 
 X2 <- .5* (A[N,N] + A[N-l,N-l] •- TEMPI); 
 LI «- ABS (A[N,N] - XI); 
 L2 «- ABS (A[N,N] - X2); 
 #LE LI < L2 #THEN #D0 #BEGIN 
 
 EE1 «- XI; 
 
 EE2 *- X2; #END #ELSE #D0 #BEGIN 
 
 EE1 «- X2; 
 
 EE2 - XI; #END; 
 #IF ABS((EE1-E1)/EE1) < .5 #THEN 
 
 El «- EE1 #ELSE 
 
 El <- 0.; 
 #IF ABS((EE2-E2)/EE2) < .5 #THEN 
 
 R0WOP: #FOR (j) #SIM(jj) #D0 #BEGIN 
 
 #LF I < N - 2 #THEN 
 
 AN <- ALP* (A[I,J] f PI* A[I+1,J] + P2* A[I+2,J]) 
 
 #ELSE 
 
 AN <- ALP* (A[I,J] + PI* A[I+1,J]); 
 
 A[I,J] <- A[I,J] - AN; 
 
 A[I+1,J] «- A[I+1,J] - PI* AN; 
 #IF I < N-2 #THEN 
 
 A[I+2,J] <- A[I+2,J] - P2* AN; 
 CLMOP: #IF I < N-2 #TKEN 
 
 AN - ALP* (A [J, I] +■ PI* A[J,I+1] + P2* A[J,I+2]) 
 
 #ELSE 
 
 AN +- ALP* (A[J,I] + PI* A[J,I+1]); 
 
 A[J,I] <- A[J,I] - AN; 
 
 A[J,I+l] «- A[J,I+1] - PI* AN; 
 #IF I < N-2 #THEN 
 
 A[J,I+2] *- A[J,I+2] - P2* AN; 
 #END; 
 
 - k9 
 
#BEGIN 
 
 #IF I < N-3 #THEN #D0 #BEGIN 
 
 AN «- ALP * P2* A[l+3,I+2]j 
 
 A[I+3,l] - - AN; 
 
 A[I+3,I+1]«--P1* AN; 
 
 A[I+3,I+2] - A[l+3 7 I+2] - P2* AN; 
 
 #END; 
 
 E2 <- EE2 #ELSE 
 
 E2 - 0.; 
 
 ESUM <- E1+E2; 
 
 EPDT <- E1*E2; 
 #0R (I) #3EQ (II) #D0 
 
 #LE 1=1 #THEN #D0 #BEGIN 
 
 Rl «- A[l,l]* A(A[l,l] - ESUM) + A[l,2] * A[2,l] + EPDT; 
 
 R2 «- A[2,l]* (A[2,2] - ESUM + A[l,l]); 
 
 R3 - A[2,l]* A[3,2]; 
 
 A[3,l] *" 0; # EN D #ELSE #D0 #BEGIN 
 
 Rl = A[I,I-l]; 
 
 R2 = A[I+1,I-1]; 
 #IF I=N-1 #THEN R3 - #ELSE R3 <- A[I+2,I-l]; 
 
 #END; 
 #IF Rl > #THEN 
 
 K <- SQRT (Rlt 2 + R2t 2 + R3t 2 ) #ELSE 
 
 K <- - SQRT (Rlt 2 + R2t 2 + R3 t 2); 
 #IF K ^ #THEN #D0 #BEGIN 
 
 ALP <- Rl/K + 1.; 
 
 PI <- R2/(RH-k); 
 
 P2 - R3/(Rl+k); #END #ELSE #D0 #BEGIN 
 
 ALP <- 2; 
 
 PI <- P2 <- 0.; #END; 
 
 JJ:: [1,1+1, ,N]; 
 
 #IF ABS(A[N,N-1]) > .000001 #THEN #D0 #BEGIN 
 
 #IF ABS(A[N-l,N-2]) > .000001 #THEN #G0T0 #L00P3; 
 
 TEMPI «_ SqRT((A[N,N] + A[N-l,N-l])t2-U* (A[N--l,N-l]* 
 A[N,N] - A[N,N-l]* A[N-1,N])); 
 
 50 
 
#END; 
 
 #END 
 
 RT1 *- .5* (A[N,N] + A[N-1,N-1] + TEMPI); 
 
 RT2 «- .5* (A[N,K] + A[N-1,N-1] - TEMPI); 
 #PRINT RT1, RT2; 
 
 N - N-2; #END #ELSE #D0 #BEGIN 
 #IF ABS(A[W-l,N-2]) < .000001 #THEN #D0 #BEGIN 
 
 #PRINT A[N,N], A[N-l,N-l]; 
 
 N <- N-2; #END #ELSE #D0 #BEGIN 
 
 #PRINT A[N,N]; 
 
 N «- N-l; #END; 
 #END; 
 
 #IF N < 1 #THEN #G0T0 L00P2; 
 #ST0P; 
 
 - 51 
 
TO GET UPPER HESSENBERG FORM 
 
 H «- 2 
 
 2 
 S = S2 
 
 2 
 
 2K 5 T 
 
 2 
 
 "FT 
 
 I [A(J,(H-1).)] ! 
 
 J=H 
 
 2K 2 *- S 2 - A(H,H-l)xS 
 
 L=1,2,...,H-1 * 
 
 K=H+l,H+2,...,N-> 
 
 UT(L) *- 
 UT(H) <- A(H,H-1)-S 
 
 UT(K) <- A(K,H-1) 
 
 I is 6kx6k 
 Identity 
 
 P <- I -TRANSPOSE (UT)xUT/2K' 
 
 A - PxA 
 A «- AxTRANSPOSE(P) 
 
 H <- H + 1 
 
 yes 
 
 52 - 
 
FIND ORIGIN SHIFT 
 
 d> 
 
 El «- « 
 
 E2 «- 
 
 N «- 6U 
 
 El <- E2 «- 
 
 Find two roots of the 
 
 last 2x2 matrix call 
 
 XI, X2 
 
 LI 
 L2 
 
 a -XI 
 n,n 
 
 a ' -X2 
 
 ESUM «- E1+E2 
 EPDT «- E1*E2 
 
 EE1 - XI 
 
 EE2 «- X2 
 
 El «- EE1 
 
 E2 «- EE2 
 
 - 53 - 
 
SET UP AND DO ROW AND COLUMN OPERATIONS 
 
 © 
 
 i=2,3,-..,N-2-» 
 
 i=N-l 
 
 Rl=a. . . 
 
 R2=a, 
 
 Rl=a (a -ESUM)+a *a +EPI 
 XfX 1,1 ±,d ^->1 
 
 !,l (a 2 r 
 
 -ESUM+a ) 
 
 R>a 2,l * a 3,2 
 
 I 2 2 2 
 
 K = -jRl +R2 +R3 
 
 yes 
 
 ALP «- Rl/K + 1 
 PI - R2/(R1+K) 
 P2 «- R3/(R1+K) 
 
 ROW OPERATION 
 
 COLUMN OPERATION 
 
 2 2 2 
 = jRl +R2 +R3 
 
 ALP 4- 2 
 PI <- P2 ♦- 
 
 5^ - 
 
END TEST 
 
 PRT.A(N,N) 
 N *- N - 1 
 
 no 
 
 III 
 
 no 
 
 yes 
 
 STOP 
 
 Find two roots of the 
 
 last 2x2 matrix call 
 
 RT1, RT2 
 
 PRT. RT1, RT2 
 N <- N - 2 
 
 - 55 - 
 
AUG ] 
 
JUN 2 1969