From 0db72a2de2259a6f2648fc2663c65613469f9291 Mon Sep 17 00:00:00 2001 From: heister Date: Wed, 3 Sep 2008 14:22:26 +0000 Subject: [PATCH] Added new class CompressedSimpleSparsityPattern and mentioned my contribution. git-svn-id: https://svn.dealii.org/trunk@16733 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/source/dofs/dof_constraints.cc | 10 + deal.II/deal.II/source/dofs/dof_tools.cc | 149 +++++- deal.II/doc/authors.html | 3 + deal.II/doc/doxygen/headers/matrices.h | 42 +- .../lac/include/lac/block_sparsity_pattern.h | 126 ++++- .../lac/compressed_set_sparsity_pattern.h | 3 + .../lac/compressed_simple_sparsity_pattern.h | 465 ++++++++++++++++++ .../include/lac/compressed_sparsity_pattern.h | 7 +- deal.II/lac/include/lac/sparsity_pattern.h | 17 +- deal.II/lac/source/block_sparsity_pattern.cc | 56 +++ .../compressed_simple_sparsity_pattern.cc | 241 +++++++++ deal.II/lac/source/sparsity_pattern.cc | 45 ++ 12 files changed, 1136 insertions(+), 28 deletions(-) create mode 100644 deal.II/lac/include/lac/compressed_simple_sparsity_pattern.h create mode 100644 deal.II/lac/source/compressed_simple_sparsity_pattern.cc diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index 748bba78ec..342d9ff716 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -1809,6 +1810,11 @@ add_entries_local_to_global (const std::vector (const std::vector &, + CompressedSimpleSparsityPattern &, + const bool) const; + template void ConstraintMatrix:: add_entries_local_to_global (const std::vector &, BlockSparsityPattern &, @@ -1822,5 +1828,9 @@ add_entries_local_to_global (const std::vecto BlockCompressedSetSparsityPattern &, const bool) const; +template void ConstraintMatrix:: +add_entries_local_to_global (const std::vector &, + BlockCompressedSimpleSparsityPattern &, + const bool) const; DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 998f451bed..fbbac19936 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -31,6 +31,7 @@ #include #include #include +#include #include #include @@ -5078,6 +5079,12 @@ DoFTools::make_sparsity_pattern, CompressedSetSparsityPattern &sparsity, const ConstraintMatrix &); template void +DoFTools::make_sparsity_pattern, + CompressedSimpleSparsityPattern> +(const DoFHandler &dof, + CompressedSimpleSparsityPattern &sparsity, + const ConstraintMatrix &); +template void DoFTools::make_sparsity_pattern, BlockSparsityPattern> (const DoFHandler &dof, @@ -5095,6 +5102,12 @@ DoFTools::make_sparsity_pattern, (const DoFHandler &dof, BlockCompressedSetSparsityPattern &sparsity, const ConstraintMatrix &); +template void +DoFTools::make_sparsity_pattern, + BlockCompressedSimpleSparsityPattern> +(const DoFHandler &dof, + BlockCompressedSimpleSparsityPattern &sparsity, + const ConstraintMatrix &); template void DoFTools::make_sparsity_pattern, @@ -5118,6 +5131,13 @@ DoFTools::make_sparsity_pattern, CompressedSetSparsityPattern &sparsity, const ConstraintMatrix &); +template void +DoFTools::make_sparsity_pattern, + CompressedSimpleSparsityPattern> +(const hp::DoFHandler &dof, + CompressedSimpleSparsityPattern &sparsity, + const ConstraintMatrix &); + template void DoFTools::make_sparsity_pattern, @@ -5137,6 +5157,12 @@ DoFTools::make_sparsity_pattern, (const hp::DoFHandler &dof, BlockCompressedSetSparsityPattern &sparsity, const ConstraintMatrix &); +template void +DoFTools::make_sparsity_pattern, + BlockCompressedSimpleSparsityPattern> +(const hp::DoFHandler &dof, + BlockCompressedSimpleSparsityPattern &sparsity, + const ConstraintMatrix &); template void @@ -5158,6 +5184,12 @@ DoFTools::make_sparsity_pattern, const Table<2,Coupling>&, CompressedSetSparsityPattern&); template void +DoFTools::make_sparsity_pattern, + CompressedSimpleSparsityPattern> +(const DoFHandler&, + const Table<2,Coupling>&, + CompressedSimpleSparsityPattern&); +template void DoFTools::make_sparsity_pattern, BlockSparsityPattern> (const DoFHandler&, @@ -5175,7 +5207,12 @@ DoFTools::make_sparsity_pattern, (const DoFHandler&, const Table<2,Coupling>&, BlockCompressedSetSparsityPattern&); - +template void +DoFTools::make_sparsity_pattern, + BlockCompressedSimpleSparsityPattern> +(const DoFHandler&, + const Table<2,Coupling>&, + BlockCompressedSimpleSparsityPattern&); template void DoFTools::make_sparsity_pattern, @@ -5196,6 +5233,12 @@ DoFTools::make_sparsity_pattern, const Table<2,Coupling>&, CompressedSetSparsityPattern&); template void +DoFTools::make_sparsity_pattern, + CompressedSimpleSparsityPattern> +(const hp::DoFHandler&, + const Table<2,Coupling>&, + CompressedSimpleSparsityPattern&); +template void DoFTools::make_sparsity_pattern, BlockSparsityPattern> (const hp::DoFHandler&, @@ -5213,6 +5256,12 @@ DoFTools::make_sparsity_pattern, (const hp::DoFHandler&, const Table<2,Coupling>&, BlockCompressedSetSparsityPattern&); +template void +DoFTools::make_sparsity_pattern, + BlockCompressedSimpleSparsityPattern> +(const hp::DoFHandler&, + const Table<2,Coupling>&, + BlockCompressedSimpleSparsityPattern&); template void @@ -5234,6 +5283,12 @@ DoFTools::make_sparsity_pattern, const DoFHandler &dof_col, CompressedSetSparsityPattern &sparsity); template void +DoFTools::make_sparsity_pattern, + CompressedSimpleSparsityPattern> +(const DoFHandler &dof_row, + const DoFHandler &dof_col, + CompressedSimpleSparsityPattern &sparsity); +template void DoFTools::make_sparsity_pattern, BlockSparsityPattern> (const DoFHandler &dof_row, @@ -5251,6 +5306,12 @@ DoFTools::make_sparsity_pattern, (const DoFHandler &dof_row, const DoFHandler &dof_col, BlockCompressedSetSparsityPattern &sparsity); +template void +DoFTools::make_sparsity_pattern, + BlockCompressedSimpleSparsityPattern> +(const DoFHandler &dof_row, + const DoFHandler &dof_col, + BlockCompressedSimpleSparsityPattern &sparsity); template void DoFTools::make_sparsity_pattern, @@ -5271,6 +5332,13 @@ DoFTools::make_sparsity_pattern, const hp::DoFHandler &dof_col, CompressedSetSparsityPattern &sparsity); template void +DoFTools::make_sparsity_pattern, + CompressedSimpleSparsityPattern> +(const hp::DoFHandler &dof_row, + const hp::DoFHandler &dof_col, + CompressedSimpleSparsityPattern &sparsity); + +template void DoFTools::make_sparsity_pattern, BlockSparsityPattern> (const hp::DoFHandler &dof_row, @@ -5288,7 +5356,12 @@ DoFTools::make_sparsity_pattern, (const hp::DoFHandler &dof_row, const hp::DoFHandler &dof_col, BlockCompressedSetSparsityPattern &sparsity); - +template void +DoFTools::make_sparsity_pattern, + BlockCompressedSimpleSparsityPattern> +(const hp::DoFHandler &dof_row, + const hp::DoFHandler &dof_col, + BlockCompressedSimpleSparsityPattern &sparsity); // #if deal_II_dimension > 1 template void @@ -5307,6 +5380,11 @@ DoFTools::make_boundary_sparsity_pattern,Compresse const std::vector &, CompressedSetSparsityPattern &); template void +DoFTools::make_boundary_sparsity_pattern,CompressedSimpleSparsityPattern> +(const DoFHandler& dof, + const std::vector &, + CompressedSimpleSparsityPattern &); +template void DoFTools::make_boundary_sparsity_pattern,BlockSparsityPattern> (const DoFHandler& dof, const std::vector &, @@ -5321,6 +5399,11 @@ DoFTools::make_boundary_sparsity_pattern,BlockComp (const DoFHandler& dof, const std::vector &, BlockCompressedSetSparsityPattern &); +template void +DoFTools::make_boundary_sparsity_pattern,BlockCompressedSimpleSparsityPattern> +(const DoFHandler& dof, + const std::vector &, + BlockCompressedSimpleSparsityPattern &); template void DoFTools::make_boundary_sparsity_pattern,SparsityPattern> @@ -5338,6 +5421,11 @@ DoFTools::make_boundary_sparsity_pattern,Compr const std::vector &, CompressedSetSparsityPattern &); template void +DoFTools::make_boundary_sparsity_pattern,CompressedSimpleSparsityPattern> +(const hp::DoFHandler& dof, + const std::vector &, + CompressedSimpleSparsityPattern &); +template void DoFTools::make_boundary_sparsity_pattern,BlockSparsityPattern> (const hp::DoFHandler& dof, const std::vector &, @@ -5352,6 +5440,11 @@ DoFTools::make_boundary_sparsity_pattern,Block (const hp::DoFHandler& dof, const std::vector &, BlockCompressedSetSparsityPattern &); +template void +DoFTools::make_boundary_sparsity_pattern,BlockCompressedSimpleSparsityPattern> +(const hp::DoFHandler& dof, + const std::vector &, + BlockCompressedSimpleSparsityPattern &); template void @@ -5373,6 +5466,12 @@ DoFTools::make_boundary_sparsity_pattern,Compresse const std::vector &dof_to_boundary_mapping, CompressedSetSparsityPattern &sparsity); template void +DoFTools::make_boundary_sparsity_pattern,CompressedSimpleSparsityPattern> +(const DoFHandler& dof, + const FunctionMap::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + CompressedSimpleSparsityPattern &sparsity); +template void DoFTools::make_boundary_sparsity_pattern,BlockSparsityPattern> (const DoFHandler& dof, const FunctionMap::type &boundary_indicators, @@ -5390,6 +5489,12 @@ DoFTools::make_boundary_sparsity_pattern,BlockComp const FunctionMap::type &boundary_indicators, const std::vector &dof_to_boundary_mapping, BlockCompressedSetSparsityPattern &sparsity); +template void +DoFTools::make_boundary_sparsity_pattern,BlockCompressedSimpleSparsityPattern> +(const DoFHandler& dof, + const FunctionMap::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + BlockCompressedSimpleSparsityPattern &sparsity); template void DoFTools::make_boundary_sparsity_pattern,SparsityPattern> @@ -5410,6 +5515,12 @@ DoFTools::make_boundary_sparsity_pattern,Compr const std::vector &dof_to_boundary_mapping, CompressedSetSparsityPattern &sparsity); template void +DoFTools::make_boundary_sparsity_pattern,CompressedSimpleSparsityPattern> +(const hp::DoFHandler& dof, + const FunctionMap::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + CompressedSimpleSparsityPattern &sparsity); +template void DoFTools::make_boundary_sparsity_pattern,BlockSparsityPattern> (const hp::DoFHandler& dof, const FunctionMap::type &boundary_indicators, @@ -5427,6 +5538,12 @@ DoFTools::make_boundary_sparsity_pattern,Block const FunctionMap::type &boundary_indicators, const std::vector &dof_to_boundary_mapping, BlockCompressedSetSparsityPattern &sparsity); +template void +DoFTools::make_boundary_sparsity_pattern,BlockCompressedSimpleSparsityPattern> +(const hp::DoFHandler& dof, + const FunctionMap::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + BlockCompressedSimpleSparsityPattern &sparsity); template void @@ -5442,6 +5559,10 @@ DoFTools::make_flux_sparsity_pattern,CompressedSet (const DoFHandler &dof, CompressedSetSparsityPattern &sparsity); template void +DoFTools::make_flux_sparsity_pattern,CompressedSimpleSparsityPattern> +(const DoFHandler &dof, + CompressedSimpleSparsityPattern &sparsity); +template void DoFTools::make_flux_sparsity_pattern,BlockSparsityPattern> (const DoFHandler &dof, BlockSparsityPattern &sparsity); @@ -5453,6 +5574,10 @@ template void DoFTools::make_flux_sparsity_pattern,BlockCompressedSetSparsityPattern> (const DoFHandler &dof, BlockCompressedSetSparsityPattern &sparsity); +template void +DoFTools::make_flux_sparsity_pattern,BlockCompressedSimpleSparsityPattern> +(const DoFHandler &dof, + BlockCompressedSimpleSparsityPattern &sparsity); template void DoFTools::make_flux_sparsity_pattern,SparsityPattern> @@ -5467,6 +5592,10 @@ DoFTools::make_flux_sparsity_pattern,Compresse (const hp::DoFHandler &dof, CompressedSetSparsityPattern &sparsity); template void +DoFTools::make_flux_sparsity_pattern,CompressedSimpleSparsityPattern> +(const hp::DoFHandler &dof, + CompressedSimpleSparsityPattern &sparsity); +template void DoFTools::make_flux_sparsity_pattern,BlockSparsityPattern> (const hp::DoFHandler &dof, BlockSparsityPattern &sparsity); @@ -5478,6 +5607,10 @@ template void DoFTools::make_flux_sparsity_pattern,BlockCompressedSetSparsityPattern> (const hp::DoFHandler &dof, BlockCompressedSetSparsityPattern &sparsity); +template void +DoFTools::make_flux_sparsity_pattern,BlockCompressedSimpleSparsityPattern> +(const hp::DoFHandler &dof, + BlockCompressedSimpleSparsityPattern &sparsity); #if deal_II_dimension > 1 @@ -5500,6 +5633,12 @@ DoFTools::make_flux_sparsity_pattern,CompressedSet const Table<2,Coupling>&, const Table<2,Coupling>&); template void +DoFTools::make_flux_sparsity_pattern,CompressedSimpleSparsityPattern> +(const DoFHandler &dof, + CompressedSimpleSparsityPattern &, + const Table<2,Coupling>&, + const Table<2,Coupling>&); +template void DoFTools::make_flux_sparsity_pattern,BlockSparsityPattern> (const DoFHandler &dof, BlockSparsityPattern &, @@ -5517,6 +5656,12 @@ DoFTools::make_flux_sparsity_pattern,BlockCompress BlockCompressedSetSparsityPattern &, const Table<2,Coupling>&, const Table<2,Coupling>&); +template void +DoFTools::make_flux_sparsity_pattern,BlockCompressedSimpleSparsityPattern> +(const DoFHandler &dof, + BlockCompressedSimpleSparsityPattern &, + const Table<2,Coupling>&, + const Table<2,Coupling>&); #endif diff --git a/deal.II/doc/authors.html b/deal.II/doc/authors.html index a3b9acb590..72d326ba7c 100644 --- a/deal.II/doc/authors.html +++ b/deal.II/doc/authors.html @@ -69,6 +69,9 @@ alphabetical order):
  • Ivan Christov: Step-25 tutorial program. +
  • Timo Heister: + CompressedSimpleSparsityPattern class. +
  • Luca Heltai: Gmsh format mesh reader and writer;   some of the meshes in the GridGenerator class;   diff --git a/deal.II/doc/doxygen/headers/matrices.h b/deal.II/doc/doxygen/headers/matrices.h index 3c1ade3c03..ca9f3ae8f5 100644 --- a/deal.II/doc/doxygen/headers/matrices.h +++ b/deal.II/doc/doxygen/headers/matrices.h @@ -143,23 +143,53 @@ * memory often also takes significant compute time. The tradeoff to avoid * excessive memory allocation cannot be avoided, however. * - * In deal.II, the CompressedSparsityPattern and CompressedSetSparsityPattern - * classes implement this "dynamic" memory scheme. They use different storage - * schemes and appear to both have use cases where they are better than the - * respective other. In either case, these classes are typically used in the - * following way: + * The following classes implement this "dynamic" memory scheme in deal.II: + * - CompressedSparsityPattern + * - CompressedSimpleSparsityPattern + * - CompressedSetSparsityPattern + * + * These classes have different performance characteristics and memory + * requirements. Which one is the "best" changes from case to case because + * it is dependent on the number of dofs, the number of couplings per dof, + * the strategy used to insert and the amount of memory available. + * + * CompressedSparsityPattern and CompressedSimpleSparsityPattern are very + * similar, where CompressedSimpleSparsityPattern trades some memory (requires + * up to twice the memory in the worst case) for additional speed which is + * noticeable in cases with many nonzero entries. CompressedSetSparsityPattern + * on the other hand uses a lot more memory but may perform better in cases with + * many nonzero entries per row. See for example the @ref step_27 "step-27" + * and @ref step_22 "step-22" tutorial programs. + * + * As a guideline you should start using CompressedSparsityPattern and try the + * other variants if you run into problems. Switching between them should be as + * simple as changing the class name because all classes have the same interface + * for adding entries. + * In either case, these classes are typically used in the following way + * (replace the class CompressedSparsityPattern with a different one from above): * @verbatim * CompressedSparsityPattern compressed_pattern (dof_handler.n_dofs()); * DoFTools::make_sparsity_pattern (dof_handler, * compressed_pattern); * constraints.condense (compressed_pattern); * - * SparsityPattern finaly_sparsity_pattern; + * SparsityPattern final_sparsity_pattern; * final_sparsity_pattern.copy_from (compressed_pattern); * @endverbatim * * The intermediate, compressed sparsity pattern is directly copied into the * "compressed" form of the final static pattern. * + *

    Dynamic block sparsity patterns

    + * + * The following classes implement an array of dynamic sparsity patterns: + * - BlockCompressedSparsityPattern + * - BlockCompressedSetSparsityPattern + * - BlockCompressedSimpleSparsityPattern + * + * These classes inherit the same tradeoffs regarding their efficiency from + * their non-block classes (see above). See their documentation and + * @ref step_22 "step-22" for more information. + * * @ingroup Matrices */ diff --git a/deal.II/lac/include/lac/block_sparsity_pattern.h b/deal.II/lac/include/lac/block_sparsity_pattern.h index 821ab634c9..2952413051 100644 --- a/deal.II/lac/include/lac/block_sparsity_pattern.h +++ b/deal.II/lac/include/lac/block_sparsity_pattern.h @@ -22,6 +22,7 @@ #include #include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -30,6 +31,7 @@ DEAL_II_NAMESPACE_OPEN template class BlockSparseMatrix; class BlockSparsityPattern; class BlockCompressedSparsityPattern; +class BlockCompressedSimpleSparsityPattern; class BlockCompressedSetSparsityPattern; /*! @addtogroup Sparsity @@ -521,6 +523,22 @@ class BlockSparsityPattern : public BlockSparsityPatternBase * of this object is lost. */ void copy_from (const BlockCompressedSetSparsityPattern &csp); + + /** + * Copy data from an object of + * type + * BlockCompressedSimpleSparsityPattern, + * i.e. resize this object to the + * size of the given argument, + * and copy over the contents of + * each of the + * subobjects. Previous content + * of this object is lost. + */ + void copy_from (const BlockCompressedSimpleSparsityPattern &csp); + + + }; @@ -532,16 +550,13 @@ class BlockSparsityPattern : public BlockSparsityPatternBase * rather acts as a @p typedef to introduce the name of this class, without * requiring the user to specify the templated name of the base class. For * information on the interface of this class refer to the base class. The - * individual blocks are based on the CompressedSparistyPattern class. + * individual blocks are based on the CompressedSparsityPattern class. * * This class is an example of the "dynamic" type of @ref Sparsity. * - * Note: This class inherits the same tradeoffs regarding its - * efficiency as CompressedSparsityPattern and as opposed to - * CompressedSetSparsity. Since the two classes (and consequently - * BlockCompressedSparsityPattern and - * BlockCompressedSetSparsityPattern) are interface-compatible it is - * frequently worth checking which version is more efficient. + * Note: + * There are several, exchangeable variations of this class, see @ref Sparsity, + * section 'Dynamic block sparsity patterns' for more information. * * Note: This class used to be called * CompressedBlockSparsityPattern. However, since it's a block wrapper around @@ -668,12 +683,9 @@ typedef BlockCompressedSparsityPattern CompressedBlockSparsityPattern; * * This class is an example of the "dynamic" type of @ref Sparsity. * - * Note: This class inherits the same tradeoffs regarding its - * efficiency as CompressedSparsityPattern and as opposed to - * CompressedSetSparsity. Since the two classes (and consequently - * BlockCompressedSparsityPattern and - * BlockCompressedSetSparsityPattern) are interface-compatible it is - * frequently worth checking which version is more efficient. + * Note: + * There are several, exchangeable variations of this class, see @ref Sparsity, + * section 'Dynamic block sparsity patterns' for more information. * * This class is used in @ref step_22 "step-22". * @@ -747,6 +759,94 @@ class BlockCompressedSetSparsityPattern : public BlockSparsityPatternBaseNote: + * There are several, exchangeable variations of this class, see @ref Sparsity, + * section 'Dynamic block sparsity patterns' for more information. + * + * @author Timo Heister, 2008 + */ +class BlockCompressedSimpleSparsityPattern : public BlockSparsityPatternBase +{ + public: + + /** + * Initialize the matrix empty, + * that is with no memory + * allocated. This is useful if + * you want such objects as + * member variables in other + * classes. You can make the + * structure usable by calling + * the reinit() function. + */ + BlockCompressedSimpleSparsityPattern (); + + /** + * Initialize the matrix with the + * given number of block rows and + * columns. The blocks themselves + * are still empty, and you have + * to call collect_sizes() after + * you assign them sizes. + */ + BlockCompressedSimpleSparsityPattern (const unsigned int n_rows, + const unsigned int n_columns); + + /** + * Initialize the pattern with + * two BlockIndices for the block + * structures of matrix rows and + * columns. This function is + * equivalent to calling the + * previous constructor with the + * length of the two index vector + * and then entering the index + * values. + */ + BlockCompressedSimpleSparsityPattern (const std::vector& row_block_sizes, + const std::vector& col_block_sizes); + + /** + * Resize the matrix to a tensor + * product of matrices with + * dimensions defined by the + * arguments. + * + * The matrix will have as many + * block rows and columns as + * there are entries in the two + * arguments. The block at + * position (i,j) will + * have the dimensions + * row_block_sizes[i] + * times col_block_sizes[j]. + */ + void reinit (const std::vector< unsigned int > &row_block_sizes, + const std::vector< unsigned int > &col_block_sizes); + + /** + * Allow the use of the reinit + * functions of the base class as + * well. + */ + using BlockSparsityPatternBase::reinit; +}; + + + /*@}*/ /*---------------------- Template functions -----------------------------------*/ diff --git a/deal.II/lac/include/lac/compressed_set_sparsity_pattern.h b/deal.II/lac/include/lac/compressed_set_sparsity_pattern.h index f1439d45b9..584f6580a4 100644 --- a/deal.II/lac/include/lac/compressed_set_sparsity_pattern.h +++ b/deal.II/lac/include/lac/compressed_set_sparsity_pattern.h @@ -82,6 +82,9 @@ template class SparseMatrix; * *

    Notes

    * + * There are several, exchangeable variations of this class, see @ref Sparsity, + * section '"Dynamic" or "compressed" sparsity patterns' for more information. + * * This class is a variation of the CompressedSparsityPattern class. * Instead of using sorted vectors together with a caching algorithm * for storing the column indices of nonzero entries, the std::set diff --git a/deal.II/lac/include/lac/compressed_simple_sparsity_pattern.h b/deal.II/lac/include/lac/compressed_simple_sparsity_pattern.h new file mode 100644 index 0000000000..86e2eb4532 --- /dev/null +++ b/deal.II/lac/include/lac/compressed_simple_sparsity_pattern.h @@ -0,0 +1,465 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 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__compressed_simple_sparsity_pattern_h +#define __deal2__compressed_simple_sparsity_pattern_h + + +#include +#include +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +template class SparseMatrix; + + +/*! @addtogroup Sparsity + *@{ + */ + + +/** + * This class acts as an intermediate form of the + * SparsityPattern class. From the interface it mostly + * represents a SparsityPattern object that is kept compressed + * at all times. However, since the final sparsity pattern is not + * known while constructing it, keeping the pattern compressed at all + * times can only be achieved at the expense of either increased + * memory or run time consumption upon use. The main purpose of this + * class is to avoid some memory bottlenecks, so we chose to implement + * it memory conservative, but the chosen data format is too unsuited + * to be used for actual matrices. It is therefore necessary to first + * copy the data of this object over to an object of type + * SparsityPattern before using it in actual matrices. + * + * Another viewpoint is that this class does not need up front allocation of a + * certain amount of memory, but grows as necessary. An extensive description + * of sparsity patterns can be found in the documentation of the @ref Sparsity + * module. + * + * This class is an example of the "dynamic" type of @ref Sparsity. + * + *

    Interface

    + * + * Since this class is intended as an intermediate replacement of the + * SparsityPattern class, it has mostly the same interface, with + * small changes where necessary. In particular, the add() + * function, and the functions inquiring properties of the sparsity + * pattern are the same. + * + * + *

    Usage

    + * + * Use this class as follows: + * @verbatim + * CompressedSimpleSparsityPattern compressed_pattern (dof_handler.n_dofs()); + * DoFTools::make_sparsity_pattern (dof_handler, + * compressed_pattern); + * constraints.condense (compressed_pattern); + * + * SparsityPattern sp; + * sp.copy_from (compressed_pattern); + * @endverbatim + * + * + *

    Notes

    + * + * There are several, exchangeable variations of this class, see @ref Sparsity, + * section '"Dynamic" or "compressed" sparsity patterns' for more information. + * + * @author Timo Heister, 2008 + */ +class CompressedSimpleSparsityPattern : public Subscriptor +{ + public: + /** + * Initialize the matrix empty, + * that is with no memory + * allocated. This is useful if + * you want such objects as + * member variables in other + * classes. You can make the + * structure usable by calling + * the reinit() function. + */ + CompressedSimpleSparsityPattern (); + + /** + * Copy constructor. This constructor is + * only allowed to be called if the + * matrix structure to be copied is + * empty. This is so in order to prevent + * involuntary copies of objects for + * temporaries, which can use large + * amounts of computing time. However, + * copy constructors are needed if yo + * want to use the STL data types on + * classes like this, e.g. to write such + * statements like v.push_back + * (CompressedSparsityPattern());, + * with @p v a vector of @p + * CompressedSparsityPattern objects. + */ + CompressedSimpleSparsityPattern (const CompressedSimpleSparsityPattern &); + + /** + * Initialize a rectangular + * matrix with @p m rows and + * @p n columns. + */ + CompressedSimpleSparsityPattern (const unsigned int m, + const unsigned int n); + + /** + * Initialize a square matrix of + * dimension @p n. + */ + CompressedSimpleSparsityPattern (const unsigned int n); + + /** + * Copy operator. For this the + * same holds as for the copy + * constructor: it is declared, + * defined and fine to be called, + * but the latter only for empty + * objects. + */ + CompressedSimpleSparsityPattern & operator = (const CompressedSimpleSparsityPattern &); + + /** + * Reallocate memory and set up + * data structures for a new + * matrix with @p m rows and + * @p n columns, with at most + * max_entries_per_row() nonzero + * entries per row. + */ + void reinit (const unsigned int m, + const unsigned int n); + + /** + * Since this object is kept + * compressed at all times anway, + * this function does nothing, + * but is declared to make the + * interface of this class as + * much alike as that of the + * SparsityPattern class. + */ + void compress (); + + /** + * Return whether the object is + * empty. It is empty if no + * memory is allocated, which is + * the same as that both + * dimensions are zero. + */ + bool empty () const; + + /** + * Return the maximum number of + * entries per row. Note that + * this number may change as + * entries are added. + */ + unsigned int max_entries_per_row () const; + + /** + * Add a nonzero entry to the + * matrix. If the entry already + * exists, nothing bad happens. + */ + void add (const unsigned int i, + const unsigned int j); + + /** + * Check if a value at a certain + * position may be non-zero. + */ + bool exists (const unsigned int i, + const unsigned int j) const; + + /** + * Make the sparsity pattern + * symmetric by adding the + * sparsity pattern of the + * transpose object. + * + * This function throws an + * exception if the sparsity + * pattern does not represent a + * square matrix. + */ + void symmetrize (); + + /** + * Print the sparsity of the + * matrix. The output consists of + * one line per row of the format + * [i,j1,j2,j3,...]. i + * is the row number and + * jn are the allocated + * columns in this row. + */ + void print (std::ostream &out) const; + + /** + * Print the sparsity of the matrix in a + * format that @p gnuplot understands and + * which can be used to plot the sparsity + * pattern in a graphical way. The format + * consists of pairs i j of + * nonzero elements, each representing + * one entry of this matrix, one per line + * of the output file. Indices are + * counted from zero on, as usual. Since + * sparsity patterns are printed in the + * same way as matrices are displayed, we + * print the negative of the column + * index, which means that the + * (0,0) element is in the top + * left rather than in the bottom left + * corner. + * + * Print the sparsity pattern in + * gnuplot by setting the data style + * to dots or points and use the + * @p plot command. + */ + void print_gnuplot (std::ostream &out) const; + + /** + * Return number of rows of this + * matrix, which equals the dimension + * of the image space. + */ + unsigned int n_rows () const; + + /** + * Return number of columns of this + * matrix, which equals the dimension + * of the range space. + */ + unsigned int n_cols () const; + + /** + * Number of entries in a specific row. + */ + unsigned int row_length (const unsigned int row) const; + + /** + * Access to column number field. + * Return the column number of + * the @p indexth entry in @p row. + */ + unsigned int column_number (const unsigned int row, + const unsigned int index) const; + + /** + * Compute the bandwidth of the matrix + * represented by this structure. The + * bandwidth is the maximum of + * $|i-j|$ for which the index pair + * $(i,j)$ represents a nonzero entry + * of the matrix. + */ + unsigned int bandwidth () const; + + /** + * Return the number of nonzero elements + * allocated through this sparsity pattern. + */ + unsigned int n_nonzero_elements () const; + + /** + * Return whether this object stores only + * those entries that have been added + * explicitly, or if the sparsity pattern + * contains elements that have been added + * through other means (implicitly) while + * building it. For the current class, + * the result is always true. + * + * This function mainly serves the + * purpose of describing the current + * class in cases where several kinds of + * sparsity patterns can be passed as + * template arguments. + */ + static + bool stores_only_added_elements (); + + + + private: + /** + * Number of rows that this sparsity + * structure shall represent. + */ + unsigned int rows; + + /** + * Number of columns that this sparsity + * structure shall represent. + */ + unsigned int cols; + + /** + * Store some data for each row + * describing which entries of this row + * are nonzero. Data is stored sorted in + * the @p entries std::vector. + * The vector per row is dynamically + * growing upon insertion doubling its + * memory each time. + */ + struct Line + { + public: + /** + * Storage for the column indices of + * this row. This array is always + * kept sorted. + */ + std::vector entries; + + /** + * Constructor. + */ + Line (); + + /** + * Add the given column number to + * this line. + */ + void add (const unsigned int col_num); + + }; + + + /** + * Actual data: store for each + * row the set of nonzero + * entries. + */ + std::vector lines; +}; + +/*@}*/ +/*---------------------- Inline functions -----------------------------------*/ + + +inline +void +CompressedSimpleSparsityPattern::Line::add (const unsigned int j) +{ + /* first check the last element (or if line is still empty) */ + if ( (entries.size()==0) || ( *(--entries.end()) < j) ) + { + entries.push_back(j); + return; + } + + /* do a binary search to find the place where to insert: */ + std::vector::iterator it = std::lower_bound(entries.begin(), + entries.end(), + j); + + /* If this entry is a duplicate, exit immediately */ + if (*it == j) + return; + + /* Insert at the right place in the vector. + Vector grows automatically to fit + elements. Always doubles its size.*/ + entries.insert(it, j); +} + + + +inline +unsigned int +CompressedSimpleSparsityPattern::n_rows () const +{ + return rows; +} + + + +inline +unsigned int +CompressedSimpleSparsityPattern::n_cols () const +{ + return cols; +} + + + +inline +void +CompressedSimpleSparsityPattern::add (const unsigned int i, + const unsigned int j) +{ + Assert (i class SparseMatrix; * *

    Notes

    * - * A variation of this class is the CompressedSetSparsityPattern class - * that appears to be more efficient in some situations, in particular - * when using hp finite elements or, more generally, in cases with - * many nonzero entries. See for example the @ref step_27 "step-27" - * and @ref step_22 "step-22" tutorial programs. + * There are several, exchangeable variations of this class, see @ref Sparsity, + * section '"Dynamic" or "compressed" sparsity patterns' for more information. * * @author Wolfgang Bangerth, 2001 */ diff --git a/deal.II/lac/include/lac/sparsity_pattern.h b/deal.II/lac/include/lac/sparsity_pattern.h index dc226b1783..6b29135671 100644 --- a/deal.II/lac/include/lac/sparsity_pattern.h +++ b/deal.II/lac/include/lac/sparsity_pattern.h @@ -31,6 +31,7 @@ template class VectorSlice; class CompressedSparsityPattern; class CompressedSetSparsityPattern; +class CompressedSimpleSparsityPattern; @@ -931,16 +932,28 @@ class SparsityPattern : public Subscriptor const bool optimize_diagonal = true); + /** + * Copy data from an object of + * type + * CompressedSetSparsityPattern. + * Previous content of this + * object is lost, and the + * sparsity pattern is in + * compressed mode afterwards. + */ + void copy_from (const CompressedSetSparsityPattern &csp, + const bool optimize_diagonal = true); + /** * Copy data from an object of * type - * CompressedSetSparsityPattern. + * CompressedSimpleSparsityPattern. * Previous content of this * object is lost, and the * sparsity pattern is in * compressed mode afterwards. */ - void copy_from (const CompressedSetSparsityPattern &csp, + void copy_from (const CompressedSimpleSparsityPattern &csp, const bool optimize_diagonal = true); diff --git a/deal.II/lac/source/block_sparsity_pattern.cc b/deal.II/lac/source/block_sparsity_pattern.cc index 1789835f05..6f6ea0be85 100644 --- a/deal.II/lac/source/block_sparsity_pattern.cc +++ b/deal.II/lac/source/block_sparsity_pattern.cc @@ -307,6 +307,7 @@ BlockSparsityPatternBase::print_gnuplot(std::ostream &out) template class BlockSparsityPatternBase; template class BlockSparsityPatternBase; +template class BlockSparsityPatternBase; template class BlockSparsityPatternBase; @@ -396,7 +397,22 @@ BlockSparsityPattern::copy_from (const BlockCompressedSparsityPattern &csp) collect_sizes(); } +void +BlockSparsityPattern::copy_from (const BlockCompressedSimpleSparsityPattern &csp) +{ + // delete old content, set block + // sizes anew + reinit (csp.n_block_rows(), csp.n_block_cols()); + // copy over blocks + for (unsigned int i=0; icollect_sizes(); } +BlockCompressedSimpleSparsityPattern::BlockCompressedSimpleSparsityPattern () +{} + + + +BlockCompressedSimpleSparsityPattern:: +BlockCompressedSimpleSparsityPattern (const unsigned int n_rows, + const unsigned int n_columns) + : + BlockSparsityPatternBase(n_rows, + n_columns) +{} + + +BlockCompressedSimpleSparsityPattern:: +BlockCompressedSimpleSparsityPattern (const std::vector& row_indices, + const std::vector& col_indices) + : + BlockSparsityPatternBase(row_indices.size(), + col_indices.size()) +{ + for (unsigned int i=0;iblock(i,j).reinit(row_indices[i],col_indices[j]); + this->collect_sizes(); +} + + +void +BlockCompressedSimpleSparsityPattern::reinit ( + const std::vector< unsigned int > &row_block_sizes, + const std::vector< unsigned int > &col_block_sizes) +{ + BlockSparsityPatternBase::reinit(row_block_sizes.size(), col_block_sizes.size()); + for (unsigned int i=0;iblock(i,j).reinit(row_block_sizes[i],col_block_sizes[j]); + this->collect_sizes(); +} + BlockCompressedSetSparsityPattern::BlockCompressedSetSparsityPattern () diff --git a/deal.II/lac/source/compressed_simple_sparsity_pattern.cc b/deal.II/lac/source/compressed_simple_sparsity_pattern.cc new file mode 100644 index 0000000000..4a3c702434 --- /dev/null +++ b/deal.II/lac/source/compressed_simple_sparsity_pattern.cc @@ -0,0 +1,241 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006 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 +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +CompressedSimpleSparsityPattern::CompressedSimpleSparsityPattern () + : + rows(0), + cols(0) +{} + + + +CompressedSimpleSparsityPattern:: +CompressedSimpleSparsityPattern (const CompressedSimpleSparsityPattern &s) + : + Subscriptor(), + rows(0), + cols(0) +{ + Assert (s.rows == 0, ExcInvalidConstructorCall()); + Assert (s.cols == 0, ExcInvalidConstructorCall()); +} + + + +CompressedSimpleSparsityPattern::CompressedSimpleSparsityPattern (const unsigned int m, + const unsigned int n) + : + rows(0), + cols(0) +{ + reinit (m,n); +} + + + +CompressedSimpleSparsityPattern::CompressedSimpleSparsityPattern (const unsigned int n) + : + rows(0), + cols(0) +{ + reinit (n,n); +} + + + +CompressedSimpleSparsityPattern & +CompressedSimpleSparsityPattern::operator = (const CompressedSimpleSparsityPattern &s) +{ + Assert (s.rows == 0, ExcInvalidConstructorCall()); + Assert (s.cols == 0, ExcInvalidConstructorCall()); + + Assert (rows == 0, ExcInvalidConstructorCall()); + Assert (cols == 0, ExcInvalidConstructorCall()); + + return *this; +} + + + +void +CompressedSimpleSparsityPattern::reinit (const unsigned int m, + const unsigned int n) +{ + rows = m; + cols = n; + + std::vector new_lines (rows); + lines.swap (new_lines); +} + + + +void +CompressedSimpleSparsityPattern::compress () +{} + + + +bool +CompressedSimpleSparsityPattern::empty () const +{ + return ((rows==0) && (cols==0)); +} + + + +unsigned int +CompressedSimpleSparsityPattern::max_entries_per_row () const +{ + unsigned int m = 0; + for (unsigned int i=0; i(lines[i].entries.size())); + } + + return m; +} + + + +bool +CompressedSimpleSparsityPattern::exists (const unsigned int i, + const unsigned int j) const +{ + Assert (i::const_iterator + j=lines[row].entries.begin(); + j != lines[row].entries.end(); + ++j) + // add the transpose entry if + // this is not the diagonal + if (row != *j) + add (*j, row); + } +} + + + +void +CompressedSimpleSparsityPattern::print (std::ostream &out) const +{ + for (unsigned int row=0; row::const_iterator + j=lines[row].entries.begin(); + j != lines[row].entries.end(); ++j) + out << ',' << *j; + + out << ']' << std::endl; + } + + AssertThrow (out, ExcIO()); +} + + + +void +CompressedSimpleSparsityPattern::print_gnuplot (std::ostream &out) const +{ + for (unsigned int row=0; row::const_iterator + j=lines[row].entries.begin(); + j != lines[row].entries.end(); ++j) + // while matrix entries are usually + // written (i,j), with i vertical and + // j horizontal, gnuplot output is + // x-y, that is we have to exchange + // the order of output + out << *j << " " << -static_cast(row) << std::endl; + } + + + AssertThrow (out, ExcIO()); +} + + + +unsigned int +CompressedSimpleSparsityPattern::bandwidth () const +{ + unsigned int b=0; + for (unsigned int row=0; row::const_iterator + j=lines[row].entries.begin(); + j != lines[row].entries.end(); ++j) + if (static_cast(std::abs(static_cast(row-*j))) > b) + b = std::abs(static_cast(row-*j)); + } + + return b; +} + + + +unsigned int +CompressedSimpleSparsityPattern::n_nonzero_elements () const +{ + unsigned int n=0; + for (unsigned int i=0; i #include #include +#include #include #include @@ -646,6 +647,50 @@ SparsityPattern::copy_from (const CompressedSetSparsityPattern &csp, compress (); } +void +SparsityPattern::copy_from (const CompressedSimpleSparsityPattern &csp, + const bool optimize_diag) +{ + // first determine row lengths for + // each row. if the matrix is + // quadratic, then we might have to + // add an additional entry for the + // diagonal, if that is not yet + // present. as we have to call + // compress anyway later on, don't + // bother to check whether that + // diagonal entry is in a certain + // row or not + const bool is_square = optimize_diag && (csp.n_rows() == csp.n_cols()); + std::vector row_lengths (csp.n_rows()); + for (unsigned int i=0; i