From 19c8bfedc2b956315a19211a91ec07eda6e45c97 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Tue, 8 Dec 2009 16:02:22 +0000 Subject: [PATCH] Implement an interface row_begin and row_end that all the non-blocked sparsity patterns share. This simplifies some copy_from functions, where we now only need one implementation and no specializations. Also use an empty IndexSet for CompressedSimpleSparsity when we work in serial. git-svn-id: https://svn.dealii.org/trunk@20212 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/chunk_sparsity_pattern.h | 34 +-- .../lac/compressed_set_sparsity_pattern.h | 2 + .../lac/compressed_simple_sparsity_pattern.h | 71 ++++++- .../include/lac/compressed_sparsity_pattern.h | 72 +++++-- .../lac/include/lac/sparse_matrix.templates.h | 2 +- deal.II/lac/include/lac/sparsity_pattern.h | 198 ++++++++++-------- deal.II/lac/source/chunk_sparsity_pattern.cc | 107 ++++------ .../compressed_simple_sparsity_pattern.cc | 59 +++--- deal.II/lac/source/sparsity_pattern.cc | 130 ++---------- .../lac/source/trilinos_sparsity_pattern.cc | 173 ++------------- 10 files changed, 353 insertions(+), 495 deletions(-) diff --git a/deal.II/lac/include/lac/chunk_sparsity_pattern.h b/deal.II/lac/include/lac/chunk_sparsity_pattern.h index d2e3252a7c..3defa0f90f 100644 --- a/deal.II/lac/include/lac/chunk_sparsity_pattern.h +++ b/deal.II/lac/include/lac/chunk_sparsity_pattern.h @@ -442,32 +442,18 @@ class ChunkSparsityPattern : public Subscriptor const bool optimize_diagonal = true); /** - * Copy data from an object of - * type - * CompressedSparsityPattern. - * Previous content of this - * object is lost, and the - * sparsity pattern is in - * compressed mode afterwards. - */ - void copy_from (const CompressedSparsityPattern &csp, - const unsigned int chunk_size, - 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 + * Copy data from an object of type + * CompressedSparsityPattern, + * CompressedSetSparsityPattern or + * CompressedSimpleSparsityPattern. + * Previous content of this object is + * lost, and the sparsity pattern is in * compressed mode afterwards. */ - void copy_from (const CompressedSetSparsityPattern &csp, - const unsigned int chunk_size, - const bool optimize_diagonal = true); - + template + void copy_from (const SparsityType &csp, + const unsigned int chunk_size, + const bool optimize_diagonal = true); /** * Take a full matrix and use its 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 946526c24e..a5991ef58c 100644 --- a/deal.II/lac/include/lac/compressed_set_sparsity_pattern.h +++ b/deal.II/lac/include/lac/compressed_set_sparsity_pattern.h @@ -492,6 +492,7 @@ CompressedSetSparsityPattern::row_length (const unsigned int row) const } + inline CompressedSetSparsityPattern::row_iterator CompressedSetSparsityPattern::row_begin (const unsigned int row) const @@ -500,6 +501,7 @@ CompressedSetSparsityPattern::row_begin (const unsigned int row) const } + inline CompressedSetSparsityPattern::row_iterator CompressedSetSparsityPattern::row_end (const unsigned int row) const diff --git a/deal.II/lac/include/lac/compressed_simple_sparsity_pattern.h b/deal.II/lac/include/lac/compressed_simple_sparsity_pattern.h index 677fdf5ec0..e3bde0a2a3 100644 --- a/deal.II/lac/include/lac/compressed_simple_sparsity_pattern.h +++ b/deal.II/lac/include/lac/compressed_simple_sparsity_pattern.h @@ -87,6 +87,14 @@ template class SparseMatrix; class CompressedSimpleSparsityPattern : public Subscriptor { public: + /** + * An iterator that can be used to + * iterate over the elements of a single + * row. The result of dereferencing such + * an iterator is a column index. + */ + typedef std::vector::const_iterator row_iterator; + /** * Initialize the matrix empty, * that is with no memory @@ -302,6 +310,18 @@ class CompressedSimpleSparsityPattern : public Subscriptor unsigned int column_number (const unsigned int row, const unsigned int index) const; + /** + * Return an iterator that can loop over + * all entries in the given + * row. Dereferencing the iterator yields + * a column index. + */ + row_iterator row_begin (const unsigned int row) const; + + /** + * Returns the end of the current row. + */ + row_iterator row_end (const unsigned int row) const; /** * Compute the bandwidth of the matrix * represented by this structure. The @@ -487,10 +507,12 @@ CompressedSimpleSparsityPattern::add (const unsigned int i, Assert (i 0 && !rowset.is_element(i)) return; - lines[rowset.index_within_set(i)].add (j); + const unsigned int rowindex = + rowset.size()==0 ? i : rowset.nth_index_in_set(i); + lines[rowindex].add (j); } @@ -505,10 +527,12 @@ CompressedSimpleSparsityPattern::add_entries (const unsigned int row, { Assert (row < rows, ExcIndexRange (row, 0, rows)); - if (!rowset.is_element(row)) + if (rowset.size() > 0 && !rowset.is_element(row)) return; - lines[rowset.index_within_set(row)].add_entries (begin, end, indices_are_sorted); + const unsigned int rowindex = + rowset.size()==0 ? row : rowset.nth_index_in_set(row); + lines[rowindex].add_entries (begin, end, indices_are_sorted); } @@ -524,9 +548,12 @@ unsigned int CompressedSimpleSparsityPattern::row_length (const unsigned int row) const { Assert (row < n_rows(), ExcIndexRange (row, 0, n_rows())); - Assert( rowset.is_element(row), ExcInternalError()); + if (rowset.size() > 0 && !rowset.is_element(row)) + return 0; - return lines[rowset.index_within_set(row)].entries.size(); + const unsigned int rowindex = + rowset.size()==0 ? row : rowset.nth_index_in_set(row); + return lines[rowindex].entries.size(); } @@ -537,14 +564,38 @@ CompressedSimpleSparsityPattern::column_number (const unsigned int row, const unsigned int index) const { Assert (row < n_rows(), ExcIndexRange (row, 0, n_rows())); - Assert (index < lines[rowset.index_within_set(row)].entries.size(), - ExcIndexRange (index, 0, lines[rowset.index_within_set(row)].entries.size())); - Assert( rowset.is_element(row), ExcInternalError()); + Assert( rowset.size() == 0 || rowset.is_element(row), ExcInternalError()); + + const unsigned int local_row = rowset.size() ? rowset.index_within_set(row) : row; + Assert (index < lines[local_row].entries.size(), + ExcIndexRange (index, 0, lines[local_row].entries.size())); + return lines[local_row].entries[index]; +} + + + +inline +CompressedSimpleSparsityPattern::row_iterator +CompressedSimpleSparsityPattern::row_begin (const unsigned int row) const +{ + Assert (row < n_rows(), ExcIndexRange (row, 0, n_rows())); + const unsigned int local_row = rowset.size() ? rowset.index_within_set(row) : row; + return lines[local_row].entries.begin(); +} + - return lines[rowset.index_within_set(row)].entries[index]; + +inline +CompressedSimpleSparsityPattern::row_iterator +CompressedSimpleSparsityPattern::row_end (const unsigned int row) const +{ + Assert (row < n_rows(), ExcIndexRange (row, 0, n_rows())); + const unsigned int local_row = rowset.size() ? rowset.index_within_set(row) : row; + return lines[local_row].entries.end(); } + inline const IndexSet & CompressedSimpleSparsityPattern::row_index_set () const diff --git a/deal.II/lac/include/lac/compressed_sparsity_pattern.h b/deal.II/lac/include/lac/compressed_sparsity_pattern.h index 4603b5fc95..d910a8352f 100644 --- a/deal.II/lac/include/lac/compressed_sparsity_pattern.h +++ b/deal.II/lac/include/lac/compressed_sparsity_pattern.h @@ -89,6 +89,14 @@ template class SparseMatrix; class CompressedSparsityPattern : public Subscriptor { public: + /** + * An iterator that can be used to + * iterate over the elements of a single + * row. The result of dereferencing such + * an iterator is a column index. + */ + typedef std::vector::const_iterator row_iterator; + /** * Initialize the matrix empty, * that is with no memory @@ -100,7 +108,7 @@ class CompressedSparsityPattern : public Subscriptor * the reinit() function. */ CompressedSparsityPattern (); - + /** * Copy constructor. This constructor is * only allowed to be called if the @@ -126,7 +134,7 @@ class CompressedSparsityPattern : public Subscriptor */ CompressedSparsityPattern (const unsigned int m, const unsigned int n); - + /** * Initialize a square matrix of * dimension @p n. @@ -142,7 +150,7 @@ class CompressedSparsityPattern : public Subscriptor * objects. */ CompressedSparsityPattern & operator = (const CompressedSparsityPattern &); - + /** * Reallocate memory and set up * data structures for a new @@ -153,7 +161,7 @@ class CompressedSparsityPattern : public Subscriptor */ void reinit (const unsigned int m, const unsigned int n); - + /** * Since this object is kept * compressed at all times anway, @@ -164,7 +172,7 @@ class CompressedSparsityPattern : public Subscriptor * SparsityPattern class. */ void compress (); - + /** * Return whether the object is * empty. It is empty if no @@ -187,7 +195,7 @@ class CompressedSparsityPattern : public Subscriptor * matrix. If the entry already * exists, nothing bad happens. */ - void add (const unsigned int i, + void add (const unsigned int i, const unsigned int j); /** @@ -197,7 +205,7 @@ class CompressedSparsityPattern : public Subscriptor * happens. */ template - void add_entries (const unsigned int row, + void add_entries (const unsigned int row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted = false); @@ -208,7 +216,7 @@ class CompressedSparsityPattern : public Subscriptor */ bool exists (const unsigned int i, const unsigned int j) const; - + /** * Make the sparsity pattern * symmetric by adding the @@ -221,7 +229,7 @@ class CompressedSparsityPattern : public Subscriptor * square matrix. */ void symmetrize (); - + /** * Print the sparsity of the * matrix. The output consists of @@ -285,6 +293,19 @@ class CompressedSparsityPattern : public Subscriptor unsigned int column_number (const unsigned int row, const unsigned int index) const; + /** + * Return an iterator that can loop over + * all entries in the given + * row. Dereferencing the iterator yields + * a column index. + */ + row_iterator row_begin (const unsigned int row) const; + + /** + * Returns the end of the current row. + */ + row_iterator row_end (const unsigned int row) const; + /** * Compute the bandwidth of the matrix * represented by this structure. The @@ -319,7 +340,7 @@ class CompressedSparsityPattern : public Subscriptor */ static bool stores_only_added_elements (); - + private: /** * Number of rows that this sparsity @@ -396,7 +417,7 @@ class CompressedSparsityPattern : public Subscriptor /** * Size of the cache. */ - static const unsigned int cache_size = 8; + static const unsigned int cache_size = 8; public: /** @@ -444,8 +465,8 @@ class CompressedSparsityPattern : public Subscriptor */ void flush_cache () const; }; - - + + /** * Actual data: store for each * row the set of nonzero @@ -474,7 +495,7 @@ CompressedSparsityPattern::Line::add (const unsigned int j) // the cache first if (cache_entries == cache_size && cache_entries != 0) flush_cache (); - + cache[cache_entries] = j; ++cache_entries; } @@ -564,6 +585,29 @@ CompressedSparsityPattern::column_number (const unsigned int row, +inline +CompressedSparsityPattern::row_iterator +CompressedSparsityPattern::row_begin (const unsigned int row) const +{ + Assert (row < n_rows(), ExcIndexRange (row, 0, n_rows())); + + if (lines[row].cache_entries != 0) + lines[row].flush_cache (); + return lines[row].entries.begin(); +} + + + +inline +CompressedSparsityPattern::row_iterator +CompressedSparsityPattern::row_end (const unsigned int row) const +{ + Assert (row < n_rows(), ExcIndexRange (row, 0, n_rows())); + return lines[row].entries.end(); +} + + + inline bool CompressedSparsityPattern::stores_only_added_elements () diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index f41edee86e..1226e4e420 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -136,7 +136,7 @@ SparseMatrix::operator = (const double d) Assert (cols->compressed || cols->empty(), SparsityPattern::ExcNotCompressed()); if (val) - std::fill_n (&val[0], cols->n_nonzero_elements(), 0.); + memset (&val[0], cols->n_nonzero_elements()*sizeof(number), 0); return *this; } diff --git a/deal.II/lac/include/lac/sparsity_pattern.h b/deal.II/lac/include/lac/sparsity_pattern.h index 27e8647348..d0cc5182f9 100644 --- a/deal.II/lac/include/lac/sparsity_pattern.h +++ b/deal.II/lac/include/lac/sparsity_pattern.h @@ -77,7 +77,7 @@ namespace internals if (len==0) return first; - + while (true) { // if length equals 8 or less, @@ -130,12 +130,12 @@ namespace internals Assert (false, ExcInternalError()); } } - - - + + + const unsigned int half = len >> 1; const unsigned int * const middle = first + half; - + // if the value is larger than // that pointed to by the // middle pointer, then the @@ -162,7 +162,7 @@ namespace internals */ unsigned int get_column_index_from_iterator (const unsigned int i); - + /** * Helper function to get the * column index from a @@ -186,17 +186,17 @@ namespace internals template unsigned int get_column_index_from_iterator (const std::pair &i); - + } } - + namespace SparsityPatternIterators { // forward declaration class Iterator; - + /** * Accessor class for iterators into * sparsity patterns. This class is also @@ -300,8 +300,8 @@ namespace SparsityPatternIterators /** @addtogroup Exceptions * @{ */ - - //@} + + //@} protected: /** * The sparsity pattern we operate on @@ -324,7 +324,7 @@ namespace SparsityPatternIterators * nonzero entry in the matrix. */ void advance (); - + /** * Grant access to iterator class. */ @@ -344,11 +344,11 @@ namespace SparsityPatternIterators * Constructor. Create an iterator * into the sparsity pattern @p sp for the * given row and the index within it. - */ + */ Iterator (const SparsityPattern *sp, const unsigned int row, const unsigned int index); - + /** * Prefix increment. */ @@ -403,7 +403,7 @@ namespace SparsityPatternIterators * accessor class. */ Accessor accessor; - }; + }; } @@ -428,7 +428,15 @@ class SparsityPattern : public Subscriptor typedef SparsityPatternIterators::Iterator const_iterator; - + + /** + * Typedef an iterator class that allows + * to walk over the nonzero elements of a + * row of a sparsity pattern. + */ + typedef + const unsigned int * row_iterator; + /** * Typedef an iterator class that allows * to walk over all nonzero elements of a @@ -442,8 +450,8 @@ class SparsityPattern : public Subscriptor typedef SparsityPatternIterators::Iterator iterator; - - + + /** * Define a value which is used * to indicate that a certain @@ -469,7 +477,7 @@ class SparsityPattern : public Subscriptor * variable may change over time. */ static const unsigned int invalid_entry = numbers::invalid_unsigned_int; - + /** * Initialize the matrix empty, * that is with no memory @@ -481,7 +489,7 @@ class SparsityPattern : public Subscriptor * the reinit() function. */ SparsityPattern (); - + /** * Copy constructor. This * constructor is only allowed to @@ -557,7 +565,7 @@ class SparsityPattern : public Subscriptor const unsigned int n, const std::vector& row_lengths, const bool optimize_diagonal = true); - + /** * Initialize a quadratic matrix * of dimension n with @@ -637,7 +645,7 @@ class SparsityPattern : public Subscriptor SparsityPattern (const SparsityPattern &original, const unsigned int max_per_row, const unsigned int extra_off_diagonals); - + /** * Destructor. */ @@ -652,7 +660,7 @@ class SparsityPattern : public Subscriptor * objects. */ SparsityPattern & operator = (const SparsityPattern &); - + /** * Reallocate memory and set up data * structures for a new matrix with @@ -707,7 +715,7 @@ class SparsityPattern : public Subscriptor const unsigned int n, const VectorSlice > &row_lengths, const bool optimize_diagonal = true); - + /** * This function compresses the sparsity * structure that this object represents. @@ -769,7 +777,34 @@ class SparsityPattern : public Subscriptor * iterator for the last row of a matrix. */ inline iterator end (const unsigned int r) const; - + + /** + * STL-like iterator with the first entry + * of row r. + * + * Note that if the given row is empty, + * i.e. does not contain any nonzero + * entries, then the iterator returned by + * this function equals + * end(r). Note also that the + * iterator may not be dereferencable in + * that case. + */ + inline row_iterator row_begin (const unsigned int r) const; + + /** + * Final iterator of row r. It + * points to the first element past the + * end of line @p r, or past the end of + * the entire sparsity pattern. + * + * Note that the end iterator is not + * necessarily dereferencable. This is in + * particular the case if it is the end + * iterator for the last row of a matrix. + */ + inline row_iterator row_end (const unsigned int r) const; + /** * This function can be used as a * replacement for reinit(), @@ -916,43 +951,18 @@ class SparsityPattern : public Subscriptor const bool optimize_diagonal = true); /** - * Copy data from an object of - * type - * CompressedSparsityPattern. - * Previous content of this - * object is lost, and the - * sparsity pattern is in - * compressed mode afterwards. - */ - 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); - - /** - * Copy data from an object of - * type + * Copy data from an object of type + * CompressedSparsityPattern, + * CompressedSetSparsityPattern or * CompressedSimpleSparsityPattern. - * Previous content of this - * object is lost, and the - * sparsity pattern is in + * Previous content of this object is + * lost, and the sparsity pattern is in * compressed mode afterwards. */ - void copy_from (const CompressedSimpleSparsityPattern &csp, + template + void copy_from (const CompressedSparsityType &csp, const bool optimize_diagonal = true); - /** * Take a full matrix and use its * nonzero entries to generate a @@ -967,7 +977,7 @@ class SparsityPattern : public Subscriptor template void copy_from (const FullMatrix &matrix, const bool optimize_diagonal = true); - + /** * Return whether the object is empty. It * is empty if no memory is allocated, @@ -1014,7 +1024,7 @@ class SparsityPattern : public Subscriptor * @deprecated Use * SparseMatrix::const_iterator */ - unsigned int operator() (const unsigned int i, + unsigned int operator() (const unsigned int i, const unsigned int j) const; /** @@ -1039,7 +1049,7 @@ class SparsityPattern : public Subscriptor */ std::pair matrix_position (const unsigned int global_index) const; - + /** * Add a nonzero entry to the matrix. * This function may only be called @@ -1048,7 +1058,7 @@ class SparsityPattern : public Subscriptor * If the entry already exists, nothing * bad happens. */ - void add (const unsigned int i, + void add (const unsigned int i, const unsigned int j); /** @@ -1061,11 +1071,11 @@ class SparsityPattern : public Subscriptor * exist, nothing bad happens. */ template - void add_entries (const unsigned int row, + void add_entries (const unsigned int row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted = false); - + /** * Make the sparsity pattern * symmetric by adding the @@ -1078,7 +1088,7 @@ class SparsityPattern : public Subscriptor * quadratic matrix. */ void symmetrize (); - + /** * Return number of rows of this * matrix, which equals the dimension @@ -1154,7 +1164,7 @@ class SparsityPattern : public Subscriptor unsigned int n_nonzero_elements () const; /** - * Return whether the structure is + * Return whether the structure is * compressed or not. */ bool is_compressed () const; @@ -1214,7 +1224,7 @@ class SparsityPattern : public Subscriptor * template arguments. */ bool stores_only_added_elements () const; - + /** * @deprecated * @@ -1278,7 +1288,7 @@ class SparsityPattern : public Subscriptor */ void partition (const unsigned int n_partitions, std::vector &partition_indices) const; - + /** * Write the data of this object * en bloc to a file. This is @@ -1325,7 +1335,7 @@ class SparsityPattern : public Subscriptor * more. */ void block_read (std::istream &in); - + /** * Print the sparsity of the * matrix. The output consists of @@ -1372,7 +1382,7 @@ class SparsityPattern : public Subscriptor unsigned int memory_consumption () const; /** - * This is kind of an expert mode. Get + * This is kind of an expert mode. Get * access to the rowstart array, but * read-only. * @@ -1624,8 +1634,8 @@ class SparsityPattern : public Subscriptor * Is special treatment of * diagonals enabled? */ - bool diagonal_optimized; - + bool diagonal_optimized; + /** * Make all sparse matrices * friends of this class. @@ -1644,7 +1654,7 @@ class SparsityPattern : public Subscriptor namespace SparsityPatternIterators -{ +{ inline Accessor:: Accessor (const SparsityPattern *sparsity_pattern, @@ -1693,7 +1703,7 @@ namespace SparsityPatternIterators Accessor::column() const { Assert (is_valid_entry() == true, ExcInvalidIterator()); - + return (sparsity_pattern ->get_column_numbers()[sparsity_pattern ->get_rowstart_indices()[a_row]+a_index]); @@ -1715,7 +1725,7 @@ namespace SparsityPatternIterators inline bool Accessor::operator == (const Accessor& other) const - { + { return (sparsity_pattern == other.sparsity_pattern && a_row == other.a_row && a_index == other.a_index); @@ -1726,22 +1736,22 @@ namespace SparsityPatternIterators inline bool Accessor::operator < (const Accessor& other) const - { + { Assert (sparsity_pattern == other.sparsity_pattern, ExcInternalError()); - + return (a_row < other.a_row || (a_row == other.a_row && a_index < other.a_index)); } - + inline void Accessor::advance () - { + { Assert (a_row < sparsity_pattern->n_rows(), ExcIteratorPastEnd()); - + ++a_index; // if at end of line: cycle until we @@ -1758,7 +1768,7 @@ namespace SparsityPatternIterators break; } } - + inline @@ -1832,7 +1842,7 @@ namespace SparsityPatternIterators { return accessor < other.accessor; } - + } @@ -1897,6 +1907,26 @@ SparsityPattern::end (const unsigned int r) const +inline +SparsityPattern::row_iterator +SparsityPattern::row_begin (const unsigned int r) const +{ + Assert (r(std::distance (begin, end)) == n_rows, ExcIteratorRange (std::distance (begin, end), n_rows)); - + // first determine row lengths for // each row. if the matrix is // quadratic, then we might have to @@ -2075,7 +2105,7 @@ SparsityPattern::copy_from (const unsigned int n_rows, const unsigned int col = internal::SparsityPatternTools::get_column_index_from_iterator(*j); Assert (col < n_cols, ExcIndexRange(col,0,n_cols)); - + if ((col!=row) || !is_square) *cols++ = col; }; diff --git a/deal.II/lac/source/chunk_sparsity_pattern.cc b/deal.II/lac/source/chunk_sparsity_pattern.cc index e6b1a1675e..e5e135391d 100644 --- a/deal.II/lac/source/chunk_sparsity_pattern.cc +++ b/deal.II/lac/source/chunk_sparsity_pattern.cc @@ -14,6 +14,7 @@ #include #include #include +#include #include @@ -35,7 +36,7 @@ ChunkSparsityPattern::ChunkSparsityPattern (const ChunkSparsityPattern &s) { Assert (s.rows == 0, ExcInvalidConstructorCall()); Assert (s.cols == 0, ExcInvalidConstructorCall()); - + reinit (0,0,0,0, false); } @@ -59,7 +60,7 @@ ChunkSparsityPattern::ChunkSparsityPattern ( const unsigned int n, const std::vector& row_lengths, const unsigned int chunk_size, - const bool optimize_diag) + const bool optimize_diag) { Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); @@ -81,7 +82,7 @@ ChunkSparsityPattern::ChunkSparsityPattern ( const unsigned int m, const std::vector& row_lengths, const unsigned int chunk_size, - const bool optimize_diag) + const bool optimize_diag) { Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); @@ -104,7 +105,7 @@ ChunkSparsityPattern::operator = (const ChunkSparsityPattern &s) // perform the checks in the underlying // object as well sparsity_pattern = s.sparsity_pattern; - + return *this; } @@ -136,8 +137,8 @@ ChunkSparsityPattern::reinit ( const bool optimize_diag) { Assert (row_lengths.size() == m, ExcInvalidNumber (m)); - Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); - + Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); + rows = m; cols = n; @@ -151,7 +152,7 @@ ChunkSparsityPattern::reinit ( // arithmetic equals // ((m+chunk_size-1)/chunk_size): const unsigned int m_chunks = (m+chunk_size-1) / chunk_size, - n_chunks = (n+chunk_size-1) / chunk_size; + n_chunks = (n+chunk_size-1) / chunk_size; // compute the maximum number of chunks in // each row. the passed array denotes the @@ -184,10 +185,11 @@ ChunkSparsityPattern::compress () +template void -ChunkSparsityPattern::copy_from (const CompressedSparsityPattern &csp, - const unsigned int chunk_size, - const bool optimize_diag) +ChunkSparsityPattern::copy_from (const SparsityType &csp, + const unsigned int chunk_size, + const bool optimize_diag) { Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); @@ -196,56 +198,21 @@ ChunkSparsityPattern::copy_from (const CompressedSparsityPattern &csp, // pattern std::vector entries_per_row (csp.n_rows(), 0); for (unsigned int row = 0; row 0, ExcInvalidNumber (chunk_size)); - - // count number of entries per row, then - // initialize the underlying sparsity - // pattern - std::vector entries_per_row (csp.n_rows(), 0); - for (unsigned int row = 0; row &matrix, reinit (matrix.m(), matrix.n(), entries_per_row, chunk_size, optimize_diag); - - // then actually fill it + + // then actually fill it for (unsigned int row=0; row> c; - AssertThrow (c == ']', ExcIO()); + AssertThrow (c == ']', ExcIO()); } @@ -594,12 +561,24 @@ ChunkSparsityPattern::memory_consumption () const // explicit instantiations template +void ChunkSparsityPattern::copy_from (const CompressedSparsityPattern &, + const unsigned int, + const bool); +template +void ChunkSparsityPattern::copy_from (const CompressedSetSparsityPattern &, + const unsigned int, + const bool); +template +void ChunkSparsityPattern::copy_from (const CompressedSimpleSparsityPattern &, + const unsigned int, + const bool); +template void ChunkSparsityPattern::copy_from (const FullMatrix &, const unsigned int, - bool); + const bool); template void ChunkSparsityPattern::copy_from (const FullMatrix &, const unsigned int, - bool); + const bool); DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/lac/source/compressed_simple_sparsity_pattern.cc b/deal.II/lac/source/compressed_simple_sparsity_pattern.cc index cb8b932238..224f760178 100644 --- a/deal.II/lac/source/compressed_simple_sparsity_pattern.cc +++ b/deal.II/lac/source/compressed_simple_sparsity_pattern.cc @@ -61,9 +61,9 @@ CompressedSimpleSparsityPattern::Line::add_entries (ForwardIterator begin, // actually doing something. ForwardIterator my_it = begin; unsigned int col = *my_it; - std::vector::iterator it = + std::vector::iterator it = std::lower_bound(entries.begin(), entries.end(), col); - while (*it == col) + while (*it == col) { ++my_it; if (my_it == end) @@ -149,10 +149,10 @@ CompressedSimpleSparsityPattern::Line::add_entries (ForwardIterator begin, entries.push_back(col); it = entries.end()-1; } - else { + else { // do a binary search to find the place // where to insert: - it2 = std::lower_bound(entries.begin(), entries.end(), col); + it2 = std::lower_bound(entries.begin(), entries.end(), col); // If this entry is a duplicate, continue // immediately Insert at the right place @@ -207,7 +207,7 @@ CompressedSimpleSparsityPattern::Line::memory_consumption () const return entries.capacity()*sizeof(unsigned int)+sizeof(Line); } - + CompressedSimpleSparsityPattern::CompressedSimpleSparsityPattern () : rows(0), @@ -234,7 +234,7 @@ CompressedSimpleSparsityPattern (const CompressedSimpleSparsityPattern &s) CompressedSimpleSparsityPattern::CompressedSimpleSparsityPattern (const unsigned int m, const unsigned int n, const IndexSet & rowset_ -) +) : rows(0), cols(0), @@ -278,15 +278,8 @@ CompressedSimpleSparsityPattern::reinit (const unsigned int m, rows = m; cols = n; rowset=rowset_; - if (rowset.size()==0) - { - rowset.set_size(m); - rowset.add_range(0,m); - } - - Assert( rowset.size()==m, ExcInternalError()); - std::vector new_lines (rowset.n_elements()); + std::vector new_lines (rowset.size()==0 ? rows : rowset.n_elements()); lines.swap (new_lines); } @@ -310,26 +303,27 @@ 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 +bool CompressedSimpleSparsityPattern::exists (const unsigned int i, const unsigned int j) const { Assert (i::const_iterator j=lines[row].entries.begin(); @@ -373,9 +368,9 @@ CompressedSimpleSparsityPattern::symmetrize () void CompressedSimpleSparsityPattern::print (std::ostream &out) const { - for (unsigned int row=0; row::const_iterator j=lines[row].entries.begin(); @@ -392,10 +387,11 @@ CompressedSimpleSparsityPattern::print (std::ostream &out) const void CompressedSimpleSparsityPattern::print_gnuplot (std::ostream &out) const -{ +{ for (unsigned int row=0; row::const_iterator j=lines[row].entries.begin(); @@ -409,7 +405,7 @@ CompressedSimpleSparsityPattern::print_gnuplot (std::ostream &out) const << -static_cast(rowindex) << std::endl; } - + AssertThrow (out, ExcIO()); } @@ -422,7 +418,8 @@ CompressedSimpleSparsityPattern::bandwidth () const unsigned int b=0; for (unsigned int row=0; row::const_iterator j=lines[row].entries.begin(); @@ -430,7 +427,7 @@ CompressedSimpleSparsityPattern::bandwidth () const if (static_cast(std::abs(static_cast(rowindex-*j))) > b) b = std::abs(static_cast(rowindex-*j)); } - + return b; } @@ -444,7 +441,7 @@ CompressedSimpleSparsityPattern::n_nonzero_elements () const { n += lines[i].entries.size(); } - + return n; } @@ -455,7 +452,7 @@ CompressedSimpleSparsityPattern::memory_consumption () const //TODO: IndexSet... unsigned int mem = sizeof(CompressedSimpleSparsityPattern); for (unsigned int i=0; i void -SparsityPattern::copy_from (const CompressedSparsityPattern &csp, +SparsityPattern::copy_from (const CSP &csp, const bool optimize_diag) { // first determine row lengths for @@ -563,74 +564,18 @@ SparsityPattern::copy_from (const CompressedSparsityPattern &csp, // 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()); + const bool do_diag_optimize = optimize_diag && (csp.n_rows() == csp.n_cols()); std::vector row_lengths (csp.n_rows()); for (unsigned int i=0; i row_lengths (csp.n_rows()); - for (unsigned int i=0; i row_lengths (csp.n_rows()); - for (unsigned int i=0; i @@ -1155,6 +1048,9 @@ partition (const unsigned int n_partitions, // explicit instantiations +template void SparsityPattern::copy_from (const CompressedSparsityPattern &, bool); +template void SparsityPattern::copy_from (const CompressedSetSparsityPattern &, bool); +template void SparsityPattern::copy_from (const CompressedSimpleSparsityPattern &, bool); template void SparsityPattern::copy_from (const FullMatrix &, bool); template void SparsityPattern::copy_from (const FullMatrix &, bool); diff --git a/deal.II/lac/source/trilinos_sparsity_pattern.cc b/deal.II/lac/source/trilinos_sparsity_pattern.cc index 4d2c07ee3e..b3ab4a11e6 100644 --- a/deal.II/lac/source/trilinos_sparsity_pattern.cc +++ b/deal.II/lac/source/trilinos_sparsity_pattern.cc @@ -294,7 +294,6 @@ namespace TrilinosWrappers static_cast(input_col_map.NumGlobalElements()), ExcDimensionMismatch (sp.n_cols(), input_col_map.NumGlobalElements())); - Assert (exchange_data == false, ExcNotImplemented()); column_space_map = std::auto_ptr (new Epetra_Map (input_col_map)); graph.reset(); @@ -326,75 +325,8 @@ namespace TrilinosWrappers ExcDimensionMismatch (graph->NumGlobalRows(), sp.n_rows())); - std::vector row_indices; const unsigned int n_rows = sp.n_rows(); - - for (unsigned int row=0; rowEpetra_CrsGraph::InsertGlobalIndices (row, row_length, - &row_indices[0]); - } - - compress(); - } - - - - template <> - void - SparsityPattern::reinit (const Epetra_Map &input_row_map, - const Epetra_Map &input_col_map, - const CompressedSimpleSparsityPattern &sp, - const bool exchange_data) - { - Assert (sp.n_rows() == - static_cast(input_row_map.NumGlobalElements()), - ExcDimensionMismatch (sp.n_rows(), - input_row_map.NumGlobalElements())); - Assert (sp.n_cols() == - static_cast(input_col_map.NumGlobalElements()), - ExcDimensionMismatch (sp.n_cols(), - input_col_map.NumGlobalElements())); - - column_space_map = std::auto_ptr (new Epetra_Map (input_col_map)); - graph.reset(); - compressed = false; - - Assert (input_row_map.LinearMap() == true, - ExcMessage ("This function is not efficient if the map is not contiguous.")); - - std::vector n_entries_per_row(input_row_map.MaxMyGID()- - input_row_map.MinMyGID() + 1); - - for (unsigned int row=input_row_map.MinMyGID(); - row(input_row_map.MaxMyGID()+1); - ++row) - n_entries_per_row[row-input_row_map.MinMyGID()] = sp.row_length(row); - - if (input_row_map.Comm().NumProc() > 1) - graph = std::auto_ptr - (new Epetra_FECrsGraph(Copy, input_row_map, - n_entries_per_row[0], - false)); - else - graph = std::auto_ptr - (new Epetra_FECrsGraph(Copy, input_row_map, input_col_map, - n_entries_per_row[0], - false)); - - Assert (graph->NumGlobalRows() == (int)sp.n_rows(), - ExcDimensionMismatch (graph->NumGlobalRows(), - sp.n_rows())); - - const unsigned int n_rows = sp.n_rows(); - std::vector row_indices; + std::vector row_indices; // Include possibility to exchange data // since CompressedSimpleSparsityPattern is @@ -403,85 +335,20 @@ namespace TrilinosWrappers if (input_row_map.MyGID(row) ) { const int row_length = sp.row_length(row); + if (row_length == 0) + continue; + row_indices.resize (row_length, -1); - for (int col=0; col < row_length; ++col) - row_indices[col] = sp.column_number (row, col); + typename SparsityType::row_iterator col_num = sp.row_begin (row), + row_end = sp.row_end(row); + for (unsigned int col = 0; col_num != row_end; ++col_num, ++col) + row_indices[col] = *col_num; graph->Epetra_CrsGraph::InsertGlobalIndices (row, row_length, &row_indices[0]); } - else if ( exchange_data && sp.row_index_set().is_element(row) ) - { - const int row_length = sp.row_length(row); - row_indices.resize (row_length, -1); - - for (int col=0; col < row_length; ++col) - row_indices[col] = sp.column_number (row, col); - - graph->InsertGlobalIndices (1, (int*)&row, row_length, &row_indices[0]); - } - - compress(); - } - - - - template<> - void - SparsityPattern::reinit (const Epetra_Map &input_row_map, - const Epetra_Map &input_col_map, - const CompressedSetSparsityPattern &sp, - const bool exchange_data) - { - Assert (exchange_data == false, ExcNotImplemented()); - Assert (sp.n_rows() == - static_cast(input_row_map.NumGlobalElements()), - ExcDimensionMismatch (sp.n_rows(), - input_row_map.NumGlobalElements())); - Assert (sp.n_cols() == - static_cast(input_col_map.NumGlobalElements()), - ExcDimensionMismatch (sp.n_cols(), - input_col_map.NumGlobalElements())); - - column_space_map = std::auto_ptr (new Epetra_Map (input_col_map)); - graph.reset(); - compressed = false; - - Assert (input_row_map.LinearMap() == true, - ExcMessage ("This function is not efficient if the map is not contiguous.")); - - std::vector n_entries_per_row(input_row_map.MaxMyGID()- - input_row_map.MinMyGID() + 1); - for (unsigned int row=input_row_map.MinMyGID(); - row(input_row_map.MaxMyGID()+1); - ++row) - { - n_entries_per_row[row-input_row_map.MinMyGID()] = sp.row_length(row); - } - - - if (input_row_map.Comm().NumProc() > 1) - graph = std::auto_ptr - (new Epetra_FECrsGraph(Copy, input_row_map, - n_entries_per_row[0], - false)); - else - graph = std::auto_ptr - (new Epetra_FECrsGraph(Copy, input_row_map, input_col_map, - n_entries_per_row[0], - false)); - - Assert (graph->NumGlobalRows() == (int)sp.n_rows(), - ExcDimensionMismatch (graph->NumGlobalRows(), - sp.n_rows())); - - - const unsigned int n_rows = sp.n_rows(); - std::vector row_indices; - - for (unsigned int row=0; rowEpetra_CrsGraph::InsertGlobalIndices (row, row_length, - &row_indices[0]); + graph->InsertGlobalIndices (1, (int*)&row, row_length, &row_indices[0]); } compress(); @@ -856,6 +719,16 @@ namespace TrilinosWrappers const Epetra_Map &, const dealii::CompressedSparsityPattern &, bool); + template void + SparsityPattern::reinit (const Epetra_Map &, + const Epetra_Map &, + const dealii::CompressedSetSparsityPattern &, + bool); + template void + SparsityPattern::reinit (const Epetra_Map &, + const Epetra_Map &, + const dealii::CompressedSimpleSparsityPattern &, + bool); } -- 2.39.5