*
* @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
{
*/
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.
*/
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.
*/
<< "uses a map that is not compatible to the one in vector "
<< arg1
<< ". Check preconditioner and matrix setup.");
+//@}
friend class SolverBase;
*
* @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.
*/
#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,
preconditioner->SetUseTranspose(false);
}
- inline
- Epetra_Operator &
- PreconditionBase::trilinos_operator () const
- {
- AssertThrow (preconditioner, ExcMessage("Trying to dereference a null pointer."));
- return (*preconditioner);
- }
-
#endif
}
}
+ 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 -------------------------- */
/* -------------------------- 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
projects / dealii.git / commitdiff