From ddc7bdfaf102c992a794f912d6be85b210033133 Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 21 Aug 2008 20:35:47 +0000 Subject: [PATCH] Add the Trilinos algebraic multigrid preconditioner. git-svn-id: https://svn.dealii.org/trunk@16641 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/lac/trilinos_preconditioner_amg.h | 110 +++++++++++ .../lac/source/trilinos_preconditioner_amg.cc | 185 ++++++++++++++++++ 2 files changed, 295 insertions(+) create mode 100755 deal.II/lac/include/lac/trilinos_preconditioner_amg.h create mode 100755 deal.II/lac/source/trilinos_preconditioner_amg.cc diff --git a/deal.II/lac/include/lac/trilinos_preconditioner_amg.h b/deal.II/lac/include/lac/trilinos_preconditioner_amg.h new file mode 100755 index 0000000000..9af75d80b2 --- /dev/null +++ b/deal.II/lac/include/lac/trilinos_preconditioner_amg.h @@ -0,0 +1,110 @@ +//--------------------------------------------------------------------------- +// $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_precondition_amg_h +#define __deal2__trilinos_precondition_amg_h + + +#include +#include +#include + +#include + + +#ifdef DEAL_II_USE_TRILINOS + +// some forward declarations +namespace ML_Epetra +{ + class MultiLevelPreconditioner; +} +class Epetra_Map; +class Epetra_CrsMatrix; + + + +DEAL_II_NAMESPACE_OPEN + +namespace TrilinosWrappers +{ + // This implements an AMG + // preconditioner based on the + // Trilinos ML implementation. + // What this class does is twofold. + // When the constructor of the class + // is invoked, a ML preconditioner + // object is created based on the + // DoFHandler and 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. + + // There are a few pecularities in + // the constructor. 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. + class PreconditionerTrilinosAmg : public Subscriptor + { + public: + PreconditionerTrilinosAmg (); + ~PreconditionerTrilinosAmg (); + + void initialize (const dealii::SparseMatrix &matrix, + const std::vector &null_space, + const unsigned int null_space_dimension, + const bool higher_order_elements, + const bool elliptic, + const bool output_details, + const double drop_tolerance = 1e-13); + + void vmult (dealii::Vector &dst, + const dealii::Vector &src) const; + + private: + + boost::shared_ptr multigrid_operator; + + Epetra_SerialComm communicator; + boost::shared_ptr Map; + boost::shared_ptr Matrix; + }; +} + + + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_USE_TRILINOS + +/*---------------------------- trilinos_precondition_amg_base.h ---------------------------*/ + +#endif +/*---------------------------- trilinos_precondition_amg_base.h ---------------------------*/ diff --git a/deal.II/lac/source/trilinos_preconditioner_amg.cc b/deal.II/lac/source/trilinos_preconditioner_amg.cc new file mode 100755 index 0000000000..72cb844b5a --- /dev/null +++ b/deal.II/lac/source/trilinos_preconditioner_amg.cc @@ -0,0 +1,185 @@ +//--------------------------------------------------------------------------- +// $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 +#include +#include +#include +#include + + + +#ifdef DEAL_II_USE_TRILINOS + +DEAL_II_NAMESPACE_OPEN + +namespace TrilinosWrappers +{ + + PreconditionerTrilinosAmg::PreconditionerTrilinosAmg () + {} + + void PreconditionerTrilinosAmg::initialize ( + const dealii::SparseMatrix &matrix, + const std::vector &null_space, + const unsigned int null_space_dimension, + const bool elliptic, + const bool higher_order_elements, + const bool output_details, + const double drop_tolerance + ) + { + Assert (drop_tolerance >= 0, + ExcMessage ("Drop tolerance must be a non-negative number.")); + + const unsigned int n_rows = matrix.m(); + const SparsityPattern *sparsity_pattern = &(matrix.get_sparsity_pattern()); + + // Init Epetra Matrix, avoid + // storing the nonzero elements. + { + Map.reset (new Epetra_Map(n_rows, 0, communicator)); + + std::vector row_lengths (n_rows); + for (dealii::SparseMatrix::const_iterator p = matrix.begin(); + p != matrix.end(); ++p) + if (std::abs(p->value()) > drop_tolerance) + ++row_lengths[p->row()]; + + Matrix.reset (new Epetra_CrsMatrix(Copy, *Map, &row_lengths[0], true)); + + const unsigned int max_nonzero_entries + = *std::max_element (row_lengths.begin(), row_lengths.end()); + + std::vector values(max_nonzero_entries, 0); + std::vector row_indices(max_nonzero_entries); + + for (unsigned int row=0; row::const_iterator p = matrix.begin(row); + p != matrix.end(row); ++p) + if (std::abs(p->value()) > drop_tolerance) + { + row_indices[index] = p->column(); + values[index] = p->value(); + ++index; + } + + Assert (index == static_cast(row_lengths[row]), + ExcMessage("Filtering out zeros could not " + "be successfully finished!")); + + Matrix->InsertGlobalValues(row, row_lengths[row], + &values[0], &row_indices[0]); + } + + Matrix->FillComplete(); + } + + // Build the AMG preconditioner. + Teuchos::ParameterList parameter_list; + + // The implementation is able + // to distinguish between + // matrices from elliptic problems + // and convection dominated + // problems. We use the standard + // options 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). + if (elliptic) + { + 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); + } + + if (output_details) + parameter_list.set("ML output", 10); + else + parameter_list.set("ML output", 0); + + if (higher_order_elements) + parameter_list.set("aggregation: type", "MIS"); + + Assert (n_rows * null_space_dimension == null_space.size(), + ExcDimensionMismatch(n_rows * null_space_dimension, + null_space.size())); + + if (null_space_dimension > 1) + { + parameter_list.set("null space: type", "pre-computed"); + parameter_list.set("null space: dimension", int(null_space_dimension)); + parameter_list.set("null space: vectors", (double *)&null_space[0]); + } + + multigrid_operator = boost::shared_ptr + (new ML_Epetra::MultiLevelPreconditioner(*Matrix, parameter_list, true)); + + if (output_details) + multigrid_operator->PrintUnused(0); + } + + // 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. Luckily, it is sufficient + // 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 PreconditionerTrilinosAmg::vmult (dealii::Vector &dst, + const dealii::Vector &src) const + { + Epetra_Vector LHS (View, *Map, dst.begin()); + Epetra_Vector RHS (View, *Map, const_cast(src.begin())); + + const int res = multigrid_operator->ApplyInverse (RHS, LHS); + + Assert (res == 0, + ExcMessage ("Trilinos AMG MultiLevel preconditioner returned " + "with an error!")); + } + +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_USE_TRILINOS -- 2.39.5