Note on the ill-conditioned eigenvalue problem in elastic vibrations

COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 6 (1975) 65-77 © NORTH-HOLLAND PUBLISHING COMPANY

NOTE ON THE ILL-CONDITIONED EIGENVALUE PROBLEM IN ELASTIC VIBRATIONS G. FRIK and Th. LUNDE JOHNSEN lnstitut fiir Statik und Dynamik der Luft- und Raum[ahrtkonstruktionen, University of Stuttgart, Germany Received 23 December 1974 Two cases which could cause ill-conditioning of the lower modes and frequencies of an undamped elastic structure are discussed. The first case is when rigid body motions are connected to a relatively small moment of inertia. The second case is when the stiffness matrix is ill-conditioned with respect to the solution of linear equations. It is shown that the first kind of ill-conditioning is not serious provided that an appropriate method of computation is used. The second kind, which is more common and more serious, is dealt with in more detail. It is shown how the results may be improved in various ways by applying the Natural Factor approach of Argyris and Bronlund and computing singular values and vectors.

1. Introduction We consider the general eigenvalue problem Mv

= KVA

0.1 )

representing elastic vibrations of a structure idealized by the finite element method [1] to [4]. We shall further restrict ourselves to statically condensed problems (Guyan's reduction) [5]. In this case the matrices M and K are positive semidefinite and usually densely populated. The general eigenvalue problem (1.1) may now be solved by some fast direct-reduction method like the Householder-QR-iteration method [6]. In order to apply some such method, the problem must be reduced to special form. As we are in general only interested in a restricted number of the largest eigenvalues with corresponding eigenvectors (the smallest eigenvalues being physically meaningless), it is important to perform the reduction to special form in such a manner that the larger eigenvalues are preserved as accurately as possible. There are two properties that may lead to numerical difficulties. Firstly, the stiffness matrix K and possibly also M could be singular. This means that the structure is supported in a kinematically undetermined way, as will be the case for a free-free structure like a rocket or an aeroplane. The difficulty arises when a rigid body mode is connected to a relatively small moment of inertia. Although the limiting case (when this moment of inertia is zero) causes some eigenvalues to become unstable (see [14]), we find that the eigenvalues of (1.1) are rather insensitive even for very small moments of inertia. The second case occurs when the stiffness matrix (or the elastic part of the stiffness matrix after removal of rigid body modes) is nonsingular but ill-conditioned with respect to solution of linear equations. In this case the larger eigenvalues of (1.1) are senstive to small perturbations of K. Hence, as soon as K is explicitly formed, they are inaccurately determined, regardless of the

66

G. Frik and Th. Lunde Johnsen, Note on the ill-conditioned eigenvalue problem in elastic vibrations

reduction method. In this case the natural factor approach of Argyris and Bronlund [7,8] should be applied. Thus we obtain the larger eigenvalues more accurately. If a similar technique is applied to the mass matrix and the singular values are computed, we show that the smaller eigenvalues also become more accurate. This result also applies to the Simultaneous Vector iteration technique. In the following we assume that K and M are scaled such that the norms IIKII = IIMII = 1 .

( 1.2)

In particular we assume that K is reasonably equilibrated by a symmetric scaling. We shall not discuss various methods of scaling here, but we would like to mention that if rigid body modes are present, the scaling procedure is not as simple as in the nonsingular case [9, 10], and the particular form of a stiffness matrix should be utilized. Reduction Irethods exist where scaling is not important (see [3 D. Here, however, we shall provide error estimates in terms of norms, and furthermore represent orthogonal complements in terms of elementary unitary Hermitians. For these reasons we make the assumptions above.

2. Removal of rigid body modes

2.1. Methods In the static condensation procedure we may always make sure that M is nonsingular. Hence, we may Cholesky-decompose Minto (2.1)

and compute DM =B- t KB- l

.

(2.2)

The eigenvalues of (1.1) are the reciprocals of those of D M , and the zero eigenvalues of D M correspond to the rigid body modes. As the larger eigenvalues of D M may become very large due to the transformation with B-t and B- 1 , we may expect high relative errors in the smaller eigenvalues of D M due to the tri-diagonalization step of the eigenvalue procedure. Hence the method (2.1), (2.2) must be discarded. Note, however, that the scaling of K is in this case of no importance. In [11-13] the following method is applied: Let X, (n, m), be some set of orthonormal vectors satisfying KX=O,

xtX=I.

(2.3)

Then we compute Y=MX

and an orthogonal matrix

(2.4)

G. Frik and Th. Lunde Johnsen, Note on the ill-conditioned eigenvalue problem in elastic vibrations

67

m

(2.5) such that Qty=

[~J}m o

,

(2.6)

}n-m

i.e. Householder's method. Partitioning (2.7) L..,-'

L..,-'

m n-m

we have

stE y=O '

(n-m,m).

(2.8)

(n, n) ,

(2.9)

If the matrix

is nonsingular, then we may congruent-transform 0.1) to get MEEVE

= KEEVEA ,

(2.10)

where

MEE =S~MSE' KEE

=S~KSE

(2.11)

.

Hence, with the Cholesky decomposition (2.12) we obtain the special eigenvalue problem (2.13)

(2.14) Another, standard technique [14, 15], is to congruent-transform (1.1) with the matrix

Q= [X(9,

X] ,

L..,-'

L..,-'

n-m

m

(2.15)

68

G. Frik and Th. Lunde Johnsen, Note on the ill-conditioned eigenvalue problem in elastic vibrations

where X"', (n, n-m), is the orthogonal complement of X. Partitioning,

[ '"

'"

'"

'"

M=QtMQ=

M ll

(2.16)

"'t M 12

and

l

o

}n-m

oj }m

(2.17)

Letting

(2.18) and making Cholesky factorization '"

K ll

_ -

t

Ull U ll

'

(2.19)

we obtain '" U- 1 - U-tMD 11II II ll'

(2.20)

D 11 having the same eigenvalues as (1.1).

2.2. Round off error analysis It is obvious from (2.18) that the nonsingularity of the rigid body mass matrix M22 is of importance for the projection method (2.16)-(2.20) (see also [14 D. In the congruent transformation method (2.4)-(2.14) th~singularity of M22 implies also that S (of (2.9)) is singular; hence both methods break down if M~2 is singular. It is easy to prove (see [11-13]) the relation (2.21) between condition numbers. It turns out, however, that only the rigid body mass is important. As M 22 is positive definite, we have, using 0.2),

(2.22) From (2.20) we obtain, after some manipulations,

(2.23) and

(2.24)

G. Frik and Th. Lunde Johnsen, Note on the ill-conditioned eigenvalue problem in elastic vibrations

69

Let us now consider a perturbed eigenvalue problem (M

+ oM)u = (K + OK)u~

.

(2.25)

Let V be the K-orthogonal set of eigenvectors of (1.1) and

v

= [V, X] ,

(2.26)

(n, n) .

Then a congruent transformation of (2.25) with

Vgives (2.27)

with perturbations (2.28) Partitioning F M' F K' A and ii as follows:

FK

= rFKll

FK12 J }n-m

(2.29)

,

Fi /2 F K22 }m

