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
],
[],
[
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`
<ol>
<li>
+ <p>
+ 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.
+ <br>
+ (Martin Kronbichler 2010/06/24)
+ </p>
+
+ <li>
+ <p>
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,
+++ /dev/null
-//---------------------------------------------------------------------------
-// $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 <base/config.h>
-
-#ifdef DEAL_II_USE_TRILINOS
-
-# include <lac/exceptions.h>
-# include <lac/trilinos_precondition.h>
-# include <lac/trilinos_sparse_matrix.h>
-
-# include <Teuchos_RCP.hpp>
-# include <Thyra_LinearOperatorDecl.hpp>
-# include <Epetra_Operator.h>
-
-// 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<Thyra::ConstLinearOperator<double> > preconditioner;
-
- friend class SolverBlockBase;
- };
-}
-
-/*@}*/
-
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif // DEAL_II_USE_TRILINOS
-
-/*-------------------------- trilinos_precondition_block.h ---------------------------*/
-
-#endif
-/*-------------------------- trilinos_precondition_block.h ---------------------------*/
+++ /dev/null
-//---------------------------------------------------------------------------
-// $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 <base/config.h>
-
-#ifdef DEAL_II_USE_TRILINOS
-
-# include <lac/exceptions.h>
-# include <lac/solver_control.h>
-
-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
- * <tt>Ax=b</tt>. 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
- * <tt>Ax=b</tt>. 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 -----------------------*/
const unsigned int blocking = 1<<BLOCK_LEVEL;
Number sum1, sum2, sum3, sum = Number();
- const Number * X = val, *X_end = X + vec_size,
+ const Number * X = val, *X_end = X + vec_size,
*X_end3 = X + ((vec_size>>(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));
// the result vector. this is necessary
// because
// operator*(complex<float>,complex<double>)
- // 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<blocking; ++j)
{
- sum2 = 0;
+ sum2 = Number();
for (unsigned int k=0; k<blocking; ++k)
{
- sum3 = 0;
+ sum3 = Number();
for (unsigned int i=0; i<blocking; ++i)
sum3 += *X++ * Number(numbers::NumberTraits<Number2>::conjugate(*Y++));
sum2 += sum3;
}
while (X != X_end2)
{
- sum2 = 0.;
+ sum2 = Number();
for (unsigned int j=0; j<blocking; ++j)
{
- sum3 = 0;
+ sum3 = Number();
for (unsigned int i=0; i<blocking; ++i)
sum3 += *X++ * Number(numbers::NumberTraits<Number2>::conjugate(*Y++));
sum2 += sum3;
}
while (X != X_end3)
{
- sum3 = 0;
+ sum3 = Number();
for (unsigned int i=0; i<blocking; ++i)
sum3 += *X++ * Number(numbers::NumberTraits<Number2>::conjugate(*Y++));
sum += sum3;
}
- sum3 = 0;
- while (X != X_end)
+ sum3 = Number();
+ while (X != X_end)
sum3 += *X++ * Number(numbers::NumberTraits<Number2>::conjugate(*Y++));
sum += sum3;
const unsigned int blocking = 1<<BLOCK_LEVEL;
real_type sum1, sum2, sum3, sum = 0.;
- const Number * X = val, *X_end = X + vec_size,
+ const Number * X = val, *X_end = X + vec_size,
*X_end3 = X + ((vec_size>>(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));
sum1 = 0.;
for (unsigned int j=0; j<blocking; ++j)
{
- sum2 = 0;
+ sum2 = 0.;
for (unsigned int k=0; k<blocking; ++k)
{
- sum3 = 0;
+ sum3 = 0.;
for (unsigned int i=0; i<blocking; ++i)
sum3 += numbers::NumberTraits<Number>::abs_square(*X++);
sum2 += sum3;
sum2 = 0.;
for (unsigned int j=0; j<blocking; ++j)
{
- sum3 = 0;
+ sum3 = 0.;
for (unsigned int i=0; i<blocking; ++i)
sum3 += numbers::NumberTraits<Number>::abs_square(*X++);
sum2 += sum3;
}
while (X != X_end3)
{
- sum3 = 0;
+ sum3 = 0.;
for (unsigned int i=0; i<blocking; ++i)
sum3 += numbers::NumberTraits<Number>::abs_square(*X++);
sum += sum3;
}
- sum3 = 0;
- while (X != X_end)
+ sum3 = 0.;
+ while (X != X_end)
sum3 += numbers::NumberTraits<Number>::abs_square(*X++);
sum += sum3;
Assert (vec_size!=0, ExcEmptyObject());
const unsigned int blocking = 1<<BLOCK_LEVEL;
- Number sum1, sum2, sum3, sum = 0.;
- const Number * X = val, *X_end = X + vec_size,
+ Number sum1, sum2, sum3, sum = Number();
+ const Number * X = val, *X_end = X + vec_size,
*X_end3 = X + ((vec_size>>(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<blocking; ++j)
{
- sum2 = 0;
+ sum2 = Number();
for (unsigned int k=0; k<blocking; ++k)
{
- sum3 = 0;
+ sum3 = Number();
for (unsigned int i=0; i<blocking; ++i)
sum3 += *X++;
sum2 += sum3;
}
while (X != X_end2)
{
- sum2 = 0.;
+ sum2 = Number();
for (unsigned int j=0; j<blocking; ++j)
{
- sum3 = 0;
+ sum3 = Number();
for (unsigned int i=0; i<blocking; ++i)
sum3 += *X++;
sum2 += sum3;
}
while (X != X_end3)
{
- sum3 = 0;
+ sum3 = Number();
for (unsigned int i=0; i<blocking; ++i)
sum3 += *X++;
sum += sum3;
}
- sum3 = 0;
- while (X != X_end)
+ sum3 = Number();
+ while (X != X_end)
sum3 += *X++;
sum += sum3;
Vector<Number>::l1_norm () const
{
Assert (vec_size!=0, ExcEmptyObject());
-
+
const unsigned int blocking = 1<<BLOCK_LEVEL;
real_type sum1, sum2, sum3, sum = 0.;
- const Number * X = val, *X_end = X + vec_size,
+ const Number * X = val, *X_end = X + vec_size,
*X_end3 = X + ((vec_size>>(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));
sum1 = 0.;
for (unsigned int j=0; j<blocking; ++j)
{
- sum2 = 0;
+ sum2 = 0.;
for (unsigned int k=0; k<blocking; ++k)
{
- sum3 = 0;
+ sum3 = 0.;
for (unsigned int i=0; i<blocking; ++i)
sum3 += numbers::NumberTraits<Number>::abs(*X++);
sum2 += sum3;
sum2 = 0.;
for (unsigned int j=0; j<blocking; ++j)
{
- sum3 = 0;
+ sum3 = 0.;
for (unsigned int i=0; i<blocking; ++i)
sum3 += numbers::NumberTraits<Number>::abs(*X++);
sum2 += sum3;
}
while (X != X_end3)
{
- sum3 = 0;
+ sum3 = 0.;
for (unsigned int i=0; i<blocking; ++i)
sum3 += numbers::NumberTraits<Number>::abs(*X++);
sum += sum3;
}
- sum3 = 0;
- while (X != X_end)
+ sum3 = 0.;
+ while (X != X_end)
sum3 += numbers::NumberTraits<Number>::abs(*X++);
sum += sum3;
const unsigned int blocking = 1<<BLOCK_LEVEL;
real_type sum1, sum2, sum3, sum = 0.;
- const Number * X = val, *X_end = X + vec_size,
+ const Number * X = val, *X_end = X + vec_size,
*X_end3 = X + ((vec_size>>(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));
sum1 = 0.;
for (unsigned int j=0; j<blocking; ++j)
{
- sum2 = 0;
+ sum2 = 0.;
for (unsigned int k=0; k<blocking; ++k)
{
- sum3 = 0;
+ sum3 = 0.;
for (unsigned int i=0; i<blocking; ++i)
sum3 += std::pow(numbers::NumberTraits<Number>::abs(*X++),p);
sum2 += sum3;
sum2 = 0.;
for (unsigned int j=0; j<blocking; ++j)
{
- sum3 = 0;
+ sum3 = 0.;
for (unsigned int i=0; i<blocking; ++i)
sum3 += std::pow(numbers::NumberTraits<Number>::abs(*X++),p);
sum2 += sum3;
}
while (X != X_end3)
{
- sum3 = 0;
+ sum3 = 0.;
for (unsigned int i=0; i<blocking; ++i)
sum3 += std::pow(numbers::NumberTraits<Number>::abs(*X++),p);
sum += sum3;
}
- sum3 = 0;
- while (X != X_end)
+ sum3 = 0.;
+ while (X != X_end)
sum3 += std::pow(numbers::NumberTraits<Number>::abs(*X++),p);
sum += sum3;
{
Assert (vec_size!=0, ExcEmptyObject());
- real_type max = 0;
+ real_type max = 0.;
for (unsigned int i=0; i<vec_size; ++i)
max = std::max (numbers::NumberTraits<Number>::abs(val[i]), max);
{
Assert(state == matrix, ExcState(state));
const int nn = this->n_cols();
- Assert(nn == this->n_rows(), ExcNotQuadratic());
+ Assert(static_cast<unsigned int>(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<unsigned int>(nn) == B.n_cols(),
+ ExcDimensionMismatch (nn, B.n_cols()));
+ Assert(eigenvectors.size() <= static_cast<unsigned int>(nn),
+ ExcMessage ("eigenvectors.size() > matrix.n_cols()"));
wr.resize(nn);
wi.resize(nn); //This is set purley for consistency reasons with the
{
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<unsigned int>(nn); ++j)
{
eigenvectors[i](j) = values_A[col_begin+j];
}
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();
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);
}
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
+++ /dev/null
-//---------------------------------------------------------------------------
-// $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 <lac/trilinos_precondition_block.h>
-
-#ifdef DEAL_II_USE_TRILINOS
-
-# include <lac/trilinos_sparse_matrix.h>
-# include <lac/trilinos_block_sparse_matrix.h>
-
-
-DEAL_II_NAMESPACE_OPEN
-
-namespace TrilinosWrappers
-{
-
- PreconditionBlockBase::PreconditionBlockBase()
- {}
-
-
-
- PreconditionBlockBase::~PreconditionBlockBase()
- {}
-
-} /* end of namespace TrilinosWrappers */
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif // DEAL_II_USE_TRILINOS
+++ /dev/null
-//---------------------------------------------------------------------------
-// $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 <lac/trilinos_solver_block.h>
-
-#ifdef DEAL_II_USE_TRILINOS
-
-# include <lac/trilinos_block_sparse_matrix.h>
-# include <lac/trilinos_block_vector.h>
-# include <lac/trilinos_precondition_block.h>
-
-# include <base/utilities.h>
-
-# include <Teuchos_RCP.hpp>
-# include <Teuchos_ParameterList.hpp>
-
-# include <Epetra_Vector.h>
-
-# include <Thyra_VectorDecl.hpp>
-# include <Thyra_VectorImpl.hpp>
-# include <Thyra_VectorSpaceImpl.hpp>
-# include <Thyra_MultiVectorBase.hpp>
-# include <Thyra_MultiVectorDefaultBase.hpp>
-# include <Thyra_LinearOperatorDecl.hpp>
-# include <Thyra_LinearOperatorImpl.hpp>
-# include <Thyra_DefaultBlockedLinearOp.hpp>
-# include <Thyra_EpetraLinearOp.hpp>
-# include <Thyra_EpetraThyraWrappers.hpp>
-# include <Thyra_LinearOpWithSolveFactoryHelpers.hpp>
-# include <Thyra_AztecOOLinearOpWithSolveFactory.hpp>
-
-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<x.n_blocks(); ++i)
- {
- tmpx.block(i).reinit (A.block(i,i).domain_partitioner(), false);
- tmpx.block(i) = x.block(i);
- tmpb.block(i).reinit (A.block(i,i).range_partitioner(), false);
- tmpb.block(i) = b.block(i);
- }
- tmpx.collect_sizes();
- tmpb.collect_sizes();
-
- solve (A, tmpx, tmpb, preconditioner);
- x = tmpx;
- }
-
-
-
- void
- SolverBlockBase::solve (const BlockSparseMatrix &input_A,
- MPI::BlockVector &input_x,
- const MPI::BlockVector &input_b,
- const PreconditionBlockBase &preconditioner)
- {
-
- // Get block information about
- // the preconditioner.
- Thyra::ConstLinearOperator<double> 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<n_rows; ++i)
- {
- AssertThrow (input_x.block(i).vector_partitioner().UniqueGIDs() == true,
- ExcOverlappingMaps("vector", "x"));
- AssertThrow (input_b.block(i).vector_partitioner().UniqueGIDs() == true,
- ExcOverlappingMaps("vector", "b"));
- for (unsigned int j=0; j<n_cols; ++j)
- {
- {
- std::ostringstream i_str;
- i_str << i;
- std::ostringstream error_component;
- error_component << "x.block(" << i_str.str() << ")";
- AssertThrow (input_x.block(j).vector_partitioner().SameAs(
- input_A.block(i,j).domain_partitioner()) == true,
- ExcNonMatchingMaps (error_component.str()));
- }
- {
- std::ostringstream i_str;
- i_str << j;
- std::ostringstream error_component;
- error_component << "b.block(" << i_str.str() << ")";
- AssertThrow (input_b.block(i).vector_partitioner().SameAs(
- input_A.block(i,j).range_partitioner()) == true,
- ExcNonMatchingMaps (error_component.str()));
- }
- }
- }
-
-
- // Wrap Epetra vectors into
- // Thyra vectors in order to
- // be able to work with the
- // abstract Thyra/Stratimikos
- // interfaces to the AztecOO
- // solver.
- //
- // As a first step, create a
- // copy of the vector spaces.
- std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > >
- epetra_vector_spaces;
- for (unsigned int i=0; i<n_rows; ++i)
- {
- Teuchos::RCP<const Thyra::VectorSpaceBase<double> > 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<Thyra::PhysicallyBlockedLinearOpBase<double> >
- tmpA = Teuchos::rcp(new Thyra::DefaultBlockedLinearOp<double>());
- tmpA->beginBlockFill(n_rows,n_cols);
- for (unsigned int i=0; i<n_rows; ++i)
- for (unsigned int j=0; j<n_cols; ++j)
- {
- Teuchos::RCP<const Thyra::LinearOpBase<double> >
- 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<const Thyra::LinearOpBase<double> > baseA = tmpA;
-
- Thyra::ConstLinearOperator<double> A = baseA;
-
- // Convert block vectors for
- // right hand side and
- // solution to Thyra objects.
- Thyra::Vector<double> rhs = A.range().createMember();
- Thyra::Vector<double> sol = A.domain().createMember();
- for (unsigned int i=0; i<n_rows; ++i)
- {
- Epetra_Vector *block_ptr = const_cast<Epetra_Vector*>
- ((input_b.block(i).trilinos_vector())(0));
- Teuchos::RCP<Thyra::VectorBase<double> > tmp_rhs_i
- = Thyra::create_Vector(Teuchos::rcp(block_ptr, false),
- epetra_vector_spaces[i]);
- const Thyra::Vector<double> rhs_i = tmp_rhs_i;
- rhs.setBlock(i, rhs_i);
-
- block_ptr = (input_x.block(i).trilinos_vector())(0);
- Teuchos::RCP<Thyra::VectorBase<double> > tmp_sol_i
- = Thyra::create_Vector(Teuchos::rcp(block_ptr, false),
- epetra_vector_spaces[i]);
- Thyra::Vector<double> 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<Teuchos::ParameterList> 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
- // <tt>solver_name</tt> 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<Thyra::LinearOpWithSolveFactoryBase<double> >
- aztecBlockLowsFactory = Teuchos::rcp
- (new Thyra::AztecOOLinearOpWithSolveFactory());
-
- aztecBlockLowsFactory->setParameterList(aztecBlockParams);
-
- Teuchos::RCP<Thyra::LinearOpWithSolveBase<double> > 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<const Thyra::LinearOpSourceBase<double> > A_ptr
- = Teuchos::rcp(new Thyra::DefaultLinearOpSource<double>(A.constPtr()),
- true);
- Teuchos::RCP<Thyra::PreconditionerBase<double> > P_ptr =
- Teuchos::rcp (new Thyra::DefaultPreconditioner<double>(), true);
- {
- Thyra::DefaultPreconditioner<double>
- *defaultPrec = &Teuchos::dyn_cast<Thyra::DefaultPreconditioner<double> >
- (*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<double>
- 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<int,unsigned int> 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
projects / dealii-svn.git / commitdiff