Re: Tensor in GSL.

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

On Thu, 18 Dec 2003, Andrea Riciputi wrote:

> Hi,
> I'm wondering if tensors (i.e. arrays with more than 2 dimensions) will 
> be implemented in GSL or if someone is already working on this. I need 
> high dimensional arrays for my own code, and I was considering to 
> encode such a data structure. Perhaps it could be included in GSL when 
> (and if) I get it working.

I have been working on it (or rather, I worked on it a few months ago
and got it to where I have a functional setup).  It works well enough
that I've implemented it in a couple of very simple programs, but isn't
by any means fully tested or debugged and might even leak (although I
don't THINK it leaks and DID test it for leaks:-).

On the good side, the stuff I've written has the following (I think
desirable) features.  It supports up up to 8th rank tensors -- although
there may be a way in C to do arbitrary rank with a single void *
pointer or some sort of recursion, I didn't find it and so my code has a
very repetitive and scalable structure for each rank one at a time.  The
upper and lower indices are passed into the allocation routine as
vectors of length <= 8.  This permits one to efficiently allocate and
dereference e.g.  angular momentum (triangular) tensor forms with
indices that go like:

  ylm[0...lmax][-lmax...lmax]

or

  3j[0..l1max][-l1max..l1max][0..l2max][-l2max..l2max][0..l3max][-l3max..l3max]

Finally, the routines transparently allocate/manipulate tensors of
structs, strings, or arbitrary data types -- not just int, double, etc.
All allocation is of void type entities with a specified stride, all
actual typing in an application is accomplished by casting the void type
pointer to a corresponding typed pointer.

On the bad side, the code to accomplish this is not exactly trivial, and
getting all this to work smoothly caused some damage to the underlying
block, vector and matrix interfaces (the latter, in particular, are
completely obsoleted).  I have not yet tackled the fairly daunting task
of re-integrating the result with the rest of the GSL.  Things that rely
on only the actual underlying data block will probably port well enough,
and I can probably modify the tensor code if necessary to accomplish
this.  However, code that relies on higher rank views will almost
certainly break, and break pretty seriously.  The GSL group will have to
decide whether or not a consistent nth rank tensor interface is worth it
as a replacement, worth it as an alternative (renaming "matrix" to e.g.
"mtxptr" everywhere and rewriting code to preserve or modify the
existing raw data block struct only), or if there is a better way to do
it altogether.

At this point it is stalled primarily for want of time.  I'm really
heavily committed.  Some would say insanely overcommitted. I write a
column for Cluster World (and am a week past deadline as I sit here
typing:-), am "active" on the beowulf list, teach physics (and JUST got
my grades in), have several GPL projects I run or participate in, do
both computer science and physics research, and then there is my family
and my wife, the physician:-). 

One other concern: I'm a self-taught C programmer (granted that I taught
myself back in the early 80's and can program in more languages than I
can remember the names of so I'm an EXPERIENCED self-taught
programmer:-).  I'm not at all certain that the stuff I'm doing with the
void * typing and casting is the only or best way of accomplishing a
typed matrix addressable tensor with arbitrary length entities via
library calls although it seems like a slick trick and works and I
didn't have any luck trying various alternatives that MIGHT have worked.

The application interface is hence a tiny bit clunky.  Here is an
example for 4th rank tensor from my test routine:

 (Note that the tensor allocated has indices that run from -3 to 1, 0 to
 8, -7 to 5, -1 to 1.)

 size_t lower[] = {-3,0,-7,-1};
 size_t upper[] = {1,8,5,1}
 int *it;
 int ****t4;

 int i;
 int block_cnt;
 new_gsl_tensor *t;

 /*
  * 4th rank tensor/4 dim matrix
  */
 t = tensor_alloc(4,sizeof(int),0,0,lower,upper);
 /*
  * Extract pointer to its data block
  */
 it = (int *) t->block->data;
 /*
  * Determine the number of elements in the block
  */
 block_cnt = t->block->size/t->block->stride;
 /*
  * Initialize them with simple counter values.  The tensor has a very
  * straightforward "vector view" within the block struct, but it is
  * necessarily type void so it can be cast to arbitrary type.
  */
 for(i=0;i<block_cnt;i++){
   it[i] = i;
 }
 /*
  * Extract and CAST pointer to actual dereferencable, typed matrix.
  */
 t4 = (int ****) t->matrix->mptr;

t4[-3][0][-1][-1] should now equal 1, etc.  I include the various
block.h, matrix.h, and tensor.h includes below so that the damage that
might be done to the block/vector/matrix API can be assessed.

The code is all fully GPL and I'd be happy to put my development
directory with test program/example up on the web if you're interested,
or of course you can start over.

   rgb

(see includes below)

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

/*
* $Id: block.h,v 1.3 2003/06/24 15:52:51 rgb Exp $
*
* See copyright in copyright.h and the accompanying file COPYING
*
*/

