Re: containers tentative design summary

[Gerard Jungman <jungman at lanl dot gov>]

On Mon, 2009-10-05 at 19:45 -0400, James Bergstra wrote:

> I'm thinking of something like
> 
> struct marray { int rank; int * dims; int * strides; dtype * base; };
> marray_add_inplace( const marray a, const marray b, marray c);  //N.B.
> the passing on the stack
> 
> Inside that function, we can check for special cases of interest and
> optimize them as necessary.

This is not so clear. You can't have a test inside the indexing
functions; that would obviously kill you. So we still need a way
to pick the correct indexing function at compile-time; I think
this means that we need the fixed-rank types. And the indexing
calculation is best done with 'dims' and 'strides' on the stack,
rather than through the pointer indirection.

Maybe you are suggesting the above as an implementation detail,
to be wrapped by the fixed-rank types? That would allow for
compile-time selection of the indexing, but it would still
probably force an indirection through those pointers.


> One disadvantage of such a structure I think you'll agree, is the
> internal pointers, which interfere with convenient stack allocation.

Yes. I rejected internal pointers for that reason.


> A somewhat ugly alternative (which I am nonetheless suggesting) would
> be the following...
> 
> { int rank; dtype * base; int dim0, dim1, dim2; int str0, str1, str2;
> int * extra_dims; int * extra_strides; };
> 
> For low-rank multiarrays, this structure can be manipulated on the
> stack, so it is possible to write expressions like
> 
> marray_add_inplace(marray_reshape2(a, 5, 50), marray_slice(b, 17),
> marray_transpose(c));  // add a reshaped view of a to the 17th
> matrix-slice of b, storing result in a transposed view of c.
> 
> At the same time, arbitrarily high ranks are supported if the user is
> willing to put up with the syntax of heap-based manipulation.

This is closer to what I have done, although I just separated the
fixed-rank and the arbitrary rank (heap-based) cases. In the end,
that seemed like the cleanest choice.


> (Of course we can haggle over how many extra ranks can be built into
> the static part of the structure.)

The way I have done it, I avoid this question. Here are some snippets
from what I have, just so you get an idea of how I am thinking about
it. I'll work downwards, from the multi-arrays, down to the
implementation details.

First, we must have separate types for the fixed rank cases, in order
to insure that we select the correct indexing function at compile-time.
So...

  typedef struct { ... } gsl_marray_1;  // rank 1
  typedef struct { ... } gsl_marray_2;  // rank 2
  typedef struct { ... } gal_marray_3;  // rank 3
  etc.


The basic ingredient for indexing is a dimension and stride
in an index place. So...

  typedef struct
  {
    size_t  dimens; // dimension in this index place
    size_t  stride; // addressing stride for this index
  }
  gsl_marray_idx;


Given this, the marray declarations can look like:

  typedef struct
  {
    int            rank;
    double *       data;
    gsl_marray_idx idx[1];
    ...
  }
  gsl_marray_1

  typedef struct
  {
    int            rank;
    double *       data;
    gsl_marray_idx idx[2];
    ...
  }
  gsl_marray_2

  etc.


Now, at the level of implementations, we write everything in
terms of (rank, gsl_marray_idx[]), avoiding an intermediate type.
We can freely mix generic versions of functions and specializations.
Generic versions handle the general case, for operations which are
not performance critical (like formatted i/o and stuff...).
Specializations handle the things that matter, like indexing.

Examples:

  // some "low-level" io implementation function;
  // handles the access generically
  gsl_marray_fread_raw(FILE * f, double * data, int rank,
                       const gsl_marray_idx idx[]);


  // indexing (offset calculation) requires specialization
  inline
  size_t
  gsl_marray_idx_1_offset(const gsl_marray_idx idx[], const int index[])
  { return index[0]*idx[0].stride; }

  inline
  size_t
  gsl_marray_idx_2_offset(const gsl_marray_idx idx[], const int index[])
  { return index[0]*idx[0].stride + index[1]*idx[1].stride; }

  etc.

  inline
  gsl_marray_1_offset(const gsl_marray_1 * ma, const int index[])
  { return gsl_marray_idx_1_offset(ma->idx, index); }


You get the idea. These are taken (almost) verbatim from the code
that I have written so far. Finally, for the generic case we have
something like:

  typedef struct
  {
    int              rank;
    gsl_marray_idx * idx;
    double *         data;
    ...
  };

with the obvious heap allocating constructor, etc. The same
low-level generic functions as above are also used to implement
the class-specific functions here.

Anyway, the whole thing seems to be layered properly, and there
are no "tricky" bits. No rocket science. Just following my nose.


I will shortly produce a complete package for people to examine. But
you know how it is. As David Byrne once sang, "you can't see it 'til
it's finished...".

In the meantime, keep suggesting stuff. Discussion always helps.

--
G. Jungman