/* newblock.c - An idea for n-dimensional double precision data
   structure. This is pseudocode to describe the idea, does not
   compile.
*/

/* Structure for n-dimensional double precision data storage */

typedef struct 
{
  size_t dim;     /* data structure dimension */
  size_t * sizes; /* dimensions of system, e.g. {3,4,5} for real 3*4*5 tensor */
  double * data;  /* continuous data storage */
  size_t layout;  /* data layout type - LAST_INDEX_CONTINUOUS =
		     "column major" for matrices, or
		     FIRST_INDEX_CONTINUOUS = "row major" (default)
		     for matrices. Note: It might be possible to
		     extend this to support other storage schemes,
		     e.g. diagonal/banded storage
		     BANDED_INDEX_CONTINUOUS. */
} base_block;

/* Structure for a generic n-dimensional view object to access base
   block data. Access to data storage is done via this view
   object. There can be multiple views referring to the same data
   storage.
*/

typedef struct 
{
  size_t * firstcell; /* indices to first access cell (minimum indices) */
  size_t * lastcell; /* indices to last access cell (maximum indices) */
  size_t * stride; /* strides for each dimension for the view */ 
  base_block * block;
  double * data; /* pointer to block->data */
  size_t owner; /* indicator for view owning the data */
} base_view;

/* vector object based on base block and view */

typedef struct 
{
  base_view * view;
} vector;

/* matrix object based on base block and view */

typedef struct 
{
  base_view * view;
} matrix;

matrix * matrix_alloc (size_t n1, size_t n2) 
{
  /* Allocates matrix object (both view and storage objects) in
     row-major layout. */
  // 1. allocate base_block with dim=2, sizes={n1,n2}
  // 2. allocate base_view with firstcell={1,1}, 
  //    lastcell={n1,n2}, stride={0,0} (=whole matrix)
  // 3. set view owner=TRUE
  // 4. return base_view;
}

matrix * matrix_alloc_with_layout (size_t n1, size_t n2, size_t layout) 
{
  /* Allocates matrix object (both view and storage objects) in the
     specified layout. */
}

int vector_alloc_view_of_matrix (vector * v, matrix * m, 
				 size_t * firstcell, size_t * lastcell, 
				 size_t * stride) 
{
  /* Creates a new vector view of the matrix m */
  // 1. allocate base_view with given firstcell, lastcell and strides
  // 2. associate view with the storage
  // 3. set owner=FALSE
  // 4. return base_view;
}

/* Example code */

int main (void) 
{
  matrix * a;
  matrix * b;
  vector * v;
  int retval;

  /* allocate matrix with 3 rows and 4 columns */
  a = matrix_alloc (3, 4);

  /* User can change the layout between row- and column-major at any
     point. After change to column-major, those GSL functions that
     require row-major storage would raise error if the matrix is
     passed to them.
  */
  retval = matrix_set_layout (a, LAST_INDEX_CONTINUOUS); 

  {
    double temp[] = { 0, 10, 20, 01, 11, 21, 02, 12, 22, 03, 13, 23 };
    size_t templen = 12;
    retval = matrix_set_doublevec (a, temp, templen);
    /* This would initialize matrix
       a = [00 01 02 03
            10 11 12 13
            20 21 22 23] */
  }

  /* Set up a vector view of the diagonal */
  const size_t b0[] = { 1, 1 };
  const size_t e0[] = { 3, 3 };
  const size_t s0[] = { 1, 1 };
  retval = vector_alloc_view_of_matrix (v, a, &b0, &e0, &s0);
  
  /* Access cells in matrix through the vector */
  double val;
  val = vector_get(v, 1); // returns 00 from matrix a
  val = vector_get(v, 2); // returns 11 from matrix a
  val = vector_get(v, 3); // returns 22 from matrix a
  retval = vector_set(v, 2, 55); // sets value 55 in place of 11 to matrix a

  size_t l;
  l = vector_length(v); // returns 3;

  double diag[3];
  retval = vector_get_doublevec (&diag, v); // copies diagonal element values to diag
}

/* complex number matrix (one possible implementation) based on base
   block and view */

typedef struct 
{
  base_view * view;
} complex_matrix;

complex_matrix * complex_matrix_alloc (size_t n1, size_t n2) {
  /* Allocates matrix object (both view and storage objects) */
  // 1. allocate base_block with dim=3, sizes={n1,n2,2}, so that 
  //    the last dimension spans the real and imaginary parts
  // 2. allocate base_view with firstcell={1,1,1} and 
  //    lastcell={n1,n2,2} (=whole matrix)
  // 3. set view owner=TRUE
  // 4. return base_view;
}