[boost] Re: Boost Mathematicians

From: David Abrahams (dave_at_[hidden])
Date: 2004-07-05 15:39:10

Michael Stevens <Michael.Stevens_at_[hidden]> writes:

> First of all I am currently making sure Boost is updated with the latest uBLAS 
> development version. This already has updated documentation which points to 
> all the correct locations for discussions, CVS and Wiki. This should make it 
> much easier for users to find the correct locations if the stable version in 
> Boost does not meet their needs.

Thanks for that.

> To facilitate the development of solutions I will place many of the points 
> into the uBLAS bug tracker on SF. 
>
>> 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.
>
> This is a big issue. I believe we need to split the documentation
> issues from the design issues. 

I don't think you can really do that.  The documentation is
reflective of the thinking that went into the design.

> Documentation is high on the priority list of many uBLAS users and
> developers. Most new users start with a requirement for a C++ matrix
> and BLAS library with nice syntax, efficiency and documentation! The
> need for clearly defined and documented concepts comes for most
> people a little latter when extensibility becomes an issue.

When evaluating a library, I look to the documentation to decide
whether it's well-thought-out or it's going to cause me more
aggravation than it saves.  Documentation that doesn't clearly
describe what I must do to use and/or can count on from a library is
a show-stopper.

> The storage concept is not as bad as it appears. The concept
> although not documented is sufficiently well defined 

defined where?

> that creating a new model is unproblematic. At present the concepts
> are documented by their models :-)

I guess from the :-) you are joking.

>> No linear algebra concepts (vector space, field, etc.)
>
> I not sure if these concepts really need to be defined. They are
> already well described in many mathematical text books and in
> Wikipedia. 

Yes, but the issue is how they map onto expressions in the programming
language.  How do you multiply two elements of a Ring?  Do you use
"prod" or operator*?

> Almost all the major reference works on BLAS and Matrices
> simply state elements are either Real or Complex and are therefore
> implicitly a field.  

We're not just concerned about the elements of matrices, but the
vectors and matrices themselves.  This becomes a real issue when you
have to write algorithms that might operate on matrices of reals or
matrices of matrices -- yes, these things do come up in real
applications.

> Of course a C++ linear algebra library should be generic enough to
> work with any type which is a numeric approximation to a field.

"Be generic" in the sense of being founded on well-defined concepts,
yeah ;-)

>> * Redundant specification of element type in matrix/vector storage.
>
> Not sure what is meant here, can you elaborate ?

See my reply to Toon.

>> * 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.
>
> I think many uBLAS users have come to appreciate the brevity of
> naming. I think the concensus in the Boost review was as long as
> they were used consistenty these names are acceptable.

An aesthetic difference.  They should all be non-members anyway so
you can inject your own brief names.

>> * prod should be named operator*; this is a linear algebra library
>>   after all.
>
> Oh, there has been much ink spent on this! Joerg argues against
> operator* for two reasons I think. a) operator* brings an
> expectation of commutativity and this is not the case in linear
> algebra 

Then I suppose you also have an expectation that

     (A)(B)

can always be rewritten as

     (B)(A)

?

> b) The are many products (inner,outer,element) and there
> should be a clarity of naming.

If element product is what I think it is, it should use operator*
also.  I don't believe outer matrix product belongs in a linear
algebra library, but if you are going to do it you should use some
other function name.  Nobody has any confusion over the meaning of

      k j x A B

where k and j are real, x is a vector, and A, B are matrices when they
see it in mathematics.  In fact the uBlas overview uses this exact
notation.  There should be no problem replacing whitespace by *.  More
importantly, that multiplication is the *same* operation in many
cases: the set of matrices forms a Ring in the same way that the set
of reals does, so there should be a generic operation for multiplying
the elements of rings.

>> * Matrix/Vector concept/class interfaces are way too "fat" and need to
>>   be minimized (e.g. rbegin/rend *member* functions should be
>>   eliminated).
>
> Interesting. Certainly the fatness is a worry. I think rbegin/rend
> need to stay as members as this is part of the Reversible Container
> concept.

Whose Reversible Container concept?  If uBlas has one it should be
changed to use nonmembers, to minimize the amount of redundant
interface required to satisfy it.

BTW does anything in uBlas actually _use_ that concept?  If not, maybe
the concept should be eliminated.  I certainly didn't notice it
mentioned anywhere as a requirement.

>> * 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.
>
> For std::slice (valarray) compatibility stride should be last. Not sure why 
> this was not picked up on a long time ago!!
> valarray (and uBLAS) use size (not end) for good reason however.

...which is...?

> I also think it would be a mistake to drop range. slices can be
> degenerate (zero stride) and this has a runtime cost which range
> does not need to share.

OK.

>> 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.
>> 
>>       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.
>
> This was changed by Joerg after the last stable release to Boost. I will check 
> there has been no regression.
> The fix propagates an expression complexity that is statically checked. This 
> check explictly prevents the creation of nested products!! 

That's really going in the wrong direction, though.  There's no
reason the user should be forced to intervene here.

> Therefore the correct code is now that with the explicit temporary
> or with one of the parameterized prod's.  The change is aimed at the
> user. They must clearly express the evaluation and storage
> complexity of the expressions.

How is that "aimed at the user" (like a weapon, maybe? ;->)

> The change of semantics between debug/release builds is
> unfortunate. The partial disabling of the expression templates in
> debug mode is currently the default and a config macro. With current
> compiler technology this prevents debuggers crashing and users
> scrolling through type names that are pages long. In a perfect
> environment the default would be the opposite.

I think we can do better.

> Interesting solution. Agreed, the current implementation is closed to other 
> objects because of the intrusive base class. Sadly uBLAS was designed before 
> enable_if. Any pointers to compile/runtime effects of the enable_if solution?

I can't really imagine what you're asking for here.

> We would however lose the static type checking that argument
> matching with the CRTP base class gives us. At this point I think we
> land in a deep philosophical argument as to the strengths and
> weaknesses of generic progreamming %-)

I don't know what static type checking you think you'd lose.  It's
certainly not evident to me.

>> * Abstraction penalty tests need to compare the library with tuned
>>   fortran BLAS and ATLAS, not naive 'C' implementation.
>
> There was some very good work on this about a year ago. Sadly the
> results are no longer accessable on the web. Since the loss to Intel
> of the KAI compiler it is hard to see any library archiving really
> exceptional results in this regard.

I know :(.  Some people still have KAI.

>> 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.
>
> Strange gcc-3.3.3 is my primary compiler. Can you specify the problem more 
> exactly. Which gcc-3.3 subversion is causing the problem?

Jeremy?

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