From: Wolfgang Bangerth Date: Mon, 21 Jun 2004 17:02:47 +0000 (+0000) Subject: Use a more sophisticated algorithm for the compressed sparsity pattern that uses... X-Git-Tag: v8.0.0~15019 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=dd61fabc1b8f0a2f092a746b7a5f212f8ed499ed;p=dealii.git Use a more sophisticated algorithm for the compressed sparsity pattern that uses caches and std::vector instead of std::set. std::set allocates 20 bytes each time we add a single integer, which can't be efficient both from memory and run-time viewpoints. git-svn-id: https://svn.dealii.org/trunk@9440 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/compressed_sparsity_pattern.h b/deal.II/lac/include/lac/compressed_sparsity_pattern.h index 4e675de5a5..3ad0248d9e 100644 --- a/deal.II/lac/include/lac/compressed_sparsity_pattern.h +++ b/deal.II/lac/include/lac/compressed_sparsity_pattern.h @@ -347,12 +347,93 @@ class CompressedSparsityPattern : public Subscriptor */ 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 @p 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 @p flush_cache function is + * marked const, and the data members are + * marked @p mutable. + */ + struct Line + { + /** + * Storage for the column indices of + * this row, unless they are still in + * the cache. This array is always + * kept sorted. + */ + mutable std::vector entries; + + /** + * Size of the cache. + */ + static const unsigned int cache_size = 8; + + /** + * Cache of entries that have not yet + * been written to @p entries; + */ + mutable unsigned int cache[cache_size]; + + /** + * Number of entries in the cache. + */ + mutable unsigned int cache_entries; + + /** + * Constructor. + */ + Line (); + + /** + * Add the given column number to + * this line. + */ + void add (const unsigned int col_num); + + /** + * Flush the cache my merging it with + * the @p entries array. + */ + void flush_cache () const; + }; + + /** * Actual data: store for each * row the set of nonzero * entries. */ - std::vector > column_indices; + std::vector lines; }; /*@}*/ @@ -367,6 +448,7 @@ CompressedSparsityPattern::n_rows () const } + inline unsigned int CompressedSparsityPattern::n_cols () const @@ -374,4 +456,226 @@ CompressedSparsityPattern::n_cols () const return cols; } + + +inline +void +CompressedSparsityPattern::add (const unsigned int i, + const unsigned int j) +{ + Assert (i new_entries; + new_entries.reserve (n_new_entries); + unsigned int cache_position = 0; + unsigned int entry_position = 0; + while ((entry_position #include + +const unsigned int CompressedSparsityPattern::Line::cache_size; + + + CompressedSparsityPattern::CompressedSparsityPattern () : rows(0), @@ -84,8 +89,8 @@ CompressedSparsityPattern::reinit (const unsigned int m, rows = m; cols = n; - std::vector > new_column_indices (rows); - column_indices.swap (new_column_indices); + std::vector new_lines (rows); + lines.swap (new_lines); } @@ -108,35 +113,28 @@ unsigned int CompressedSparsityPattern::max_entries_per_row () const { unsigned int m = 0; - for (unsigned int i=1; i(column_indices[i].size())); - + for (unsigned int i=0; i(lines[i].entries.size())); + } + return m; } -void -CompressedSparsityPattern::add (const unsigned int i, - const unsigned int j) -{ - Assert (i::const_iterator i=column_indices[row].begin(); - i!=column_indices[row].end(); ++i) + { + lines[row].flush_cache (); + for (std::vector::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 != *i) - add (*i, row); + if (row != *j) + add (*j, row); + } } void CompressedSparsityPattern::print_gnuplot (std::ostream &out) const -{ +{ for (unsigned int row=0; row::const_iterator i=column_indices[row].begin(); - i!=column_indices[row].end(); ++i) - // 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 << *i << " " << -static_cast(row) << std::endl; + { + lines[row].flush_cache (); + for (std::vector::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 -CompressedSparsityPattern::row_length (const unsigned int row) const -{ - return column_indices[row].size(); -} - - - -unsigned int -CompressedSparsityPattern::column_number (const unsigned int row, - const unsigned int index) const -{ - Assert (index < column_indices[row].size(), - ExcIndexRange (index, 0, column_indices[row].size())); - std::set::const_iterator p = column_indices[row].begin(); - std::advance (p, index); - return *p; -} - - - unsigned int CompressedSparsityPattern::bandwidth () const { unsigned int b=0; for (unsigned int row=0; row::const_iterator i=column_indices[row].begin(); - i!=column_indices[row].end(); ++i) - if (static_cast(std::abs(static_cast(row-*i))) > b) - b = std::abs(static_cast(row-*i)); - + { + lines[row].flush_cache (); + + for (std::vector::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; } @@ -229,6 +219,10 @@ CompressedSparsityPattern::n_nonzero_elements () const { unsigned int n=0; for (unsigned int i=0; i