<h3>Specific improvements</h3>
<ol>
+ <li> New: RelaxationBlock classes for geometric multigrid now support parallel
+ computations using Trilinos.
+ <br>
+ (Timo Heister, Guido Kanschat, 2016/08/08)
+ </li>
+
<li> Improved: The regular and hp versions of
DoFTools::make_flux_sparsity_pattern() no longer use the user flags of the
underlying triangulation to determine if entries along a certain face have been
* implementation relies on a straight forward implementation of the Gauss-
* Seidel process.
*
+ * Parallel computations require you to specify an initialized
+ * ghost vector in AdditionalData::temp_ghost_vector.
+ *
* @ingroup Preconditioners
* @author Guido Kanschat
* @date 2010
*/
-template <typename MatrixType, typename inverse_type=typename MatrixType::value_type>
+template <typename MatrixType,
+ typename InverseNumberType = typename MatrixType::value_type,
+ typename VectorType = Vector<double> >
class RelaxationBlock :
- protected PreconditionBlockBase<inverse_type>
+ protected PreconditionBlockBase<InverseNumberType>
{
private:
/**
/**
* Value type for inverse matrices.
*/
- typedef inverse_type value_type;
+ typedef InverseNumberType value_type;
public:
/**
*/
AdditionalData (const double relaxation = 1.,
const bool invert_diagonal = true,
- const bool same_diagonal = false);
+ const bool same_diagonal = false,
+ const typename PreconditionBlockBase<InverseNumberType>::Inversion inversion
+ = PreconditionBlockBase<InverseNumberType>::gauss_jordan,
+ const double threshold = 0.,
+ VectorType *temp_ghost_vector = NULL);
/**
* The mapping from indices to blocks. Each row of this pattern enumerates
* particular if their sizes differ.
*/
bool same_diagonal;
+
/**
* Choose the inversion method for the blocks.
*/
- typename PreconditionBlockBase<inverse_type>::Inversion inversion;
+ typename PreconditionBlockBase<InverseNumberType>::Inversion inversion;
/**
* If #inversion is SVD, we can compute the Penrose-Moore inverse of the
* </ol>
*/
std::vector<std::vector<unsigned int> > order;
+
+ /**
+ * Temporary ghost vector that is used in the relaxation method when
+ * performing parallel MPI computations. The user is required to have this
+ * point to an initialized vector that contains all indices
+ * that appear in the @p block_list sa ghost values. Typically, this the
+ * set of locally active level DoFs. Unused when VectorType is a serial
+ * vector type like Vector<double>.
+ */
+ mutable VectorType *temp_ghost_vector;
+
/**
* Return the memory allocated in this object.
*/
* vectors). For the Jacobi step, the calling function must copy @p dst to
* @p pref after this.
*/
- template <class VECTOR>
void do_step (
- VECTOR &dst,
- const VECTOR &prev,
- const VECTOR &src,
+ VectorType &dst,
+ const VectorType &prev,
+ const VectorType &src,
const bool backward) const;
+
/**
* Pointer to the matrix. Make sure that the matrix exists as long as this
* class needs it, i.e. until calling @p invert_diagblocks, or (if the
* inverse matrices should not be stored) until the last call of the
* preconditioning @p vmult function of the derived classes.
*/
- SmartPointer<const MatrixType,RelaxationBlock<MatrixType,inverse_type> > A;
+ SmartPointer<const MatrixType,RelaxationBlock<MatrixType, InverseNumberType, VectorType> > A;
/**
* Control information.
*/
- SmartPointer<const AdditionalData, RelaxationBlock<MatrixType,inverse_type> > additional_data;
+ SmartPointer<const AdditionalData, RelaxationBlock<MatrixType, InverseNumberType, VectorType> > additional_data;
};
* @author Guido Kanschat
* @date 2010
*/
-template<typename MatrixType, typename inverse_type = typename MatrixType::value_type>
+template<typename MatrixType,
+ typename InverseNumberType = typename MatrixType::value_type,
+ typename VectorType = Vector<double> >
class RelaxationBlockJacobi : public virtual Subscriptor,
- protected RelaxationBlock<MatrixType, inverse_type>
+ protected RelaxationBlock<MatrixType, InverseNumberType, VectorType>
{
public:
/**
/**
* Make type publicly available.
*/
- using typename RelaxationBlock<MatrixType,inverse_type>::AdditionalData;
+ using typename RelaxationBlock<MatrixType, InverseNumberType, VectorType>::AdditionalData;
/**
* Make initialization function publicly available.
*/
- using RelaxationBlock<MatrixType, inverse_type>::initialize;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::initialize;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::clear;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::clear;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::empty;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::empty;
+ /**
+ * Make function of base class public again.
+ */
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::size;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::size;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::inverse;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::inverse;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::inverse_householder;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::inverse_householder;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::inverse_svd;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::inverse_svd;
- using PreconditionBlockBase<inverse_type>::log_statistics;
+ using PreconditionBlockBase<InverseNumberType>::log_statistics;
/**
* Perform one step of the Jacobi iteration.
*/
- template <class VECTOR>
- void step (VECTOR &dst, const VECTOR &rhs) const;
+ void step (VectorType &dst, const VectorType &rhs) const;
/**
* Perform one step of the Jacobi iteration.
*/
- template <class VECTOR>
- void Tstep (VECTOR &dst, const VECTOR &rhs) const;
+ void Tstep (VectorType &dst, const VectorType &rhs) const;
/**
* Return the memory allocated in this object.
* @author Guido Kanschat
* @date 2010
*/
-template<typename MatrixType, typename inverse_type = typename MatrixType::value_type>
+template<typename MatrixType,
+ typename InverseNumberType = typename MatrixType::value_type,
+ typename VectorType = Vector<double> >
class RelaxationBlockSOR : public virtual Subscriptor,
- protected RelaxationBlock<MatrixType, inverse_type>
+ protected RelaxationBlock<MatrixType, InverseNumberType, VectorType>
{
public:
/**
/**
* Make type publicly available.
*/
- using typename RelaxationBlock<MatrixType,inverse_type>::AdditionalData;
+ using typename RelaxationBlock<MatrixType, InverseNumberType, VectorType>::AdditionalData;
/**
* Make initialization function publicly available.
*/
- using RelaxationBlock<MatrixType, inverse_type>::initialize;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::initialize;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::clear;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::clear;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::empty;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::empty;
+ /**
+ * Make function of base class public again.
+ */
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::size;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::size;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::inverse;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::inverse;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::inverse_householder;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::inverse_householder;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::inverse_svd;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::inverse_svd;
- using PreconditionBlockBase<inverse_type>::log_statistics;
+ using PreconditionBlockBase<InverseNumberType>::log_statistics;
/**
* Perform one step of the SOR iteration.
*/
- template <class VECTOR>
- void step (VECTOR &dst, const VECTOR &rhs) const;
+ void step (VectorType &dst, const VectorType &rhs) const;
/**
* Perform one step of the transposed SOR iteration.
*/
- template <class VECTOR>
- void Tstep (VECTOR &dst, const VECTOR &rhs) const;
+ void Tstep (VectorType &dst, const VectorType &rhs) const;
};
* @author Guido Kanschat
* @date 2010
*/
-template<typename MatrixType, typename inverse_type = typename MatrixType::value_type>
+template<typename MatrixType,
+ typename InverseNumberType = typename MatrixType::value_type,
+ typename VectorType = Vector<double> >
class RelaxationBlockSSOR : public virtual Subscriptor,
- protected RelaxationBlock<MatrixType, inverse_type>
+ protected RelaxationBlock<MatrixType, InverseNumberType, VectorType>
{
public:
/**
/**
* Make type publicly available.
*/
- using typename RelaxationBlock<MatrixType,inverse_type>::AdditionalData;
+ using typename RelaxationBlock<MatrixType, InverseNumberType, VectorType>::AdditionalData;
/**
* Make initialization function publicly available.
*/
- using RelaxationBlock<MatrixType, inverse_type>::initialize;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::initialize;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::clear;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::clear;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::empty;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::empty;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::size;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::size;
+ /**
+ * Make function of base class public again.
+ */
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::inverse;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::inverse;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::inverse_householder;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::inverse_householder;
+ using RelaxationBlock<MatrixType, InverseNumberType, VectorType>::inverse_svd;
/**
* Make function of base class public again.
*/
- using RelaxationBlock<MatrixType, inverse_type>::inverse_svd;
- using PreconditionBlockBase<inverse_type>::log_statistics;
+ using PreconditionBlockBase<InverseNumberType>::log_statistics;
/**
* Perform one step of the SSOR iteration.
*/
- template <class VECTOR>
- void step (VECTOR &dst, const VECTOR &rhs) const;
+ void step (VectorType &dst, const VectorType &rhs) const;
/**
* Perform one step of the transposed SSOR iteration.
*/
- template <class VECTOR>
- void Tstep (VECTOR &dst, const VECTOR &rhs) const;
+ void Tstep (VectorType &dst, const VectorType &rhs) const;
};
#include <deal.II/lac/relaxation_block.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/vector_memory.h>
+#include <deal.II/lac/trilinos_vector.h>
DEAL_II_NAMESPACE_OPEN
-template <typename MatrixType, typename inverse_type>
+template <typename MatrixType, typename InverseNumberType, typename VectorType>
inline
-RelaxationBlock<MatrixType,inverse_type>::AdditionalData::AdditionalData
+RelaxationBlock<MatrixType, InverseNumberType, VectorType>::AdditionalData::AdditionalData
(const double relaxation,
const bool invert_diagonal,
- const bool same_diagonal)
+ const bool same_diagonal,
+ const typename PreconditionBlockBase<InverseNumberType>::Inversion inversion,
+ const double threshold,
+ VectorType *temp_ghost_vector)
:
relaxation(relaxation),
invert_diagonal(invert_diagonal),
same_diagonal(same_diagonal),
- inversion(PreconditionBlockBase<inverse_type>::gauss_jordan),
- threshold(0.)
+ inversion(inversion),
+ threshold(threshold),
+ temp_ghost_vector (temp_ghost_vector)
{}
-template <typename MatrixType, typename inverse_type>
+template <typename MatrixType, typename InverseNumberType, typename VectorType>
inline
std::size_t
-RelaxationBlock<MatrixType,inverse_type>::AdditionalData::memory_consumption() const
+RelaxationBlock<MatrixType, InverseNumberType, VectorType>::AdditionalData::memory_consumption() const
{
std::size_t result = sizeof(*this)
+ block_list.memory_consumption()
}
-template <typename MatrixType, typename inverse_type>
+template <typename MatrixType, typename InverseNumberType, typename VectorType>
inline
void
-RelaxationBlock<MatrixType,inverse_type>::initialize (const MatrixType &M,
- const AdditionalData ¶meters)
+RelaxationBlock<MatrixType, InverseNumberType, VectorType>::initialize (const MatrixType &M,
+ const AdditionalData ¶meters)
{
Assert (parameters.invert_diagonal, ExcNotImplemented());
}
-template <typename MatrixType, typename inverse_type>
+template <typename MatrixType, typename InverseNumberType, typename VectorType>
inline
void
-RelaxationBlock<MatrixType,inverse_type>::clear ()
+RelaxationBlock<MatrixType, InverseNumberType, VectorType>::clear ()
{
A = 0;
additional_data = 0;
- PreconditionBlockBase<inverse_type>::clear ();
+ PreconditionBlockBase<InverseNumberType>::clear ();
}
-template <typename MatrixType, typename inverse_type>
+template <typename MatrixType, typename InverseNumberType, typename VectorType>
inline
void
-RelaxationBlock<MatrixType,inverse_type>::invert_diagblocks ()
+RelaxationBlock<MatrixType, InverseNumberType, VectorType>::invert_diagblocks ()
{
const MatrixType &M=*A;
- FullMatrix<inverse_type> M_cell;
+ FullMatrix<InverseNumberType> M_cell;
if (this->same_diagonal())
{
}
switch (this->inversion)
{
- case PreconditionBlockBase<inverse_type>::gauss_jordan:
+ case PreconditionBlockBase<InverseNumberType>::gauss_jordan:
this->inverse(block).reinit(bs, bs);
this->inverse(block).invert(M_cell);
break;
- case PreconditionBlockBase<inverse_type>::householder:
+ case PreconditionBlockBase<InverseNumberType>::householder:
this->inverse_householder(block).initialize(M_cell);
break;
- case PreconditionBlockBase<inverse_type>::svd:
+ case PreconditionBlockBase<InverseNumberType>::svd:
this->inverse_svd(block).reinit(bs, bs);
this->inverse_svd(block) = M_cell;
this->inverse_svd(block).compute_inverse_svd(additional_data->threshold);
this->inverses_computed(true);
}
-
-template <typename MatrixType, typename inverse_type>
-template <class VECTOR>
+namespace internal
+{
+ /**
+ * Default implementation for serial vectors. Here we don't need to make a
+ * copy into a ghosted vector, so just return a reference to @p prev.
+ */
+ template <class VectorType>
+ const VectorType &
+ prepare_ghost_vector(
+ const VectorType &prev,
+ VectorType *other)
+ {
+ // If the following Assertion triggers, you either set temp_ghost_vector
+ // for a serial computation (don't!), or nobody implemented, instantiated, and
+ // tested the parallel version for your vector type.
+ Assert(other==NULL, ExcNotImplemented());
+ (void)other;
+ return prev;
+ }
+
+ /**
+ * Specialization for Trilinos. Use the ghosted vector.
+ */
+ template <>
+ const TrilinosWrappers::MPI::Vector &
+ prepare_ghost_vector(
+ const TrilinosWrappers::MPI::Vector &prev,
+ TrilinosWrappers::MPI::Vector *other)
+ {
+ Assert(other!=NULL,
+ ExcMessage("You need to provide a ghosted vector in RelaxationBlock::AdditionalData::temp_trilinos_ghost_vector."));
+ Assert(other->size()==prev.size(), ExcInternalError());
+
+ // import ghost values:
+ *other = prev;
+ return *other;
+ }
+} // end namespace internal
+
+template <typename MatrixType, typename InverseNumberType, typename VectorType>
inline
void
-RelaxationBlock<MatrixType,inverse_type>::do_step (VECTOR &dst,
- const VECTOR &prev,
- const VECTOR &src,
- const bool backward) const
+RelaxationBlock<MatrixType, InverseNumberType, VectorType>::do_step (VectorType &dst,
+ const VectorType &prev,
+ const VectorType &src,
+ const bool backward) const
{
Assert (additional_data->invert_diagonal, ExcNotImplemented());
+ const VectorType &ghosted_prev = internal::prepare_ghost_vector(prev, additional_data->temp_ghost_vector);
+
const MatrixType &M=*this->A;
- Vector<typename VECTOR::value_type> b_cell, x_cell;
+ Vector<typename VectorType::value_type> b_cell, x_cell;
const bool permutation_empty = additional_data->order.size() == 0;
const unsigned int n_permutations = (permutation_empty)
b_cell(row_cell) = src(row->column());
for (typename MatrixType::const_iterator entry = M.begin(row->column());
entry != M.end(row->column()); ++entry)
- b_cell(row_cell) -= entry->value() * prev(entry->column());
+ b_cell(row_cell) -= entry->value() * ghosted_prev(entry->column());
}
// Apply inverse diagonal
this->inverse_vmult(block, x_cell, b_cell);
}
#endif
// Store in result vector
- row=additional_data->block_list.begin(block);
+ row = additional_data->block_list.begin(block);
for (size_type row_cell=0; row_cell<bs; ++row_cell, ++row)
dst(row->column()) += additional_data->relaxation * x_cell(row_cell);
}
//----------------------------------------------------------------------//
-template <typename MatrixType, typename inverse_type>
-template <class VECTOR>
-void RelaxationBlockJacobi<MatrixType,inverse_type>::step
-(VECTOR &dst,
- const VECTOR &src) const
+template <typename MatrixType, typename InverseNumberType, typename VectorType>
+void RelaxationBlockJacobi<MatrixType, InverseNumberType, VectorType>::step
+(VectorType &dst,
+ const VectorType &src) const
{
- GrowingVectorMemory<VECTOR > mem;
- typename VectorMemory<VECTOR >::Pointer aux = mem;
+ GrowingVectorMemory<VectorType> mem;
+ typename VectorMemory<VectorType>::Pointer aux = mem;
aux->reinit(dst, false);
*aux = dst;
this->do_step(dst, *aux, src, false);
}
-template <typename MatrixType, typename inverse_type>
-template <class VECTOR>
-void RelaxationBlockJacobi<MatrixType,inverse_type>::Tstep
-(VECTOR &dst,
- const VECTOR &src) const
+template <typename MatrixType, typename InverseNumberType, typename VectorType>
+void RelaxationBlockJacobi<MatrixType, InverseNumberType, VectorType>::Tstep
+(VectorType &dst,
+ const VectorType &src) const
{
- GrowingVectorMemory<VECTOR > mem;
- typename VectorMemory<VECTOR >::Pointer aux = mem;
+ GrowingVectorMemory<VectorType> mem;
+ typename VectorMemory<VectorType>::Pointer aux = mem;
aux->reinit(dst, false);
*aux = dst;
this->do_step(dst, *aux, src, true);
//----------------------------------------------------------------------//
-template <typename MatrixType, typename inverse_type>
-template <class VECTOR>
-void RelaxationBlockSOR<MatrixType,inverse_type>::step
-(VECTOR &dst,
- const VECTOR &src) const
+template <typename MatrixType, typename InverseNumberType, typename VectorType>
+void RelaxationBlockSOR<MatrixType, InverseNumberType, VectorType>::step
+(VectorType &dst,
+ const VectorType &src) const
{
this->do_step(dst, dst, src, false);
}
-template <typename MatrixType, typename inverse_type>
-template <class VECTOR>
-void RelaxationBlockSOR<MatrixType,inverse_type>::Tstep
-(VECTOR &dst,
- const VECTOR &src) const
+template <typename MatrixType, typename InverseNumberType, typename VectorType>
+void RelaxationBlockSOR<MatrixType, InverseNumberType, VectorType>::Tstep
+(VectorType &dst,
+ const VectorType &src) const
{
this->do_step(dst, dst, src, true);
}
//----------------------------------------------------------------------//
-template <typename MatrixType, typename inverse_type>
-template <class VECTOR>
-void RelaxationBlockSSOR<MatrixType,inverse_type>::step
-(VECTOR &dst,
- const VECTOR &src) const
+template <typename MatrixType, typename InverseNumberType, typename VectorType>
+void RelaxationBlockSSOR<MatrixType, InverseNumberType, VectorType>::step
+(VectorType &dst,
+ const VectorType &src) const
{
this->do_step(dst, dst, src, false);
this->do_step(dst, dst, src, true);
}
-template <typename MatrixType, typename inverse_type>
-template <class VECTOR>
-void RelaxationBlockSSOR<MatrixType,inverse_type>::Tstep
-(VECTOR &dst,
- const VECTOR &src) const
+template <typename MatrixType, typename InverseNumberType, typename VectorType>
+void RelaxationBlockSSOR<MatrixType, InverseNumberType, VectorType>::Tstep
+(VectorType &dst,
+ const VectorType &src) const
{
this->do_step(dst, dst, src, true);
this->do_step(dst, dst, src, false);
-for (S1, S2 : REAL_SCALARS)
+for (S1, S2, S3 : REAL_SCALARS)
{
- template class RelaxationBlock<SparseMatrix<S1>, S2>;
- template class RelaxationBlockJacobi<SparseMatrix<S1>, S2>;
- template class RelaxationBlockSOR<SparseMatrix<S1>, S2>;
- template class RelaxationBlockSSOR<SparseMatrix<S1>, S2>;
+ template class RelaxationBlock<SparseMatrix<S1>, S2, Vector<S3> >;
+ template class RelaxationBlockJacobi<SparseMatrix<S1>, S2, Vector<S3> >;
+ template class RelaxationBlockSOR<SparseMatrix<S1>, S2, Vector<S3> >;
+ template class RelaxationBlockSSOR<SparseMatrix<S1>, S2, Vector<S3> >;
}
-for (S1, S2, S3: REAL_SCALARS)
-{
-// ------------ RelaxationBlockJacobi -----------------
- template
- void RelaxationBlockJacobi<SparseMatrix<S1>, S2>::step<S3>
- (Vector<S3> &, const Vector<S3> &) const;
- template
- void RelaxationBlockJacobi<SparseMatrix<S1>, S2>::Tstep<S3>
- (Vector<S3> &, const Vector<S3> &) const;
-
-// ------------ RelaxationBlockSOR -----------------
- template
- void RelaxationBlockSOR<SparseMatrix<S1>, S2>::step<S3>
- (Vector<S3> &, const Vector<S3> &) const;
- template
- void RelaxationBlockSOR<SparseMatrix<S1>, S2>::Tstep<S3>
- (Vector<S3> &, const Vector<S3> &) const;
-
-// ------------ RelaxationBlockSSOR -----------------
- template
- void RelaxationBlockSSOR<SparseMatrix<S1>, S2>::step<S3>
- (Vector<S3> &, const Vector<S3> &) const;
- template
- void RelaxationBlockSSOR<SparseMatrix<S1>, S2>::Tstep<S3>
- (Vector<S3> &, const Vector<S3> &) const;
-}
-
-
-for (S2 : REAL_SCALARS)
+for (S1 : REAL_SCALARS)
{
#ifdef DEAL_II_WITH_TRILINOS
- template class RelaxationBlock<TrilinosWrappers::SparseMatrix, S2>;
- template class RelaxationBlockJacobi<TrilinosWrappers::SparseMatrix, S2>;
- template class RelaxationBlockSOR<TrilinosWrappers::SparseMatrix, S2>;
- template class RelaxationBlockSSOR<TrilinosWrappers::SparseMatrix, S2>;
-
- template void RelaxationBlockJacobi<TrilinosWrappers::SparseMatrix, S2>::step
- (TrilinosWrappers::MPI::Vector &, const TrilinosWrappers::MPI::Vector &) const;
- template void RelaxationBlockJacobi<TrilinosWrappers::SparseMatrix, S2>::Tstep
- (TrilinosWrappers::MPI::Vector &, const TrilinosWrappers::MPI::Vector &) const;
-
- template void RelaxationBlockSOR<TrilinosWrappers::SparseMatrix, S2>::step
- (TrilinosWrappers::MPI::Vector &, const TrilinosWrappers::MPI::Vector &) const;
- template void RelaxationBlockSOR<TrilinosWrappers::SparseMatrix, S2>::Tstep
- (TrilinosWrappers::MPI::Vector &, const TrilinosWrappers::MPI::Vector &) const;
-
-// ------------ RelaxationBlockSSOR -----------------
- template void RelaxationBlockSSOR<TrilinosWrappers::SparseMatrix, S2>::step
- (TrilinosWrappers::MPI::Vector &, const TrilinosWrappers::MPI::Vector &) const;
- template void RelaxationBlockSSOR<TrilinosWrappers::SparseMatrix, S2>::Tstep
- (TrilinosWrappers::MPI::Vector &, const TrilinosWrappers::MPI::Vector &) const;
+ template class RelaxationBlock<TrilinosWrappers::SparseMatrix, S1, TrilinosWrappers::MPI::Vector>;
+ template class RelaxationBlockJacobi<TrilinosWrappers::SparseMatrix, S1, TrilinosWrappers::MPI::Vector>;
+ template class RelaxationBlockSOR<TrilinosWrappers::SparseMatrix, S1, TrilinosWrappers::MPI::Vector>;
+ template class RelaxationBlockSSOR<TrilinosWrappers::SparseMatrix, S1, TrilinosWrappers::MPI::Vector>;
#endif
}
// ---------------------------------------------------------------------
/*
- * Author: Guido Kanschat, Texas A&M University, 2009
+ * Test RelaxationBlockJacobi in parallel
*/
#include "../tests.h"
MGCoarseGridLACIteration<SolverCG<TrilinosWrappers::MPI::Vector>,TrilinosWrappers::MPI::Vector>
coarse_grid_solver(coarse_solver, coarse_matrix, identity);
- typedef RelaxationBlockJacobi<TrilinosWrappers::SparseMatrix, double>
+ typedef RelaxationBlockJacobi<TrilinosWrappers::SparseMatrix, double, TrilinosWrappers::MPI::Vector>
Smoother;
MGLevelObject<typename Smoother::AdditionalData> smoother_data;
smoother_data.resize(0, triangulation.n_levels() - 1);
mg::SmootherRelaxation<Smoother, TrilinosWrappers::MPI::Vector> mg_smoother;
+ MGLevelObject<TrilinosWrappers::MPI::Vector> temp_vectors(0, triangulation.n_levels() - 1);
+
+
for (unsigned int l = smoother_data.min_level() + 1; l <= smoother_data.max_level(); ++l)
{
DoFTools::make_cell_patches(smoother_data[l].block_list, this->dof_handler, l);
- smoother_data[l].block_list.compress();
+ if (smoother_data[l].block_list.n_rows()>0)
+ smoother_data[l].block_list.compress();
smoother_data[l].relaxation = 0.7;
smoother_data[l].inversion = PreconditionBlockBase<double>::svd;
+ TrilinosWrappers::MPI::Vector *ghost = &(temp_vectors[l]);
+ IndexSet relevant_dofs;
+ DoFTools::extract_locally_relevant_level_dofs(dof_handler, l, relevant_dofs);
+ ghost->reinit(dof_handler.locally_owned_mg_dofs(l),
+ relevant_dofs,
+ MPI_COMM_WORLD);
+ smoother_data[l].temp_ghost_vector = ghost;
}
mg_smoother.initialize(mg_matrix, smoother_data);
- mg_smoother.set_steps(1);
- mg_smoother.set_variable(true);
+ mg_smoother.set_steps(2);
+ mg_smoother.set_variable(false);
mg::Matrix<TrilinosWrappers::MPI::Vector> mgmatrix(mg_matrix);
mg::Matrix<TrilinosWrappers::MPI::Vector> mgdown(mg_matrix_dg_down);
Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads());
MPILogInitAll log;
- try
- {
- FE_DGQ<2> fe1(2);
- InteriorPenaltyProblem<2> test1(fe1);
- test1.run(6);
- }
- catch (std::exception &exc)
- {
- std::cerr << std::endl << std::endl
- << "----------------------------------------------------"
- << std::endl;
- std::cerr << "Exception on processing: " << std::endl
- << exc.what() << std::endl
- << "Aborting!" << std::endl
- << "----------------------------------------------------"
- << std::endl;
- return 1;
- }
- catch (...)
- {
- std::cerr << std::endl << std::endl
- << "----------------------------------------------------"
- << std::endl;
- std::cerr << "Unknown exception!" << std::endl
- << "Aborting!" << std::endl
- << "----------------------------------------------------"
- << std::endl;
- return 1;
- }
-
+ FE_DGQ<2> fe1(2);
+ InteriorPenaltyProblem<2> test1(fe1);
+ test1.run(6);
return 0;
}
-
\ No newline at end of file
+
+DEAL:0::Element: FE_DGQ<2>(2)
+DEAL:0::Step 0
+DEAL:0::Triangulation 12 cells, 2 levels
+DEAL:0::DoFHandler 108 dofs, level dofs 27 108
+DEAL:0::Assemble matrix
+DEAL:0::Assemble multilevel matrix
+DEAL:0::Assemble right hand side
+DEAL:0::Solve
+DEAL:0:cg::Starting value 20.6576
+DEAL:0:cg::Convergence step 15 value 4.82452e-13
+DEAL:0::Estimate 1.10045
+DEAL:0::Step 1
+DEAL:0::Triangulation 18 cells, 3 levels
+DEAL:0::DoFHandler 162 dofs, level dofs 27 108 72
+DEAL:0::Assemble matrix
+DEAL:0::Assemble multilevel matrix
+DEAL:0::Assemble right hand side
+DEAL:0::Solve
+DEAL:0:cg::Starting value 21.1045
+DEAL:0:cg::Convergence step 15 value 9.35135e-13
+DEAL:0::Estimate 1.03046
+DEAL:0::Step 2
+DEAL:0::Triangulation 27 cells, 4 levels
+DEAL:0::DoFHandler 243 dofs, level dofs 27 108 108 72
+DEAL:0::Assemble matrix
+DEAL:0::Assemble multilevel matrix
+DEAL:0::Assemble right hand side
+DEAL:0::Solve
+DEAL:0:cg::Starting value 21.3244
+DEAL:0:cg::Convergence step 16 value 7.03391e-13
+DEAL:0::Estimate 0.836764
+DEAL:0::Step 3
+DEAL:0::Triangulation 39 cells, 5 levels
+DEAL:0::DoFHandler 351 dofs, level dofs 27 108 144 108 72
+DEAL:0::Assemble matrix
+DEAL:0::Assemble multilevel matrix
+DEAL:0::Assemble right hand side
+DEAL:0::Solve
+DEAL:0:cg::Starting value 21.4584
+DEAL:0:cg::Convergence step 17 value 2.88305e-13
+DEAL:0::Estimate 0.644552
+DEAL:0::Step 4
+DEAL:0::Triangulation 54 cells, 6 levels
+DEAL:0::DoFHandler 486 dofs, level dofs 27 108 180 144 108 72
+DEAL:0::Assemble matrix
+DEAL:0::Assemble multilevel matrix
+DEAL:0::Assemble right hand side
+DEAL:0::Solve
+DEAL:0:cg::Starting value 21.5437
+DEAL:0:cg::Convergence step 16 value 6.40351e-13
+DEAL:0::Estimate 0.474987
+DEAL:0::Step 5
+DEAL:0::Triangulation 84 cells, 7 levels
+DEAL:0::DoFHandler 756 dofs, level dofs 27 108 324 216 144 108 72
+DEAL:0::Assemble matrix
+DEAL:0::Assemble multilevel matrix
+DEAL:0::Assemble right hand side
+DEAL:0::Solve
+DEAL:0:cg::Starting value 25.5435
+DEAL:0:cg::Convergence step 16 value 4.60775e-13
+DEAL:0::Estimate 0.356044
+
+DEAL:1::Element: FE_DGQ<2>(2)
+DEAL:1::Step 0
+DEAL:1::Triangulation 12 cells, 2 levels
+DEAL:1::DoFHandler 108 dofs, level dofs 27 108
+DEAL:1::Assemble matrix
+DEAL:1::Assemble multilevel matrix
+DEAL:1::Assemble right hand side
+DEAL:1::Solve
+DEAL:1:cg::Starting value 20.6576
+DEAL:1:cg::Convergence step 15 value 4.82452e-13
+DEAL:1::Estimate 1.10045
+DEAL:1::Step 1
+DEAL:1::Triangulation 18 cells, 3 levels
+DEAL:1::DoFHandler 162 dofs, level dofs 27 108 72
+DEAL:1::Assemble matrix
+DEAL:1::Assemble multilevel matrix
+DEAL:1::Assemble right hand side
+DEAL:1::Solve
+DEAL:1:cg::Starting value 21.1045
+DEAL:1:cg::Convergence step 15 value 9.35135e-13
+DEAL:1::Estimate 1.03046
+DEAL:1::Step 2
+DEAL:1::Triangulation 27 cells, 4 levels
+DEAL:1::DoFHandler 243 dofs, level dofs 27 108 108 72
+DEAL:1::Assemble matrix
+DEAL:1::Assemble multilevel matrix
+DEAL:1::Assemble right hand side
+DEAL:1::Solve
+DEAL:1:cg::Starting value 21.3244
+DEAL:1:cg::Convergence step 16 value 7.03391e-13
+DEAL:1::Estimate 0.836764
+DEAL:1::Step 3
+DEAL:1::Triangulation 39 cells, 5 levels
+DEAL:1::DoFHandler 351 dofs, level dofs 27 108 144 108 72
+DEAL:1::Assemble matrix
+DEAL:1::Assemble multilevel matrix
+DEAL:1::Assemble right hand side
+DEAL:1::Solve
+DEAL:1:cg::Starting value 21.4584
+DEAL:1:cg::Convergence step 17 value 2.88305e-13
+DEAL:1::Estimate 0.644552
+DEAL:1::Step 4
+DEAL:1::Triangulation 54 cells, 6 levels
+DEAL:1::DoFHandler 486 dofs, level dofs 27 108 180 144 108 72
+DEAL:1::Assemble matrix
+DEAL:1::Assemble multilevel matrix
+DEAL:1::Assemble right hand side
+DEAL:1::Solve
+DEAL:1:cg::Starting value 21.5437
+DEAL:1:cg::Convergence step 16 value 6.40351e-13
+DEAL:1::Estimate 0.474987
+DEAL:1::Step 5
+DEAL:1::Triangulation 84 cells, 7 levels
+DEAL:1::DoFHandler 756 dofs, level dofs 27 108 324 216 144 108 72
+DEAL:1::Assemble matrix
+DEAL:1::Assemble multilevel matrix
+DEAL:1::Assemble right hand side
+DEAL:1::Solve
+DEAL:1:cg::Starting value 25.5435
+DEAL:1:cg::Convergence step 16 value 4.60775e-13
+DEAL:1::Estimate 0.356044
+
+
+DEAL:2::Element: FE_DGQ<2>(2)
+DEAL:2::Step 0
+DEAL:2::Triangulation 12 cells, 2 levels
+DEAL:2::DoFHandler 108 dofs, level dofs 27 108
+DEAL:2::Assemble matrix
+DEAL:2::Assemble multilevel matrix
+DEAL:2::Assemble right hand side
+DEAL:2::Solve
+DEAL:2:cg::Starting value 20.6576
+DEAL:2:cg::Convergence step 15 value 4.82452e-13
+DEAL:2::Estimate 1.10045
+DEAL:2::Step 1
+DEAL:2::Triangulation 18 cells, 3 levels
+DEAL:2::DoFHandler 162 dofs, level dofs 27 108 72
+DEAL:2::Assemble matrix
+DEAL:2::Assemble multilevel matrix
+DEAL:2::Assemble right hand side
+DEAL:2::Solve
+DEAL:2:cg::Starting value 21.1045
+DEAL:2:cg::Convergence step 15 value 9.35135e-13
+DEAL:2::Estimate 1.03046
+DEAL:2::Step 2
+DEAL:2::Triangulation 27 cells, 4 levels
+DEAL:2::DoFHandler 243 dofs, level dofs 27 108 108 72
+DEAL:2::Assemble matrix
+DEAL:2::Assemble multilevel matrix
+DEAL:2::Assemble right hand side
+DEAL:2::Solve
+DEAL:2:cg::Starting value 21.3244
+DEAL:2:cg::Convergence step 16 value 7.03391e-13
+DEAL:2::Estimate 0.836764
+DEAL:2::Step 3
+DEAL:2::Triangulation 39 cells, 5 levels
+DEAL:2::DoFHandler 351 dofs, level dofs 27 108 144 108 72
+DEAL:2::Assemble matrix
+DEAL:2::Assemble multilevel matrix
+DEAL:2::Assemble right hand side
+DEAL:2::Solve
+DEAL:2:cg::Starting value 21.4584
+DEAL:2:cg::Convergence step 17 value 2.88305e-13
+DEAL:2::Estimate 0.644552
+DEAL:2::Step 4
+DEAL:2::Triangulation 54 cells, 6 levels
+DEAL:2::DoFHandler 486 dofs, level dofs 27 108 180 144 108 72
+DEAL:2::Assemble matrix
+DEAL:2::Assemble multilevel matrix
+DEAL:2::Assemble right hand side
+DEAL:2::Solve
+DEAL:2:cg::Starting value 21.5437
+DEAL:2:cg::Convergence step 16 value 6.40351e-13
+DEAL:2::Estimate 0.474987
+DEAL:2::Step 5
+DEAL:2::Triangulation 84 cells, 7 levels
+DEAL:2::DoFHandler 756 dofs, level dofs 27 108 324 216 144 108 72
+DEAL:2::Assemble matrix
+DEAL:2::Assemble multilevel matrix
+DEAL:2::Assemble right hand side
+DEAL:2::Solve
+DEAL:2:cg::Starting value 25.5435
+DEAL:2:cg::Convergence step 16 value 4.60775e-13
+DEAL:2::Estimate 0.356044
+