Hi,

Say you had to solve A * x = lambda B * x. Following the book's instructions one could write a code involving choldc, tred2, tqli and other minor bits and pieces to connect the routines together.

I have done so, actually for the slightly more complicated case of complex Hermitian matrices, and it works perfectly in my tests with big (200x200+) random matrices. But in the actual problem I want to apply this to, where A and B are 5x5 real near-singular matrices, the process fails for one of the five eigenvalues. It correctly gets the triple one, and one of the two single eigenvalues, but fails for the other one.

In particular, the Cholesky decomposition of B returns a lower triangular matrix with a diagonal of 1 1 1 1 1d-7 (i.e., det(B)~1d-14 from this). I cannot but associate the four 1's with the correct eigenvalues and the 1e-7 with the one that goes wrong.

I've also tried linear combinations of A and B (e.g. A'=A+kB and B'=B+A'/w, which has the same eigenvectors and related eigenvalues), but this doesn't improve things. I'm sure the problem is solvable, though.

My question is, is it possible to tackle this problem as I've outlined, or do I need to use a different method when A and/or B are nearly singular?

Thanks,
Pablo Lopez