[boost] Re: multidimensional arrays (was Boost Mathematicians)

From: David Abrahams (dave_at_[hidden])
Date: 2004-07-09 17:14:53

Matthias Schabel <boost_at_[hidden]> writes:

> I've been following this thread with some interest as I have had a
> number of similar complaints with the design and implementation of
> uBlas and have been frustrated for years by the lack of a
> comprehensive and flexible framework for handling multidimensional
> array data in a manner which enables the average end user to ignore
> the details of implementation of storage scheme, etc... while
> providing a consistent interface with sufficient comprehensiveness to
> handle as many reasonable use cases as possible.  It seems to me that
> a Boost matrix/linear algebra library could very profitably be built
> on top of a well designed policy-based multidimensional array class.

Maybe so, but I wouldn't start there.  The way any particular class
template is defined should be one of the last concerns in the design
of any generic framework.  It's all about concepts and algorithms.

> BTW, I am aware of the existence of Blitz++, which provides a number
> of interesting capabilities and is to be lauded for its demonstration
> that C++ multidimensional array performance need not be inferior to
> FORTRAN.  However, it seems like the support of that library is weak
> at present and there appears to be no effort to drive toward
> standardization.
>
> A couple of thoughts on such a D-dimensional array class :
>
> 1) it should implement as much of the STL one-dimensional container
> interface as possible 

If by that you mean driving towards something based on the STL
(reversible) container requirements in tables 65 and 66, then I
strongly disagree.  Those are bloated concepts and little or no
generic code has been written to take advantage of them.

> so that algorithms which are unconcerned with the order of traversal
> can use an efficient 1D iterator interface

That, I agree, is important.

>  - for example, if I want
> to replace every element of a matrix with its square, I should be
> able to do this :
>
> matrix_type mat;
>
> for (matrix_type::iterator it=mat.begin();it!=mat.end();++it) *it =
> (*it)*(*it);

That's potentially rather inefficient with all those .end()
invocations.  Anyway, I don't think low-level iteration over elements
should be required by the concepts to be exposed as the primary
begin() and end() iterators of a matrix, but through some kind of
indirection.  Since I am predisposed to using a free-function
interface:

   std::for_each(begin(elements(mat)), end(elements(mat)), _1 = _1*_1);

> The advantage here is that the algorithm doesn't need to know if the
> matrix is sparse or dense, as that logic is in the iterator, and
> applies equally well to higher rank arrays.  Indexed element access
> can just be viewed as a mapping of a D-dimensional index to the
> corresponding 1D offset.

There's precedent for that in the MTL already.

<snip other details of array interface>

I have no problem with any of these details being part of any specific
array class.  Whether or not they should be part of the library's
concept requirements is another question (on which I haven't formed
an opinion yet).

> If there is interest, I have written some starting code that I would
> be happy to share which addresses most of these points to some extent
> and provides some other nice properties such as allowing transposition
> of arbitrary indices with no need to reorder the underlying storage
> (that is mat.transpose(0,1) causes all subsequent calls to mat(i,j) to
> provide the transposed element without moving the elements
> themselves).  

I'm not sure that's a nice property, since it means that the
traversal order of the matrix is not a compile-time property.

> I have functioning dense, sparse, and diagonal
> implementations (with a proxy reference type) as well as indirection,
> subrange, and mask slicing.  The design is uses a single policy
> template for the underlying implementation :
>
> Array<T,R,Imp> - T is the value type of elements, R is the rank, and
> Imp contains the functional implementation details.
>
> I'll see about shaping the code up this weekend and posting it on the
> Yahoo files section if there's interest....

Your matrix implementation might be very nice, but from my point of
view, using any specific implementation as the basis of a new linear
algebra library design is starting from the wrong end of things... so
for me, at least, looking at it in detail would be premature.

-- 
Dave Abrahams
Boost Consulting
http://www.boost-consulting.com