Dear all,
a few weeks ago, I have sent you this bug report already without giving
you enough information. Now, I would like to provide the missing
information and correct some statements.
After using the procedure gsl_eigen_hermv to diagonalize a matrix, I
suspect that the sign convention for the eigenvectors should be
reconsidered. More specifically, I diagonalize a matrix which has only a
few non-zero entries, namely at m(2i,2i+1) and m(2i+1,2i) for all i. In
my case, I diagonalize a 8x8 matrix with real entries, e.g. 1.0.
The eigenvectors should be of the following form (not normalized):
v1=(1,1,0,0,0....)
v2=(-1,1,0,0,0,...),
other eigenvectors can be obtained by shifting the non-zero coefficients
of v1 and v2 by two places to the right.
However, diagonalizing the 8x8 matrix described above, I obtain the
eigenvectors as described above except the last one which is
v8=(0,0,0,0,0,0,1,-1) instead of (0,0,0,0,0,0,-1,1), i.e., multiplied
with (-1). Therefore, I think that the sign convention is either not
implemented correctly or the convention used is not broad enough.
Please find in the attachment to this e-mail a simple compilable code
which demonstrates this problem. I have checked the described results
for GSL versions 1.15 and 1.16.
Kind regards,
Walter
#include <iostream>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
int main()
{
int Ns = 3;
int Nwf = (int) pow(2,Ns);
//initialisation of the matrix
double Hamil[] = { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
};
gsl_matrix_complex_view Ham = gsl_matrix_complex_view_array(Hamil,Nwf,Nwf);
//print the matrix
{
std::cout << std::endl << "matrix" << std::endl;
for (int mp1=0;mp1<Nwf;mp1++)
{
for (int mp2=0;mp2<Nwf;mp2++)
{
printf("%5.2f",GSL_REAL(gsl_matrix_complex_get(&Ham.matrix,mp1,mp2)));
printf("%5.2fi ",GSL_IMAG(gsl_matrix_complex_get(&Ham.matrix,mp1,mp2)));
}
std::cout << std::endl;
}
}
gsl_vector *eval = gsl_vector_calloc (Nwf);
gsl_matrix_complex *evec = gsl_matrix_complex_calloc (Nwf,Nwf);
gsl_eigen_hermv_workspace * w = gsl_eigen_hermv_alloc (Nwf);
gsl_eigen_hermv (&Ham.matrix, eval, evec, w);
gsl_eigen_hermv_free (w);
//print the eigenvalues and eigenvectors of the matrix
{
for (int i = 0; i < Nwf; i++)
{
double eval_i = gsl_vector_get (eval, i);
gsl_vector_complex_view evec_i = gsl_matrix_complex_column (evec, i);
printf ("eigenvalue = %g\n", eval_i);
printf ("eigenvector = \n");
for (int j = 0; j < Nwf; ++j)
{
gsl_complex z = gsl_vector_complex_get(&evec_i.vector, j);
printf("%g + %gi\n", GSL_REAL(z), GSL_IMAG(z));
}
}
}
}
