PETScWrappers::Solver Bring the class into the new century

diff --git a/include/deal.II/lac/petsc_compatibility.h b/include/deal.II/lac/petsc_compatibility.h
index 96bd2306294baf257445057fac1539bd52b4bb4d..46d61f78221f437f72861f02b473518dd56eb5de 100644 (file)
--- a/include/deal.II/lac/petsc_compatibility.h
+++ b/include/deal.II/lac/petsc_compatibility.h
@@ -53,25 +53,6 @@ namespace PETScWrappers
 
 
 
-  /**
-   * Destroy a Krylov Subspace (KSP) PETSc solver. This function wraps
-   * KSPDestroy with a version check (the signature of this function changed
-   * in PETSc 3.2.0).
-   *
-   * @warning Since the primary intent of this function is to enable RAII
-   * semantics in the PETSc wrappers, this function will not throw an
-   * exception if an error occurs, but instead just returns the error code
-   * given by MatDestroy.
-   */
-  inline PetscErrorCode
-  destroy_krylov_solver(KSP &krylov_solver)
-  {
-    // PETSc will check whether or not matrix is nullptr.
-    return KSPDestroy(&krylov_solver);
-  }
-
-
-
   /**
    * Set a PETSc matrix option. This function wraps MatSetOption with a
    * version check.

diff --git a/include/deal.II/lac/petsc_solver.h b/include/deal.II/lac/petsc_solver.h
index a68796670920d9798633dee61bbecf1a6094b193..eaee6fca717bb0abe6e434d76892d52f5beca387 100644 (file)
--- a/include/deal.II/lac/petsc_solver.h
+++ b/include/deal.II/lac/petsc_solver.h
@@ -21,13 +21,13 @@
 
 #ifdef DEAL_II_WITH_PETSC
 
+#  include <deal.II/base/smartpointer.h>
+
 #  include <deal.II/lac/exceptions.h>
 #  include <deal.II/lac/solver_control.h>
 
 #  include <petscksp.h>
 
-#  include <memory>
-
 #  ifdef DEAL_II_WITH_SLEPC
 #    include <deal.II/lac/slepc_spectral_transformation.h>
 #  endif
@@ -77,11 +77,15 @@ namespace PETScWrappers
    * linear solver and can be obtained from the <a
    * href="http://www.mcs.anl.gov/petsc">documentation and manual pages</a>.
    *
+   * @note Command line options relative to convergence tolerances are ignored
+   * when the solver is instantiated with a SolverControl.
+   *
    * @note Repeated calls to solve() on a solver object with a Preconditioner
    * must be used with care. The preconditioner is initialized in the first
    * call to solve() and subsequent calls reuse the solver and preconditioner
    * object. This is done for performance reasons. The solver and
-   * preconditioner can be reset by calling reset().
+   * preconditioner can be reset by calling reset() or by using the command
+   * line option "-ksp_reuse_preconditioner false".
    *
    * @ingroup PETScWrappers
    */
