Re: new double precision data structure?

[Tuomo Keskitalo <Tuomo dot Keskitalo at iki dot fi>]

On 09/16/2009 03:49 AM, Gerard Jungman wrote:

> It smells like a view, and indeed it is. But it is a notion of
> view which is more central to the design than the current one.
> Everything should be a view.

I think we really need to sort the data structure out now.

I'm not sure what you have in mind, but to get some discussion going on, 
I've attached an (incomplete and untested) idea for a general 
n-dimensional double precision data structure. Is this anything like 
what you pictured?

I'd like to hear people's comments on this. Would something like this be 
of use? Feel free to give your full critique.

--
Tuomo.Keskitalo@iki.fi
http://iki.fi/tuomo.keskitalo

/* 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;
}