Re: [boost] uBLAS formal review

From: Joerg Walter (jhr.walter_at_[hidden])
Date: 2002-06-29 08:04:14

----- Original Message -----
From: "Martin Weiser" <weiser_at_[hidden]>
To: <boost_at_[hidden]>
Sent: Friday, June 28, 2002 11:20 AM
Subject: Re: [boost] uBLAS formal review

> On Freitag, 28. Juni 2002 00:24, Joerg Walter wrote:
> > ----- Original Message -----
> > From: "Martin Weiser" <weiser_at_[hidden]>
> > > (i) Given the performance drop for large matrices, there should be a
> > > possibility to supply high performance, specialized algorithms for
> > > core algorithms such as _gemm. I tried to follow the instantiation
> > > and execution path for C.assign(prod(A,B)) and got the impression
> > > that it would be a hard job intercepting the expression template
> > > chain in order to provide a specialized algorithm (correct me if I'm
> > > wrong).
> >
> > We'd have to specialize vector_assign<> and matrix_assign<>.
>
> Ok, if you say it's possible I'll be satisfied. After all you know the
> code better than me ;)

Yep. But I shouldn't only state it, but also show. So I wrote a little
sample:

----------
#include <iostream>
#include <boost/numeric/ublas/vector.h>
#include <boost/numeric/ublas/matrix.h>
#include <boost/numeric/ublas/io.h>

namespace boost { namespace numerics {

    // vector assignment_operation scalar
    template<>
    struct vector_assign_scalar<class scalar_multiplies_assign<double,
double> > {
        void operator () (vector<double> &v, const double &t) {
            std::cout << "here we are" << std::endl;
            for (unsigned i = 0; i < v.size (); ++ i)
                v (i) *= t;
        }
    };

} }

int main () {
    numerics::vector<double> v (1);
    double d = 0;
    v.clear ();
    v *= d;
    return 0;
}

----------

which should show, how to achieve this (works for me under GCC 3.1). And
yes, I realize, that vector_assign<> and matrix_assign<> don't have the
optimal interface: they should have been free functions in fact, but this is
the way, we got the thing working on MSVC 6.0 ;-(.

> It's just that I considered specializing matrix_assign<> for dgemm and,
> reading the code, thought the necessary expression type information to be
> lost somewhere. Probably I'm just to dumb. Preferring a design that
> allows such specializations, I just wanted to draw your attention to this
> point.

Regards

Joerg