Hi all,

I'm using boost::proto to evaluate expressions involving Eigen matrices, in the context of a Finite Element Method (FEM) solver. The goal is to be able to write high-level math expressions over a FEM mesh, that are then evaluated using boost::proto.

However, whenever I have an expression involving more than 2 Eigen matrices, I run into trouble because temporary variables cease to exist. The following (stand-alone) sample mimics what happens inside proto:

Code: Select all

#include <iostream>

#include <Eigen/Dense>

using namespace Eigen;

/// Multiply left and right
template<typename T1, typename T2>
  typename ProductReturnType<T1, T2 >::Type
multb(const T1& left, const T2& right)
{
  return (left * right);
}

/// multiply left with a nested call to multb(right, left)
template<typename T1, typename T2>
  typename ProductReturnType
  <
    T1,
    typename ProductReturnType<T2, T1>::Type
  >::Type
multa(const T1& left, const T2& right)
{
  return left * multb(right, left);
}

int main(void)
{
  Matrix<double, 1, 2> A(2., 2.);
  Matrix<double, 2, 1> B(2., 2.);
 
  std::cout << multa(A, B) << std::endl; // should be 16 16
 
  return 0;
};


Basically, whenever there is a nested function call returning an expression, the temporary that is created gets destroyed before it is used in the evaluation. In the example above, multb produces a CoeffBasedProduct expression, and multa takes this by reference. Since the result of multb gets destoyed as soon as multb returns, the evaluation at the cout call uses a dangling reference.

In debug mode, this results in a segfault, while in optimized mode everything works, presumably because the function calls are inlined. The same thing happens with boost::proto, which basicly traverses the proto expressions using ultimately a set of nested calls to functors. I also run into the problem with GeneralProduct expressions, when there is a reference to an intermediate temp matrix.

Is there any way around this problem, without sacrificing performance? I tried this with Eigen 3.0 beta 2 and with the latest from hg.

Kind regards,

Bart

Ah, classical question

Your analysis is correct. Of course we've had to solve this before in Eigen itself (everytime you use chained matrix products with Eigen), I'll try to make you some example code tomorrow if nobody beats me to it (Gael, please go ahead ).

There is indeed an issue whenever you try to return a complex expression containing a product expression because this later gets evaluated when you nest it. It is like doing:

Code: Select all

Transpose<Mat> foo() {
  Mat temp = ...;
  return temp.transpose();
}

Here temp will be nested by reference, because this is a Matrix, and then the returned object will refer to a dead object.

The only solution is to enforce nesting by value:

Code: Select all

Transpose<NestByValue<Mat> > foo() {
  Mat temp = ...;
  return temp.nestByValue().transpose();
}

Unfortunately, this NestByValue wrapper is not compatible with general expressions yet, and it requires a specialization:

Code: Select all

template<typename Lhs, typename Rhs, int Mode> class NestByValue<GeneralProduct<Lhs,Rhs,Mode> >
  : public ProductBase<NestByValue<GeneralProduct<Lhs,Rhs,Mode> >,
                       typename GeneralProduct<Lhs,Rhs,Mode>::_LhsNested,
                       typename GeneralProduct<Lhs,Rhs,Mode>::_RhsNested>
{
  public:

    typedef GeneralProduct<Lhs,Rhs,Mode> NestedProduct;
   
    typedef ProductBase<NestByValue,
                       typename NestedProduct::_LhsNested,
                       typename NestedProduct::_RhsNested> Base;
    typedef typename Base::Scalar Scalar;
    typedef typename Base::PlainObject PlainObject;

    NestByValue(const NestedProduct& prod)
      : Base(prod.lhs(),prod.rhs()), m_prod(prod) { }
   
    template<typename Dest>
    inline void evalTo(Dest& dst) const { dst.setZero(); scaleAndAddTo(dst,Scalar(1)); }

    template<typename Dest>
    inline void addTo(Dest& dst) const { scaleAndAddTo(dst,Scalar(1)); }

    template<typename Dest>
    inline void subTo(Dest& dst) const { scaleAndAddTo(dst,-Scalar(1)); }

    template<typename Dest>
    inline void scaleAndAddTo(Dest& dst,Scalar alpha) const { m_prod.derived().scaleAndAddTo(dst,Scalar(1)); }

    operator const NestedProduct&() const { return m_prod; }
   
  protected:
    const NestedProduct m_prod;
};

