Re: [ublas] Aliasing for += and -=

From: Ian McCulloch (ianmcc_at_[hidden])
Date: 2005-06-17 06:43:03

Gunter Winkler wrote:

> On Friday 17 June 2005 10:43, Michael Stevens wrote:
>> Guys,
>>
>> On Thursday 16 June 2005 21:00, Gunter Winkler wrote:
>> > Am Donnerstag, 16. Juni 2005 13:10 schrieb Toon Knapen:
>> > > IMHO the way to have both is as I have described below: we know that
>> > > in uBLAS different containers are never aliased. Aliasing only occurs
>> > > when one is working with views.
>> > >
>> > > So from the moment only containers
>> > > are used in an expression, the expression engine can safely assume
>> > > there is no aliasing.
>>
>> ??? I'm still not sure we are on the same wavelength with regard to
>> aliases. Hopefully the following is correct.
> 
> in order to give an example:
> 
> double alpha;
> Vector x;
> Matrix A;
> x += alpha * inner_prod( x, prod(A, x) );
> 
> Now ask yourself: Do we need a temporary or not? So it is really hard to
> find out wether an expression aliases the lvalue or not.

Well, that expression does not make sense (adding a scalar to a vector), but
anyway, itsn't the analysis rather straightforward?  The issue is whether
updating an element of the left hand side can affect subsequent
computations on the right hand side.  For the case of, eg

x += prod(A, x);

this is clearly "yes".  It should instead be x += copy(prod(A,x)); or
something similar.  Interestingly, this expression maps onto a single BLAS
call, without a temporary, (in the presence of aliasing the result is
unspecified [or undefined?], so you are also allowed to 'accidentally' get
the correct answer too ;) so the copy() function needs to be fairly smart.

For a different case,

Vector x,y;
Matrix A;

x += y * inner_prod(x, prod(A, x));

the result of inner_prod() is just a number, so the outer loop does not
involve any aliasing.

The difficult cases are where you do not know the origin of the arguments:

void foo(Matrix& M, vector& v1, vector const& v2)
{
   v1 += prod(M, v2);  // do we need a copy here?
}

but these problems always exist:

void foo(std::vector<double>& v, double const& x)
{
   for (int i = 0; i < v.size(); ++i) v[i] += x;  // potential aliasing!
}

Cheers,
Ian