0

I try to write a custom preconditioner, but do not understand the concept of Eigen's preconditioner so far. My current state looks like this, but there is something obvious wrong...

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();
    }
};

I want to set some null-space information and remove it from the rhs in every cg-step. Somehow I have the feeling this is somehow not possible with the preconditioner concept. At least I can't find the right way to get it implmented.

ref.:
- solving a singular matrix
- https://bitbucket.org/eigen/eigen/src/1a24287c6c133b46f8929cf5a4550e270ab66025/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h?at=default&fileviewer=file-view-default#BasicPreconditioners.h-185

further information:

the nullspace is constructed this way:

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

So the question is maybe not so much related to eigen but more a basic c++ question: How can I manipulate the rhs in a const function by a projection which depends on the rhs. So I guess the const concept doesn't fit my needs.

But removing the 'const' from the function results in errors like this:

.../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]
Lo L
  • 21
  • 3
  • What behavior do you observe/expect? How does `null_space` look, typically? – chtz Mar 03 '18 at 08:49
  • The problem is the programming interface. If I use this modification to the rhs the programm crashes. I do not have good possibilities to debug the problem and I guess the matrix multiplication should be right. The null-space is simple a matrix consisting of vectors which have no impact on the system: A * N = 0 in the precontitioner I want to create an orthonormal matrix and project it on the rhs. Therefor I used [this equation:](http://www.math.lsa.umich.edu/~speyer/417/OrthoProj.pdf) This should remove rigid-body movements from the solution. – Lo L Mar 03 '18 at 17:45
  • btw. after checking again I found some cases where my implemetation indeed worked. In this cases the result of the system was 0. So I guess the rhs was also 0. So it's possible the precontitioner is not called at all or it is called but as the projection of the nullspace 0 it has no impact. – Lo L Mar 03 '18 at 18:08

0 Answers0