Re: [eigen] aliasing system

```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.

>
> 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.

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/