]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
new classes BlockList, RelaxationBlockSOR and RelaxationBlockSSOR
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 19 Oct 2010 19:20:22 +0000 (19:20 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 19 Oct 2010 19:20:22 +0000 (19:20 +0000)
git-svn-id: https://svn.dealii.org/trunk@22387 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/lac/include/lac/block_list.h [new file with mode: 0644]
deal.II/lac/include/lac/precondition_block.h
deal.II/lac/include/lac/relaxation_block.h [new file with mode: 0644]
deal.II/lac/include/lac/relaxation_block.templates.h [new file with mode: 0644]
deal.II/lac/source/relaxation_block.cc [new file with mode: 0644]
deal.II/lac/source/relaxation_block.inst.in [new file with mode: 0644]

index 31c376625cc1f814bb99ded482194a9d073e74bd..f049643af7b4c29299f697144f8090fef1a28434 100644 (file)
@@ -252,7 +252,16 @@ through DoFHandler::get_tria() and DoFHandler::get_fe(), respectively.
 <h3>lac</h3>
 
 <ol>
-  <li><p>New: The ConstraintMatrix::merge function now takes a second
+
+  <li><p>New: The classes RelaxationBlockSOR and RelaxationBlockSSOR
+  implement overlapping Schwarz relaxation methods. Additionally,
+  their base class RelaxationBlock and the helper class BlockList have
+  been added to the library.
+  <br>
+  (GK 2010/10/19)
+  </p></li>
+  
+  <li><p>Improved: The ConstraintMatrix::merge function now takes a second
   argument that indicates what should happen if the two objects to be
   merged have constraints on the very same degree of freedom.
   <br>
diff --git a/deal.II/lac/include/lac/block_list.h b/deal.II/lac/include/lac/block_list.h
new file mode 100644 (file)
index 0000000..e0cc124
--- /dev/null
@@ -0,0 +1,338 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 2010 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__block_list_h
+#define __deal2__block_list_h
+
+
+#include <base/subscriptor.h>
+
+#include <vector>
+#include <utility>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * A vector of index sets listing the indices of small blocks of a
+ * linear system. For each block, the indices in that block are
+ * listed.
+ *
+ * The focus of this class is on small blocks of degrees of freedom
+ * associated with a single mesh cell or a small patch. These indices
+ * may be contiguous or not and we do not optimize for the case they
+ * are. For larger sets, the use of a vector of IndexSet objects might
+ * be advisable.
+ *
+ * BlockList objects can conveniently be initialized with iterator
+ * ranges from DoFHandler or MGDoFHandler, using either their cell or
+ * multigrid indices. Other initializations are possible to be
+ * implemented.
+ *
+ * @author Guido Kanschat
+ * @date 2010
+ */
+class BlockList :
+  public Subscriptor
+{
+  public:
+                                    /// The container for each index set
+    typedef std::vector<unsigned int> block_container;
+                                    /// The iterator for individual indices
+    typedef block_container::const_iterator const_iterator;
+    
+                                    /**
+                                     * Set up all index sets using an
+                                     * DoF iterator range. This
+                                     * function will call
+                                     * <tt>begin->get_dof_indices()</tt>
+                                     * with a signature like
+                                     * DoFCellAccessor::get_dof_indices().
+                                     * Typically, the iterators will
+                                     * loop over active cells of a
+                                     * triangulation.
+                                     *
+                                     * In addition, the function
+                                     * needs the total number of
+                                     * blocks as its first argument.
+                                     */
+    template <typename ITERATOR>
+    void initialize(unsigned int n_blocks,
+                   const ITERATOR begin,
+                   const typename identity<ITERATOR>::type end);
+    
+                                    /**
+                                     * Set up all index sets using an
+                                     * DoF iterator range. This
+                                     * function will call
+                                     * <tt>begin->get_mg_dof_indices()</tt>
+                                     * with a signature like
+                                     * MGDoFCellAccessor::get_mg_dof_indices().
+                                     * Typically, the iterators loop
+                                     * over the cells of a single
+                                     * level or a Triangulation.
+                                     *
+                                     * In addition, the function
+                                     * needs the total number of
+                                     * blocks as its first argument.
+                                     */
+    template <typename ITERATOR>
+    void initialize_mg(unsigned int n_blocks,
+                      const ITERATOR begin,
+                      const typename identity<ITERATOR>::type end);
+    
+                                    /**
+                                     * Set up all index sets using an
+                                     * DoF iterator range. This
+                                     * function will call
+                                     * <tt>begin->get_dof_indices()</tt>
+                                     * with a signature like
+                                     * DoFCellAccessor::get_dof_indices().
+                                     * Typically, the iterators will
+                                     * loop over active cells of a
+                                     * triangulation.
+                                     *
+                                     * The argument vector
+                                     * <tt>selected_dofs</tt> should
+                                     * have the length  of dofs per
+                                     * cell (thus, this function is
+                                     * not suitable for hp), and a
+                                     * true value for each degree of
+                                     * freedom which should be added
+                                     * to the index set of this
+                                     * cell. If you are working on a
+                                     * single block of a block
+                                     * system, the <tt>offset</tt> is
+                                     * the start index of this block.
+                                     *
+                                     * In addition, the function
+                                     * needs the total number of
+                                     * blocks as its first argument.
+                                     */
+    template <typename ITERATOR>
+    void initialize(unsigned int n_blocks,
+                   const ITERATOR begin,
+                   const typename identity<ITERATOR>::type end,
+                   const std::vector<bool>& selected_dofs,
+                   unsigned int offset = 0);
+                                    /**
+                                     * Set up all index sets using an
+                                     * DoF iterator range. This
+                                     * function will call
+                                     * <tt>begin->get_mg_dof_indices()</tt>
+                                     * with a signature like
+                                     * MGDoFCellAccessor::get_mg_dof_indices().
+                                     * Typically, the iterators will
+                                     * loop over cells on a single
+                                     * level of a triangulation.
+                                     *
+                                     * The argument vector
+                                     * <tt>selected_dofs</tt> should
+                                     * have the length  of dofs per
+                                     * cell (thus, this function is
+                                     * not suitable for hp), and a
+                                     * true value for each degree of
+                                     * freedom which should be added
+                                     * to the index set of this
+                                     * cell. If you are working on a
+                                     * single block of a block
+                                     * system, the <tt>offset</tt> is
+                                     * the start index of this block.
+                                     *
+                                     * In addition, the function
+                                     * needs the total number of
+                                     * blocks as its first argument.
+                                     */
+    template <typename ITERATOR>
+    void initialize_mg(unsigned int n_blocks,
+                      const ITERATOR begin,
+                      const typename identity<ITERATOR>::type end,
+                      const std::vector<bool>& selected_dofs,
+                      unsigned int offset = 0);
+
+                                    /**
+                                     * The number of blocks.
+                                     */
+    unsigned int size() const;
+    
+                                    /**
+                                     * The size of a single block.
+                                     */
+    unsigned int block_size(unsigned int block) const;
+
+                                    /**
+                                     * Iterator to the first index in block.
+                                     */
+    const_iterator begin(unsigned int block) const;
+                                    /**
+                                     * End iterator for a single block.
+                                     */
+    const_iterator end(unsigned int block) const;
+                                    /**
+                                     * Return the position of
+                                     * <tt>index</tt> in
+                                     * <tt>block</tt>, or
+                                     * numbers::invalid_unsigned_int,
+                                     * if the index is not in the block.
+                                     */
+    unsigned int local_index(unsigned int block, unsigned int index) const;
+    
+  private:
+                                    /**
+                                     * The container for t he index sets.
+                                     */
+    std::vector<block_container> index_sets;
+};
+
+
+template <typename ITERATOR>
+inline
+void
+BlockList::initialize(unsigned int n_blocks, const ITERATOR begin, const typename identity<ITERATOR>::type end)
+{
+  index_sets.resize(n_blocks);
+  std::vector<unsigned int> indices;
+  unsigned int k = 0;
+  for (ITERATOR cell = begin; cell != end; ++cell, ++k)
+    {
+      indices.resize(cell->get_fe().dofs_per_cell);
+      cell->get_dof_indices(indices);
+
+      for (std::vector<unsigned int>::const_iterator i=indices.begin();
+          i != indices.end(); ++i)
+       index_sets[k].push_back(*i);
+    }
+}
+
+
+template <typename ITERATOR>
+inline
+void
+BlockList::initialize_mg(unsigned int n_blocks, const ITERATOR begin, const typename identity<ITERATOR>::type end)
+{
+  index_sets.resize(n_blocks);
+  std::vector<unsigned int> indices;
+  unsigned int k = 0;
+  for (ITERATOR cell = begin; cell != end; ++cell, ++k)
+    {
+      indices.resize(cell->get_fe().dofs_per_cell);
+      cell->get_mg_dof_indices(indices);
+
+      for (std::vector<unsigned int>::const_iterator i=indices.begin();
+          i != indices.end(); ++i)
+       index_sets[k].push_back(*i);
+    }
+}
+
+
+template <typename ITERATOR>
+inline
+void
+BlockList::initialize(
+  unsigned int n_blocks,
+  const ITERATOR begin,
+  const typename identity<ITERATOR>::type end,
+  const std::vector<bool>& selected_dofs,
+  unsigned int offset)
+{
+  index_sets.resize(n_blocks);
+  std::vector<unsigned int> indices;
+  unsigned int k = 0;
+  for (ITERATOR cell = begin; cell != end; ++cell, ++k)
+    {
+      indices.resize(cell->get_fe().dofs_per_cell);
+      AssertDimension(selected_dofs.size(), cell->get_fe().dofs_per_cell);
+      
+      cell->get_dof_indices(indices);
+      
+      for (unsigned int i=0; i<indices.size(); ++i)
+       if (selected_dofs[i])
+         index_sets[k].push_back(indices[i]-offset);
+    }
+}
+
+
+template <typename ITERATOR>
+inline
+void
+BlockList::initialize_mg(
+  unsigned int n_blocks,
+  const ITERATOR begin,
+  const typename identity<ITERATOR>::type end,
+  const std::vector<bool>& selected_dofs,
+  unsigned int offset)
+{
+  index_sets.resize(n_blocks);
+  std::vector<unsigned int> indices;
+  unsigned int k = 0;
+  for (ITERATOR cell = begin; cell != end; ++cell, ++k)
+    {
+      indices.resize(cell->get_fe().dofs_per_cell);
+      AssertDimension(selected_dofs.size(), cell->get_fe().dofs_per_cell);
+      
+      cell->get_mg_dof_indices(indices);
+      
+      for (unsigned int i=0; i<indices.size(); ++i)
+       if (selected_dofs[i])
+         index_sets[k].push_back(indices[i]-offset);
+    }
+}
+
+
+inline
+unsigned int
+BlockList::size() const
+{
+  return index_sets.size();
+}
+
+
+inline
+unsigned int
+BlockList::block_size(unsigned int block) const
+{
+  return index_sets[block].size();
+}
+
+
+inline
+BlockList::const_iterator
+BlockList::begin(unsigned int block) const
+{
+  AssertIndexRange(block, index_sets.size());
+  return index_sets[block].begin();
+}
+
+
+inline
+BlockList::const_iterator
+BlockList::end(unsigned int block) const
+{
+  AssertIndexRange(block, index_sets.size());
+  return index_sets[block].end();
+}
+
+
+inline
+unsigned int
+BlockList::local_index(unsigned int block, unsigned int index) const
+{
+  AssertIndexRange(block, index_sets.size());
+  const block_container& b = index_sets[block];
+  for (unsigned i=0;i<b.size();++i)
+    if (b[i] == index)
+      return i;
+  return numbers::invalid_unsigned_int;
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 100be2e1d824473ccf5301a5cef0e43788977ac1..631e718bc117b73f1348f2b1b5a48abc9d145bfa 100644 (file)
@@ -699,7 +699,13 @@ class PreconditionBlockJacobi : public virtual Subscriptor,
  * elements outside the diagonal blocks may be distributed
  * arbitrarily.
  *
- * See PreconditionBlock for requirements on the matrix.
+ * See PreconditionBlock for requirements on the matrix. The blocks
+ * used in this class must be contiguous and non-overlapping. An
+ * overlapping Schwarz relaxation method can be found in
+ * RelaxationBlockSOR; that class does not offer preconditioning,
+ * though.
+ *
+ * <h3>Permutations</h3>
  *
  * Optionally, the entries of the source vector can be treated in the
  * order of the indices in the permutation vector set by
@@ -713,6 +719,8 @@ class PreconditionBlockJacobi : public virtual Subscriptor,
  * may not change the order inside blocks or swap single indices
  * between blocks.
  *
+ * <h3>Instantiations</h3>
+ *
  * @note Instantiations for this template are provided for <tt>@<float@> and
  * @<double@></tt>; others can be generated in application programs (see the
  * section on @ref Instantiations in the manual).
@@ -918,7 +926,11 @@ class PreconditionBlockSOR : public virtual Subscriptor,
  * based on the implementation in PreconditionBlockSOR.  This
  * class requires storage of the diagonal blocks and their inverses.
  *
- * See PreconditionBlock for requirements on the matrix.
+ * See PreconditionBlock for requirements on the matrix. The blocks
+ * used in this class must be contiguous and non-overlapping. An
+ * overlapping Schwarz relaxation method can be found in
+ * RelaxationBlockSSOR; that class does not offer preconditioning,
+ * though.
  *
  * @note Instantiations for this template are provided for <tt>@<float@> and
  * @<double@></tt>; others can be generated in application programs (see the
diff --git a/deal.II/lac/include/lac/relaxation_block.h b/deal.II/lac/include/lac/relaxation_block.h
new file mode 100644 (file)
index 0000000..32c9b39
--- /dev/null
@@ -0,0 +1,351 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2010 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__relaxation_block_h
+#define __deal2__relaxation_block_h
+
+#include <base/subscriptor.h>
+#include <base/smartpointer.h>
+#include <lac/vector.h>
+#include <lac/precondition_block_base.h>
+#include <lac/block_list.h>
+
+#include <vector>
+#include <set>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * Base class for the implementation of overlapping, multiplicative
+ * Schwarz relaxation methods and smoothers.
+ *
+ * This class uses the infrastructure provided by
+ * PreconditionBlockBase. It adds functions to initialize with a block
+ * list and to do the relaxation step. The actual relaxation method
+ * with the interface expected by SolverRelaxation and
+ * MGSmootherRelaxation is in the derived classes.
+ *
+ * This class allows for more general relaxation methods than
+ * PreconditionBlock, since the index sets may be arbitrary and
+ * overlapping, while there only contiguous, disjoint sets of equal
+ * size are allowed. As a drawback, this class cannot be used as a
+ * preconditioner, since its implementation relies on a straight
+ * forward implementation of the Gauss-Seidel process.
+ *
+ * @author Guido Kanschat
+ * @date 2010
+ */
+template <class MATRIX, typename inverse_type=typename MATRIX::value_type>
+class RelaxationBlock :
+  protected PreconditionBlockBase<inverse_type>
+{
+  private:
+                                    /**
+                                     * Define number type of matrix.
+                                     */
+    typedef typename MATRIX::value_type number;
+
+                                    /**
+                                     * Value type for inverse matrices.
+                                     */
+    typedef inverse_type value_type;
+
+  public:
+                                    /**
+                                     * Parameters for block
+                                     * relaxation methods.
+                                     */
+    class AdditionalData
+    {
+      public:
+                                        /**
+                                         * Default constructor
+                                         */
+       AdditionalData ();
+       
+                                        /**
+                                         * Constructor.
+                                         */
+       AdditionalData (const BlockList& block_list,
+                       const double relaxation = 1.,
+                       const bool invert_diagonal = true,
+                       const bool same_diagonal = false);
+       
+                                        /**
+                                         * Pointer to the DoFHandler
+                                         * providing the cells.
+                                         */
+       SmartPointer<const BlockList, typename RelaxationBlock<MATRIX, inverse_type>::AdditionalData> block_list;
+       
+                                        /**
+                                         * Relaxation parameter.
+                                         */
+       double relaxation;
+
+                                        /**
+                                         * Invert diagonal during initialization.
+                                         */
+       bool invert_diagonal;
+
+                                        /**
+                                         * Assume all diagonal blocks
+                                         * are equal to save memory.
+                                         */
+       bool same_diagonal;
+    };
+    
+                                    /**
+                                     * Initialize matrix and block
+                                     * size.  We store the matrix and
+                                     * the block size in the
+                                     * preconditioner object. In a
+                                     * second step, the inverses of
+                                     * the diagonal blocks may be
+                                     * computed.
+                                     *
+                                     * Additionally, a relaxation
+                                     * parameter for derived classes
+                                     * may be provided.
+                                     */
+    void initialize (const MATRIX& A,
+                    const AdditionalData parameters);
+
+                                    /**
+                                     * Deletes the inverse diagonal
+                                     * block matrices if existent,
+                                     * sets the blocksize to 0, hence
+                                     * leaves the class in the state
+                                     * that it had directly after
+                                     * calling the constructor.
+                                     */
+    void clear();
+
+                                    /**
+                                     * Checks whether the object is empty.
+                                     */
+    bool empty () const;
+
+                                    /**
+                                     * Read-only access to entries.
+                                     * This function is only possible
+                                     * if the inverse diagonal blocks
+                                     * are stored.
+                                     */
+    value_type el(unsigned int i,
+                 unsigned int j) const;
+
+                                    /**
+                                     * Stores the inverse of the
+                                     * diagonal blocks in
+                                     * @p inverse. This costs some
+                                     * additional memory - for DG
+                                     * methods about 1/3 (for double
+                                     * inverses) or 1/6 (for float
+                                     * inverses) of that used for the
+                                     * matrix - but it makes the
+                                     * preconditioning much faster.
+                                     *
+                                     * It is not allowed to call this
+                                     * function twice (will produce
+                                     * an error) before a call of
+                                     * <tt>clear(...)</tt>  because at the
+                                     * second time there already
+                                     * exist the inverse matrices.
+                                     *
+                                     * After this function is called,
+                                     * the lock on the matrix given
+                                     * through the @p use_matrix
+                                     * function is released, i.e. you
+                                     * may overwrite of delete it.
+                                     * You may want to do this in
+                                     * case you use this matrix to
+                                     * precondition another matrix.
+                                     */
+    void invert_diagblocks();
+
+  protected:
+                                    /**
+                                     * Perform one block relaxation
+                                     * step.
+                                     *
+                                     * Depending on the arguments @p
+                                     * dst and @p pref, this performs
+                                     * an SOR step (both reference
+                                     * the same vector) of a Jacobi
+                                     * step (both are different
+                                     * vectors). For the Jacobi step,
+                                     * the calling function must copy
+                                     * @p dst to @p pref after this.
+                                     */
+    template <typename number2>
+    void do_step (
+      Vector<number2>       &dst,
+      const Vector<number2> &prev,
+      const Vector<number2> &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 preconditoining
+                                     * @p vmult function of the
+                                     * derived classes.
+                                     */
+    SmartPointer<const MATRIX,RelaxationBlock<MATRIX,inverse_type> > A;
+
+                                    /**
+                                     * Control information.
+                                     */
+    AdditionalData additional_data;
+};
+
+
+/**
+ * Block Gauss-Seidel method with possibly overlapping blocks.
+ *
+ * This class implements the step() and Tstep() functions expected by
+ * SolverRelaxation and MGSmootherRelaxation. They perform a
+ * multiplicative Schwarz method on the blocks provided in the
+ * BlockList of AdditionalData. Differing from PreconditionBlockSOR,
+ * these blocks may be of varying size, non-contiguous, and
+ * overlapping. On the other hand, this class does not implement the
+ * preconditioner interface expected by Solver objects.
+ *
+ * @author Guido Kanschat
+ * @data 2010
+ */
+template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
+class RelaxationBlockSOR : public virtual Subscriptor,
+                            protected RelaxationBlock<MATRIX, inverse_type>
+{
+  public:
+                                    /**
+                                     * Default constructor.
+                                     */
+//    RelaxationBlockSOR();
+    
+                                    /**
+                                     * Define number type of matrix.
+                                     */
+    typedef typename MATRIX::value_type number;
+
+                                    /**
+                                     * Make type publicly available.
+                                     */
+    using RelaxationBlock<MATRIX,inverse_type>::AdditionalData;
+
+                                    /**
+                                     * Make initialization function
+                                     * publicly available.
+                                     */
+    using RelaxationBlock<MATRIX, inverse_type>::initialize;
+
+                                    /**
+                                     * Make function of base class public again.
+                                     */
+    using RelaxationBlock<MATRIX, inverse_type>::clear;
+
+                                    /**
+                                     * Make function of base class public again.
+                                     */
+    using RelaxationBlock<MATRIX, inverse_type>::empty;
+                                    /**
+                                     * Perform one step of the SOR
+                                     * iteration.
+                                     */
+    template <typename number2>
+    void step (Vector<number2>& dst, const Vector<number2>& rhs) const;
+
+                                    /**
+                                     * Perform one step of the
+                                     * transposed SOR iteration.
+                                     */
+    template <typename number2>
+    void Tstep (Vector<number2>& dst, const Vector<number2>& rhs) const;
+
+  protected:
+                                    /**
+                                     * Constructor to be used by
+                                     * RelaxationBlockSSOR.
+                                     */
+//    RelaxationBlockSOR(bool store);
+};
+
+
+/**
+ * Symmetric block Gauss-Seidel method with possibly overlapping blocks.
+ *
+ * This class implements the step() and Tstep() functions expected by
+ * SolverRelaxation and MGSmootherRelaxation. They perform a
+ * multiplicative Schwarz method on the blocks provided in the
+ * BlockList of AdditionalData in symmetric fashion. Differing from
+ * PreconditionBlockSSOR, these blocks may be of varying size,
+ * non-contiguous, and overlapping. On the other hand, this class does
+ * not implement the preconditioner interface expected by Solver
+ * objects.
+ *
+ * @author Guido Kanschat
+ * @data 2010
+ */
+template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
+class RelaxationBlockSSOR : public virtual Subscriptor,
+                           protected RelaxationBlock<MATRIX, inverse_type>
+{
+  public:    
+                                    /**
+                                     * Define number type of matrix.
+                                     */
+    typedef typename MATRIX::value_type number;
+
+                                    /**
+                                     * Make type publicly available.
+                                     */
+    using RelaxationBlock<MATRIX,inverse_type>::AdditionalData;
+
+                                    /**
+                                     * Make initialization function
+                                     * publicly available.
+                                     */
+    using RelaxationBlock<MATRIX, inverse_type>::initialize;
+
+                                    /**
+                                     * Make function of base class public again.
+                                     */
+    using RelaxationBlock<MATRIX, inverse_type>::clear;
+
+                                    /**
+                                     * Make function of base class public again.
+                                     */
+    using RelaxationBlock<MATRIX, inverse_type>::empty;
+                                    /**
+                                     * Perform one step of the SOR
+                                     * iteration.
+                                     */
+    template <typename number2>
+    void step (Vector<number2>& dst, const Vector<number2>& rhs) const;
+
+                                    /**
+                                     * Perform one step of the
+                                     * transposed SOR iteration.
+                                     */
+    template <typename number2>
+    void Tstep (Vector<number2>& dst, const Vector<number2>& rhs) const;
+};
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/lac/include/lac/relaxation_block.templates.h b/deal.II/lac/include/lac/relaxation_block.templates.h
new file mode 100644 (file)
index 0000000..4206ff1
--- /dev/null
@@ -0,0 +1,222 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 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__relaxation_block_templates_h
+#define __deal2__relaxation_block_templates_h
+
+#include <lac/relaxation_block.h>
+#include <lac/full_matrix.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+template <class MATRIX, typename inverse_type>
+inline
+RelaxationBlock<MATRIX,inverse_type>::AdditionalData::AdditionalData ()
+               :
+               relaxation(1.),
+               invert_diagonal(true),
+               same_diagonal(false)
+{}
+
+
+template <class MATRIX, typename inverse_type>
+inline
+RelaxationBlock<MATRIX,inverse_type>::AdditionalData::AdditionalData (
+  const BlockList& bl,
+  const double relaxation,
+  const bool invert_diagonal,
+  const bool same_diagonal)
+               :
+               block_list(&bl),
+               relaxation(relaxation),
+               invert_diagonal(invert_diagonal),
+               same_diagonal(same_diagonal)
+{}
+
+
+template <class MATRIX, typename inverse_type>
+inline
+void
+RelaxationBlock<MATRIX,inverse_type>::initialize (
+  const MATRIX& M,
+  const AdditionalData parameters)
+{
+  Assert (parameters.invert_diagonal, ExcNotImplemented());
+  
+  clear();
+//  Assert (M.m() == M.n(), ExcNotQuadratic());
+  A = &M;
+  additional_data = parameters;
+  
+  this->reinit(additional_data.block_list->size(), 0, additional_data.same_diagonal);
+  
+  if (parameters.invert_diagonal)
+    invert_diagblocks();
+}
+
+
+template <class MATRIX, typename inverse_type>
+inline
+void
+RelaxationBlock<MATRIX,inverse_type>::clear ()
+{
+  A = 0;
+  additional_data.block_list = 0;
+  PreconditionBlockBase<inverse_type>::clear ();
+}
+
+
+template <class MATRIX, typename inverse_type>
+inline
+void
+RelaxationBlock<MATRIX,inverse_type>::invert_diagblocks ()
+{
+  const MATRIX &M=*A;
+  FullMatrix<inverse_type> M_cell;
+
+  if (this->same_diagonal())
+    {
+      Assert(false, ExcNotImplemented());
+    }
+  else
+    {
+      for (unsigned int block=0;block<additional_data.block_list->size();++block)
+       {
+         const unsigned int bs = additional_data.block_list->block_size(block);
+         M_cell.reinit(bs, bs);
+
+                                          // Copy rows for this block
+                                          // into the matrix for the
+                                          // diagonal block
+         BlockList::const_iterator row = additional_data.block_list->begin(block);
+         BlockList::const_iterator end = additional_data.block_list->end(block);
+         for (unsigned int row_cell=0; row_cell<bs; ++row_cell, ++row)
+           {
+//TODO:[GK] Optimize here
+             for (typename MATRIX::const_iterator entry = M.begin(*row);
+                  entry != M.end(*row); ++entry)
+               {
+                 const unsigned int column = entry->column();
+                 const unsigned int col_cell = additional_data.block_list->local_index(block, column);
+                 if (col_cell != numbers::invalid_unsigned_int)
+                   M_cell(row_cell, col_cell) = entry->value();
+               }
+           }
+                                          // Now M_cell contains the
+                                          // diagonal block. Now
+                                          // store it and its
+                                          // inverse, if so requested.
+         if (this->store_diagonals())
+           {
+             this->diagonal(block).reinit(bs, bs);
+             this->diagonal(block) = M_cell;
+           }
+         this->inverse(block).reinit(bs, bs);
+         this->inverse(block).invert(M_cell);
+       }
+    }
+  this->inverses_computed(true);
+}
+
+
+template <class MATRIX, typename inverse_type>
+template <typename number2>
+inline
+void
+RelaxationBlock<MATRIX,inverse_type>::do_step (
+  Vector<number2>       &dst,
+  const Vector<number2> &prev,
+  const Vector<number2> &src,
+  const bool backward) const
+{
+  Assert (additional_data.invert_diagonal, ExcNotImplemented());
+  
+  const MATRIX &M=*this->A;
+  Vector<number2> b_cell, x_cell;
+  
+  const unsigned int n_blocks = additional_data.block_list->size();
+  for (unsigned int bi=0;bi<n_blocks;++bi)
+    {
+      const unsigned int block = backward ? (n_blocks - bi - 1) : bi;
+      const unsigned int bs = additional_data.block_list->block_size(block);
+      
+      b_cell.reinit(bs);
+      x_cell.reinit(bs);
+                                      // Collect off-diagonal parts
+      BlockList::const_iterator row = additional_data.block_list->begin(block);
+      BlockList::const_iterator end = additional_data.block_list->end(block);
+      for (unsigned int row_cell=0; row_cell<bs; ++row_cell, ++row)
+       {
+         b_cell(row_cell) = src(*row);
+         for (typename MATRIX::const_iterator entry = M.begin(*row);
+              entry != M.end(*row); ++entry)
+           b_cell(row_cell) -= entry->value() * prev(entry->column());
+       }
+                                      // Apply inverse diagonal
+      this->inverse(block).vmult(x_cell, b_cell);
+                                      // Store in result vector
+      row=additional_data.block_list->begin(block);
+      for (unsigned int row_cell=0; row_cell<bs; ++row_cell, ++row)
+       dst(*row) = prev(*row) + additional_data.relaxation * x_cell(row_cell);
+    }
+}
+
+
+//----------------------------------------------------------------------//
+
+template <class MATRIX, typename inverse_type>
+template <typename number2>
+void RelaxationBlockSOR<MATRIX,inverse_type>::step (
+  Vector<number2>       &dst,
+  const Vector<number2> &src) const
+{
+  this->do_step(dst, dst, src, false);
+}
+
+
+template <class MATRIX, typename inverse_type>
+template <typename number2>
+void RelaxationBlockSOR<MATRIX,inverse_type>::Tstep (
+  Vector<number2>       &dst,
+  const Vector<number2> &src) const
+{
+  this->do_step(dst, dst, src, true);
+}
+
+
+template <class MATRIX, typename inverse_type>
+template <typename number2>
+void RelaxationBlockSSOR<MATRIX,inverse_type>::step (
+  Vector<number2>       &dst,
+  const Vector<number2> &src) const
+{
+  this->do_step(dst, dst, src, true);
+  this->do_step(dst, dst, src, false);
+}
+
+
+template <class MATRIX, typename inverse_type>
+template <typename number2>
+void RelaxationBlockSSOR<MATRIX,inverse_type>::Tstep (
+  Vector<number2>       &dst,
+  const Vector<number2> &src) const
+{
+  this->do_step(dst, dst, src, false);
+  this->do_step(dst, dst, src, true);
+}
+
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+#endif
diff --git a/deal.II/lac/source/relaxation_block.cc b/deal.II/lac/source/relaxation_block.cc
new file mode 100644 (file)
index 0000000..26906dc
--- /dev/null
@@ -0,0 +1,20 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2010 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/relaxation_block.templates.h>
+#include <lac/sparse_matrix.h>
+
+DEAL_II_NAMESPACE_OPEN
+#include "relaxation_block.inst"
+DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/lac/source/relaxation_block.inst.in b/deal.II/lac/source/relaxation_block.inst.in
new file mode 100644 (file)
index 0000000..5e09c10
--- /dev/null
@@ -0,0 +1,40 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2005, 2006, 2007, 2009, 2010 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.
+//
+//---------------------------------------------------------------------------
+
+
+for (S1, S2 : REAL_SCALARS)
+  {
+    template class RelaxationBlock<SparseMatrix<S1>, S2>;
+    template class RelaxationBlockSOR<SparseMatrix<S1>, S2>;
+    template class RelaxationBlockSSOR<SparseMatrix<S1>, S2>;
+  }
+
+
+for (S1, S2, S3 : REAL_SCALARS)
+  {
+// ------------ 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;
+  }
+

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.