// classes. In particular, we will need
// interfaces to the matrix and vector
// classes based on Trilinos as well as
- // generic preconditioners and the Trilinos
- // Algebraic Multigrid (AMG) preconditioner
- // that we will use for the $A$ block of the
- // Stokes matrix:
+ // Trilinos preconditioners:
#include <lac/trilinos_block_vector.h>
#include <lac/trilinos_sparse_matrix.h>
#include <lac/trilinos_block_sparse_matrix.h>
#include <lac/trilinos_precondition.h>
-#include <lac/trilinos_precondition_amg.h>
// Finally, here are two C++ headers that
// haven't been included yet by one of the
DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components,
null_space);
Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0,0),
- true, true, 5e-2, null_space, false);
+ TrilinosWrappers::PreconditionAMG::AdditionalData
+ (true, true, 5e-2, null_space, 0, false));
// TODO: we could throw away the (0,0)
// block here since things have been
#include <lac/trilinos_sparse_matrix.h>
#include <lac/trilinos_block_sparse_matrix.h>
#include <lac/trilinos_precondition.h>
-#include <lac/trilinos_precondition_amg.h>
#include <grid/tria.h>
#include <grid/grid_generator.h>
null_space);
Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0,0),
- true, true, 5e-2, null_space, false);
+ TrilinosWrappers::PreconditionAMG::AdditionalData
+ (true, true, 5e-2, null_space, 0, false));
Mp_preconditioner = boost::shared_ptr<TrilinosWrappers::PreconditionIC>
(new TrilinosWrappers::PreconditionIC());
#include <base/config.h>
#include <base/subscriptor.h>
+#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
+
+#include <boost/shared_ptr.hpp>
#ifdef DEAL_II_USE_TRILINOS
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+# include <Epetra_MpiComm.h>
+#else
+# include <Epetra_SerialComm.h>
+#endif
+#include <Epetra_Map.h>
+
#include <Teuchos_RefCountPtr.hpp>
+#include <Epetra_Operator.h>
// forward declarations
class Ifpack_Preconditioner;
+namespace ML_Epetra
+{
+ class MultiLevelPreconditioner;
+}
+
DEAL_II_NAMESPACE_OPEN
* sparse matrix.
*/
PreconditionBase ();
+
+ /**
+ * Constructor. Does not do
+ * anything. The
+ * <tt>initialize</tt> function
+ * of the derived classes will
+ * have to create the
+ * preconditioner from a given
+ * sparse matrix.
+ */
+ ~PreconditionBase ();
/**
* Apply the preconditioner.
void vmult (VectorBase &dst,
const VectorBase &src) const;
+ /**
+ * Apply the preconditioner on
+ * deal.II data structures
+ * instead of the ones provided
+ * in the Trilinos wrapper
+ * class.
+ */
+ void vmult (dealii::Vector<double> &dst,
+ const dealii::Vector<double> &src) const;
+
/**
* Exception.
*/
friend class SolverBase;
friend class SolverSaddlePoint;
- //protected:
+ protected:
/**
* This is a pointer to the
- * preconditioner object.
+ * preconditioner object that
+ * is used when applying the
+ * preconditioner.
+ */
+ Teuchos::RCP<const Epetra_Operator> preconditioner;
+
+ /**
+ * Internal communication
+ * pattern in case the matrix
+ * needs to be copied from
+ * deal.II format.
*/
- Teuchos::RefCountPtr<Ifpack_Preconditioner> preconditioner;
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ Epetra_MpiComm communicator;
+#else
+ Epetra_SerialComm communicator;
+#endif
+ /**
+ * Internal Trilinos map in
+ * case the matrix needs to be
+ * copied from deal.II format.
+ */
+ std::auto_ptr<Epetra_Map> map;
};
const double min_diagonal = 0);
/**
- * Relaxation parameter and
- * minimal diagonal value.
+ * This specifies the
+ * relaxation parameter in the
+ * Jacobi preconditioner.
*/
- double omega, min_diagonal;
+ double omega;
+
+ /**
+ * This specifies the minimum
+ * value the diagonal elements
+ * should have. This might be
+ * necessary when the Jacobi
+ * preconditioner is used on
+ * matrices with zero diagonal
+ * elements. In that case, a
+ * straight-forward application
+ * would not be possible since
+ * we would divide by zero.
+ */
+ double min_diagonal;
};
/**
*/
void initialize (const SparseMatrix &matrix,
const AdditionalData &additional_data = AdditionalData());
+
+ private:
+ /**
+ * This is a pointer to the
+ * Ifpack data contained in
+ * this preconditioner.
+ */
+ Teuchos::RCP<Ifpack_Preconditioner> ifpack;
};
const unsigned int overlap = 0);
/**
- * Relaxation parameter,
- * minimal diagonal element,
- * and overlap.
+ * This specifies the (over-)
+ * relaxation parameter in the
+ * SSOR preconditioner.
+ */
+ double omega;
+
+ /**
+ * This specifies the minimum
+ * value the diagonal elements
+ * should have. This might be
+ * necessary when the SSOR
+ * preconditioner is used on
+ * matrices with zero diagonal
+ * elements. In that case, a
+ * straight-forward application
+ * would not be possible since
+ * we divide by the diagonal
+ * element.
+ */
+ double min_diagonal;
+
+ /**
+ * This determines how large
+ * the overlap of the local
+ * matrix portions on each
+ * processor in a parallel
+ * application should be.
*/
- double omega, min_diagonal;
unsigned int overlap;
};
*/
void initialize (const SparseMatrix &matrix,
const AdditionalData &additional_data = AdditionalData());
+
+ private:
+ /**
+ * This is a pointer to the
+ * Ifpack data contained in
+ * this preconditioner.
+ */
+ Teuchos::RCP<Ifpack_Preconditioner> ifpack;
};
AdditionalData (const double omega = 1,
const double min_diagonal = 0,
const unsigned int overlap = 0);
-
+
/**
- * Relaxation parameter,
- * minimal diagonal element,
- * and overlap.
+ * This specifies the (over-)
+ * relaxation parameter in the
+ * SOR preconditioner.
+ */
+ double omega;
+
+ /**
+ * This specifies the minimum
+ * value the diagonal elements
+ * should have. This might be
+ * necessary when the SOR
+ * preconditioner is used on
+ * matrices with zero diagonal
+ * elements. In that case, a
+ * straight-forward application
+ * would not be possible since
+ * we divide by the diagonal
+ * element.
+ */
+ double min_diagonal;
+
+ /**
+ * This determines how large
+ * the overlap of the local
+ * matrix portions on each
+ * processor in a parallel
+ * application should be.
*/
- double omega, min_diagonal;
unsigned int overlap;
};
* are any.
*/
void initialize (const SparseMatrix &matrix,
- const AdditionalData &additional_data = AdditionalData());
+ const AdditionalData &additional_data = AdditionalData());
+
+ private:
+ /**
+ * This is a pointer to the
+ * Ifpack data contained in
+ * this preconditioner.
+ */
+ Teuchos::RCP<Ifpack_Preconditioner> ifpack;
};
const double ic_atol = 0.,
const double ic_rtol = 1.,
const unsigned int overlap = 0);
-
+
/**
- * IC parameters and overlap.
+ * This specifies the amount of
+ * additional fill-in elements
+ * besides the sparse matrix
+ * structure. When
+ * <tt>ic_fill</tt> is large,
+ * this means that many
+ * fill-ins will be added, so
+ * that the IC preconditioner
+ * comes closer to a direct
+ * sparse Cholesky
+ * decomposition. Note,
+ * however, that this will
+ * drastically increase the
+ * memory requirement,
+ * especially when the
+ * preconditioner is used in
+ * 3D.
*/
unsigned int ic_fill;
- double ic_atol, ic_rtol;
- unsigned int overlap;
+
+ /**
+ * This specifies the amount of
+ * an absolute perturbation
+ * that will be added to the
+ * diagonal of the matrix,
+ * which sometimes can help to
+ * get better preconditioners.
+ */
+ double ic_atol;
+
+ /**
+ * This specifies the factor by
+ * which the diagonal of the
+ * matrix will be scaled, which
+ * sometimes can help to get
+ * better preconditioners.
+ */
+ double ic_rtol;
+
+ /**
+ * This determines how large
+ * the overlap of the local
+ * matrix portions on each
+ * processor in a parallel
+ * application should be.
+ */
+ unsigned int overlap;
};
/**
*/
void initialize (const SparseMatrix &matrix,
const AdditionalData &additional_data = AdditionalData());
+
+ private:
+ /**
+ * This is a pointer to the
+ * Ifpack data contained in
+ * this preconditioner.
+ */
+ Teuchos::RCP<Ifpack_Preconditioner> ifpack;
};
const double ilu_atol = 0.,
const double ilu_rtol = 1.,
const unsigned int overlap = 0);
-
+
/**
- * ILU parameters and overlap.
+ * This specifies the amount of
+ * additional fill-in elements
+ * besides the sparse matrix
+ * structure. When
+ * <tt>ilu_fill</tt> is large,
+ * this means that many
+ * fill-ins will be added, so
+ * that the ILU preconditioner
+ * comes closer to a (direct)
+ * sparse LU
+ * decomposition. Note,
+ * however, that this will
+ * drastically increase the
+ * memory requirement,
+ * especially when the
+ * preconditioner is used in
+ * 3D.
*/
unsigned int ilu_fill;
- double ilu_atol, ilu_rtol;
- unsigned int overlap;
+
+ /**
+ * This specifies the amount of
+ * an absolute perturbation
+ * that will be added to the
+ * diagonal of the matrix,
+ * which sometimes can help to
+ * get better preconditioners.
+ */
+ double ilu_atol;
+
+ /**
+ * This specifies the factor by
+ * which the diagonal of the
+ * matrix will be scaled, which
+ * sometimes can help to get
+ * better preconditioners.
+ */
+ double ilu_rtol;
+
+ /**
+ * This determines how large
+ * the overlap of the local
+ * matrix portions on each
+ * processor in a parallel
+ * application should be.
+ */
+ unsigned int overlap;
};
/**
*/
void initialize (const SparseMatrix &matrix,
const AdditionalData &additional_data = AdditionalData());
+
+ private:
+ /**
+ * This is a pointer to the
+ * Ifpack data contained in
+ * this preconditioner.
+ */
+ Teuchos::RCP<Ifpack_Preconditioner> ifpack;
};
*/
AdditionalData (const unsigned int overlap = 0);
+
/**
- * Block direct parameters and overlap.
+ * This determines how large
+ * the overlap of the local
+ * matrix portions on each
+ * processor in a parallel
+ * application should be.
*/
unsigned int overlap;
};
*/
void initialize (const SparseMatrix &matrix,
const AdditionalData &additional_data = AdditionalData());
+
+ private:
+ /**
+ * This is a pointer to the
+ * Ifpack data contained in
+ * this preconditioner.
+ */
+ Teuchos::RCP<Ifpack_Preconditioner> ifpack;
};
+
+/**
+ * This class implements an algebraic multigrid (AMG) preconditioner
+ * based on the Trilinos ML implementation, which is a black-box
+ * preconditioner that works well for many PDE-based linear problems.
+ * What this class does is twofold. When the initialize() function is
+ * invoked, a ML preconditioner object is created based on the matrix
+ * that we want the preconditioner to be based on. A call of the
+ * respective <code>vmult</code> function does call the respective
+ * operation in the Trilinos package, where it is called
+ * <code>ApplyInverse</code>. Use of this class is explained in the
+ * @ref step_31 "step-31" tutorial program.
+ *
+ * There are a few pecularities in initialize(). Since the Trilinos
+ * objects we want to use are heavily dependent on Epetra objects, the
+ * fundamental construction routines for vectors and matrices in
+ * Trilinos, we do a copy of our deal.II preconditioner matrix to a
+ * Epetra matrix. This is of course not optimal, but for the time
+ * being there is no direct support for our data interface. When
+ * doing this time-consuming operation, we can still profit from the
+ * fact that some of the entries in the preconditioner matrix are zero
+ * and hence can be neglected.
+ *
+ * The implementation is able to distinguish between matrices from
+ * elliptic problems and convection dominated problems. We use the
+ * standard options provided by Trilinos ML for elliptic problems,
+ * except that we use a Chebyshev smoother instead of a symmetric
+ * Gauss-Seidel smoother. For most elliptic problems, Chebyshev
+ * provides a better damping of high frequencies (in the algebraic
+ * sense) than Gauss-Seidel (SSOR).
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2008
+ */
+ class PreconditionAMG : public PreconditionBase
+ {
+ public:
+
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, we
+ * pretend to work on elliptic
+ * problems with linear finite
+ * elements on a scalar
+ * equation.
+ */
+ AdditionalData (const bool elliptic = true,
+ const bool higher_order_elements = false,
+ const double aggregation_threshold = 1e-4,
+ const std::vector<std::vector<bool> > &null_space = std::vector<std::vector<bool> > (1),
+ const unsigned int smoother_overlap = 0,
+ const bool output_details = false);
+
+ /**
+ * Determines whether the AMG
+ * preconditioner should be
+ * optimized for elliptic
+ * problems (ML option smoothed
+ * aggregation SA, using a
+ * Chebyshev smoother) or for
+ * non-elliptic problems (ML
+ * option non-symmetric
+ * smoothed aggregation NSSA,
+ * smoother is SSOR with
+ * underrelaxation).
+ */
+ bool elliptic;
+
+ /**
+ * Determines whether the
+ * matrix that the
+ * preconditioner is built upon
+ * is generated from linear or
+ * higher-order elements.
+ */
+ bool higher_order_elements;
+
+ /**
+ * This threshold tells the AMG
+ * setup how the coarsening
+ * should be performed. In the
+ * AMG used by ML, all points
+ * that strongly couple with
+ * the tentative coarse-level
+ * point form one
+ * aggregate. The term
+ * <em>strong coupling</em> is
+ * controlled by the variable
+ * <tt>aggregation_threshold</tt>,
+ * meaning that all elements
+ * that are not smaller than
+ * <tt>aggregation_threshold</tt>
+ * times the diagonal element
+ * do couple strongly.
+ */
+ double aggregation_threshold;
+
+ /**
+ * Specifies the constant modes
+ * (near null space) of the
+ * matrix. This parameter tells
+ * AMG whether we work on a
+ * scalar equation (where the
+ * near null space only
+ * consists of ones) or on a
+ * vector-valued equation.
+ */
+ std::vector<std::vector<bool> > null_space;
+
+ /**
+ * Determines the overlap in
+ * the SSOR/Chebyshev error
+ * smoother when run in
+ * parallel.
+ */
+ unsigned int smoother_overlap;
+
+ /**
+ * If this flag is set to
+ * <tt>true</tt>, then internal
+ * information from the ML
+ * preconditioner is printed to
+ * screen. This can be useful
+ * when debugging the
+ * preconditioner.
+ */
+ bool output_details;
+ };
+
+ /**
+ * Let Trilinos compute a
+ * multilevel hierarchy for the
+ * solution of a linear system
+ * with the given matrix. The
+ * function uses the matrix
+ * format specified in
+ * TrilinosWrappers::SparseMatrix.
+ */
+ void initialize (const SparseMatrix &matrix,
+ const AdditionalData &additional_data = AdditionalData());
+
+ /**
+ * Let Trilinos compute a
+ * multilevel hierarchy for the
+ * solution of a linear system
+ * with the given matrix. This
+ * function takes a deal.ii
+ * matrix and copies the
+ * content into a Trilinos
+ * matrix, so the function can
+ * be considered rather
+ * inefficient.
+ */
+ void initialize (const ::dealii::SparseMatrix<double> &deal_ii_sparse_matrix,
+ const AdditionalData &additional_data = AdditionalData(),
+ const double drop_tolerance = 1e-13);
+
+ /**
+ * This function can be used
+ * for a faster recalculation
+ * of the preconditioner
+ * construction when the matrix
+ * entries underlying the
+ * preconditioner have changed,
+ * but the matrix sparsity
+ * pattern has remained the
+ * same. What this function
+ * does is taking the already
+ * generated coarsening
+ * structure, computing the AMG
+ * prolongation and restriction
+ * according to a smoothed
+ * aggregation strategy and
+ * then building the whole
+ * multilevel hiearchy. This
+ * function can be considerably
+ * faster than the initialize
+ * function, since the
+ * coarsening pattern is
+ * usually the most difficult
+ * thing to do when setting up
+ * the AMG ML preconditioner.
+ */
+ void reinit ();
+
+ private:
+
+ /**
+ * A pointer to the
+ * preconditioner object.
+ */
+ Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> multilevel_operator;
+
+ /**
+ * A copy of the deal.II matrix
+ * into Trilinos format.
+ */
+ boost::shared_ptr<SparseMatrix> Matrix;
+ };
+
}
/*@}*/
#include <Ifpack.h>
#include <Teuchos_ParameterList.hpp>
+#include <Epetra_Vector.h>
+#include <Epetra_MultiVector.h>
+#include <Epetra_Vector.h>
+#include <ml_include.h>
+#include <ml_MultiLevelPreconditioner.h>
DEAL_II_NAMESPACE_OPEN
{
PreconditionBase::PreconditionBase()
+ :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ communicator (MPI_COMM_WORLD)
+#endif
{}
+ PreconditionBase::~PreconditionBase()
+ {
+ preconditioner.release();
+ }
+
+
+
void
PreconditionBase::vmult (VectorBase &dst,
const VectorBase &src) const
+ // For the implementation of
+ // the <code>vmult</code>
+ // function we note that
+ // invoking a call of the
+ // Trilinos preconditioner
+ // requires us to use Epetra
+ // vectors as well. It is
+ // faster to provide a view,
+ // i.e., feed Trilinos with a
+ // pointer to the data, so we
+ // avoid copying the content
+ // of the vectors during the
+ // iteration. In the
+ // declaration of the right
+ // hand side, we need to cast
+ // the source vector (that is
+ // <code>const</code> in all
+ // deal.II calls) to
+ // non-constant value, as this
+ // is the way Trilinos wants
+ // to have them.
+ void PreconditionBase::vmult (dealii::Vector<double> &dst,
+ const dealii::Vector<double> &src) const
+ {
+
+ Epetra_Vector LHS (View, preconditioner->OperatorDomainMap(),
+ dst.begin());
+ Epetra_Vector RHS (View, preconditioner->OperatorRangeMap(),
+ const_cast<double*>(src.begin()));
+
+ const int res = preconditioner->ApplyInverse (RHS, LHS);
+
+ Assert (res == 0,
+ ExcMessage ("Trilinos AMG MultiLevel preconditioner returned "
+ "with an error!"));
+ }
+
+
+
/* -------------------------- PreconditionJacobi -------------------------- */
PreconditionJacobi::AdditionalData::
const AdditionalData &additional_data)
{
preconditioner.release();
+ ifpack.release();
- preconditioner = Teuchos::rcp (Ifpack().Create ("point relaxation", &*matrix.matrix, 0));
- Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
- "preconditioner"));
+ ifpack = Teuchos::rcp (Ifpack().Create ("point relaxation",
+ &*matrix.matrix, 0));
+
+ Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
+ "preconditioner"));
int ierr;
parameter_list.set ("relaxation: min diagonal value",
additional_data.min_diagonal);
- ierr = preconditioner->SetParameters(parameter_list);
+ ierr = ifpack->SetParameters(parameter_list);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- ierr = preconditioner->Initialize();
+ ierr = ifpack->Initialize();
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- ierr = preconditioner->Compute();
+ ierr = ifpack->Compute();
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ preconditioner = Teuchos::rcp (&*ifpack, false);
}
const AdditionalData &additional_data)
{
preconditioner.release();
+ ifpack.release();
- preconditioner = Teuchos::rcp (Ifpack().Create ("point relaxation",
- &*matrix.matrix,
- additional_data.overlap));
- Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
- "preconditioner"));
+ ifpack = Teuchos::rcp (Ifpack().Create ("point relaxation",
+ &*matrix.matrix,
+ additional_data.overlap));
+
+ Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
+ "preconditioner"));
int ierr;
additional_data.min_diagonal);
parameter_list.set ("schwarz: combine mode", "Add");
- ierr = preconditioner->SetParameters(parameter_list);
+ ierr = ifpack->SetParameters(parameter_list);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- ierr = preconditioner->Initialize();
+ ierr = ifpack->Initialize();
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- ierr = preconditioner->Compute();
+ ierr = ifpack->Compute();
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ preconditioner = Teuchos::rcp (&*ifpack, false);
}
const AdditionalData &additional_data)
{
preconditioner.release();
+ ifpack.release();
- preconditioner = Teuchos::rcp (Ifpack().Create ("point relaxation",
+ ifpack = Teuchos::rcp (Ifpack().Create ("point relaxation",
&*matrix.matrix,
additional_data.overlap));
- Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
- "preconditioner"));
+
+ Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
+ "preconditioner"));
int ierr;
additional_data.min_diagonal);
parameter_list.set ("schwarz: combine mode", "Add");
- ierr = preconditioner->SetParameters(parameter_list);
+ ierr = ifpack->SetParameters(parameter_list);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- ierr = preconditioner->Initialize();
+ ierr = ifpack->Initialize();
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- ierr = preconditioner->Compute();
+ ierr = ifpack->Compute();
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ preconditioner = Teuchos::rcp (&*ifpack, false);
}
const AdditionalData &additional_data)
{
preconditioner.release();
+ ifpack.release();
- preconditioner = Teuchos::rcp (Ifpack().Create ("IC", &*matrix.matrix,
- additional_data.overlap));
- Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
+ ifpack = Teuchos::rcp (Ifpack().Create ("IC", &*matrix.matrix,
+ additional_data.overlap));
+
+ Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
"preconditioner"));
int ierr;
parameter_list.set ("fact: relative threshold",additional_data.ic_rtol);
parameter_list.set ("schwarz: combine mode", "Add");
- ierr = preconditioner->SetParameters(parameter_list);
+ ierr = ifpack->SetParameters(parameter_list);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- ierr = preconditioner->Initialize();
+ ierr = ifpack->Initialize();
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- ierr = preconditioner->Compute();
+ ierr = ifpack->Compute();
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ preconditioner = Teuchos::rcp (&*ifpack, false);
}
const AdditionalData &additional_data)
{
preconditioner.release();
+ ifpack.release();
- preconditioner = Teuchos::rcp (Ifpack().Create ("ILU", &*matrix.matrix,
- additional_data.overlap));
- Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
+ ifpack = Teuchos::rcp (Ifpack().Create ("ILU", &*matrix.matrix,
+ additional_data.overlap));
+
+ Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
"preconditioner"));
int ierr;
parameter_list.set ("fact: relative threshold",additional_data.ilu_rtol);
parameter_list.set ("schwarz: combine mode", "Add");
- ierr = preconditioner->SetParameters(parameter_list);
+ ierr = ifpack->SetParameters(parameter_list);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- ierr = preconditioner->Initialize();
+ ierr = ifpack->Initialize();
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- ierr = preconditioner->Compute();
+ ierr = ifpack->Compute();
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ preconditioner = Teuchos::rcp (&*ifpack, false);
}
const AdditionalData &additional_data)
{
preconditioner.release();
+ ifpack.release();
- preconditioner = Teuchos::rcp (Ifpack().Create ("Amesos", &*matrix.matrix,
- additional_data.overlap));
- Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this "
- "preconditioner"));
+ ifpack = Teuchos::rcp (Ifpack().Create ("Amesos", &*matrix.matrix,
+ additional_data.overlap));
+ Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
+ "preconditioner"));
int ierr;
Teuchos::ParameterList parameter_list;
parameter_list.set ("schwarz: combine mode", "Add");
- ierr = preconditioner->SetParameters(parameter_list);
+ ierr = ifpack->SetParameters(parameter_list);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- ierr = preconditioner->Initialize();
+ ierr = ifpack->Initialize();
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- ierr = preconditioner->Compute();
+ ierr = ifpack->Compute();
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ preconditioner = Teuchos::rcp (&*ifpack, false);
+ }
+
+
+
+/* -------------------------- PreconditionAMG -------------------------- */
+
+ PreconditionAMG::AdditionalData::
+ AdditionalData (const bool elliptic,
+ const bool higher_order_elements,
+ const double aggregation_threshold,
+ const std::vector<std::vector<bool> > &null_space,
+ const unsigned int smoother_overlap,
+ const bool output_details)
+ :
+ elliptic (elliptic),
+ higher_order_elements (higher_order_elements),
+ aggregation_threshold (aggregation_threshold),
+ null_space (null_space),
+ smoother_overlap (smoother_overlap),
+ output_details (output_details)
+ {}
+
+
+
+ void
+ PreconditionAMG:: initialize (const SparseMatrix &matrix,
+ const AdditionalData &additional_data)
+ {
+ preconditioner.release();
+ multilevel_operator.release();
+
+ const unsigned int n_rows = matrix.m();
+ const unsigned int null_space_dimension = additional_data.null_space.size();
+
+ // Build the AMG preconditioner.
+ Teuchos::ParameterList parameter_list;
+
+ if (additional_data.elliptic == true)
+ {
+ ML_Epetra::SetDefaults("SA",parameter_list);
+ parameter_list.set("smoother: type", "Chebyshev");
+ parameter_list.set("smoother: sweeps", 4);
+ }
+ else
+ {
+ ML_Epetra::SetDefaults("NSSA",parameter_list);
+ parameter_list.set("aggregation: type", "Uncoupled");
+ parameter_list.set("aggregation: block scaling", true);
+ }
+
+ parameter_list.set("smoother: ifpack overlap",
+ (int)additional_data.smoother_overlap);
+ parameter_list.set("aggregation: threshold",
+ additional_data.aggregation_threshold);
+
+ if (additional_data.output_details)
+ parameter_list.set("ML output", 10);
+ else
+ parameter_list.set("ML output", 0);
+
+ if (additional_data.higher_order_elements)
+ parameter_list.set("aggregation: type", "MIS");
+
+ const Epetra_Map * domain_map = &(matrix.matrix->DomainMap());
+
+ Epetra_MultiVector null_space_modes (*domain_map, null_space_dimension);
+
+ if (null_space_dimension > 1)
+ {
+ Assert (n_rows == additional_data.null_space[0].size(),
+ ExcDimensionMismatch(n_rows,
+ additional_data.null_space[0].size()));
+ Assert (n_rows == (unsigned int)null_space_modes.GlobalLength(),
+ ExcDimensionMismatch(n_rows,
+ null_space_modes.GlobalLength()));
+
+ const unsigned int my_size = domain_map->NumMyElements();
+ Assert (my_size == (unsigned int)domain_map->MaxLID()+1,
+ ExcDimensionMismatch (my_size, domain_map->MaxLID()+1));
+
+ // Reshape null space as a
+ // contiguous vector of
+ // doubles so that Trilinos
+ // can read from it.
+ for (unsigned int d=0; d<null_space_dimension; ++d)
+ for (unsigned int row=0; row<my_size; ++row)
+ {
+ int global_row_id = domain_map->GID(row);
+ null_space_modes.ReplaceMyValue(row, d,
+ (double)additional_data.null_space[d][global_row_id]);
+ }
+
+ parameter_list.set("null space: type", "pre-computed");
+ parameter_list.set("null space: dimension", null_space_modes.NumVectors());
+ parameter_list.set("null space: vectors", null_space_modes.Values());
+ }
+
+ multilevel_operator = Teuchos::rcp (new ML_Epetra::MultiLevelPreconditioner(
+ *matrix.matrix, parameter_list, true));
+
+ if (additional_data.output_details)
+ multilevel_operator->PrintUnused(0);
+
+ preconditioner = Teuchos::rcp (&*multilevel_operator);
+ }
+
+
+
+ void
+ PreconditionAMG::
+ initialize (const ::dealii::SparseMatrix<double> &deal_ii_sparse_matrix,
+ const AdditionalData &additional_data,
+ const double drop_tolerance)
+ {
+ const unsigned int n_rows = deal_ii_sparse_matrix.m();
+
+ // Init Epetra Matrix, avoid
+ // storing the nonzero
+ // elements.
+
+ map.reset (new Epetra_Map(n_rows, 0, communicator));
+
+ Matrix.reset();
+ Matrix = boost::shared_ptr<SparseMatrix> (new SparseMatrix());
+
+ Matrix->reinit (*map, deal_ii_sparse_matrix, drop_tolerance);
+ Matrix->compress();
+
+ initialize (*Matrix, additional_data);
+ }
+
+
+ void PreconditionAMG::
+ reinit ()
+ {
+ multilevel_operator->ReComputePreconditioner();
}
}
#include <lac/trilinos_sparse_matrix.h>
#include <lac/trilinos_vector_base.h>
#include <lac/trilinos_precondition.h>
-#include <lac/trilinos_precondition_amg.h>
#include <cmath>
#ifdef DEAL_II_USE_TRILINOS
-#include <Ifpack.h>
-#include <ml_MultiLevelPreconditioner.h>
+#include <Epetra_Operator.h>
DEAL_II_NAMESPACE_OPEN
- void
- SolverBase::solve (const SparseMatrix &A,
- VectorBase &x,
- const VectorBase &b,
- const PreconditionBase &preconditioner)
- {
- solve (A, x, b, &*preconditioner.preconditioner);
- }
-
-
-
- void
- SolverBase::solve (const SparseMatrix &A,
- VectorBase &x,
- const VectorBase &b,
- const PreconditionAMG &preconditioner)
- {
- solve (A, x, b, &*preconditioner.multilevel_operator);
- }
-
-
-
SolverControl &
SolverBase::control() const
{
SolverBase::solve (const SparseMatrix &A,
VectorBase &x,
const VectorBase &b,
- const Epetra_Operator* preconditioner)
+ const PreconditionBase &preconditioner)
{
int ierr;
// Introduce the
// preconditioner, ...
- ierr = solver.SetPrecOperator (const_cast<Epetra_Operator*>(preconditioner));
+ ierr = solver.SetPrecOperator (const_cast<Epetra_Operator*>
+ (&*preconditioner.preconditioner));
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
// ... set some options, ...