From 28b22c1df83cf3f14b9bf60df5c363f72187b9a9 Mon Sep 17 00:00:00 2001 From: kanschat Date: Tue, 19 Oct 2010 19:20:22 +0000 Subject: [PATCH] new classes BlockList, RelaxationBlockSOR and RelaxationBlockSSOR git-svn-id: https://svn.dealii.org/trunk@22387 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 11 +- deal.II/lac/include/lac/block_list.h | 338 +++++++++++++++++ deal.II/lac/include/lac/precondition_block.h | 16 +- deal.II/lac/include/lac/relaxation_block.h | 351 ++++++++++++++++++ .../include/lac/relaxation_block.templates.h | 222 +++++++++++ deal.II/lac/source/relaxation_block.cc | 20 + deal.II/lac/source/relaxation_block.inst.in | 40 ++ 7 files changed, 995 insertions(+), 3 deletions(-) create mode 100644 deal.II/lac/include/lac/block_list.h create mode 100644 deal.II/lac/include/lac/relaxation_block.h create mode 100644 deal.II/lac/include/lac/relaxation_block.templates.h create mode 100644 deal.II/lac/source/relaxation_block.cc create mode 100644 deal.II/lac/source/relaxation_block.inst.in diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 31c376625c..f049643af7 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -252,7 +252,16 @@ through DoFHandler::get_tria() and DoFHandler::get_fe(), respectively.

