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

Jacobi’s original algorithm of 1846 searchedthe whole upper triangle at each stage and set the largest off-diagonal element to zero.This is a reasonable strategy for hand calculation, but it is prohibitive on a computersince the search alone makes each Jacobi rotation a process of order N 2 instead of N .A better strategy for our purposes is the cyclic Jacobi method, where oneannihilates elements in strict order.

For example, one can simply proceed downthe rows: P12 , P13 , ..., P1n ; then P23 , P24 , etc. One can show that convergenceis generally quadratic for both the original or the cyclic Jacobi methods, fornondegenerate eigenvalues. One such set of n(n − 1)/2 Jacobi rotations is calleda sweep.The program below, based on the implementations in [1,2], uses two furtherrefinements:• In the first three sweeps, we carry out the pq rotation only if |apq | > for some threshold value=1 S05 n2(11.1.25)where S0 is the sum of the off-diagonal moduli,S0 =X|ars |(11.1.26)r<s• After four sweeps, if |apq | |app | and |apq | |aqq |, we set |apq | = 0and skip the rotation.

The criterion used in the comparison is |apq | <10−(D+2)|app |, where D is the number of significant decimal digits on themachine, and similarly for |aqq |.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).V0 = V · Pi11.1 Jacobi Transformations of a Symmetric Matrix467#include <math.h>#include "nrutil.h"#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\a[k][l]=h+s*(g-h*tau);void jacobi(float **a, int n, float d[], float **v, int *nrot)Computes all eigenvalues and eigenvectors of a real symmetric matrix a[1..n][1..n].

Onoutput, elements of a above the diagonal are destroyed. d[1..n] returns the eigenvalues of a.v[1..n][1..n] is a matrix whose columns contain, on output, the normalized eigenvectors ofa. nrot returns the number of Jacobi rotations that were required.{int j,iq,ip,i;float tresh,theta,tau,t,sm,s,h,g,c,*b,*z;b=vector(1,n);z=vector(1,n);for (ip=1;ip<=n;ip++) {Initialize to the identity matrix.for (iq=1;iq<=n;iq++) v[ip][iq]=0.0;v[ip][ip]=1.0;}for (ip=1;ip<=n;ip++) {Initialize b and d to the diagonalb[ip]=d[ip]=a[ip][ip];of a.z[ip]=0.0;This vector will accumulate termsof the form tapq as in equa}*nrot=0;tion (11.1.14).for (i=1;i<=50;i++) {sm=0.0;for (ip=1;ip<=n-1;ip++) {Sum off-diagonal elements.for (iq=ip+1;iq<=n;iq++)sm += fabs(a[ip][iq]);}if (sm == 0.0) {The normal return, which reliesfree_vector(z,1,n);on quadratic convergence tofree_vector(b,1,n);machine underflow.return;}if (i < 4)tresh=0.2*sm/(n*n);...on the first three sweeps.elsetresh=0.0;...thereafter.for (ip=1;ip<=n-1;ip++) {for (iq=ip+1;iq<=n;iq++) {g=100.0*fabs(a[ip][iq]);After four sweeps, skip the rotation if the off-diagonal element is small.if (i > 4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip])&& (float)(fabs(d[iq])+g) == (float)fabs(d[iq]))a[ip][iq]=0.0;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).In the following routine the n×n symmetric matrix a is stored as a[1..n][1..n]. On output, the superdiagonal elements of a are destroyed, but the diagonaland subdiagonal are unchanged and give full information on the original symmetricmatrix a. The vector d[1..n] returns the eigenvalues of a.

During the computation,it contains the current diagonal of a. The matrix v[1..n][1..n] outputs thenormalized eigenvector belonging to d[k] in its kth column. The parameter nrot isthe number of Jacobi rotations that were needed to achieve convergence.Typical matrices require 6 to 10 sweeps to achieve convergence, or 3n2 to25n Jacobi rotations. Each rotation requires of order 4n operations, each consistingof a multiply and an add, so the total labor is of order 12n3 to 20n3 operations.Calculation of the eigenvectors as well as the eigenvalues changes the operationcount from 4n to 6n per rotation, which is only a 50 percent overhead.468Chapter 11.Eigensystems}}for (ip=1;ip<=n;ip++) {b[ip] += z[ip];d[ip]=b[ip];z[ip]=0.0;}Update d with the sum of tapq ,and reinitialize z.}nrerror("Too many iterations in routine jacobi");}Note that the above routine assumes that underflows are set to zero.

Onmachines where this is not true, the program must be modified.The eigenvalues are not ordered on output. If sorting is desired, the followingroutine can be invoked to reorder the output of jacobi or of later routines in thischapter. (The method, straight insertion, is N 2 rather than N log N ; but since youhave just done an N 3 procedure to get the eigenvalues, you can afford yourselfthis little indulgence.)void eigsrt(float d[], float **v, int n)Given the eigenvalues d[1..n] and eigenvectors v[1..n][1..n] as output from jacobi(§11.1) or tqli (§11.3), this routine sorts the eigenvalues into descending order, and rearrangesthe columns of v correspondingly.

The method is straight insertion.{int k,j,i;float p;for (i=1;i<n;i++) {p=d[k=i];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).else if (fabs(a[ip][iq]) > tresh) {h=d[iq]-d[ip];if ((float)(fabs(h)+g) == (float)fabs(h))t=(a[ip][iq])/h;t = 1/(2θ)else {theta=0.5*h/(a[ip][iq]);Equation (11.1.10).t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));if (theta < 0.0) t = -t;}c=1.0/sqrt(1+t*t);s=t*c;tau=s/(1.0+c);h=t*a[ip][iq];z[ip] -= h;z[iq] += h;d[ip] -= h;d[iq] += h;a[ip][iq]=0.0;for (j=1;j<=ip-1;j++) {Case of rotations 1 ≤ j < p.ROTATE(a,j,ip,j,iq)}for (j=ip+1;j<=iq-1;j++) {Case of rotations p < j < q.ROTATE(a,ip,j,j,iq)}for (j=iq+1;j<=n;j++) {Case of rotations q < j ≤ n.ROTATE(a,ip,j,iq,j)}for (j=1;j<=n;j++) {ROTATE(v,j,ip,j,iq)}++(*nrot);}11.2 Reduction of a Symmetric Matrix to Tridiagonal Form469}}CITED REFERENCES AND FURTHER READING:Golub, G.H., and Van Loan, C.F.

1989, Matrix Computations, 2nd ed. (Baltimore: Johns HopkinsUniversity Press), §8.4.Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide, 2nd ed., vol. 6 ofLecture Notes in Computer Science (New York: Springer-Verlag). [1]Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag). [2]11.2 Reduction of a Symmetric Matrixto Tridiagonal Form: Givens andHouseholder ReductionsAs already mentioned, the optimum strategy for finding eigenvalues andeigenvectors is, first, to reduce the matrix to a simple form, only then beginning aniterative procedure.