[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: Re: [eigen]
- From: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
- Date: Tue, 1 Jun 2010 08:38:30 -0400
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:received:in-reply-to :references:date:message-id:subject:from:to:content-type :content-transfer-encoding; bh=wqVGhnzvA3IDV0v9eVQRHYPun1vAi/ZTmZw1ExY3FkE=; b=yIUOmOU/T2kTae/vFdNz52P20ecadngvLA9hZCc/QdDS28Zn4UWS4JBdouEulYE5PR OwHrpdAeTgF8LltJOQSX6Btx+IBf/trLGSqwz2YDo7v/sUelOAC2CcmRVlov2XAbLhE/ vMspj2K/UMupcDuOgr8MWLtDRqC89LvYJBWGg=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type:content-transfer-encoding; b=RSY7gWYqTpfWMAohI5ZQUEYyo+qk5Xvfg9cdyg9YAMtsyU6eMQp1L0ecjbswXu/Yhv QQN90xdI43x1diuz+ShydBe8XKDLY/csRQfR6mntZwNVUhRUCd1c9Rc2Z8kuJ+P7WKgD 27DGMT9Pavi15ibKNdqDJ0EkWSM07O5ItQvcc=
2010/6/1 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.
Oops. I had forgotten about that!
> 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.
Oops (bis).
>
> 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.
Oops (ter). Good point! So, documentation issue.
>
> 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).
Ah yes, an assertion would help there.
> For the solvers, we could enable the EvalBeforeAssiging on a per solver basis.
OK. This will require some changes in ReturnByValue as, currently, it
doesn't care about the flags of Derived.
>
> What do you think?
I revert my change! Sorry!
Benoit
>
> 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.
>>
>