/* 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 #include #include #include #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; }