This is the mail archive of the
gsl-discuss@sourceware.org
mailing list for the GSL project.
Re: generalized eigensystems
> Here is an interesting test case, it aborts with a "shouldn't get
> here" error in extended precision but works in double precision (via
> GSL_IEEE_MODE).
Well, I've found out that the algorithm is failing under -O2 when
computing the eigenvalues of a 2-by-2 generalized block. To
do this I had copied the equations directly from Moler's paper
which is the same thing EISPACK did, but LAPACK has a fancier
routine which computes scaling factors and checks for overflow/underflow
etc. Apparantly the bug shows up when looking at a block with
"almost" complex eigenvalues - ie: a real eigenvalue block with
a determinant very close to 0. In the no-optimization case it
computes the determinant as slightly positive, and with optimization
on it gets computed as slightly negative which must be the result
of some sort of underflow condition.
Looks like I'll have to port over the lapack routine...bleh. I should
have a patch ready within the next week or so.
Patrick Alken