Eigen3 Sparse Solver noncopyable

Question

iam working on a numerical code and want to evaluate how Sparse and Dense Matrix-LU decomposition (and later others as well) differ for the usecase of the code. Eigens Dense Decomposition Objects can be copyable, and that is used to cache these, using a boost::variant, for more flexibility later.

I want to achieve the same with the Sparse Solvers, but have my difficulties doing that. The minimal example below should illustrate how my approach is.

The question is, why are sparse-solvers non copyable? Can i just write my own copy operations, or are they definitly illformed. How could i work around that problem?

Thank you :)

/// -------------------------------- 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;
}

Plod