#include <stdio.h>
#include <assert.h>
#include <gsl_complex.h>
#include <gsl_complex_math.h>
#include <gsl_math.h>
#include <gsl_eigen.h>
#include <gsl_cblas.h>

int N = 5;
gsl_complex Q_complex_couple = {{0.6,0.0}};
double      Q_real_couple    = 0.6;

#define MAT_SET(matrix,i,j,real,imag) 			\
do {							\
  gsl_complex *_c = gsl_matrix_complex_ptr(matrix,i,j);	\
  GSL_REAL(*_c) = real;					\
  GSL_IMAG(*_c) = imag;					\
}while(0)


//only works for square matrix
void write_matrix_complex(FILE *f, gsl_matrix_complex *m, char *form){
  int i,j,sz;
  char buffer[100];
  gsl_complex c;
  
  sprintf(buffer, "(%s,%s)", form, form);
  
  sz = m->size1;
  
  for(i=0; i<sz; i++){
    for(j=0; j<sz; j++){
      c = gsl_matrix_complex_get(m, i, j);
      fprintf(f, buffer, GSL_REAL(c), GSL_IMAG(c));
    }
    fprintf(f, "\n");
  }
}

void write_matrix(FILE *f, gsl_matrix *m, char *form){
  int i,j,sz;
  sz = m->size1;

  for(i=0; i<sz; i++){
    for(j=0; j<sz; j++)
      fprintf(f, form, gsl_matrix_get(m, i, j));
    fprintf(f, "\n");
  }
}




void fill_H_complex(gsl_matrix_complex *H){
  int i,j;
  
  for(i = 0; i < N; i++)
    for(j = 0; j < N; j++)
      gsl_matrix_complex_set(H, i, j, Q_complex_couple);
  
  for(i = 0; i < N; i++)
    MAT_SET(H,i,i, -3.0 - 1.0/(1.0+i), 0.0);
  
  
}

void fill_H(gsl_matrix *H){
  int i,j;
  for(i = 0; i < N; i++)
    for(j = 0; j < N; j++)
      gsl_matrix_set(H, i, j, Q_real_couple);
  
  for(i = 0; i < N; i++)
    gsl_matrix_set(H,i,i,-3.0 - 1.0/(1.0+i));
}


int main_complex(void){
  int i, j, res;
  gsl_matrix_complex *H, *H_copy;
  gsl_complex c;
  gsl_vector *eval;
  gsl_matrix_complex *evec;
  gsl_eigen_hermv_workspace *w;

  H      = gsl_matrix_complex_alloc(N,N);
  H_copy = gsl_matrix_complex_alloc(N,N);
  fill_H_complex(H);
  gsl_matrix_complex_memcpy(H_copy, H);
  
  
  eval = gsl_vector_alloc (N);
  evec = gsl_matrix_complex_alloc (N, N);
  w    = gsl_eigen_hermv_alloc (N);
  
  res = gsl_eigen_hermv(H_copy, eval, evec, w);
  if(res){
    printf("error: %s\n", gsl_strerror(res));
    exit(1);
  }
  gsl_eigen_hermv_free(w);
  gsl_eigen_hermv_sort (eval, evec,
			GSL_EIGEN_SORT_VAL_ASC);
 

  write_matrix_complex(stdout, H, "%8.4f");fprintf(stdout, "\n\n");

  for(i=0; i < N; i++)
    fprintf(stdout, "%15d%15e%15e\n", 
	    i,
	    GSL_REAL(gsl_matrix_complex_get(H,i,i)),
	    gsl_vector_get(eval,i));

  fprintf(stdout, "\n\n");
  write_matrix_complex(stdout, evec, "%8.4e");fprintf(stdout, "\n\n");
}

int main_real(void){
  int i, j, res;
  gsl_matrix *H, *H_copy;
  gsl_vector *eval;
  gsl_matrix *evec;
  gsl_eigen_symmv_workspace *w;
  
  H      = gsl_matrix_alloc(N,N);
  H_copy = gsl_matrix_alloc(N,N);
  fill_H(H);
  gsl_matrix_memcpy(H_copy, H);
  
  
  eval = gsl_vector_alloc (N);
  evec = gsl_matrix_alloc (N, N);
  w    = gsl_eigen_symmv_alloc (N);
  
  //  write_matrix_complex(stdout, H, "%10f");exit(1);
  
  res = gsl_eigen_symmv(H_copy, eval, evec, w);
  if(res){
    printf("error: %s\n", gsl_strerror(res));
    exit(1);
  }
  gsl_eigen_symmv_free(w);
  gsl_eigen_symmv_sort (eval, evec,
			GSL_EIGEN_SORT_VAL_ASC);
  
  
  write_matrix(stdout, H, "%12.4f");fprintf(stdout, "\n\n");
  
  for(i=0; i < N; i++)
    fprintf(stdout, "%15d%15e%15e\n", 
	    i,
	    gsl_matrix_get(H,i,i),
	    gsl_vector_get(eval,i));

  fprintf(stdout, "\n\n");
    
  write_matrix(stdout, evec, "%12.4e");fprintf(stdout, "\n\n");
}

int main(void){
  fprintf(stdout, "COMPLEX VERSION FAILS!!\n");
  fprintf(stdout, "=======================\n\n");
  main_complex();
  fprintf(stdout, "REAL VERSION WORKS!!\n");
  fprintf(stdout, "========================\n\n");
  main_real();
}