@@ -89,14 +93,20 @@ namespace PETScWrappers
   {
   public:
     /**
-     * Constructor.
+     * Default constructor without SolverControl. The resulting solver will
+     * use PETSc's default convergence testing routines.
+     */
+    SolverBase();
+
+    /**
+     * Constructor with a SolverControl instance.
      */
     SolverBase(SolverControl &cn);
 
     /**
      * Destructor.
      */
-    virtual ~SolverBase() = default;
+    virtual ~SolverBase();
 
     /**
      * Solve the linear system <tt>Ax=b</tt>. Depending on the information
@@ -111,7 +121,6 @@ namespace PETScWrappers
           const VectorBase &      b,
           const PreconditionBase &preconditioner);
 
-
     /**
      * Resets the contained preconditioner and solver object. See class
      * description for more details.
@@ -119,7 +128,6 @@ namespace PETScWrappers
     virtual void
     reset();
 
-
     /**
      * Sets a prefix name for the solver object. Useful when customizing the
      * PETSc KSP object with command-line options.
@@ -127,7 +135,6 @@ namespace PETScWrappers
     void
     set_prefix(const std::string &prefix);
 
-
     /**
      * Access to object that controls convergence.
      */
@@ -141,21 +148,54 @@ namespace PETScWrappers
     void
     initialize(const PreconditionBase &preconditioner);
 
+    /**
+     * Return the PETSc KSP object.
+     */
+    KSP
+    petsc_ksp();
+
+    /**
+     * Conversion operator to gain access to the underlying PETSc type. If you
+     * do this, you cut this class off some information it may need, so this
+     * conversion operator should only be used if you know what you do.
+     */
+    operator KSP() const;
+
   protected:
+    /**
+     * The PETSc KSP object.
+     */
+    KSP ksp;
+
     /**
      * Reference to the object that controls convergence of the iterative
      * solver. In fact, for these PETSc wrappers, PETSc does so itself, but we
      * copy the data from this object before starting the solution process,
      * and copy the data back into it afterwards.
      */
-    SolverControl &solver_control;
+    SmartPointer<SolverControl, SolverBase> solver_control;
+
+    /**
+     * Utility to create the KSP object and attach convergence test.
+     */
+    void
+    initialize_ksp_with_comm(const MPI_Comm &comm);
 
     /**
      * %Function that takes a Krylov Subspace Solver context object, and sets
      * the type of solver that is requested by the derived class.
      */
     virtual void
-    set_solver_type(KSP &ksp) const = 0;
+    set_solver_type(KSP &ksp) const;
+
+    /**
+     * Utility to use deal.II convergence testing.
+     *
+     * This call changes the convergence criterion when the instance of the
+     * class has a SolverControl object associated.
+     */
+    void
+    perhaps_set_convergence_test() const;
 
     /**
      * Solver prefix name to qualify options specific to the PETSc KSP object
@@ -179,41 +219,6 @@ namespace PETScWrappers
                      KSPConvergedReason *reason,
                      void *              solver_control);
 
-    /**
-     * A structure that contains the PETSc solver and preconditioner objects.
-     * This object is preserved between subsequent calls to the solver if the
-     * same preconditioner is used as in the previous solver step. This may
-     * save some computation time, if setting up a preconditioner is
-     * expensive, such as in the case of an ILU for example.
-     *
-     * The actual declaration of this class is complicated by the fact that
-     * PETSc changed its solver interface completely and incompatibly between
-     * versions 2.1.6 and 2.2.0 :-(
-     *
-     * Objects of this type are explicitly created, but are destroyed when the
-     * surrounding solver object goes out of scope, or when we assign a new
-     * value to the pointer to this object. The respective *Destroy functions
-     * are therefore written into the destructor of this object, even though
-     * the object does not have a constructor.
-     */
-    struct SolverData
-    {
-      /**
-       * Destructor
-       */
-      ~SolverData();
-
-      /**
-       * Object for Krylov subspace solvers.
-       */
-      KSP ksp;
-    };
-
-    /**
-     * Pointer to an object that stores the solver context. This is recreated
-     * in the main solver routine if necessary.
-     */
-    std::unique_ptr<SolverData> solver_data;
 
 #  ifdef DEAL_II_WITH_SLEPC
     // Make the transformation class a friend, since it needs to set the KSP
@@ -945,38 +950,6 @@ namespace PETScWrappers
     virtual void
     set_solver_type(KSP &ksp) const override;
 
-  private:
-    /**
-     * A function that is used in PETSc as a callback to check convergence. It
-     * takes the information provided from PETSc and checks it against
-     * deal.II's own SolverControl objects to see if convergence has been
-     * reached.
-     */
-    static PetscErrorCode
-    convergence_test(KSP                 ksp,
-                     const PetscInt      iteration,
-                     const PetscReal     residual_norm,
-                     KSPConvergedReason *reason,
-                     void *              solver_control);
-
-    /**
-     * A structure that contains the PETSc solver and preconditioner objects.
-     * Since the solve member function in the base is not used here, the
-     * private SolverData struct located in the base could not be used either.
-     */
-    struct SolverDataMUMPS
-    {
-      /**
-       * Destructor
-       */
-      ~SolverDataMUMPS();
-
-      KSP ksp;
-      PC  pc;
-    };
-
-    std::unique_ptr<SolverDataMUMPS> solver_data;
-
     /**
      * Flag specifies whether matrix being factorized is symmetric or not. It
      * influences the type of the used preconditioner (PCLU or PCCHOLESKY)

diff --git a/source/lac/petsc_solver.cc b/source/lac/petsc_solver.cc
index 834a3b0c9f79ad462c2c6924f6ba093f9f8b17d8..060525f3e07ae4618466d57a10082e9310cad2fc 100644 (file)
--- a/source/lac/petsc_solver.cc
+++ b/source/lac/petsc_solver.cc
@@ -26,25 +26,57 @@
 #  include <deal.II/lac/petsc_precondition.h>
 #  include <deal.II/lac/petsc_vector_base.h>
 
-#  include <petscversion.h>
-
-#  include <cmath>
-#  include <memory>
+// Shorthand notation for PETSc error codes.
+#  define AssertPETSc(code)                          \
+    do                                               \
+      {                                              \
+        PetscErrorCode ierr = (code);                \
+        AssertThrow(ierr == 0, ExcPETScError(ierr)); \
+      }                                              \
+    while (0)
 
 DEAL_II_NAMESPACE_OPEN
 
 namespace PETScWrappers
 {
-  SolverBase::SolverData::~SolverData()
+  SolverBase::SolverBase()
+    : ksp(nullptr)
+    , solver_control(nullptr)
+  {}
+
+
+  SolverBase::SolverBase(SolverControl &cn)
+    : ksp(nullptr)
+    , solver_control(&cn)
+  {}
+
+
+
+  void
+  SolverBase::set_solver_type(KSP &) const
+  {}
+
+
+
+  SolverBase::~SolverBase()
   {
-    destroy_krylov_solver(ksp);
+    AssertPETSc(KSPDestroy(&ksp));
   }
 
 
 
-  SolverBase::SolverBase(SolverControl &cn)
-    : solver_control(cn)
-  {}
+  KSP
+  SolverBase::petsc_ksp()
+  {
+    return ksp;
+  }
+
+
+
+  SolverBase::operator KSP() const
+  {
+    return ksp;
+  }
 
 
 
@@ -54,86 +86,54 @@ namespace PETScWrappers
                     const VectorBase &      b,
                     const PreconditionBase &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)
+    if (ksp == nullptr)
       {
-        solver_data = std::make_unique<SolverData>();
-
-        PetscErrorCode ierr =
-          KSPCreate(A.get_mpi_communicator(), &solver_data->ksp);
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
+        initialize_ksp_with_comm(A.get_mpi_communicator());
 
         // 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));
-
-        // setting the preconditioner overwrites the used matrices.
-        // hence, we need to set the matrices after the preconditioner.
-        Mat B;
-        ierr = PCGetOperators(preconditioner.get_pc(), nullptr, &B);
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
-        ierr = KSPSetOperators(solver_data->ksp, A, B);
-        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_solver_type(ksp);
+
+        AssertPETSc(KSPSetPC(ksp, preconditioner.get_pc()));
+
+        /*
+         * by default we set up the preconditioner only once.
+         * this can be overriden by command line.
+         */
+        AssertPETSc(KSPSetReusePreconditioner(ksp, PETSC_TRUE));
       }
 
