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

[ Thread Index | Date Index | More Archives ]

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 =;
>> Eigen automatically performs the solve operation inplace, just like it
>> would do with :
>> Advantages:
>> - simpler API: only one function
>> - works fine with writable temporary objects, e.g., 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 =;
C =;

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):
> 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 (
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

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 ?


> Is someone currently working on the 'Direct-access type unification'-TODO?
> (or on my issue from above's mail?)
> 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
    Reference(T& obj) : m_object(obj) { }
    operator T&() { return m_object; }
    T& object() { return m_object; }
    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)

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

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

  cout << b << endl;
  return 0;

Mail converted by MHonArc 2.6.19+