and for some mysterious reasons this still does not work when you nest a outer product as in your example (but it does work for other cases), and surprisingly this does not generate extra temporaries... This part of Eigen became so complex... but the problem *is* complex. A more elegant and general solution would be give up evaluating at the expression construction time, but defer it at the evaluation time through an analysis of the expression more like NT²/boost::proto does....

What I don't get is, why would NestByValue still be needed now that we nest expressions by value by default? Is it that matrix products are special in this respect because they may hold temporary matrices which one doesn't want to nest by value?

I actually like our approach (evaluate at construction time) better than NT2's because our approach means that we can take advantage of these evaluations to limit the expression tree depth, so we scale better to deep chained products.

We just need to make sure that we understand our own code and document it, and are able to answer questions such as his.

What I don't get is that our own operator* in Product.h seems to behave exactly like his function, so why is it that our operator* works and his function doesn't ?

ggael wrote:and for some mysterious reasons this still does not work when you nest a outer product as in your example (but it does work for other cases), and surprisingly this does not generate extra temporaries... This part of Eigen became so complex... but the problem *is* complex. A more elegant and general solution would be give up evaluating at the expression construction time, but defer it at the evaluation time through an analysis of the expression more like NT²/boost::proto does....

Thanks for the explanation about NestByValue, I had seen this and wondered if it would solve my problem, but could not figure out how to use it. Regarding the exact case of my example, I had actually managed to fix it, by adding the following specialization at line 235 of Core/products/CoeffBasedProduct.h:

Code: Select all

template<typename Lhs, typename Rhs, int N, typename PlainObject>
struct ei_nested<CoeffBasedProduct<Lhs,Rhs,NestByRefBit>, N, PlainObject>
{
  typedef CoeffBasedProduct<Lhs,Rhs,NestByRefBit> const type;
};

Not sure if this is a good solution or not (performance wise), but I think the CoeffBasedProduct itself is light-weight so it should be OK to pass by value. I thought it was irrelevant, because I still ran into trouble with more complex expressions. Maybe this combined with the NestByValue for GeneralProduct would work in all cases, I'll try this ASAP.

bjacob wrote:What I don't get is, why would NestByValue still be needed now that we nest expressions by value by default? Is it that matrix products are special in this respect because they may hold temporary matrices which one doesn't want to nest by value?

Well, only pure expressions are nested by value, not the ones with storage, i.e., not Matrix, Array, *Product*, etc. So yes NestByValue is still needed for them.

What I don't get is that our own operator* in Product.h seems to behave exactly like his function, so why is it that our operator* works and his function doesn't ?

There is a big difference between:

D = A * (B * A);

and

D = multiply_on_both_side(A,B);

where multiply_on_both_side tries to return a chained product expression.

In the former case the temporary expression holding the result of (B * A) is valid during the whole lifetime of the expression...


Now, perhaps there exist a simpler solution to always nest by value, e.g., by relying on the fact the compiler is smart enough to remove copies of temporaries, or, in the case of dynamic sized object, we could add a "temporary" runtime flag enabling shallow copies, or in ProductBase, we could try to explicitly do a shallow copy in the copy ctor using swap... For statically allocated objects, however, I'm stuck. I don't see any alternative than relying on the smartness of the compiler.


I'll open a bug about that topic.

ggael wrote:There is a big difference between:

D = A * (B * A);

and

D = multiply_on_both_side(A,B);

where multiply_on_both_side tries to return a chained product expression.

In the former case the temporary expression holding the result of (B * A) is valid during the whole lifetime of the expression...

Oh, right, because it's created in the caller's stack frame. Ok, thanks for the explanation!

Now, perhaps there exist a simpler solution to always nest by value, e.g., by relying on the fact the compiler is smart enough to remove copies of temporaries, or, in the case of dynamic sized object, we could add a "temporary" runtime flag enabling shallow copies, or in ProductBase, we could try to explicitly do a shallow copy in the copy ctor using swap... For statically allocated objects, however, I'm stuck. I don't see any alternative than relying on the smartness of the compiler.

Right now, is it completely impossible, or just very hard to write a function returning a Eigen matrix product expression?

Is it enough to put a NestByValue?

At least one thing that we need, is good documentation. If it's hard or impossible to write functions returning product expressions, that needs to be well documented. We should have a special topic page on functions returning eigen expressions, mentioning the auto keyword, and explaining how matrix products are a special case.

I'll open a bug about that topic.[/quote]

I see, http://eigen.tuxfamily.org/bz/show_bug.cgi?id=99, let's continue the discussion there.