+    // setting the preconditioner overwrites the used matrices.
+    // hence, we need to set the matrices after the preconditioner.
+    Mat B;
+    AssertPETSc(KSPGetOperators(ksp, nullptr, &B));
+    AssertPETSc(KSPSetOperators(ksp, A, B));
+
     // set the command line option prefix name
-    PetscErrorCode ierr =
-      KSPSetOptionsPrefix(solver_data->ksp, prefix_name.c_str());
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetOptionsPrefix(ksp, prefix_name.c_str()));
 
     // set the command line options provided
     // by the user to override the defaults
-    ierr = KSPSetFromOptions(solver_data->ksp);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetFromOptions(ksp));
 
     // 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));
+    AssertPETSc(KSPSetUp(ksp));
 
-    // do not destroy solver object
-    //    solver_data.reset ();
+    AssertPETSc(KSPSolve(ksp, b, x));
 
     // in case of failure: throw
     // exception
-    if (solver_control.last_check() != SolverControl::success)
+    if (solver_control &&
+        solver_control->last_check() != SolverControl::success)
       AssertThrow(false,
-                  SolverControl::NoConvergence(solver_control.last_step(),
-                                               solver_control.last_value()));
+                  SolverControl::NoConvergence(solver_control->last_step(),
+                                               solver_control->last_value()));
     // otherwise exit as normal
   }
 