lac

    -
  1. New: The ConstraintMatrix::merge function now takes a second + +

  2. 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. +
    + (GK 2010/10/19) +

  3. + +
  4. 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.
    diff --git a/deal.II/lac/include/lac/block_list.h b/deal.II/lac/include/lac/block_list.h new file mode 100644 index 0000000000..e0cc1248a0 --- /dev/null +++ b/deal.II/lac/include/lac/block_list.h @@ -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 + +#include +#include + +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 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 + * begin->get_dof_indices() + * 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 + void initialize(unsigned int n_blocks, + const ITERATOR begin, + const typename identity::type end); + + /** + * Set up all index sets using an + * DoF iterator range. This + * function will call + * begin->get_mg_dof_indices() + * 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 + void initialize_mg(unsigned int n_blocks, + const ITERATOR begin, + const typename identity::type end); + + /** + * Set up all index sets using an + * DoF iterator range. This + * function will call + * begin->get_dof_indices() + * with a signature like + * DoFCellAccessor::get_dof_indices(). + * Typically, the iterators will + * loop over active cells of a + * triangulation. + * + * The argument vector + * selected_dofs 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 offset is + * the start index of this block. + * + * In addition, the function + * needs the total number of + * blocks as its first argument. + */ + template + void initialize(unsigned int n_blocks, + const ITERATOR begin, + const typename identity::type end, + const std::vector& selected_dofs, + unsigned int offset = 0); + /** + * Set up all index sets using an + * DoF iterator range. This + * function will call + * begin->get_mg_dof_indices() + * 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 + * selected_dofs 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 offset is + * the start index of this block. + * + * In addition, the function + * needs the total number of + * blocks as its first argument. + */ + template + void initialize_mg(unsigned int n_blocks, + const ITERATOR begin, + const typename identity::type end, + const std::vector& 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 + * index in + * block, 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 index_sets; +}; + + +template +inline +void +BlockList::initialize(unsigned int n_blocks, const ITERATOR begin, const typename identity::type end) +{ + index_sets.resize(n_blocks); + std::vector 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::const_iterator i=indices.begin(); + i != indices.end(); ++i) + index_sets[k].push_back(*i); + } +} + + +template +inline +void +BlockList::initialize_mg(unsigned int n_blocks, const ITERATOR begin, const typename identity::type end) +{ + index_sets.resize(n_blocks); + std::vector 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::const_iterator i=indices.begin(); + i != indices.end(); ++i) + index_sets[k].push_back(*i); + } +} + + +template +inline +void +BlockList::initialize( + unsigned int n_blocks, + const ITERATOR begin, + const typename identity::type end, + const std::vector& selected_dofs, + unsigned int offset) +{ + index_sets.resize(n_blocks); + std::vector 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 +inline +void +BlockList::initialize_mg( + unsigned int n_blocks, + const ITERATOR begin, + const typename identity::type end, + const std::vector& selected_dofs, + unsigned int offset) +{ + index_sets.resize(n_blocks); + std::vector 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; iPermutations * * 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. * + *

    Instantiations

    + * * @note Instantiations for this template are provided for @ and * @; 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 @ and * @; 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 index 0000000000..32c9b39bd3 --- /dev/null +++ b/deal.II/lac/include/lac/relaxation_block.h @@ -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 +#include +#include +#include +#include + +#include +#include + +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 RelaxationBlock : + protected PreconditionBlockBase +{ + 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::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 + * clear(...) 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 + void do_step ( + Vector &dst, + const Vector &prev, + const Vector &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 > 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 RelaxationBlockSOR : public virtual Subscriptor, + protected RelaxationBlock +{ + public: + /** + * Default constructor. + */ +// RelaxationBlockSOR(); + + /** + * Define number type of matrix. + */ + typedef typename MATRIX::value_type number; + + /** + * Make type publicly available. + */ + using RelaxationBlock::AdditionalData; + + /** + * Make initialization function + * publicly available. + */ + using RelaxationBlock::initialize; + + /** + * Make function of base class public again. + */ + using RelaxationBlock::clear; + + /** + * Make function of base class public again. + */ + using RelaxationBlock::empty; + /** + * Perform one step of the SOR + * iteration. + */ + template + void step (Vector& dst, const Vector& rhs) const; + + /** + * Perform one step of the + * transposed SOR iteration. + */ + template + void Tstep (Vector& dst, const Vector& 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 RelaxationBlockSSOR : public virtual Subscriptor, + protected RelaxationBlock +{ + public: + /** + * Define number type of matrix. + */ + typedef typename MATRIX::value_type number; + + /** + * Make type publicly available. + */ + using RelaxationBlock::AdditionalData; + + /** + * Make initialization function + * publicly available. + */ + using RelaxationBlock::initialize; + + /** + * Make function of base class public again. + */ + using RelaxationBlock::clear; + + /** + * Make function of base class public again. + */ + using RelaxationBlock::empty; + /** + * Perform one step of the SOR + * iteration. + */ + template + void step (Vector& dst, const Vector& rhs) const; + + /** + * Perform one step of the + * transposed SOR iteration. + */ + template + void Tstep (Vector& dst, const Vector& 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 index 0000000000..4206ff170b --- /dev/null +++ b/deal.II/lac/include/lac/relaxation_block.templates.h @@ -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 +#include + +DEAL_II_NAMESPACE_OPEN + +template +inline +RelaxationBlock::AdditionalData::AdditionalData () + : + relaxation(1.), + invert_diagonal(true), + same_diagonal(false) +{} + + +template +inline +RelaxationBlock::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 +inline +void +RelaxationBlock::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 +inline +void +RelaxationBlock::clear () +{ + A = 0; + additional_data.block_list = 0; + PreconditionBlockBase::clear (); +} + + +template +inline +void +RelaxationBlock::invert_diagblocks () +{ + const MATRIX &M=*A; + FullMatrix M_cell; + + if (this->same_diagonal()) + { + Assert(false, ExcNotImplemented()); + } + else + { + for (unsigned int block=0;blocksize();++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_cellcolumn(); + 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 +template +inline +void +RelaxationBlock::do_step ( + Vector &dst, + const Vector &prev, + const Vector &src, + const bool backward) const +{ + Assert (additional_data.invert_diagonal, ExcNotImplemented()); + + const MATRIX &M=*this->A; + Vector b_cell, x_cell; + + const unsigned int n_blocks = additional_data.block_list->size(); + for (unsigned int bi=0;biblock_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_cellvalue() * 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 +template +void RelaxationBlockSOR::step ( + Vector &dst, + const Vector &src) const +{ + this->do_step(dst, dst, src, false); +} + + +template +template +void RelaxationBlockSOR::Tstep ( + Vector &dst, + const Vector &src) const +{ + this->do_step(dst, dst, src, true); +} + + +template +template +void RelaxationBlockSSOR::step ( + Vector &dst, + const Vector &src) const +{ + this->do_step(dst, dst, src, true); + this->do_step(dst, dst, src, false); +} + + +template +template +void RelaxationBlockSSOR::Tstep ( + Vector &dst, + const Vector &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 index 0000000000..26906dcde6 --- /dev/null +++ b/deal.II/lac/source/relaxation_block.cc @@ -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 +#include + +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 index 0000000000..5e09c1050c --- /dev/null +++ b/deal.II/lac/source/relaxation_block.inst.in @@ -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, S2>; + template class RelaxationBlockSOR, S2>; + template class RelaxationBlockSSOR, S2>; + } + + +for (S1, S2, S3 : REAL_SCALARS) + { +// ------------ RelaxationBlockSOR ----------------- + template + void RelaxationBlockSOR, S2>::step + (Vector &, const Vector &) const; + template + void RelaxationBlockSOR, S2>::Tstep + (Vector &, const Vector &) const; + +// ------------ RelaxationBlockSSOR ----------------- + template + void RelaxationBlockSSOR, S2>::step + (Vector &, const Vector &) const; + template + void RelaxationBlockSSOR, S2>::Tstep + (Vector &, const Vector &) const; + } + -- 2.39.5