and assuming II FK II

~

1, we may put FK /2 = 0 without lack of generality and obtain

Al - FMl1 - FM/2(A2-FM22

-FK X)-I F1/2)ii l = iii X.

(2.30)

The error in Xmay now be bounded by the Wieland-Hoffman theorem or the Gerschgorin circle theorem. It follows from (2.30) that the error term influenced by small rigid body masses is quadratic in 110M II II V1I 2 • If V, therefore, is well-conditioned, so is K 1l' In this case (2.30) shows that small rigid body masses have little influence on the accuracy of the eigenvalues when applying the projection method (2.16) - (2.20). Using a backward round off error analysis, we may show that the computed matrices /(11 and Mil of(2.18) and (2.19) may be obtained by exact computation on slightly perturbed matrices [16, p.152]. Furthermore the error in the eigenvalue based on the computed D ll of (2.20) are of the order IIKilIl2-t, where t is the number of binary digits in the mantissa of a single precision computer word. For the congruent transformation method (2.4)-(2.14) a similar statement is true, but here the error in the eigenvalues based on the computed DEE are of the order IIKE~112-t. Comparing with (2.24), we see that this bound is not as good as the one for the projection method. As pointed out by Rosanoff [22], the rigid body modes as computed from the geometry are more accurately determined than would be the case if they were computed from the stiffness matrix. It follows that although this has little influence on the accuracy of the eigenvalues, the

70

