From: kronbichler Date: Thu, 24 Jun 2010 18:10:33 +0000 (+0000) Subject: Remove SolverBlock and PreconditionBlock from TrilinosWrappers to make deal.II compat... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e18b94fc2ef8214027733a58bf630fb57b0ef7b0;p=dealii-svn.git Remove SolverBlock and PreconditionBlock from TrilinosWrappers to make deal.II compatible with Trilinos v. 10.2. Use deal.II solvers with block matrices in that case. git-svn-id: https://svn.dealii.org/trunk@21329 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/aclocal.m4 b/deal.II/aclocal.m4 index 6d10cfb290..a407c939b2 100644 --- a/deal.II/aclocal.m4 +++ b/deal.II/aclocal.m4 @@ -6124,21 +6124,7 @@ AC_DEFUN(DEAL_II_CHECK_TRILINOS_HEADER_FILES, dnl Sacado.hpp \ Teuchos_ParameterList.hpp \ Teuchos_RCP.hpp \ - Teuchos_RefCountPtr.hpp \ - Thyra_AztecOOLinearOpWithSolveFactory.hpp \ - Thyra_DefaultBlockedLinearOp.hpp \ - Thyra_DefaultInverseLinearOp.hpp \ - Thyra_EpetraLinearOp.hpp \ - Thyra_EpetraThyraWrappers.hpp \ - Thyra_InverseLinearOperator.hpp \ - Thyra_LinearOperatorDecl.hpp \ - Thyra_LinearOperatorImpl.hpp \ - Thyra_LinearOpWithSolveFactoryHelpers.hpp \ - Thyra_MultiVectorBase.hpp \ - Thyra_MultiVectorDefaultBase.hpp \ - Thyra_VectorDecl.hpp \ - Thyra_VectorImpl.hpp \ - Thyra_VectorSpaceImpl.hpp + Teuchos_RefCountPtr.hpp ], [], [ diff --git a/deal.II/configure b/deal.II/configure index 756e4db6f8..79e39f2cdb 100755 --- a/deal.II/configure +++ b/deal.II/configure @@ -13825,21 +13825,7 @@ for ac_header in Amesos.h \ Sacado.hpp \ Teuchos_ParameterList.hpp \ Teuchos_RCP.hpp \ - Teuchos_RefCountPtr.hpp \ - Thyra_AztecOOLinearOpWithSolveFactory.hpp \ - Thyra_DefaultBlockedLinearOp.hpp \ - Thyra_DefaultInverseLinearOp.hpp \ - Thyra_EpetraLinearOp.hpp \ - Thyra_EpetraThyraWrappers.hpp \ - Thyra_InverseLinearOperator.hpp \ - Thyra_LinearOperatorDecl.hpp \ - Thyra_LinearOperatorImpl.hpp \ - Thyra_LinearOpWithSolveFactoryHelpers.hpp \ - Thyra_MultiVectorBase.hpp \ - Thyra_MultiVectorDefaultBase.hpp \ - Thyra_VectorDecl.hpp \ - Thyra_VectorImpl.hpp \ - Thyra_VectorSpaceImpl.hpp + Teuchos_RefCountPtr.hpp do as_ac_Header=`$as_echo "ac_cv_header_$ac_header" | $as_tr_sh` diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index b65e10b51c..b54602945a 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -31,6 +31,18 @@ inconvenience this causes.
  1. +

    + Removed: The class TrilinosWrappers::SolverBlock and + TrilinosWrappers::PreconditionBlock have been deleted for compatibility + reasons with Trilinos 10.2. These solvers can be exchanged by deal.II's + own iterative solvers with block matrices and vectors without loss in + performance. +
    + (Martin Kronbichler 2010/06/24) +

    + +
  2. +

    Changed: The SparseMatrix::operator() used to always check that an entry exists when accessing it. If it didn't, it would throw an exception that could be caught in a user program. Since these accesses are very frequent, diff --git a/deal.II/lac/include/lac/trilinos_precondition_block.h b/deal.II/lac/include/lac/trilinos_precondition_block.h deleted file mode 100644 index 65e6070455..0000000000 --- a/deal.II/lac/include/lac/trilinos_precondition_block.h +++ /dev/null @@ -1,89 +0,0 @@ -//--------------------------------------------------------------------------- -// $Id$ -// Version: $Name$ -// -// Copyright (C) 2008, 2009 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_block_h -#define __deal2__trilinos_precondition_block_h - - -#include - -#ifdef DEAL_II_USE_TRILINOS - -# include -# include -# include - -# include -# include -# include - -// some forward declarations -class Ifpack_Preconditioner; - - -DEAL_II_NAMESPACE_OPEN - -/*! @addtogroup TrilinosWrappers - *@{ - */ - -namespace TrilinosWrappers -{ - // forward declarations - class BlockSparseMatrix; - -/** - * This implements an interface class of preconditioners for block - * matrices based on Trilinos, which are intended to be used with the - * solver classes SolverBlockCG and SolverBlockGMRES. This is based on - * the Trilinos Thyra abstract interface of linear operators. This - * class does not implement any method, and derived classes will need - * to do that job. - * - * @ingroup TrilinosWrappers - * @ingroup Preconditioners - * @author Martin Kronbichler, 2008, 2009 - */ - class PreconditionBlockBase : public Subscriptor - { - public: - /** - * Empty constructor. - */ - PreconditionBlockBase(); - - /** - * Destructor. - */ - ~PreconditionBlockBase(); - - protected: - /** - * Preconditioner object. - */ - Teuchos::RCP > preconditioner; - - friend class SolverBlockBase; - }; -} - -/*@}*/ - - -DEAL_II_NAMESPACE_CLOSE - -#endif // DEAL_II_USE_TRILINOS - -/*-------------------------- trilinos_precondition_block.h ---------------------------*/ - -#endif -/*-------------------------- trilinos_precondition_block.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/trilinos_solver_block.h b/deal.II/lac/include/lac/trilinos_solver_block.h deleted file mode 100644 index 76cd35ec9a..0000000000 --- a/deal.II/lac/include/lac/trilinos_solver_block.h +++ /dev/null @@ -1,381 +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_block_solver_h -#define __deal2__trilinos_block_solver_h - - -#include - -#ifdef DEAL_II_USE_TRILINOS - -# include -# include - -DEAL_II_NAMESPACE_OPEN - - -namespace TrilinosWrappers -{ - // forward declarations - class BlockSparseMatrix; - class BlockVector; - namespace MPI - { - class BlockVector; - } - class SparseMatrix; - class PreconditionBase; - class PreconditionBlockBase; - - -/** - * Base class for solver classes using the Trilinos solvers on block - * matrices based on the Trilinos abstract solver interface Thyra. - * - * @ingroup TrilinosWrappers - * @author Martin Kronbichler, 2008 - */ - class SolverBlockBase - { - public: - - /** - * Enumeration object that is - * set in the constructor of - * the derived classes and - * tells Trilinos which solver - * to use. This option can also - * be set in the user program, - * so one might use this base - * class instead of one of the - * specialized derived classes - * when the solver should be - * set at runtime. Currently - * enabled options are: - */ - enum SolverBlockName {cg, gmres} solver_name; - - /** - * 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 bool right_preconditioning = false, - const bool output_details = false); - - /** - * The number of vectors that - * should be used in the GMRES - * algorithm before the solver - * is restarted. - */ - unsigned int gmres_restart_parameter; - - /** - * Flag that turns right - * preconditioning on. - */ - bool right_preconditioning; - - /** - * A flag to determine whether - * solver details should be - * written to screen. - */ - bool output_details; - }; - - /** - * Constructor. Takes the - * solver control object and - * creates the solver. - */ - SolverBlockBase (SolverControl &cn); - - /** - * Second constructor. This - * constructor takes an enum - * object that specifies the - * solver name and sets the - * appropriate Krylov - * method. - */ - SolverBlockBase (const enum SolverBlockName solver_name, - SolverControl &cn); - - /** - * Destructor. - */ - virtual ~SolverBlockBase (); - - /** - * 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. This interface is - * designed for Trilinos - * running in parallel. - */ - void - solve (const BlockSparseMatrix &A, - MPI::BlockVector &x, - const MPI::BlockVector &b, - const PreconditionBlockBase &preconditioner); - - /** - * 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. This interface is - * designed for Trilinos - * running in serial. - */ - void - solve (const BlockSparseMatrix &A, - BlockVector &x, - const BlockVector &b, - const PreconditionBlockBase &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"); - - /** - * Exception - */ - DeclException2 (ExcOverlappingMaps, - std::string, std::string, - << "The Epetra map in the " << arg1 << " " << arg2 - << " is overlapping between the processors, which" - << " is forbidden when doing a solve."); - - /** - * Exception. - */ - DeclException1 (ExcNonMatchingMaps, - std::string, - << "The Epetra map in the vector " - << arg1 - << " does not match the one used in the matrix. " - << "Check vector and matrix setup."); - - 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; - - protected: - - /** - * Store a copy of the flags for this - * particular solver. - */ - AdditionalData additional_data; - - }; - - - -/** - * An implementation of the solver interface using the Trilinos CG - * solver for block matrices and vectors. - * - * @ingroup TrilinosWrappers - * @author Martin Kronbichler, 2008 - */ - class SolverBlockCG : public SolverBlockBase - { - 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 bool right_preconditioning = false, - const bool output_details = false); - - /** - * Flag that turns right - * preconditioning on. - */ - bool right_preconditioning; - - /** - * A flag to determine whether - * solver details should be - * written to screen. - */ - bool output_details; - }; - - /** - * 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. - */ - SolverBlockCG (SolverControl &cn, - const AdditionalData &data = AdditionalData()); - - protected: - - /** - * Store a copy of the flags for this - * particular solver. - */ - const AdditionalData data; - }; - - - -/** - * An implementation of the solver interface using the Trilinos GMRES - * solver for block matrices and vectors. - * - * @author Martin Kronbichler, 2008 - */ - class SolverBlockGMRES : public SolverBlockBase - { - 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 bool right_preconditioning = false, - const bool output_details = false); - - /** - * The number of vectors that - * should be used in the GMRES - * algorithm before the solver - * is restarted. - */ - unsigned int gmres_restart_parameter; - - /** - * Flag that turns right - * preconditioning on. - */ - bool right_preconditioning; - - /** - * A flag to determine whether - * solver details should be - * written to screen. - */ - bool output_details; - }; - - /** - * 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. - */ - SolverBlockGMRES (SolverControl &cn, - const AdditionalData &data = AdditionalData()); - - protected: - - /** - * Store a copy of the flags for this - * particular solver. - */ - const AdditionalData data; - }; - - -} - -DEAL_II_NAMESPACE_CLOSE - -#endif // DEAL_II_USE_TRILINOS - -/*---------------------- trilinos_solver_block.h -----------------------*/ - -#endif -/*---------------------- trilinos_solver_block.h -----------------------*/ diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h index 397f344b0d..773acec024 100644 --- a/deal.II/lac/include/lac/vector.templates.h +++ b/deal.II/lac/include/lac/vector.templates.h @@ -330,7 +330,7 @@ Number Vector::operator * (const Vector& v) const const unsigned int blocking = 1<>(BLOCK_LEVEL))<<(BLOCK_LEVEL)), *X_end2 = X + ((vec_size>>(2*BLOCK_LEVEL))<<(2*BLOCK_LEVEL)), *X_end1 = X + ((vec_size>>(3*BLOCK_LEVEL))<<(3*BLOCK_LEVEL)); @@ -341,25 +341,25 @@ Number Vector::operator * (const Vector& v) const // the result vector. this is necessary // because // operator*(complex,complex) - // is not defined by default. do the + // is not defined by default. do the // operations block-wise with post-update. - // use three nested loops, which should make - // the roundoff error very small up to + // use three nested loops, which should make + // the roundoff error very small up to // about 20m entries. in the - // end do extra loops with the remainders. + // end do extra loops with the remainders. // this blocked algorithm has been proposed // by Castaldo, Whaley and Chronopoulos - // (SIAM J. Sci. Comput. 31, 1156-1174, + // (SIAM J. Sci. Comput. 31, 1156-1174, // 2008) while (X != X_end1) { - sum1 = 0.; + sum1 = Number(); for (unsigned int j=0; j::conjugate(*Y++)); sum2 += sum3; @@ -370,10 +370,10 @@ Number Vector::operator * (const Vector& v) const } while (X != X_end2) { - sum2 = 0.; + sum2 = Number(); for (unsigned int j=0; j::conjugate(*Y++)); sum2 += sum3; @@ -382,14 +382,14 @@ Number Vector::operator * (const Vector& v) const } while (X != X_end3) { - sum3 = 0; + sum3 = Number(); for (unsigned int i=0; i::conjugate(*Y++)); sum += sum3; } - sum3 = 0; - while (X != X_end) + sum3 = Number(); + while (X != X_end) sum3 += *X++ * Number(numbers::NumberTraits::conjugate(*Y++)); sum += sum3; @@ -405,7 +405,7 @@ Vector::norm_sqr () const const unsigned int blocking = 1<>(BLOCK_LEVEL))<<(BLOCK_LEVEL)), *X_end2 = X + ((vec_size>>(2*BLOCK_LEVEL))<<(2*BLOCK_LEVEL)), *X_end1 = X + ((vec_size>>(3*BLOCK_LEVEL))<<(3*BLOCK_LEVEL)); @@ -415,10 +415,10 @@ Vector::norm_sqr () const sum1 = 0.; for (unsigned int j=0; j::abs_square(*X++); sum2 += sum3; @@ -432,7 +432,7 @@ Vector::norm_sqr () const sum2 = 0.; for (unsigned int j=0; j::abs_square(*X++); sum2 += sum3; @@ -441,14 +441,14 @@ Vector::norm_sqr () const } while (X != X_end3) { - sum3 = 0; + sum3 = 0.; for (unsigned int i=0; i::abs_square(*X++); sum += sum3; } - sum3 = 0; - while (X != X_end) + sum3 = 0.; + while (X != X_end) sum3 += numbers::NumberTraits::abs_square(*X++); sum += sum3; @@ -462,21 +462,21 @@ Number Vector::mean_value () const Assert (vec_size!=0, ExcEmptyObject()); const unsigned int blocking = 1<>(BLOCK_LEVEL))<<(BLOCK_LEVEL)), *X_end2 = X + ((vec_size>>(2*BLOCK_LEVEL))<<(2*BLOCK_LEVEL)), *X_end1 = X + ((vec_size>>(3*BLOCK_LEVEL))<<(3*BLOCK_LEVEL)); while (X != X_end1) { - sum1 = 0.; + sum1 = Number(); for (unsigned int j=0; j::mean_value () const } while (X != X_end2) { - sum2 = 0.; + sum2 = Number(); for (unsigned int j=0; j::mean_value () const } while (X != X_end3) { - sum3 = 0; + sum3 = Number(); for (unsigned int i=0; i::real_type Vector::l1_norm () const { Assert (vec_size!=0, ExcEmptyObject()); - + const unsigned int blocking = 1<>(BLOCK_LEVEL))<<(BLOCK_LEVEL)), *X_end2 = X + ((vec_size>>(2*BLOCK_LEVEL))<<(2*BLOCK_LEVEL)), *X_end1 = X + ((vec_size>>(3*BLOCK_LEVEL))<<(3*BLOCK_LEVEL)); @@ -533,10 +533,10 @@ Vector::l1_norm () const sum1 = 0.; for (unsigned int j=0; j::abs(*X++); sum2 += sum3; @@ -550,7 +550,7 @@ Vector::l1_norm () const sum2 = 0.; for (unsigned int j=0; j::abs(*X++); sum2 += sum3; @@ -559,14 +559,14 @@ Vector::l1_norm () const } while (X != X_end3) { - sum3 = 0; + sum3 = 0.; for (unsigned int i=0; i::abs(*X++); sum += sum3; } - sum3 = 0; - while (X != X_end) + sum3 = 0.; + while (X != X_end) sum3 += numbers::NumberTraits::abs(*X++); sum += sum3; @@ -597,7 +597,7 @@ Vector::lp_norm (const real_type p) const const unsigned int blocking = 1<>(BLOCK_LEVEL))<<(BLOCK_LEVEL)), *X_end2 = X + ((vec_size>>(2*BLOCK_LEVEL))<<(2*BLOCK_LEVEL)), *X_end1 = X + ((vec_size>>(3*BLOCK_LEVEL))<<(3*BLOCK_LEVEL)); @@ -607,10 +607,10 @@ Vector::lp_norm (const real_type p) const sum1 = 0.; for (unsigned int j=0; j::abs(*X++),p); sum2 += sum3; @@ -624,7 +624,7 @@ Vector::lp_norm (const real_type p) const sum2 = 0.; for (unsigned int j=0; j::abs(*X++),p); sum2 += sum3; @@ -633,14 +633,14 @@ Vector::lp_norm (const real_type p) const } while (X != X_end3) { - sum3 = 0; + sum3 = 0.; for (unsigned int i=0; i::abs(*X++),p); sum += sum3; } - sum3 = 0; - while (X != X_end) + sum3 = 0.; + while (X != X_end) sum3 += std::pow(numbers::NumberTraits::abs(*X++),p); sum += sum3; @@ -655,7 +655,7 @@ Vector::linfty_norm () const { Assert (vec_size!=0, ExcEmptyObject()); - real_type max = 0; + real_type max = 0.; for (unsigned int i=0; i::abs(val[i]), max); diff --git a/deal.II/lac/source/lapack_full_matrix.cc b/deal.II/lac/source/lapack_full_matrix.cc index 5d37b9a90c..74bede304e 100644 --- a/deal.II/lac/source/lapack_full_matrix.cc +++ b/deal.II/lac/source/lapack_full_matrix.cc @@ -318,10 +318,12 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric ( { Assert(state == matrix, ExcState(state)); const int nn = this->n_cols(); - Assert(nn == this->n_rows(), ExcNotQuadratic()); + Assert(static_cast(nn) == this->n_rows(), ExcNotQuadratic()); Assert(B.n_rows() == B.n_cols(), ExcNotQuadratic()); - Assert(nn == B.n_cols(), ExcDimensionMismatch (nn, B.n_cols())); - Assert(eigenvectors.size() <= nn, ExcMessage ("eigenvectors.size() > matrix.n_cols()")); + Assert(static_cast(nn) == B.n_cols(), + ExcDimensionMismatch (nn, B.n_cols())); + Assert(eigenvectors.size() <= static_cast(nn), + ExcMessage ("eigenvectors.size() > matrix.n_cols()")); wr.resize(nn); wi.resize(nn); //This is set purley for consistency reasons with the @@ -393,7 +395,7 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric ( { unsigned int col_begin(i*nn); eigenvectors[i].reinit(nn, true); - for (unsigned int j=0; j < nn; ++j) + for (unsigned int j=0; j < static_cast(nn); ++j) { eigenvectors[i](j) = values_A[col_begin+j]; } diff --git a/deal.II/lac/source/trilinos_precondition.cc b/deal.II/lac/source/trilinos_precondition.cc index 8aa2df9364..d654172101 100644 --- a/deal.II/lac/source/trilinos_precondition.cc +++ b/deal.II/lac/source/trilinos_precondition.cc @@ -513,8 +513,8 @@ namespace TrilinosWrappers PreconditionAMG:: initialize (const SparseMatrix &matrix, const AdditionalData &additional_data) { - preconditioner.release(); - multilevel_operator.release(); + preconditioner.reset (); + multilevel_operator.reset (); const unsigned int n_rows = matrix.m(); @@ -632,10 +632,9 @@ namespace TrilinosWrappers PreconditionAMG::initialize (const SparseMatrix &matrix, const Teuchos::ParameterList &ml_parameters) { - multilevel_operator = Teuchos::rcp (new ML_Epetra::MultiLevelPreconditioner( - matrix.trilinos_matrix(), ml_parameters, - true)); - preconditioner = Teuchos::rcp (multilevel_operator.get(), false); + multilevel_operator.reset (new ML_Epetra::MultiLevelPreconditioner + (matrix.trilinos_matrix(), ml_parameters)); + preconditioner.reset (multilevel_operator.get(), false); } @@ -648,6 +647,8 @@ namespace TrilinosWrappers const double drop_tolerance, const ::dealii::SparsityPattern *use_this_sparsity) { + preconditioner.reset(); + multilevel_operator.reset (); const unsigned int n_rows = deal_ii_sparse_matrix.m(); // Init Epetra Matrix using an diff --git a/deal.II/lac/source/trilinos_precondition_block.cc b/deal.II/lac/source/trilinos_precondition_block.cc deleted file mode 100644 index ece2a4a611..0000000000 --- a/deal.II/lac/source/trilinos_precondition_block.cc +++ /dev/null @@ -1,40 +0,0 @@ -//--------------------------------------------------------------------------- -// $Id$ -// Version: $Name$ -// -// Copyright (C) 2008, 2009 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 - - -DEAL_II_NAMESPACE_OPEN - -namespace TrilinosWrappers -{ - - PreconditionBlockBase::PreconditionBlockBase() - {} - - - - PreconditionBlockBase::~PreconditionBlockBase() - {} - -} /* end of namespace TrilinosWrappers */ - -DEAL_II_NAMESPACE_CLOSE - -#endif // DEAL_II_USE_TRILINOS diff --git a/deal.II/lac/source/trilinos_solver_block.cc b/deal.II/lac/source/trilinos_solver_block.cc deleted file mode 100644 index cbae827f56..0000000000 --- a/deal.II/lac/source/trilinos_solver_block.cc +++ /dev/null @@ -1,429 +0,0 @@ -//--------------------------------------------------------------------------- -// $Id$ -// Version: $Name$ -// -// Copyright (C) 2008, 2009 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 -# include - -# include - -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include - -DEAL_II_NAMESPACE_OPEN - -namespace TrilinosWrappers -{ - - SolverBlockBase::AdditionalData:: - AdditionalData (const unsigned int gmres_restart_parameter, - const bool right_preconditioning, - const bool output_details) - : - gmres_restart_parameter (gmres_restart_parameter), - right_preconditioning (right_preconditioning), - output_details (output_details) - {} - - - - SolverBlockBase::SolverBlockBase (SolverControl &cn) - : - solver_name (gmres), - solver_control (cn) - {} - - - - SolverBlockBase::SolverBlockBase (SolverBlockBase::SolverBlockName solver_name, - SolverControl &cn) - : - solver_name (solver_name), - solver_control (cn) - {} - - - - SolverBlockBase::~SolverBlockBase () - {} - - - - SolverControl & - SolverBlockBase::control() const - { - return solver_control; - } - - - - void - SolverBlockBase::solve (const BlockSparseMatrix &A, - BlockVector &x, - const BlockVector &b, - const PreconditionBlockBase &preconditioner) - { - - Assert (x.n_blocks() == b.n_blocks(), - ExcDimensionMismatch(x.n_blocks(), b.n_blocks())); - - // Just copy everything over - // to a distributed vector... - MPI::BlockVector tmpx (x.n_blocks()); - MPI::BlockVector tmpb (b.n_blocks()); - for (unsigned int i=0; i P = *preconditioner.preconditioner; - const unsigned int n_rows = P.numBlockRows(); - const unsigned int n_cols = P.numBlockCols(); - - // Check the dimensions of the - // block vectors and matrix - // and check if the parallel - // distribution of the - // matrices and vectors is the - // same. - AssertThrow (n_rows == n_cols, - ExcDimensionMismatch (n_rows, n_cols)); - - AssertThrow (input_x.n_blocks() == input_A.n_block_cols(), - ExcDimensionMismatch(input_x.n_blocks(), - input_A.n_block_cols())); - AssertThrow (input_b.n_blocks() == input_A.n_block_rows(), - ExcDimensionMismatch(input_b.n_blocks(), - input_A.n_block_rows())); - AssertThrow (input_A.n_block_rows() == n_rows, - ExcDimensionMismatch(input_A.n_block_rows(), n_rows)); - AssertThrow (input_A.n_block_cols() == n_cols, - ExcDimensionMismatch(input_A.n_block_cols(), n_cols)); - - for (unsigned int i=0; i > > - epetra_vector_spaces; - for (unsigned int i=0; i > tmp_space - = Thyra::create_VectorSpace( - Teuchos::rcp(&input_A.block(i,i).domain_partitioner(), - false)); - - epetra_vector_spaces.push_back(tmp_space); - } - - // Convert the block matrix to - // a Thyra linear operator - // that acts on blocks. - Teuchos::RCP > - tmpA = Teuchos::rcp(new Thyra::DefaultBlockedLinearOp()); - tmpA->beginBlockFill(n_rows,n_cols); - for (unsigned int i=0; i > - A_ij = Thyra::epetraLinearOp(Teuchos::rcp(&input_A.block(i,j).trilinos_matrix(), - false)); - tmpA->setBlock(i, j, A_ij); - } - tmpA->endBlockFill(); - tmpA->setObjectLabel("Blocked system matrix"); - Thyra::RCP > baseA = tmpA; - - Thyra::ConstLinearOperator A = baseA; - - // Convert block vectors for - // right hand side and - // solution to Thyra objects. - Thyra::Vector rhs = A.range().createMember(); - Thyra::Vector sol = A.domain().createMember(); - for (unsigned int i=0; i - ((input_b.block(i).trilinos_vector())(0)); - Teuchos::RCP > tmp_rhs_i - = Thyra::create_Vector(Teuchos::rcp(block_ptr, false), - epetra_vector_spaces[i]); - const Thyra::Vector rhs_i = tmp_rhs_i; - rhs.setBlock(i, rhs_i); - - block_ptr = (input_x.block(i).trilinos_vector())(0); - Teuchos::RCP > tmp_sol_i - = Thyra::create_Vector(Teuchos::rcp(block_ptr, false), - epetra_vector_spaces[i]); - Thyra::Vector sol_i = tmp_sol_i; - sol.setBlock(i, sol_i); - } - - // ------- Now build a solver factory for the block problem ------- - - // Set up a parameter list for - // the AztecOO solver. - Teuchos::RCP aztecBlockParams - = rcp(new Teuchos::ParameterList("AztecOO block solver parameters")); - - // Set the number of - // iterations and the solver - // tolerance from the solver - // control object. - aztecBlockParams->sublist("Forward Solve").set("Max Iterations", - (int)solver_control.max_steps()); - aztecBlockParams->sublist("Forward Solve").set("Tolerance", - solver_control.tolerance()); - aztecBlockParams->sublist("Forward Solve") - .sublist("AztecOO Settings").set("Convergence Test", "no scaling"); - - // The settings about the - // actual solver that is going - // to be used - let this be - // determined by the - // solver_name set in - // a derived class. - switch (solver_name) - { - case cg: - aztecBlockParams->sublist("Forward Solve") - .sublist("AztecOO Settings").set("Aztec Solver", "CG"); - break; - case gmres: - aztecBlockParams->sublist("Forward Solve") - .sublist("AztecOO Settings").set("Aztec Solver", "GMRES"); - aztecBlockParams->sublist("Forward Solve") - .sublist("AztecOO Settings").set("Size of Krylov Subspace", - (int)additional_data.gmres_restart_parameter); - break; - default: - ExcNotImplemented(); - } - - // We use an externally - // provided preconditioner, so - // do not set anything here. - aztecBlockParams->sublist("Forward Solve") - .sublist("AztecOO Settings").set("Aztec Preconditioner", "none"); - - if (additional_data.output_details == true) - aztecBlockParams->sublist("Forward Solve") - .sublist("AztecOO Settings").set("Output Frequency", 1); - else - aztecBlockParams->sublist("VerboseObject").set("Verbosity Level", - "none"); - - // Create solve factory from - // parameter list. - Teuchos::RCP > - aztecBlockLowsFactory = Teuchos::rcp - (new Thyra::AztecOOLinearOpWithSolveFactory()); - - aztecBlockLowsFactory->setParameterList(aztecBlockParams); - - Teuchos::RCP > rcpBlockAztec - = aztecBlockLowsFactory->createOp(); - - - // ---- Set up the preconditioned inverse object and do the solve ---- - - // Extract base operators from - // block matrix A and block - // preconditioner P. - Teuchos::RCP > A_ptr - = Teuchos::rcp(new Thyra::DefaultLinearOpSource(A.constPtr()), - true); - Teuchos::RCP > P_ptr = - Teuchos::rcp (new Thyra::DefaultPreconditioner(), true); - { - Thyra::DefaultPreconditioner - *defaultPrec = &Teuchos::dyn_cast > - (*P_ptr); - if (additional_data.right_preconditioning == true) - (*defaultPrec).initializeRight(P.constPtr()); - else - (*defaultPrec).initializeUnspecified(P.constPtr()); - } - - // Attach the matrix and the - // preconditioner to the - // solver object. - aztecBlockLowsFactory->initializePreconditionedOp(A_ptr, P_ptr, - &*rcpBlockAztec, - Thyra::SUPPORT_SOLVE_FORWARD_ONLY); - - // Now do the solve. - Thyra::SolveStatus - solveStatus = Thyra::solve( *rcpBlockAztec, Thyra::NOTRANS, - *rhs.constPtr().get(), &*sol.ptr().get() ); - - // Extract the number of - // needed iterations from the - // solver message. - unsigned int n_iterations; - { - std::string output_text = solveStatus.message; - unsigned int position = output_text.find ("using "); - const std::pair tmp - = Utilities::get_integer_at_position (output_text, position+6); - n_iterations = tmp.first; - } - - // Do a dummy check of - // convergence to tell the - // SolverControl object the - // number of iterations. We do - // not impose the achieved - // tolerance here since that - // would require additional - // work (computing a - // residual). - solver_control.check (n_iterations, 0); - - } - - - - - -/* ------------------------ SolverCG -------------------------- */ - - SolverBlockCG::AdditionalData:: - AdditionalData (const bool right_preconditioning, - const bool output_details) - : - right_preconditioning (right_preconditioning), - output_details (output_details) - {} - - - - SolverBlockCG::SolverBlockCG (SolverControl &cn, - const AdditionalData &data) - : - SolverBlockBase (cn), - data (data) - { - solver_name = cg; - additional_data.right_preconditioning = data.right_preconditioning; - additional_data.output_details = data.output_details; - } - - -/* ---------------------- SolverGMRES ------------------------- */ - - SolverBlockGMRES::AdditionalData:: - AdditionalData (const unsigned int gmres_restart_parameter, - const bool right_preconditioning, - const bool output_details) - : - gmres_restart_parameter (gmres_restart_parameter), - right_preconditioning (right_preconditioning), - output_details (output_details) - {} - - - - SolverBlockGMRES::SolverBlockGMRES (SolverControl &cn, - const AdditionalData &data) - : - SolverBlockBase (cn), - data (data) - { - solver_name = gmres; - additional_data.gmres_restart_parameter = data.gmres_restart_parameter; - additional_data.right_preconditioning = data.right_preconditioning; - additional_data.output_details = data.output_details; - } - - -} - -DEAL_II_NAMESPACE_CLOSE - -#endif // DEAL_II_USE_PETSC