petsc_solver.cc

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2019 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------


#include <deal.II/base/logstream.h>
#include <deal.II/base/std_cxx14/memory.h>

#include <deal.II/lac/petsc_solver.h>

#ifdef DEAL_II_WITH_PETSC

#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/petsc_compatibility.h>
#  include <deal.II/lac/petsc_matrix_base.h>
#  include <deal.II/lac/petsc_precondition.h>
#  include <deal.II/lac/petsc_vector_base.h>

#  include <petscversion.h>

#  include <cmath>

DEAL_II_NAMESPACE_OPEN

namespace PETScWrappers
{
  SolverBase::SolverData::~SolverData()
  {
    destroy_krylov_solver(ksp);
  }



  SolverBase::SolverBase(SolverControl &cn, const MPI_Comm &mpi_communicator)
    : solver_control(cn)
    , mpi_communicator(mpi_communicator)
  {}



  void
  SolverBase::solve(const MatrixBase &        A,
                    VectorBase &              x,
                    const VectorBase &        b,
                    const PreconditionerBase &preconditioner)
  {
    /*
      TODO: PETSc duplicates communicators, so this does not work (you put
    MPI_COMM_SELF in, but get something other out when you ask PETSc for the
    communicator. This mainly fails due to the MatrixFree classes, that can not
    ask PETSc for a communicator. //Timo Heister
    Assert(A.get_mpi_communicator()==mpi_communicator, ExcMessage("PETSc Solver
    and Matrix need to use the same MPI_Comm."));
    Assert(x.get_mpi_communicator()==mpi_communicator, ExcMessage("PETSc Solver
    and Vector need to use the same MPI_Comm."));
    Assert(b.get_mpi_communicator()==mpi_communicator, ExcMessage("PETSc Solver
    and Vector need to use the same MPI_Comm."));
    */

    // first create a solver object if this
    // is necessary
    if (solver_data.get() == nullptr)
      {
        solver_data = std_cxx14::make_unique<SolverData>();

        PetscErrorCode ierr = KSPCreate(mpi_communicator, &solver_data->ksp);
        AssertThrow(ierr == 0, ExcPETScError(ierr));

        // let derived classes set the solver
        // type, and the preconditioning
        // object set the type of
        // preconditioner
        set_solver_type(solver_data->ksp);

        ierr = KSPSetPC(solver_data->ksp, preconditioner.get_pc());
        AssertThrow(ierr == 0, ExcPETScError(ierr));

        // make sure the preconditioner has an associated matrix set
        const Mat B = preconditioner;
        AssertThrow(B != nullptr,
                    ExcMessage("PETSc preconditioner should have an "
                               "associated matrix set to be used in solver."));

        // setting the preconditioner overwrites the used matrices.
        // hence, we need to set the matrices after the preconditioner.
#  if DEAL_II_PETSC_VERSION_LT(3, 5, 0)
        // the last argument is irrelevant here,
        // since we use the solver only once anyway
        ierr = KSPSetOperators(solver_data->ksp,
                               A,
                               preconditioner,
                               SAME_PRECONDITIONER);
#  else
        ierr = KSPSetOperators(solver_data->ksp, A, preconditioner);
#  endif
        AssertThrow(ierr == 0, ExcPETScError(ierr));

        // then a convergence monitor
        // function. that function simply
        // checks with the solver_control
        // object we have in this object for
        // convergence
        ierr = KSPSetConvergenceTest(solver_data->ksp,
                                     &convergence_test,
                                     reinterpret_cast<void *>(&solver_control),
                                     nullptr);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
      }

    // set the command line option prefix name
    PetscErrorCode ierr =
      KSPSetOptionsPrefix(solver_data->ksp, prefix_name.c_str());
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // set the command line options provided
    // by the user to override the defaults
    ierr = KSPSetFromOptions(solver_data->ksp);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // then do the real work: set up solver
    // internal data and solve the
    // system.
    ierr = KSPSetUp(solver_data->ksp);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    ierr = KSPSolve(solver_data->ksp, b, x);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // do not destroy solver object
    //    solver_data.reset ();

    // in case of failure: throw
    // exception
    if (solver_control.last_check() != SolverControl::success)
      AssertThrow(false,
                  SolverControl::NoConvergence(solver_control.last_step(),
                                               solver_control.last_value()));
    // otherwise exit as normal
  }


  void
  SolverBase::set_prefix(const std::string &prefix)
  {
    prefix_name = prefix;
  }


