Re: gslclapack

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

On Wed, 29 Mar 2006, Gerard Jungman wrote:

> One more small point (thanks to Jim Amundson for pointing this out):
> ATLAS not only has blas functions, but also has a few lapack routines.
> I don't know what their intention is/was in this direction, but this
> may add something to the stew.

If we're spicing up the stew, this might be a good time to reconsider
the way matrices are handled in GSL in general.  Some time ago there was
a discussion of multidimensional (high rank) tensors -- 3-9 dimensions,
basically, of use in e.g. relativistic theories, lattice gauge theories,
a whole bunch of physics -- which in C are typically dynamically
allocated via looped pointers and offsets.  There is also a strong
motivation to allocate the actual physical memory as a single contiguous
vector (perhaps inside a struct) so it can be passed easily to
vector-based e.g. ode solvers.  Since such multidimensional matrices may
themselves be matrices of structs as easily as matrices of e.g. ints of
doubles, in the idealest of cases the constructors and destructors would
need to do their work with void typing and permit a cast to determine
what exactly is in the matrix and how large it is.

At that time I wrote such a matrix allocator in something like the GSL
style -- the style itself had to be slightly modified as the basic
objects in use pretty much presupposed 1-2d simple matrices.  It could
manage up to (IIRC, I haven't looked at it or used it for a couple of
years now) 9th rank tensors or thereabouts with more or less arbitrary
index ranges on each slot.  If there is any interest in this I'd be
happy to revive it and/or contribute it for consideration.

The point being that it sounds like doing lapack "right" in GSL may
involve messing some with the GSL matrix definition anyway, which in
turn may involve a "major" number bump -- routines that use matrices may
need to be rewritten to accomodate this.  If this happens, it would be a
really good time to rethink tensors altogether and implement SOME sort
of high-rank tensor in the GSL, ideally making matrices 2d object in a
smooth hierarchy.

The other thing that would need to be considered at that time, BTW, is
that whatever matrix schema is used, it needs to minimize redirection
and "object" like access in favor of efficiency.  One reason that people
still use Fortran for a lot of linear algebra-based arithmetic is that
its fairly rigorous matrix typing permits compiler-level optimizations
that object-like implementations of C and/or C++ do not.  This is the
reason that I altered the way matrices were packed in my code -- the
actual storage is always a single contiguous vector *data whose offsets
are repacked into the usual **...*matrix objects.  One lets the compiler
manage the offset arithmetic during access in matrix[i][j]...[n] form,
but one HAS the data in data[i] form under a single pointer, where a
cast permits either one to manage the right offsets.  I'm not certain
that this is the most efficient way of doing it, BTW -- it is just one
that makes it very easy to code and to be able to visualize e.g.
blocking and stride while one does so.

This is the kind of thing that it would be good for a Real Computer
Scientist (and C God) to look at and pronounce upon.  I'm at best a
possibly gifted amateur here for all of the C code I've written over the
years -- I haven't focussed on operational efficiency the way the
compiler guys have.  Maybe Greg Lindahl or somebody else from pathscale
or the gcc project could give us guidance on the most efficient ways of
implementing this, or we could even work hand in hand with the gcc folks
to get "blessed" optimization of whatever way tensor forms are
implemented in pointers so that we could eventually make the GSL as
efficient as a fortran library in this regard at least with gcc.

   rgb

--
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