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

From: Matthias Schabel (boost_at_[hidden])
Date: 2004-07-09 16:24:27

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.  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 so that
algorithms which are unconcerned with the order of traversal can use an  
efficient 1D iterator interface -
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);

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.

2) it should support two common syntaxes :

        mat[i][j] - zero offset C-style access
        mat(i,j) - possibly non-zero offset access

3) implementations for both statically allocated and dynamically  
allocated memory should be feasible.

4) for genericity (in algorithms which work for various values of D)  
there should also be an index type Index<D>
which may be passed to operator() as an index :

        mat(Index<2>(i,j)) - equivalent to mat(i,j)

5) should support various types of slicing : indirection, subranges,  
masking, etc...

6) should facilitate range checking if desired.

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 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....

Matthias

> Jeremy and I have just completed a re-evaluation of uBlas based on
> what's in uBlas' own CVS repository, having not discovered that until
> recently (you should have that info on uBlas' Boost page!)  We have
> some major complaints with uBlas' design.  The list is long, and the
> issues run deep enough that we don't believe that uBlas is a suitable
> foundation for the work we want to do.
>
> Here is a partial list of things we take issue with:
>
> Interface Design
> ----------------
>
> * Not grounded in Generic Programming.  The concept taxonomies, to the
>   extent they exist, are weak, and poorly/incorrectly documented.
>   Aspects of the design that should be generic are not (e.g. only
>   certain storage containers are supported, rather than supplying
>   storage concept requirements).  No linear algebra concepts (vector
>   space, field, etc.)  The library is so out-of-conformance with our
>   expectations for Generic Programming that this one item by itself
>   is probably enough to make it unsuitable for us.
>
> * Redundant specification of element type in matrix/vector storage.
>
> * size1 and size2 should be named num_rows and num_columns or
>   something memnonic
>
> * iterator1 and iterator2 should be named column_iterator and
>   row_iterator or something memnonic.
>
> * prod should be named operator*; this is a linear algebra library
>   after all.
>
> * begin() and end() should never violate O(1) complexity expectations.
>
> * insert(i,x) and erase(i) names used inconsistently with standard  
> library.
>
> * Matrix/Vector concept/class interfaces are way too "fat" and need to
>   be minimized (e.g. rbegin/rend *member* functions should be
>   eliminated).
>
> * The slice interface is wrong; stride should come last and be
>   optional; 2nd argument should be end and not size; then a separate
>   range interface could be eliminated.
>
> * No support for unorderd sparse formats -- it can't be made to fit
>   into the uBlas framework.
>
> Implementation
> --------------
>
> * Expressions that require temporaries are not supported by uBLAS
>   under release mode. They are supported under debug mode.  For
>   example, the following program compiles under debug mode, but not
>   under release mode.
>
>   #include <boost/numeric/ublas/matrix.hpp>
>   #include <boost/numeric/ublas/io.hpp>
>
>   int main () {
>       using namespace boost::numeric::ublas;
>       matrix<double> m (3, 3);
>       vector<double> v (3);
>       for (unsigned i = 0; i < std::min (m.size1 (), v.size ()); ++ i)  
> {
> 	  for (unsigned j = 0; j < m.size2 (); ++ j)
> 	      m (i, j) = 3 * i + j;
> 	  v (i) = i;
>       }
>       std::cout << prod (prod(v, m), m) << std::endl;
>   }
>
>   The workaround to make it compile under release mode is to
>   explicitly insert the creation of a temporary:
>
>       std::cout << prod (vector<double>(prod(v, m)), m) << std::endl;
>
>   There should be no such surprises when moving from debug to
>   release.  Debug mode should use expression templates, too, as the
>   differences can cause other surprises.
>
> * Should use iterator_adaptor. There is a ton of boilerplate iterator
>   code in the uBLAS that needs to be deleted.
>
> * Should use enable_if instead of CRTP to implement operators. uBLAS
>   avoids the ambiguity problem by only using operator* for
>   vector-scalar, matrix-scalar ops, but that's only a partial
>   solution.  Its expressions can't interact with objects from other
>   libraries (e.g. multi-array) because they require the intrusive CRTP
>   base class.
>
> * Certain operations, especially on sparse matrices and vectors, and
>   when dealing with matrix_row and matrix_column proxies have the
>   wrong complexity.  Functions such as begin() and end() are suppose
>   to be constant time. I see calls to find1 and find2, which look like
>   expensive functions (they each contain a loop).
>
> Testing
> -------
> * There should be a readme describing the organization of the tests.
>
> * Tests should not print stuff that must be manually inspected for
>   correctness.
>
> * Test programs should instead either complete successfully
>   (with exit code 0) or not (and print why it failed).
>
> * Abstraction penalty tests need to compare the library with tuned
>   fortran BLAS and ATLAS, not naive 'C' implementation.
>
> Documentation
> -------------
>
> * In really bad shape.  Redundant boilerplate tables make the head
>   spin rather than providing useful information.
>
> * Needs to be user-centered, not implementation centered.
>
> * Need a good set of tutorials.
>
> Compilers
> ---------
>
> * Simple example doesn't compile with gcc 3.3. Got it to compile by
>   #if'ing out operator value_type() in vector_expression.hpp.
>
>> On Saturday 05 June 2004 13:53, David Abrahams wrote:
>>>> What do you plan for MTL?  How is it different than ublas?
>>>
>>> MTL is aimed at linear algebra, whereas IIUC ublas is not.  
>>
>> Well the L and A in uBLAS certainly stand for Linear Algebra! Of  
>> course the B
>> stands for Basic and uBLAS's primary aim is to provide the standard  
>> set of
>> BLAS functions in a modern C++ environment. Of course as it stands the
>> complete uBLAS library is more then just the BLAS functions and  
>> includes some
>> common Linear Algebra algorithms and many useful types.
>
> Whoops.  Of course that's correct.
>
>> That said I think it is important to separate BLAS functions from
>> domain specific linear algebra algorithm development. This is
>> something that proved itself since the seventies.
>>
>>> There's a lot more to what's in the current plan than I can lay out
>>> here, but the focus will be on support for different kinds of  
>>> matrices
>>> with combinations of these aspects
>>>
>>>       o Shape
>>>       o Orientation
>>>       o Symmetry
>>>       o Sparsity
>>>       o Blocking
>>>       o Upper / Lower storage
>>>       o Unit diagonal
>>
>> Other then the many forms of blocking (other then banded) uBLAS  
>> supports all
>> these in its design.
>
> I believe that a design with really good support for blocking can't be
> easily grafted onto an existing design that doesn't have it.
>
>> This really is its strength! To a large extent they can even be
>> combine these properties where it makes mathematical sense. For
>> example you can wrap up one of a number of sparse matrix types in a
>> symmetric adaptor.
>
> This stuff was all present in MTL2 IIRC.
>
>>> and operations like:
>>>
>>>       o scalar-vector (vector-scalar) multiplication
>>>       o vector addition (and subtraction)
>>>       o apply linear operator (left)
>>>       o norm
>>>       o inner product
>>>       o triangular solve
>> Other then 'apply linear operator' these are all in uBLAS!
>>
>>> with expression templates, and of course, zero abstraction penalty  
>>> ;-)
>> Of course uBLAS does this all with ET, but the abstraction penalty  
>> may not be
>> zero :-)
>>
>> Other then the lack of ET in the current MTL the big difference
>> between the two libraries is the definition of iterators. Neither
>> design seems to be perfect with regard to efficiency.
>
> No, and I have some ideas for addressing that.
>
>> Since uBLAS is already in Boost and has a well established and clean  
>> user
>> syntax it would seem strange to ignore it.
>
> Yeah, I stopped ignoring it long enough to determine for sure that we
> should probably ignore it :(.
>
>> For the perspective of building further Linear Algebra algorithms it
>> would not be too hard to use the syntax sufficiently portably so
>> that a future MTL with expression templates could not be used
>> interchangeably.
>
> We have some problems with the syntax too, as you can see from the
> above.  That said, if the design of MTL makes sparing use of members
> and instead relies on free functions, you should be able to make
> uBlas syntax adapters ;-)
>
> --  
> Dave Abrahams
> Boost Consulting
> http://www.boost-consulting.com
>
> _______________________________________________
> Unsubscribe & other changes:  
> http://lists.boost.org/mailman/listinfo.cgi/boost
>
>
------------------------------------------------------------------------ 
---------------------------
Matthias Schabel, Ph.D.                           	
Utah Center for Advanced Imaging Research
729 Arapeen Drive
Salt Lake City, UT  84108
801-587-9413 (work)
801-585-3592 (fax)
801-706-5760 (cell)
801-484-0811 (home)
mschabel at ucair med utah edu