  void
  SolverBase::reset()
  {
    solver_data.reset();
  }


  SolverControl &
  SolverBase::control() const
  {
    return solver_control;
  }


  int
  SolverBase::convergence_test(KSP /*ksp*/,
                               const PetscInt      iteration,
                               const PetscReal     residual_norm,
                               KSPConvergedReason *reason,
                               void *              solver_control_x)
  {
    SolverControl &solver_control =
      *reinterpret_cast<SolverControl *>(solver_control_x);

    const SolverControl::State state =
      solver_control.check(iteration, residual_norm);

    switch (state)
      {
        case ::SolverControl::iterate:
          *reason = KSP_CONVERGED_ITERATING;
          break;

        case ::SolverControl::success:
          *reason = static_cast<KSPConvergedReason>(1);
          break;

        case ::SolverControl::failure:
          if (solver_control.last_step() > solver_control.max_steps())
            *reason = KSP_DIVERGED_ITS;
          else
            *reason = KSP_DIVERGED_DTOL;
          break;

        default:
          Assert(false, ExcNotImplemented());
      }

    // return without failure
    return 0;
  }

  void
  SolverBase::initialize(const PreconditionerBase &preconditioner)
  {
    PetscErrorCode ierr;

    solver_data = std_cxx14::make_unique<SolverData>();

    ierr = KSPCreate(mpi_communicator, &solver_data->ksp);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // let derived classes set the solver
    // type, and the preconditioning
    // object set the type of
    // preconditioner
    set_solver_type(solver_data->ksp);

    ierr = KSPSetPC(solver_data->ksp, preconditioner.get_pc());
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // then a convergence monitor
    // function. that function simply
    // checks with the solver_control
    // object we have in this object for
    // convergence
    ierr = KSPSetConvergenceTest(solver_data->ksp,
                                 &convergence_test,
                                 reinterpret_cast<void *>(&solver_control),
                                 nullptr);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // set the command line options provided
    // by the user to override the defaults
    ierr = KSPSetFromOptions(solver_data->ksp);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }



  /* ---------------------- SolverRichardson ------------------------ */

  SolverRichardson::AdditionalData::AdditionalData(const double omega)
    : omega(omega)
  {}



  SolverRichardson::SolverRichardson(SolverControl &       cn,
                                     const MPI_Comm &      mpi_communicator,
                                     const AdditionalData &data)
    : SolverBase(cn, mpi_communicator)
    , additional_data(data)
  {}


  void
  SolverRichardson::set_solver_type(KSP &ksp) const
  {
    PetscErrorCode ierr = KSPSetType(ksp, KSPRICHARDSON);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // set the damping factor from the data
    ierr = KSPRichardsonSetScale(ksp, additional_data.omega);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // in the deal.II solvers, we always
    // honor the initial guess in the
    // solution vector. do so here as well:
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // Hand over the absolute
    // tolerance and the maximum
    // iteration number to the PETSc
    // convergence criterion. The
    // custom deal.II SolverControl
    // object is ignored by the PETSc
    // Richardson method (when no
    // PETSc monitoring is present),
    // since in this case PETSc
    // uses a faster version of
    // the Richardson iteration,
    // where no residual is
    // available.
    ierr = KSPSetTolerances(ksp,
                            PETSC_DEFAULT,
                            this->solver_control.tolerance(),
                            PETSC_DEFAULT,
                            this->solver_control.max_steps() + 1);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }


  /* ---------------------- SolverChebychev ------------------------ */

  SolverChebychev::SolverChebychev(SolverControl &       cn,
                                   const MPI_Comm &      mpi_communicator,
                                   const AdditionalData &data)
    : SolverBase(cn, mpi_communicator)
    , additional_data(data)
  {}


  void
  SolverChebychev::set_solver_type(KSP &ksp) const
  {
    PetscErrorCode ierr = KSPSetType(ksp, KSPCHEBYSHEV);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // in the deal.II solvers, we always
    // honor the initial guess in the
    // solution vector. do so here as well:
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }


  /* ---------------------- SolverCG ------------------------ */

  SolverCG::SolverCG(SolverControl &       cn,
                     const MPI_Comm &      mpi_communicator,
                     const AdditionalData &data)
    : SolverBase(cn, mpi_communicator)
    , additional_data(data)
  {}


