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

* problem with gsl_eigen_hermv(), possibly due to bug in gsl_linalg_complex_householder_transform()
@ 2002-10-15 15:26 Steve Martin
  2002-10-17  0:11 ` Gert Van den Eynde
  0 siblings, 1 reply; 5+ messages in thread
From: Steve Martin @ 2002-10-15 15:26 UTC (permalink / raw)
  To: gsl-discuss

[-- Attachment #1: Type: TEXT/PLAIN, Size: 1771 bytes --]

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 #2: Type: TEXT/PLAIN, Size: 3865 bytes --]

/* 
This mini-program purports to show that gsl_eigen_hermv() fails when the
Hermitian matrix contains some 0 entries. In particular, this occurs
whenever the matrix is diagonal. This is illustrated here with the matrix
{{2,0,0},{0,3,0},{0,0,4}}. Adding an off-diagonal part to the 1,2 and 2,1
entries, or to the 1,3 and 3,1 entries in the matrix does NOT avoid the
problem. However, adding even a very small non-zero part to the 2,3 and
3,2 entries does avoid the problem.

The same problem occurs for diagonal 2x2 matrices (not shown).

To compile: gcc testbug.c -lgsl -lm -lgslcblas
*/

#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <gsl/gsl_eigen.h>

#define size 3
#define eps 0.000000000000001

int findeigenvecs(complex double hmatrix[size][size]);

int main(void){

complex double hmatrixdiagdata[size][size] = 
                  {{2., 0., 0.},
                   {0., 3., 0.},
                   {0., 0., 4.}};

complex double hmatrix12data[size][size] = 
                     {{2., 1., 0.},
                      {1., 3., 0.},
                      {0., 0., 4.}};

complex double hmatrix13data[size][size] = 
                     {{2., 0., 1.},
                      {0., 3., 0.},
                      {1., 0., 4.}};

complex double hmatrixeps23data[size][size] = 
                     {{2., 0., 0.},
                      {0., 3., eps},
                      {0., eps, 4.}};

    printf("\n");
    printf("First find the eigenvectors of the diagonal matrix\n"); 
    printf("{{2.,0.,0.},{0.,3.,0.},{0.,0.,4.}}:\n");

    findeigenvecs(hmatrixdiagdata);

    printf("Uh-oh, nan appears in the last entry of the last eigenvector!\n");
    printf("That should not have happened.\n");

    printf("\n");
    printf("Now, try the same matrix but with non-zero 1,2 and 2,1 entries:\n");
    printf("{{2.,1.,0.},{1.,3.,0.},{0.,0.,4.}}:\n");

    findeigenvecs(hmatrix12data);

    printf("And we see that the problem still occurs.\n");
    printf("\n");

    printf("Now, try the original matrix but with non-zero 1,3 and 3,1 entries:\n");
    printf("{{2.,0.,1.},{0,3.,0.},{1.,0.,4.}}:\n");

    findeigenvecs(hmatrix13data);

    printf("And we see that the problem still occurs, in fact twice.\n");

    printf("\n");

    printf("Now, include a very small part in the 2,3 and 3,2 positions:\n");

    findeigenvecs(hmatrixeps23data);

    printf("And we see that the problem does not occur in this case.\n");
    printf("\n");

    return 0;
}				/* End of main */

int findeigenvecs(complex double hmatrix[size][size])
{
    int i, j;
    gsl_complex tempconvert;
    gsl_matrix_complex *hermmatrix = gsl_matrix_complex_alloc(size, size);
    gsl_matrix_complex *eigenvecs = gsl_matrix_complex_alloc(size, size);
    gsl_vector *eigenvals = gsl_vector_alloc(size);
    gsl_eigen_hermv_workspace *w = gsl_eigen_hermv_alloc(size);

/* Translate the Hermitian matrix hmatrix into gsl form as hermmatrix */
    for (i = 0; i < size; i++)
	for (j = 0; j < size; j++) {
	    GSL_SET_COMPLEX(&tempconvert, creal(hmatrix[i][j]),
			    cimag(hmatrix[i][j]));
	    gsl_matrix_complex_set(hermmatrix, i, j, tempconvert);
	};

/* Here we go. */
    gsl_eigen_hermv(hermmatrix, eigenvals, eigenvecs, w);

/* Now print out the eigenvectors we have found. */
    for (i = 0; i < size; i++) {
        printf("eigenvector %d =\n",i+1);
	for (j = 0; j < size; j++) {
	    printf("%f+I*%f    ",
		   GSL_REAL(gsl_matrix_complex_get(eigenvecs, j, i)),
		   GSL_IMAG(gsl_matrix_complex_get(eigenvecs, j, i)));
	}
	printf("\n");
    }

/* Free up, just to be good. */
    gsl_eigen_hermv_free(w);
    gsl_matrix_complex_free(eigenvecs);
    gsl_vector_free(eigenvals);
    gsl_matrix_complex_free(hermmatrix);

    return 0;
}

^ permalink raw reply	[flat|nested] 5+ messages in thread