From: kronbichler Date: Wed, 8 Oct 2008 13:22:13 +0000 (+0000) Subject: Update lac directory X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=30497265ad1bf23f3b6dd0e8e75638916b43ddbf;p=dealii-svn.git Update lac directory git-svn-id: https://svn.dealii.org/trunk@17147 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/trilinos_precondition_amg.h b/deal.II/lac/include/lac/trilinos_precondition_amg.h deleted file mode 100755 index c361615ba1..0000000000 --- a/deal.II/lac/include/lac/trilinos_precondition_amg.h +++ /dev/null @@ -1,235 +0,0 @@ -//--------------------------------------------------------------------------- -// $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 -#include -#include - -#include - - -#ifdef DEAL_II_USE_TRILINOS -#ifdef DEAL_II_COMPILER_SUPPORTS_MPI -# include -#else -# include -#endif -#include - -// some forward declarations -namespace ML_Epetra -{ - class MultiLevelPreconditioner; -} -class Epetra_Map; -class Epetra_CrsMatrix; - - - -DEAL_II_NAMESPACE_OPEN - -/*! @addtogroup TrilinosWrappers - *@{ - */ - -namespace TrilinosWrappers -{ - - // Forward declarations. - class SolverBase; - -/** - * 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 Subscriptor - { - public: - /** - * Constructor. - */ - PreconditionAMG (); - - /** - * 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 bool elliptic = true, - const bool higher_order_elements = false, - const double aggregation_threshold = 1e-4, - const std::vector > &null_space = std::vector > (1), - const bool output_details = false); - - /** - * 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 bool elliptic = true, - const bool higher_order_elements = false, - const double aggregation_threshold = 1e-4, - const std::vector > &null_space = std::vector > (), - const bool output_details = false); - - /** - * 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 to take the already - * generated coarsening - * structure, compute the AMG - * prolongation and restriction - * according to a smoothed - * aggregation strategy and - * then builds the whole - * multilevel hiearchy. This - * function can be considerably - * faster in that case, since - * the coarsening pattern is - * usually the most difficult - * thing to do when setting up - * the AMG ML preconditioner. - */ - void reinit (); - - /** - * Apply the preconditioner. - */ - void vmult (VectorBase &dst, - const VectorBase &src) const; - - /** - * Do the same as before, but - * now use deal.II vectors instead - * of the ones provided in the - * Trilinos wrapper class. - */ - void vmult (dealii::Vector &dst, - const dealii::Vector &src) const; - - /** - * Exception. - */ - DeclException1 (ExcNonMatchingMaps, - std::string, - << "The sparse matrix the preconditioner is based on " - << "uses a map that is not compatible to the one in vector" - << arg1 - << ". Check preconditioner and matrix setup."); - - private: - - /** - * A pointer to the - * preconditioner object. - */ - Teuchos::RefCountPtr - multilevel_operator; - - /** - * Internal communication - * pattern in case the matrix - * needs to be copied from - * deal.II format. - */ -#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. - */ - boost::shared_ptr Map; - - /** - * A copy of the deal.II matrix - * into Trilinos format. - */ - boost::shared_ptr Matrix; - - friend class SolverBase; - friend class SolverSaddlePoint; - }; -} - -/*@}*/ - - -DEAL_II_NAMESPACE_CLOSE - -#endif // DEAL_II_USE_TRILINOS - -/*---------------------------- trilinos_precondition_amg.h ---------------------------*/ - -#endif -/*---------------------------- trilinos_precondition_amg.h ---------------------------*/ diff --git a/deal.II/lac/source/trilinos_precondition_amg.cc b/deal.II/lac/source/trilinos_precondition_amg.cc deleted file mode 100755 index 29c33f0fe1..0000000000 --- a/deal.II/lac/source/trilinos_precondition_amg.cc +++ /dev/null @@ -1,217 +0,0 @@ -//--------------------------------------------------------------------------- -// $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 - -#ifdef DEAL_II_USE_TRILINOS - -#include -#include -#include -#include -#include - -DEAL_II_NAMESPACE_OPEN - -namespace TrilinosWrappers -{ - - PreconditionAMG::PreconditionAMG () -#ifdef DEAL_II_COMPILER_SUPPORTS_MPI - : - communicator (MPI_COMM_WORLD) -#endif - {} - - - - void - PreconditionAMG:: - initialize (const SparseMatrix &matrix, - const bool elliptic, - const bool higher_order_elements, - const double aggregation_threshold, - const std::vector > &null_space, - const bool output_details) - { - multilevel_operator.release(); - - const unsigned int n_rows = matrix.m(); - const unsigned int null_space_dimension = null_space.size(); - - // Build the AMG preconditioner. - Teuchos::ParameterList parameter_list; - - 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); - } - - parameter_list.set("aggregation: threshold", aggregation_threshold); - - 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"); - - 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 == null_space[0].size(), - ExcDimensionMismatch(n_rows, - 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)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 (output_details) - multilevel_operator->PrintUnused(0); - } - - - - void - PreconditionAMG:: - initialize (const ::dealii::SparseMatrix &deal_ii_sparse_matrix, - const bool elliptic, - const bool higher_order_elements, - const double aggregation_threshold, - const std::vector > &null_space, - const bool output_details) - { - 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); - Matrix->compress(); - - initialize (*Matrix, elliptic, higher_order_elements, - aggregation_threshold, null_space, output_details); - } - - - void PreconditionAMG:: - reinit () - { - multilevel_operator->ReComputePreconditioner(); - } - - - - void PreconditionAMG::vmult (VectorBase &dst, - const VectorBase &src) const - { - 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 - // 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 PreconditionAMG::vmult (dealii::Vector &dst, - const dealii::Vector &src) const - { - Assert (Map->SameAs(Matrix->matrix->RowMap()), - ExcNonMatchingMaps("dst")); - Assert (Map->SameAs(Matrix->matrix->RowMap()), - ExcNonMatchingMaps("src")); - - Epetra_Vector LHS (View, multilevel_operator->OperatorDomainMap(), - dst.begin()); - Epetra_Vector RHS (View, multilevel_operator->OperatorRangeMap(), - const_cast(src.begin())); - - const int res = multilevel_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