  void
  SolverCG::set_solver_type(KSP &ksp) const
  {
    PetscErrorCode ierr = KSPSetType(ksp, KSPCG);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // in the deal.II solvers, we always
    // honor the initial guess in the
    // solution vector. do so here as well:
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }


  /* ---------------------- SolverBiCG ------------------------ */

  SolverBiCG::SolverBiCG(SolverControl &       cn,
                         const MPI_Comm &      mpi_communicator,
                         const AdditionalData &data)
    : SolverBase(cn, mpi_communicator)
    , additional_data(data)
  {}


  void
  SolverBiCG::set_solver_type(KSP &ksp) const
  {
    PetscErrorCode ierr = KSPSetType(ksp, KSPBICG);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // in the deal.II solvers, we always
    // honor the initial guess in the
    // solution vector. do so here as well:
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }


  /* ---------------------- SolverGMRES ------------------------ */

  SolverGMRES::AdditionalData::AdditionalData(
    const unsigned int restart_parameter,
    const bool         right_preconditioning)
    : restart_parameter(restart_parameter)
    , right_preconditioning(right_preconditioning)
  {}



  SolverGMRES::SolverGMRES(SolverControl &       cn,
                           const MPI_Comm &      mpi_communicator,
                           const AdditionalData &data)
    : SolverBase(cn, mpi_communicator)
    , additional_data(data)
  {}


  void
  SolverGMRES::set_solver_type(KSP &ksp) const
  {
    PetscErrorCode ierr = KSPSetType(ksp, KSPGMRES);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // set the restart parameter from the
    // data. we would like to use the simple
    // code that is commented out, but this
    // leads to nasty warning and error
    // messages due to some stupidity on
    // PETSc's side: KSPGMRESSetRestart is
    // implemented as a macro in which return
    // statements are hidden. This may work
    // if people strictly follow the PETSc
    // coding style of always having
    // functions return an integer error
    // code, but the present function isn't
    // like this.
    /*
        ierr = KSPGMRESSetRestart (ksp, additional_data.restart_parameter);
        AssertThrow (ierr == 0, ExcPETScError(ierr));
    */
    // so rather expand their macros by hand,
    // and do some equally nasty stuff that at
    // least doesn't yield warnings...
    int (*fun_ptr)(KSP, int);
    ierr = PetscObjectQueryFunction(reinterpret_cast<PetscObject>(ksp),
                                    "KSPGMRESSetRestart_C",
                                    reinterpret_cast<void (**)()>(&fun_ptr));
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    ierr = (*fun_ptr)(ksp, additional_data.restart_parameter);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // Set preconditioning side to
    // right
    if (additional_data.right_preconditioning)
      {
        ierr = KSPSetPCSide(ksp, PC_RIGHT);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
      }

    // in the deal.II solvers, we always
    // honor the initial guess in the
    // solution vector. do so here as well:
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }


  /* ---------------------- SolverBicgstab ------------------------ */

  SolverBicgstab::SolverBicgstab(SolverControl &       cn,
                                 const MPI_Comm &      mpi_communicator,
                                 const AdditionalData &data)
    : SolverBase(cn, mpi_communicator)
    , additional_data(data)
  {}


  void
  SolverBicgstab::set_solver_type(KSP &ksp) const
  {
    PetscErrorCode ierr = KSPSetType(ksp, KSPBCGS);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // in the deal.II solvers, we always
    // honor the initial guess in the
    // solution vector. do so here as well:
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }


  /* ---------------------- SolverCGS ------------------------ */

  SolverCGS::SolverCGS(SolverControl &       cn,
                       const MPI_Comm &      mpi_communicator,
                       const AdditionalData &data)
    : SolverBase(cn, mpi_communicator)
    , additional_data(data)
  {}


  void
  SolverCGS::set_solver_type(KSP &ksp) const
  {
    PetscErrorCode ierr = KSPSetType(ksp, KSPCGS);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // in the deal.II solvers, we always
    // honor the initial guess in the
    // solution vector. do so here as well:
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }


  /* ---------------------- SolverTFQMR ------------------------ */

  SolverTFQMR::SolverTFQMR(SolverControl &       cn,
                           const MPI_Comm &      mpi_communicator,
                           const AdditionalData &data)
    : SolverBase(cn, mpi_communicator)
    , additional_data(data)
  {}


