Hi folks,

during development I have run into a problem with Eigen expressions. I have tried to create a minimal example (see below).

The classes MockEquationSystem and MockTimeDeriv are mock classes for unit tests and should represent components of an equation system. I have some additional glue code that adds them together and should produce code as in the lambda g. I have also created the other lambda f which should in theory produce the same result as g. However when I run it, then the result of f is printed and as soon as g() is called, the program segfaults.

I compiled it with:

Code: Select all

g++ -Wall -Wextra -g --std=c++14 filename.cpp


with gdb I can reproduce the following back trace:

Code: Select all

#0  0x0000000000404010 in Eigen::internal::scalar_difference_op<double>::operator() (this=0x7fffffffe270, a=@0x14: <error reading variable>, b=@0x62c330: -0.50407537515464962) at /usr/include/eigen3/Eigen/src/Core/Functors.h:190
#1  0x00000000004042d7 in Eigen::CwiseBinaryOpImpl<Eigen::internal::scalar_difference_op<double>, Eigen::GeneralProduct<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 4> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::Dense>::coeff (this=0x7fffffffe260, index=0) at /usr/include/eigen3/Eigen/src/Core/CwiseBinaryOp.h:189
#2  0x0000000000403686 in Eigen::CwiseBinaryOpImpl<Eigen::internal::scalar_sum_op<double>, Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<double>, Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<double>, Eigen::CwiseUnaryOp<Eigen::internal::scalar_quotient1_op<double>, Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> > const> const> const, Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double>, Eigen::GeneralProduct<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 4> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const> const, Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double>, Eigen::GeneralProduct<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 4> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const> const, Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double>, Eigen::GeneralProduct<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 4> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const, Eigen::Dense>::coeff (this=0x7fffffffe1d0, index=0) at /usr/include/eigen3/Eigen/src/Core/CwiseBinaryOp.h:188
#3  0x0000000000402622 in Eigen::DenseCoeffsBase<Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<double>, Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<double>, Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<double>, Eigen::CwiseUnaryOp<Eigen::internal::scalar_quotient1_op<double>, Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> > const> const> const, Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double>, Eigen::GeneralProduct<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 4> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const> const, Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double>, Eigen::GeneralProduct<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 4> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const> const, Eigen::CwiseBinaryOp<Eigen::internal::scalar_difference_op<double>, Eigen::GeneralProduct<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 4> const, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const>, 0>::operator() (this=0x7fffffffe1d0, index=0) at /usr/include/eigen3/Eigen/src/Core/DenseCoeffsBase.h:174
#4  0x000000000040125e in main () at segfaul.cpp:64


Interestingly, if I let the compiler optimize (additional -O3), then there is no segfault. But the output of f() and g() is not equal.

Also with experimenting, I found out that the problem is in MockTimeDeriv, specifically in the part that eval returns the parameter x. When it doesn't return it, everything is fine. However I don't really get this, as MockEquationSystem also returns a function of x.

Code: Select all

#include <eigen3/Eigen/Core>
#include <iostream>

using Eigen::MatrixXd;
using Eigen::VectorXd;

class MockEquationSystem {
public:
  explicit MockEquationSystem(const std::size_t N)
      : m_A(MatrixXd::Random(N, N)), m_b(VectorXd::Random(N)) {}

  MockEquationSystem(const MockEquationSystem &other) = delete;

  auto eval(const VectorXd &x) const { return A() * x - b(); }

  const MatrixXd &A() const { return m_A; }

  const VectorXd &b() const { return m_b; }

private:
  const MatrixXd m_A;
  const VectorXd m_b;
};

class MockTimeDeriv {
public:
  MockTimeDeriv(const std::size_t N, const double dt) : m_dt(dt), m_N(N) {}

  MockTimeDeriv(const MockTimeDeriv &other) = delete;

  double time_step() const { return m_dt; }

  auto eval(const VectorXd &x) const {
    return (x - VectorXd::Ones(N())) / time_step();
  }

private:
  const double m_dt;
  std::size_t N() const { return m_N; }
  const std::size_t m_N;
};

int main() {
  const int N = 20;

  const MockEquationSystem m1(N);
  const MockEquationSystem m2(N);
  const MockEquationSystem m3(N);
  const MockTimeDeriv mTD(N, 0.1);

  const VectorXd x = VectorXd::Random(N);

  auto f = [&mTD, &m1, &m2, &m3](auto x) {
    return (x - VectorXd::Ones(N)) / mTD.time_step() + m1.A() * x - m1.b() +
           m2.A() * x - m2.b() + m3.A() * x - m3.b();
  };

  auto g = [&mTD, &m1, &m2, &m3](auto x) {
    return mTD.eval(x) + m1.eval(x) + m2.eval(x) + m3.eval(x);
  };

  for (std::size_t i = 0; i < N; i++) {
    std::cout << f(x)(i) << std::endl;
    std::cout << g(x)(i) << std::endl;
  }

  return 0;
}


Thank you in advance for any suggestions since I am stuck.

Please don't abuse of the "auto" keyword when using expression template libraries. Your code creates a lot of expression referencing dead objects, see: http://eigen.tuxfamily.org/dox-devel/TopicPitfalls.html

Hi,

I was puzzled by your answer at first, since I can evaluate the single components without any problems, so I thought that shouldn't be a problem. So I forbid Eigen to allocate memory for temporaries and found out that MockEquationSystem seems to create a temporary during a call to eval(). I guess that could be probably one of the issues that leads to the segfault.

This leads me to another question. Provided I should not use auto, how can I create expressions otherwise in an elegant way? I would like to be able to create something as the function f() does for different cases (e.g. only 1 or 2 MockEquationSystem or completely different components). How could I accomplish that without having to write every single case by hand (when I should not create functions return expressions)?

Thanks for any suggestions

Yes, .eval() evaluate the expression and returns a temporary Matrix<> object (or Array or SparseMatrix, etc.).

Ok, auto is not always a bad idea, if you are 100% sure that 1) you really want an expression and not a Matrix-like object type, and 2) that no temporary has been created on the way. Regarding this second aspect, I strongly recommend you to move to Eigen3.3 (beta1 or devel branch) because in 3.2, some subexpressions are directly evaluated within a temporary (e.g., matrix products).

So in your case this concerns matrix-vector products like "m1.A() * x". You can either move to 3.3, or if N is typically very small (say <20), move to lazy-products: m1.A().lazyProduct(x) which lazily evaluates the product in a coefficient-wise manner.

I have tried your suggestion with using 3.3 and it fixed the segfault. Thank you for your help!