From 54459a823c39779511f0e22bde416ec96575185c Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Mon, 9 Jan 2017 11:11:58 +0100 Subject: [PATCH] Make Trilinos Preconditioner classes compatible with LinearOperators. This patch adds some core functionality to the Trilinos Preconditioner that is later necessary for them to be compatible with LinearOperators. --- .../changes/minor/20170104Jean-PaulPelteret_1 | 4 + include/deal.II/lac/trilinos_precondition.h | 111 ++++++++++++++++-- source/lac/trilinos_precondition.cc | 70 +++++++++++ 3 files changed, 172 insertions(+), 13 deletions(-) create mode 100644 doc/news/changes/minor/20170104Jean-PaulPelteret_1 diff --git a/doc/news/changes/minor/20170104Jean-PaulPelteret_1 b/doc/news/changes/minor/20170104Jean-PaulPelteret_1 new file mode 100644 index 0000000000..7b80fcf638 --- /dev/null +++ b/doc/news/changes/minor/20170104Jean-PaulPelteret_1 @@ -0,0 +1,4 @@ +New: The Trilinos preconditioner classes have been extended to be compatible +with LinearOperators. +
+(Jean-Paul Pelteret, 2017/01/04) diff --git a/include/deal.II/lac/trilinos_precondition.h b/include/deal.II/lac/trilinos_precondition.h index e3a6110359..d056b4cfb9 100644 --- a/include/deal.II/lac/trilinos_precondition.h +++ b/include/deal.II/lac/trilinos_precondition.h @@ -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 &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 } diff --git a/source/lac/trilinos_precondition.cc b/source/lac/trilinos_precondition.cc index 8cd96849a2..d3e3d23973 100644 --- a/source/lac/trilinos_precondition.cc +++ b/source/lac/trilinos_precondition.cc @@ -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(&matrix.trilinos_matrix()), + 0)); + + Ifpack_Preconditioner *ifpack = static_cast + (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 -- 2.39.5