From 776167a22c92d7185bf5f4631b509cfbb91b97c7 Mon Sep 17 00:00:00 2001 From: kayser-herold Date: Wed, 6 Jun 2007 15:22:40 +0000 Subject: [PATCH] Added an alternative implementation of a compressed sparsity pattern, which is based on the std::set class. Due to unknown reasons, this class seems to perform better than the original CSP in the context of hp-adaptivity. git-svn-id: https://svn.dealii.org/trunk@14758 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/compressed_set_sparsity_pattern.h | 510 ++++++++++++++++++ deal.II/lac/include/lac/sparsity_pattern.h | 15 + .../source/compressed_set_sparsity_pattern.cc | 220 ++++++++ deal.II/lac/source/sparsity_pattern.cc | 52 ++ 4 files changed, 797 insertions(+) create mode 100644 deal.II/lac/include/lac/compressed_set_sparsity_pattern.h create mode 100644 deal.II/lac/source/compressed_set_sparsity_pattern.cc diff --git a/deal.II/lac/include/lac/compressed_set_sparsity_pattern.h b/deal.II/lac/include/lac/compressed_set_sparsity_pattern.h new file mode 100644 index 0000000000..63cb2fbaad --- /dev/null +++ b/deal.II/lac/include/lac/compressed_set_sparsity_pattern.h @@ -0,0 +1,510 @@ +//--------------------------------------------------------------------------- +// $Id: compressed_sparsity_pattern.h 14038 2006-10-23 02:46:34Z bangerth $ +// 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. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__compressed_set_sparsity_pattern_h +#define __deal2__compressed_set_sparsity_pattern_h + + +#include +#include +#include + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +template class SparseMatrix; + +/*! @addtogroup Matrix1 + *@{ + */ + + +/** + * 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. + * + * + *

Rationale

+ * + * When constructing the sparsity pattern of a matrix, you usually + * first have to provide an empty sparsity pattern object with a fixed + * maximal number of entries per row. To find out about this maximal + * row length, one usually calls the function + * DoFHandler::max_couplings_between_dofs() which returns an estimate for + * that quantity or DoFTools::compute_row_length_vector(), which gives + * a better estimate for each row. While this estimate is usually + * quite good in 2d and exact in 1d, it is often significantly too + * large in 3d and especially for higher order elements. Furthermore, + * normally only a small fraction of the rows of a matrix will end up + * having the maximal number of nonzero entries per row (usually those + * nodes adjacent to hanging nodes), most have much less. In effect, + * the empty SparsityPattern object has allocated much too much + * memory. Although this unnecessarily allocated memory is later freed + * when SparsityPattern::compress() is called, this overallocation + * has, with higher order elements and in 3d, sometimes been so large + * that the program aborted due to lack of memory. + * + * This class therefore provides an alternative representation of a + * sparsity pattern: we don't specify a maximal row length initially, + * but store a set of column indices indicating possible nonzero + * entries in the sparsity pattern for each row. This is very much + * like the final "compressed" format used in the + * SparsityPattern object after compression, but uses a less + * compact memory storage format, since the exact number of entries + * per row is only known a posteriori and since it may change (for the + * SparsityPattern class, no more changes are allowed after + * compressing it). We can therefore not store all the column indices + * in a big array, but have to use a vector of sets. This can later be + * used to actually initialize a SparsityPattern object with the + * then final set of necessary indices. + * + * + *

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. + * + *

Notes

+ * + * 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 NZ entries, the std::set + * container is used. This solution might not be the fastest in + * all situations, but seems to work much better than the + * CompressedSparsityPattern in the context of hp-adaptivity. + * On the other hand, a benchmark where NZ entries were randomly inserted + * into the sparsity pattern revealed that this class is slower + * by a factor 4-6 in this situation. Hence, currently the suggestion + * is to carefully analyse which of the CompressedSparsityPattern + * classes works best in a certain setting. An algorithm which + * performs equally well in all situations still has to be found. + * + *

Usage