@@ -148,14 +148,18 @@ namespace PETScWrappers
   void
   SolverBase::reset()
   {
-    solver_data.reset();
+    AssertPETSc(KSPDestroy(&ksp));
   }
 
 
   SolverControl &
   SolverBase::control() const
   {
-    return solver_control;
+    AssertThrow(
+      solver_control,
+      ExcMessage(
+        "You need to create the solver with a SolverControl object if you want to call the function that returns it."));
+    return *solver_control;
   }
 
 
@@ -179,7 +183,7 @@ namespace PETScWrappers
           break;
 
         case ::dealii::SolverControl::success:
-          *reason = static_cast<KSPConvergedReason>(1);
+          *reason = KSP_CONVERGED_RTOL;
           break;
 
         case ::dealii::SolverControl::failure:
@@ -200,39 +204,44 @@ namespace PETScWrappers
 
 
   void
-  SolverBase::initialize(const PreconditionBase &preconditioner)
+  SolverBase::initialize_ksp_with_comm(const MPI_Comm &comm)
+  {
+    // Create the PETSc KSP object
+    AssertPETSc(KSPCreate(comm, &ksp));
+
+    // then a convergence monitor
+    // function that simply
+    // checks with the solver_control
+    // object we have in this object for
+    // convergence
+    perhaps_set_convergence_test();
+  }
+
+
+
+  void
+  SolverBase::perhaps_set_convergence_test() const
   {
-    PetscErrorCode ierr;
+    if (ksp && solver_control)
+      AssertPETSc(
+        KSPSetConvergenceTest(ksp, &convergence_test, solver_control, nullptr));
+  }
 
-    solver_data = std::make_unique<SolverData>();
 
-    ierr = KSPCreate(preconditioner.get_mpi_communicator(), &solver_data->ksp);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+  void
+  SolverBase::initialize(const PreconditionBase &preconditioner)
+  {
+    initialize_ksp_with_comm(preconditioner.get_mpi_communicator());
 
     // 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_solver_type(ksp);
 
     // set the command line options provided
     // by the user to override the defaults
-    ierr = KSPSetFromOptions(solver_data->ksp);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetFromOptions(ksp));
   }
 
 
@@ -263,18 +272,15 @@ namespace PETScWrappers
   void
   SolverRichardson::set_solver_type(KSP &ksp) const
   {
-    PetscErrorCode ierr = KSPSetType(ksp, KSPRICHARDSON);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetType(ksp, KSPRICHARDSON));
 
     // set the damping factor from the data
-    ierr = KSPRichardsonSetScale(ksp, additional_data.omega);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPRichardsonSetScale(ksp, additional_data.omega));
 
     // 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));
+    AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
 
     // Hand over the absolute
     // tolerance and the maximum
@@ -289,12 +295,11 @@ namespace PETScWrappers
     // 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));
+    AssertPETSc(KSPSetTolerances(ksp,
+                                 PETSC_DEFAULT,
+                                 this->solver_control->tolerance(),
+                                 PETSC_DEFAULT,
+                                 this->solver_control->max_steps() + 1));
   }
 
 
@@ -317,14 +322,12 @@ namespace PETScWrappers
   void
   SolverChebychev::set_solver_type(KSP &ksp) const
   {
-    PetscErrorCode ierr = KSPSetType(ksp, KSPCHEBYSHEV);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetType(ksp, KSPCHEBYSHEV));
 
     // 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));
+    AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
   }
 
 
@@ -346,14 +349,12 @@ namespace PETScWrappers
   void
   SolverCG::set_solver_type(KSP &ksp) const
   {
-    PetscErrorCode ierr = KSPSetType(ksp, KSPCG);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetType(ksp, KSPCG));
 
     // 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));
+    AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
   }
 
 
@@ -375,14 +376,12 @@ namespace PETScWrappers
   void
   SolverBiCG::set_solver_type(KSP &ksp) const
   {
-    PetscErrorCode ierr = KSPSetType(ksp, KSPBICG);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetType(ksp, KSPBICG));
 
     // 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));
