Algorithm 297: Eigenvalues and Eigenvectors of the symmetric system

p r o c e d u r e eigensolve (a, b, x, n, nondef); value n; i n t e g e r n; l a b e l nondef; a r r a y a, b, x; c o m m e n t solves the equation (A --XB)X = 0 for symmetric A, B in a, b[1 :n, l:n] provided one of these is either positive or negative-definite.B is decomposed symmetrically so that B = LL' and the equation transformed to (C--X/)Y = 0 wherethe original equation is rearranged as ( B -( 1 / h ) A ) X = 0 and solved as such for positive-definite A or as ( B --( -l / X ) ( --A ) ) X = 0 for negative-definite A. Failure to achieve one useful transformation from the four possibilities leads to exit via the label nondef.The procedure calls procedure symmetric QR 2 [P. A. Businger, Algorithm 254, Eigenvalues and eigenvectors of a real symmetric matrix by the QR method.Comm.A C M 8 (April, 1965), 218-219] to evaluate the roots and vectors of ( C -h i ) Y = 0. T h a t procedure leaves untouched the strictly upper triangle of C. In conformity with this, eigensolve preserves the strictly upper triangles of A and B.If, before entry to eigensolve, the user saves the diagonals of A, B, both these arrays may, if necessary, be fully restored after exit.On exit from the procedure the eigenvalues occupy the diagonal elements a [i, i] with the eigenvectors in corresponding columns of x[l:n, l:n]; b e g i n i n t e g e r i, j, k, jlessl, ilessl, adi, bdi; r e a l t, sum, xij, length; B o o l e a n recip ; p r o c e d u r e LCHOLESKI(a, n, fail) ; value n; i n t e g e r n; label fail; a r r a y a; c o m m e n t performs the symmetric decomposition A = LL' for positive definite A in all :n, 1 :n].The lower triangle of A is overwritten by L. The strictly upper triangle of A is intact.For nonpositive-definite A the procedure exits via label parameter fail; b e g i n i n t e g e r i, j, k, jlessl; r e a l ajj, ajk, aij; jlessl := O;for j := 1 s t e p 1 u n t i l n do b e g i n ajj := a[j, j];for k := 1 s t e p 1 u n t i l jlessl do b e g i n ajk := a[j, k]; ajj := ajj -ajk × ajk e n d ;if ajj < 0.0 t h e n go to fail; ajj := a[j, j] := sqrl(ajj);for i := j + 1 s t e p 1 u n t i l n do b e g i n aij := a[i, j]; for k := 1 s t e p 1 u n t i l jlessl do aij := aij --a[i, k] X a[j, k]; a[i, j] := aij/ajj