G. Frik and Th. Lunde Johnsen, Note on the ill-conditioned eigenvalue problem in elastic vibrations

lower eigenmodes computed become remarkably more accurate than if they were computed solely from the computed matrices M and K.

3. The ill-conditioned stiffness matrix Many structures result in an ill-conditioned stiffness matrix. Ill-conditioning may be due to highly varying material properties as well as to the geometry - e.g. a cantilever. In both cases lower mode information is poorly represented in the stiffness matrix. The Natural Factor method [7, 8] (see also [19]) may be used to obtain accurate eigenvalues in such cases [17, 18]. Let A be the natural factor of K satisfying (3.1)

Then A is formed directly from the physical data instead of from K. The natural factor is decomposed as (3.2)

A=QU,

where Q is orthogonal and' U upper triangular [20, p. 233]. The dynamic matrix is given as D= U-tMU.

(3.3)

The Natural Factor method is clearly compatible with static condensation and the methods of removing rigid body modes as described in section 2. The improvement of accuracy, however, is lost if the matrix D M of (2.2) is computed; hence that should be avoided. As the computation (3.3) is subjected to round off errors, we obtain, using the techniques described in [2l,p.30l] and [6,p.135], IID e -D c li ~ 0 (IIU-II1I1U-Il1l1u-tMIl 2- t )

II D bll

"'"

IIDbll

'

(3.4)

where De is the exact and Dc the computed dynamic matrix. In general the factor (3.5) will lie between 1 and II U- 1 11, but as the lower modes will frequently be connected to large masses, we have normally F(U, M) = 0(1).

(3.6)

As we also have F(U, V) <; 0(1< (M)) ,

(3.7)

(3.6) is always true if M is well conditioned. From the standard stiffness matrix - Cholesky approach

G. Frik and Th. Lunde Johnsen, Note on the ill-conditioned eigenvalue problem in elastic vibrations

71

we have the realistic bound

liD

-DII

e",

e

~

O(llU- 1 11 2 2- t ) ,

(3.8)

IIDeII which is clearly inferior to the bound (3.4). If the natural factor of M is also computed, M =BtB,

(3.9)

such that E =BU- 1

(3.10)

satisfies (3.11 ) we need form only E and compute its singular values a l' a 2 , exact E e and the computed E e we obtain

liE -E II e e ~ 0(1/ (K)1/2 2-tI2) liEe II "'" no

••• ,

an' For the difference between the

(3.12)

•

The eigenvalue or singular value methods based on the QR iteration will normally produce errors of the order (3.13) or (3.14) where Aie , aie denote the computed values as compared to the exact corresponding values AiD' aiE of the computed D, E. The bounds (3.13) and (3.14) should be compared with)(3.4) and (3.12) respectively. As the bounds (3.4) and (3.12) overestimate the error for the smaller values, (3.13) or (3.14) will dominate for those. Since = Ai' the bound (3.14) will give lower relative errors for the smaller eigenvalues, indicating that the method of singular values is more accurate. Some product form like (3.10) is used in the Simultaneous Vector Iteration; hence this is comparable with the errors in the singular value decomposition, although in general it is more !It:;curate since the product form is kept during the iteration procedure.

a;

4. Examples The method described have been tested on several examples to check the validity of the error analysis. We describe the results of the following tests.

72

G. Frik and Th. Lunde Johnsen, Note on the ill-conditioned eigenvalue problem in elastic vibrations

Example ,. Removing rigid body modes ofan example with well-conditioned eigenvectors using the two lttethods of section 2 - Congruent Transformation and Projection. In this case the data were generated as follows. Given two diagonal matrices A=

]}m

AI [

(4.1)

A 2 }n-m '

we determine a matrix X, (n, n), using a random number generator. Defining the matrices K= X-tjX- I

we obtain the

,

M =X-tAX- I

(4.2)

,

eigenprobl~'TI

(4.3)

MX=KXA.

If we partition X into X

= [XI '-,o-J

m

X2]

(4.4)

,

'-,o-J

n-m

the problem (4.3) splits up into (4.5)

The problem (4.3) is now solved by the two methods of section 2 and the results are compared. By letting Al of (4.1) contain one or more very small entries, we obtain a correspondingly illconditioned matrix M22 (see (2.16)). Table 1 Results for example I Assumed eigenvalues

