// forward declarations
class Ifpack_Preconditioner;
-
DEAL_II_NAMESPACE_OPEN
/*! @addtogroup TrilinosWrappers
class VectorBase;
class SparseMatrix;
-
-/**
+ class SolverBase;
+
+/**
+ * The base class for all preconditioners based on Trilinos sparse
+ * matrices.
+ *
* @ingroup TrilinosWrappers
* @ingroup Preconditioners
* @author Martin Kronbichler, 2008
*/
- class PreconditionJacobi : public Subscriptor
+ class PreconditionBase : public Subscriptor
+ {
+ public:
+
+ /**
+ * Standardized data struct to
+ * pipe additional flags to the
+ * preconditioner.
+ */
+ struct AdditionalData
+ {};
+
+ /**
+ * 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;
+
+ /**
+ * Exception.
+ */
+ DeclException1 (ExcNonMatchingMaps,
+ std::string,
+ << "The sparse matrix that the preconditioner is based "
+ << "on a map that is not compatible to the one in vector "
+ << arg1
+ << ". Check preconditioner and matrix setup.");
+
+ friend class SolverBase;
+
+ protected:
+ /**
+ * This is a pointer to the
+ * preconditioner object.
+ */
+ Teuchos::RefCountPtr<Ifpack_Preconditioner> preconditioner;
+
+ };
+
+
+/**
+ * A wrapper class for a (pointwise) Jacobi preconditioner for
+ * Trilinos matrices. This preconditioner works both in serial and in
+ * parallel, depending on the matrix it is based on.
+ *
+ * The AdditionalData data structure allows to set preconditioner
+ * options. For the Jacobi preconditioner, these options are the
+ * damping parameter <tt>omega</tt> and a <tt>min_diagonal</tt>
+ * argument that can be used to make the preconditioner work even if
+ * the matrix contains some zero elements on the diagonal. The default
+ * settings are 1 for the damping parameter and zero for the diagonal
+ * augmentation.
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2008
+ */
+ class PreconditionJacobi : public PreconditionBase
{
public:
};
/**
- * Constructor. Does not do
- * anything. The
- * <tt>initialize</tt> function
- * will have to set create the
- * partitioner.
- */
- PreconditionJacobi ();
-
- /**
- * Take the matrix which is
- * used to form the
- * preconditioner, and
- * additional flags if there
- * are any.
+ * Take the sparse matrix the
+ * preconditioner object should
+ * be built of, and additional
+ * flags (damping parameter,
+ * etc.) if there are any.
*/
void initialize (const SparseMatrix &matrix,
const AdditionalData &additional_data = AdditionalData());
-
- /**
- * Apply the preconditioner.
- */
- void vmult (VectorBase &dst,
- const VectorBase &src) const;
-
- private:
- Teuchos::RefCountPtr<Ifpack_Preconditioner> preconditioner;
- //std::auto_ptr<Ifpack_Preconditioner> preconditioner;
- //boost::shared_ptr<Ifpack_Preconditioner> preconditioner;
};
-/**
+/**
+ * A wrapper class for a (pointwise) SSOR preconditioner for Trilinos
+ * matrices. This preconditioner works both in serial and in parallel,
+ * depending on the matrix it is based on.
+ *
+ * The AdditionalData data structure allows to set preconditioner
+ * options. For the SSOR preconditioner, these options are the
+ * damping/relaxation parameter <tt>omega</tt>, a
+ * <tt>min_diagonal</tt> argument that can be used to make the
+ * preconditioner work even if the matrix contains some zero elements
+ * on the diagonal, and a parameter <tt>overlap</tt> that determines
+ * if and how much overlap there should be between the matrix
+ * partitions on the various MPI processes. The default settings are 1
+ * for the relaxation parameter, 0 for the diagonal augmentation and 0
+ * for the overlap.
+ *
+ * Note that a parallel applicatoin of the SSOR preconditioner is
+ * actually a block-Jacobi preconditioner with block size equal to the
+ * local matrix size. Spoken more technically, this parallel operation
+ * is an <a
+ * href="http://en.wikipedia.org/wiki/Additive_Schwarz_method">additive
+ * Schwarz method</a> with an SSOR <em>approximate solve</em> as inner
+ * solver, based on the outer parallel partitioning.
+ *
* @ingroup TrilinosWrappers
* @ingroup Preconditioners
* @author Wolfgang Bangerth, 2008
*/
- class PreconditionSSOR : public Subscriptor
+ class PreconditionSSOR : public PreconditionBase
{
public:
};
/**
- * Constructor. Does not do
- * anything. The
- * <tt>initialize</tt> function
- * will have to set create the
- * partitioner.
- */
- PreconditionSSOR ();
-
- /**
- * Take the matrix which is
- * used to form the
- * preconditioner, and
- * additional flags if there
+ * Take the sparse matrix the
+ * preconditioner object should
+ * be built of, and additional
+ * flags (damping parameter,
+ * overlap in parallel
+ * computations, etc.) if there
* are any.
*/
void initialize (const SparseMatrix &matrix,
const AdditionalData &additional_data = AdditionalData());
-
- /**
- * Apply the preconditioner.
- */
- void vmult (VectorBase &dst,
- const VectorBase &src) const;
-
- private:
- Teuchos::RefCountPtr<Ifpack_Preconditioner> preconditioner;
};
-/**
+/**
+ * A wrapper class for a (pointwise) SOR preconditioner for Trilinos
+ * matrices. This preconditioner works both in serial and in parallel,
+ * depending on the matrix it is based on.
+ *
+ * The AdditionalData data structure allows to set preconditioner
+ * options. For the SOR preconditioner, these options are the
+ * damping/relaxation parameter <tt>omega</tt>, a
+ * <tt>min_diagonal</tt> argument that can be used to make the
+ * preconditioner work even if the matrix contains some zero elements
+ * on the diagonal, and a parameter <tt>overlap</tt> that determines
+ * if and how much overlap there should be between the matrix
+ * partitions on the various MPI processes. The default settings are 1
+ * for the relaxation parameter, 0 for the diagonal augmentation and 0
+ * for the overlap.
+ *
+ * Note that a parallel applicatoin of the SOR preconditioner is
+ * actually a block-Jacobi preconditioner with block size equal to the
+ * local matrix size. Spoken more technically, this parallel operation
+ * is an <a
+ * href="http://en.wikipedia.org/wiki/Additive_Schwarz_method">additive
+ * Schwarz method</a> with an SOR <em>approximate solve</em> as inner
+ * solver, based on the outer parallel partitioning.
+ *
* @ingroup TrilinosWrappers
* @ingroup Preconditioners
* @author Martin Kronbichler, 2008
*/
- class PreconditionSOR : public Subscriptor
+ class PreconditionSOR : public PreconditionBase
{
public:
};
/**
- * Constructor. Does not do
- * anything. The
- * <tt>initialize</tt> function
- * will have to set create the
- * partitioner.
- */
- PreconditionSOR ();
-
- /**
- * Take the matrix which is
- * used to form the
- * preconditioner, and
- * additional flags if there
+ * Take the sparse matrix the
+ * preconditioner object should
+ * be built of, and additional
+ * flags (damping parameter,
+ * overlap in parallel
+ * computations etc.) if there
* are any.
*/
void initialize (const SparseMatrix &matrix,
const AdditionalData &additional_data = AdditionalData());
-
- /**
- * Apply the preconditioner.
- */
- void vmult (VectorBase &dst,
- const VectorBase &src) const;
-
- private:
- Teuchos::RefCountPtr<Ifpack_Preconditioner> preconditioner;
};
-/**
+/**
+ * A wrapper class for an incomplete Cholesky factorization (IC)
+ * preconditioner for @em symmetric Trilinos matrices. This
+ * preconditioner works both in serial and in parallel, depending on
+ * the matrix it is based on. In general, an incomplete factorization
+ * does not take all fill-in elements that would appear in a full
+ * factorization (that is the basis for a direct solve). Trilinos
+ * allows to set the amount of fill-in elements, governed by the
+ * additional data argument <tt>ic_fill</tt>, so one can gradually
+ * choose between a factorization on the sparse matrix structure only
+ * (<tt>ic_fill=0</tt>) to a full factorization (<tt>ic_fill</tt> in
+ * the range of 10 to 50, depending on the spatial dimension of the
+ * PDE problem and the degree of the finite element basis functions;
+ * generally, more required fill-in elements require this parameter to
+ * be set to a higher integer value).
+ *
+ * The AdditionalData data structure allows to set preconditioner
+ * options. Besides the fill-in argument, these options are some
+ * options for perturbations (see the documentation of the
+ * AdditionalData structure for details), and a parameter
+ * <tt>overlap</tt> that determines if and how much overlap there
+ * should be between the matrix partitions on the various MPI
+ * processes. The default settings are 1 for the relaxation parameter,
+ * 0 for the diagonal augmentation and 0 for the overlap.
+ *
+ * Note that a parallel applicatoin of the IC preconditioner is
+ * actually a block-Jacobi preconditioner with block size equal to the
+ * local matrix size. Spoken more technically, this parallel operation
+ * is an <a
+ * href="http://en.wikipedia.org/wiki/Additive_Schwarz_method">additive
+ * Schwarz method</a> with an IC <em>approximate solve</em> as inner
+ * solver, based on the (outer) parallel partitioning.
+ *
* @ingroup TrilinosWrappers
* @ingroup Preconditioners
* @author Martin Kronbichler, 2008
*/
- class PreconditionIC : public Subscriptor
+ class PreconditionIC : public PreconditionBase
{
public:
/**
unsigned int overlap;
};
- /**
- * (Empty) constructor.
- */
- PreconditionIC ();
-
/**
* Initialize function. Takes
- * the matrix which is used to
- * form the preconditioner, and
- * additional flags if there
- * are any.
+ * the matrix the
+ * preconditioner should be
+ * computed of, and additional
+ * flags if there are any.
*/
void initialize (const SparseMatrix &matrix,
const AdditionalData &additional_data = AdditionalData());
-
- /**
- * Apply the preconditioner.
- */
- void vmult (VectorBase &dst,
- const VectorBase &src) const;
-
- private:
- Teuchos::RefCountPtr<Ifpack_Preconditioner> preconditioner;
};
-/**
+/**
+ * A wrapper class for an incomplete LU factorization (ILU)
+ * preconditioner for Trilinos matrices. This preconditioner works
+ * both in serial and in parallel, depending on the matrix it is based
+ * on. In general, an incomplete factorization does not take all
+ * fill-in elements that would appear in a full factorization (that is
+ * the basis for a direct solve). Trilinos allows to set the amount of
+ * fill-in elements, governed by the additional data argument
+ * <tt>ilu_fill</tt>, so one can gradually choose between a
+ * factorization on the sparse matrix structure only
+ * (<tt>ilu_fill=0</tt>) to a full factorization (<tt>ilu_fill</tt> in
+ * the range of 10 to 50, depending on the spatial dimension of the
+ * PDE problem and the degree of the finite element basis functions;
+ * generally, more required fill-in elements require this parameter to
+ * be set to a higher integer value).
+ *
+ * The AdditionalData data structure allows to set preconditioner
+ * options. Besides the fill-in argument, these options are some
+ * options for perturbations (see the documentation of the
+ * AdditionalData structure for details), and a parameter
+ * <tt>overlap</tt> that determines if and how much overlap there
+ * should be between the matrix partitions on the various MPI
+ * processes. The default settings are 1 for the relaxation parameter,
+ * 0 for the diagonal augmentation and 0 for the overlap.
+ *
+ * Note that a parallel applicatoin of the ILU preconditioner is
+ * actually a block-Jacobi preconditioner with block size equal to the
+ * local matrix size. Spoken more technically, this parallel operation
+ * is an <a
+ * href="http://en.wikipedia.org/wiki/Additive_Schwarz_method">additive
+ * Schwarz method</a> with an ILU <em>approximate solve</em> as inner
+ * solver, based on the (outer) parallel partitioning.
+ *
* @ingroup TrilinosWrappers
* @ingroup Preconditioners
* @author Martin Kronbichler, 2008
*/
- class PreconditionILU : public Subscriptor
+ class PreconditionILU : public PreconditionBase
{
public:
/**
unsigned int overlap;
};
- /**
- * (Empty) constructor.
- */
- PreconditionILU ();
-
/**
* Initialize function. Takes
* the matrix which is used to
*/
void initialize (const SparseMatrix &matrix,
const AdditionalData &additional_data = AdditionalData());
-
- /**
- * Apply the preconditioner.
- */
- void vmult (VectorBase &dst,
- const VectorBase &src) const;
-
- private:
- Teuchos::RefCountPtr<Ifpack_Preconditioner> preconditioner;
};
#else
# include <Epetra_SerialComm.h>
#endif
+#include <Teuchos_RefCountPtr.hpp>
// some forward declarations
namespace ML_Epetra
namespace TrilinosWrappers
{
-
+
+ // Forward declarations.
+ class SolverBase;
+
/**
* This class implements an algebraic multigrid (AMG) preconditioner
- * based on the Trilinos ML implementation. 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
+ * 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.
*
* 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 is
- * better than Gauss-Seidel (SSOR).
+ * 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
*/
PreconditionAMG ();
- /**
- * Destructor.
- */
- ~PreconditionAMG ();
-
/**
* Let Trilinos compute a
* multilevel hierarchy for the
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> > (),
+ const std::vector<std::vector<bool> > &null_space = std::vector<std::vector<bool> > (1),
const bool output_details = false);
/**
* the AMG ML preconditioner.
*/
void reinit ();
-
+
/**
- * Multiply the source vector
- * with the preconditioner
- * represented by this object,
- * and return it in the
- * destination vector.
+ * Apply the preconditioner.
*/
- void vmult (VectorBase &dst,
- const VectorBase &src) const;
+ void vmult (VectorBase &dst,
+ const VectorBase &src) const;
/**
* Do the same as before, but
void vmult (dealii::Vector<double> &dst,
const dealii::Vector<double> &src) const;
+ /**
+ * Exception.
+ */
+ DeclException1 (ExcNonMatchingMaps,
+ std::string,
+ << "The sparse matrix that the preconditioner is based "
+ << "on a map that is not compatible to the one in vector"
+ << arg1
+ << ". Check preconditioner and matrix setup.");
+
private:
/**
- * Pointer to the Trilinos AMG object.
+ * A pointer to the
+ * preconditioner object.
*/
- boost::shared_ptr<ML_Epetra::MultiLevelPreconditioner> multigrid_operator;
+ Teuchos::RefCountPtr<ML_Epetra::MultiLevelPreconditioner>
+ multilevel_operator;
/**
* Internal communication
* case the matrix needs to be
* copied from deal.II format.
*/
- boost::shared_ptr<Epetra_Map> Map;
+ boost::shared_ptr<Epetra_Map> Map;
/**
* A copy of the deal.II matrix
* into Trilinos format.
*/
boost::shared_ptr<SparseMatrix> Matrix;
+
+ friend class SolverBase;
};
}
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2008 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__trilinos_solver_h
+#define __deal2__trilinos_solver_h
+
+#include <base/config.h>
+#include <lac/exceptions.h>
+#include <lac/solver_control.h>
+#include <boost/shared_ptr.hpp>
+
+
+#ifdef DEAL_II_USE_TRILINOS
+
+#include <Epetra_LinearProblem.h>
+#include <AztecOO.h>
+#include <Epetra_Operator.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace TrilinosWrappers
+{
+ // forward declarations
+ class SparseMatrix;
+ class VectorBase;
+ class PreconditionBase;
+ class PreconditionAMG;
+
+
+/**
+ * Base class for solver classes using the Trilinos solvers. Since
+ * solvers in Trilinos are selected based on flags passed to a generic
+ * solver object, basically all the actual solver calls happen in this
+ * class, and derived classes simply set the right flags to select one
+ * solver or another, or to set certain parameters for individual
+ * solvers. For a general discussion on the Trilinos solver package
+ * AztecOO, we refer to the <a href =
+ * "http://trilinos.sandia.gov/packages/aztecoo/AztecOOUserGuide.pdf">AztecOO
+ * user guide</a>.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, 2008
+ */
+ class SolverBase
+ {
+ public:
+
+ /**
+ * Standardized data struct to
+ * pipe additional data to the
+ * solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Restart parameter in case
+ * the derived class is
+ * GMRES.
+ *
+ * TODO: Find a better way for
+ * doing this.
+ */
+ AdditionalData (const unsigned int gmres_restart_parameter = 30);
+
+ const unsigned int gmres_restart_parameter;
+ };
+
+ /**
+ * Constructor. Takes the
+ * solver control object and
+ * the MPI communicator over
+ * which parallel computations
+ * are to happen.
+ */
+ SolverBase (SolverControl &cn);
+
+ /**
+ * Destructor.
+ */
+ virtual ~SolverBase ();
+
+ /**
+ * Solve the linear system
+ * <tt>Ax=b</tt>. Depending on
+ * the information provided by
+ * derived classes and the
+ * object passed as a
+ * preconditioner, one of the
+ * linear solvers and
+ * preconditioners of Trilinos
+ * is chosen.
+ */
+ void
+ solve (const SparseMatrix &A,
+ VectorBase &x,
+ const VectorBase &b,
+ const PreconditionBase &preconditioner);
+
+ /**
+ * Solve the linear system
+ * <tt>Ax=b</tt>. In contrast
+ * to the other <tt>vmult</tt>
+ * function, this solves the
+ * system with an AMG
+ * preconditioner.
+ */
+ void
+ solve (const SparseMatrix &A,
+ VectorBase &x,
+ const VectorBase &b,
+ const PreconditionAMG &preconditioner);
+
+ /**
+ * Access to object that controls
+ * convergence.
+ */
+ SolverControl & control() const;
+
+ /**
+ * Exception
+ */
+ DeclException1 (ExcTrilinosError,
+ int,
+ << "An error with error number " << arg1
+ << " occured while calling a Trilinos function");
+
+ protected:
+
+ /**
+ * Reference to the object that
+ * controls convergence of the
+ * iterative solver. In fact,
+ * for these Trilinos wrappers,
+ * Trilinos does so itself, but
+ * we copy the data from this
+ * object before starting the
+ * solution process, and copy
+ * the data back into it
+ * afterwards.
+ */
+ SolverControl &solver_control;
+
+ /**
+ * String object that is set in
+ * the constructor of the
+ * derived classes and tells
+ * Trilinos which solver to
+ * use. Currently enabled
+ * options are <tt>"cg", "cgs",
+ * "gmres", "bicgstab"</tt>,
+ * and <tt>"tfqmr"</tt>.
+ */
+ enum SolverName {cg, cgs, gmres, bicgstab, tfqmr} solver_name;
+
+ private:
+
+ /**
+ * Since we might need to call
+ * the solver with different
+ * types of preconditioners
+ * (Ifpack, ML), we write an
+ * additional function that
+ * does the actual job, whereas
+ * the public <tt>vmult</tt>
+ * operations do call this
+ * function with the respective
+ * argument.
+ */
+ void solve (const SparseMatrix &A,
+ VectorBase &x,
+ const VectorBase &b,
+ const Epetra_Operator* preconditioner);
+
+ /**
+ * A structure that collects
+ * the Trilinos sparse matrix,
+ * the right hand side vector
+ * and the solution vector,
+ * which is passed down to the
+ * Trilinos solver.
+ */
+ std::auto_ptr<Epetra_LinearProblem> linear_problem;
+
+ /**
+ * A structure that contains
+ * the Trilinos solver and
+ * preconditioner objects.
+ */
+ AztecOO solver;
+
+ /**
+ * Store a copy of the flags for this
+ * particular solver.
+ */
+ const AdditionalData additional_data;
+
+ };
+
+
+
+/**
+ * An implementation of the solver interface using the Trilinos CG
+ * solver.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, 2008
+ */
+ class SolverCG : public SolverBase
+ {
+ public:
+ /**
+ * Standardized data struct to
+ * pipe additional data to the
+ * solver.
+ */
+ struct AdditionalData
+ {};
+
+ /**
+ * Constructor. In contrast to
+ * deal.II's own solvers, there is no
+ * need to give a vector memory
+ * object.
+ *
+ * The last argument takes a structure
+ * with additional, solver dependent
+ * flags for tuning.
+ */
+ SolverCG (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ protected:
+ /**
+ * Store a copy of the flags for this
+ * particular solver.
+ */
+ const AdditionalData additional_data;
+ };
+
+
+
+/**
+ * An implementation of the solver interface using the Trilinos CGS
+ * solver.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, 2008
+ */
+ class SolverCGS : public SolverBase
+ {
+ public:
+ /**
+ * Standardized data struct to
+ * pipe additional data to the
+ * solver.
+ */
+ struct AdditionalData
+ {};
+
+ /**
+ * Constructor. In contrast to
+ * deal.II's own solvers, there is no
+ * need to give a vector memory
+ * object.
+ *
+ * The last argument takes a structure
+ * with additional, solver dependent
+ * flags for tuning.
+ */
+ SolverCGS (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ protected:
+ /**
+ * Store a copy of the flags for this
+ * particular solver.
+ */
+ const AdditionalData additional_data;
+ };
+
+
+
+/**
+ * An implementation of the solver interface using the Trilinos GMRES
+ * solver.
+ *
+ * @author Martin Kronbichler, 2008
+ */
+ class SolverGMRES : public SolverBase
+ {
+ public:
+ /**
+ * Standardized data struct to
+ * pipe additional data to the
+ * solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, set the
+ * number of temporary vectors to
+ * 30, i.e. do a restart every 30
+ * iterations.
+ */
+ AdditionalData (const unsigned int restart_parameter = 30);
+
+ /**
+ * Maximum number of
+ * tmp vectors.
+ */
+ unsigned int restart_parameter;
+ };
+
+ /**
+ * Constructor. In contrast to
+ * deal.II's own solvers, there is no
+ * need to give a vector memory
+ * object.
+ *
+ * The last argument takes a structure
+ * with additional, solver dependent
+ * flags for tuning.
+ */
+ SolverGMRES (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ protected:
+ /**
+ * Store a copy of the flags for this
+ * particular solver.
+ */
+ const AdditionalData additional_data;
+ };
+
+
+
+/**
+ * An implementation of the solver interface using the Trilinos BiCGStab
+ * solver.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, 2008
+ */
+ class SolverBicgstab : public SolverBase
+ {
+ public:
+ /**
+ * Standardized data struct to
+ * pipe additional data to the
+ * solver.
+ */
+ struct AdditionalData
+ {};
+
+ /**
+ * Constructor. In contrast to
+ * deal.II's own solvers, there is no
+ * need to give a vector memory
+ * object.
+ *
+ * The last argument takes a structure
+ * with additional, solver dependent
+ * flags for tuning.
+ */
+ SolverBicgstab (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ protected:
+ /**
+ * Store a copy of the flags for this
+ * particular solver.
+ */
+ const AdditionalData additional_data;
+ };
+
+
+
+/**
+ * An implementation of the solver interface using the Trilinos TFQMR
+ * solver.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, 2008
+ */
+ class SolverTFQMR : public SolverBase
+ {
+ public:
+ /**
+ * Standardized data struct to
+ * pipe additional data to the
+ * solver.
+ */
+ struct AdditionalData
+ {};
+
+ /**
+ * Constructor. In contrast to
+ * deal.II's own solvers, there is no
+ * need to give a vector memory
+ * object.
+ *
+ * The last argument takes a structure
+ * with additional, solver dependent
+ * flags for tuning.
+ */
+ SolverTFQMR (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ protected:
+ /**
+ * Store a copy of the flags for this
+ * particular solver.
+ */
+ const AdditionalData additional_data;
+ };
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS
+
+/*---------------------------- trilinos_solver.h ---------------------------*/
+
+#endif
+/*---------------------------- trilinos_solver.h ---------------------------*/
// $Id$
// Version: $Name$
//
-// Copyright (C) 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+// Copyright (C) 2008 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
#endif // DEAL_II_USE_TRILINOS
-/*---------------------------- trilinos_matrix_base.h ---------------------------*/
+/*---------------------------- trilinos_sparse_matrix.h ---------------------------*/
#endif
-/*---------------------------- trilinos_matrix_base.h ---------------------------*/
+/*---------------------------- trilinos_sparse_matrix.h ---------------------------*/
//
// Copyright (C) 2008 by the deal.II authors
//
-// This file is subject to QPL and may not be distributed
+// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
// to the file deal.II/doc/license.html for the text and
// further information on this license.
#ifdef DEAL_II_USE_TRILINOS
#include <Ifpack.h>
+#include <Teuchos_ParameterList.hpp>
DEAL_II_NAMESPACE_OPEN
namespace TrilinosWrappers
{
+ PreconditionBase::PreconditionBase()
+ {}
+
+
+
+ void
+ PreconditionBase::vmult (VectorBase &dst,
+ const VectorBase &src) const
+ {
+ Assert (dst.vector->Map().SameAs(preconditioner->OperatorRangeMap()),
+ ExcNonMatchingMaps("dst"));
+ Assert (src.vector->Map().SameAs(preconditioner->OperatorDomainMap()),
+ ExcNonMatchingMaps("src"));
+
+ const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+/* -------------------------- PreconditionJacobi -------------------------- */
+
PreconditionJacobi::AdditionalData::
AdditionalData (const double omega,
const double min_diagonal)
min_diagonal (min_diagonal)
{}
-
-
- PreconditionJacobi::PreconditionJacobi ()
- {}
-
void
- void
- PreconditionJacobi::vmult (VectorBase &dst,
- const VectorBase &src) const
- {
- const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
+/* -------------------------- PreconditionSSOR -------------------------- */
PreconditionSSOR::AdditionalData::
AdditionalData (const double omega,
overlap (overlap)
{}
-
-
- PreconditionSSOR::PreconditionSSOR ()
- {}
-
void
- void
- PreconditionSSOR::vmult (VectorBase &dst,
- const VectorBase &src) const
- {
- const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
+/* -------------------------- PreconditionSOR -------------------------- */
PreconditionSOR::AdditionalData::
AdditionalData (const double omega,
overlap (overlap)
{}
-
-
- PreconditionSOR::PreconditionSOR ()
- {}
-
void
- void
- PreconditionSOR::vmult (VectorBase &dst,
- const VectorBase &src) const
- {
- const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
+/* -------------------------- PreconditionIC -------------------------- */
PreconditionIC::AdditionalData::
AdditionalData (const unsigned int ic_fill,
- PreconditionIC::PreconditionIC ()
- {}
-
-
-
void
PreconditionIC::initialize (const SparseMatrix &matrix,
const AdditionalData &additional_data)
- void
- PreconditionIC::vmult (VectorBase &dst,
- const VectorBase &src) const
- {
- const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
+/* -------------------------- PreconditionILU -------------------------- */
PreconditionILU::AdditionalData::
AdditionalData (const unsigned int ilu_fill,
- PreconditionILU::PreconditionILU ()
- {}
-
-
-
void
PreconditionILU::initialize (const SparseMatrix &matrix,
const AdditionalData &additional_data)
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
-
-
- void
- PreconditionILU::vmult (VectorBase &dst,
- const VectorBase &src) const
- {
- const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
}
DEAL_II_NAMESPACE_CLOSE
#include <Epetra_Map.h>
#include <Epetra_MultiVector.h>
#include <Epetra_Vector.h>
-#include <Teuchos_ParameterList.hpp>
#include <ml_include.h>
#include <ml_MultiLevelPreconditioner.h>
PreconditionAMG::PreconditionAMG ()
#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
- :
- communicator (MPI_COMM_WORLD)
+ :
+ communicator (MPI_COMM_WORLD)
#endif
{}
-
- PreconditionAMG::~PreconditionAMG ()
- {}
-
-
void
PreconditionAMG::
parameter_list.set("null space: vectors", null_space_modes.Values());
}
- multigrid_operator = boost::shared_ptr<ML_Epetra::MultiLevelPreconditioner>
- (new ML_Epetra::MultiLevelPreconditioner(
+ multilevel_operator = Teuchos::rcp (new ML_Epetra::MultiLevelPreconditioner(
*matrix.matrix, parameter_list, true));
if (output_details)
- multigrid_operator->PrintUnused(0);
+ multilevel_operator->PrintUnused(0);
}
void PreconditionAMG::
reinit ()
{
- multigrid_operator->ReComputePreconditioner();
+ multilevel_operator->ReComputePreconditioner();
}
-
-
-
+
+
+
void PreconditionAMG::vmult (VectorBase &dst,
const VectorBase &src) const
{
- const int ierr = multigrid_operator->ApplyInverse (*src.vector,
- *dst.vector);
-
- Assert (ierr == 0, ExcTrilinosError(ierr));
+ Assert (dst.vector->Map().SameAs(multilevel_operator->OperatorRangeMap()),
+ ExcNonMatchingMaps("dst"));
+ Assert (src.vector->Map().SameAs(multilevel_operator->OperatorDomainMap()),
+ ExcNonMatchingMaps("src"));
+
+ const int ierr = multilevel_operator->ApplyInverse (*src.vector, *dst.vector);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
}
+
// For the implementation of
// the <code>vmult</code>
// function we note that
const dealii::Vector<double> &src) const
{
Assert (Map->SameAs(Matrix->matrix->RowMap()),
- ExcMessage("The sparse matrix given to the preconditioner uses "
- "a map that is not compatible. Check ML preconditioner "
- "setup."));
+ ExcNonMatchingMaps("dst"));
Assert (Map->SameAs(Matrix->matrix->RowMap()),
- ExcMessage("The sparse matrix given to the preconditioner uses "
- "a map that is not compatible. Check ML preconditioner "
- "setup."));
+ ExcNonMatchingMaps("src"));
- Epetra_Vector LHS (View, multigrid_operator->OperatorDomainMap(),
+ Epetra_Vector LHS (View, multilevel_operator->OperatorDomainMap(),
dst.begin());
- Epetra_Vector RHS (View, multigrid_operator->OperatorRangeMap(),
+ Epetra_Vector RHS (View, multilevel_operator->OperatorRangeMap(),
const_cast<double*>(src.begin()));
- const int res = multigrid_operator->ApplyInverse (RHS, LHS);
+ const int res = multilevel_operator->ApplyInverse (RHS, LHS);
Assert (res == 0,
ExcMessage ("Trilinos AMG MultiLevel preconditioner returned "
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2008 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+
+
+#include <lac/trilinos_solver.h>
+#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>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+
+ SolverBase::AdditionalData::AdditionalData (const unsigned int gmres_restart_parameter)
+ :
+ gmres_restart_parameter (gmres_restart_parameter)
+ {}
+
+
+
+ SolverBase::SolverBase (SolverControl &cn)
+ :
+ solver_control (cn),
+ solver_name (gmres)
+ {}
+
+
+
+ SolverBase::~SolverBase ()
+ {}
+
+
+
+ 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
+ {
+ return solver_control;
+ }
+
+
+
+ void
+ SolverBase::solve (const SparseMatrix &A,
+ VectorBase &x,
+ const VectorBase &b,
+ const Epetra_Operator* preconditioner)
+ {
+ int ierr;
+
+ linear_problem.reset();
+
+ // We need an
+ // Epetra_LinearProblem object
+ // to let the AztecOO solver
+ // know about the matrix and
+ // vectors.
+ linear_problem = std::auto_ptr<Epetra_LinearProblem> (
+ new Epetra_LinearProblem(&*(A.matrix), &*x.vector,
+ &*b.vector));
+
+ // Next we can allocate the
+ // AztecOO solver...
+ solver.SetProblem(*linear_problem);
+
+ // ... and we can specify the
+ // solver to be used.
+ switch (solver_name)
+ {
+ case cg:
+ solver.SetAztecOption(AZ_solver, AZ_cg);
+ break;
+ case cgs:
+ solver.SetAztecOption(AZ_solver, AZ_cgs);
+ break;
+ case gmres:
+ solver.SetAztecOption(AZ_solver, AZ_gmres);
+ solver.SetAztecOption(AZ_kspace, additional_data.gmres_restart_parameter);
+ break;
+ case bicgstab:
+ solver.SetAztecOption(AZ_solver, AZ_bicgstab);
+ break;
+ case tfqmr:
+ solver.SetAztecOption(AZ_solver, AZ_tfqmr);
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+
+ // Introduce the
+ // preconditioner, ...
+ ierr = solver.SetPrecOperator (const_cast<Epetra_Operator*>(preconditioner));
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ // ... set some options, ...
+ solver.SetAztecOption (AZ_output, AZ_none);
+ solver.SetAztecOption (AZ_conv, AZ_noscaled);
+
+ // ... and then solve!
+ ierr = solver.Iterate (solver_control.max_steps(),
+ solver_control.tolerance());
+ AssertThrow (ierr >= 0, ExcTrilinosError(ierr));
+
+ // Finally, let the deal.II
+ // SolverControl object know
+ // what has happened. If the
+ // solve succeeded, the status
+ // of the solver control will
+ // turn into
+ // SolverControl::success.
+ solver_control.check (solver.NumIters(), solver.TrueResidual());
+
+ if (solver_control.last_check() != SolverControl::success)
+ throw SolverControl::NoConvergence (solver_control.last_step(),
+ solver_control.last_value());
+ }
+
+
+
+
+
+/* ---------------------- SolverCG ------------------------ */
+
+ SolverCG::SolverCG (SolverControl &cn,
+ const AdditionalData &data)
+ :
+ SolverBase (cn),
+ additional_data (data)
+ {
+ solver_name = cg;
+ }
+
+
+/* ---------------------- SolverGMRES ------------------------ */
+
+ SolverGMRES::AdditionalData::
+ AdditionalData (const unsigned int restart_parameter)
+ :
+ restart_parameter (restart_parameter)
+ {}
+
+
+
+ SolverGMRES::SolverGMRES (SolverControl &cn,
+ const AdditionalData &data)
+ :
+ SolverBase (cn),
+ additional_data (data.restart_parameter)
+ {
+ solver_name = gmres;
+ }
+
+
+/* ---------------------- SolverBicgstab ------------------------ */
+
+ SolverBicgstab::SolverBicgstab (SolverControl &cn,
+ const AdditionalData &data)
+ :
+ SolverBase (cn),
+ additional_data (data)
+ {
+ solver_name = bicgstab;
+ }
+
+
+/* ---------------------- SolverCGS ------------------------ */
+
+ SolverCGS::SolverCGS (SolverControl &cn,
+ const AdditionalData &data)
+ :
+ SolverBase (cn),
+ additional_data (data)
+ {
+ solver_name = cgs;
+ }
+
+
+/* ---------------------- SolverTFQMR ------------------------ */
+
+ SolverTFQMR::SolverTFQMR (SolverControl &cn,
+ const AdditionalData &data)
+ :
+ SolverBase (cn),
+ additional_data (data)
+ {
+ solver_name = tfqmr;
+ }
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_PETSC
// $Id$
// Version: $Name$
//
-// Copyright (C) 2004, 2005, 2006, 2008 by the deal.II authors
+// Copyright (C) 2008 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer