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

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

#define N 20

/**
 * Create two 2N-by-2N matrices A and B whose elements are given by:
 *   A(i,j) = i*N + j
 *   B(i,j) = i + j*N
 * Suppose we want to operate on A_ul and B_lr, the upper-left N-by-N submatrix
 * of A and the lower-right N-by-N submatrix of B, respectively, and suppose we
 * want to use submatrices instead of writing some loop-like code.
 */
int
main (int argc, char ** argv)
{
  size_t i, j;
  gsl_matrix *A, *A_ul, *B, *B_lr;

  fprintf (stderr, "GSL version: %s\n\n", gsl_version);

  /* Initialize A and B */
  A = gsl_matrix_alloc (2*N, 2*N);
  for (i = 0; i < 2*N; i++)
    {
      for (j = 0; j < 2*N; j++)
        gsl_matrix_set (A, i, j, i * 2*N + j);
    }

  B = gsl_matrix_alloc (2*N, 2*N);
  for (i = 0; i < 2*N; i++)
    {
      for (j = 0; j < 2*N; j++)
        gsl_matrix_set (B, i, j, i + j * 2*N);
    }

  /* Dump gsl_matrix structs. */
  fprintf (stderr, "A = %p\n", A);
  fprintf (stderr, "A->data = %p\n", A->data);
  fprintf (stderr, "A->nrows = %lu\n", A->size1);
  fprintf (stderr, "A->ncols = %lu\n", A->size2);
  fputc ('\n', stderr);

  fprintf (stderr, "B = %p\n", B);
  fprintf (stderr, "B->data = %p\n", B->data);
  fprintf (stderr, "B->nrows = %lu\n", B->size1);
  fprintf (stderr, "B->ncols = %lu\n", B->size2);
  fputc ('\n', stderr);

  /* Get upper-left and lower-right submatrices. */
  {
    gsl_matrix_view view1 = gsl_matrix_submatrix (A, 0, 0, N, N);
    A_ul = &(view1.matrix);
  }

  {
    gsl_matrix_view view2 = gsl_matrix_submatrix (B, N, N, N, N);
    B_lr = &(view2.matrix);
  }

  /* Dump gsl_matrix structs of submatrices. */
  fprintf (stderr, "A_ul = %p\n", A_ul);
  fprintf (stderr, "A_ul->data = %p\n", A_ul->data);
  fprintf (stderr, "A_ul->nrows = %lu\n", A_ul->size1);
  fprintf (stderr, "A_ul->ncols = %lu\n", A_ul->size2);
  fputc ('\n', stderr);

  fprintf (stderr, "B_lr = %p\n", B_lr);
  fprintf (stderr, "B_lr->data = %p\n", B_lr->data);
  fprintf (stderr, "B_lr->nrows = %lu\n", B_lr->size1);
  fprintf (stderr, "B_lr->ncols = %lu\n", B_lr->size2);
  fputc ('\n', stderr);

  gsl_matrix_free (A);
  gsl_matrix_free (B);

  return 0;
}