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

Subject: [proto] Using proto with expressions containing matrices from the EIgen library
From: Bart Janssens (bart.janssens_at_[hidden])
Date: 2010-10-27 06:50:17

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.

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.

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

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

Kind regards,

-- 
Bart