]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make Trilinos Solver classes compatible with LinearOperators.
authorJean-Paul Pelteret <jppelteret@gmail.com>
Mon, 9 Jan 2017 10:14:39 +0000 (11:14 +0100)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Tue, 17 Jan 2017 08:58:19 +0000 (09:58 +0100)
This patch adds some core functionality to the Trilinos Solver that is
later necessary for them to be compatible with LinearOperators.
Specifically, this is two additional solve() functions that accept
Trilinos EpetraOperators as the input Matrix (with one accomodating
wrapped Trilinos preconditioners, and the other another EpetraOperator
as a preconditioner.).

doc/news/changes/minor/20170104Jean-PaulPelteret_2 [new file with mode: 0644]
include/deal.II/lac/trilinos_solver.h
source/lac/trilinos_solver.cc

diff --git a/doc/news/changes/minor/20170104Jean-PaulPelteret_2 b/doc/news/changes/minor/20170104Jean-PaulPelteret_2
new file mode 100644 (file)
index 0000000..2b2e6b0
--- /dev/null
@@ -0,0 +1,5 @@
+New: Additional solve() functions have been added to the Trilinos solver
+classes, which provide the necessary extension for them to be compatible
+with LinearOperators.
+<br>
+(Jean-Paul Pelteret, 2017/01/04)
index 44798f1a98bd537bad70d461cad7b8c3b8331fde..2ec95a4a8513034b97ff32ceef600a084ce22ec5 100644 (file)
@@ -61,7 +61,8 @@ namespace TrilinosWrappers
    * ParameterList) and is similar to the deal.II class SolverSelector.
    *
    * @ingroup TrilinosWrappers
