6

我正在研究数字代码，并希望评估稀疏和密集矩阵-LU 分解（以及后来的其他分解）对于代码的用例有何不同。Eigens Dense Decomposition Objects 可以是可复制的，用于缓存这些，使用 boost::variant，以便以后获得更大的灵活性。

我想用稀疏求解器达到同样的效果，但这样做有困难。下面的最小示例应该说明我的方法。

问题是，为什么稀疏求解器不可复制？我可以只编写自己的复制操作，还是它们绝对不正确。我该如何解决这个问题？

谢谢 ：）

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

0 回答 0