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

#include <gsl/gsl_complex.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_version.h>

int
main (int argc, char ** argv)
{
  fprintf (stdout, "GSL version: %s\n\n", gsl_version);

  int i, j;
  FILE * os;
  gsl_complex z;

  const int tda = 5;
  const int size = 3;

  gsl_matrix_complex * m = gsl_matrix_complex_alloc (tda, tda);

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

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

  fputs ("original complex-valued matrix:\n", stdout);
  //gsl_matrix_complex_fprintf (stdout, m, "%g");
  for (i = 0; i < tda; i++) // Use this version; it provides low-level access.
    {
      for (j = 0; j < tda; j++)
        {
          fprintf (stdout, "{%g, %g}",
                   m->data[2 * (i * tda + j)],
                   m->data[2 * (i * tda + j) + 1]);
          fputc ((j < tda - 1) ? '\t' : '\n', stdout);
        }
    }
  fputc ('\n', stdout);

  os = fopen ("matrix.bin", "w");
  gsl_matrix_complex_fwrite (os, m);
  fclose (os);

  gsl_matrix_complex_set_zero (m);

  os = fopen ("matrix.bin", "r");
  gsl_matrix_complex_fread (os, m);
  fclose (os);

  fputs ("matrix after writing to and reading from disk:\n", stdout);
  //gsl_matrix_complex_fprintf (stdout, m, "%g");
  for (i = 0; i < tda; i++) // Use this version; it provides low-level access.
    {
      for (j = 0; j < tda; j++)
        {
          fprintf (stdout, "{%g, %g}",
                   m->data[2 * (i * tda + j)],
                   m->data[2 * (i * tda + j) + 1]);
          fputc ((j < tda - 1) ? '\t' : '\n', stdout);
        }
    }
  fputc ('\n', stdout);

  gsl_matrix_complex_free (m);

  return EXIT_SUCCESS;
}