Proto for Boost mailing list page: Re: [proto] Using proto with expressions containing matrices from the EIgen library

Subject: Re: [proto] Using proto with expressions containing matrices from the EIgen library
From: Eric Niebler (eric_at_[hidden])
Date: 2010-10-27 12:57:29

On 10/27/2010 3:50 AM, Bart Janssens wrote:
> Hi all,
> 
> I am trying to use proto for the evaluation of expressions involved in
> a Finite Element Method. The idea is that developers can write
> expressions that resemble the mathematical formulation of the problem,
> but that get processed efficiently to use matrix operations on
> fixed-size matrices. The matrix library that is used as backend is
> Eigen (http://eigen.tuxfamily.org).
> 
> In my first attempt, I used proto::eval with a context that
> understands some custom terminals and functions, and delegates other
> stuff (like the math operators) to the proto default context.
> Unfortunately, this approach works only on expressions involving at
> most 2 matrices, since in more complex cases Eigen constructs its own
> nested expression templates and these may refer to their operands by
> reference. The operands may be temporaries, however, (as generated by
> a nested call to proto::eval), so the program either crashes or gives
> the wrong result. In a release compile, everything works, presumably
> because all the nested proto calls get inlined.

Are you saying that proto::eval is generating temporary matrices? I
don't understand why, but there's always a way to get Proto to store
intermediate nodes by value instead of by reference. Without knowing
what you're doing, I can't say more.

> In my second attempt, I tried to use an object transform, hoping the
> temporaries would persist until the evaluation happens. The code looks
> like this:
> 
> template<typename Left, typename Right>
> struct MultExpr
> {
>   typedef typename Eigen::ProductReturnType<Left, Right>::Type type;
> };
> 
> struct EvaluateExpr :
>   boost::proto::or_
>   <
>     boost::proto::when
>     <
>       boost::proto::multiplies<boost::proto::_, boost::proto::_>,
>       MultExpr
>       <
>         EvaluateExpr(boost::proto::_left), EvaluateExpr(boost::proto::_right)
>       >(EvaluateExpr(boost::proto::_left), EvaluateExpr(boost::proto::_right))
>     >,
>     boost::proto::when
>     <
>       boost::proto::_,
>       ContextEvaluator
>     >
>   >
> {};
> 
> An example of the expression to evaluate:
> transpose(mapped_gradient(nodes)) * transpose(jacobian_adjoint(nodes))
> * jacobian_adjoint(nodes)
> 
> where nodes is a terminal that can be evaluated using a context. In
> the EvaluateExpr grammar, ContextEvaluator is a primitive transform
> that simply delegates evaluation to proto::eval, using the data
> argument as context. The context stores the matrices involved, so any
> references to matrices returned are safe to use.

If you've gone to the trouble writing a primitive transform, why bother
using proto::eval at all? Do it all in the transform. You might find
your bug that way.

> The Eigen product expression is now constructed using object transform
> MultExpr, but it turns out to exhibit the same problem with
> temporaries.

Can ContextEvaluator ever return a reference to a temporary object? Is
Eigen::ProductReturnType<Left, Right>::Type ever a reference? Does that
type hold Left and Right by reference?

> Does anyone see a solution to this? At the moment I am thinking of
> having the EvaluateExpr transform build a fusion::cons list to store
> all the intermediate expressions, and refer to those by reference in
> subsequent expressions, but I'm not sure how to pull this off yet, or
> if there is a simpler way.
> 
> I have also asked the question on the Eigen forum, together with a
> small sample that exposes the temporary variable problem without using
> proto:
> http://forum.kde.org/viewtopic.php?f=74&t=91119

Looking at that code, I strongly suspect that
Eigen::ProductReturnType<Left, Right>::Type is a type that holds Right
by reference. In fact, I'm 99% sure it does. Print the name of the
resulting type of Eigen::ProductReturnType<Left, Right>::Type and see if
that's the case. If so, you've found your bug. Can you see why?

-- 
Eric Niebler
BoostPro Computing
http://www.boostpro.com