]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make Trilinos Preconditioner classes compatible with LinearOperators.
authorJean-Paul Pelteret <jppelteret@gmail.com>
Mon, 9 Jan 2017 10:11:58 +0000 (11:11 +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 Preconditioner
that is later necessary for them to be compatible with LinearOperators.

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

diff --git a/doc/news/changes/minor/20170104Jean-PaulPelteret_1 b/doc/news/changes/minor/20170104Jean-PaulPelteret_1
new file mode 100644 (file)
index 0000000..7b80fcf
--- /dev/null
@@ -0,0 +1,4 @@
+New: The Trilinos preconditioner classes have been extended to be compatible
+with LinearOperators.
+<br>
+(Jean-Paul Pelteret, 2017/01/04)
index e3a6110359cd48591e45c6e0f7ab5d83de49d3ec..d056b4cfb9a9c100e82b2ba61f89c2d88d91e29e 100644 (file)
@@ -72,7 +72,8 @@ namespace TrilinosWrappers
    *
    * @ingroup TrilinosWrappers
    * @ingroup Preconditioners
-   * @author Martin Kronbichler, 2008
+   * @author Martin Kronbichler, 2008; extension for full compatibility with
+   * LinearOperator class: Jean-Paul Pelteret, 2015
    */
   class PreconditionBase : public Subscriptor
   {
@@ -112,6 +113,22 @@ namespace TrilinosWrappers
      */
     void clear ();
 
+    /**
+     * Return the MPI communicator object in use with this matrix.
+     */
+    MPI_Comm get_mpi_communicator () const;
+
+    /**
+     * Sets an internal flag so that all operations performed by the matrix,
+     * i.e., multiplications, are done in transposed order. However, this does
+     * not reshape the matrix to transposed form directly, so care should be
+     * taken when using this flag.
+     *
+     * @note Calling this function any even number of times in succession will
+     * return the object to its original state.
+     */
+    void transpose ();
+
     /**
      * Apply the preconditioner.
      */
@@ -153,14 +170,42 @@ namespace TrilinosWrappers
                          const dealii::LinearAlgebra::distributed::Vector<double> &src) const;
 
     /**
-     * Return a reference to the underlaying Trilinos Epetra_Operator. So you
-     * can use the preconditioner with unwrapped Trilinos solver.
+     * @name Access to underlying Trilinos data
+     */
+//@{
+    /**
      *
      * Calling this function from an uninitialized object will cause an
      * exception.
      */
-    Epetra_Operator &trilinos_operator() const;
+    Epetra_Operator &trilinos_operator () const;
+    //@}
+
+    /**
+     * @name Partitioners
+     */
+//@{
+
+    /**
+     * Return the partitioning of the domain space of this matrix, i.e., the
+     * partitioning of the vectors this matrix has to be multiplied with.
+     */
+    IndexSet locally_owned_domain_indices() const;
+
+    /**
+     * Return the partitioning of the range space of this matrix, i.e., the
+     * partitioning of the vectors that are result from matrix-vector
+     * products.
+     */
+    IndexSet locally_owned_range_indices() const;
+
+//@}
 
+    /**
+     * @addtogroup Exceptions
+     *
+     */
+//@{
     /**
      * Exception.
      */
@@ -170,6 +215,7 @@ namespace TrilinosWrappers
                     << "uses a map that is not compatible to the one in vector "
                     << arg1
                     << ". Check preconditioner and matrix setup.");
+//@}
 
     friend class SolverBase;
 
@@ -1798,12 +1844,34 @@ namespace TrilinosWrappers
    *
    * @ingroup TrilinosWrappers
    * @ingroup Preconditioners
-   * @author Bruno Turcksin, 2013
+   * @author Bruno Turcksin, 2013; extension for full compatibility with
+   * LinearOperator class: Jean-Paul Pelteret, 2016
    */
   class PreconditionIdentity : public PreconditionBase
   {
   public:
 
+    /**
+     * This function is only present to provide the interface of a
+     * preconditioner to be handed to a smoother.  This does nothing.
+     */
+    struct AdditionalData
+    {
+      /**
+       * Constructor.
+       */
+      AdditionalData () {}
+    };
+
+    /**
+     * The matrix argument is ignored and here just for compatibility with more
+     * complex preconditioners.
+     * @note This function must be called when this preconditioner is to be
+     * wrapped in a LinearOperator without an exemplar materix.
+     */
+    void initialize (const SparseMatrix   &matrix,
+                     const AdditionalData &additional_data = AdditionalData());
+
     /**
      * Apply the preconditioner, i.e., dst = src.
      */
@@ -1853,6 +1921,31 @@ namespace TrilinosWrappers
 
 #ifndef DOXYGEN
 
+
+  inline
+  void
+  PreconditionBase::transpose ()
+  {
+    // This only flips a flag that tells
+    // Trilinos that any vmult operation
+    // should be done with the
+    // transpose. However, the matrix
+    // structure is not reset.
+    int ierr;
+
+    if (!preconditioner->UseTranspose())
+      {
+        ierr = preconditioner->SetUseTranspose (true);
+        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+      }
+    else
+      {
+        ierr = preconditioner->SetUseTranspose (false);
+        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+      }
+  }
+
+
   inline
   void
   PreconditionBase::vmult (VectorBase       &dst,
@@ -1972,14 +2065,6 @@ namespace TrilinosWrappers
     preconditioner->SetUseTranspose(false);
   }
 
-  inline
-  Epetra_Operator &
-  PreconditionBase::trilinos_operator () const
-  {
-    AssertThrow (preconditioner, ExcMessage("Trying to dereference a null pointer."));
-    return (*preconditioner);
-  }
-
 #endif
 
 }
index 8cd96849a2cea820829b3046783e939129577220..d3e3d23973a2097d7d6f965297665690807518a3 100644 (file)
@@ -70,6 +70,33 @@ namespace TrilinosWrappers
   }
 
 
+  MPI_Comm
+  PreconditionBase::get_mpi_communicator () const
+  {
+    return communicator.GetMpiComm();
+  }
+
+
+  Epetra_Operator &
+  PreconditionBase::trilinos_operator () const
+  {
+    AssertThrow (preconditioner, ExcMessage("Trying to dereference a null pointer."));
+    return (*preconditioner);
+  }
+
+
+  IndexSet
+  PreconditionBase::locally_owned_domain_indices() const
+  {
+    return IndexSet(preconditioner->OperatorDomainMap());
+  }
+
+
+  IndexSet
+  PreconditionBase::locally_owned_range_indices() const
+  {
+    return IndexSet(preconditioner->OperatorRangeMap());
+  }
 
   /* -------------------------- PreconditionJacobi -------------------------- */
 
@@ -674,6 +701,49 @@ namespace TrilinosWrappers
 
   /* -------------------------- PreconditionIdentity --------------------- */
 
+  void
+  PreconditionIdentity::initialize (const SparseMatrix   &matrix,
+                                    const AdditionalData &)
+  {
+    // What follows just configures a dummy preconditioner that
+    // sets up the domain and range maps, as well as the communicator.
+    // It is never used as the vmult, Tvmult operations are
+    // given a custom defintion.
+    // Note: This is only required in order to wrap this
+    // preconditioner in a LinearOperator without an exemplar
+    // matrix.
+
+    // From PreconditionJacobi:
+    // release memory before reallocation
+    preconditioner.reset ();
+    preconditioner.reset (Ifpack().Create
+                          ("point relaxation",
+                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
+                           0));
+
+    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
+                                    (preconditioner.get());
+    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
+                                     "preconditioner"));
+
+    int ierr;
+
+    Teuchos::ParameterList parameter_list;
+    parameter_list.set ("relaxation: sweeps", 1);
+    parameter_list.set ("relaxation: type", "Jacobi");
+    parameter_list.set ("relaxation: damping factor", 1.0);
+    parameter_list.set ("relaxation: min diagonal value", 0.0);
+
+    ierr = ifpack->SetParameters(parameter_list);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    ierr = ifpack->Initialize();
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    ierr = ifpack->Compute();
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
   void
   PreconditionIdentity::vmult(VectorBase       &dst,
                               const VectorBase &src) const

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.