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

[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());

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