  void
  SolverTFQMR::set_solver_type(KSP &ksp) const
  {
    PetscErrorCode ierr = KSPSetType(ksp, KSPTFQMR);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // in the deal.II solvers, we always
    // honor the initial guess in the
    // solution vector. do so here as well:
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }


  /* ---------------------- SolverTCQMR ------------------------ */

  SolverTCQMR::SolverTCQMR(SolverControl &       cn,
                           const MPI_Comm &      mpi_communicator,
                           const AdditionalData &data)
    : SolverBase(cn, mpi_communicator)
    , additional_data(data)
  {}


  void
  SolverTCQMR::set_solver_type(KSP &ksp) const
  {
    PetscErrorCode ierr = KSPSetType(ksp, KSPTCQMR);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // in the deal.II solvers, we always
    // honor the initial guess in the
    // solution vector. do so here as well:
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }


  /* ---------------------- SolverCR ------------------------ */

  SolverCR::SolverCR(SolverControl &       cn,
                     const MPI_Comm &      mpi_communicator,
                     const AdditionalData &data)
    : SolverBase(cn, mpi_communicator)
    , additional_data(data)
  {}


  void
  SolverCR::set_solver_type(KSP &ksp) const
  {
    PetscErrorCode ierr = KSPSetType(ksp, KSPCR);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // in the deal.II solvers, we always
    // honor the initial guess in the
    // solution vector. do so here as well:
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }


  /* ---------------------- SolverLSQR ------------------------ */

  SolverLSQR::SolverLSQR(SolverControl &       cn,
                         const MPI_Comm &      mpi_communicator,
                         const AdditionalData &data)
    : SolverBase(cn, mpi_communicator)
    , additional_data(data)
  {}


  void
  SolverLSQR::set_solver_type(KSP &ksp) const
  {
    PetscErrorCode ierr = KSPSetType(ksp, KSPLSQR);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // in the deal.II solvers, we always
    // honor the initial guess in the
    // solution vector. do so here as well:
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }


  /* ---------------------- SolverPreOnly ------------------------ */

  SolverPreOnly::SolverPreOnly(SolverControl &       cn,
                               const MPI_Comm &      mpi_communicator,
                               const AdditionalData &data)
    : SolverBase(cn, mpi_communicator)
    , additional_data(data)
  {}


