Re: containers tentative design summary

[James Bergstra <bergstrj@iro.umontreal.ca>]

On Mon, Oct 5, 2009 at 6:56 PM, Gerard Jungman <jungman@lanl.gov> wrote:
> On Mon, 2009-10-05 at 10:50 -0400, James Bergstra wrote:
>> Two comments:
>>
>> I'm a bit rusty with my C structs... but you would need two distinct
>> static classes to have const and non-const data pointers for your view
>> right?
>
> Yes, sorry. I left that out.
>
>
>> Also, it sounds like writing code that will work for a tensor of any
>> rank (e.g. add two tensors together) might be either tedious or
>> impossible.  I recognize that part of the problem is the lack of
>> templating and polymorphism, but it would at least be comforting to
>> see just how bad the situation is via a use case or two in the design
>> documentation.   I (naively?) fear that to get good performance will
>> require a whole library of functions for even the most basic of
>> operations.
>>
>> gsl_marray_add_1_0( gsl_marray_1, double );
>> gsl_marray_add_1_1( gsl_marray_1, gsl_marray_1);
>> gsl_marray_add_1_2( gsl_marray_1, gsl_marray_2);
>> gsl_marray_add_2_2(... )
>> ...
>>
>> gsl_marray_sub_1_0( ... )
>>
>> Maybe a system of macros could be designed to help here, but it sounds
>> like it will never be as easy as writing a couple of for-loops.
>
>
> First, I want to be sure that we distinguish the multi-array
> case from the vector/matrix/tensor case. Algebra for vectors and
> matrices is well-defined and already on place; we only need to
> re-write some wrappers, etc. It is handled by BLAS calls. As I
> mentioned, this places a constraint on matrices, that the
> fast traversals be of unit stride. It seems like we just
> have to live with that constraint, for the matrix type.
>
> Addition, subtraction, scaling of multi-arrays is not hard
> because it is only defined within the same rank. So there
> is only a linear complexity, and no combinatorial explosion
> in the interfaces, for these operations. Of course, there
> are issues of shape conformance at run-time, but that is
> also true for vectors and matrices.
>
> That leaves multi-array algebra as an open question. By algebra,
> I mean technically the "cross-rank" operations, which form
> some kind of general multi-linear algebra. Sounds pretty
> hairy, as you suspect.
>

Another idea here would be to implement general-cases in terms of
structures with a dynamic rank right?  The rank does not need to be
part of the struct, but can be a variable.

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.
One disadvantage of such a structure I think you'll agree, is the
internal pointers, which interfere with convenient stack allocation.

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.

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

James
-- 
http://www-etud.iro.umontreal.ca/~bergstrj