Matthias Schabel <boost@schabel-family.org> 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