  void
  SolverPreOnly::set_solver_type(KSP &ksp) const
  {
    PetscErrorCode ierr = KSPSetType(ksp, KSPPREONLY);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    // The KSPPREONLY solver of
    // PETSc never calls the convergence
    // monitor, which leads to failure
    // even when everything was ok.
    // Therefore the SolverControl status
    // is set to some nice values, which
    // guarantee a nice result at the end
    // of the solution process.
    solver_control.check(1, 0.0);

    // Using the PREONLY solver with
    // a nonzero initial guess leads
    // PETSc to produce some error messages.
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_FALSE);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }


  /* ---------------------- SparseDirectMUMPS------------------------ */

  SparseDirectMUMPS::SolverDataMUMPS::~SolverDataMUMPS()
  {
    destroy_krylov_solver(ksp);
    // the 'pc' object is owned by the 'ksp' object, and consequently
    // does not have to be destroyed explicitly here
  }


  SparseDirectMUMPS::SparseDirectMUMPS(SolverControl &       cn,
                                       const MPI_Comm &      mpi_communicator,
                                       const AdditionalData &data)
    : SolverBase(cn, mpi_communicator)
    , additional_data(data)
    , symmetric_mode(false)
  {}


  void
  SparseDirectMUMPS::set_solver_type(KSP &ksp) const
  {
    PetscErrorCode ierr = KSPSetType(ksp, KSPPREONLY);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    solver_control.check(1, 0.0);

    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_FALSE);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }

  void
  SparseDirectMUMPS::solve(const MatrixBase &A,
                           VectorBase &      x,
                           const VectorBase &b)
  {
#  ifdef DEAL_II_PETSC_WITH_MUMPS

    Mat F;

    PetscInt ival = 2, icntl = 7;
    PetscInt its;
    PetscReal rnorm;

    if (solver_data == nullptr)
      {
        solver_data = std_cxx14::make_unique<SolverDataMUMPS>();

        PetscErrorCode ierr = KSPCreate(mpi_communicator, &solver_data->ksp);
        AssertThrow(ierr == 0, ExcPETScError(ierr));

#    if DEAL_II_PETSC_VERSION_LT(3, 5, 0)
        ierr =
          KSPSetOperators(solver_data->ksp, A, A, DIFFERENT_NONZERO_PATTERN);
#    else
        ierr = KSPSetOperators(solver_data->ksp, A, A);
#    endif
        AssertThrow(ierr == 0, ExcPETScError(ierr));

        set_solver_type(solver_data->ksp);

        ierr = KSPGetPC(solver_data->ksp, &solver_data->pc);
        AssertThrow(ierr == 0, ExcPETScError(ierr));

        if (symmetric_mode)
          ierr = PCSetType(solver_data->pc, PCCHOLESKY);
        else
          ierr = PCSetType(solver_data->pc, PCLU);
        AssertThrow(ierr == 0, ExcPETScError(ierr));

        ierr = KSPSetConvergenceTest(solver_data->ksp,
                                     &convergence_test,
                                     reinterpret_cast<void *>(&solver_control),
                                     PETSC_NULL);
        AssertThrow(ierr == 0, ExcPETScError(ierr));

#    if DEAL_II_PETSC_VERSION_LT(3, 9, 0)
        ierr = PCFactorSetMatSolverPackage(solver_data->pc, MATSOLVERMUMPS);
#    else
        ierr = PCFactorSetMatSolverType(solver_data->pc, MATSOLVERMUMPS);
#    endif
        AssertThrow(ierr == 0, ExcPETScError(ierr));

#    if DEAL_II_PETSC_VERSION_LT(3, 9, 0)
        ierr = PCFactorSetUpMatSolverPackage(solver_data->pc);
#    else
        ierr = PCFactorSetUpMatSolverType(solver_data->pc);
#    endif
        AssertThrow(ierr == 0, ExcPETScError(ierr));

        ierr = PCFactorGetMatrix(solver_data->pc, &F);
        AssertThrow(ierr == 0, ExcPETScError(ierr));

        ierr = MatMumpsSetIcntl(F, icntl, ival);
        AssertThrow(ierr == 0, ExcPETScError(ierr));

        ierr = KSPSetOptionsPrefix(solver_data->ksp, prefix_name.c_str());
        AssertThrow(ierr == 0, ExcPETScError(ierr));

        ierr = KSPSetFromOptions(solver_data->ksp);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
      }

    PetscErrorCode ierr = KSPSolve(solver_data->ksp, b, x);
    AssertThrow(ierr == 0, ExcPETScError(ierr));

    if (solver_control.last_check() != SolverControl::success)
      {
        AssertThrow(false,
                    SolverControl::NoConvergence(solver_control.last_step(),
                                                 solver_control.last_value()));
      }
    else
      {
        ierr = KSPGetIterationNumber(solver_data->ksp, &its);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
        ierr = KSPGetResidualNorm(solver_data->ksp, &rnorm);
        AssertThrow(ierr == 0, ExcPETScError(ierr));
      }

#  else // DEAL_II_PETSC_WITH_MUMPS
    Assert(
      false,
      ExcMessage(
        "Your PETSc installation does not include a copy of "
        "the MUMPS package necessary for this solver. You will need to configure "
        "PETSc so that it includes MUMPS, recompile it, and then re-configure "
        "and recompile deal.II as well."));

    // Cast to void to silence compiler warnings
    (void)A;
    (void)x;
    (void)b;
#  endif
  }

  PetscErrorCode
  SparseDirectMUMPS::convergence_test(KSP /*ksp*/,
                                      const PetscInt      iteration,
                                      const PetscReal     residual_norm,
                                      KSPConvergedReason *reason,
                                      void *              solver_control_x)
  {
    SolverControl &solver_control =
      *reinterpret_cast<SolverControl *>(solver_control_x);

    const SolverControl::State state =
      solver_control.check(iteration, residual_norm);

    switch (state)
      {
        case ::SolverControl::iterate:
          *reason = KSP_CONVERGED_ITERATING;
          break;

        case ::SolverControl::success:
          *reason = static_cast<KSPConvergedReason>(1);
          break;

        case ::SolverControl::failure:
          if (solver_control.last_step() > solver_control.max_steps())
            *reason = KSP_DIVERGED_ITS;
          else
            *reason = KSP_DIVERGED_DTOL;
          break;

        default:
          Assert(false, ExcNotImplemented());
      }

    return 0;
  }

  void
  SparseDirectMUMPS::set_symmetric_mode(const bool flag)
  {
    symmetric_mode = flag;
  }

} // namespace PETScWrappers

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_PETSC