Как отслеживать объекты Eigen посредством решения()?

Question

Этот вопрос относится к cast from Eigen::CwiseBinaryOp to MatrixXd causes segfault. Вероятно, это будет простое решение, как первое.

В этом минимальном примере я определяю Holder, который содержит и собственную матрицу, и возвращает его через функцию-член get(). Точно так же, Decomp это шаблон выражение для разложения LDLT этой матрицы, а Solve решает для AX = B, получая 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; 

Следующий тест терпит неудачу на 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; 
} 

, но если использовать прокомментированные строки в файле заголовка, все работает нормально. К сожалению, это переводит оценку разложения на конструкцию решателя, что не подходит для моего использования. Как правило, я хочу построить сложное выражение expr с участием этого Solve и позвонить по телефону expr.get() позже.

Как я могу решить эту проблему? Существует ли общее правило, чтобы избежать дальнейших связанных вопросов?

Answer 1

Чтобы избежать бесполезных и дорогих копий, внутренняя структура solve_retval хранит разложение и правую сторону по ссылкам const. Однако объект LDLT, созданный в функции Decomp::get, удаляется одновременно с возвратом этой функции, и поэтому объект solve_retval ссылается на мертвые объекты.

Возможное обходное решение состоит в том, чтобы добавить Decomp::result_type объект в Solve и инициализировать его в Solve::get. Кроме того, чтобы избежать многочисленных глубоких копий, я предлагаю использовать константные ссылки на несколько атрибутов следующим образом:

#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: 
    const Hold& mat_; 
    mutable result_type result_; 
    mutable bool init_; 

public: 
    Decomp(const Hold& mat) : mat_(mat), init_(false) {} 

    const result_type& get() const { 
    if(!init_) { 
     init_ = true; 
     result_.compute(mat_.get()); 
     return result_; 
    } 
    } 
}; 

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: 
    const Derived& decomp_; 
    const OtherDerived& mat_; 

public: 
    Solve(const Derived& decomp, const OtherDerived& mat) 
     : decomp_(decomp), mat_(mat) {} 

    result_type get() const { 
    return decomp_.get().solve(mat_.get()); 
    } 
}; 

Общее правило для хранения тяжелых предметов по константной ссылке (чтобы избежать глубоких копий) и легких выражений по стоимости (для сокращения временных сроков жизни).