This question is related to cast from Eigen::CwiseBinaryOp to MatrixXd causes segfault . It will probably have as simple a solution as the former.
In this minimal example, I define Holder
, which holds and Eigen matrix, and returns it via its get()
member function. Similarly, Decomp
is an expression template for the LDLT decomposition of this matrix, and Solve
solves for AX=B, yielding X.
#include <Eigen/Dense>
#include <Eigen/Cholesky>
template <class EigenType> class Holder {
public:
typedef EigenType result_type;
private:
result_type in_;
public:
Holder(const EigenType& in) : in_(in) {}
const result_type& get() const { return in_; }
};
template <class Hold> class Decomp {
public:
typedef typename Eigen::LDLT
<typename Eigen::MatrixBase<typename Hold::result_type>::PlainObject>
result_type;
private:
Hold mat_;
public:
Decomp(const Hold& mat) : mat_(mat) {}
result_type get() const { return mat_.get().ldlt(); }
};
template <class Derived, class OtherDerived> class Solve {
public:
typedef typename Eigen::internal::solve_retval
<typename Derived::result_type, typename OtherDerived::result_type>
result_type;
private:
Derived decomp_;
// typename Derived::result_type decomp_;
OtherDerived mat_;
public:
Solve(const Derived& decomp, const OtherDerived& mat)
: decomp_(decomp), mat_(mat) {}
//: decomp_(decomp.get()), mat_(mat) {}
result_type get() const { return decomp_.get().solve(mat_.get()); }
// result_type get() const { return decomp_.solve(mat_.get()); }
};
typedef Holder<Eigen::MatrixXd> MatrixHolder;
typedef Decomp<MatrixHolder> MatrixDecomp;
typedef Solve<MatrixDecomp, MatrixHolder> SimpleSolve;
The following test fails on X.get()
#include "Simple.h"
#include <Eigen/Dense>
#include <iostream>
int main(int, char * []) {
MatrixHolder A(Eigen::MatrixXd::Identity(3, 3));
MatrixHolder B(Eigen::MatrixXd::Random(3, 2));
MatrixDecomp ldlt(A);
SimpleSolve X(ldlt, B);
std::cout << X.get() << std::endl;
return 0;
}
but if you use the commented out lines in the header file, everything works fine. Unfortunately, this moves the evaluation of the decomposition to the construction of the solver, which is not suited to my use. Typically, I want to build a complicated expression expr
involving this Solve
, and call expr.get()
later on.
How can I solve this problem? Is there a general rule to follow, to avoid further related questions?
To avoid useless and expensive copies, the internal
solve_retval
structure stores the decomposition and right-hand-side by const references. However, theLDLT
object created in theDecomp::get
function is deleted at the same time this function returns, and thus thesolve_retval
object references dead objects.One possible workaround, is to add a
Decomp::result_type
object inSolve
, and initialize it inSolve::get
. Moreover, to avoid multiple deep copies, I suggest to use const references for a few attributes as follow:The general rule is to store heavy objects by const reference (to avoid deep copies) and lightweight expressions by value (to reduce temporary life issues).