Re: [boost] Re: slicing and ranges of matrices

From: Toon Knapen (toon.knapen_at_[hidden])
Date: 2002-01-18 14:28:21

Ralf W. Grosse-Kunstleve wrote:

> FWIW: a few theoretical remarks (born out of very pragmatic needs):
> 
> I came to the conclusion that multi_array will not
> immediately work for us (see my other messages). Therefore
> I thought some more about multi dimensional arrays. IMHO
> Toon's questions only show that a more general design is
> needed. There are three major concepts that I can
> distinguish, and which I believe should be presented and
> parameterized as such:
> 
> - Memory management
> - Access schemes
> - Iterators
> 
> Points:
> 
> Memory management:
>   Ultimately any array is a one-dimensional array. Here I
>   see six important memory management policies that are
>   useful in different contexts:
>     - automatic allocation, fixed size() (aka boost::array)
>     - automatic allocation (fixed capacity), variable size()
>     - dynamic allocation, fixed size()
>         deep copy semantics
>         reference counted
>     - dynamic allocation, variable size() (e.g. std::vector)
>         deep copy semantics
>         reference counted

is'nt the difference between e.g. 'automatic allocation, fixed size' and 
'dynamic allocation fixed size' not more a question of the way you 
instantiate it instead of being a policy ? (just to sync terminology to 
make sure that I know what you're saying)

>   The automatically allocated (i.e. stack) arrays will tend
>   to be small, therefore deepcopy semantics is always
>   appropriate.
>   For the dynamically allocated arrays, both deepcopy and
>   handle-based (referenced counting) copy semantics are
>   useful, depending on the application. However, an array
>   that is reference counted can easily provide a deep copy,
>   but not vice versa. Therefore in practice four memory
>   management varieties should be sufficient.
>   A comprehensive package should give access to each of
>   these.
> 
> Access schemes:
>   An access scheme is a transformation from a general index
>   object to a 1-dimensional array index. A simple index
>   object is a tuple of integers. A simple access scheme is
>   that of a conventional n-dimensional grid of points (as
>   represented by multi_array). In real-world applications
>   there are many other useful types of access schemes,
>   mostly related to some type of symmetry (e.g.  Hermitian
>   symmetry of complex-valued arrays).
>   A comprehensive package should parameterize the access
>   scheme (as a template argument) and provide the most
>   common schemes as part of the package.
>   Toon would provide a class implementing his custom
>   access scheme, and then instantiate, e.g.
>   multi_array<int, toon_accessor>.
>   
>   A slice is just a special access scheme, or a generalized
>   simple access scheme.

are'nt you mixing storage scheme and access schemes. After all, we want 
all matrices (sparse, skyline, dense, ....) to have the same interface 
(and thus 'access'-eable in the same way) but depending on the storage 
one way provides faster access compared to another access method (or 
complexity guarantees or dependent on the storage scheme underneath, 
which is actually similar to the STL)

thus an array library should be templated on its storage (dense, sparse, 
skyline, HB, ...) and on the other hand provide different access methods 
like subarray views (continuous, strided, indexed(non-strided)), 
iterators, ...

BTW talking about access methods. After the performance discussion about 
multi_array, I'm inclined to say that the generic access methods for a 
multi-dimensional array lib are probably not performant enough for 
numerics without the proper specialisations for the 1D and 2D case.

I think the linear-algebra world is big enough to justify this 
specialisation.

> 
> Iterators:
>   In the context of multi-dimensional arrays, iterators are
>   a somewhat artificial concept. An iterator that does not
>   keep track of the index object corresponding to the position
>   in the underlying 1-dimensional array (e.g. the mapping from
>   the 1-dimensional index to the corresponding three-tuple)
>   is either fast but of limited use, or general but slow.
>   A general multi-dimensional iterator must have some sort
>   of nested loop inside. Therefore an iterator and a nested
>   loop become almost indistinguishable.  However, I find it
>   much more natural to think of a nested loop rather than
>   an iterator. One could even argue that the STL iterators
>   are the trivial form of a nested loop, and the nested
>   loop is the proper abstraction for higher dimensions.

What I remember from looking at the MTL version 2 implementation is that 
  Jeremy succeeded in very intelligently implemting blocking based on 
iterators (or am I mistaking ?)

> 
> A partial implementation of these ideas can be found in the
> cvs tree at cctbx.sf.net (e.g. shared_storage.h, ndim.h,
> loops.h; maps/ndim.h implements a specialized accessor;
> naming is still a bit disorganized). We are currently
> refining the design. I would be happy to share more details
> if there is an interest.

I'm interested. Does this mean that you have implemented another 
matrix-library (no blitz, mtl, ublas) ?