Re: problem with gsl_eigen_hermv(), possibly due to bug in gsl_linalg_complex_householder_transform()

[Gert Van den Eynde <gvdeynde at sckcen dot be>]

Hi,

Just to let you know I can reproduce this behaviour. I think we can classify this as a bug. I'll try to have a look at it later (I don't have my Golub and Van Loan at hand here), but if someone else has a solution, feel free....


Bye,

Gert


On Tue, 15 Oct 2002 17:26:32 -0500 (CDT)
Steve Martin <spm@zippy.physics.niu.edu> wrote:

> Hi,
> 
> I think there may be a bug in gsl 1.2 which affects the computation of
> eigenvectors of Hermitian matrices using gsl_eigen_hermv(), but quite
> possibly other things too.
> 
> Everything works fine except when the Hermitian matrix has some zero
> entries. Then the eigenvectors returned are NaN in one or more entries.  
> In particular, the bug always seems to arise when the original matrix is
> diagonal.
> 
> Attached is a small program which illustrates the problem for the 3x3
> matrices {{2,0,0},{0,3,0},{0,0,4}} and {{2,1,0},{1,2,0},{0,0,4}} and
> {{2,0,1},{0,3,0},{1,0,4}}. In contrast, I have checked that Hermitian
> matrices with generic non-zero entries do not see the problem. Also, not
> all matrices with some zero entries have a problem, for example
> {{2,0,0},{0,3,x},{0,x,4}} is handled correctly for any non-zero value of
> x.
> 
> I do not have a patch, but I think the problem may be as follows:
> gsl_eigen_hermv() in eigen/hermv.c calls the function
> gsl_linalg_hermtd_decomp() in linalg/hermtd.c, which in turn calls
> gsl_linalg_complex_householder_transform() in linalg/householdercomplex.c 
> In the latter function, the variable beta_r is divided by, but it seems 
> to be 0 in the bad cases mentioned above. 
> 
> Information on my setup:
> I am using gsl 1.2, compiled without optimization. To be specific,
> I changed -O2 to -O0 on lines 1268 and 1274 of configure. Otherwise,
> the install was the generic one: make check produced no failures.
> My hardware is a dual pentium III running Linux kernel 2.4.9-34smp.
> The output of gcc -v is:
> gcc version 2.96 20000731 (Red Hat Linux 7.2 2.96-108.7.2)
> (I'm running standard Red Hat 7.2 with all updates applied.)
> 
> I'll be happy to try to provide more info if needed. 
> 
> Thanks in advance for whatever you can tell me about this,
> Steve Martin
>