1 3.552713678800 5E-14} 23.5527136788005E-14 A 3 3.000000000000 OE+OO 1 4 4.000000000000 OE+OO 5 5.000000000000 OE+O 6 6.000000000000 OE+OO 7 7.0000000000000E+00 8 8.0000000000000E+00 A 2 9 9.0000000000000E+00 10 \.000000000000 OE+OI 11 1.100000000000 OE+OI

Eigenvalues after removal of 4 rigid body modes, computed by a) Congrnenttransformanon

b) Projection method

1 \.099999999997 3E+OI 2 9.999999999737 2E+00 38.9999999999138E+00 47.9999999999516E+00 5 6.999999999889 4E+00 65.9999999999955E+00 7 4.999999999996 2E+00

1 2 3 4 5 6 7

\.099999999999 7E+OI \.0000000000004E+OI 8.999999999998 OE+OO 8.000000000000 OE+OO 7.000000000002 5E+00 5.999 999 999 999 OE +00 4.999999999998 8E+00

Table 1 displays results obtained on a CDC 6600 computer with a mantissa of about 14.4 decimal digits.

G. Frik and Th. Lunde Johnsen, Note on the ill-conditioned eigenvalue problem in elastic vibrations

73

Interpretation of the results. Although some of the rigid body masses Ai are only slightly above the working accuracy of the computer, we obtain very accurate eigenvalues for both methods. The better accuracy observed for the projection method seems to be generally true for all examples run so far, although frequently less outstanding than in this example. If A 1 (see (4.1» contains some negative entries, then the matrix M22 of (2.16) will not be positive definite and the computation (2.19) breaks down. In this case the following approach works well: 1. Compute the eigenvalues A2 and the eigenvectors ii of M22 • 2. Compute the diagonal matrices At) and At) such that the diagonal elements ofA'~+) are the positive square roots of the absolute values of the diagonal elements of A2 and also :i'(-)~(+) 2

u2

= ~A 2

(4.6)

.

3. Compute B(+) = 12

