490 Chapter 11.Eigensystems of the matrix and look for the next eigenvalue.Alternatively,an-1.n-2 may become negligible.In this case the eigenvalues of the 2 x 2 matrix in the lower right-hand corner may be taken to be eigenvalues.We delete the nth and(n-1)st rows and columns of the matrix and continue. The test for convergence to an eigenvalue is combined with a test for negligible subdiagonal elements that allows splitting of the matrix into submatrices.We find the largest i such that ai.is negligible.Ifi=n,we have found a single eigenvalue.If i=n-1,we have found two eigenvalues.Otherwise we continue the iteration on the submatrix in rows i to n(i being set to unity if there is no small subdiagonal element). 三 After determining i,the submatrix in rows i to n is examined to see if the product 81 of any two consecutive subdiagonal elements is small enough that we can work with an even smaller submatrix,starting say in row m.We start with m =n-2 and decrement it down to i+1,computing p,q,and r according to equations(11.6.23) with 1 replaced by m and 2 by m+1.Ifthese were indeed the elements of the special "first"Householder matrix in a double QR step,then applying the Householder matrix would lead to nonzero elements in positions(m+1,m-1),(m+2,m-1). and(m +2,m).We require that the first two of these elements be small compared with the local diagonal elements am-1.m-1,amm and am+1.m+1.A satisfactory approximate criterion is 9 lam,m-1(lg+r)<pl(lam+1.m+1+lamml+lam-1.m-1) (11.6.26) Very rarely,the procedure described so far will fail to converge.On such matrices,experience shows that if one double step is performed with any shifts 9起N that are of order the norm of the matrix,convergence is subsequently very rapid. Accordingly,if ten iterations occur without determining an eigenvalue,the usual shifts are replaced for the next iteration by shifts defined by kg+kg+1=1.5×(0an,n-1+lan-1,n-2l0 (11.6.27) ksks+1=(lan,n-1+lan-1.n-2)2 The factor 1.5 was arbitrarily chosen to lessen the likelihood of an"unfortunate" choice of shifts.This strategy is repeated after 20 unsuccessful iterations.After 30 unsuccessful iterations,the routine reports failure. 61 Numerica 10621 4310 The operation count for the QR algorithm described here is~5k2 per iteration, where k is the current size of the matrix.The typical average number of iterations per (outside Recipes eigenvalue is~1.8,so the total operation count for all the eigenvalues is~3n3.This estimate neglects any possible efficiency due to splitting or sparseness of the matrix. North Software. The following routine hor is based algorithmically on the above description, in turn following the implementations in [1,2]. #include <math.h> #include "nrutil.h" void hgr(float **a,int n,float wr[],float wi[]) Finds all eigenvalues of an upper Hessenberg matrix a[1..n][1..n].On input a can be exactly as output from elmhes 811.5;on output it is destroyed.The real and imaginary parts of the eigenvalues are returned in wr[1..n]and wi[1..n],respectively.490 Chapter 11. Eigensystems Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). of the matrix and look for the next eigenvalue. Alternatively, a n−1,n−2 may become negligible. In this case the eigenvalues of the 2 × 2 matrix in the lower right-hand corner may be taken to be eigenvalues. We delete the nth and (n − 1)st rows and columns of the matrix and continue. The test for convergence to an eigenvalue is combined with a test for negligible subdiagonal elements that allows splitting of the matrix into submatrices. We find the largest i such that ai,i−1 is negligible. If i = n, we have found a single eigenvalue. If i = n−1, we have found two eigenvalues. Otherwise we continue the iteration on the submatrix in rowsi to n (i being set to unity if there is no small subdiagonal element). After determining i, the submatrix in rowsi to n is examined to see if the product of any two consecutive subdiagonal elements is small enough that we can work with an even smaller submatrix, starting say in row m. We start with m = n − 2 and decrement it down to i + 1, computing p, q, and r according to equations (11.6.23) with 1 replaced by m and 2 by m+ 1. If these were indeed the elements of the special “first” Householder matrix in a double QR step, then applying the Householder matrix would lead to nonzero elements in positions (m + 1, m − 1),(m + 2, m − 1), and (m + 2, m). We require that the first two of these elements be small compared with the local diagonal elements am−1,m−1, amm and am+1,m+1. A satisfactory approximate criterion is |am,m−1|(|q| + |r|)  |p|(|am+1,m+1| + |amm| + |am−1,m−1|) (11.6.26) Very rarely, the procedure described so far will fail to converge. On such matrices, experience shows that if one double step is performed with any shifts that are of order the norm of the matrix, convergence is subsequently very rapid. Accordingly, if ten iterations occur without determining an eigenvalue, the usual shifts are replaced for the next iteration by shifts defined by ks + ks+1 = 1.5 × (|an,n−1| + |an−1,n−2|) ksks+1 = (|an,n−1| + |an−1,n−2|) 2 (11.6.27) The factor 1.5 was arbitrarily chosen to lessen the likelihood of an “unfortunate” choice of shifts. This strategy is repeated after 20 unsuccessful iterations. After 30 unsuccessful iterations, the routine reports failure. The operation count for the QR algorithm described here is ∼ 5k 2 per iteration, where k is the current size of the matrix. The typical average number of iterations per eigenvalue is ∼ 1.8, so the total operation count for all the eigenvalues is ∼ 3n3. This estimate neglects any possible efficiency due to splitting or sparseness of the matrix. The following routine hqr is based algorithmically on the above description, in turn following the implementations in [1,2]. #include <math.h> #include "nrutil.h" void hqr(float **a, int n, float wr[], float wi[]) Finds all eigenvalues of an upper Hessenberg matrix a[1..n][1..n]. On input a can be exactly as output from elmhes §11.5; on output it is destroyed. The real and imaginary parts of the eigenvalues are returned in wr[1..n] and wi[1..n], respectively. {