This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: Two suggestions...


On 11 Jun 2003, Gerard Jungman wrote:

> I think his main point (please correct me if I'm missing it)
> is that this way you have a simple "view" scheme for switching
> between vectors and matrices. The example of calling the ODE solver
> with an underlying matrix argument is a good one, and it highlights
> the issue.

Well, a vector view works for arbitrary tensors as well, as long as the
original data object is stored at a single *vector pointer.  So the ODE
solver is no problem -- linear algebra is (or can be).

The trick is in the embedded tensors of lower rank and the problem
persists even for vectors and (2D) matrices, really.  A row,column
ordered (Unix/C) matrix can easily be represented as separate columns,
because it is easy to extract the pointers to the column heads
(otherwise sequentially addressable).  However, a view of rows as a
"vector" is not so simple, as each row is a single element from each
column, separated by a large stride from one another in the parent data
*vector object.

A tensor of rank N thus has a single coordinate ordering in which rank
N-1 tensors can be found (recursively down to rank 1 vectors), and lots
of coordinate orderings in which the entire data object must be used,
with appropriate stride, in order to construct a view.  One is
computationally efficient (pass a single pointer cast to the appropriate
object time) and one is very INefficient, effectly getting the
components one at a time with indexing arithmetic.

So one can get >>an<< ordered-coordinate view of a tensor that can be
passed to BLAS etc, but not every view can be passed without pain or an
intermediate reordered tensor or view being constructed.  I think.

> Of course a tensor mechanism is more general than this, but
> if the main motivation is to have a uniform interface available
> for vectors and matrices, then it is an important idea.
> 
> We should consider doing something about this. I think there are
> two choices:
> 
>   1) create some hack which allows switching objects between
>      vectors and matrices, as a kind of special case
> 
>   2) create a tensor interface, as Robert suggests
> 
> I don't see any reason to prefer the first over the second, so why
> not consider tensors? Anyway, I don't think they would create any
> problems or require any changes to the current interfaces.
> Probably requires some fiddling under the hood if we want to
> do it right though...

Or accept that only certain tensor mappings are ever going to make
numerical sense for purely computational/toplogigical reasons.  A tensor
as a *vector is always ok, as long as it is allocated all at once so it
is a single contiguous data object.  Index-ordered views are ok, as
unique pointers to embedded tensors of lower rank exist in index order.
However, non-index-order views (or repacked representations in
***pointers) are possible only via a routine that basically uses the
*vector view with appropriate index resolution logic and will be very
difficult to pass to e.g. BLAS.

> Also, I thought the comment about vectors of random numbers
> was probably correct; it might be a bit tedious to add a
> vector version of the rng functions, but it seems like
> that's the least we can do, given all the good work that
> has already gone into that library component.

In time I may try to tackle this one myself as well, at least to see how
hard it is to implement.  However not for weeks to months.

   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




Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]