M12 H(~(+»)-l 2 ,

B(-) 12

=M12 H(~(-»)-l 2

(4.7)

and finally

M12 M-221 M12t =B(+)(B(-»)t 12 12

(4.8)

(see (2.18». Note that the congruent transformation method (2.6) -(2.Il) also works well in this case. In addition, note that if we try to compute B 12 =M12 M221 and M12 M221 M: 2 =M12 B~2 instead of (4.8), then the errors due to the inversion of M22 will produce meaningless results when applied to examples as ill-conditioned as example 1. This has been confirmed in several cases. Example 2. Vibrating frame structure with ill-conditioned stiffness matrix.

y

Fig. 1. Vibrating frame structure.

The ill-conditioned stiffness matrix is achieved by letting some of the bars be very flexible. The mass matrix is well-conditioned. Eigenvalues are computed with the classical stiffness matrix by using a) Cholesky approach, b) Natural Factor approach with dynamic matrix, and c) Natural Factor approach with singular value computation. Table 2 displays results obtained on a UNIYAC 1108 computer with a mantissa of about 8 decimal digits.

74

G. Frik and Th. Lunde Johnsen, Note on the ill-conditioned eigenvalue problem in elastic vibrations Table 2 Results for example 2 Exact eigenvalues

Cholesky approach

1 2 3 4 5 6 7 8 9 10 11 12 13 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14

1.299695216590346£+11 5.701685651328843£+09 6.027254702976073£+08 6.710478809179807£+07 3.555888713938353£+07 7~54719369610128£+06

1.024559321256470£+06 9.512 988 603 804167£+04 2.105771 021631336£+00 1.718926687755162£-01 1.419780091414438£-02 5.307496994420958£-03 6.206774 886466903£-04 4.797411 659080374£-04

8.892 010 9£+09 9.7451468£+08 1.0844303£+08 3.7108927£+07 1.2127104£+07 1.368001 7£+06 1.7408519£+05 9.0102889£+04 3.5196365£+02 1.7598183£+02 8.799092 2£+01 4.3995463£+01 2.1997733£+01 0.0000000

Natural Factor approach with dynamic matrix

Natural Factor approach with singular value computation

I 2 3 4 5 6 7 8 9 10 11 12 13 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14

1.2994823£+11 5.6997852£+09 6.0231248£+08 6.7100274£+07 3.5556921£+07 7.647 863 6£+06 1.0232006£+06 8.7931302£+04 1.0991413£+04 5.4957070£+03 2.7478536£+03 1.3739269£+03 6.8696349£+02 0.0000000

1.2994621£+11 5.6997898£+09 6.0231450£+08 6.7104476£+07 3.5558658£+07 7.6552843£+06 1.0245620£+06 9.5129987£+04 2.1057554£+00 1.4191075£-02 1.7177740£-01 6.2133084£-04 5.302 156 1£-03 4.7976021£-04

Interpretation of the results. By inspecting the exact eigenvalues, keeping in mind that the mass matrix is well-conditioned, we observe that the stiffness matrix in this case is so ill-conditioned that the traditional approach of constructing the stiffness matrix explicitly and computing its Cholesky factor in order to compute the dynamic matrix is doomed to failure. The larger eigenvalues are therefore not even correct in order of magnitude, as correctly foreseen by (3.8). For the smaller eigenvalues, on the other hand, some useful information is still present in the dynamic matrix, but due to errors in the QR iteration method this is lost for the smallest eigenvalues as indicated by (3.13). In the case of the remaining eigenvalues, i.e. those in the centre of the eigenvalue spectrum, some useful information did indeed pass through the two processes, as indicated by the fact that eigenvalue no. 8 contains one good digit. Looking now at the Natural Factor approach with explicit computation of the dynamic matrix, we note that the larger eigenvalues are fairly accurate, as indicated by the error bound (3.4). Note that (3.6) will hold because we have a very well-conditioned mass matrix (see (3.7)). The inaccuracy of the smaller eigenvalues is due only to errors in the QR iteration process (see (3.13)). Finally, we note that the Natural

G. Frik and Th. Lunde Johnsen, Note on the ill-conditioned eigenvalue problem in elastic vibrations

75

Factor approach with singular value decomposition is similar to the previous approach and, in addition, evaluates the smaller eigenvalues to a greater accuracy. This is because the difference in the magnitude between the eigenvalues is now represented in terms of the singular values which are the square roots of the eigenvalues so that the inaccuracy is much less serious (see (3.14))

5. Conclusions The rigid body transformation which separates the elastic modes and the rigid body modes of a free vibrating structure is usually very stable. It breaks down when zero masses are connected to rigid body modes, but in this case the system is physically meaningless. Even for extremely small masses connected to the rigid body modes the rigid transformations of section 2 are very stable, at least when the matrix of eigenvectors of the problem is well-conditioned. If the stiffness matrix is ill-conditioned, the method utilizing the natural factor produces remarkably better eigenvalues, and the cost of the extra computation required is minimal. If in addition the mass matrix is also factorized and the singular values are computed directly, the smaller eigenvalues will also be more accurately computed. One should, however, take into account that the singular value computation costs about 30 - 40% more computation time than the eigenvalue computation. We would like to conclude these considerations with the following Remark. If the eigenvalues and eigenvectors are computed by Simultaneous Vector Iteration using the Natural Factor approach, their accuracy will be at least as good as the singular value approach. This is due to the fact that the dynamic matrix is never formed and that the methods of computing the singular values and the eigenvalues are numerically equivalent. The following application of the natural factors is also worth noting: Let X, (n, k), be a set of computed eigenvectors of (1.1) corresponding to the k largest eigenvalues, and let us assume that the computation of X is based upon K and M as explicitly formed in the computer. Now, if K is very illconditioned, as is frequently the case, the algorithm will not be able to decouple the first few lower modes. This may be seen by looking at the residual Kx -M

1

I

x

=11.

(5.1 )

If Ais large, we have K x ~ 0, and hence only the invariant subspace corresponding to all large Ais well determined. Hence, if the (k+ 1)-th eigenvalue Ak+l of (1.1) satisfiesiAk+1 ~ AI' we would expect to obtain a very accurate representation of the invariant subspace corresponding to the first k eigenvalues. Now suppose that A and B are the natural factors of K and M respectively. Then, if we compute YK =AX

(5.2)

and YM=BX,

(5.3)

the columns of YK will be approximately orthonormal, and in any case Yk YK would be wellconditioned. This means that the information lost in forming K is well preserved in the smaller

76

G. Frik and Th. Lunde Johnsen, Note on the ill·conditioned eigenvalue problem in elastic vibrations

eigenproblem

- -

Mx=KxA,

(5.4)

where

M= ytM Y.M

K= ytK Y.K

(5.5)

The eigenvalues of (5.4) and therefore of an accuracy comparable to those which would have been found by the Natural Factor approach using the factor U of (3.2). This remark indicates a further important aspect of the Natural Factor method, for instance in algorithms which are not compatible with the Natural Factor method (cf. the sectioning method in [23]).

Acknowledgement This work is a part of a larger project guided by Professor J.H. Argyris whose support is gratefully acknowledged.

References [1) J.H. Argyris, Energy theorems and structural analysis, Aircraft Eng. 26 (1954), Aircraft Eng. 27 (955). [2) J.H. Argyris, Energy theorems and analysis, Part I (Butterworths, London, 1960). [3] J.H. Argyris, Some results on the free·free oscillations of aircraft type structures, Research Report to the Boeing Company, Airplace Division (Nov. 1964). Also read at the International Symposium of the IUTAM on the Mechanics of Linear Vibrations, Societe Fran\(aise des Mecaniciens and British National Committee for Theoretical and Applied Mechanics, Paris, April 13-15, 1965, VI Session. Revue Fran\(aise de Mecanique, 3eme trimestre (1965) 59-73. [4] J.H. Argyris, Continua and discontinua, Opening paper to the Air Force Conference on Matrix Methods in Structural Mechanics at Wright Patterson Air Force Base, Dayton/Ohio (Oct. 1965) 1-198. [5) R.J. Guyan, Reduction of stiffness and mass matrix, AIAA Journal 3 (1965) 380. [6] J.H. Wilkinson, The algebraic eigenvalue problem (Clarendon Press, Oxford, 1965). [7] O.E. Briinlund, Computation of the Cholesky factor of a stiffness matrix direct from the factor of its initial quadratic form, ISD Report No. 142 (Stuttgart, May 1973). [8) J .H. Argyris and O.E. Briinlund, The natural factor formulation of the matrix displacement method, Computer Meth. Appl. Mech. Eng. 5 (1975) 97-119. [9) R.A. Rosanoff, J.F. Gloudeman and S. Levy, Numerical conditions of stiffness matrix formulations for frame structures, Proc. of Conf. on Matrix Methods in Structural Mechanics, AFFDL-TR-68-150, Wright-Patterson Air Force Base, Ohio, (1968) 1029-1060. [10] J.R. Roy, Numerical error in structural solutions, Journal of Structural Division, ASCE 97, No. ST4 (1971) 1039-1053. [11) DYNAN user's reference manual, ISD Report No. 97 (1971). [12] DYNAN·lecture notes with computational examples, ISD Report No. 109 (1971). [13) DYNAN Programmer's reference manual (1972). [14) G. Fix and R. Heiberger, An algorithm for the ill-conditioned general eigenvalue problem, SIAM J. Numer. Anal. 9 (1972) 78-88. [15) G. Peters and J.H. Wilkinson, Ax = BXA and the generalized eigenproblem, SIAM Num. Anal. 7 (1970) 479-492. [16) J.H. Wilkinson, Rundungsfehler (Springer, Berlin, Heidelberg, New York, 1969). [17) Th. Lunde Johnsen and J.R. Roy, On systems of linear equations of the form A tAx =b, error analysis and certain consequences for structural applications, Compo Meth. Appl. Mech. Eng. 3 (1974) 357 -374. [18] J.H. Argyris, Th. Lunde Johnsen, R.A. Rosanoff and J.R. Roy, On numerical error in the finite element method, Compo Meth AppI. Mech. Eng., to be published.

G. Frik and Th. Lunde Johnsen. Note on the ill-conditioned eigenvalue problem in elastic vibrations (19) J.H. Argyris and D.W. Scharpf, Some general considerations of the natural mode technique, Aero. J. Royal Aero. Soc. 73, Part 1, Small displacements, pp. 218-226 and Part 2, Large displacements, pp. 361-368 (1969). (20) R. Zurmiihl, Matrizen und ihre technische Anwendung, 4th ed. (Springer, 1964). [21) J.H. Wilkinson and C. Reinsch, Handbook for automatic computation 2, Linear algebra (Springer, 1971). (22) R.A. Rosanoff, Private communication (1972). (23) P.S. Jensen, The solution oflarge symmetric eigenproblems by sectioning, SIAM J. Numer. Anal. 9 (1972) 534 -545.

77