#include <stdlib.h>
#include <stdio.h>

#include <gsl/gsl_blas.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_version.h>

// Sets matrix elements to predefined values depending on row and column
// indices. The resulting matrix is complex Hermitian.
int
set_matrix (gsl_matrix_complex * m, const int size)
{
  int i, j;

  for (i = 0; i < size; i++)
    {
      gsl_matrix_complex_set (m, i, i, (gsl_complex) {{ (i+1) * (i+1), 0}});

      for (j = 0; j < i; j++)
        {
          gsl_matrix_complex_set (m, i, j, (gsl_complex) {{ (i+1) * (j+1),  j-5 }});
          gsl_matrix_complex_set (m, j, i, (gsl_complex) {{ (i+1) * (j+1), -j+5 }});
      }
    }

  return 0;
}

// Prints the matrix line-by-line to stdout.
int
print_matrix (const gsl_matrix_complex * m, const int nrows, const int ncols,
              const char * label)
{
  int i, j;
  gsl_complex z;

  fprintf (stdout, "%s\n\n", label);

  fprintf (stdout,
           "tda = %i\nsize1 = %i\nsize2 = %i\n\n",
           m->tda, m->size1, m->size2);

//  gsl_matrix_complex_fprintf (stdout, m, "%f");
//  fputc ('\n', stdout);

  for (i = 0; i < nrows; i++)
    {
      for (j = 0; j < ncols; j++)
        {
          z = gsl_matrix_complex_get (m, i, j);
          fprintf (stdout, "{% f, % f}", GSL_REAL(z), GSL_IMAG(z));
          fputc ((j < ncols - 1) ? '\t' : '\n', stdout);
        }
    }
  fputc ('\n', stdout);

  return 0;
}

// Prints the vector on single line to stdout.
int
print_vector (const gsl_vector * v, const int nsize, const char * label)
{
  int i;

  fprintf (stdout, "%s\n\n", label);

//  gsl_vector_fprintf (stdout, v, "%f");
//  fputc ('\n', stdout);

  for (i = 0; i < nsize; i++)
    {
      fprintf (stdout, "%f", gsl_vector_get (v, i));
      fputc ((i < nsize - 1) ? '\t' : '\n', stdout);
    }
  fputc ('\n', stdout);

  return 0;
}

// Prints the data block line-by-line assuming that matrix is complex-valued,
// but using double variables.
int
dump_matrix (const double * data, const int nrows, const int ncols,
             const int tda, const char * label)
{
  int i, j;

  fprintf (stdout, "%s (dump)\n\n", label);

  fprintf (stdout, "tda = %i\nsize1 = %i\nsize2 = %i\n\n", tda, nrows, ncols);

  for (i = 0; i < 2 * nrows; i += 2)
    {
      for (j = 0; j < 2 * ncols; j += 2)
        {
          fprintf (stdout, "{% f, % f}",
                   data[i * tda + j], data[i * tda + j + 1]);
          fputc ((j < 2 * (ncols - 1)) ? '\t' : '\n', stdout);
        }
    }
  fputc ('\n', stdout);

  return 0;
}

int
main (int argc, char ** argv)
{
  int i, j;
  gsl_complex z;

  const int tda = 7;
  const int size = 5;

  // Print GSL version.
  fprintf (stdout, "GSL v%s\n", gsl_version);
  fputc ('\n', stdout);

  gsl_matrix_complex * m = gsl_matrix_complex_calloc (tda, tda);
  gsl_matrix_complex * evec = gsl_matrix_complex_calloc (tda, tda);
  gsl_vector * eval = gsl_vector_calloc (tda);

  m->size1 = m->size2 = size;
  evec->size1 = evec->size2 = size;
  eval->size = size;

  // Fill matrix.
  set_matrix (m, size);

  // Print matrix.
  print_matrix (m, size, size, "matrix");

  // Diagonalize matrix.
  gsl_eigen_hermv_workspace * ws = gsl_eigen_hermv_alloc (tda);
  gsl_eigen_hermv (m, eval, evec, ws);
  gsl_eigen_hermv_free (ws);

  // Print diagonalized matrix.
  print_matrix (m, size, size, "matrix after diagonalization routine");

  // Print vector of eigenvalues.
  print_vector (eval, size, "vector of eigenvalues");

  // Print matrix of eigenvectors.
  print_matrix (evec, size, size, "matrix of eigenvectors");

  // Print full matrix of eigenvectors.
  evec->size1 = evec->size2 = tda;
  print_matrix (evec, tda, tda, "FULL matrix of eigenvectors");
  evec->size1 = evec->size2 = size;

  // Print matrix of eigenvector using doubles (just checking).
  dump_matrix (evec->data, size, size, tda, "matrix of eigenvectors");
  dump_matrix (evec->data, tda, tda, tda, "FULL matrix of eigenvectors");

  // Restore matrix.
  set_matrix (m, size);

  // Similarity transformation: conjg(trans(evec)).m.evec
  gsl_matrix_complex * tmp = gsl_matrix_complex_calloc (tda, tda);
  tmp->size1 = tmp->size2 = size;

  gsl_complex alpha = GSL_COMPLEX_ONE;
  gsl_complex beta = GSL_COMPLEX_ZERO;

  gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, alpha, m, evec, beta, tmp);
  gsl_blas_zgemm (CblasConjTrans, CblasNoTrans, alpha, evec, tmp, beta, m);

  gsl_matrix_complex_free (tmp);

  // Print transformed matrix.
  print_matrix (m, size, size, "matrix m after similarity transformation");

  // Write eigenvectors to disk.
  const char * filename = "evec.bin";
  FILE * stream = fopen (filename, "w");
  gsl_matrix_complex_fwrite (stream, evec);
  fclose (stream);

  // Check file size.
  int filesize_est = size * size * sizeof (gsl_complex);
  fprintf (stdout, "estimated size of evec file: %i byte\n", filesize_est);

  stream = fopen (filename, "r");
  fseek (stream, 0L, SEEK_END);
  int filesize_real = ftell (stream);
  fclose (stream);
  fprintf (stdout, "real size of evec file: %i byte\n", filesize_real);
  fputc ('\n', stdout);

  // Read eigenvectors from disk.
  gsl_matrix_complex_set_zero (evec);

  stream = fopen (filename, "r");
  gsl_matrix_complex_fread (stream, evec);
  fclose (stream);

  // Print matrix of eigenvectors.
  print_matrix (evec, size, size, "matrix of eigenvectors after reading");

  // Print full matrix of eigenvectors.
  evec->size1 = evec->size2 = tda;
  print_matrix (evec, tda, tda, "FULL matrix of eigenvectors after reading");
  evec->size1 = evec->size2 = size;

  // Brute-force check of the evec matrix on disk.
  double * evec_bf = (double *) calloc (size * size, sizeof (gsl_complex));

  stream = fopen (filename, "r");
  int nread = fread (evec_bf, sizeof (gsl_complex), size * size, stream);
  fclose (stream);

  fprintf (stdout, "number of elements read: %i\n\n", nread);

  dump_matrix (evec_bf, size, size, size, "matrix of eigenvectors after reading");

  free (evec_bf);


  gsl_matrix_complex_free (m);
  gsl_matrix_complex_free (evec);
  gsl_vector_free (eval);

  return EXIT_SUCCESS;
}