Re: [eigen]

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

wow ! that's a rather radical change killing half of the purpose of the ReturnByValue. Moreover, at least the partial piv LU and LLT solvers works in place, so:

x = m.selfadjointView<Lower>().solve(x);

and

x = m.lu().solve(x);

already worked fine, without any memory alloc. So enforcing the use of noalias here is well, a nonsense.

Moreover (bis), for dynamic sized matrices, m = m.inverse() use partial piv LU, and should works fine inplace.

Moreover (ter), while I agree that the solve functions are more like products, inverse() is more like transpose() for which aliasing has to be manually solved.

So overall, I'm not a big fan of this change. If you want we can do automatic aliasing detection for the fixed sized inverse (in debug mode only, reusing the mechanism of the transpose function). For the solvers, we could enable the EvalBeforeAssiging on a per solver basis.

What do you think?

gael

On Sun, May 30, 2010 at 7:46 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

This is now fixed..

The problem was not specific to inverse(), but existed for all ReturnByValue-enabled function. For example, doing x = y.solve(x) would probably have exhibited the same aliasing issues.

The fix was to make ReturnByValue evaluate-before-assigning. If you want the old behavior (slightly more efficient), you can use .noalias() to tell Eigen not to care about aliasing, like this:

  x.noalias() = y.inverse();  /// don't do this if y aliases x !!

Since this is the behavior for products, it felt natural, API-wise, to do the same for other "multiplicative" operations such as inverse(), and more generally, solve() and friends.

Benoit

2010/5/30 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>

aaahh... this had me puzzled for a while, but what you really have here is an aliasing problem. With the devel branch, when you do S = S.inverse(), S starts being overwritten with its inverse and then read again to compute the rest of the inverse... i need to have a look at the code to see if we're getting any speed gain from that, a priori i'd say let's work like in 2.0, taking care of aliasing, and if later there is demand for a super optimized noalias variant, we can then add it with an explicit name like a.noalias()  = a.inverse(), like for products.

Benoit

2010/5/28 SHIROKOBROD Oleg <Oleg.Shirokobrod@xxxxxxxxxx>

Hello,

 

I’ve met a problem with development branch. Here you are the case.

                                   { -2.0                                     -1.0 }

I have a 2x2 matrix S = {                                                   }

                                   { 0.00026130287408377900     -2.0 }

 

and I want to inverse it. This is the code I compiled with MSVS 2008

 

#include <Eigen\Core>

#include <Eigen\LU>

 

Eigen::Matrix<double, 2, 2> S;

S << -2.0 , -1.0 , 0.00026130287408377900 , -2.0;

S = S.inverse();

 

Eigen2 gives the right answer

 

{-0.499967339	0.24998367}
S.inverse  =  {
{-6.53215E-05	-0.499967339}

 

 

Eigen3 (development branch I got on 26.05) gives a wrong answer

 

 

{-0.499967339	0.24998367}
S.inverse  =  {
{-6.53215E-05	-0.124983670170524}

 

 

What’s wrong?

 

Thank you,

 

Oleg Shirokobrod

 

___________________________________________________________________________
This e-mail is confidential and is for the addressee only.   Please refer to
www.oxinst.com/email-statement for regulatory information.