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

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

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

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 ?

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

inline Reference<DenseBase> ref() { return Reference<DenseBase>(*this); }

#define EIGEN_DENSEBASE_PLUGIN "inoutplugin.h"

#include <typeinfo>
#include <iostream>

class ReferenceBase {};
template<typename T> class Reference;


#include <Eigen/Eigen>

using namespace Eigen;
using namespace std;

template<typename T> class Reference
{
  public:
    Reference(T& obj) : m_object(obj) { }
    operator T&() { return m_object; }
    T& object() { return m_object; }
  protected:
    T& m_object;
};

template<typename Derived>
void foo(Reference<DenseBase<Derived> > _bar)
{
  Derived& bar(_bar.object().derived());
  bar.x() -= 1;
}

template<typename Derived>
void foo(DenseBase<Derived>& bar)
{
  foo(bar.ref());
}

int main()
{
  int x = 2;
  RowVector4i b(1,2,3,4);

  foo(b);
  foo(b.ref());
  foo(b.head<2>().ref());
//   foo(b.head<2>()); // does not compile

  cout << b << endl;
  
  return 0;
}