0

我尝试编写自定义预处理器，但到目前为止不了解 Eigen 预处理器的概念。我现在的状态看起来是这样的，但是有什么明显的错误……

class NullSpaceProjector: public Eigen::IdentityPreconditioner
{
  public:
    Eigen::MatrixXd null_space_1;
    Eigen::MatrixXd null_space_2;

    template<typename Rhs>
    inline Rhs solve(Rhs& b) const { 
    return b - this->null_space_1 * (this->null_space_2 * b);
    }

    void setNullSpace(Eigen::MatrixXd null_space) {
    // normalize + orthogonalize the nullspace
    this->null_space_1 = null_space * ((null_space.transpose() * null_space).inverse());
    this->null_space_2 = null_space.transpose();
    }
};

我想设置一些零空间信息并将其从每个 cg 步骤的 rhs 中删除。不知何故，我觉得这在某种程度上是不可能的预处理器概念。至少我找不到正确的方法来实现它。

参考：
-求解奇异矩阵
- https://bitbucket.org/eigen/eigen/src/1a24287c6c133b46f8929cf5a4550e270ab66025/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h?at=default&fileviewer=file-view-default#BasicPreconditioners.h-185

更多信息：

零空间是这样构造的：

Eigen::MatrixXd LscmRelax::get_nullspace()
{
    Eigen::MatrixXd null_space;
    null_space.setZero(this->flat_vertices.size() * 2, 3);

    for (int i=0; i<this->flat_vertices.cols(); i++)
    {
        null_space(i * 2, 0) =  1;  // x-translation
        null_space(i * 2 + 1, 1) =  1;  // y-translation
        null_space(i * 2, 2) =  -  this->flat_vertices(1, i);  // z-rotation
        null_space(i * 2 + 1, 2) = this->flat_vertices(0, i);  // z-rotation

    }
    return null_space;
}

所以这个问题可能与 eigen 没有太大关系，而是一个基本的 c++ 问题：如何通过依赖于 rhs 的投影来操作 const 函数中的 rhs。所以我猜 const 概念不符合我的需要。

但是从函数中删除“const”会导致如下错误：

.../include/eigen3/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h:63:5: error: passing ‘const lscmrelax::NullSpaceProjector’ as ‘this’ argument of ‘Eigen::VectorXd lscmrelax::NullSpaceProjector::solve(const VectorXd&)’ discards qualifiers [-fpermissive]
4

0 回答 0