#include <gsl/gsl.h>
#include <marray.h>

#define RANK 3
#define DIM1 5
#define DIM2 2
#define DIM3 4


int main (void)
 {
  size_t dimensions[3] = { DIM1, DIM2, DIM3 };  /* to pass marray dimensions */
  size_t i, j, k;        /* marray indices */
  size_t indices[3];     /* to pass marray indices */

  /* Allocate space for the marrays */
  marray * a = marray_calloc (RANK, dimensions);
  marray * b = marray_calloc (RANK, dimensions);
  marray * t = marray_calloc (RANK, dimensions);

  /* Fill them */
  for (i = 0; i < DIM1; i++)
    {
      for (j = 0; j < DIM2; j++)
        {
          for (k = 0; k < DIM3; k++)
            {
              counter++;
              indices[0] = i;  indices[1] = j;  indices[2] = k;
              marray_set (a, indices, 3 + i + 5 * j + 2 * k);
              marray_set (b, indices, 3 + 2 * i + 4 * j + k);
            }
        }
    }

  /*
   * Do some operations
   */

  /* Sum */
  marray_memcpy (t, a);
  marray_add (t, b);     /* t = a + b */
  marray_fprintf (stdout, t, "%g");

  /* Element division */
  marray_memcpy (t, a);
  marray_div_elements (t, b);  /* t = a ./ b */
  marray_fprintf (stdout, t, "%g");

  marray_free (a);
  marray_free (b);

  marray_free (t);

  return 0;
}