11.5 Reduction of a General Matrix to Hessenberg Form 485 required elements in a column of the matrix.Householder reduction to Hessenberg form is in fact an accepted technique.An alternative,however,is a procedure analogous to Gaussian elimination with pivoting.We will use this elimination procedure since it is about a factor of 2 more efficient than the Householder method, and also since we want to teach you the method.It is possible to construct matrices for which the Householder reduction,being orthogonal,is stable and elimination is not,but such matrices are extremely rare in practice. Straight Gaussian elimination is not a similarity transformation of the matrix Accordingly,the actual elimination procedure used is slightly different.Before the rth stage,the original matrix A =A1 has become Ar,which is upper Hessenberg in its first r-1 rows and columns.The rth stage then consists of the following sequence of operations: Find the element of maximum magnitude in the rth column below the diagonal.If it is zero,skip the next two"bullets"and the stage is done. Cam Otherwise,suppose the maximum element was in row r'. Interchange rows r'and r+1.This is the pivoting procedure.To make the permutation a similarity transformation,also interchange columns r' and r+1. 三d 令 For i=r+2,r+3....,N,compute the multiplier dir Press. ni,r+1三 0r+1, Subtract ni.+times row r+1 from row i.To make the elimination a similarity transformation,also add ni+times column i to columnr+1. A total of N-2 such stages are required. IENTIFIC When the magnitudes of the matrix elements vary over many orders,you should try to rearrange the matrix so that the largest elements are in the top left-hand corner. This reduces the roundofferror,since the reduction proceeds from left to right. Since we are concerned only with eigenvalues,the routine elmhes does not 最管 keep track of the accumulated similarity transformation.The operation count is about 5N3/6 for large N. 10621 #include <math.h> #define SWAP(g,h){y=(g);(g)=(h);(h)=y;} Numerical Recipes 43108 void elmhes(float a,int n) Reduction to Hessenberg form by the elimination method.The real,nonsymmetric matrix (outside a[1..n][1..n]is replaced by an upper Hessenberg matrix with identical eigenvalues.Rec- ommended,but not required,is that this routine be preceded by balanc.On output,the North Software. Hessenberg matrix is in elements a[i][j]with i<j+1.Elements with i>j+1 are to be thought of as zero,but are returned with random values. Ame int m,j,i; float y,xi for (m=2;m<n;m++){ m is called r+1 in the text. x=0.0; 1m; for (j-m;j<en;j++){ Find the pivot. if (fabs(a[j][m-1])>fabs(x)){ x=a[j][m-1]; 1=j;11.5 Reduction of a General Matrix to Hessenberg Form 485 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). required elements in a column of the matrix. Householder reduction to Hessenberg form is in fact an accepted technique. An alternative, however, is a procedure analogous to Gaussian elimination with pivoting. We will use this elimination procedure since it is about a factor of 2 more efficient than the Householder method, and also since we want to teach you the method. It is possible to construct matrices for which the Householder reduction, being orthogonal, is stable and elimination is not, but such matrices are extremely rare in practice. Straight Gaussian elimination is not a similarity transformation of the matrix. Accordingly, the actual elimination procedure used is slightly different. Before the rth stage, the original matrix A ≡ A1 has become Ar, which is upper Hessenberg in its first r − 1 rows and columns. The rth stage then consists of the following sequence of operations: • Find the element of maximum magnitude in the rth column below the diagonal. If it is zero, skip the next two “bullets” and the stage is done. Otherwise, suppose the maximum element was in row r
. • Interchange rows r
and r + 1. This is the pivoting procedure. To make the permutation a similarity transformation, also interchange columns r
and r + 1. • For i = r + 2, r + 3,...,N, compute the multiplier ni,r+1 ≡ air ar+1,r Subtract ni,r+1 times row r + 1 from row i. To make the elimination a similarity transformation, also add ni,r+1 times column i to column r + 1. A total of N − 2 such stages are required. When the magnitudes of the matrix elements vary over many orders, you should try to rearrange the matrix so that the largest elements are in the top left-hand corner. This reduces the roundoff error, since the reduction proceeds from left to right. Since we are concerned only with eigenvalues, the routine elmhes does not keep track of the accumulated similarity transformation. The operation count is about 5N3/6 for large N. #include <math.h> #define SWAP(g,h) {y=(g);(g)=(h);(h)=y;} void elmhes(float **a, int n) Reduction to Hessenberg form by the elimination method. The real, nonsymmetric matrix a[1..n][1..n] is replaced by an upper Hessenberg matrix with identical eigenvalues. Rec￾ommended, but not required, is that this routine be preceded by balanc. On output, the Hessenberg matrix is in elements a[i][j] with i ≤ j+1. Elements with i > j+1 are to be thought of as zero, but are returned with random values. { int m,j,i; float y,x; for (m=2;m<n;m++) { m is called r + 1 in the text. x=0.0; i=m; for (j=m;j<=n;j++) { Find the pivot. if (fabs(a[j][m-1]) > fabs(x)) { x=a[j][m-1]; i=j;