-   * @author Martin Kronbichler, 2008, 2009
+   * @author Martin Kronbichler, 2008, 2009; extension for full compatibility
+   * with LinearOperator class: Jean-Paul Pelteret, 2015
    */
   class SolverBase
   {
@@ -164,11 +165,59 @@ namespace TrilinosWrappers
      * Trilinos is chosen.
      */
     void
-    solve (Epetra_Operator        &A,
+    solve (const Epetra_Operator  &A,
            VectorBase             &x,
            const VectorBase       &b,
            const PreconditionBase &preconditioner);
 
+    /**
+     * Solve the linear system <tt>Ax=b</tt> where both <tt>A</tt> and its
+     * @p precondtioner are an operator.
+     * This function can be used when both <tt>A</tt> and the @p preconditioner
+     * are LinearOperators derived from a TrilinosPayload.
+     * Depending on the information provided by derived classes and the object
+     * passed as a preconditioner, one of the linear solvers and preconditioners
+     * of Trilinos is chosen.
+     */
+    void
+    solve (const Epetra_Operator  &A,
+           VectorBase             &x,
+           const VectorBase       &b,
+           const Epetra_Operator  &preconditioner);
+
+    /**
+     * Solve the linear system <tt>Ax=b</tt> where <tt>A</tt> is an operator,
+     * and the vectors @p x and @p b are native Trilinos vector types.
+     * This function can be used when <tt>A</tt> is a LinearOperators derived
+     * from a TrilinosPayload.
+     * Depending on the information provided by derived classes and the object
+     * passed as a preconditioner, one of the linear solvers and preconditioners
+     * of Trilinos is chosen.
+     */
+    void
+    solve (const Epetra_Operator    &A,
+           Epetra_MultiVector       &x,
+           const Epetra_MultiVector &b,
+           const PreconditionBase   &preconditioner);
+
+    /**
+     * Solve the linear system <tt>Ax=b</tt> where both <tt>A</tt> and its
+     * @p precondtioner are an operator, and the vectors @p x and @p b are
+     * native Trilinos vector types.
+     * This function can be used when both <tt>A</tt> and the @p preconditioner
+     * are LinearOperators derived from a TrilinosPayload.
+     * Depending on the information provided by derived classes and the object
+     * passed as a preconditioner, one of the linear solvers and preconditioners
+     * of Trilinos is chosen.
+     */
+    void
+    solve (const Epetra_Operator    &A,
+           Epetra_MultiVector       &x,
+           const Epetra_MultiVector &b,
+           const Epetra_Operator    &preconditioner);
+
+
+
     /**
      * Solve the linear system <tt>Ax=b</tt>. Depending on the information
      * provided by derived classes and the object passed as a preconditioner,
@@ -257,7 +306,15 @@ namespace TrilinosWrappers
      * The solve function is used to set properly the Epetra_LinearProblem,
      * once it is done this function solves the linear problem.
      */
-    void do_solve(const PreconditionBase &preconditioner);
+    template<typename Preconditioner>
+    void do_solve(const Preconditioner &preconditioner);
+
+    /**
+    * A function that sets the preconditioner that the solver will apply
+    */
+    template<typename Preconditioner>
+    void set_preconditioner (AztecOO              &solver,
+                             const Preconditioner &preconditioner);
 
     /**
      * A structure that collects the Trilinos sparse matrix, the right hand
index 3c343844f78e9dc8a66629f53f4ea8f0668e5141..400b7df1eeac79d40154264d94f5c6c69285306c 100644 (file)
@@ -88,8 +88,10 @@ namespace TrilinosWrappers
 
 
 
+  // Note: "A" is set as a constant reference so that all patterns for ::solve
+  //       can be used by the inverse_operator of LinearOperator
   void
-  SolverBase::solve (Epetra_Operator        &A,
+  SolverBase::solve (const Epetra_Operator  &A,
                      VectorBase             &x,
                      const VectorBase       &b,
                      const PreconditionBase &preconditioner)
@@ -99,7 +101,7 @@ namespace TrilinosWrappers
     // We need an Epetra_LinearProblem object to let the AztecOO solver know
     // about the matrix and vectors.
     linear_problem.reset
-    (new Epetra_LinearProblem(&A,
+    (new Epetra_LinearProblem(const_cast<Epetra_Operator *>(&A),
                               &x.trilinos_vector(),
                               const_cast<Epetra_MultiVector *>(&b.trilinos_vector())));
 
@@ -108,6 +110,72 @@ namespace TrilinosWrappers
 
 
 
+  // Note: "A" is set as a constant reference so that all patterns for ::solve
+  //       can be used by the inverse_operator of LinearOperator
+  void
+  SolverBase::solve (const Epetra_Operator  &A,
+                     VectorBase             &x,
+                     const VectorBase       &b,
+                     const Epetra_Operator  &preconditioner)
+  {
+    linear_problem.reset();
+
+    // We need an Epetra_LinearProblem object to let the AztecOO solver know
+    // about the matrix and vectors.
+    linear_problem.reset
+    (new Epetra_LinearProblem(const_cast<Epetra_Operator *>(&A),
+                              &x.trilinos_vector(),
+                              const_cast<Epetra_MultiVector *>(&b.trilinos_vector())));
+
+    do_solve(preconditioner);
+  }
+
+
+
+  // Note: "A" is set as a constant reference so that all patterns for ::solve
+  //       can be used by the inverse_operator of LinearOperator
+  void
+  SolverBase::solve (const Epetra_Operator    &A,
+                     Epetra_MultiVector       &x,
+                     const Epetra_MultiVector &b,
+                     const PreconditionBase   &preconditioner)
+  {
+    linear_problem.reset();
+
+    // We need an Epetra_LinearProblem object to let the AztecOO solver know
+    // about the matrix and vectors.
+    linear_problem.reset
+    (new Epetra_LinearProblem(const_cast<Epetra_Operator *>(&A),
+                              &x,
+                              const_cast<Epetra_MultiVector *>(&b)));
+
+    do_solve(preconditioner);
+  }
+
+
+
+  // Note: "A" is set as a constant reference so that all patterns for ::solve
+  //       can be used by the inverse_operator of LinearOperator
+  void
+  SolverBase::solve (const Epetra_Operator    &A,
+                     Epetra_MultiVector       &x,
+                     const Epetra_MultiVector &b,
+                     const Epetra_Operator    &preconditioner)
+  {
+    linear_problem.reset();
+
+    // We need an Epetra_LinearProblem object to let the AztecOO solver know
+    // about the matrix and vectors.
+    linear_problem.reset
+    (new Epetra_LinearProblem(const_cast<Epetra_Operator *>(&A),
+                              &x,
+                              const_cast<Epetra_MultiVector *>(&b)));
+
+    do_solve(preconditioner);
+  }
+
+
+
   void
   SolverBase::solve (const SparseMatrix           &A,
                      dealii::Vector<double>       &x,
@@ -214,9 +282,9 @@ namespace TrilinosWrappers
   }
 
 
-
+  template<typename Preconditioner>
   void
-  SolverBase::do_solve(const PreconditionBase &preconditioner)
+  SolverBase::do_solve(const Preconditioner &preconditioner)
   {
     int ierr;
 
@@ -246,16 +314,8 @@ namespace TrilinosWrappers
         Assert (false, ExcNotImplemented());
       }
 
-    // Introduce the preconditioner, if the identity preconditioner is used,
-    // the precondioner is set to none, ...
-    if (preconditioner.preconditioner.use_count()!=0)
-      {
-        ierr = solver.SetPrecOperator (const_cast<Epetra_Operator *>
-                                       (preconditioner.preconditioner.get()));
-        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-      }
-    else
-      solver.SetAztecOption(AZ_precond,AZ_none);
+    // Set the preconditioner
+    set_preconditioner(solver, preconditioner);
 
     // ... set some options, ...
     solver.SetAztecOption (AZ_output, additional_data.output_solver_details ?
@@ -299,6 +359,32 @@ namespace TrilinosWrappers
 
 
 
+  template<>
+  void
+  SolverBase::set_preconditioner(AztecOO                &solver,
+                                 const PreconditionBase &preconditioner)
+  {
+    // Introduce the preconditioner, if the identity preconditioner is used,
+    // the precondioner is set to none, ...
+    if (preconditioner.preconditioner.use_count()!=0)
+      {
+        const int ierr = solver.SetPrecOperator (const_cast<Epetra_Operator *>
+                                                 (preconditioner.preconditioner.get()));
+        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+      }
+    else
+      solver.SetAztecOption(AZ_precond,AZ_none);
+  }
+
+
+  template<>
+  void
+  SolverBase::set_preconditioner(AztecOO               &solver,
+                                 const Epetra_Operator &preconditioner)
+  {
+    const int ierr = solver.SetPrecOperator (const_cast<Epetra_Operator *>(&preconditioner));
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
 
 
   /* ---------------------- SolverCG ------------------------ */
@@ -637,6 +723,26 @@ namespace TrilinosWrappers
 
 }
 
+
+// explicit instantiations
+// TODO: put these instantiations into generic file
+namespace TrilinosWrappers
+{
+  template void
+  SolverBase::do_solve(const PreconditionBase &preconditioner);
+
+  template void
+  SolverBase::do_solve(const Epetra_Operator &preconditioner);
+
+  template void
+  SolverBase::set_preconditioner(AztecOO                &solver,
+                                 const PreconditionBase &preconditioner);
+
+  template void
+  SolverBase::set_preconditioner(AztecOO               &solver,
+                                 const Epetra_Operator &preconditioner);
+}
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif // DEAL_II_WITH_PETSC

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.