Hello

If I do this simple Matrix Computation . I should exipit aliasing, but why do I not? Because G.block(0,0,10,12)*P is a temporary?, Because of vectorization? I turned it of with EIGEN_DONT_VECTORIZE before header include....

Code: Select all

   Matrix<double, 12,12> G;
   G.setConstant(1);
   Matrix<double,12,1> P;
   P.setConstant(3);

   // here we should get aliasing because we update sequential
   P.segment<10>(0).noalias() = (P.segment<10>(0) + G.block(0,0,10,12)*P);
   cout <<"With Aliasing:" << P <<endl;

   P.setConstant(3);
   P.segment<10>(0).noalias() = (P.segment<10>(0) + G.block(0,0,10,12)*P).eval();
   cout <<"With out Aliasing:" << P <<endl;


Output

Code: Select all

With Aliasing:
39
39
39
39
39
39
39
39
39
39
3
3
With out Aliasing:
39
39
39
39
39
39
39
39
39
39
3
3

The other thing is this statement: (P_back, P_front, m_c are Vectors, G, R are Matrices).
How do I know where in this statement are temporaries involved?
Nowhere because I write .noalias(), and everything gets evaluated into P_front, so no temporaries for G*P_back...?

Code: Select all

P_front.noalias() = P_back - R*(G *P_back + m_c);


Thanks!!

this is because in your example:

P.segment<10>(0).noalias() = (P.segment<10>(0) + G.block(0,0,10,12)*P);

the nested product G.block(0,0,10,12)*P has to be evaluated into a temporary. To avoid it you should use operator+=:

P.segment<10>(0).noalias() += G.block(0,0,10,12)*P;

Then you get a wrong result.

In your second example you have two temporaries, one for each products. Again, you can avoid one by writing:

P_front = P_back;
P_front.noalias() -= R*(G *P_back + m_c);

In the near future, Eigen will be able to this kind of high level optimizations for you.

there are some hints there: http://eigen.tuxfamily.org/dox-devel/To ... ssion.html