
Ublas : 
Subject: [ublas] Patch/proposal for overloading operator * in ublas
From: Jesse Perla (jesseperla_at_[hidden])
Date: 20090817 17:25:39
As discussed in previous posts, people really can't get past the lack of *
overloading in ublas. I know this was partially a philosophical/design
decision, but it is extremely important to "joe user". Can this be easily
remedied?
Here is a simple patch for discussion that implements:
matrix<double> A;
vector<double> y;
cout << 2.0 * y << endl;
cout << y * 2 << endl;
cout << A * 2 << endl;
cout << 2 * A << endl;
cout << A * y << endl;
cout << y * A << endl;
cout << A * A << endl;
The patch can be applied in the ublas directory. All it does is:
 For the scalar * vector and scalar * matrix operator overloads, it uses
boost::enable_if to remove anything from those overloads that doesn't fit
boost::is_arithmetic. This is to prevent a clash with the new operator *
 Then, new * operators are added for matrix*vector, vector*matrix, and
matrix*matrix, which forward directly to the appropriate prod() dispatcher.
What is missing:
 Is this as simple as I think it is? Is there any chance of overhead
here?
 What about y * y for two vectors? I think we should hold off on this
because of the inner_prod vs. outer_prod confusion. Perhaps with the
extended numeric_traits created for the bindings this could be used to
emulate some kind of matlab approach with explicitly stated orientation (row
vs. column vectors), but it may be more confusing than not. I personally
hate the conformance bugs that come out of * meaning both inner and outer
product in matlab.
 What about b *= A; and A *= B? Am I right to think this kind of
aliased multiplication should be using axpy_prod(b, A, b, false); and
axpy_prod(A, B, A, false); ? Should we just call this in the overload?
This operation comes up all over the place in the statistical work that i
use it for.
The big one, of course, is nested matrix multiplications.. I don't entirely
know how to do it, but here is my proposal:
 Operator * will do a (relatively) inefficient nested matrix
multiplication using a heap allocated, rowmajor matrix<T> as the temporary
object. What type to use for T? Use some kind of type promotion on the
value_type inside of the two argument types....
 The docs tell people what it does, and how to have finer control with
the prod(A, B, temp); approach.
 In line 4834, etc. of matrix_expression it tests for the complexity at
compile time. Instead, in the operator* overload, we could detect the
complexity at compile time and tagdispatch to a function that would create
the appropriate temporary object based on the rules I mentioned above. It
would then call the prod() dispatcher. We probably would have to do this
for the vector*matrix operator overloads as well.
 Would this work? Any overhead or gotchas? I don't think I am capable
of implementing it myself, but can help experts out with testing, etc?
Jesse