+    AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
   }
 
 
@@ -413,24 +412,20 @@ namespace PETScWrappers
   void
   SolverGMRES::set_solver_type(KSP &ksp) const
   {
-    PetscErrorCode ierr = KSPSetType(ksp, KSPGMRES);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetType(ksp, KSPGMRES));
 
-    ierr = KSPGMRESSetRestart(ksp, additional_data.restart_parameter);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPGMRESSetRestart(ksp, additional_data.restart_parameter));
 
     // Set preconditioning side to right
     if (additional_data.right_preconditioning)
       {
-        ierr = KSPSetPCSide(ksp, PC_RIGHT);
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
+        AssertPETSc(KSPSetPCSide(ksp, PC_RIGHT));
       }
 
     // 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));
+    AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
   }
 
 
@@ -452,14 +447,12 @@ namespace PETScWrappers
   void
   SolverBicgstab::set_solver_type(KSP &ksp) const
   {
-    PetscErrorCode ierr = KSPSetType(ksp, KSPBCGS);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetType(ksp, KSPBCGS));
 
     // 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));
+    AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
   }
 
 
@@ -481,14 +474,12 @@ namespace PETScWrappers
   void
   SolverCGS::set_solver_type(KSP &ksp) const
   {
-    PetscErrorCode ierr = KSPSetType(ksp, KSPCGS);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetType(ksp, KSPCGS));
 
     // 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));
+    AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
   }
 
 
@@ -510,14 +501,12 @@ namespace PETScWrappers
   void
   SolverTFQMR::set_solver_type(KSP &ksp) const
   {
-    PetscErrorCode ierr = KSPSetType(ksp, KSPTFQMR);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetType(ksp, KSPTFQMR));
 
     // 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));
+    AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
   }
 
 
@@ -539,14 +528,12 @@ namespace PETScWrappers
   void
   SolverTCQMR::set_solver_type(KSP &ksp) const
   {
-    PetscErrorCode ierr = KSPSetType(ksp, KSPTCQMR);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetType(ksp, KSPTCQMR));
 
     // 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));
+    AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
   }
 
 
@@ -568,14 +555,12 @@ namespace PETScWrappers
   void
   SolverCR::set_solver_type(KSP &ksp) const
   {
-    PetscErrorCode ierr = KSPSetType(ksp, KSPCR);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetType(ksp, KSPCR));
 
     // 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));
+    AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
   }
 
 
@@ -599,14 +584,19 @@ namespace PETScWrappers
   void
   SolverLSQR::set_solver_type(KSP &ksp) const
   {
-    PetscErrorCode ierr = KSPSetType(ksp, KSPLSQR);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetType(ksp, KSPLSQR));
 
     // 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));
+    AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
+
+    // The KSPLSQR implementation overwrites the user-defined
+    // convergence test at creation (i.e. KSPSetType) time.
+    // This is probably a bad design decision in PETSc.
+    // Anyway, here we make sure we use our own convergence
+    // test.
+    perhaps_set_convergence_test();
   }
 
 
@@ -628,8 +618,7 @@ namespace PETScWrappers
   void
   SolverPreOnly::set_solver_type(KSP &ksp) const
   {
-    PetscErrorCode ierr = KSPSetType(ksp, KSPPREONLY);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetType(ksp, KSPPREONLY));
 
     // The KSPPREONLY solver of
     // PETSc never calls the convergence
@@ -639,13 +628,12 @@ namespace PETScWrappers
     // is set to some nice values, which
     // guarantee a nice result at the end
     // of the solution process.
-    solver_control.check(1, 0.0);
+    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));
+    AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_FALSE));
   }
 
 
@@ -668,14 +656,6 @@ namespace PETScWrappers
 
 
 
-  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
-  }
-
-
   void
   SparseDirectMUMPS::set_solver_type(KSP &ksp) const
   {
@@ -684,8 +664,7 @@ namespace PETScWrappers
      * preconditioner.  Its use is due to SparseDirectMUMPS being a direct
      * (rather than iterative) solver
      */
-    PetscErrorCode ierr = KSPSetType(ksp, KSPPREONLY);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetType(ksp, KSPPREONLY));
 
     /*
      * The KSPPREONLY solver of PETSc never calls the convergence monitor,
@@ -693,14 +672,13 @@ namespace PETScWrappers
      * 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);
+    solver_control->check(1, 0.0);
 
     /*
      * Using a PREONLY solver with a nonzero initial guess leads PETSc to
      * produce some error messages.
      */
