# include <petscpc.h>
+# include <functional>
+
DEAL_II_NAMESPACE_OPEN
# ifndef DOXYGEN
class MatrixBase;
class VectorBase;
- class SolverBase;
# endif
/**
explicit PreconditionBase(const MPI_Comm &mpi_communicator);
/**
- * Constructor. This constructor is deprecated.
+ * Constructor.
*
- * @deprecated
*/
- DEAL_II_DEPRECATED
PreconditionBase();
/**
/**
* Destroys the preconditioner, leaving an object like just after having
- * called the constructor.
+ * called the default constructor.
*/
void
clear();
void
Tvmult(VectorBase &dst, const VectorBase &src) const;
+ /**
+ * Explictly call setup. This is usually not needed since PETSc will
+ * automatically call the setup function when needed.
+ */
+ void
+ setup();
/**
* Give access to the underlying PETSc object.
/**
* Return the MPI communicator object used by this preconditioner.
*/
- MPI_Comm
+ const MPI_Comm &
get_mpi_communicator() const;
protected:
- /**
- * The communicator to be used for this preconditioner.
- */
- MPI_Comm mpi_communicator;
-
/**
* The PETSc preconditioner object
*/
PC pc;
- /**
- * A pointer to the matrix that acts as a preconditioner.
- */
- Mat matrix;
-
/**
* Internal function to create the PETSc preconditioner object. Fails if
* called twice.
*/
void
- create_pc();
+ create_pc_with_mat(const MatrixBase &);
/**
- * Conversion operator to get a representation of the matrix that
- * represents this preconditioner. We use this inside the actual solver,
- * where we need to pass this matrix to the PETSc solvers.
+ * Internal function to create the PETSc preconditioner object.
*/
- operator Mat() const;
-
- // Make the solver class a friend, since it needs to call the conversion
- // operator.
- friend class SolverBase;
+ void
+ create_pc_with_comm(const MPI_Comm &);
};
initialize();
};
+ class PreconditionShell : public PreconditionBase
+ {
+ public:
+ /**
+ * Empty Constructor. You need to call initialize() before using this
+ * object.
+ */
+ PreconditionShell() = default;
+
+ /**
+ * Constructor. Take the matrix which is used to form the preconditioner.
+ */
+ PreconditionShell(const MatrixBase &matrix);
+
+ /**
+ * Same as above but without setting a matrix to form the preconditioner.
+ */
+ PreconditionShell(const MPI_Comm &communicator);
+
+
+ /**
+ * The callback for the application of the preconditioner. Defaults to
+ * a copy operation if not provided.
+ */
+ std::function<int(VectorBase &dst, const VectorBase &src)> apply;
+
+ /**
+ * The callback for the application of the transposed preconditioner.
+ * Defaults to the non-transpose operation if not provided.
+ */
+ std::function<int(VectorBase &dst, const VectorBase &src)> applyT;
+
+ protected:
+ /**
+ * Initialize the preconditioner object without knowing a particular
+ * matrix. This function sets up the PCSHELL preconditioner
+ */
+ void
+ initialize(const MPI_Comm &comm);
+
+ /**
+ * Initialize the preconditioner object with a particular
+ * matrix. This function sets up the PCSHELL preconditioner
+ */
+ void
+ initialize(const MatrixBase &matrix);
+
+ private:
+ /**
+ * Callback-function invoked by PCApply
+ */
+ static int
+ pcapply(PC pc, Vec src, Vec dst);
+
+ /**
+ * Callback-function invoked by PCApplyTranspose
+ */
+ static int
+ pcapply_transpose(PC pc, Vec src, Vec dst);
+ };
+
/**
* Alias for backwards-compatibility.
* @deprecated Use PETScWrappers::PreconditionBase instead.
namespace PETScWrappers
{
PreconditionBase::PreconditionBase(const MPI_Comm &comm)
- : mpi_communicator(comm)
- , pc(nullptr)
- , matrix(nullptr)
- {}
-
-
+ : pc(nullptr)
+ {
+ create_pc_with_comm(comm);
+ }
PreconditionBase::PreconditionBase()
- : PreconditionBase(MPI_COMM_NULL)
+ : pc(nullptr)
{}
-
-
PreconditionBase::~PreconditionBase()
{
try
{}
}
-
void
PreconditionBase::clear()
{
- mpi_communicator = MPI_COMM_NULL;
-
- matrix = nullptr;
-
- if (pc != nullptr)
+ if (pc)
{
PetscErrorCode ierr = PCDestroy(&pc);
- pc = nullptr;
AssertThrow(ierr == 0, ExcPETScError(ierr));
}
}
-
void
PreconditionBase::vmult(VectorBase &dst, const VectorBase &src) const
{
AssertThrow(ierr == 0, ExcPETScError(ierr));
}
-
-
void
PreconditionBase::Tvmult(VectorBase &dst, const VectorBase &src) const
{
AssertThrow(ierr == 0, ExcPETScError(ierr));
}
+ void
+ PreconditionBase::setup()
+ {
+ AssertThrow(pc != nullptr, StandardExceptions::ExcInvalidState());
+ PetscErrorCode ierr = PCSetUp(pc);
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+ }
- MPI_Comm
+ const MPI_Comm &
PreconditionBase::get_mpi_communicator() const
{
- return mpi_communicator;
+ static MPI_Comm comm = PETSC_COMM_SELF;
+ MPI_Comm pcomm = PetscObjectComm(reinterpret_cast<PetscObject>(pc));
+ if (pcomm != MPI_COMM_NULL)
+ comm = pcomm;
+ return comm;
}
-
-
void
- PreconditionBase::create_pc()
+ PreconditionBase::create_pc_with_mat(const MatrixBase &matrix)
{
// only allow the creation of the
// preconditioner once
AssertThrow(pc == nullptr, StandardExceptions::ExcInvalidState());
- MPI_Comm comm = get_mpi_communicator();
-
- PetscErrorCode ierr = PCCreate(comm, &pc);
+ MPI_Comm comm;
+ PetscErrorCode ierr = PetscObjectGetComm(
+ reinterpret_cast<PetscObject>(static_cast<const Mat &>(matrix)), &comm);
AssertThrow(ierr == 0, ExcPETScError(ierr));
+ create_pc_with_comm(comm);
+
# if DEAL_II_PETSC_VERSION_LT(3, 5, 0)
ierr = PCSetOperators(pc, matrix, matrix, SAME_PRECONDITIONER);
# else
AssertThrow(ierr == 0, ExcPETScError(ierr));
}
+ void
+ PreconditionBase::create_pc_with_comm(const MPI_Comm &comm)
+ {
+ clear();
+ PetscErrorCode ierr = PCCreate(comm, &pc);
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+ }
const PC &
PreconditionBase::get_pc() const
}
- PreconditionBase::operator Mat() const
- {
- return matrix;
- }
-
-
/* ----------------- PreconditionJacobi -------------------- */
PreconditionJacobi::PreconditionJacobi()
- : PreconditionBase(MPI_COMM_NULL)
+ : PreconditionBase()
{}
{
clear();
- mpi_communicator = matrix_.get_mpi_communicator();
-
- matrix = static_cast<Mat>(matrix_);
additional_data = additional_data_;
- create_pc();
+ create_pc_with_mat(matrix_);
initialize();
-
- PetscErrorCode ierr = PCSetUp(pc);
- AssertThrow(ierr == 0, ExcPETScError(ierr));
}
/* ----------------- PreconditionBlockJacobi -------------------- */
PreconditionBlockJacobi::PreconditionBlockJacobi()
- : PreconditionBase(MPI_COMM_NULL)
+ : PreconditionBase()
{}
PreconditionBlockJacobi::PreconditionBlockJacobi(
{
clear();
- mpi_communicator = matrix_.get_mpi_communicator();
-
- matrix = static_cast<Mat>(matrix_);
additional_data = additional_data_;
- create_pc();
+ create_pc_with_mat(matrix_);
initialize();
-
- PetscErrorCode ierr = PCSetUp(pc);
- AssertThrow(ierr == 0, ExcPETScError(ierr));
}
/* ----------------- PreconditionSOR -------------------- */
PreconditionSOR::PreconditionSOR()
- : PreconditionBase(MPI_COMM_NULL)
+ : PreconditionBase()
{}
{
clear();
- mpi_communicator = matrix_.get_mpi_communicator();
-
- matrix = static_cast<Mat>(matrix_);
additional_data = additional_data_;
- create_pc();
+ create_pc_with_mat(matrix_);
PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
AssertThrow(ierr == 0, ExcPETScError(ierr));
ierr = PCSetFromOptions(pc);
AssertThrow(ierr == 0, ExcPETScError(ierr));
-
- ierr = PCSetUp(pc);
- AssertThrow(ierr == 0, ExcPETScError(ierr));
}
/* ----------------- PreconditionSSOR -------------------- */
PreconditionSSOR::PreconditionSSOR()
- : PreconditionBase(MPI_COMM_NULL)
+ : PreconditionBase()
{}
{
clear();
- mpi_communicator = matrix_.get_mpi_communicator();
-
- matrix = static_cast<Mat>(matrix_);
additional_data = additional_data_;
- create_pc();
+ create_pc_with_mat(matrix_);
PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
AssertThrow(ierr == 0, ExcPETScError(ierr));
ierr = PCSetFromOptions(pc);
AssertThrow(ierr == 0, ExcPETScError(ierr));
-
- ierr = PCSetUp(pc);
- AssertThrow(ierr == 0, ExcPETScError(ierr));
}
/* ----------------- PreconditionICC -------------------- */
PreconditionICC::PreconditionICC()
- : PreconditionBase(MPI_COMM_NULL)
+ : PreconditionBase()
{}
{
clear();
- mpi_communicator = matrix_.get_mpi_communicator();
-
- matrix = static_cast<Mat>(matrix_);
additional_data = additional_data_;
- create_pc();
+ create_pc_with_mat(matrix_);
PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCICC));
AssertThrow(ierr == 0, ExcPETScError(ierr));
ierr = PCSetFromOptions(pc);
AssertThrow(ierr == 0, ExcPETScError(ierr));
-
- ierr = PCSetUp(pc);
- AssertThrow(ierr == 0, ExcPETScError(ierr));
}
/* ----------------- PreconditionILU -------------------- */
PreconditionILU::PreconditionILU()
- : PreconditionBase(MPI_COMM_NULL)
+ : PreconditionBase()
{}
{
clear();
- mpi_communicator = matrix_.get_mpi_communicator();
-
- matrix = static_cast<Mat>(matrix_);
additional_data = additional_data_;
- create_pc();
+ create_pc_with_mat(matrix_);
PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCILU));
AssertThrow(ierr == 0, ExcPETScError(ierr));
ierr = PCSetFromOptions(pc);
AssertThrow(ierr == 0, ExcPETScError(ierr));
-
- ierr = PCSetUp(pc);
- AssertThrow(ierr == 0, ExcPETScError(ierr));
}
PreconditionBoomerAMG::PreconditionBoomerAMG()
- : PreconditionBase(MPI_COMM_NULL)
+ : PreconditionBase()
{}
# ifdef DEAL_II_PETSC_WITH_HYPRE
clear();
- mpi_communicator = matrix_.get_mpi_communicator();
-
- matrix = static_cast<Mat>(matrix_);
additional_data = additional_data_;
- create_pc();
+ create_pc_with_mat(matrix_);
initialize();
- PetscErrorCode ierr = PCSetUp(pc);
- AssertThrow(ierr == 0, ExcPETScError(ierr));
-
# else // DEAL_II_PETSC_WITH_HYPRE
(void)matrix_;
(void)additional_data_;
PreconditionParaSails::PreconditionParaSails()
- : PreconditionBase(MPI_COMM_NULL)
+ : PreconditionBase()
{}
{
clear();
- mpi_communicator = matrix_.get_mpi_communicator();
-
- matrix = static_cast<Mat>(matrix_);
additional_data = additional_data_;
# ifdef DEAL_II_PETSC_WITH_HYPRE
- create_pc();
+ create_pc_with_mat(matrix_);
PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
AssertThrow(ierr == 0, ExcPETScError(ierr));
ierr = PCSetFromOptions(pc);
AssertThrow(ierr == 0, ExcPETScError(ierr));
- ierr = PCSetUp(pc);
- AssertThrow(ierr == 0, ExcPETScError(ierr));
-
# else // DEAL_II_PETSC_WITH_HYPRE
(void)pc;
Assert(false,
/* ----------------- PreconditionNone ------------------------- */
PreconditionNone::PreconditionNone()
- : PreconditionBase(MPI_COMM_NULL)
+ : PreconditionBase()
{}
{
clear();
- mpi_communicator = matrix_.get_mpi_communicator();
-
- matrix = static_cast<Mat>(matrix_);
additional_data = additional_data_;
- create_pc();
+ create_pc_with_mat(matrix_);
PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCNONE));
AssertThrow(ierr == 0, ExcPETScError(ierr));
ierr = PCSetFromOptions(pc);
AssertThrow(ierr == 0, ExcPETScError(ierr));
-
- ierr = PCSetUp(pc);
- AssertThrow(ierr == 0, ExcPETScError(ierr));
}
PreconditionLU::PreconditionLU()
- : PreconditionBase(MPI_COMM_NULL)
+ : PreconditionBase()
{}
{
clear();
- mpi_communicator = matrix_.get_mpi_communicator();
-
- matrix = static_cast<Mat>(matrix_);
additional_data = additional_data_;
- create_pc();
+ create_pc_with_mat(matrix_);
PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCLU));
AssertThrow(ierr == 0, ExcPETScError(ierr));
ierr = PCSetFromOptions(pc);
AssertThrow(ierr == 0, ExcPETScError(ierr));
-
- ierr = PCSetUp(pc);
- AssertThrow(ierr == 0, ExcPETScError(ierr));
}
/* ----------------- PreconditionBDDC -------------------- */
template <int dim>
PreconditionBDDC<dim>::PreconditionBDDC()
- : PreconditionBase(MPI_COMM_NULL)
+ : PreconditionBase()
{}
// The matrix must be of IS type. We check for this to avoid the PETSc error
// in order to suggest the correct matrix reinit method.
{
- MatType current_type;
- ierr = MatGetType(matrix, ¤t_type);
+ MatType current_type;
+ Mat A, P;
+ PetscBool flg;
+
+ ierr = PCGetOperators(pc, &A, &P);
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+ ierr = PCGetUseAmat(pc, &flg);
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+
+ ierr = MatGetType(flg ? A : P, ¤t_type);
AssertThrow(ierr == 0, ExcPETScError(ierr));
AssertThrow(
strcmp(current_type, MATIS) == 0,
{
clear();
- mpi_communicator = matrix_.get_mpi_communicator();
-
- matrix = static_cast<Mat>(matrix_);
additional_data = additional_data_;
- create_pc();
+ create_pc_with_mat(matrix_);
initialize();
+ }
- PetscErrorCode ierr = PCSetUp(pc);
+ /* ----------------- PreconditionShell -------------------- */
+
+ PreconditionShell::PreconditionShell(const MatrixBase &matrix)
+ {
+ initialize(matrix);
+ }
+
+ PreconditionShell::PreconditionShell(const MPI_Comm &comm)
+ {
+ initialize(comm);
+ }
+
+ void
+ PreconditionShell::initialize(const MPI_Comm &comm)
+ {
+ PetscErrorCode ierr;
+ if (pc)
+ {
+ ierr = PCDestroy(&pc);
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+ }
+ create_pc_with_comm(comm);
+
+ ierr = PCSetType(pc, PCSHELL);
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+ ierr = PCShellSetContext(pc, static_cast<void *>(this));
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+ ierr = PCShellSetApply(pc, PreconditionShell::pcapply);
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+ ierr = PCShellSetApplyTranspose(pc, PreconditionShell::pcapply_transpose);
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+ ierr = PCShellSetName(pc, "deal.II user solve");
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+ }
+
+ void
+ PreconditionShell::initialize(const MatrixBase &matrix)
+ {
+ initialize(matrix.get_mpi_communicator());
+ PetscErrorCode ierr;
+# if DEAL_II_PETSC_VERSION_LT(3, 5, 0)
+ ierr = PCSetOperators(pc, matrix, matrix, DIFFERENT_NONZERO_PATTERN);
+# else
+ ierr = PCSetOperators(pc, matrix, matrix);
+# endif
AssertThrow(ierr == 0, ExcPETScError(ierr));
}
+ int
+ PreconditionShell::pcapply(PC ppc, Vec x, Vec y)
+ {
+ void *ctx;
+
+ PetscFunctionBeginUser;
+ PetscErrorCode ierr = PCShellGetContext(ppc, &ctx);
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+
+ VectorBase src(x);
+ VectorBase dst(y);
+ auto user = static_cast<PreconditionShell *>(ctx);
+ if (user->apply)
+ user->apply(dst, src);
+ else
+ dst.sadd(0, src);
+ PetscFunctionReturn(0);
+ }
+
+ int
+ PreconditionShell::pcapply_transpose(PC ppc, Vec x, Vec y)
+ {
+ void *ctx;
+
+ PetscFunctionBeginUser;
+ PetscErrorCode ierr = PCShellGetContext(ppc, &ctx);
+ AssertThrow(ierr == 0, ExcPETScError(ierr));
+
+ auto user = static_cast<PreconditionShell *>(ctx);
+ VectorBase src(x);
+ VectorBase dst(y);
+ if (user->applyT)
+ user->applyT(dst, src);
+ else if (user->apply)
+ user->apply(dst, src);
+ else
+ dst.sadd(0, src);
+ PetscFunctionReturn(0);
+ }
+
} // namespace PETScWrappers