c11-5 (Numerical Recipes in C), страница 2

Recommended, but not required, is that this routine be preceded by balanc. On output, theHessenberg matrix is in elements a[i][j] with i ≤ j+1. Elements with i > j+1 are to bethought 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;Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.

Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).required elements in a column of the matrix. Householder reduction to Hessenbergform is in fact an accepted technique. An alternative, however, is a procedureanalogous to Gaussian elimination with pivoting. We will use this eliminationprocedure 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 matricesfor which the Householder reduction, being orthogonal, is stable and elimination isnot, 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 therth stage, the original matrix A ≡ A1 has become Ar , which is upper Hessenbergin its first r − 1 rows and columns.

The rth stage then consists of the followingsequence of operations:• Find the element of maximum magnitude in the rth column below thediagonal. If it is zero, skip the next two “bullets” and the stage is done.Otherwise, suppose the maximum element was in row r0 .• Interchange rows r 0 and r + 1. This is the pivoting procedure. To makethe permutation a similarity transformation, also interchange columns r 0and r + 1.• For i = r + 2, r + 3, . .

. , N , compute the multiplier486Chapter 11.Eigensystems}}CITED REFERENCES AND FURTHER READING:Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag). [1]Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide, 2nd ed., vol.

6 ofLecture Notes in Computer Science (New York: Springer-Verlag). [2]Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§6.5.4. [3]11.6 The QR Algorithm for Real HessenbergMatricesRecall the following relations for the QR algorithm with shifts:Qs · (As − ks 1) = Rs(11.6.1)where Q is orthogonal and R is upper triangular, andAs+1 = Rs · QTs + ks 1= Qs · As · QTs(11.6.2)The QR transformation preserves the upper Hessenberg form of the original matrixA ≡ A1 , and the workload on such a matrix is O(n2 ) per iteration as opposedto O(n3 ) on a general matrix. As s → ∞, As converges to a form wherethe eigenvalues are either isolated on the diagonal or are eigenvalues of a 2 × 2submatrix on the diagonal.As we pointed out in §11.3, shifting is essential for rapid convergence. A keydifference here is that a nonsymmetric real matrix can have complex eigenvalues.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.

Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).}}if (i != m) {Interchange rows and columns.for (j=m-1;j<=n;j++) SWAP(a[i][j],a[m][j])for (j=1;j<=n;j++) SWAP(a[j][i],a[j][m])}if (x) {Carry out the elimination.for (i=m+1;i<=n;i++) {if ((y=a[i][m-1]) != 0.0) {y /= x;a[i][m-1]=y;for (j=m;j<=n;j++)a[i][j] -= y*a[m][j];for (j=1;j<=n;j++)a[j][m] += y*a[j][i];}}}.