From: kronbichler Date: Fri, 19 Sep 2008 11:45:20 +0000 (+0000) Subject: Implemented a few more wrappers to Trilinos Ifpack preconditioners (Jacobi, SOR,... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b4801d2a76851b7d2736ad124c131b61b6416f1c;p=dealii-svn.git Implemented a few more wrappers to Trilinos Ifpack preconditioners (Jacobi, SOR, IC, ILU). Use the same interface as for all the other interfaces, i.e., let an initialize() function do the actual job. git-svn-id: https://svn.dealii.org/trunk@16855 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 818bf64b99..db3be88046 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -1144,8 +1144,8 @@ BoussinesqFlowProblem::build_stokes_preconditioner () // keep the (1,1) block, though Mp_preconditioner = boost::shared_ptr - (new TrilinosWrappers::PreconditionSSOR( - stokes_preconditioner_matrix.block(1,1),1.2)); + (new TrilinosWrappers::PreconditionSSOR()); + Mp_preconditioner->initialize(stokes_preconditioner_matrix.block(1,1),1.2); std::cout << std::endl; @@ -1841,11 +1841,11 @@ void BoussinesqFlowProblem::solve () 1e-8*temperature_rhs.l2_norm()); SolverCG cg (solver_control); - TrilinosWrappers::PreconditionSSOR preconditioner (temperature_matrix, - 1.2); + TrilinosWrappers::PreconditionSSOR preconditioner; + preconditioner.initialize (temperature_matrix, 1.2); + cg.solve (temperature_matrix, temperature_solution, - temperature_rhs, - preconditioner); + temperature_rhs, preconditioner); // produce a consistent temperature field temperature_constraints.distribute (temperature_solution); diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index b72f043873..ec451be4a2 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -800,8 +800,9 @@ BoussinesqFlowProblem::build_stokes_preconditioner () true, true, 5e-2, null_space, false); Mp_preconditioner = boost::shared_ptr - (new TrilinosWrappers::PreconditionSSOR( - stokes_preconditioner_matrix.block(1,1),1.2)); + (new TrilinosWrappers::PreconditionSSOR()); + + Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1),1.2); pcout << std::endl; @@ -1286,11 +1287,11 @@ void BoussinesqFlowProblem::solve () 1e-8*temperature_rhs.l2_norm()); SolverCG cg (solver_control); - TrilinosWrappers::PreconditionSSOR preconditioner (temperature_matrix, - 1.2); + TrilinosWrappers::PreconditionSSOR preconditioner; + preconditioner.initialize (temperature_matrix, 1.2); + cg.solve (temperature_matrix, temperature_solution, - temperature_rhs, - preconditioner); + temperature_rhs, preconditioner); TrilinosWrappers::Vector localized_temperature_solution (temperature_solution); temperature_constraints.distribute (localized_temperature_solution); diff --git a/deal.II/lac/include/lac/trilinos_precondition.h b/deal.II/lac/include/lac/trilinos_precondition.h index aebf1abc02..280d2f4362 100755 --- a/deal.II/lac/include/lac/trilinos_precondition.h +++ b/deal.II/lac/include/lac/trilinos_precondition.h @@ -15,11 +15,13 @@ #include -#include +#include #ifdef DEAL_II_USE_TRILINOS +#include + // forward declarations class Ifpack_Preconditioner; @@ -36,43 +38,487 @@ namespace TrilinosWrappers class VectorBase; class SparseMatrix; +/** + * @ingroup TrilinosWrappers + * @ingroup Preconditioners + * @author Martin Kronbichler, 2008 + */ + class PreconditionJacobi : public Subscriptor + { + public: + + /** + * Standardized data struct to + * pipe additional flags to the + * preconditioner. The + * parameter omega + * specifies the relaxation + * parameter in the Jacobi + * preconditioner. The + * parameter + * min_diagonal can be + * used to make the application + * of the preconditioner also + * possible when some diagonal + * elements are zero. In a + * default application this + * would mean that we divide by + * zero, so by setting the + * parameter + * min_diagonal to a + * small nonzero value the SOR + * will work on a matrix that + * is not too far away from the + * one we want to + * treat. + */ + struct AdditionalData + { + /** + * Constructor. By default, set + * the damping parameter to + * one, and do not modify the + * diagonal. + */ + AdditionalData (const double omega = 1, + const double min_diagonal = 0); + + /** + * Relaxation parameter and + * minimal diagonal value. + */ + double omega, min_diagonal; + }; + + /** + * 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. + */ + 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; + }; + + + + /** * @ingroup TrilinosWrappers * @ingroup Preconditioners * @author Wolfgang Bangerth, 2008 */ - class PreconditionSSOR + class PreconditionSSOR : public Subscriptor { public: + + /** + * Standardized data struct to + * pipe additional flags to the + * preconditioner. The + * parameter omega + * specifies the relaxation + * parameter in the SSOR + * preconditioner. The + * parameter + * min_diagonal can be + * used to make the application + * of the preconditioner also + * possible when some diagonal + * elements are zero. In a + * default application this + * would mean that we divide by + * zero, so by setting the + * parameter + * min_diagonal to a + * small nonzero value the SOR + * will work on a matrix that + * is not too far away from the + * one we want to + * treat. Finally, + * overlap governs the + * overlap of the partitions + * when the preconditioner runs + * in parallel, forming a + * so-called additive Schwarz + * preconditioner. + */ + struct AdditionalData + { + /** + * Constructor. By default, set + * the damping parameter to + * one, we do not modify the + * diagonal, and there is no + * overlap (i.e. in parallel, + * we run a BlockJacobi + * preconditioner, where each + * block is inverted + * approximately by an SSOR. + */ + AdditionalData (const double omega = 1, + const double min_diagonal = 0, + const unsigned int overlap = 0); + + /** + * Relaxation parameter, + * minimal diagonal element, + * and overlap. + */ + double omega, min_diagonal; + unsigned int overlap; + }; + + /** + * 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 + * 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; + }; + + + + +/** + * @ingroup TrilinosWrappers + * @ingroup Preconditioners + * @author Martin Kronbichler, 2008 + */ + class PreconditionSOR : public Subscriptor + { public: + /** * Standardized data struct to * pipe additional flags to the + * preconditioner. The + * parameter omega + * specifies the relaxation + * parameter in the SOR + * preconditioner. The + * parameter + * min_diagonal can be + * used to make the application + * of the preconditioner also + * possible when some diagonal + * elements are zero. In a + * default application this + * would mean that we divide by + * zero, so by setting the + * parameter + * min_diagonal to a + * small nonzero value the SOR + * will work on a matrix that + * is not too far away from the + * one we want to + * treat. Finally, + * overlap governs the + * overlap of the partitions + * when the preconditioner runs + * in parallel, forming a + * so-called additive Schwarz * preconditioner. */ struct AdditionalData { - /** - * Constructor. By default, - * set the damping parameter - * to one. - */ - AdditionalData (const double omega = 1); + /** + * Constructor. By default, set + * the damping parameter to + * one, we do not modify the + * diagonal, and there is no + * overlap (i.e. in parallel, + * we run a BlockJacobi + * preconditioner, where each + * block is inverted + * approximately by an SOR. + */ + AdditionalData (const double omega = 1, + const double min_diagonal = 0, + const unsigned int overlap = 0); - /** - * Relaxation parameter. - */ - double omega; + /** + * Relaxation parameter, + * minimal diagonal element, + * and overlap. + */ + double omega, min_diagonal; + unsigned int overlap; + }; + + /** + * 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 + * 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; + }; + + + +/** + * @ingroup TrilinosWrappers + * @ingroup Preconditioners + * @author Martin Kronbichler, 2008 + */ + class PreconditionIC : public Subscriptor + { + public: + /** + * Standardized data struct to + * pipe additional parameters + * to the preconditioner. The + * Trilinos IC decomposition + * allows for some fill-in, so + * it actually is a threshold + * incomplete Cholesky + * factorization. The amount of + * fill-in, and hence, the + * amount of memory used by + * this preconditioner, is + * controlled by the parameter + * ic_fill, which + * specifies this as a + * double. When forming the + * preconditioner, for certain + * problems bad conditioning + * (or just bad luck) can cause + * the preconditioner to be + * very poorly + * conditioned. Hence it can + * help to add diagonal + * perturbations to the + * original matrix and form the + * preconditioner for this + * slightly better + * matrix. ic_atol is + * an absolute perturbation + * that is added to the + * diagonal before forming the + * prec, and ic_rtol + * is a scaling factor $rtol + * \geq 1$. The last parameter + * specifies the overlap of the + * partitions when the + * preconditioner runs in + * parallel. + */ + struct AdditionalData + { + /** + * Constructor. By default, set + * the drop tolerance to 0, the + * level of extra fill-ins is + * set to be zero (just use the + * matrix structure, do not + * generate any additional + * fill-in), the tolerance + * level are 0 and 1, + * respectively, and the + * overlap in case of a + * parallel execution is + * zero. This overlap in a + * block-application of the IC + * in the parallel case makes + * the preconditioner a + * so-called additive Schwarz + * preconditioner. + */ + AdditionalData (const double ic_fill = 0., + const double ic_atol = 0., + const double ic_rtol = 1., + const unsigned int overlap = 0); + + /** + * IC parameters and overlap. + */ + double ic_fill, ic_atol, ic_rtol; + unsigned int overlap; }; /** - * Constructor. Take the matrix which - * is used to form the preconditioner, - * and additional flags if there are - * any. + * (Empty) constructor. + */ + PreconditionIC (); + + /** + * Initialize function. Takes + * the matrix which is used to + * form the preconditioner, 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; + }; + + + +/** + * @ingroup TrilinosWrappers + * @ingroup Preconditioners + * @author Martin Kronbichler, 2008 + */ + class PreconditionILU : public Subscriptor + { + public: + /** + * Standardized data struct to + * pipe additional parameters + * to the preconditioner. The + * Trilinos ILU decomposition + * allows for some fill-in, so + * it actually is a threshold + * incomplete LU + * factorization. The amount of + * fill-in, and hence, the + * amount of memory used by + * this preconditioner, is + * controlled by the parameter + * ilu_fill, which + * specifies this as a + * double. When forming the + * preconditioner, for certain + * problems bad conditioning + * (or just bad luck) can cause + * the preconditioner to be + * very poorly + * conditioned. Hence it can + * help to add diagonal + * perturbations to the + * original matrix and form the + * preconditioner for this + * slightly better + * matrix. ilu_atol is + * an absolute perturbation + * that is added to the + * diagonal before forming the + * prec, and ilu_rtol + * is a scaling factor $rtol + * \geq 1$. The last parameter + * specifies the overlap of the + * partitions when the + * preconditioner runs in + * parallel. + */ + struct AdditionalData + { + /** + * Constructor. By default, the + * level of extra fill-ins is + * set to be zero (just use the + * matrix structure, do not + * generate any additional + * fill-in), the tolerance + * level are 0 and 1, + * respectively, and the + * overlap in case of a + * parallel execution is + * zero. This overlap in a + * block-application of the ILU + * in the parallel case makes + * the preconditioner a + * so-called additive Schwarz + * preconditioner. + */ + AdditionalData (const double ilu_fill = 0., + const double ilu_atol = 0., + const double ilu_rtol = 1., + const unsigned int overlap = 0); + + /** + * ILU parameters and overlap. + */ + double ilu_drop, ilu_fill, ilu_atol, ilu_rtol; + unsigned int overlap; + }; + + /** + * (Empty) constructor. + */ + PreconditionILU (); + + /** + * Initialize function. Takes + * the matrix which is used to + * form the preconditioner, and + * additional flags if there + * are any. */ - PreconditionSSOR (const SparseMatrix &matrix, - const AdditionalData &additional_data = AdditionalData()); + void initialize (const SparseMatrix &matrix, + const AdditionalData &additional_data = AdditionalData()); /** * Apply the preconditioner. @@ -81,7 +527,7 @@ namespace TrilinosWrappers const VectorBase &src) const; private: - boost::shared_ptr preconditioner; + Teuchos::RefCountPtr preconditioner; }; @@ -94,7 +540,7 @@ DEAL_II_NAMESPACE_CLOSE #endif // DEAL_II_USE_TRILINOS -/*---------------------------- trilinos_precondition_base.h ---------------------------*/ +/*---------------------------- trilinos_precondition.h ---------------------------*/ #endif -/*---------------------------- trilinos_precondition_base.h ---------------------------*/ +/*---------------------------- trilinos_precondition.h ---------------------------*/ diff --git a/deal.II/lac/source/trilinos_precondition.cc b/deal.II/lac/source/trilinos_precondition.cc index 0d9eb8768f..7a96d6bc3a 100755 --- a/deal.II/lac/source/trilinos_precondition.cc +++ b/deal.II/lac/source/trilinos_precondition.cc @@ -26,26 +26,96 @@ DEAL_II_NAMESPACE_OPEN namespace TrilinosWrappers { + PreconditionJacobi::AdditionalData:: + AdditionalData (const double omega, + const double min_diagonal) + : + omega (omega), + min_diagonal (min_diagonal) + {} + + + + PreconditionJacobi::PreconditionJacobi () + {} + + + + void + PreconditionJacobi::initialize (const SparseMatrix &matrix, + const AdditionalData &additional_data) + { + + preconditioner = Teuchos::rcp (Ifpack().Create ("point relaxation", &*matrix.matrix, 0)); + Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this " + "preconditioner")); + + int ierr; + + Teuchos::ParameterList parameter_list; + parameter_list.set ("relaxation: type", "Jacobi"); + parameter_list.set ("relaxation: damping factor", additional_data.omega); + parameter_list.set ("relaxation: min diagonal value", + additional_data.min_diagonal); + + ierr = preconditioner->SetParameters(parameter_list); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = preconditioner->Initialize(); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = preconditioner->Compute(); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + PreconditionJacobi::vmult (VectorBase &dst, + const VectorBase &src) const + { + const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + PreconditionSSOR::AdditionalData:: - AdditionalData (const double omega) + AdditionalData (const double omega, + const double min_diagonal, + const unsigned int overlap) : - omega (omega) + omega (omega), + min_diagonal (min_diagonal), + overlap (overlap) {} - PreconditionSSOR::PreconditionSSOR (const SparseMatrix &matrix, - const AdditionalData &additional_data) - : - preconditioner (Ifpack().Create ("point relaxation", - &*matrix.matrix, 0)) + + PreconditionSSOR::PreconditionSSOR () + {} + + + + void + PreconditionSSOR::initialize (const SparseMatrix &matrix, + const AdditionalData &additional_data) { - Assert (preconditioner != 0, ExcMessage ("Trilinos could not create this preconditioner")); + + preconditioner = Teuchos::rcp (Ifpack().Create ("point relaxation", + &*matrix.matrix, + additional_data.overlap)); + Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this " + "preconditioner")); int ierr; Teuchos::ParameterList parameter_list; parameter_list.set ("relaxation: type", "symmetric Gauss-Seidel"); parameter_list.set ("relaxation: damping factor", additional_data.omega); + parameter_list.set ("relaxation: min diagonal value", + additional_data.min_diagonal); + parameter_list.set ("schwarz: combine mode", "Add"); ierr = preconditioner->SetParameters(parameter_list); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); @@ -58,11 +128,188 @@ namespace TrilinosWrappers } + void PreconditionSSOR::vmult (VectorBase &dst, const VectorBase &src) const { - preconditioner->ApplyInverse (*src.vector, *dst.vector); + const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + PreconditionSOR::AdditionalData:: + AdditionalData (const double omega, + const double min_diagonal, + const unsigned int overlap) + : + omega (omega), + min_diagonal (min_diagonal), + overlap (overlap) + {} + + + + PreconditionSOR::PreconditionSOR () + {} + + + + void + PreconditionSOR::initialize (const SparseMatrix &matrix, + const AdditionalData &additional_data) + { + + preconditioner = Teuchos::rcp (Ifpack().Create ("point relaxation", + &*matrix.matrix, + additional_data.overlap)); + Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this " + "preconditioner")); + + int ierr; + + Teuchos::ParameterList parameter_list; + parameter_list.set ("relaxation: type", "Gauss-Seidel"); + parameter_list.set ("relaxation: damping factor", additional_data.omega); + parameter_list.set ("relaxation: min diagonal value", + additional_data.min_diagonal); + parameter_list.set ("schwarz: combine mode", "Add"); + + ierr = preconditioner->SetParameters(parameter_list); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = preconditioner->Initialize(); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = preconditioner->Compute(); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + PreconditionSOR::vmult (VectorBase &dst, + const VectorBase &src) const + { + const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + PreconditionIC::AdditionalData:: + AdditionalData (const double ic_fill, + const double ic_atol, + const double ic_rtol, + const unsigned int overlap) + : + ic_fill (ic_fill), + ic_atol (ic_atol), + ic_rtol (ic_rtol), + overlap (overlap) + {} + + + + PreconditionIC::PreconditionIC () + {} + + + + void + PreconditionIC::initialize (const SparseMatrix &matrix, + const AdditionalData &additional_data) + { + preconditioner = Teuchos::rcp (Ifpack().Create ("IC", &*matrix.matrix, + additional_data.overlap)); + Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this " + "preconditioner")); + + int ierr; + + Teuchos::ParameterList parameter_list; + parameter_list.set ("fact: level-of-fill",additional_data.ic_fill); + parameter_list.set ("fact: absolute threshold",additional_data.ic_atol); + parameter_list.set ("fact: relative threshold",additional_data.ic_rtol); + parameter_list.set ("schwarz: combine mode", "Add"); + + ierr = preconditioner->SetParameters(parameter_list); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = preconditioner->Initialize(); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = preconditioner->Compute(); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + PreconditionIC::vmult (VectorBase &dst, + const VectorBase &src) const + { + const int ierr = preconditioner->ApplyInverse (*src.vector, *dst.vector); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + PreconditionILU::AdditionalData:: + AdditionalData (const double ilu_fill, + const double ilu_atol, + const double ilu_rtol, + const unsigned int overlap) + : + ilu_fill (ilu_fill), + ilu_atol (ilu_atol), + ilu_rtol (ilu_rtol), + overlap (overlap) + {} + + + + PreconditionILU::PreconditionILU () + {} + + + + void + PreconditionILU::initialize (const SparseMatrix &matrix, + const AdditionalData &additional_data) + { + preconditioner = Teuchos::rcp (Ifpack().Create ("ILU", &*matrix.matrix, + additional_data.overlap)); + Assert (&*preconditioner != 0, ExcMessage ("Trilinos could not create this " + "preconditioner")); + + int ierr; + + Teuchos::ParameterList parameter_list; + parameter_list.set ("fact: level-of-fill",additional_data.ilu_fill); + parameter_list.set ("fact: absolute threshold",additional_data.ilu_atol); + parameter_list.set ("fact: relative threshold",additional_data.ilu_rtol); + parameter_list.set ("schwarz: combine mode", "Add"); + + ierr = preconditioner->SetParameters(parameter_list); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = preconditioner->Initialize(); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = preconditioner->Compute(); + 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)); } }