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

Question

Я работаю над числовым кодом и хочу оценить, как различаются 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;
}