Matrix-free linear solver implementation cpp

Question

I need to solve the linear system associated with the Schur complement of the Stokes problem to compute the pressure field. For this purpose I would like to use a Krylov solver with a matrix-free approach in order to avoid the storage of the dense matrix associated with the Schur complement.

I did not find any demo about the usage of a Krylov solver with a matrix-free linear operator and hence I tried by myself. Unfortunately I can't succeed in making it working. In the following the cpp code for a simple matrix-free test case derived from the 2D Poisson demo. Where I'm wrong?

class MyLinearOperator : public LinearOperator 
{
  private:

  MyLinearOperator(const GenericMatrix& M): _M(M) {};

  void mult(const GenericVector& x, GenericVector& y) const {
    _M.mult(x,y);
  };

  std::size_t size(std::size_t dim) const {
    return _M.size(dim);
  };

  private:

   const GenericMatrix& _M;

};

...

Poisson::BilinearForm a(V, V);
Poisson::LinearForm L(V);
Source f;
dUdN g;
L.f = f;
L.g = g;

Function u(V);
boost::shared_ptr<Matrix> A(new Matrix);
Vector b;  
assemble_system(*A, b, a, L, bc);

I have tested two different ways (case #) to setup the matrix-free solver based on the documented "KrylovSolver" class using either a reference (case 1) or a boost pointer
(case 2) to the matrix-free linear operator.

Case 1

MyLinearOperator MyA(*A);
KrylovSolver KMF_solver("gmres", "none");
KMF_solver.solve(MyA, *u.vector(), b);

Case 2

boost::shared_ptr<const MyLinearOperator> MyAp( new MyLinearOperator(*A));     
KrylovSolver KMF_solver("gmres", "none");
KMF_solver.set_operator(MyAp);
KMF_solver.solve(*u.vector(), b);

I have tested the code using both "PETSc" and "uBLAS". With PETSc I get a segmentation fault error for "case 1" while for "case 2" the following error is reported:

krylov-matrix-free: /usr/include/boost/smart_ptr/shared_ptr.hpp:418: T* boost::shared_ptr::operator->() const [with T = const dolfin::PETScBaseMatrix]: Assertion `px != 0' failed.
[pc-carini:22641] *** Process received signal ***
[pc-carini:22641] Signal: Aborted (6)
[pc-carini:22641] Signal code: (-6)
[pc-carini:22641] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x364a0) [0x7fcaaa09c4a0]
[pc-carini:22641] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0x35) [0x7fcaaa09c425]
[pc-carini:22641] [ 2] /lib/x86_64-linux-gnu/libc.so.6(abort+0x17b) [0x7fcaaa09fb8b]
[pc-carini:22641] [ 3] /lib/x86_64-linux-gnu/libc.so.6(+0x2f0ee) [0x7fcaaa0950ee]
[pc-carini:22641] [ 4] /lib/x86_64-linux-gnu/libc.so.6(+0x2f192) [0x7fcaaa095192]
[pc-carini:22641] [ 5] /usr/lib/libdolfin.so.1.2(+0x432240) [0x7fcaab06b240]
[pc-carini:22641] [ 6] /usr/lib/libdolfin.so.1.2(_ZN6dolfin17PETScKrylovSolver5solveERNS_11PETScVectorERKS1_+0xb5d) [0x7fcaab0cf6dd]
[pc-carini:22641] [ 7] /usr/lib/libdolfin.so.1.2(_ZN6dolfin12KrylovSolver5solveERNS_13GenericVectorERKS1_+0x148) [0x7fcaab0cba58]
[pc-carini:22641] [ 8] ./krylov-matrix-free(main+0x912) [0x40ec42]
[pc-carini:22641] [ 9] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xed) [0x7fcaaa08776d]
[pc-carini:22641] [10] ./krylov-matrix-free() [0x40f6dd]
[pc-carini:22641] *** End of error message ***
Aborted (core dumped)

With uBLAS no error is reported but the solver seems not work.

Thank you!

Comments

Please, post complete codes (case 1, 2) so we can try it. Report version of FEniCS also.

---

Here you are the complete code (CASE 2 is commented)

    #include <dolfin.h>
    #include "Poisson.h"

    using namespace dolfin;

    // Source term (right-hand side)
    class Source : public Expression
    {
      void eval(Array<double>& values, const Array<double>& x) const
      {
        double dx = x[0] - 0.5;
        double dy = x[1] - 0.5;
        values[0] = 10*exp(-(dx*dx + dy*dy) / 0.02);
      }
    };

    // Normal derivative (Neumann boundary condition)
    class dUdN : public Expression
    {
      void eval(Array<double>& values, const Array<double>& x) const
      {
        values[0] = sin(5*x[0]);
      }
    };

    // Sub domain for Dirichlet boundary condition
    class DirichletBoundary : public SubDomain
    {
      bool inside(const Array<double>& x, bool on_boundary) const
      {
        return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS;
      }
    };

    // --------------------------------------------------------

    class MyLinearOperator : public LinearOperator
    {

     public: 

     // Constructor

     MyLinearOperator(const GenericMatrix& M) _M(M) {};

     // Overloading matrix-vector multiplication

     void mult(const GenericVector& x, GenericVector& y) const {

      _M.mult(x,y);

     };

     // Overloading size function

     std::size_t size(std::size_t dim) const {

      return _M.size(dim);

     };

     std::string str(bool verbose) const {

      return "My Linear Operator";

     };

     private:

      const GenericMatrix& _M;

    };

    // --------------------------------------------------------

    int main()
    {

      set_log_level(PROGRESS);

      // Linear algebra backend

      parameters["linear_algebra_backend"] = "PETSc";

      // Create mesh and function space
      UnitSquareMesh mesh(32, 32);
      Poisson::FunctionSpace V(mesh);

      // Define boundary condition
      Constant u0(0.0);
      DirichletBoundary boundary;
      DirichletBC bc(V, u0, boundary);

      // Define variational forms
      Poisson::BilinearForm a(V, V);
      Poisson::LinearForm L(V);
      Source f;
      dUdN g;
      L.f = f;
      L.g = g;

      Function u(V);

      // Assembling linear system

      boost::shared_ptr<Matrix> A(new Matrix);
      Vector b,  res;  

      assemble_system(*A, b, a, L, bc);

      // Computing solution using Krylov matrix-free solver (no preconditioner)

      MyLinearOperator MyA(*A);

      KrylovSolver KMF_solver("gmres", "none");
      KMF_solver.parameters["relative_tolerance"] = 1.0e-10;
      KMF_solver.parameters["report"] = true;
      //KMF_solver.parameters["monitor_convergence"] = true;
      KMF_solver.solve(MyA, *u.vector(), b);
      MyA.mult(*u.vector(), res);
      res.axpy(-1.0, b);


      // CASE 2 (no preconditioner - pointer version)

      /*

      boost::shared_ptr<const MyLinearOperator> MyAp( new MyLinearOperator(*A));  

      KrylovSolver KMF_solver("gmres", "none");
      KMF_solver.set_operator(MyAp);
      KMF_solver.parameters["relative_tolerance"] = 1.0e-10;
      KMF_solver.parameters["report"] = true;
      KMF_solver.solve(*u.vector(), b);
      MyAp->mult(*u.vector(), res);
      res.axpy(-1.0, b);

      */

      info("Residual Krylov sol. matrix-free inf-norm %1.5e", res.norm("linf"));


     // Plot solution
      plot(u);
      interactive();

      return 0;

    }

---

I forgot: FEniCS version is 1.2.0 from ubuntu repository

Answer 1

Finally browsing the library source code dolfin/test/unit/la/cpp/LinearOperator.cpp I succeed in making the matrix-free implementation working by simply modifying the constructor of MyLinearOperator class as follows

MyLinearOperator(const GenericMatrix& M, const GenericVector& x): LinearOperator(x, x), _M(M) {};

where "x" is a dummy vector of the same size of the vector on which the linear opeartor acts which is employed within the initialization list by "LinearOperator(x, x)". I don't understand the meaning of this instruction since it is not reported within the LinearOperator class documentation: it seems you need to specify the input/output dimensions of the linear mapping.

Comments

Thanks for leaving the comment, it saved me a bad headache!
It looks like some methods are missing in the documentation [they are commented with two slashes instead of three in the code...].