#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));
          }
      }
  }
}