From: kronbichler Date: Wed, 8 Oct 2008 13:09:28 +0000 (+0000) Subject: Made the interface to Trilinos preconditioners more generic. Now AMG can be controlle... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fa83ab5cd4ba18c251dfcfed0a921c80500cb0a7;p=dealii-svn.git Made the interface to Trilinos preconditioners more generic. Now AMG can be controlled in the same way as Ifpack preconditioners. Moreover, AMG now has an AdditionalData() field which takes precondition parameters. git-svn-id: https://svn.dealii.org/trunk@17145 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index b17f68bf53..6ef4629c47 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -54,15 +54,11 @@ // 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 #include #include #include -#include // Finally, here are two C++ headers that // haven't been included yet by one of the @@ -1231,7 +1227,8 @@ BoussinesqFlowProblem::build_stokes_preconditioner () 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 diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 853aecf245..d1b125f44e 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -28,7 +28,6 @@ #include #include #include -#include #include #include @@ -805,7 +804,8 @@ BoussinesqFlowProblem::build_stokes_preconditioner () 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 (new TrilinosWrappers::PreconditionIC()); diff --git a/deal.II/lac/include/lac/trilinos_precondition.h b/deal.II/lac/include/lac/trilinos_precondition.h index 2093233864..54db649b4d 100755 --- a/deal.II/lac/include/lac/trilinos_precondition.h +++ b/deal.II/lac/include/lac/trilinos_precondition.h @@ -16,14 +16,31 @@ #include #include +#include +#include + +#include #ifdef DEAL_II_USE_TRILINOS +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI +# include +#else +# include +#endif +#include + #include +#include // forward declarations class Ifpack_Preconditioner; +namespace ML_Epetra +{ + class MultiLevelPreconditioner; +} + DEAL_II_NAMESPACE_OPEN @@ -68,6 +85,17 @@ namespace TrilinosWrappers * sparse matrix. */ PreconditionBase (); + + /** + * 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. @@ -75,6 +103,16 @@ namespace TrilinosWrappers 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 &dst, + const dealii::Vector &src) const; + /** * Exception. */ @@ -88,13 +126,33 @@ namespace TrilinosWrappers 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 preconditioner; + + /** + * Internal communication + * pattern in case the matrix + * needs to be copied from + * deal.II format. */ - Teuchos::RefCountPtr 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 map; }; @@ -156,10 +214,25 @@ namespace TrilinosWrappers 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; }; /** @@ -171,6 +244,14 @@ namespace TrilinosWrappers */ 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; }; @@ -257,11 +338,34 @@ namespace TrilinosWrappers 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; }; @@ -276,6 +380,14 @@ namespace TrilinosWrappers */ 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; }; @@ -360,13 +472,36 @@ namespace TrilinosWrappers 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; }; @@ -380,7 +515,15 @@ namespace TrilinosWrappers * 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; }; @@ -490,13 +633,55 @@ namespace TrilinosWrappers 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 + * ic_fill 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; }; /** @@ -508,6 +693,14 @@ namespace TrilinosWrappers */ 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; }; @@ -616,13 +809,55 @@ namespace TrilinosWrappers 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 + * ilu_fill 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; }; /** @@ -634,6 +869,14 @@ namespace TrilinosWrappers */ 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; }; @@ -673,8 +916,13 @@ namespace TrilinosWrappers */ 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; }; @@ -688,9 +936,219 @@ namespace TrilinosWrappers */ 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; }; + +/** + * 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 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. + * + * 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 > &null_space = std::vector > (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 + * strong coupling is + * controlled by the variable + * aggregation_threshold, + * meaning that all elements + * that are not smaller than + * aggregation_threshold + * 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 > 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 + * true, 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 &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 multilevel_operator; + + /** + * A copy of the deal.II matrix + * into Trilinos format. + */ + boost::shared_ptr Matrix; + }; + } /*@}*/ diff --git a/deal.II/lac/source/trilinos_precondition.cc b/deal.II/lac/source/trilinos_precondition.cc index 9b4c4ae03c..b09cae7295 100755 --- a/deal.II/lac/source/trilinos_precondition.cc +++ b/deal.II/lac/source/trilinos_precondition.cc @@ -20,6 +20,11 @@ #include #include +#include +#include +#include +#include +#include DEAL_II_NAMESPACE_OPEN @@ -28,10 +33,21 @@ namespace TrilinosWrappers { 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 @@ -47,6 +63,45 @@ namespace TrilinosWrappers + // For the implementation of + // the vmult + // 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 + // const in all + // deal.II calls) to + // non-constant value, as this + // is the way Trilinos wants + // to have them. + void PreconditionBase::vmult (dealii::Vector &dst, + const dealii::Vector &src) const + { + + Epetra_Vector LHS (View, preconditioner->OperatorDomainMap(), + dst.begin()); + Epetra_Vector RHS (View, preconditioner->OperatorRangeMap(), + const_cast(src.begin())); + + const int res = preconditioner->ApplyInverse (RHS, LHS); + + Assert (res == 0, + ExcMessage ("Trilinos AMG MultiLevel preconditioner returned " + "with an error!")); + } + + + /* -------------------------- PreconditionJacobi -------------------------- */ PreconditionJacobi::AdditionalData:: @@ -64,10 +119,13 @@ namespace TrilinosWrappers 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; @@ -78,14 +136,16 @@ namespace TrilinosWrappers 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); } @@ -109,12 +169,14 @@ namespace TrilinosWrappers 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; @@ -126,14 +188,16 @@ namespace TrilinosWrappers 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); } @@ -157,12 +221,14 @@ namespace TrilinosWrappers 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; @@ -174,14 +240,16 @@ namespace TrilinosWrappers 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); } @@ -207,10 +275,12 @@ namespace TrilinosWrappers 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; @@ -221,14 +291,16 @@ namespace TrilinosWrappers 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); } @@ -254,10 +326,12 @@ namespace TrilinosWrappers 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; @@ -268,14 +342,16 @@ namespace TrilinosWrappers 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); } @@ -295,25 +371,164 @@ namespace TrilinosWrappers 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 > &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; dGID(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 &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 (new SparseMatrix()); + + Matrix->reinit (*map, deal_ii_sparse_matrix, drop_tolerance); + Matrix->compress(); + + initialize (*Matrix, additional_data); + } + + + void PreconditionAMG:: + reinit () + { + multilevel_operator->ReComputePreconditioner(); } } diff --git a/deal.II/lac/source/trilinos_solver.cc b/deal.II/lac/source/trilinos_solver.cc index 7f865051f3..57ba88ac89 100644 --- a/deal.II/lac/source/trilinos_solver.cc +++ b/deal.II/lac/source/trilinos_solver.cc @@ -16,14 +16,12 @@ #include #include #include -#include #include #ifdef DEAL_II_USE_TRILINOS -#include -#include +#include DEAL_II_NAMESPACE_OPEN @@ -50,28 +48,6 @@ namespace TrilinosWrappers - 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 { @@ -84,7 +60,7 @@ namespace TrilinosWrappers SolverBase::solve (const SparseMatrix &A, VectorBase &x, const VectorBase &b, - const Epetra_Operator* preconditioner) + const PreconditionBase &preconditioner) { int ierr; @@ -129,7 +105,8 @@ namespace TrilinosWrappers // Introduce the // preconditioner, ... - ierr = solver.SetPrecOperator (const_cast(preconditioner)); + ierr = solver.SetPrecOperator (const_cast + (&*preconditioner.preconditioner)); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); // ... set some options, ...