Re: [eigen] Why no LU::solveInPlace()?

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2010/9/3 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> 2010/9/3 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>> On Sat, Aug 21, 2010 at 1:21 AM, Christoph Hertzberg
>> <chtz@xxxxxxxxxxxxxxxxxxxxxxxx> wrote:
>>> On 17.08.2010 09:36, Gael Guennebaud wrote:
>>>>
>>>> We definitely have to uniformize that, but actually we could also
>>>> consider removing all the solveInPlace functions...
>>>>
>>>> Indeed, when you write:
>>>>
>>>> X = A.lu().solve(X);
>>>>
>>>> Eigen automatically performs the solve operation inplace, just like it
>>>> would do with :
>>>>
>>>> A.lu().solveInPlace(X);
>>>>
>>>> Advantages:
>>>> - simpler API: only one function
>>>> - works fine with writable temporary objects, e.g., X.col(j) =
>>>> A.lu().solve(X.col(j));
>>>>
>>>> Drawbacks:
>>>> - need to write the "source" and "destination" object twice
>>>> - requires one runtime test (one pointer comparison which is disabled
>>>> at compile time if the types do not match).
>>>
>>> I really hope that this comparison will be optimized away in most cases by
>>> any recent compiler ...
>>
>> I tried the following:
>>
>> MatrixXf A, B, C;
>> B = A.lu().solve(B);
>> C = A.lu().solve(B);
>>
>> in the first case the if is optimized away, and in the second it is
>> not. But anyway here we are talking about a single if among O(N^2)
>> comparisons so do we really care ?
>> <SNIP>
>> Anyway, I would to have more opinion about whether there really is a
>> strong need for an additional solveInPlace(X) function when X =
>> solve(X) can already work in place ?
>
> For MatrixXf we don't. For Matrix2f we might (as the branching's cost
> can be non-negligible)... I really can't make up my mind. On the one
> hand it is unusual for us in Eigen to do this kind of if(), and it is
> unusual for us to tolerate even small / potential overheads when the
> matrices at hand may be fixed-size. On the other hand, since we don't
> currently have fixed-size specializations for these decompositions, we
> do currently have overhead in these decompositions. Yet again, this
> new if() would be a kind of overhead inherent to the API. Also, what
> if the B is an expression that takes longer to write. Then having to
> write it twice will be annoying. For example:
>
> B.foo().bar().long().expression() = A.solve(B.foo().bar().long().expression());

Wait, there are 2 more problems here that I didn't anticipate before.

The first problem is that  B.foo() will return a new (lightweight)
expression object everytime, so pointer comparisons will fail.

The second problem is that just taking a pointer to an expression
object could easily prevent the compiler from "optimizing away" this
object. I know we are already doing that in transpose() aliasing
detection but that is just in debug mode, right? Here we're now
talking about adding that also in non-debug mode and that could be a
bigger issue. Although I guess that the cost of constructing the
expression is totally negligible in this particular case (solving
being expensive enough) I thought that was worth bringing up. In the
case of complex enough expressions and small fixed sizes, this could
easily become non-negligible.

Benoit

>
> Finally, a solveInPlace() method is easy to look up and understand for
> people who just want to write code quickly and be sure that stuff is
> optimized. It's explicit.
>
> Maybe a good compromise is to keep solveInPlace() and also do the if()
> in solve(). Should make everyone happy: people who care about the cost
> of the if() will typically be able to use in-place solving anyway.
>
> Benoit
>
>>
>>>> Recall that to make the solveInPlace() function to work with writable
>>>> temporary objects we either have to declare the argument as a const
>>>> reference or enforce the user to name the temporary that can be a real
>>>> pain without the C++1x auto keyword.
>>>>
>>>> What other people think about it?
>>>
>>> I think the const correctness issues are somewhat related to something I
>>> complained about last month (2nd issue of that mail):
>>>
>>> http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2010/07/msg00095.html
>>>
>>> I think the clean way would be to pass some wrapper as const reference,
>>> which internally holds a non-const reference/pointer.
>>> Actually the only objects which can be passed to solveInPlace are
>>> direct-access types (unless I'm wrong) so technically you only need to pass
>>> a pointer to the start of the data together with the sizes and strides
>>> (where the latter can be templated).
>>
>> yes this is similar to what I propose a long time ago when we did not
>> have the ReturnByValue mechanism and having in/out arguments was
>> mandatory (http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2009/08/msg00073.html).
>> But since our functions are templated, you are enforced to explicitely
>> request the conversion to the wrapper type. See the attached example:
>>
>> MatrixXf m;
>> solveInPlace(m); // ok
>> solveInPlace(m.ref()); // ok
>> solveInPlace(m.col(i).ref()); // ok
>> solveInPlace(m.col(i)); // does not compile
>>
>> hm.... maybe a better solution would be to have a
>> CanBePassedByValueBase class that would be inherited by all
>> lightweight expressions, and then we would have two overloads for
>> solveInPlace:
>>
>> template<typename Derived> inline void solveInPlace(MatrixBase<Derived>& mat);
>> template<typename Derived> inline void
>> solveInPlace(CanBePassedByValueBase<Derived> mat);
>>
>>
>> gael.
>>
>>>
>>> Is someone currently working on the 'Direct-access type unification'-TODO?
>>> (or on my issue from above's mail?)
>>> http://eigen.tuxfamily.org/index.php?title=Todo_for_3.0
>>>
>>> I'd volunteer to work out some suggestions/concepts, but currently I am
>>> quite busy and you couldn't expect something before mid-September from me :(
>>>
>>> Christoph
>>>
>>>
>>> --
>>> ----------------------------------------------
>>> Dipl.-Inf. Christoph Hertzberg
>>> Cartesium 0.051
>>> Universität Bremen
>>> Enrique-Schmidt-Straße 5
>>> 28359 Bremen
>>>
>>> Tel: (+49) 421-218-64252
>>> ----------------------------------------------
>>>
>>>
>>>
>>
>