Nash - Compact Numerical Methods for Computers (523163), страница 25
Текст из файла (страница 25)
Standardisation of a complex vector (cont.)begin {STEP 3}b := T[m,i]*T[m,i]+U[m,i]*U[m,i]; {the magnitude of element m}if b>g then {STEP 4}begin {STEP 5}k := m; {save the index of the largest element}g := b; {and its size}end; {if b>g}end; {loop on m -- STEP 6}end; {if n>1}e := T[k,i]/g; {STEP 7}s := -U[k,i]/g; {e & s establish the rotation constant in Eq.
9.29}fork:= 1 to n do {STEP 8}begin {the rotation of the elements}g := T[k,i]*e-U[k,i]*s; U[k,i] := U[k,i]*e+T[k,i]*s; T[k,i] := g;end; {loop on k}end, {loop on i -- over the eigensolutions}end; {alg11.pas == stdceigv}Algorithm 12. Residuals of a complex eigensolutionprocedure comres( i, n: integer; {eigensolution index for whichresiduals wanted, and order of problem}A, Z, T, U, Acopy, Zcopy : rmatrix);{outputfrom comeig (alg26). A and Z store theeigenvalues, T and U the eigenvectors, andAcopy and Zcopy provide a copy of theoriginal complex matrix.}(alg12pas == Residuals for complex eigenvalues and eigenvectors.This is slightly different in form from the step-and-description algorithmgiven in the first edition of Compact Numerical Methods; we work with thei’th eigensolution as produced by comeig.Copyright 1988 J.C.Nash}varj, k: integer;g, s, ss : real;beginwriteln(‘alg12.pas -- complex eigensolution residuals’);ss := 0.0; {sum of squares accumulator}for j := 1 to n do {STEP 1}begin {computation of the residuals, noting that theeigenvalue is located on the diagonal of A and Z}s := -A[i,i]*T[j,i]+Z[i,i]*U[j,i]; g := -Z[i,i]q[j,i]-A[i,i]*U[j,i];{s + sqrt(-1) g = -eigenvalue * vector_element_j}for k := 1 to n dobegins := s+Acopy[j,k]*T[k,i]-Zcopy[j,k]*U[k,i];g := g+Acopy[j,k]*U[k,i]+Zcopy[j,k]*T[k,i];end; {loop on k}writeln(‘(‘,s,’,‘,g,’)’);ss := ss+s*s+g*g;The algebraic eigenvalue problem113Algorithm 12.
Residuals of a complex eigensolution (cont.)end; {loop on j}writeln(‘Sum of squares = ’,ss);writeln;end; {alg12.pas == comres}Algorithm 26. Eigensolutions of a complex matrix by Eberlein’s methodprocedure comeig( n : integer; {order of problem}var itcount: integer; {on entry, a limit to the iterationcount; on output, the iteration count to convergencewhich is set negative if the limit is exceeded}var A, Z, T, U : rmatrix); {the matrix for which theeigensolutions are to be found is A + sqrt(-1)*Z,and this will be transformed so that the diagonalelements become the eigenvalues; on exit,T + sqrt(-1)*U has the eigenvectors in its columns}{alg26.pas == Pascal version of Eberlein’s comeig.Translated from comeig.bas, comeig.for and the original Algolversion in Eberlein (1971).Copyright J C Nash 1988}varRvec : rvector;{Rvec was called en by Eberlein, but we use upper case vector names.}i, itlimit, j, k, kl, m, m9, n1 : integer;aki, ami, bv, br, bi : real;c, cli, clr, c2i, c2r, ca, cb, ch, cos2a, cot2x, cotx, cx : real;d, de, di, diag, dr, e, ei, er, eps, eta, g, hi, hj, hr : real;isw, max, nc, nd, root1, root2, root : real;s, s1i, s1r, s2i, s2r, sa, sb, sh, sig, sin2a, sx : real;tanh, tau, te, tee, tern, tep, tse, zki, zmi : real;mark : boolean;begin {comeig}{Commentary in this version has been limited to notes on differencesbetween this implementation and the Algol original.}writeln(‘alg26.pas -- comeig’);eps := Calceps; {calculate machine precision}{NOTE: we have not scaled eps here, but probably should do so to avoidunnecessary computations.}mark := false; n1 := n-1;for i := l to n dobeginfor j := 1 to n dobegin {initialize eigenvector arrays}T[i,j] := 0.0; U[i,j] := 0.0; if i=j then T[i,i] := 1.0;end; {loop on j}end, {loop on i}itlimit := itcount; {use value on entry as a limit}itcount := 0; {and then set counter to zero}while (itcount<=itlimit) and (not mark) doI114Compact numerical methods for computersAlgorithm 26.
Eigensolutions of a complex matrix by Eberlein’s method (cont.)beginitcount := itcount+1; {safety loop counter}tau := 0.0; {convergence criteria}diag := 0.0; {to accumulate diagonal norm}for k := l to n dobegintem := 0;for i := 1 to n do if i<>k then tern := tem+ABS(A[i,k])+ABS(Z[i,k]);tau := tau+tem; tep := abs(A{k,k])+abs(Z[k,k]);diag := diag+tep; {rem accumulate diagonal norm}Rvec[k] := tem+tep;end; {loop on k}writeln(‘TAU=’,tau,‘ AT ITN ’,itcount);for k := 1 to n1 do {interchange rows and columns}beginmax := Rvec[k]; i := k; k1 := k+1;for j := k1 to n dobeginif max<Rvec[j] thenbeginmax := Rvec[j]; i := j;end; {if max<Rvec[j]}end; {loop on j}if i<>k thenbeginRvec[i] := Rvec[k];for j := 1 to n dobegintep := A[k,j]; A[k,j] := A[i,j]; A[i,j] := tep; tep := Z[k,j];Z[k,j] := Z[i,j]; Z[i,j] := tep;end; {loop on j}for j := l to n dobegintep := A[j,k]; A[j,k] := A[j,i]; A[j,i] := tep; tep := Z[j,k];Z[j,k] := Z[j,i]; Z[j,i] := tep; tep := T[j,k]; T[j,k] := T[j,i];T[j,i] := tep; tep := U[j,k]; U[j,k] := U[j,i]; U[j,i] := tep;end; {loop on j}end; {if i<>k}end; {loop on k}if tau>=l00.0*eps then {note possible change in convergence test fromform of Eberlein to one which uses size of diagonal elements}begin {sweep}mark := true;for k := 1 to n1 do {main outer loop}begink1 := k+1;for m := k1 to n do {main inner loop}beginhj := 0.0; hr := 0.0; hi := 0.0; g := 0.0;for i := l to n dobeginif (i<>k) and (i<>m) thenThe algebraic eigenvalue problem115Algorithm 26.
Eigensolutions of a complex matrix by Eberlein’s method (cont.)beginhr := hr+A[k,i]*A[m,i]+Z[k,i]*Z[m,i];hr := hr-A[i,k]*A[i,m]-Z[i,k]*Z[i,m];hi := hi+Z[k,i]*A[m,i]-A[k,i]*Z[m,i];hi := hi-A[i,k]*Z[i,m]+Z[i,k]*A[i,m];te := A[i,k]*A[i,k]+Z[i,k]*Z[i,k]+A[m,i]*A[m,i]+Z[m,i]*Z[m,i];tee := A[i,m]*A[i,m]+Z[i,m]*Z[i,m]+A[k,i]*A[k,i]+Z[k,i]*Z[k,i];g := g+te+tee; hj := hj-te+tee;end; {if i<>k and i<>m}end; {loop on i}br := A[k,m]+A[m,k]; bi := Z[k,m]+Z[m,k]; er := A[k,m]-A[m,k];ei := Z[k,m]-Z[m,k]; dr := A[k,k]-A[m,m]; di := Z[k,k]-Z[m,m];te := br*br+ei*ei+dr*dr; tee := bi*bi+er*er+di*di;if te>=tee thenbeginisw := 1.0; c := br; s := ei; d := dr, de := di;root2 := sqrt(te);endelse {te<tee}beginisw := -1.0; c := bi; s := -er; d := di; de := dr,root2 := sqrt(tee);end;root1 := sqrt(s*s+c*c); sig := -1.0; if d>=0.0 then sig := 1.0;sa := 0.0; ca := -1.0; if c>=0.0 then ca := 1.0;if root1<=eps thenbeginsx := 0.0; sa := 0.0; cx := 1.0; ca := 1.0;if isw<=0.0 thenbegine := ei; bv := -br;endelsebegine := er; bv := bi;end; {if isw<=0.0}nd := d*d+de*de;endelse {root1>eps}beginif abs(s)>eps thenbeginca := c/root1; sa := s/root1;end; {abs(s)>eps}cot2x := d/rootl; cotx := cot2x+(sig*sqrt(1.0+cot2x*cot2x));sx := sig/sqrt(1.0+cotx*cotx); cx := sx*cotx;{find rotated elements}eta := (er*br+ei*bi)/root1; tse := (br*bi-er*ei)/root1;te := sig*(tse*d-de*root1)/root2; tee := (d*de+root1*tse)/root2;nd := root2*root2+tee*tee; tee := hj*cx*sx; cos2a := ca*ca-sa*sa;sin2a := 2.0*ca*sa; tem := hr*cos2a+hi*sin2a;tep := hi*cos2a-hr*sin2a; hr := hr*cx*cx-tem*sx*sx-ca*tee;116Compact numerical methods for computersAlgorithm 26.
Eigensolutions of a complex matrix by Eberlein’s method (cont.)hi := hi*cx*cx+tep*sx*sx-sa*tee;bv := isw*te*ca+eta*sa; e := ca*eta-isw*te*sa;end; {root1>eps}{label ‘enter1’ is here in Algol version}s := hr-sig*root2*e; c := hi-sig*root2*bv; root := sqrt(c*c+s*s);if root<eps thenbegincb := 1.0; ch := 1.0; sb := 0.0; sh := 0.0;end {if root<eps}else {root>=eps}begincb := -c/root; sb := s/root; tee := cb*bv-e*sb; nc := tee*tee;tanh :=root/(g+2.0*(nc+nd)); ch := 1.0/sqrt(1.0-tanh*tanh);sh := ch*tanh;end; {root>=eps}tem := sx*sh*(sa*cb-sb*ca); c1r := cx*ch-tem; c2r := cx*ch+tem;c1i := -sx*sh*(ca*cb+sa*sb); c2i := c1i; tep := sx*ch*ca;tem := cx*sh*sb* s1r := tep-tem; s2r := -tep-tem; tep := sx*ch*sa;tem := cx*sh*cb: s1i := tep+tem; s2i := tep-tem;tem := sqrt(s1r*s1r+s1i*s1i); tep := sqrt(s2r*s2r+s2i*s2i);if tep>eps then mark := false;if (tep>eps) and (tem>eps) thenbeginfor i := 1 to n dobeginaki := A[k,i]; ami := A[m,i]; zki := Z[k,i]; zmi := Z[m,i];A[k,i] := c1r*aki-c1i*zki+s1r*ami-s1i*zmi;Z[k,i] := c1r*zki+c1i*aki+s1r*zmi+s1i*ami;A[m,i] := s2r*aki-s2i*zki+c2r*ami-c2i*zmi;Z[m,i] := s2r*zki+s2i*aki+c2r*zmi+c2i*ami;end; {loop on i}for i := l to n dobeginaki := A[i,k]; ami := A[i,m]; zki := Z[i,k]; zmi := Z[i,m];A[i,k] := c2r*aki-c2i*zki-s2r*ami+s2i*zmi;Z[i,k] := c2r*zki+c2i*aki-s2r*zmi-s2i*ami;A[i,m] := -s1r*aki+s1i*zki+c1r*ami-c1i*zmi;Z[i,m] := -s1r*zki-s1i*aki+c1r*zmi+c1i*ami;aki := T[i,k]; ami := T[i,m]; zki := U[i,k]; zmi := U[i,m];T[i,k] := c2r*aki-c2i*zki-s2r*ami+s2i*zmi;U[i,k] := c2r*zki+c2i*aki-s2r*zmi-s2i*ami;T[i,m] := -s1r*aki+s1i*zki+c1r*ami-c1i*zmi;U[i,m] := -s1r*zki-s1i*aki+c1r*zmi+c1i*ami;end; {loop on i}end; {if tep and tern>eps}end; {loop on m}end; {loop on k}end {if tau>=l00*eps}else mark := true; {to indicate convergence}end; {while itcount<=itlimit}if itcount>itlimit then itcount := -itcount; {negative iteration countmeans process has not converged properly}end; {alg26.pas == comeig}The algebraic eigenvalue problem117Example 9.2.
Eigensolutions of a complex matrixThe following output from a Data General NOVA operating in 23-bit binaryarithmetic shows the computation of the eigensolutions of a complex matrix dueto Eberlein from the test set published by Gregory and Karney (1969). Thenotation ( , ) is used to indicate a complex number, real part followed byimaginary part. Note that the residuals computed are quite large by comparisonwith those for a real symmetric matrix. This is due to the increased difficulty ofthe problem, to the extra operations needed to take account of the complexnumbers and to the standardisation of the eigenvectors, which will introduce someadditional errors (for instance, in the first eigenvector, 5·96046E–8 for zero).This comment must be tempered by the observation that the norm of the matrix isquite large, so that the residuals divided by this norm are still only a reasonablysmall multiple of the machine precision.RUNENHCMG - COMEIGORDER? 3ELEMENT ( 1 , 1ELEMENT ( 1 , 2ELEMENT ( 1 , 3ELEMENT ( 2 , 1ELEMENT ( 2 , 2ELEMENT ( 2 , 3ELEMENT ( 3 , 1ELEMENT ( 3 , 2ELEMENT ( 3 , 3AT SEPT 3 74);REAL=?);REAL=?);REAL=?);REAL=?);REAL=?);REAL=?);REAL=?);REAL=?);REAL=?1 IMAGINARY? 23 IMAGINARY? 421 IMAGINARY? 2243 IMAGINARY? 4413 IMAGINARY? 1415 IMAGINARY? 165 IMAGINARY? 67 IMAGINARY? 825 IMAGINARY? 26TAU= 194 AT ITN 1TAU= 99,7552 AT ITN 2TAU= 64,3109 AT ITN 3TAU= 25,0133 AT ITN 4TAU= 7,45953 AT ITN 5TAU= .507665 AT ITN 6TAU= 6.23797E-4 AT ITN 7TAU= 1.05392E-7 AT ITN 8EIGENSOLUTIONSRAW VECTOR 1( .371175 ,-.114606 )( .873341 ,-.29618 )( .541304 ,-.178142 )EIGENVALUE 1 =( 39,7761 , 42,9951 )VECTOR( .42108 , 1.15757E-2 )( 1 , 5.96046E-8 )( .617916 , 5.57855E-3 )RESIDUALS( 2.2918E-4 , 2.34604E-4 )( 5.16415E-4 , 5.11169E-4 )( 3.70204E-4 , 3.77655E-4 )RAW VECTOR 2(-9.52917E-2 ,-.491205 )( 1.19177 , .98026 )(-.342159 ,-9.71221E-2 )118Compact numerical methods for computersEIGENVALUES 2 =( 6.7008 ,-7.87591 )VECTOR(-.249902 ,-.206613 )( 1 , 1.19209E-7 )(-.211227 , 9.22453E-2 )RESIDUALS(-3.8147E-5 , 3.8147E-6 )( 7.55787E-5 ,-7.48634E-5 )(-1.52588E-5 , 2.57492E-5 )RAW VECTOR 3( .408368 , .229301 )(-.547153 ,-1.39186 )(-4.06002E-2 , .347927 )EIGENVALUE 3 =(-7.47744 , 6.88024 )VECTOR(-.242592 , .198032 )( 1 , 0 )(-.206582 ,-.110379 )RESIDUALS( 5.24521E-6 ,-4.00543E-5 )(-7.9155E-5 , 7.82013E-5 )( 2.81334E-5 ,-1.04904E-5 )PreviousChapter 10REAL SYMMETRIC MATRICES10.1.
THE EIGENSOLUTIONS OF A REAL SYMMETRIC MATRIXThe most common practical algebraic eigenvalue problem is that of determiningall the eigensolutions of a real symmetric matrix. Fortunately, this problem hasthe most agreeable properties with respect to its solution (see, for instance,Wilkinson 1965).(i) All the eigenvalues of a real symmetric matrix are real.(ii) It is possible to find a complete set of n eigenvectors for an order-n realsymmetric matrix and these can be made mutually orthogonal. Usually they arenormalised so that the Euclidean norm (sum of squares of the elements) is unity.Thus, the total eigenproblem can be expressed(10.1)AX = XEwhere the matrix X has column j such thatAxi = ejxj(10.2)where ej is the j th eigenvalue.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
