Solving a linear system using ISTL

Here we show how to solve a linear system of equations originating from a PDE directly using ISTL instead of PDELab's abstractions.

Note that generally, we recommend going the all-PDELab route as explained in Solving a linear system within PDELab.

First, we assemble our residual as in Assembling a linear system from a PDE

  X d(gfs,0.0);
  go.residual(x,d);

as well as the system matrix we want to solve:

  typedef GO::Jacobian M;
  M A(go);
  go.jacobian(x,A);

So far, these are data types provided by PDELab. We can access the underlying ISTL data types via

  typedef Dune::PDELab::Backend::Native<M> ISTLM;
  typedef Dune::PDELab::Backend::Native<X> ISTLX;

and the ISTL vectors or matrices behind a PDELab object via:

  using Dune::PDELab::Backend::native;

Now we can wrap our matrix in a LinearOperator object

  typedef Dune::MatrixAdapter<ISTLM,ISTLX,ISTLX> Operator;
  Operator matrixop(native(A));

and pass that into our desired solver along with a preconditioner of our choice.

  Dune::SeqJac<ISTLM,ISTLX,ISTLX> preconditioner(native(A), 1, .5);
  Dune::CGSolver<ISTLX> solver(matrixop, preconditioner, 1e-3, 10, 2);

Applying the solver now returns the solution of the defect problem, i.e. A*v=d.

  X v(gfs,0.0);
  Dune::InverseOperatorResult res;
  solver.apply(native(v), native(d), res);

Finally, we can subtract that from our initial guess to obtain the solution of A*x=b:

  x -= v;

In the end, let's print the solution to console.

  Dune::printvector(std::cout, native(x), "Solution", "");

Full example code: recipe-linear-system-solution-istl.cc