From: Gerard Jungman <jungman@cybermesa.com>
To: gsl-discuss@sourceware.org
Subject: containers tentative design summary
Date: Mon, 05 Oct 2009 10:12:00 -0000 [thread overview]
Message-ID: <1254708349.18519.4.camel@ForbiddenPlanet> (raw)
[-- Attachment #1: Type: text/plain, Size: 184 bytes --]
A summary of the tentative design for containers. Code will follow,
as soon as a few things are sorted out. (See previous post on
questions about block/vector/matrix).
--
G. Jungman
[-- Attachment #2: gsl-containers-design-summary.txt --]
[-- Type: text/plain, Size: 11270 bytes --]
** Multi-array.
A multi-array is a data structure managing a reference to
a data segment (think double * as an origin) and enabling
indexed access to the data. Indexing is controlled by
an indexing map, which translates indices into an
offset in the underlying data segment. The number
of indices (dimension of the domain of the indexing map),
is called the rank of the indexing map.
The data segment may be externally controlled, or internally
controlled by the multi-array. In essence, "control" means
control of the allocation/deallocation process for the segment.
No notion of view is needed. When the data is internally
controlled, we can think of the multi-array as "original".
When the data is externally controlled, we can think of
the multi-array as a view. But this is not a precise way
of thinking. More precisely, the life-cycle of the data
is simply independent of the indexing semantics, and the
need for a separate notion of view evaporates.
This lack of a separate view type is not shared by some other
designs. For example, boost introduces a distinct reference
type boost::multi_array_ref, in addition to boost:multi_array.
The reason for this choice is not clear. It may be for absolutely
optimal efficiency. Certain choices are easier to make if one
separates these types; however, my current feeling is that
the extra type is not justified.
If we decide to introduce separate view types later, we must
insure that it is done in a way which does not force other
interfaces in the library to grow combinatorially. And
without doubt, the extra types must have the same generic
interfaces as the non-view types; this is one area where
GSL fails miserably.
** Indexing map.
An indexing map translates an index vector to an offset. It is
specified by a shape. At a minimum, shapes include the dimension
of each index (number of indices addressed). When a shape allows
an index in place i to have only one value, that index place
is called "invariate". If it can have no values (null set), it is
called "null". Otherwise (the normal case) it is called "variate".
The number of indices (dimension of the domain of the indexing map),
is called the "static rank" or simply "rank" of the map. The
number of variate indices in a given map is called the
"dynamic rank" or "effective rank" of the map.
Examples using C-array notation:
double x[3][4][5]; // static rank = 3, dynamic rank = 3
double x[3][1][5]; // static rank = 3, dynamic rank = 2
// [3] and [5] are variate, [1] is invariate
** Structure of indexing maps.
We choose the most general affine form for indexing maps.
Therefore, the offset is calculated as
offset = sum(n, index[n] * stride[n])
The data segment is stored as a pointer to the origin of
the data. Therefore, the most general map calculation
takes the affine form
location = origin + sum(n, index[n] * stride[n]).
** Slicing as composition of indexing maps.
Let i[N] be an element of the index set for a rank N indexing map.
Let j[M] be an element of the index set for a rank M indexing map.
Let delta[N][M] be an integer-valued matrix and b[N] an
integer-valued vector. Then the following formula which
maps M-dimensional indices to N-dimensional indices
defines a rank M index map from any rank-N index map.
i[n] = b[n] + sum(m, delta[n][m] * j[m])
Notice the functoriality. "Old" indices are determined in
terms of "new" indices. In other language, the new indices
are pulled-back to the old context. This is the correct
functoriality to allow arbitrary chained composition of
index maps, and therefore slices of slices, to any order.
It should be clear that the index map generated this way
is always affine, and all affine indexing maps can be
generated in this way.
Therefore, a general slice of an affine indexing map
is specified by a vector b[] and matrix delta[][], of
the indicated dimensions.
The reader can convince himself that any form of strided
access is supported in this way. For example, it is equally
easy to slice out the super-diagonal of a matrix as it is to
slice out a row or column, or for that matter to slice out
any regularly strided sub-matrix, etc.
The underlying direct data access, by computation of an
offset, is itself an affine indexing map, specialized to
the case where the "old" index is one-dimensional. So
the design space of affine indexing maps is closed.
** Multi-array implementation requirements.
* Indexing calculations must be as fast as possible. The indexing
arithmetic should be inlined when possible. The shape data
structure should enable the fastest possible computation,
consistent with general strided access.
This suggests that different rank objects be of different static
types, since this is the most straight-forward way to insure that
the correct indexing calculation is chosen at compile-time. This
design corresponds to the usual designs (such as in GSL), where
vectors and matrices (for example) are different static types.
This choice also corresponds to the design of boost::multi_array,
where the rank is a template parameter and therefore a
compile-time constant which generates a separate type.
* The shape of a multi-array should be dynamic. This corresponds
to the usual design, such as in GSL, where the dimensions of
a matrix (for example) are run-time constants, not compile-time
constants.
* It should be possible to construct multi-arrays without
heap allocations. For example:
double data[64];
gsl_marray_2 ma_8x8 = gsl_marray_2_build(data, 8, 8);
Placing such objects on the stack has good implications
for efficiency and is very natural from the client's point
of view. For example, in C++, one simply writes
std::vector<double> v;
which places the object on the stack as an automatic variable.
** Multi-array implementation choices. These are open for debate,
but these are the choices I have made. No loss of functional
generality is implied by these choices, only some loss of
syntactic generality.
* Multi-array indexing conforms precisely to C-style array
indexing. Shapes are completely specified by dimensions and
strides. This means that indices always "start" at 0, and
end at dimension-1 (inclusive).
* The notion of storage order does not enter. The "fast" index
is always to the right, conforming to C array semantics. See
the discussion of matrices below.
** Vector.
A vector manages a data segment and supports indexed access
to the data through a rank one indexing map. However,
vectors may be required to satisfy extra constraints, not
compatible with rank one multi-arrays. Also, they may provide
other capability, such as indexing which "starts" at some
non-zero value. Personally, I have no desire for this extra
sort of functionality. It would only add overhead to an
otherwise clean design. My first inlination is to
implement vectors as follows:
typedef gsl_marray_1 gsl_vector; // DONE
However, the reality may not be quite as simple as this in
the end. However it is done, I think it is a good idea
to preserve the name "vector" in the design, although
it will probably imply some significant boiler-plate
in function declarations, since C does not allow us
to simply inherit the member function namespace
of another class. SIGH.
** Matrix.
Matrices are clearly different from rank two multi-arrays.
There are necessary extra capabilities, and there are also
constraints, in order to make them conform to the notions
used in linear algebra packages.
The extra capability corresponds to the identification of "row"
and "column" indices. Multi-array semantics are governed by data
locality, so the notion of "fast" index is immutable. But a matrix
may either identify its rows with fast index traversal (C order)
or its columns with fast index traversal (Fortran order). Matrices
will need to remember which and to adjust their indexing computations
accordingly.
The extra constraint on matrices, in comparison to multi-arrays,
is the need for the fast index to traverse the data with stride 1.
For example, BLAS interfaces allow for vectors of arbitrary
stride, but do not allow for a matrix "fast index" stride.
Example:
Let A be a 3x3 matrix. Let B1 be a 2x2 matrix defined as a
view (sub-matrix) of A defined by the x's in the following
diagram:
x x -
- - -
x x -
Then B1 is a matrix, in the BLAS sense. It can be passed
to BLAS functions without further reshaping. To be exact,
it would be passed with origin equal to the given origin
and with tda = 6 and dimensions (2,2).
Let B2 be the 2x2 sub-matrix view defined by the following
diagram.
x - x
- - -
x - x
Then B2 is _not_ a BLAS compatible matrix, since there is no
parameter in the BLAS interfaces which allows for the non-unit
stride in the fast index. Some reshaping, involving data
motion/copy, would be required to pass this view to BLAS
functions.
Therefore, I do not think that generic rank two multi-arrays can
be identified with matrices. Rather, a lightweight conversion
from rank-two multi-arrays to matrices must be provided, and this
conversion must fail at run-time if the multi-array fast index
stride is not unity.
** Tensor
No explicit notion of tensor currently exists, beyond that
expressed in the tensor/ extension package for GSL, by
Jordi Burguet-Castell.
The Burguet-Castell package does not meet the requirements
layed out above, mainly because it looks and works too much
like the GSL types, having been modelled on them for obvious
reasons. In particular, it allows the rank to be determined
at run-time. This may or may not be desirable, depending on
how we want tensors to interact withy multi-arrays.
However, we draw some motivation from the Burguet-Castell package.
A tensor is defined to be a multi-array where the shape is constrained
to have all dimensions equal.
Like in the case of matrices, we encounter some small technical
problems in the identification with multi-arrays. Since the
shape of a multi-array is dynamic, the compatibility of
a given multi-array with the tensor concept is question
that must be deferred to run-time.
The question of the role of multi-arrays in the implementation
of tensors remains open, as it does for matrices. But we do
require a certain level of interoperability, perhaps expressed
as light-weight conversion functions which fail to convert a
multi-array to a tensor if the shape is not consistent. This
seems like the most likely solution. It also allows the tensor
implementation to be independent of the high-level multi-array
interfaces, which strikes me as the right way to do it, from
both a general design point of view and an efficiency point of
view.
next reply other threads:[~2009-10-05 10:12 UTC|newest]
Thread overview: 30+ messages / expand[flat|nested] mbox.gz Atom feed top
2009-10-05 10:12 Gerard Jungman [this message]
2009-10-05 14:50 ` James Bergstra
2009-10-05 23:00 ` Gerard Jungman
2009-10-05 23:45 ` James Bergstra
2009-10-06 19:59 ` Gerard Jungman
[not found] ` <645d17210910060537s762d6323pfd2bec8590ad28e9@mail.gmail.com>
2009-10-06 20:02 ` Gerard Jungman
2009-10-23 21:28 ` Brian Gough
2009-10-27 23:06 ` Gerard Jungman
[not found] ` <7f1eaee30910271628h70785125m68e47c7a7b5c25b7@mail.gmail.com>
2009-10-27 23:49 ` Gerard Jungman
2009-10-29 18:06 ` Brian Gough
2009-10-29 20:41 ` Gerard Jungman
2009-10-29 21:40 ` James Bergstra
2009-10-30 16:54 ` Brian Gough
2009-10-30 16:54 ` Brian Gough
2009-10-05 23:04 ` some general questions Gerard Jungman
2009-10-06 16:21 ` Brian Gough
2009-10-06 20:32 ` Gerard Jungman
2009-10-10 14:45 ` Tuomo Keskitalo
2009-10-24 11:16 ` Brian Gough
2009-10-14 20:00 ` Brian Gough
2009-10-15 18:39 ` Tuomo Keskitalo
2009-11-03 19:44 containers tentative design summary Gerard Jungman
2009-11-09 20:41 ` Brian Gough
2009-11-09 23:06 ` Gerard Jungman
2009-11-14 15:25 ` Brian Gough
2009-11-15 9:13 ` Tuomo Keskitalo
2009-11-15 16:44 ` Jonathan Underwood
2009-11-15 18:41 ` Robert G. Brown
2009-11-16 11:56 ` Brian Gough
2009-11-20 9:37 Justin Lenzo
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=1254708349.18519.4.camel@ForbiddenPlanet \
--to=jungman@cybermesa.com \
--cc=gsl-discuss@sourceware.org \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).