From 284b4b0206fa02d957c12fc5c18ce2c15cd95495 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Tue, 23 Sep 2008 07:23:14 +0000 Subject: [PATCH] Trilinos Ifpack preconditioner wrapper classes are now derived from a PreconditionBase class. Add a wrapper class to the Trilinos AztecOO solver for basic (non-block) systems. git-svn-id: https://svn.dealii.org/trunk@16897 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/trilinos_precondition.h | 327 ++++++++----- .../include/lac/trilinos_precondition_amg.h | 63 +-- deal.II/lac/include/lac/trilinos_solver.h | 436 ++++++++++++++++++ .../lac/include/lac/trilinos_sparse_matrix.h | 6 +- deal.II/lac/source/trilinos_precondition.cc | 99 ++-- .../lac/source/trilinos_precondition_amg.cc | 49 +- deal.II/lac/source/trilinos_solver.cc | 233 ++++++++++ deal.II/lac/source/trilinos_sparse_matrix.cc | 2 +- 8 files changed, 973 insertions(+), 242 deletions(-) create mode 100644 deal.II/lac/include/lac/trilinos_solver.h create mode 100644 deal.II/lac/source/trilinos_solver.cc diff --git a/deal.II/lac/include/lac/trilinos_precondition.h b/deal.II/lac/include/lac/trilinos_precondition.h index 64414efab7..dafff8725b 100755 --- a/deal.II/lac/include/lac/trilinos_precondition.h +++ b/deal.II/lac/include/lac/trilinos_precondition.h @@ -25,7 +25,6 @@ // forward declarations class Ifpack_Preconditioner; - DEAL_II_NAMESPACE_OPEN /*! @addtogroup TrilinosWrappers @@ -37,13 +36,85 @@ namespace 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 + * initialize 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 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 omega and a min_diagonal + * 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: @@ -91,45 +162,48 @@ namespace TrilinosWrappers }; /** - * Constructor. Does not do - * anything. The - * initialize 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 preconditioner; - //std::auto_ptr preconditioner; - //boost::shared_ptr 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 omega, a + * min_diagonal argument that can be used to make the + * preconditioner work even if the matrix contains some zero elements + * on the diagonal, and a parameter overlap 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 additive + * Schwarz method with an SSOR approximate solve 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: @@ -191,43 +265,50 @@ namespace TrilinosWrappers }; /** - * Constructor. Does not do - * anything. The - * initialize 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 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 omega, a + * min_diagonal argument that can be used to make the + * preconditioner work even if the matrix contains some zero elements + * on the diagonal, and a parameter overlap 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 additive + * Schwarz method with an SOR approximate solve 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: @@ -289,42 +370,58 @@ namespace TrilinosWrappers }; /** - * Constructor. Does not do - * anything. The - * initialize 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 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 ic_fill, so one can gradually + * choose between a factorization on the sparse matrix structure only + * (ic_fill=0) to a full factorization (ic_fill 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 + * overlap 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 additive + * Schwarz method with an IC approximate solve 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: /** @@ -400,39 +497,57 @@ namespace TrilinosWrappers 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 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 + * ilu_fill, so one can gradually choose between a + * factorization on the sparse matrix structure only + * (ilu_fill=0) to a full factorization (ilu_fill 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 + * overlap 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 additive + * Schwarz method with an ILU approximate solve 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: /** @@ -507,11 +622,6 @@ namespace TrilinosWrappers unsigned int overlap; }; - /** - * (Empty) constructor. - */ - PreconditionILU (); - /** * Initialize function. Takes * the matrix which is used to @@ -521,15 +631,6 @@ namespace TrilinosWrappers */ void initialize (const SparseMatrix &matrix, const AdditionalData &additional_data = AdditionalData()); - - /** - * Apply the preconditioner. - */ - void vmult (VectorBase &dst, - const VectorBase &src) const; - - private: - Teuchos::RefCountPtr preconditioner; }; diff --git a/deal.II/lac/include/lac/trilinos_precondition_amg.h b/deal.II/lac/include/lac/trilinos_precondition_amg.h index 6c034baf32..936d80da30 100755 --- a/deal.II/lac/include/lac/trilinos_precondition_amg.h +++ b/deal.II/lac/include/lac/trilinos_precondition_amg.h @@ -30,6 +30,7 @@ #else # include #endif +#include // some forward declarations namespace ML_Epetra @@ -49,15 +50,19 @@ DEAL_II_NAMESPACE_OPEN 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 - * vmult 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 vmult function does call the respective + * operation in the Trilinos package, where it is called * ApplyInverse. Use of this class is explained in the * @ref step_31 "step-31" tutorial program. * @@ -75,8 +80,9 @@ namespace TrilinosWrappers * 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 @@ -90,11 +96,6 @@ namespace TrilinosWrappers */ PreconditionAMG (); - /** - * Destructor. - */ - ~PreconditionAMG (); - /** * Let Trilinos compute a * multilevel hierarchy for the @@ -108,7 +109,7 @@ namespace TrilinosWrappers const bool elliptic = true, const bool higher_order_elements = false, const double aggregation_threshold = 1e-4, - const std::vector > &null_space = std::vector > (), + const std::vector > &null_space = std::vector > (1), const bool output_details = false); /** @@ -156,16 +157,12 @@ namespace TrilinosWrappers * 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 @@ -176,12 +173,24 @@ namespace TrilinosWrappers void vmult (dealii::Vector &dst, const dealii::Vector &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 multigrid_operator; + Teuchos::RefCountPtr + multilevel_operator; /** * Internal communication @@ -200,13 +209,15 @@ namespace TrilinosWrappers * case the matrix needs to be * copied from deal.II format. */ - boost::shared_ptr Map; + boost::shared_ptr Map; /** * A copy of the deal.II matrix * into Trilinos format. */ boost::shared_ptr Matrix; + + friend class SolverBase; }; } diff --git a/deal.II/lac/include/lac/trilinos_solver.h b/deal.II/lac/include/lac/trilinos_solver.h new file mode 100644 index 0000000000..a2294e6fa9 --- /dev/null +++ b/deal.II/lac/include/lac/trilinos_solver.h @@ -0,0 +1,436 @@ +//--------------------------------------------------------------------------- +// $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 +#include +#include +#include + + +#ifdef DEAL_II_USE_TRILINOS + +#include +#include +#include + +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 AztecOO + * user guide. + * + * @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 + * Ax=b. 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 + * Ax=b. In contrast + * to the other vmult + * 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 "cg", "cgs", + * "gmres", "bicgstab", + * and "tfqmr". + */ + 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 vmult + * 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 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 ---------------------------*/ diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index 2363ca35fd..6639ae82c2 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -2,7 +2,7 @@ // $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 @@ -1674,7 +1674,7 @@ DEAL_II_NAMESPACE_CLOSE #endif // DEAL_II_USE_TRILINOS -/*---------------------------- trilinos_matrix_base.h ---------------------------*/ +/*---------------------------- trilinos_sparse_matrix.h ---------------------------*/ #endif -/*---------------------------- trilinos_matrix_base.h ---------------------------*/ +/*---------------------------- trilinos_sparse_matrix.h ---------------------------*/ diff --git a/deal.II/lac/source/trilinos_precondition.cc b/deal.II/lac/source/trilinos_precondition.cc index 33f50b52b3..5d9c7af690 100755 --- a/deal.II/lac/source/trilinos_precondition.cc +++ b/deal.II/lac/source/trilinos_precondition.cc @@ -4,7 +4,7 @@ // // 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. @@ -19,6 +19,7 @@ #ifdef DEAL_II_USE_TRILINOS #include +#include DEAL_II_NAMESPACE_OPEN @@ -26,6 +27,28 @@ 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) @@ -34,11 +57,6 @@ namespace TrilinosWrappers min_diagonal (min_diagonal) {} - - - PreconditionJacobi::PreconditionJacobi () - {} - void @@ -70,15 +88,7 @@ namespace TrilinosWrappers - 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, @@ -90,11 +100,6 @@ namespace TrilinosWrappers overlap (overlap) {} - - - PreconditionSSOR::PreconditionSSOR () - {} - void @@ -129,15 +134,7 @@ namespace TrilinosWrappers - 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, @@ -149,11 +146,6 @@ namespace TrilinosWrappers overlap (overlap) {} - - - PreconditionSOR::PreconditionSOR () - {} - void @@ -188,15 +180,7 @@ namespace TrilinosWrappers - 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, @@ -212,11 +196,6 @@ namespace TrilinosWrappers - PreconditionIC::PreconditionIC () - {} - - - void PreconditionIC::initialize (const SparseMatrix &matrix, const AdditionalData &additional_data) @@ -246,15 +225,7 @@ namespace TrilinosWrappers - 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, @@ -270,11 +241,6 @@ namespace TrilinosWrappers - PreconditionILU::PreconditionILU () - {} - - - void PreconditionILU::initialize (const SparseMatrix &matrix, const AdditionalData &additional_data) @@ -302,15 +268,6 @@ namespace TrilinosWrappers 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 diff --git a/deal.II/lac/source/trilinos_precondition_amg.cc b/deal.II/lac/source/trilinos_precondition_amg.cc index 765849e684..4f702efc57 100755 --- a/deal.II/lac/source/trilinos_precondition_amg.cc +++ b/deal.II/lac/source/trilinos_precondition_amg.cc @@ -19,7 +19,6 @@ #include #include #include -#include #include #include @@ -30,17 +29,12 @@ namespace TrilinosWrappers PreconditionAMG::PreconditionAMG () #ifdef DEAL_II_COMPILER_SUPPORTS_MPI - : - communicator (MPI_COMM_WORLD) + : + communicator (MPI_COMM_WORLD) #endif {} - - PreconditionAMG::~PreconditionAMG () - {} - - void PreconditionAMG:: @@ -114,12 +108,11 @@ namespace TrilinosWrappers parameter_list.set("null space: vectors", null_space_modes.Values()); } - multigrid_operator = boost::shared_ptr - (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); } @@ -155,21 +148,25 @@ namespace TrilinosWrappers 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 vmult // function we note that @@ -195,20 +192,16 @@ namespace TrilinosWrappers const dealii::Vector &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(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 " diff --git a/deal.II/lac/source/trilinos_solver.cc b/deal.II/lac/source/trilinos_solver.cc new file mode 100644 index 0000000000..7f865051f3 --- /dev/null +++ b/deal.II/lac/source/trilinos_solver.cc @@ -0,0 +1,233 @@ +//--------------------------------------------------------------------------- +// $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 +#include +#include +#include +#include + +#include + +#ifdef DEAL_II_USE_TRILINOS + +#include +#include + +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 ( + 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(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 diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 80e3f2351e..a4da44b294 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -2,7 +2,7 @@ // $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 -- 2.39.5