这个问题与从 Eigen::CwiseBinaryOp 到 MatrixXd 的转换导致段错误有关。它可能会有像前者一样简单的解决方案。
在这个最小的示例中,我定义了 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;
}
但是,如果您在头文件中使用注释掉的行,则一切正常。不幸的是,这会将分解的评估转移到求解器的构造上,这不适合我的使用。通常,我想构建一个涉及此Solve
的复杂表达式expr
,并在稍后调用expr.get()
。
如何解决这个问题?是否有一般规则可以遵循,以避免进一步的相关问题?
为了避免无用和昂贵的副本,内部solve_retval
结构存储分解和右侧的常量引用。但是,在 Decomp::get
函数中创建的LDLT
对象将在此函数返回的同时被删除,因此solve_retval
对象引用死对象。
一种可能的解决方法是在 Solve
中添加一个 Decomp::result_type
对象,并在 Solve::get
中对其进行初始化。此外,为了避免多个深度拷贝,我建议对一些属性使用 const 引用,如下所示:
#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());
}
};
一般规则是通过 const 引用(以避免深度复制)和按值存储轻量级表达式(以减少临时生命问题)存储重对象。