Eigen3 Sparse Solver не копируемый

Я работаю над числовым кодом и хочу оценить, как различаются Sparse и Dense Matrix-LU (а также и другие, более поздние) для варианта использования кода.
Eigens Dense Decomposition Objects могут быть копируемыми, и это используется для их кэширования, используя boost :: Вариант, для большей гибкости позже.

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

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

Спасибо 🙂

/// -------------------------------- DENSE APPROACH, WORKS -------------------------------

using CacheType = boost::variant<Eigen::FullPivLU<Eigen::MatrixXd>,
Eigen::PartialPivLU<Eigen::MatrixXd>>;

// visit the variant, and solve with the correct decomposition
struct DenseVisitor : boost::static_visitor<Eigen::MatrixXd> {
DenseVisitor(Eigen::MatrixXd const& r) : rhs{r} {}

template <class Decomposition>
Eigen::MatrixXd operator()(Decomposition const& d) const
{
Eigen::MatrixXd res = d.solve(rhs);
return res;
}

private:
Eigen::MatrixXd const& rhs; // reference to rhs, since () will take only one argument
};

// part of a class, having a cachetype as member
Eigen::MatrixXd solve(Eigen::MatrixXd const& A, Eigen::MatrixXd const& b)
{
// decompose if we now we changed A, and save the decomposition of A
if(cache_dirty) {
cache_ = A.partialPivLU();
cache_dirty = false;
}

// solve rhs with cached decomposition
auto result = boost::apply_visitor(DenseVisitor(b), cache_);
return result;
}/// ------------------------- SPARSE APPROACH, WORKS NOT ---------------------------------

// will be extended later on, but for now thats enough
using CacheType = boost::variant<Eigen::SparseLU<Eigen::SparseMatrix<double>>>;

// visit the variant, and solve with the correct decomposition
struct SparseVisitor : boost::static_visitor<Eigen::MatrixXd> {
SparseVisitor(Eigen::MatrixXd const& r) : rhs{r} {}

template <class Decomposition>
Eigen::MatrixXd operator()(Decomposition const& d) const
{
Eigen::MatrixXd res = d.solve(rhs);
if (d.info() != Eigen::Success)
throw std::runtime_error{"Sparse solve failed!"};

return res;
}

private:
Eigen::MatrixXd const& rhs; // reference to rhs, since () will take only one argument
};

// part of a class, having a cachetype as member, and a Pointer to A
// so the cache will only solve for b, and if necessary recompute the decomposition
Eigen::MatrixXd solve(Eigen::SparseMatrix<double>& A, Eigen::MatrixXd const& b)
{
// get decomposition, this will be extended by a visitor as well!
auto* decomp = boost::get<Eigen::SparseLU<Eigen::SparseMatrix<double>>>(cache_);
// decompose if we now we changed A, and save the decomposition of A
if(cache_dirty) {

// reanalyze the pattern
if (reanalyze) {
A.makeCompressed();
decomp->analyzePattern(A);
}

// factorize
decomp->factorize(A);

if(decomp->info() != Eigen::Success)
throw std::runtime_error{"Sparse decomposition failed"};

cache_dirty = false;
}

// solve rhs with cached decomposition
auto result = boost::apply_visitor(SparseVisitor(b), cache_);
return result;
}