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