/*
 * BLOCK
 * The block is the basic data object.  It is a wrapper for a pointer
 * to a block of memory and a small amount of attached header data.
 *
 * Note that the data block has void (no) type.  It is a raw pointer to
 * the address of a block of memory of size size (in bytes) and must thus 
 * be handled with care.  Presumably malloc/free can still manage.
 *
 * A gsl_vector was just a block with a stride and owner.  By putting
 * stride and owner into the block itself, tremendous simplification
 * is possible.  If stride=owner=NULL, obviously the block is just
 * a block.  If stride is non-NULL but owner is NULL, it is a typed
 * block.  If stride is NULL but the owner is non-NULL it is an owned
 * (but untyped) block (ownership is as likely to be useful for blocks
 * for anything else, I suppose).  If both are non-NULL, the block is
 * an (old-style) vector, but without the redundant block pointer.
 */
typedef struct
{
 size_t size;
 size_t stride;
 int owner;
 void *data;
} new_gsl_block;


/*
 * Keep it very, very simple.
 *  new_gsl_block_alloc allocates a basic block using malloc.
 *  new_gsl_block_calloc allocates AND zeros the block.
 *  free frees the block.
 */
new_gsl_block *block_alloc(size_t size, size_t stride, int owner);
new_gsl_block *block_calloc(size_t size, size_t stride, int owner);
int block_free(new_gsl_block *block);


/*
* $Id: matrix.h,v 1.2 2003/06/24 15:52:51 rgb Exp $
*
* See copyright in copyright.h and the accompanying file COPYING
*
*/

#define MATRIX_MAX_DIM 8

/*
 * MATRIX The matrix is a struct that contains the generalized, packed set
 * of pointers and parameters required to create the tensor itself and
 * pack a "matrix view" of it according to those parameters.  The matrix
 * does NOT contain a direct reference to the data block itself as there
 * may be a number of matrices packed from a single block -- for example,
 * in relativity a 2nd rank 4x4 tensor might have the embedded 3x3
 * submatrix for the space-like piece accessible by a different name.  A
 * matrix can also be allocated independent of a tensor (straight from an
 * already existing data block).  Freeing a matrix does NOT free the
 * underlying tensor data (unless done with a routine that does both).
 *
 * As noted above, the matrices themselves are of void type and hence
 * CANNOT be used to directly address the data.  They MUST be cast to a
 * pointer of the matching rank and appropriate type.  There is, however,
 * no reason not to make a bunch of shell subroutines to do the requisite
 * cast for you.
 */
typedef struct
{
  int dim;                /* Must be less than TENSOR_MAX_DIM */
  int signature;          /* 0=covariant; 1=contravariant */
  /*
   * This permits several ways of generating a matrix out of a block.
   * Note that all of these are basically vectors with indices that
   * run from 0-dim.  At least two of them (typically lower/upper
   * or lower/length) are used to generate the matrix view.
   */
  int *lower;
  int *upper;
  int *length;
  /*
   * the number of *'s here MUST be greater than or equal to
   * MATRIX_MAX_DIM.
   */
  void ********mptr;
} new_gsl_matrix;

/*
 * Very, very simple.
 *  matrix_alloc allocates a matrix struct and fills it in (a complicated process!)
 *  matrix_free frees it (but NOT the data block it references!)
 */
new_gsl_matrix *matrix_alloc(size_t dim, int signature, int *lower, int *upper);
int matrix_free(new_gsl_matrix *matrix);

/*
 * matrix_row_alloc allocates space for a matrix row of POINTERS.
 * The space is returned with void pointer type, and is presumed
 * to be filled only with memory addresses when the matrix is fully
 * packed.  As noted above, the entire matrix is cast to an appropriate
 * type before referencing it, so all pointers dereference correctly.
 */
void *matrix_row_alloc(size_t length, int lower_index);
int matrix_row_free(void *matrix_row);


/*
* $Id: tensor.h,v 1.8 2003/06/24 15:52:51 rgb Exp $
*
* See copyright in copyright.h and the accompanying file COPYING
*
*/

#include "block.h"
#include "matrix.h"


/*
 * The tensor struct.  A tensor is a block and matrix.  Nuttin' to it.
 */
typedef struct
{

 new_gsl_block *block;
 new_gsl_matrix *matrix;

} new_gsl_tensor;


/*
 * Very, very simple.
 *  alloc allocates everything (block and matrix).
 *  calloc allocates everything AND zeros the block (calls alloc).
 *  free frees everything.
 * Can't get much simpler.
 *  In almost ALL CASES once one has allocated the tensor one will want
 * extract and cast the actual packed matrix to a ***...**ptr of whatever
 * type you are trying to create.  This ptr is never freed, but the tensor
 * (or its underlying block and matrix) must be.
 *  Once this cast is done, the matrix/tensor can be referenced via the ptr
 * like ptr[i][j][k]...[m][n]
 * as long as i,j,k...m,n remain in bounds relative to the lower and upper
 * tensor index vectors passed in the tensor_alloc() call.
 */
new_gsl_tensor *tensor_alloc(size_t dim, size_t stride, int owner, int signature, int *lower, int *upper);
new_gsl_tensor *tensor_calloc(size_t dim, size_t stride, int owner, int signature, int *lower, int *upper);
int tensor_free(new_gsl_tensor *tensor);