Re: [eigen] aliasing system

[Benoît Jacob <jacob@xxxxxxxxxxxxxxx>]

On Thursday 06 September 2007 16:02:15 Andre Krause wrote:
> so then if there is no way to avoid this problem - what about a sort of
> #pragma warning or assertion that gets thrown / printed in DEBUG mode
> only, if a user tries to do
>
> m = m*m;
>
> is there some template magic available to detect this?

Very good idea. Yes, I think that can be done with partial template 
specialization. Will try.

> it shouldnt not only be #1 in the FAQ, but #1 in the very beginning of
> an "about" or Introduction to Eigen2. You should explicitly write down
> all cases that could be dangerous.

writing down all cases is impossible as there are infinitely many. But we 
cansure write down the most common cases. And THE most common one by very far 
is a = a * b and a = b * a.

> for example, what about
>
> v = 2*v; // <-- doubling elements in a vector
>
> are all expressions problematic, where the same variable name (here v)
> is on both sides of the equation?

No, v = 2*v is not problematic, because eigen2 evaluates it as
for(...) v[i] = 2*v[i];
which is good.

The only reason why matrix-matrix product is problematic, is that the 
computation of (a * b)(i, j) involves not only a(i,j) and b(i,j) but also 
other entries which may have already been overwritten earlier in this 
matrix-matrix computation.

Here's an example with 2x2 matrices. Take m to be the following matrix:
(1 2)
(3 4)
And let's see what Eigen2 does when you do m = m*m.
It loops over all indices (i, j) and does
  m(i,j) = m(i,0) * m(0,j) + m(i,1) * m(1,j)
It starts with index (i=0,j=0). So it does:
m(0,0) = m(0,0) * m(0,0) + m(0,1) * m(1,0) = 1*1+2*3 = 7
So the matrix m is now:
(7 2)
(3 4)
Next, it does index (i=1,j=0). So it does:
m(1,0) = m(1,0)* m(0,0) + m(1,1) * m(1,0) = 3 * 7 + 4 * 3 = 33
which is incorrect! The problem was that it used the _new_ value of m(1,0), 
which is 7, whereas it should have use the _old_ value, which was 1.

Hope that helps,
Benoit

Attachment: signature.asc
Description: This is a digitally signed message part.