]> https://gitweb.dealii.org/ - dealii.git/commitdiff
implement parallel RelaxationBlock 2908/head
authorTimo Heister <timo.heister@gmail.com>
Tue, 2 Aug 2016 19:04:01 +0000 (21:04 +0200)
committerTimo Heister <timo.heister@gmail.com>
Mon, 8 Aug 2016 20:50:38 +0000 (16:50 -0400)
- use a ghosted vector in RelaxationBlock
- introduce VectorType in base class RelaxationBlock
- update test for parallel RelaxationBlock
- add documentation

doc/news/changes.h
include/deal.II/lac/relaxation_block.h
include/deal.II/lac/relaxation_block.templates.h
source/lac/relaxation_block.inst.in
tests/mpi/step-39-block.cc
tests/mpi/step-39-block.with_trilinos=true.mpirun=3.output

index 2e59d52a6bc0661df807cde51cddc2af5b9ea8a4..6444cdfac0a1114d3ad2849e7eda2593f24cb5c9 100644 (file)
@@ -428,6 +428,12 @@ SparsityPattern.
 <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
index ce5597255d3372133483187be8efebcb5c409f54..696fea920631223016cc22af2fd146f7e8e38c7a 100644 (file)
@@ -43,13 +43,18 @@ DEAL_II_NAMESPACE_OPEN
  * 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:
   /**
@@ -60,7 +65,7 @@ private:
   /**
    * Value type for inverse matrices.
    */
-  typedef inverse_type value_type;
+  typedef InverseNumberType value_type;
 
 public:
   /**
@@ -82,7 +87,11 @@ 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
@@ -112,10 +121,11 @@ public:
      * 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
@@ -147,6 +157,17 @@ public:
      * </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.
      */
@@ -210,24 +231,24 @@ protected:
    * 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;
 };
 
 
@@ -246,9 +267,11 @@ protected:
  * @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:
   /**
@@ -264,50 +287,51 @@ 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.
@@ -331,9 +355,11 @@ public:
  * @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:
   /**
@@ -349,50 +375,51 @@ 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;
 };
 
 
@@ -411,9 +438,11 @@ public:
  * @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:
   /**
@@ -424,51 +453,52 @@ 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;
 };
 
 
index 70a99df7a706cabbad357a59bd75b419bc1f80d1..84eafb84ae69409d646da21ab7c3e7f2402b4e7c 100644 (file)
 #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()
@@ -51,11 +56,11 @@ RelaxationBlock<MatrixType,inverse_type>::AdditionalData::memory_consumption() c
 }
 
 
-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 &parameters)
+RelaxationBlock<MatrixType, InverseNumberType, VectorType>::initialize (const MatrixType     &M,
+    const AdditionalData &parameters)
 {
   Assert (parameters.invert_diagonal, ExcNotImplemented());
 
@@ -73,24 +78,24 @@ RelaxationBlock<MatrixType,inverse_type>::initialize (const MatrixType     &M,
 }
 
 
-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())
     {
@@ -131,14 +136,14 @@ RelaxationBlock<MatrixType,inverse_type>::invert_diagblocks ()
             }
           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);
@@ -151,20 +156,59 @@ RelaxationBlock<MatrixType,inverse_type>::invert_diagblocks ()
   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)
@@ -197,7 +241,7 @@ RelaxationBlock<MatrixType,inverse_type>::do_step (VECTOR       &dst,
               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);
@@ -208,7 +252,7 @@ RelaxationBlock<MatrixType,inverse_type>::do_step (VECTOR       &dst,
             }
 #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);
         }
@@ -219,28 +263,26 @@ RelaxationBlock<MatrixType,inverse_type>::do_step (VECTOR       &dst,
 
 //----------------------------------------------------------------------//
 
-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);
@@ -249,21 +291,19 @@ void RelaxationBlockJacobi<MatrixType,inverse_type>::Tstep
 
 //----------------------------------------------------------------------//
 
-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);
 }
@@ -271,22 +311,20 @@ void RelaxationBlockSOR<MatrixType,inverse_type>::Tstep
 
 //----------------------------------------------------------------------//
 
-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);
index 47cf6b86fbb424ff5c658e3504b13904793c293e..2cc7accd1833062b9f450a8b5d76674a970d3569 100644 (file)
 
 
 
-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
 }
 
index e24b17e9a3957ad49a522b5631917fd3f356081f..c82f3ae28b3ecce3e0533f07c6bcfad74a615051 100644 (file)
@@ -14,7 +14,7 @@
 // ---------------------------------------------------------------------
 
 /*
- * Author: Guido Kanschat, Texas A&M University, 2009
+ * Test RelaxationBlockJacobi in parallel
  */
 
 #include "../tests.h"
@@ -565,24 +565,35 @@ namespace Step39
     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);
@@ -712,35 +723,8 @@ int main(int argc, char *argv[])
   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;
 }
index 0519ecba6ea913e21689ec692e81e9e4973fbf73..75cd43ea5179eb2c8bc8c45c0f8b4a31ae9025eb 100644 (file)
@@ -1 +1,188 @@
\ 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
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.