Re: gsl container designs

["Robert G. Brown" <rgb at phy dot duke dot edu>]

On Wed, 6 Jan 2010, Tuomo Keskitalo wrote:

> Here is a crazy idea: If we take Gerards design 1, would it be too absurd to 
> use a horrible big wrapper struct like
>
> typedef struct {
>  gsl_marray *m;
>  gsl_const_marray *cm;
>  gsl_marray_1 *m1;
>  gsl_const_marray_1 *cm1;
>  gsl_marray_2 *m2;
>  gsl_const_marray_2 *cm2;
>  gsl_marray_3 *m3;
>  gsl_const_marray_3 *cm3;
>  gsl_vector *vec;
>  gsl_const_vector *cvec;
>  gsl_matrix *mat;
>  gsl_const_matrix *cmat;
>  ...
> } gsl_container;
>
> and make everything a gsl_container, and then always use them like
> &a.mat or &b.vec? Maybe the const types could be in another 
> gsl_const_container struct..
>
> One downside is that this would effectively mask the type of a or b in 
> program code (gsl_container a; // is it a vector, matrix, or marray?), which 
> might make reading code a pain..

IMO the primary reason for "fixing" vector/matrix objects in the GSL is
to enable the single vector that represents the actual data in the
object to be presented to the user in a matrix-intuitive way, even at
some cost in code efficiency (and to try to minimize that cost or
provide further benefits to counterbalance it by enabling efficient
linear algebra routines to be matched to it, for example).

Any of us who have coded for a decade or two know perfectly well how to
address a vector via explicit displacement arithmetic, burying horrible
constructs in one's code to compute the displacement of an i,j,k element
in a vector containin aa "3d matrix".  Or, one can use C pointer
constructs that are equally unreadable for the same general purpose.  We
have also had to DEBUG this sort of thing, which is extremely
unpleasant, or understand somebody else's displacement code, which is
often unpleasant to the point of being nightmarish.

What I want is extremely simple:

* A way to construct/allocate a matrix of arbitrary dimension (out to some
(un)reasonable max, e.g. 10, likely to be higher than anybody needs for
coding purposes at least at this point, and deconstruct/free the matrix
when I'm done with it, ideally with a single call in each case, ideally
for any type of variable.  I'm perfectly comfortable with requiring a
pointer and cast to establish the type, but I want to be able to use it
just like m[i][j][k]...[r][s] when all is said and done.

* The specific dimension bounds need to be user specified.  If I want to
make a 6 dimensional triple triangular array to hold e.g. 3j symbols, I
don't want to have to allocate a rectangular array most of which is
empty wasted space.  I'm not talking about sparse matrix here, just
non-rectangular arrays.  Angular momentum indices are perfect (and
common) examples.

* Access to the underlying actual data vector (non-opaque, in other
words) so it can be passed to e.g. ODE solvers and the like that can
only be sensibly written in a general way for a vector.  I want to be
able to pass the ODE solver the vector and its dimension (and perhaps
an auxiliary similar vector for the derivatives with similar packing
and indexing), but I want to be able to write the ODE deriv routine 
using the readable matrix form.

None of this is particularly difficult to accomplish.  I do it all the
time in my own code.  In the end, one can integrate sets of ODEs where
the derivatives are indexed by angular momentum l,m, where the deriviative
expression contains sums over l',m', l",m" with precomputed tables of
e.g. Gaunt numbers, and where the derivative expression is HUMAN
READABLE and directly comparable to an expression in a physics paper or
textbook.

My problem in this regard with the GSL is that so far, the matrix and
vector constructs are pretty much useless -- they are far away from the
standard way one usually does this sort of thing in C using pointers or
just allocating a vector or matrix.

One of the ADVANTAGES of C is that it permits one to work pointer magic
to improve the readability (and sometimes execution efficiency) of code.
Opaque types and set/get calls are dark evil in both regards.  I just
want to be able to do (double) I[l][m][l'][m'][l"][m"]aand/or (int)
epsilon[i][j][k] transparently and reasonably efficiently, because both
of these occur all the time in real numerical code tht I've written.  I
can do it myself and still use the GSL, and have in the past done so.
For ODE solvers this isn't much of a problem.  For linear algebra,
however, it starts being a problem.  Then issues like column-major
become important and have to be solved in a standardized way.  This
would make even my own code more efficient, as well, as most of what I do
with the I() matrix above is basically matrix multiplication or
contraction -- summing out shared indices -- and linear algebra routines
that operate on correctly allocated structures can be 2x-3x more
efficient than ones that sum across big steps in the underlying vectors.

   rgb

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

Robert G. Brown	                       http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb@phy.duke.edu