[boost] Re: Re: uBLAS and expression templates

From: Benedikt Weber (weber_at_[hidden])
Date: 2002-07-01 04:30:07

Joerg and Mathias,
Thanks you both for the fast response!

"Joerg Walter" <jhr.walter_at_[hidden]> wrote

> Ok. After discussing this issue intensely with Mathias (thanks, Mathias!),
> I've added overloads like the following:
>
>     template<class M, class E1, class E2>
>     NUMERICS_INLINE
>     M
>     prod (const matrix_expression<E1> &e1,
>           const vector_expression<E2> &e2) {
>         return M (prod (e1, e2));
>     }
>
> With these new overloads, we'd recommend to write your original sample in
> the following way:
>
>     t.restart ();
>     for (int run = 0; run < 1000; ++ run)
>         r.assign (prod (prod<matrix<double> > (A, B), v));
>     std::cout << t.elapsed () << std::endl;
>     t.restart ();
>     for (int run = 0; run < 1000; ++ run)
>         r.assign (prod (A, prod<vector<double> > (B, v)));
>     std::cout << t.elapsed () << std::endl;

This looks like a resonable pragmatic solution. Although still much writing,
it is more readable than
 r.assign (prod (A, vector<double> (prod (B, v))));

We have to see, if this covers all the cases. My first encounter to the
problem was actually when I timed
prod(s*A,B), with "s" being a scalar. I was really surprised, that it took
much longer than s*prod(A,B) until I found the reason. However, this case
can easily be avoided if one is aware of it, and it's not likely that people
write it anyway. A paragraph in the documentation would be helpful though.
There might be other cases too, which I am not aware of now.

This expression template business is quite exciting, especially for me being
exposed to it for the first time. I went through the introduction by Todd
Veldhuizen and tried a few things to better understand. My vision for a
matrix library is still that expression templates could take care of the
above optimization transparently to the user, and that they could also do
blocked operations. However, I see that this no easy task and I have to
study the subject for some time to see all the problems involved, especially
when combined with different matrix classes and iterators.

Benedikt