TrilinosWrappers::SolverDirect: adjust interface

dealii · Mar 13, 2024 · 1f07687 · 1f07687
1 parent cae4071
commit 1f07687
Show file tree

Hide file tree

Showing 2 changed files with 78 additions and 7 deletions.
diff --git a/include/deal.II/lac/trilinos_solver.h b/include/deal.II/lac/trilinos_solver.h
@@ -480,7 +480,7 @@ namespace TrilinosWrappers
    *
    * @ingroup TrilinosWrappers
    */
-  class SolverDirect
+  class SolverDirect : public Subscriptor
   {
   public:
     /**
@@ -522,6 +522,11 @@ namespace TrilinosWrappers
       std::string solver_type;
     };
 
+    /**
+     * Constructor. Creates the solver without solver control object.
+     */
+    explicit SolverDirect(const AdditionalData &data = AdditionalData());
+
     /**
      * Constructor. Takes the solver control object and creates the solver.
      */
@@ -542,23 +547,50 @@ namespace TrilinosWrappers
     void
     initialize(const SparseMatrix &A);
 
+    /**
+     * Initializes the direct solver for the matrix <tt>A</tt> and creates a
+     * factorization for it with the package chosen from the additional
+     * data structure. Note that there is no need for a preconditioner
+     * here and solve() is not called. Furthermore, @p data replaces the
+     * data stored in this instance.
+     */
+    void
+    initialize(const SparseMatrix &A, const AdditionalData &data);
+
     /**
      * Solve the linear system <tt>Ax=b</tt> based on the
-     * package set in initialize(). Note the matrix is not refactorized during
-     * this call.
+     * package set in the constructor on initialize(). Note the matrix is not
+     * refactored during this call.
      */
     void
     solve(MPI::Vector &x, const MPI::Vector &b);
 
     /**
      * Solve the linear system <tt>Ax=b</tt> based on the package set in
      * initialize() for deal.II's own parallel vectors. Note the matrix is not
-     * refactorized during this call.
+     * refactored during this call.
      */
     void
     solve(dealii::LinearAlgebra::distributed::Vector<double>       &x,
           const dealii::LinearAlgebra::distributed::Vector<double> &b);
 
+    /**
+     * Solve the linear system <tt>Ax=b</tt> based on the
+     * package set in the constructor or initialize(). Note the matrix is not
+     * refactored during this call.
+     */
+    void
+    vmult(MPI::Vector &x, const MPI::Vector &b) const;
+
+    /**
+     * Solve the linear system <tt>Ax=b</tt> based on the package set in
+     * initialize() for deal.II's own parallel vectors. Note the matrix is not
+     * refactored during this call.
+     */
+    void
+    vmult(dealii::LinearAlgebra::distributed::Vector<double>       &x,
+          const dealii::LinearAlgebra::distributed::Vector<double> &b) const;
+
     /**
      * Solve the linear system <tt>Ax=b</tt>. Creates a factorization of the
      * matrix with the package chosen from the additional data structure and
@@ -613,6 +645,11 @@ namespace TrilinosWrappers
     void
     do_solve();
 
+    /**
+     * Local dummy solver control object.
+     */
+    SolverControl solver_control_own;
+
     /**
      * Reference to the object that controls convergence of the iterative
      * solver. In fact, for these Trilinos wrappers, Trilinos does so itself,
@@ -637,7 +674,7 @@ namespace TrilinosWrappers
     /**
      * Store a copy of the flags for this particular solver.
      */
-    const AdditionalData additional_data;
+    AdditionalData additional_data;
   };
 
 

diff --git a/source/lac/trilinos_solver.cc b/source/lac/trilinos_solver.cc
@@ -636,6 +636,13 @@ namespace TrilinosWrappers
 
 
 
+  SolverDirect::SolverDirect(const AdditionalData &data)
+    : solver_control(solver_control_own)
+    , additional_data(data.output_solver_details, data.solver_type)
+  {}
+
+
+
   SolverDirect::SolverDirect(SolverControl &cn, const AdditionalData &data)
     : solver_control(cn)
     , additional_data(data.output_solver_details, data.solver_type)
@@ -696,8 +703,35 @@ namespace TrilinosWrappers
   }
 
 
+
+  void
+  SolverDirect::initialize(const SparseMatrix &A, const AdditionalData &data)
+  {
+    this->additional_data = data;
+
+    this->initialize(A);
+  }
+
+
   void
   SolverDirect::solve(MPI::Vector &x, const MPI::Vector &b)
+  {
+    this->vmult(x, b);
+  }
+
+
+
+  void
+  SolverDirect::solve(
+    dealii::LinearAlgebra::distributed::Vector<double>       &x,
+    const dealii::LinearAlgebra::distributed::Vector<double> &b)
+  {
+    this->vmult(x, b);
+  }
+
+
+  void
+  SolverDirect::vmult(MPI::Vector &x, const MPI::Vector &b) const
   {
     // Assign the empty LHS vector to the Epetra_LinearProblem object
     linear_problem->SetLHS(&x.trilinos_vector());
@@ -725,9 +759,9 @@ namespace TrilinosWrappers
 
 
   void
-  SolverDirect::solve(
+  SolverDirect::vmult(
     dealii::LinearAlgebra::distributed::Vector<double>       &x,
-    const dealii::LinearAlgebra::distributed::Vector<double> &b)
+    const dealii::LinearAlgebra::distributed::Vector<double> &b) const
   {
     Epetra_Vector ep_x(View,
                        linear_problem->GetOperator()->OperatorDomainMap(),