<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>
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
* 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
* 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).
* 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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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;
+ }
+