On using block diagonal preconditioner for the mixed Poisson problem

[Charlie Z]

Dear all,

I'm trying to use a block-diagonal preconditioner for the mixed Poisson problem from step-20, but haven't made it yet.

The diagonal preconditioner of interest is from Powell & Silvester,2003. Specifically, for the matrix system after discretization with the RT0-DG0 element pair

  [A  BT

   B  0],

they suggested a block diagonal preconditioner

  P = 

  [DA  0

   0     PS],

where DA is diag(A), and PS is AMG applied to the approximated Schur complement SD = B(DA^-1)BT.

As the above idea is quite similar to that in step-55, so I borrowed quite a few code blocks from there, e.g. InverseMatrix and BlockDiagonalPreconditioner. The resulting solve() with a SolverMinRes solver is as follows. (The rest of the code is mainly from step-20)

  template <int dim>

 void MixedLaplaceProblem<dim>::solve()

 {

   TrilinosWrappers::PreconditionJacobi diagA;

   diagA.initialize(system_matrix.block(0, 0));//DA

    ApproximateSchurComplement approximate_schur(system_matrix, solution);//SD

    InverseMatrix<ApproximateSchurComplement> approximate_inverse(approximate_schur);//PS

    const BlockDiagonalPreconditioner<TrilinosWrappers::PreconditionJacobi,

         InverseMatrix<ApproximateSchurComplement> >

         preconditioner(diagA, approximate_inverse);

    SolverControl solver_control(system_matrix.m(),

                                1e-6 * system_rhs.l2_norm());

    SolverMinRes<TrilinosWrappers::MPI::BlockVector> solver(solver_control);

   solver.solve(system_matrix, solution, system_rhs,

      preconditioner

     );

    std::cout<<"last step = "<<solver_control.last_step()<<std::endl;

 }

However, the code can't compile, and gives the following error message:

  /home/dealii/test.cc:336:5: error: no matching function for call to ‘dealii::TrilinosWrappers::PreconditionAMG::initialize(const Step20::ApproximateSchurComplement&, dealii::TrilinosWrappers::PreconditionAMG::AdditionalData&)’

      preconditioner.initialize(*matrix, data);

      ^

 In file included from /home/dealii/dealii-v9.0.0/include/deal.II/lac/generic_linear_algebra.h:141:0,

                  from /home/dealii/test.cc:31:

 /home/dealii/dealii-v9.0.0/include/deal.II/lac/trilinos_precondition.h:1516:10: note: candidate: void dealii::TrilinosWrappers::PreconditionAMG::initialize(const dealii::TrilinosWrappers::SparseMatrix&, const dealii::TrilinosWrappers::PreconditionAMG::AdditionalData&)

      void initialize (const SparseMatrix   &matrix,

The error is related to the use of PreconditionAMG in my implementation of InverseMatrix<MatrixType>::vmult:

  template <class MatrixType>

 void InverseMatrix<MatrixType>::vmult (LinearAlgebraTrilinos::MPI::Vector       &dst,

                                        const LinearAlgebraTrilinos::MPI::Vector &src) const

 {

   SolverControl solver_control(src.size(),

                            1e-6 * src.l2_norm());

   SolverCG<TrilinosWrappers::MPI::Vector>    cg (solver_control);

   dst = 0;

    LinearAlgebraTrilinos::MPI::PreconditionAMG preconditioner;

   LinearAlgebraTrilinos::MPI::PreconditionAMG::AdditionalData data;

   preconditioner.initialize(*matrix, data);

    // PreconditionIdentity preconditioner;//A Identity preconditioner works fine here

    cg.solve (*matrix, dst, src, preconditioner);

    std::cout<<"inner last step = "<<solver_control.last_step()<<std::endl;

 }

I find a conflict here is that::

1) TrilinosWrappers::PreconditionAMG can only be initialized with a matrix, but not a ApproximateSchurComplement;

2) step-20 uses ApproximateSchurComplement because we actually don't want to/can't assemble the Schur matrix or its inverse.

Could you please give me a hint on how to solve this conflict in deal.II? Is it actually possible to assemble this B(DA^-1)BT directly?

Attached is a minimal (but not working) code.

Best regards,

Charlie

[Wolfgang Bangerth]

On 8/27/19 2:23 PM, Charlie Z wrote:
>
> I find a conflict here is that::
> 1) TrilinosWrappers::PreconditionAMG can only be initialized with a
> matrix, but not a ApproximateSchurComplement;
> 2) step-20 uses ApproximateSchurComplement because we actually don't
> want to/can't assemble the Schur matrix or its inverse.

Yes, you correctly identified the issue. step-20 gets away with using
ApproximateSchurComplement as a preconditioner because, from the
perspective of GMRES, a preconditioner only has be able to provide a
matrix-vector product (i.e., apply itself to a vector).

But if you want to use an object to initialize a algebraic multigrid,
the AMG implementation needs to be able to do much more than just apply
the object to a vector -- it needs to build a multilevel hierarchy, and
for that it needs access to the elements that make up that matrix. The
problem is of course that ApproximateSchurComplement does not offer
access to the elements of the matrix it represents (because it never
explicitly computes them).

> Could you please give me a hint on how to solve this conflict in
> deal.II? Is it actually possible to assemble this B(DA^-1)BT directly?

Yes, of course. You can just multiply it out as the matrix-matrix
product of three matrices. The result is a square matrix of the right
size, with a rather large number of nonzero entries per row.

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Charlie Z]

Hi Wolfgang,

Thanks very much for your clarification.

The code works fine now with the help of TrilinosWrappers::SparseMatrix::mmult(), though, I understand, mmult is expensive. 

Regards,

Charlie

[Wolfgang Bangerth]

On 8/28/19 6:23 AM, Charlie Z wrote:
>
> The code works fine now with the help of
> TrilinosWrappers::SparseMatrix::mmult(), though, I understand, mmult is
> expensive.

Glad to hear! Yes, mmult is expensive, but then you're only doing it once.
Unless you have concrete evidence to the contrary, performance problems are
not worth chasing.