problem with gsl_eigen_hermv(), possibly due to bug ingsl_linalg_complex_householder_transform()

[Steve Martin <spm at zippy dot physics dot niu dot edu>]

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

Attachment: testbug.c
Description: Text document