+ * + * Use this class as follows: + * @verbatim + * CompressedSetSparsityPattern 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 + * + * + * @author Wolfgang Bangerth, 2001 + */ +class CompressedSetSparsityPattern : public Subscriptor +{ + public: + typedef std::set::iterator CSSPIterator; + + + /** + * 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. + */ + CompressedSetSparsityPattern (); + + /** + * 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 + * (CompressedSetSparsityPattern());, + * with @p v a vector of @p + * CompressedSetSparsityPattern objects. + * + * Usually, it is sufficient to + * use the explicit keyword to + * disallow unwanted temporaries, + * but for the STL vectors, this + * does not work. Since copying a + * structure like this is not + * useful anyway because multiple + * matrices can use the same + * sparsity structure, copies are + * only allowed for empty + * objects, as described above. + */ + CompressedSetSparsityPattern (const CompressedSetSparsityPattern &); + + /** + * Initialize a rectangular + * matrix with @p m rows and + * @p n columns. + */ + CompressedSetSparsityPattern (const unsigned int m, + const unsigned int n); + + /** + * Initialize a square matrix of + * dimension @p n. + */ + CompressedSetSparsityPattern (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. + */ + CompressedSetSparsityPattern & operator = (const CompressedSetSparsityPattern &); + + /** + * 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 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. + */ + CSSPIterator row_begin (const unsigned int row) const; + CSSPIterator row_end (const unsigned int row) 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 + * of this matrix. Actually, it returns + * the number of entries in the sparsity + * pattern; if any of the entries should + * happen to be zero, it is counted + * anyway. + * + * This function may only be called if + * the matrix struct is compressed. It + * does not make too much sense otherwise + * anyway. + */ + unsigned int n_nonzero_elements () const; + + 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 organized as + * follows: if an entry is added to a + * row, it is first added to the #cache + * variable, irrespective of whether an + * entry with same column number has + * already been added. Only if the cache + * is full do we flush it by removing + * duplicates, removing entries that are + * already stored in the @p entries + * array, sorting everything, and merging + * the two arrays. + * + * The reasoning behind this scheme is + * that memory allocation is expensive, + * and we only want to do it when really + * necessary. Previously (in deal.II + * versions up to 5.0), we used to store + * the column indices inside a std::set, + * but this would allocate 20 bytes each + * time we added an entry. Using the + * present scheme, we only need to + * allocate memory once for every 8 added + * entries, and we waste a lot less + * memory by not using a balanced tree + * for storing column indices. + * + * Since some functions that are @p const + * need to access the data of this + * object, but need to flush caches + * before, the flush_cache() function is + * marked const, and the data members are + * marked @p mutable. + * + * A small testseries about the size of + * the cache showed that the run time of + * a small program just testing the + * compressed sparsity pattern element + * insertion routine ran for 3.6 seconds + * with a cache size of 8, and 4.2 + * seconds with a cache size of 16. We + * deem even smaller cache sizes + * undesirable, since they lead to more + * memory allocations, while larger cache + * sizes lead to waste of memory. The + * original version of this class, with + * one std::set per row took 8.2 seconds + * on the same program. + */ + struct Line + { + std::set 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 +CompressedSetSparsityPattern::Line::add (const unsigned int j) +{ + entries.insert (j); +} + + + +inline +unsigned int +CompressedSetSparsityPattern::n_rows () const +{ + return rows; +} + + + +inline +unsigned int +CompressedSetSparsityPattern::n_cols () const +{ + return cols; +} + + + +inline +void +CompressedSetSparsityPattern::add (const unsigned int i, + const unsigned int j) +{ + Assert (i class SparseMatrix; template class VectorSlice; class CompressedSparsityPattern; +class CompressedSetSparsityPattern; @@ -968,6 +969,20 @@ class SparsityPattern : public Subscriptor void copy_from (const CompressedSparsityPattern &csp, 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); + + /** * Take a full matrix and use its * nonzero entries to generate a diff --git a/deal.II/lac/source/compressed_set_sparsity_pattern.cc b/deal.II/lac/source/compressed_set_sparsity_pattern.cc new file mode 100644 index 0000000000..fd50f31fca --- /dev/null +++ b/deal.II/lac/source/compressed_set_sparsity_pattern.cc @@ -0,0 +1,220 @@ +//--------------------------------------------------------------------------- +// $Id: compressed_sparsity_pattern.cc 14038 2006-10-23 02:46:34Z bangerth $ +// 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 + + + +CompressedSetSparsityPattern::CompressedSetSparsityPattern () + : + rows(0), + cols(0) +{} + + + +CompressedSetSparsityPattern:: +CompressedSetSparsityPattern (const CompressedSetSparsityPattern &s) + : + Subscriptor(), + rows(0), + cols(0) +{ + Assert (s.rows == 0, ExcInvalidConstructorCall()); + Assert (s.cols == 0, ExcInvalidConstructorCall()); +} + + + +CompressedSetSparsityPattern::CompressedSetSparsityPattern (const unsigned int m, + const unsigned int n) + : + rows(0), + cols(0) +{ + reinit (m,n); +} + + + +CompressedSetSparsityPattern::CompressedSetSparsityPattern (const unsigned int n) + : + rows(0), + cols(0) +{ + reinit (n,n); +} + + + +CompressedSetSparsityPattern & +CompressedSetSparsityPattern::operator = (const CompressedSetSparsityPattern &s) +{ + Assert (s.rows == 0, ExcInvalidConstructorCall()); + Assert (s.cols == 0, ExcInvalidConstructorCall()); + + Assert (rows == 0, ExcInvalidConstructorCall()); + Assert (cols == 0, ExcInvalidConstructorCall()); + + return *this; +} + + + +void +CompressedSetSparsityPattern::reinit (const unsigned int m, + const unsigned int n) +{ + rows = m; + cols = n; + + std::vector new_lines (rows); + lines.swap (new_lines); +} + + + +void +CompressedSetSparsityPattern::compress () +{} + + + +bool +CompressedSetSparsityPattern::empty () const +{ + return ((rows==0) && (cols==0)); +} + + + +unsigned int +CompressedSetSparsityPattern::max_entries_per_row () const +{ + unsigned int m = 0; + for (unsigned int i=0; i(lines[i].entries.size())); + } + + return m; +} + + + +bool +CompressedSetSparsityPattern::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 +CompressedSetSparsityPattern::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 +CompressedSetSparsityPattern::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 +CompressedSetSparsityPattern::n_nonzero_elements () const +{ + unsigned int n=0; + for (unsigned int i=0; i #include #include +#include #include #include @@ -597,6 +598,57 @@ SparsityPattern::copy_from (const CompressedSparsityPattern &csp, } + +void +SparsityPattern::copy_from (const CompressedSetSparsityPattern &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 void SparsityPattern::copy_from (const FullMatrix &matrix, const bool optimize_diag) -- 2.39.5