Nash - Compact Numerical Methods for Computers (523163), страница 24
Текст из файла (страница 24)
Mixed modearithmetic is avoided by the use of the constant reltest. Thevariable msame is used in place of m to avoid confusion during thescope of the ‘while’ loop.}end; {loop on i}{STEP 15 -- now part of while loop control}if msame<nRow thenbegin {STEP 16 -- multiplication of vector by transformed B matrix}for i:=1 to nRow dobegins:=0.0;for j:=1 to nRow do s:=s+A[i,j+nRow]*Y[j];X[i]:=s;end; {loop on i}end; {if msame < nRow}end; {while loop -- STEP 17 is now part of the while loop control}if itcount>=itlimit then itcount:=-itcount; {set negative toindicate failure to converge within iteration limit}shift:=ev; {to pass eigenvalue to calling program}end; {alg10.pas == gii}107108Compact numerical methods for computersExample 9.1. Inverse iterationThe following output, from a Data General NOVA operating in 23-bit binaryarithmetic, shows the application of algorithm 10 to the algebraic eigenproblem ofthe order-4 Hilbert segment (see appendix 1).RUNENHGII OCT 21 76GAUSS ELIMINATION FOR INVERSE: ITERATIONORDER=? 4HILBERT SEGMENTSHIFT=? 0APPROX EV= 0APPROX EV= 9.67397E-5APPROX EV= 9.66973E-5APPROX EV= 9.66948E-5APPROX EV= 9.66948E-5CONVERGED TO EV= 9.66948E-5 IN 5 ITNS4 EQUAL CPNTSHILBERT SEGMENTVECTOR-2.91938E-2.328714-.791412.514551RESIDUALS-2.98023E-8-7.45058E-8-2.98023E-8-2.98023E-89.3.
SOME NOTES ON THE BEHAVIOUR OF INVERSE ITERATIONThe algorithm just presented may in some details appear complicated.(i) The convergence test uses a comparison of all elements in the vector. Formany applications the norm of y or some similar measure may suffice; however, itis not foolproof, particularly when the starting vector is set to some simple choicesuch as a column of ones. In this case, inverse iteration with a diagonal matrixAii = i, a unit matrix B and a shift of zero will ‘converge’ at iteration 2, which is itsearliest opportunity to stop.
However, the vector is left very much in error,though the dominant component has converged.(ii) The eigenvalue is given byshift + xi/ yi = k + xi/ yi(9.20)from the analysis of equations (9.16)-(9.18). When the element yi is zero, ofcourse, this is not suitable for determining the eigenvalue. Therefore the programmust save the vectors x and y, search for the largest element in y, then divide itinto the corresponding element of x in order to get the eigenvalue.
It is temptingto suggest that the expression should be simplyshift + 1/yi = k + 1/y i(9.21)since a normalisation is performed at each stage. Alas, too many matrices haveThe algebraic eigenvalue problem109symmetries which cause the dominant ‘component’ of a vector to be a pair ofelements with equal magnitude but opposite sign. A program may thereforechoose first one, then the other, as dominant component, thereby calculating acorrect eigenvector but an eigenvalue having the wrong sign!(iii) The limit on the number of iterations is a necessary precaution since theconvergence test may be too stringent.
Unfortunately, the rate of convergencemay be very slow since it is governed by the eigenvalue distribution. If theprogram can be controlled manually, one may prefer to allow the process tocontinue until the approximate eigenvalue has settled on a value, then interveneto halt execution if the convergence test does not become satisfied after a fewmore iterations. In fully automatic computation, some limit, say 100 iterations,does not seem unreasonable.A comparison of Gauss elimination and Givens’ reduction algorithms forinverse iteration using the nine test matrices of appendix 1 suggests that theformer is slightly preferable. Using matrices of order 5, 10, 15 and 20 and a shiftof k = 0, I timed the execution of inverse iteration including the triangularisation.Of the 36 problems, two were theoretically impossible since the shift was exactlybetween two eigenvalues.
Both methods failed to converge in 101 iterations forthese two problems (the Wilkinson W – matrices of even order) and the vectors atthe point of termination had large residuals so these two cases were dropped fromthe test. On the remaining 34, Gauss elimination yielded a smaller residual thanGivens’ reduction on 20, had a lower execution time on 21 for a total time (a sumof the separate times) of 3489·3 seconds compared to 4160·9 seconds, and failedto satisfy the convergence requirements in 101 iterations a total of 8 times, ascompared to 11 times out of 34 for Givens’ reduction. The above tests wereperformed in BASIC on a Data General NOVA operating interpretively in sixhexadecimal (base 16) digit arithmetic.Together with the overall timing, we should consider the fact that the methodshave very different relative performances depending on the matrix in question; inparticular, how many interchanges must be performed in Gauss elimination.Furthermore, as we might expect, the Givens’ method is faster for sparsematrices.Note that algorithm 10 is not the most compact form for the ordinary algebraiceigenvalue problem (2.62), since the Gauss elimination algorithm 5 givesPA' x = LRx = Pex(9.22)by virtue of (6.27), where P is the permutation matrix defined by the interchanges resulting from pivoting.
These can be stored, as discussed in §6.3 in asingle integer vector of indices q. Then to perform inverse iteration, it is onlynecessary to store this vector plus a working array n by n instead of the n by 2narray used in algorithm 10. Two vectors x and y are still needed. The elements ofthe lower-triangular matrix L are1for i = j(9.23)Lij = 0for j > ifor j < i.mij110Compact numerical methods for computersThe subdiagonal elements mij are left in place after the Gauss elimination step.and the upper-triangular matrix R forms the upper triangle of the working array.Then the inverse iteration step (9.10a) involves the forward-substitutionL v = Px i(9.23)Ryi = v.(9.25)and back-substitutionThe latter substitution step has been treated before.
but the former is simplified bythe ones on the diagonal of L so that(9.26)(9.27)The calculation can be arranged so that u is not needed. that is so x and y are theonly working vectors needed.9.4. EIGENSOLUTIONS OF NON-SYMETRIC ANDCOMPLEX MATRICESThe algebraic eigenvalue problem is considerably more difficult to solve when thematrix A is non-cymmetric or complex. If the matrix is Hermitian. that is.
if thecomplex-conjugate transpose of A is equal to A, the eigenvalues are then real andmethods for real symmetric matrices can easily be extended to solve this case(Nash 1974). However, in general, it is possible for the matrix to be defective,that is, have less than n eigenvectors. and the problem may be ill conditioned inthat the eigenvalues may be highly sensitive to small changes in the matrixelements.
The procedures which have been published for this problem are, by andlarge, long. They must, after all, contend with the possibility of complex eigenvalues and eigenvectors. On a small machine, the space which must be allocatedfor these very rapidly exhausts the memory available, since where previously onespace was required, two must now be reserved and the corresponding programcode is more than doubled in length to handle the arithmetic.Fortunately, such matrices occur rarely in practice.
For both real and complexcases I have used a direct translation of Eberlein’s ALGOL program comeig(Wilkinson and Reinsch 1971). This uses a generalisation of the Jacobi algorithmand I have found it to function well. One point which is not treated in Eberlein’sdiscussion is that the eigenvectors can all be multiplied by any complex number, c,and still be eigenvectors. In the real symmetric case all arithmetic is in the realdomain and by normalisation of the eigenvectors the results of a computation canbe standardised with the exception of eigenvectorc corresponding to multipleeigenvalues. In the complex case. however, the eigenvectorx + iy(9.28)The algebraic eigenvalue problemwill be completely unrecognisable after multiplication byc = re iφ = r cos φ + ir sin φthat is, we obtainx' + iy' = (xr cos φ – yr sin φ) + i(xr sin φ + yr cos φ).111(9.29)(9.30)Therefore, it is useful to standardise every computed eigenvector so that thelargest.
component is(9.31)1 + i0.Furthermore it is worthwhile to compute residuals for each vector. While this maybe a trivial programming tack for those familiar with complex numbers and linearalgebra, algorithms for both standardisation and residual calculation follow as anaid to those less familiar with the eigenproblem of general square matrices.We now present three algorithms which are intended to be used together:algorithm 11, to standardise a complex (eigen)vector;algorithm 12, to compute residuals for a purported complex eigenvector;andalgorithm 26, Eberlein’s Jacobi-like method for eigensolutions of a complexsquare matrix.The driver program DR26.PAS on the software diskette is designed to use these threealgorithms in computing and presenting eigensolutions of a general square matrix,that is, a square matrix which may have real or complex elements and need not haveany symmetries in its elements.Algorithm 11.
Standardisation of a complex vectorprocedure stdceigv(n: integer; {number of elements in vector}var T, U: x-matrix); {eigenvector k is given(column k of T) + sqrt(-1) * (column k of U)= T[.,k] + sqrt(-1) * U[.,k] }{algll.pas == Standardisation of complex eigensolutions.Copyright 1988 J.C.Nash}vari, k, m : integer;b, e, g, s : real;beginwriteln(‘algll.pas -- standardized eigensolutions’);for i := 1 to n do {loop over the eigensolutions}begin {the standardization of the solution so that the largestelement of each solution is set to 1 + 0 sqrt(-1)}g := T[1,i]*T[1,i]+U[1,i]*U[1,i]; {STEP 1}{the magnitude of element 1 of eigensolution i}k := 1; {to set index for the largest element so far}if n>1 thenbegin {STEP 2}for m := 2 to n do {loop over the other elements of the solution}112Compact numerical methods for computersAlgorithm 11.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