-    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_FALSE);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_FALSE));
   }
 
   void
@@ -709,154 +687,110 @@ namespace PETScWrappers
                            const VectorBase &b)
   {
 #  ifdef DEAL_II_PETSC_WITH_MUMPS
-    /*
-     * factorization matrix to be obtained from MUMPS
-     */
-    Mat F;
-
-    /*
-     * setting MUMPS integer control parameters ICNTL to be passed to
-     * MUMPS.  Setting entry 7 of MUMPS ICNTL array (of size 40) to a value
-     * of 2. This sets use of Approximate Minimum Fill (AMF)
-     */
-    PetscInt ival = 2, icntl = 7;
-    /*
-     * number of iterations to solution (should be 1) for a direct solver
-     */
-    PetscInt its;
-    /*
-     * norm of residual
-     */
-    PetscReal rnorm;
-
     /*
      * creating a solver object if this is necessary
      */
-    if (solver_data == nullptr)
+    if (ksp == nullptr)
       {
-        solver_data = std::make_unique<SolverDataMUMPS>();
+        initialize_ksp_with_comm(A.get_mpi_communicator());
 
         /*
-         * creates the default KSP context and puts it in the location
-         * solver_data->ksp
+         * setting the solver type
          */
-        PetscErrorCode ierr =
-          KSPCreate(A.get_mpi_communicator(), &solver_data->ksp);
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
+        set_solver_type(ksp);
 
         /*
          * set the matrices involved. the last argument is irrelevant here,
          * since we use the solver only once anyway
          */
-        ierr = KSPSetOperators(solver_data->ksp, A, A);
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
-
-        /*
-         * setting the solver type
-         */
-        set_solver_type(solver_data->ksp);
+        AssertPETSc(KSPSetOperators(ksp, A, A));
 
         /*
          * getting the associated preconditioner context
          */
-        ierr = KSPGetPC(solver_data->ksp, &solver_data->pc);
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
+        PC pc;
+        AssertPETSc(KSPGetPC(ksp, &pc));
 
         /*
          * build PETSc PC for particular PCLU or PCCHOLESKY preconditioner
          * depending on whether the symmetric mode has been set
          */
         if (symmetric_mode)
-          ierr = PCSetType(solver_data->pc, PCCHOLESKY);
+          AssertPETSc(PCSetType(pc, PCCHOLESKY));
         else
-          ierr = PCSetType(solver_data->pc, PCLU);
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
-
-        /*
-         * convergence monitor function that checks with the solver_control
-         * object for convergence
-         */
-        ierr = KSPSetConvergenceTest(solver_data->ksp,
-                                     &convergence_test,
-                                     reinterpret_cast<void *>(&solver_control),
-                                     nullptr);
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
+          AssertPETSc(PCSetType(pc, PCLU));
 
-        /*
-         * set the software that is to be used to perform the lu
-         * factorization here we start to see differences with the base
-         * class solve function
-         */
+          /*
+           * set the software that is to be used to perform the lu
+           * factorization here we start to see differences with the base
+           * class solve function
+           */
 #    if DEAL_II_PETSC_VERSION_LT(3, 9, 0)
-        ierr = PCFactorSetMatSolverPackage(solver_data->pc, MATSOLVERMUMPS);
+        AssertPETSc(PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS));
 #    else
-        ierr = PCFactorSetMatSolverType(solver_data->pc, MATSOLVERMUMPS);
+        AssertPETSc(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS));
 #    endif
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
 
         /*
          * set up the package to call for the factorization
          */
 #    if DEAL_II_PETSC_VERSION_LT(3, 9, 0)
-        ierr = PCFactorSetUpMatSolverPackage(solver_data->pc);
+        AssertPETSc(PCFactorSetUpMatSolverPackage(pc));
 #    else
