This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: Two suggestions...
- From: "Robert G. Brown" <rgb at phy dot duke dot edu>
- To: Gerard Jungman <jungman at lanl dot gov>
- Cc: GSL Discussion list <gsl-discuss at sources dot redhat dot com>
- Date: Wed, 11 Jun 2003 17:45:28 -0400 (EDT)
- Subject: 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