[Bug-gsl] bug in gsl_eigen_hermv

[Matthew B. Hastings]

Hi,

I think there is a bug in gsl_eigen_hermv which causes it to give nan if
there ar many very small entries or entries equal to zero.  I checked the
bug archives, and I found a similar bug report,
http://lists.gnu.org/archive/html/bug-gsl/2002-10/msg00002.html , for
gsl-1.2.  However, I am using gsl-1.10.  The bug does not seem to appear
on gsl-1.4.  I am using Mac OS X 10.4.9 on a PowerPC G4.

Here is some code the produces the error:
/*--------------------------------------------------------*/

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>

main()
{
gsl_matrix_complex *H;
gsl_eigen_hermv_workspace *w;
gsl_vector *eigs;
gsl_matrix_complex *U;
int i,j,size;
double x;
size=7;

H=gsl_matrix_complex_alloc(size,size);

for (i=0;i<size;i++)
for (j=0;j<size;j++)
{
gsl_matrix_complex_set(H,i,j,gsl_complex_rect(0.0,0.0));
}

gsl_matrix_complex_set(H,0,2,gsl_complex_rect(pow(.1,280.0),pow(.1,119.0)));
gsl_matrix_complex_set(H,2,0,gsl_complex_rect(pow(.1,280.0),pow(.1,119.0)));

gsl_matrix_complex_set(H,2,4,gsl_complex_rect(pow(.1,320.0),pow(.1,119.0)));
gsl_matrix_complex_set(H,4,2,gsl_complex_rect(pow(.1,320.0),pow(.1,119.0)));

gsl_matrix_complex_set(H,2,6,gsl_complex_rect(0.0,pow(.1,119.0)));
gsl_matrix_complex_set(H,6,2,gsl_complex_rect(0.0,pow(.1,119.0)));

gsl_matrix_complex_set(H,4,6,gsl_complex_rect(0.0,pow(.1,119.0)));
gsl_matrix_complex_set(H,6,4,gsl_complex_rect(0.0,pow(.1,119.0)));

gsl_matrix_complex_set(H,0,4,gsl_complex_rect(pow(.1,300.0),pow(.1,119.0)));
gsl_matrix_complex_set(H,4,0,gsl_complex_rect(pow(.1,300.0),pow(.1,119.0)));

printf("Here is the matrix H\n");
for (i=0;i<size;i++) {for (j=0;j<size;j++)
printf("(%lg,%lg)
",GSL_REAL(gsl_matrix_complex_get(H,i,j)),GSL_IMAG(gsl_matrix_complex_get(H,i,j)));
printf("\n");}

w=gsl_eigen_hermv_alloc(size);
eigs=gsl_vector_alloc(size);
U=gsl_matrix_complex_alloc(size,size);

gsl_eigen_hermv(H,eigs,U,w);

printf("\n");
printf("Here are the eigenvalues of H, note the nan\n");
for (i=0;i<size;i++) printf("%lg \n",gsl_vector_get(eigs,i));
}
/*--------------------------------------------------------*/

When I compile it with   gcc problem.c -lm -lgsl -lgslcblas
-Lgsl-1.10/.libs -Lgsl-1.10/cblas/.libs -Igsl-1.10 I get the output:

Here is the matrix H
(0,0) (0,0) (1e-280,1e-119) (0,0) (1e-300,1e-119) (0,0) (0,0)
(0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
(1e-280,1e-119) (0,0) (0,0) (0,0) (9.99989e-321,1e-119) (0,0) (0,1e-119)
(0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
(1e-300,1e-119) (0,0) (9.99989e-321,1e-119) (0,0) (0,0) (0,0) (0,1e-119)
(0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
(0,0) (0,0) (0,1e-119) (0,0) (0,1e-119) (0,0) (0,0)

Here are the eigenvalues of H, note the nan
nan
nan
nan
nan
nan
nan
nan

On the other hand, when I compile it with gsl-1.4 I get

Here is the matrix H
(0,0) (0,0) (1e-280,1e-119) (0,0) (1e-300,1e-119) (0,0) (0,0)
(0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
(1e-280,1e-119) (0,0) (0,0) (0,0) (9.99989e-321,1e-119) (0,0) (0,1e-119)
(0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
(1e-300,1e-119) (0,0) (9.99989e-321,1e-119) (0,0) (0,0) (0,0) (0,1e-119)
(0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
(0,0) (0,0) (0,1e-119) (0,0) (0,1e-119) (0,0) (0,0)

Here are the eigenvalues of H, note the nan
2.23607e-119
-2.23607e-119
1.50714e-135
-1.50714e-135
3.87304e-198
0
0


Thanks,
Matt Hastings