-        ierr = PCFactorSetUpMatSolverType(solver_data->pc);
+        AssertPETSc(PCFactorSetUpMatSolverType(pc));
 #    endif
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
-
-        /*
-         * get the factored matrix F from the preconditioner context.  This
-         * routine is valid only for LU, ILU, Cholesky, and incomplete
-         * Cholesky
-         */
-        ierr = PCFactorGetMatrix(solver_data->pc, &F);
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
 
         /*
-         * Passing the control parameters to MUMPS
+         * get the factored matrix F from the preconditioner context.
          */
-        ierr = MatMumpsSetIcntl(F, icntl, ival);
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
+        Mat F;
+        AssertPETSc(PCFactorGetMatrix(pc, &F));
 
         /*
-         * set the command line option prefix name
+         * pass control parameters to MUMPS.
+         * Setting entry 7 of MUMPS ICNTL array to a value
+         * of 2. This sets use of Approximate Minimum Fill (AMF)
          */
-        ierr = KSPSetOptionsPrefix(solver_data->ksp, prefix_name.c_str());
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
+        AssertPETSc(MatMumpsSetIcntl(F, 7, 2));
 
         /*
-         * set the command line options provided by the user to override
-         * the defaults
+         * by default we set up the preconditioner only once.
+         * this can be overriden by command line.
          */
-        ierr = KSPSetFromOptions(solver_data->ksp);
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
+        AssertPETSc(KSPSetReusePreconditioner(ksp, PETSC_TRUE));
       }
 
+    /*
+     * set the matrices involved. the last argument is irrelevant here,
+     * since we use the solver only once anyway
+     */
+    AssertPETSc(KSPSetOperators(ksp, A, A));
+
+    /*
+     * set the command line option prefix name
+     */
+    AssertPETSc(KSPSetOptionsPrefix(ksp, prefix_name.c_str()));
+
+    /*
+     * set the command line options provided by the user to override
+     * the defaults
+     */
+    AssertPETSc(KSPSetFromOptions(ksp));
+
     /*
      * solve the linear system
      */
-    PetscErrorCode ierr = KSPSolve(solver_data->ksp, b, x);
-    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    AssertPETSc(KSPSolve(ksp, b, x));
 
     /*
      * in case of failure throw exception
      */
-    if (solver_control.last_check() != SolverControl::success)
+    if (solver_control &&
+        solver_control->last_check() != SolverControl::success)
       {
         AssertThrow(false,
-                    SolverControl::NoConvergence(solver_control.last_step(),
-                                                 solver_control.last_value()));
-      }
-    else
-      {
-        /*
-         * obtain convergence information. obtain the number of iterations
-         * and residual norm
-         */
-        ierr = KSPGetIterationNumber(solver_data->ksp, &its);
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
-        ierr = KSPGetResidualNorm(solver_data->ksp, &rnorm);
-        AssertThrow(ierr == 0, ExcPETScError(ierr));
+                    SolverControl::NoConvergence(solver_control->last_step(),
+                                                 solver_control->last_value()));
       }
 
 #  else // DEAL_II_PETSC_WITH_MUMPS
@@ -877,45 +811,6 @@ namespace PETScWrappers
 
 
 
-  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 ::dealii::SolverControl::iterate:
-          *reason = KSP_CONVERGED_ITERATING;
-          break;
-
-        case ::dealii::SolverControl::success:
-          *reason = static_cast<KSPConvergedReason>(1);
-          break;
-
-        case ::dealii::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)
   {

diff --git a/source/lac/slepc_spectral_transformation.cc b/source/lac/slepc_spectral_transformation.cc
index dc68c3c8c71de09541f9dd145ad76a3e3007074d..3a8cf0321be8f4964a613e6f52b4e04fc1067878 100644 (file)
--- a/source/lac/slepc_spectral_transformation.cc
+++ b/source/lac/slepc_spectral_transformation.cc
@@ -55,7 +55,7 @@ namespace SLEPcWrappers
   void
   TransformationBase::set_solver(const PETScWrappers::SolverBase &solver)
   {
-    PetscErrorCode ierr = STSetKSP(st, solver.solver_data->ksp);
+    PetscErrorCode ierr = STSetKSP(st, solver);
     AssertThrow(ierr == 0, SolverBase::ExcSLEPcError(ierr));
   }