From: David Wells Date: Tue, 20 Sep 2022 17:01:05 +0000 (-0400) Subject: Remove SparsityPatternBase. X-Git-Tag: v9.5.0-rc1~914^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F14304%2Fhead;p=dealii.git Remove SparsityPatternBase. We only have one class which inherits from this (SparsityPattern) and, since we instantiate templates for each sparsity pattern type, it isn't possible to do something useful by inheriting from this class. This was split off from SparsityPattern in 5213023ca165bfbd420423551c184a3e31886223 with the intent of writing some more general storage classes based on SparsityPatternBase, but that has never left the proof-of-concept stage. --- diff --git a/doc/news/changes/incompatibilities/20220922DavidWells b/doc/news/changes/incompatibilities/20220922DavidWells new file mode 100644 index 0000000000..cb13857e94 --- /dev/null +++ b/doc/news/changes/incompatibilities/20220922DavidWells @@ -0,0 +1,3 @@ +Removed: The base class of SparsityPattern, SparsityPatternBase, has been removed. +
+(David Wells, 2022/09/22) diff --git a/include/deal.II/lac/sparsity_pattern.h b/include/deal.II/lac/sparsity_pattern.h index 65d1034084..adc2459e6b 100644 --- a/include/deal.II/lac/sparsity_pattern.h +++ b/include/deal.II/lac/sparsity_pattern.h @@ -44,7 +44,6 @@ DEAL_II_NAMESPACE_OPEN // Forward declarations #ifndef DOXYGEN class SparsityPattern; -class SparsityPatternBase; class DynamicSparsityPattern; class ChunkSparsityPattern; template @@ -140,12 +139,12 @@ namespace SparsityPatternIterators /** * Constructor. */ - Accessor(const SparsityPatternBase *matrix, const std::size_t linear_index); + Accessor(const SparsityPattern *matrix, const std::size_t linear_index); /** * Constructor. Construct the end accessor for the given sparsity pattern. */ - Accessor(const SparsityPatternBase *matrix); + Accessor(const SparsityPattern *matrix); /** * Default constructor creating a dummy accessor. This constructor is here @@ -226,7 +225,7 @@ namespace SparsityPatternIterators /** * The sparsity pattern we operate on accessed. */ - const SparsityPatternBase *container; + const SparsityPattern *container; /** * Index in global sparsity pattern. This index represents the location @@ -264,7 +263,7 @@ namespace SparsityPatternIterators * SparsityPattern class for more information. * * @note This class operates directly on the internal data structures of the - * SparsityPatternBase class. As a consequence, some operations are cheap and + * SparsityPattern class. As a consequence, some operations are cheap and * some are not. In particular, it is cheap to access the column index of * the sparsity pattern entry pointed to. On the other hand, it is expensive * to access the row index (this requires $O(\log(N))$ operations for a @@ -287,14 +286,14 @@ namespace SparsityPatternIterators /** * Type of the stored pointer. */ - using container_pointer_type = SparsityPatternBase *; + using container_pointer_type = SparsityPattern *; /** * Constructor. Create an iterator into the sparsity pattern @p sp for the * given global index (i.e., the index of the given element counting from * the zeroth row). */ - Iterator(const SparsityPatternBase *sp, const std::size_t linear_index); + Iterator(const SparsityPattern *sp, const std::size_t linear_index); /** * Constructor. Create an iterator into the sparsity pattern @p sp for @@ -304,553 +303,6 @@ namespace SparsityPatternIterators }; } // namespace SparsityPatternIterators - - -/** - * A class that can store which elements of a matrix are nonzero (or, to be - * precise, may be nonzero) and for which we have to allocate memory to - * store their values. This class is an example of the "static" type of sparsity - * patterns (see - * @ref Sparsity). - * It uses the compressed row storage - * (CSR) format to store data, and is used as the basis for the - * derived SparsityPattern class, which is in turn used by the SparseMatrix - * class. - * - * The elements of a SparsityPatternBase, corresponding to the places where - * SparseMatrix objects can store nonzero entries, are stored row-by-row. - * The ordering of non-zero elements within each row (i.e. increasing - * column index order) depends on the derived classes. - */ -class SparsityPatternBase : public Subscriptor -{ -public: - /** - * Declare type for container size. - */ - using size_type = types::global_dof_index; - - /** - * Typedef an iterator class that allows to walk over all nonzero elements - * of a sparsity pattern. - */ - using const_iterator = SparsityPatternIterators::Iterator; - - /** - * Typedef an iterator class that allows to walk over all nonzero elements - * of a sparsity pattern. - * - * Since the iterator does not allow to modify the sparsity pattern, this - * type is the same as that for @p const_iterator. - */ - using iterator = SparsityPatternIterators::Iterator; - - /** - * Define a value which is used to indicate that a certain value in the - * #colnums array is unused, i.e. does not represent a certain column number - * index. - * - * Indices with this invalid value are used to insert new entries to the - * sparsity pattern using the add() member function, and are removed when - * calling compress(). - * - * You should not assume that the variable declared here has a certain - * value. The initialization is given here only to enable the compiler to - * perform some optimizations, but the actual value of the variable may - * change over time. - */ - static const size_type invalid_entry = numbers::invalid_size_type; - - /** - * @name Construction and Initialization - * - * Constructors, destructor, functions initializing, copying and filling an - * object. - * - * @{ - */ - - /** - * 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. - */ - SparsityPatternBase(); - - /** - * Destructor. - */ - ~SparsityPatternBase() override = default; - - /** - * Reallocate memory and set up data structures for a new matrix with @p m - * rows and @p n columns, with at most @p max_per_row - * nonzero entries per row. - * - * This function simply maps its operations to the other reinit() - * function. - */ - void - reinit(const size_type m, const size_type n, const unsigned int max_per_row); - - /** - * Reallocate memory for a matrix of size @p m times @p n. The number of - * entries for each row is taken from the array @p row_lengths which - * has to give this number of each row $i=1\ldots m$. - * - * If m*n==0 all memory is freed, resulting in a total - * reinitialization of the object. If it is nonzero, new memory is only - * allocated if the new size extends the old one. This is done to save time - * and to avoid fragmentation of the heap. - */ - void - reinit(const size_type m, - const size_type n, - const std::vector &row_lengths); - - /** - * Same as above, but with an ArrayView argument instead. - * - * The derived classes are responsible for implementation of this function. - */ - virtual void - reinit(const size_type m, - const size_type n, - const ArrayView &row_lengths) = 0; - - /** - * 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 quadratic matrix. - */ - void - symmetrize(); - - /** - * Add a nonzero entry to the matrix. This function may only be called for - * non-compressed sparsity patterns. - * - * If the entry already exists, nothing bad happens. - */ - void - add(const size_type i, const size_type j); - - /** - * @} - */ - - /** - * @name Iterators - * - * @{ - */ - - /** - * Iterator starting at the first entry of the matrix. The resulting - * iterator can be used to walk over all nonzero entries of the sparsity - * pattern. - * - * The order in which elements are accessed depends on the storage scheme - * implemented by derived classes. - */ - iterator - begin() const; - - /** - * Final iterator. - */ - iterator - end() const; - - /** - * Iterator starting at 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 dereferenceable in - * that case. - * - * The order in which elements are accessed depends on the storage scheme - * implemented by derived classes. - */ - iterator - begin(const size_type 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 dereferenceable. This is in - * particular the case if it is the end iterator for the last row of a - * matrix. - */ - iterator - end(const size_type r) const; - - - /** - * @} - */ - - /** - * @name Querying information - * - * @{ - */ - - /** - * Test for equality of two SparsityPatterns. - */ - bool - operator==(const SparsityPatternBase &) const; - - /** - * 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; - - /** - * Check if a value at a certain position may be non-zero. - */ - bool - exists(const size_type i, const size_type j) const; - - /** - * Return the maximum number of entries per row. Before compression, this - * equals the number given to the constructor, while after compression, it - * equals the maximum number of entries actually allocated by the user. - */ - size_type - max_entries_per_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. Consequently, the maximum - * bandwidth a $n\times m$ matrix can have is $\max\{n-1,m-1\}$, a diagonal - * matrix has bandwidth 0, and there are at most $2*q+1$ entries per row if - * the bandwidth is $q$. The returned quantity is sometimes called "half - * bandwidth" in the literature. - */ - size_type - 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. - */ - std::size_t - n_nonzero_elements() const; - - /** - * Return whether the structure is compressed or not. - */ - bool - is_compressed() const; - - /** - * Return number of rows of this matrix, which equals the dimension of the - * image space. - */ - size_type - n_rows() const; - - /** - * Return number of columns of this matrix, which equals the dimension of - * the range space. - */ - size_type - n_cols() const; - - /** - * Number of entries in a specific row. - */ - unsigned int - row_length(const size_type row) const; - - /** - * Determine an estimate for the memory consumption (in bytes) of this - * object. See MemoryConsumption. - */ - std::size_t - memory_consumption() const; - - /** - * @} - */ - - /** - * @name Accessing entries - * - * @{ - */ - - /** - * Access to column number field. Return the column number of the - * indexth entry in row. Note that if diagonal elements - * are optimized, the first element in each row is the diagonal element, - * i.e. column_number(row,0)==row. - * - * If the sparsity pattern is already compressed, then (except for the - * diagonal element), the entries are sorted by columns, i.e. - * column_number(row,i) < column_number(row,i+1). - */ - size_type - column_number(const size_type row, const unsigned int index) const; - - /** - * The index of a global matrix entry in its row. - * - * This function is analogous to operator(), but it computes the index not - * with respect to the total field, but only with respect to the row - * j. - */ - size_type - row_position(const size_type i, const size_type j) const; - - /** - * This is the inverse operation to operator()(): given a global index, find - * out row and column of the matrix entry to which it belongs. The returned - * value is the pair composed of row and column index. - * - * This function may only be called if the sparsity pattern is closed. The - * global index must then be between zero and n_nonzero_elements(). - * - * If N is the number of rows of this matrix, then the complexity - * of this function is log(N). - */ - std::pair - matrix_position(const std::size_t global_index) const; - - /** - * @} - */ - - /** - * @name Input/Output - * - * @{ - */ - - /** - * 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 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 plot command. - */ - void - print_gnuplot(std::ostream &out) const; - - /** - * Prints the sparsity of the matrix in a .svg file which can be opened in a - * web browser. The .svg file contains squares which correspond to the - * entries in the matrix. An entry in the matrix which contains a non-zero - * value corresponds with a red square while a zero-valued entry in the - * matrix correspond with a white square. - */ - void - print_svg(std::ostream &out) const; - - /** - * Write the data of this object to a stream for the purpose of - * serialization using the [BOOST serialization - * library](https://www.boost.org/doc/libs/1_74_0/libs/serialization/doc/index.html). - */ - template - void - save(Archive &ar, const unsigned int version) const; - - /** - * Read the data of this object from a stream for the purpose of - * serialization using the [BOOST serialization - * library](https://www.boost.org/doc/libs/1_74_0/libs/serialization/doc/index.html). - */ - template - void - load(Archive &ar, const unsigned int version); - -#ifdef DOXYGEN - /** - * Write and read the data of this object from a stream for the purpose - * of serialization using the [BOOST serialization - * library](https://www.boost.org/doc/libs/1_74_0/libs/serialization/doc/index.html). - */ - template - void - serialize(Archive &archive, const unsigned int version); -#else - // This macro defines the serialize() method that is compatible with - // the templated save() and load() method that have been implemented. - BOOST_SERIALIZATION_SPLIT_MEMBER() -#endif - - /** - * @} - */ - - /** - * @addtogroup Exceptions - * - * @{ - */ - - /** - * The operation is only allowed after the SparsityPattern has been set up - * and compress() was called. - */ - DeclExceptionMsg( - ExcNotCompressed, - "The operation you attempted is only allowed after the SparsityPattern " - "has been set up and compress() was called."); - - /** - * You tried to add an element to a row, but there was no space left. - */ - DeclException2(ExcNotEnoughSpace, - int, - int, - << "Upon entering a new entry to row " << arg1 - << ": there was no free entry any more. " << std::endl - << "(Maximum number of entries for this row: " << arg2 - << "; maybe the matrix is already compressed?)"); - - /** - * This operation changes the structure of the SparsityPattern and is not - * possible after compress() has been called. - */ - DeclExceptionMsg( - ExcMatrixIsCompressed, - "The operation you attempted changes the structure of the SparsityPattern " - "and is not possible after compress() has been called."); - - /** - * @} - */ - - -protected: - /** - * Maximum number of rows that can be stored in the #rowstart array. Since - * reallocation of that array only happens if the present one is too small, - * but never when the size of this matrix structure shrinks, #max_dim might - * be larger than #rows and in this case #rowstart has more elements than - * are used. - */ - size_type max_dim; - - /** - * Number of rows that this sparsity structure shall represent. - */ - size_type rows; - - /** - * Number of columns that this sparsity structure shall represent. - */ - size_type cols; - - /** - * Size of the actually allocated array #colnums. Here, the same applies as - * for the #rowstart array, i.e. it may be larger than the actually used - * part of the array. - */ - std::size_t max_vec_len; - - /** - * Maximum number of elements per row. This is set to the value given to the - * reinit() function (or to the constructor), or to the maximum row length - * computed from the vectors in case the more flexible constructors or - * reinit versions are called. Its value is more or less meaningless after - * compress() has been called. - */ - unsigned int max_row_length; - - /** - * Array which hold for each row which is the first element in #colnums - * belonging to that row. Note that the size of the array is one larger than - * the number of rows, because the last element is used for - * row=#rows, i.e. the row past the last used one. The value of - * #rowstart[#rows]} equals the index of the element past the end in - * #colnums; this way, we are able to write loops like for - * (i=rowstart[k]; i also for the last row. - * - * Note that the actual size of the allocated memory may be larger than the - * region that is used. The actual number of elements that was allocated is - * stored in #max_dim. - */ - std::unique_ptr rowstart; - - /** - * Array of column numbers. In this array, we store for each non-zero - * element its column number. The column numbers for the elements in row - * r are stored within the index range - * #rowstart[r]...#rowstart[r+1]. Therefore to find out - * whether a given element (r,c) exists, we have to check whether the - * column number c exists in the above-mentioned range within this - * array. If it exists, say at position p within this array, the - * value of the respective element in the sparse matrix will also be at - * position p of the values array of that class. - * - * At the beginning, all elements of this array are set to @p -1 indicating - * invalid (unused) column numbers (diagonal elements are preset if - * optimized storage is requested, though). Now, if nonzero elements are - * added, one column number in the row's respective range after the other is - * set to the column number of the added element. When compress is called, - * unused elements (indicated by column numbers @p -1) are eliminated by - * copying the column number of subsequent rows and the column numbers - * within each row (with possible exception of the diagonal element) are - * sorted, such that finding whether an element exists and determining its - * position can be done by a binary search. - */ - std::unique_ptr colnums; - - /** - * Store whether the compress() function was called for this object. - */ - bool compressed; - - // Make all sparse matrices friends of this class. - template - friend class SparseMatrix; - template - friend class SparseLUDecomposition; - template - friend class SparseILU; - template - friend class ChunkSparseMatrix; - - friend class ChunkSparsityPattern; - friend class DynamicSparsityPattern; - - // Also give access to internal details to the iterator/accessor classes. - friend class SparsityPatternIterators::Iterator; - friend class SparsityPatternIterators::Accessor; - friend class ChunkSparsityPatternIterators::Accessor; -}; - /** * This class stores a sparsity pattern in * the it->column() and * it->row(). */ -class SparsityPattern : public SparsityPatternBase +class SparsityPattern : public Subscriptor { public: /** * Declare type for container size. */ - using size_type = SparsityPatternBase::size_type; + using size_type = types::global_dof_index; + + /** + * Define a value which is used to indicate that a certain value in the + * #colnums array is unused, i.e. does not represent a certain column number + * index. + * + * Indices with this invalid value are used to insert new entries to the + * sparsity pattern using the add() member function, and are removed when + * calling compress(). + * + * You should not assume that the variable declared here has a certain + * value. The initialization is given here only to enable the compiler to + * perform some optimizations, but the actual value of the variable may + * change over time. + */ + static constexpr size_type invalid_entry = numbers::invalid_size_type; /** * Typedef an iterator class that allows to walk over all nonzero elements * of a sparsity pattern. */ - using const_iterator = SparsityPatternBase::const_iterator; + using const_iterator = SparsityPatternIterators::Iterator; /** * Typedef an iterator class that allows to walk over all nonzero elements @@ -907,14 +375,59 @@ public: * Since the iterator does not allow to modify the sparsity pattern, this * type is the same as that for @p const_iterator. */ - using iterator = SparsityPatternBase::iterator; + using iterator = SparsityPatternIterators::Iterator; + + /** + * @name Iterators + * + * @{ + */ + + /** + * Iterator starting at the first entry of the matrix. The resulting + * iterator can be used to walk over all nonzero entries of the sparsity + * pattern. + * + * The order in which elements are accessed depends on the storage scheme + * implemented by derived classes. + */ + iterator + begin() const; + + /** + * Final iterator. + */ + iterator + end() const; + + /** + * Iterator starting at 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 dereferenceable in + * that case. + * + * The order in which elements are accessed depends on the storage scheme + * implemented by derived classes. + */ + iterator + begin(const size_type 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 dereferenceable. This is in + * particular the case if it is the end iterator for the last row of a + * matrix. + */ + iterator + end(const size_type r) const; /** - * Since this class has to implement only one reinit() function, we need to - * bring all base reinit() functions into the scope so that the compiler can - * find them. + * @} */ - using SparsityPatternBase::reinit; /** * @name Construction and setup @@ -1020,11 +533,6 @@ public: const unsigned int max_per_row, const size_type extra_off_diagonals); - /** - * Destructor. - */ - ~SparsityPattern() override = default; - /** * 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 @@ -1034,14 +542,38 @@ public: operator=(const SparsityPattern &); /** - * Reallocate memory for a matrix of size @p m times @p n. The number of - * entries for each row is taken from the ArrayView @p row_lengths which - * has to give this number of each row $i=0\ldots m-1$. + * Reallocate memory for a matrix of size @p m times @p n. The number of + * entries for each row is taken from the ArrayView @p row_lengths which + * has to give this number of each row $i=0\ldots m-1$. + * + * If m*n==0 all memory is freed, resulting in a total + * reinitialization of the object. If it is nonzero, new memory is only + * allocated if the new size extends the old one. This is done to save time + * and to avoid fragmentation of the heap. + */ + void + reinit(const size_type m, + const size_type n, + const ArrayView &row_lengths); + + /** + * Same as the other reinit(), but uses a std::vector instead of an ArrayView. + */ + void + reinit(const size_type m, + const size_type n, + const std::vector &row_lengths); + + /** + * Reallocate memory and set up data structures for a new matrix with @p m + * rows and @p n columns, with at most @p max_per_row + * nonzero entries per row. + * + * This function simply maps its operations to the other reinit() + * function. */ - virtual void - reinit(const size_type m, - const size_type n, - const ArrayView &row_lengths) override; + void + reinit(const size_type m, const size_type n, const unsigned int max_per_row); /** * This function compresses the sparsity structure that this object @@ -1174,6 +706,15 @@ public: void copy_from(const FullMatrix &matrix); + /** + * Add a nonzero entry to the matrix. This function may only be called for + * non-compressed sparsity patterns. + * + * If the entry already exists, nothing bad happens. + */ + void + add(const size_type i, const size_type j); + /** * Add several nonzero entries to the specified matrix row. This function * may only be called for non-compressed sparsity patterns. @@ -1187,6 +728,16 @@ public: ForwardIterator end, const bool indices_are_sorted = false); + /** + * 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 quadratic matrix. + */ + void + symmetrize(); + /** * @} */ @@ -1199,6 +750,84 @@ public: * @{ */ + /** + * 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. + */ + std::size_t + n_nonzero_elements() const; + + /** + * Return number of rows of this matrix, which equals the dimension of the + * image space. + */ + size_type + n_rows() const; + + /** + * Return number of columns of this matrix, which equals the dimension of + * the range space. + */ + size_type + n_cols() const; + + /** + * 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; + + /** + * Check if a value at a certain position may be non-zero. + */ + bool + exists(const size_type i, const size_type j) const; + + /** + * This is the inverse operation to operator()(): given a global index, find + * out row and column of the matrix entry to which it belongs. The returned + * value is the pair composed of row and column index. + * + * This function may only be called if the sparsity pattern is closed. The + * global index must then be between zero and n_nonzero_elements(). + * + * If N is the number of rows of this matrix, then the complexity + * of this function is log(N). + */ + std::pair + matrix_position(const std::size_t global_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. Consequently, the maximum + * bandwidth a $n\times m$ matrix can have is $\max\{n-1,m-1\}$, a diagonal + * matrix has bandwidth 0, and there are at most $2*q+1$ entries per row if + * the bandwidth is $q$. The returned quantity is sometimes called "half + * bandwidth" in the literature. + */ + size_type + bandwidth() const; + + /** + * Return whether the structure is compressed or not. + */ + bool + is_compressed() const; + + /** + * Return the maximum number of entries per row. Before compression, this + * equals the number given to the constructor, while after compression, it + * equals the maximum number of entries actually allocated by the user. + */ + size_type + max_entries_per_row() const; + /** * Test for equality of two SparsityPatterns. */ @@ -1220,6 +849,12 @@ public: bool stores_only_added_elements() const; + /** + * Number of entries in a specific row. + */ + unsigned int + row_length(const size_type row) const; + /** * Determine an estimate for the memory consumption (in bytes) of this * object. See MemoryConsumption. @@ -1262,6 +897,29 @@ public: size_type operator()(const size_type i, const size_type j) const; + /** + * Access to column number field. Return the column number of the + * indexth entry in row. Note that if diagonal elements + * are optimized, the first element in each row is the diagonal element, + * i.e. column_number(row,0)==row. + * + * If the sparsity pattern is already compressed, then (except for the + * diagonal element), the entries are sorted by columns, i.e. + * column_number(row,i) < column_number(row,i+1). + */ + size_type + column_number(const size_type row, const unsigned int index) const; + + /** + * The index of a global matrix entry in its row. + * + * This function is analogous to operator(), but it computes the index not + * with respect to the total field, but only with respect to the row + * j. + */ + size_type + row_position(const size_type i, const size_type j) const; + /** * @} */ @@ -1273,6 +931,40 @@ public: * @{ */ + /** + * 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 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 plot command. + */ + void + print_gnuplot(std::ostream &out) const; + + /** + * Prints the sparsity of the matrix in a .svg file which can be opened in a + * web browser. The .svg file contains squares which correspond to the + * entries in the matrix. An entry in the matrix which contains a non-zero + * value corresponds with a red square while a zero-valued entry in the + * matrix correspond with a white square. + */ + void + print_svg(std::ostream &out) const; + /** * Write the data of this object en bloc to a file. This is done in a binary * mode, so the output is neither readable by humans nor (probably) by other @@ -1314,56 +1006,165 @@ public: * Read the data of this object from a stream for the purpose of * serialization */ - template - void - load(Archive &ar, const unsigned int version); + template + void + load(Archive &ar, const unsigned int version); + +#ifdef DOXYGEN + /** + * Write and read the data of this object from a stream for the purpose + * of serialization. + */ + template + void + serialize(Archive &archive, const unsigned int version); +#else + // This macro defines the serialize() method that is compatible with + // the templated save() and load() method that have been implemented. + BOOST_SERIALIZATION_SPLIT_MEMBER() +#endif + + /** + * @} + */ + + /** + * @addtogroup Exceptions + * + * @{ + */ + /** + * Exception + */ + DeclException2(ExcIteratorRange, + int, + int, + << "The iterators denote a range of " << arg1 + << " elements, but the given number of rows was " << arg2); + /** + * Exception + */ + DeclException1(ExcInvalidNumberOfPartitions, + int, + << "The number of partitions you gave is " << arg1 + << ", but must be greater than zero."); + + /** + * You tried to add an element to a row, but there was no space left. + */ + DeclException2(ExcNotEnoughSpace, + int, + int, + << "Upon entering a new entry to row " << arg1 + << ": there was no free entry any more. " << std::endl + << "(Maximum number of entries for this row: " << arg2 + << "; maybe the matrix is already compressed?)"); + + /** + * This operation changes the structure of the SparsityPattern and is not + * possible after compress() has been called. + */ + DeclExceptionMsg( + ExcMatrixIsCompressed, + "The operation you attempted changes the structure of the SparsityPattern " + "and is not possible after compress() has been called."); + + /** + * The operation is only allowed after the SparsityPattern has been set up + * and compress() was called. + */ + DeclExceptionMsg( + ExcNotCompressed, + "The operation you attempted is only allowed after the SparsityPattern " + "has been set up and compress() was called."); + + /** + * @} + */ +private: + /** + * Is special treatment of diagonals enabled? + */ + bool store_diagonal_first_in_row; + + /** + * Maximum number of rows that can be stored in the #rowstart array. Since + * reallocation of that array only happens if the present one is too small, + * but never when the size of this matrix structure shrinks, #max_dim might + * be larger than #rows and in this case #rowstart has more elements than + * are used. + */ + size_type max_dim; -#ifdef DOXYGEN /** - * Write and read the data of this object from a stream for the purpose - * of serialization. + * Number of rows that this sparsity structure shall represent. */ - template - void - serialize(Archive &archive, const unsigned int version); -#else - // This macro defines the serialize() method that is compatible with - // the templated save() and load() method that have been implemented. - BOOST_SERIALIZATION_SPLIT_MEMBER() -#endif + size_type rows; /** - * @} + * Number of columns that this sparsity structure shall represent. */ + size_type cols; /** - * @addtogroup Exceptions - * - * @{ + * Size of the actually allocated array #colnums. Here, the same applies as + * for the #rowstart array, i.e. it may be larger than the actually used + * part of the array. */ + std::size_t max_vec_len; + /** - * Exception + * Maximum number of elements per row. This is set to the value given to the + * reinit() function (or to the constructor), or to the maximum row length + * computed from the vectors in case the more flexible constructors or + * reinit versions are called. Its value is more or less meaningless after + * compress() has been called. */ - DeclException2(ExcIteratorRange, - int, - int, - << "The iterators denote a range of " << arg1 - << " elements, but the given number of rows was " << arg2); + unsigned int max_row_length; + /** - * Exception + * Array which hold for each row which is the first element in #colnums + * belonging to that row. Note that the size of the array is one larger than + * the number of rows, because the last element is used for + * row=#rows, i.e. the row past the last used one. The value of + * #rowstart[#rows]} equals the index of the element past the end in + * #colnums; this way, we are able to write loops like for + * (i=rowstart[k]; i also for the last row. + * + * Note that the actual size of the allocated memory may be larger than the + * region that is used. The actual number of elements that was allocated is + * stored in #max_dim. */ - DeclException1(ExcInvalidNumberOfPartitions, - int, - << "The number of partitions you gave is " << arg1 - << ", but must be greater than zero."); + std::unique_ptr rowstart; + /** - * @} + * Array of column numbers. In this array, we store for each non-zero + * element its column number. The column numbers for the elements in row + * r are stored within the index range + * #rowstart[r]...#rowstart[r+1]. Therefore to find out + * whether a given element (r,c) exists, we have to check whether the + * column number c exists in the above-mentioned range within this + * array. If it exists, say at position p within this array, the + * value of the respective element in the sparse matrix will also be at + * position p of the values array of that class. + * + * At the beginning, all elements of this array are set to @p -1 indicating + * invalid (unused) column numbers (diagonal elements are preset if + * optimized storage is requested, though). Now, if nonzero elements are + * added, one column number in the row's respective range after the other is + * set to the column number of the added element. When compress is called, + * unused elements (indicated by column numbers @p -1) are eliminated by + * copying the column number of subsequent rows and the column numbers + * within each row (with possible exception of the diagonal element) are + * sorted, such that finding whether an element exists and determining its + * position can be done by a binary search. */ -private: + std::unique_ptr colnums; + /** - * Is special treatment of diagonals enabled? + * Store whether the compress() function was called for this object. */ - bool store_diagonal_first_in_row; + bool compressed; // Make all sparse matrices friends of this class. template @@ -1395,15 +1196,15 @@ private: namespace SparsityPatternIterators { - inline Accessor::Accessor(const SparsityPatternBase *sparsity_pattern, - const std::size_t i) + inline Accessor::Accessor(const SparsityPattern *sparsity_pattern, + const std::size_t i) : container(sparsity_pattern) , linear_index(i) {} - inline Accessor::Accessor(const SparsityPatternBase *sparsity_pattern) + inline Accessor::Accessor(const SparsityPattern *sparsity_pattern) : container(sparsity_pattern) , linear_index(container->rowstart[container->rows]) {} @@ -1501,8 +1302,8 @@ namespace SparsityPatternIterators } - inline Iterator::Iterator(const SparsityPatternBase *sp, - const std::size_t linear_index) + inline Iterator::Iterator(const SparsityPattern *sp, + const std::size_t linear_index) : LinearIndexIterator(Accessor(sp, linear_index)) {} @@ -1516,58 +1317,30 @@ namespace SparsityPatternIterators -inline SparsityPatternBase::iterator -SparsityPatternBase::begin() const +inline std::size_t +SparsityPattern::n_nonzero_elements() const { - if (n_rows() > 0) - return {this, rowstart[0]}; - else - return end(); -} - - + Assert(compressed, ExcNotCompressed()); -inline SparsityPatternBase::iterator -SparsityPatternBase::end() const -{ - if (n_rows() > 0) - return {this, rowstart[rows]}; + if ((rowstart != nullptr) && (colnums != nullptr)) + return rowstart[rows] - rowstart[0]; else - return {nullptr, 0}; -} - - - -inline SparsityPatternBase::iterator -SparsityPatternBase::begin(const size_type r) const -{ - AssertIndexRange(r, n_rows()); - - return {this, rowstart[r]}; -} - - - -inline SparsityPatternBase::iterator -SparsityPatternBase::end(const size_type r) const -{ - AssertIndexRange(r, n_rows()); - - return {this, rowstart[r + 1]}; + // the object is empty or has zero size + return 0; } -inline SparsityPatternBase::size_type -SparsityPatternBase::n_rows() const +inline SparsityPattern::size_type +SparsityPattern::n_rows() const { return rows; } -inline SparsityPatternBase::size_type -SparsityPatternBase::n_cols() const +inline SparsityPattern::size_type +SparsityPattern::n_cols() const { return cols; } @@ -1575,7 +1348,7 @@ SparsityPatternBase::n_cols() const inline bool -SparsityPatternBase::is_compressed() const +SparsityPattern::is_compressed() const { return compressed; } @@ -1591,7 +1364,7 @@ SparsityPattern::stores_only_added_elements() const inline unsigned int -SparsityPatternBase::row_length(const size_type row) const +SparsityPattern::row_length(const size_type row) const { AssertIndexRange(row, rows); return rowstart[row + 1] - rowstart[row]; @@ -1600,8 +1373,8 @@ SparsityPatternBase::row_length(const size_type row) const inline SparsityPattern::size_type -SparsityPatternBase::column_number(const size_type row, - const unsigned int index) const +SparsityPattern::column_number(const size_type row, + const unsigned int index) const { AssertIndexRange(row, rows); AssertIndexRange(index, row_length(row)); @@ -1611,95 +1384,54 @@ SparsityPatternBase::column_number(const size_type row, -inline std::size_t -SparsityPatternBase::n_nonzero_elements() const +inline SparsityPattern::iterator +SparsityPattern::begin() const { - Assert(compressed, ExcNotCompressed()); - - if ((rowstart != nullptr) && (colnums != nullptr)) - return rowstart[rows] - rowstart[0]; + if (n_rows() > 0) + return {this, rowstart[0]}; else - // the object is empty or has zero size - return 0; + return end(); } -template -inline void -SparsityPatternBase::save(Archive &ar, const unsigned int) const +inline SparsityPattern::iterator +SparsityPattern::end() const { - // forward to serialization function in the base class. - ar &boost::serialization::base_object(*this); - - ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed; - - if (max_dim != 0) - ar &boost::serialization::make_array(rowstart.get(), max_dim + 1); - else - Assert(rowstart.get() == nullptr, ExcInternalError()); - - if (max_vec_len != 0) - ar &boost::serialization::make_array(colnums.get(), max_vec_len); + if (n_rows() > 0) + return {this, rowstart[rows]}; else - Assert(colnums.get() == nullptr, ExcInternalError()); + return {nullptr, 0}; } -template -inline void -SparsityPatternBase::load(Archive &ar, const unsigned int) +inline SparsityPattern::iterator +SparsityPattern::begin(const size_type r) const { - // forward to serialization function in the base class. - ar &boost::serialization::base_object(*this); - - ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed; - - if (max_dim != 0) - { - rowstart = std::make_unique(max_dim + 1); - ar &boost::serialization::make_array(rowstart.get(), max_dim + 1); - } - else - rowstart.reset(); + AssertIndexRange(r, n_rows()); - if (max_vec_len != 0) - { - colnums = std::make_unique(max_vec_len); - ar &boost::serialization::make_array(colnums.get(), max_vec_len); - } - else - colnums.reset(); + return {this, rowstart[r]}; } -template -inline void -SparsityPattern::save(Archive &ar, const unsigned int) const +inline SparsityPattern::iterator +SparsityPattern::end(const size_type r) const { - // forward to serialization function in the base class. - ar &boost::serialization::base_object(*this); - ar &store_diagonal_first_in_row; -} - - + AssertIndexRange(r, n_rows()); -template -inline void -SparsityPattern::load(Archive &ar, const unsigned int) -{ - // forward to serialization function in the base class. - ar &boost::serialization::base_object(*this); - ar &store_diagonal_first_in_row; + return {this, rowstart[r + 1]}; } inline bool -SparsityPatternBase::operator==(const SparsityPatternBase &sp2) const +SparsityPattern::operator==(const SparsityPattern &sp2) const { + if (store_diagonal_first_in_row != sp2.store_diagonal_first_in_row) + return false; + // it isn't quite necessary to compare *all* member variables. by only // comparing the essential ones, we can say that two sparsity patterns are // equal even if one is compressed and the other is not (in which case some @@ -1722,15 +1454,6 @@ SparsityPatternBase::operator==(const SparsityPatternBase &sp2) const -inline bool -SparsityPattern::operator==(const SparsityPattern &sp2) const -{ - return (static_cast(*this) == sp2) && - (store_diagonal_first_in_row == sp2.store_diagonal_first_in_row); -} - - - namespace internal { namespace SparsityPatternTools @@ -1818,6 +1541,58 @@ SparsityPattern::copy_from(const size_type n_rows, } + +template +inline void +SparsityPattern::save(Archive &ar, const unsigned int) const +{ + // forward to serialization function in the base class. + ar &boost::serialization::base_object(*this); + + ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed; + + if (max_dim != 0) + ar &boost::serialization::make_array(rowstart.get(), max_dim + 1); + else + Assert(rowstart.get() == nullptr, ExcInternalError()); + + if (max_vec_len != 0) + ar &boost::serialization::make_array(colnums.get(), max_vec_len); + else + Assert(colnums.get() == nullptr, ExcInternalError()); + ar &store_diagonal_first_in_row; +} + + + +template +inline void +SparsityPattern::load(Archive &ar, const unsigned int) +{ + // forward to serialization function in the base class. + ar &boost::serialization::base_object(*this); + + ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed; + + if (max_dim != 0) + { + rowstart = std::make_unique(max_dim + 1); + ar &boost::serialization::make_array(rowstart.get(), max_dim + 1); + } + else + rowstart.reset(); + + if (max_vec_len != 0) + { + colnums = std::make_unique(max_vec_len); + ar &boost::serialization::make_array(colnums.get(), max_vec_len); + } + else + colnums.reset(); + ar &store_diagonal_first_in_row; +} + + #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE diff --git a/source/lac/sparsity_pattern.cc b/source/lac/sparsity_pattern.cc index b59e976ed6..d6958cbbfd 100644 --- a/source/lac/sparsity_pattern.cc +++ b/source/lac/sparsity_pattern.cc @@ -31,28 +31,18 @@ DEAL_II_NAMESPACE_OPEN -#ifdef DEAL_II_MSVC -__declspec(selectany) // Weak extern binding due to multiple link error -#endif - const SparsityPatternBase::size_type SparsityPatternBase::invalid_entry; +constexpr SparsityPattern::size_type SparsityPattern::invalid_entry; + -SparsityPatternBase::SparsityPatternBase() - : max_dim(0) +SparsityPattern::SparsityPattern() + : store_diagonal_first_in_row(false) + , max_dim(0) , rows(0) , cols(0) , max_vec_len(0) , max_row_length(0) - , rowstart(nullptr) - , colnums(nullptr) , compressed(false) -{} - - - -SparsityPattern::SparsityPattern() - : SparsityPatternBase() - , store_diagonal_first_in_row(false) { reinit(0, 0, 0); } @@ -60,8 +50,7 @@ SparsityPattern::SparsityPattern() SparsityPattern::SparsityPattern(const SparsityPattern &s) - : SparsityPatternBase() - , store_diagonal_first_in_row(false) + : SparsityPattern() { (void)s; Assert(s.empty(), @@ -78,8 +67,7 @@ SparsityPattern::SparsityPattern(const SparsityPattern &s) SparsityPattern::SparsityPattern(const size_type m, const size_type n, const unsigned int max_per_row) - : SparsityPatternBase() - , store_diagonal_first_in_row(m == n) + : SparsityPattern() { reinit(m, n, max_per_row); } @@ -89,8 +77,7 @@ SparsityPattern::SparsityPattern(const size_type m, SparsityPattern::SparsityPattern(const size_type m, const size_type n, const std::vector &row_lengths) - : SparsityPatternBase() - , store_diagonal_first_in_row(m == n) + : SparsityPattern() { reinit(m, n, row_lengths); } @@ -99,7 +86,7 @@ SparsityPattern::SparsityPattern(const size_type m, SparsityPattern::SparsityPattern(const size_type m, const unsigned int max_per_row) - : SparsityPatternBase() + : SparsityPattern() { reinit(m, m, max_per_row); } @@ -108,7 +95,7 @@ SparsityPattern::SparsityPattern(const size_type m, SparsityPattern::SparsityPattern(const size_type m, const std::vector &row_lengths) - : SparsityPatternBase() + : SparsityPattern() { reinit(m, m, row_lengths); } @@ -213,18 +200,6 @@ SparsityPattern::operator=(const SparsityPattern &s) -void -SparsityPatternBase::reinit(const size_type m, - const size_type n, - const unsigned int max_per_row) -{ - // simply map this function to the other @p{reinit} function - const std::vector row_lengths(m, max_per_row); - reinit(m, n, row_lengths); -} - - - void SparsityPattern::reinit(const size_type m, const size_type n, @@ -334,6 +309,28 @@ SparsityPattern::reinit(const size_type m, +void +SparsityPattern::reinit(const size_type m, + const size_type n, + const std::vector &row_lengths) +{ + reinit(m, n, make_array_view(row_lengths)); +} + + + +void +SparsityPattern::reinit(const size_type m, + const size_type n, + const unsigned int max_per_row) +{ + // simply map this function to the other @p{reinit} function + const std::vector row_lengths(m, max_per_row); + reinit(m, n, row_lengths); +} + + + void SparsityPattern::compress() { @@ -610,18 +607,9 @@ SparsityPattern::copy_from(const FullMatrix &matrix) } -void -SparsityPatternBase::reinit(const size_type m, - const size_type n, - const std::vector &row_lengths) -{ - reinit(m, n, make_array_view(row_lengths)); -} - - bool -SparsityPatternBase::empty() const +SparsityPattern::empty() const { // let's try to be on the safe side of life by using multiple possibilities // in the check for emptiness... (sorry for this kludge -- emptying matrices @@ -643,8 +631,71 @@ SparsityPatternBase::empty() const -SparsityPatternBase::size_type -SparsityPatternBase::max_entries_per_row() const +bool +SparsityPattern::exists(const size_type i, const size_type j) const +{ + Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); + AssertIndexRange(i, rows); + AssertIndexRange(j, cols); + + for (size_type k = rowstart[i]; k < rowstart[i + 1]; ++k) + { + // entry already exists + if (colnums[k] == j) + return true; + } + return false; +} + + + +std::pair +SparsityPattern::matrix_position(const std::size_t global_index) const +{ + Assert(compressed == true, ExcNotCompressed()); + AssertIndexRange(global_index, n_nonzero_elements()); + + // first find the row in which the entry is located. for this note that the + // rowstart array indexes the global indices at which each row starts. since + // it is sorted, and since there is an element for the one-past-last row, we + // can simply use a bisection search on it + const size_type row = + (std::upper_bound(rowstart.get(), rowstart.get() + rows, global_index) - + rowstart.get() - 1); + + // now, the column index is simple since that is what the colnums array + // stores: + const size_type col = colnums[global_index]; + + // so return the respective pair + return std::make_pair(row, col); +} + + + +SparsityPattern::size_type +SparsityPattern::bandwidth() const +{ + Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); + size_type b = 0; + for (size_type i = 0; i < rows; ++i) + for (size_type j = rowstart[i]; j < rowstart[i + 1]; ++j) + if (colnums[j] != invalid_entry) + { + if (static_cast( + std::abs(static_cast(i - colnums[j]))) > b) + b = std::abs(static_cast(i - colnums[j])); + } + else + // leave if at the end of the entries of this line + break; + return b; +} + + + +SparsityPattern::size_type +SparsityPattern::max_entries_per_row() const { // if compress() has not yet been called, we can get the maximum number of // elements per row using the stored value @@ -701,7 +752,7 @@ SparsityPattern::operator()(const size_type i, const size_type j) const void -SparsityPatternBase::add(const size_type i, const size_type j) +SparsityPattern::add(const size_type i, const size_type j) { Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); AssertIndexRange(i, rows); @@ -777,68 +828,9 @@ SparsityPattern::add_entries(const size_type row, } -bool -SparsityPatternBase::exists(const size_type i, const size_type j) const -{ - Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); - AssertIndexRange(i, rows); - AssertIndexRange(j, cols); - - for (size_type k = rowstart[i]; k < rowstart[i + 1]; ++k) - { - // entry already exists - if (colnums[k] == j) - return true; - } - return false; -} - - - -SparsityPatternBase::size_type -SparsityPatternBase::row_position(const size_type i, const size_type j) const -{ - Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); - AssertIndexRange(i, rows); - AssertIndexRange(j, cols); - - for (size_type k = rowstart[i]; k < rowstart[i + 1]; ++k) - { - // entry exists - if (colnums[k] == j) - return k - rowstart[i]; - } - return numbers::invalid_size_type; -} - - - -std::pair -SparsityPatternBase::matrix_position(const std::size_t global_index) const -{ - Assert(compressed == true, ExcNotCompressed()); - AssertIndexRange(global_index, n_nonzero_elements()); - - // first find the row in which the entry is located. for this note that the - // rowstart array indexes the global indices at which each row starts. since - // it is sorted, and since there is an element for the one-past-last row, we - // can simply use a bisection search on it - const size_type row = - (std::upper_bound(rowstart.get(), rowstart.get() + rows, global_index) - - rowstart.get() - 1); - - // now, the column index is simple since that is what the colnums array - // stores: - const size_type col = colnums[global_index]; - - // so return the respective pair - return std::make_pair(row, col); -} - - void -SparsityPatternBase::symmetrize() +SparsityPattern::symmetrize() { Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); Assert(compressed == false, ExcMatrixIsCompressed()); @@ -871,8 +863,26 @@ SparsityPatternBase::symmetrize() +SparsityPattern::size_type +SparsityPattern::row_position(const size_type i, const size_type j) const +{ + Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); + AssertIndexRange(i, rows); + AssertIndexRange(j, cols); + + for (size_type k = rowstart[i]; k < rowstart[i + 1]; ++k) + { + // entry exists + if (colnums[k] == j) + return k - rowstart[i]; + } + return numbers::invalid_size_type; +} + + + void -SparsityPatternBase::print(std::ostream &out) const +SparsityPattern::print(std::ostream &out) const { Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); @@ -893,7 +903,7 @@ SparsityPatternBase::print(std::ostream &out) const void -SparsityPatternBase::print_gnuplot(std::ostream &out) const +SparsityPattern::print_gnuplot(std::ostream &out) const { Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); @@ -910,8 +920,10 @@ SparsityPatternBase::print_gnuplot(std::ostream &out) const AssertThrow(out.fail() == false, ExcIO()); } + + void -SparsityPatternBase::print_svg(std::ostream &out) const +SparsityPattern::print_svg(std::ostream &out) const { const unsigned int m = this->n_rows(); const unsigned int n = this->n_cols(); @@ -943,26 +955,6 @@ SparsityPatternBase::print_svg(std::ostream &out) const -SparsityPatternBase::size_type -SparsityPatternBase::bandwidth() const -{ - Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); - size_type b = 0; - for (size_type i = 0; i < rows; ++i) - for (size_type j = rowstart[i]; j < rowstart[i + 1]; ++j) - if (colnums[j] != invalid_entry) - { - if (static_cast( - std::abs(static_cast(i - colnums[j]))) > b) - b = std::abs(static_cast(i - colnums[j])); - } - else - // leave if at the end of the entries of this line - break; - return b; -} - - void SparsityPattern::block_write(std::ostream &out) const { @@ -1027,20 +1019,13 @@ SparsityPattern::block_read(std::istream &in) std::size_t -SparsityPatternBase::memory_consumption() const +SparsityPattern::memory_consumption() const { return (max_dim * sizeof(size_type) + sizeof(*this) + max_vec_len * sizeof(size_type)); } -std::size_t -SparsityPattern::memory_consumption() const -{ - return sizeof(*this) + SparsityPatternBase::memory_consumption(); -} - - #ifndef DOXYGEN // explicit instantiations diff --git a/tests/serialization/sparsity_pattern.output b/tests/serialization/sparsity_pattern.output index abc3bde154..cbc87c3039 100644 --- a/tests/serialization/sparsity_pattern.output +++ b/tests/serialization/sparsity_pattern.output @@ -1,6 +1,6 @@ -DEAL::0 0 0 0 0 0 16 16 16 64 5 1 0 3 7 11 14 18 23 28 32 36 41 46 50 53 57 61 64 0 1 4 1 0 2 5 2 1 3 6 3 2 7 4 0 5 8 5 1 4 6 9 6 2 5 7 10 7 3 6 11 8 4 9 12 9 5 8 10 13 10 6 9 11 14 11 7 10 15 12 8 13 13 9 12 14 14 10 13 15 15 11 14 1 +DEAL::0 0 0 0 16 16 16 64 5 1 0 3 7 11 14 18 23 28 32 36 41 46 50 53 57 61 64 0 1 4 1 0 2 5 2 1 3 6 3 2 7 4 0 5 8 5 1 4 6 9 6 2 5 7 10 7 3 6 11 8 4 9 12 9 5 8 10 13 10 6 9 11 14 11 7 10 15 12 8 13 13 9 12 14 14 10 13 15 15 11 14 1 -DEAL::0 0 0 0 0 0 16 16 16 64 5 1 0 3 7 11 14 18 23 28 32 36 41 46 50 53 57 61 64 0 1 4 1 0 2 5 2 1 3 6 3 2 7 4 0 5 8 5 1 4 6 9 6 2 5 7 10 7 3 6 11 8 4 9 12 9 5 8 10 13 10 6 9 11 14 11 7 10 15 12 8 13 13 9 12 14 14 10 13 15 15 11 14 1 +DEAL::0 0 0 0 16 16 16 64 5 1 0 3 7 11 14 18 23 28 32 36 41 46 50 53 57 61 64 0 1 4 1 0 2 5 2 1 3 6 3 2 7 4 0 5 8 5 1 4 6 9 6 2 5 7 10 7 3 6 11 8 4 9 12 9 5 8 10 13 10 6 9 11 14 11 7 10 15 12 8 13 13 9 12 14 14 10 13 15 15 11 14 1 DEAL::OK diff --git a/tests/serialization/sparsity_pattern_empty.output b/tests/serialization/sparsity_pattern_empty.output index aa452f6475..9f62171624 100644 --- a/tests/serialization/sparsity_pattern_empty.output +++ b/tests/serialization/sparsity_pattern_empty.output @@ -1,6 +1,6 @@ -DEAL::0 0 0 0 0 0 16 16 16 16 5 1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 +DEAL::0 0 0 0 16 16 16 16 5 1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 -DEAL::0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL::0 0 0 0 0 0 0 0 0 0 0 DEAL::OK diff --git a/tests/sparsity/sparsity_pattern_13.cc b/tests/sparsity/sparsity_pattern_13.cc deleted file mode 100644 index 7d1a542f53..0000000000 --- a/tests/sparsity/sparsity_pattern_13.cc +++ /dev/null @@ -1,169 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2018 - 2020 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE.md at -// the top level directory of deal.II. -// -// --------------------------------------------------------------------- - - - -// proof of concept for DynamicSparsityPattern-like CSR object without -// special treatment of diagonals - -#include "sparsity_pattern_common.h" - -class SparsityPatternStandard : public SparsityPatternBase -{ -public: - using size_type = SparsityPatternBase::size_type; - - using SparsityPatternBase::reinit; - - SparsityPatternStandard() - : SparsityPatternBase() - { - reinit(0, 0, 0); - }; - - virtual void - reinit(const size_type m, - const size_type n, - const ArrayView &row_lengths) override - { - AssertDimension(row_lengths.size(), m); - - rows = m; - cols = n; - - // delete empty matrices - if ((m == 0) || (n == 0)) - { - rowstart.reset(); - colnums.reset(); - - max_vec_len = max_dim = rows = cols = 0; - // if dimension is zero: ignore max_per_row - max_row_length = 0; - compressed = false; - - return; - } - - // find out how many entries we need in the @p{colnums} array. if this - // number is larger than @p{max_vec_len}, then we will need to reallocate - // memory - // - // note that the number of elements per row is bounded by the number of - // columns - // - std::size_t vec_len = 0; - for (size_type i = 0; i < m; ++i) - vec_len += std::min(static_cast(row_lengths[i]), n); - - // sometimes, no entries are requested in the matrix (this most often - // happens when blocks in a block matrix are simply zero). in that case, - // allocate exactly one element, to have a valid pointer to some memory - if (vec_len == 0) - { - vec_len = 1; - max_vec_len = vec_len; - colnums = std::make_unique(max_vec_len); - } - - max_row_length = - (row_lengths.size() == 0 ? - 0 : - std::min(static_cast( - *std::max_element(row_lengths.begin(), row_lengths.end())), - n)); - - // allocate memory for the rowstart values, if necessary. even though we - // re-set the pointers again immediately after deleting their old content, - // set them to zero in between because the allocation might fail, in which - // case we get an exception and the destructor of this object will be called - // -- where we look at the non-nullness of the (now invalid) pointer again - // and try to delete the memory a second time. - if (rows > max_dim) - { - max_dim = rows; - rowstart = std::make_unique(max_dim + 1); - } - - // allocate memory for the column numbers if necessary - if (vec_len > max_vec_len) - { - max_vec_len = vec_len; - colnums = std::make_unique(max_vec_len); - } - - // set the rowstart array - rowstart[0] = 0; - for (size_type i = 1; i <= rows; ++i) - rowstart[i] = rowstart[i - 1] + - std::min(static_cast(row_lengths[i - 1]), n); - Assert((rowstart[rows] == vec_len) || - ((vec_len == 1) && (rowstart[rows] == 0)), - ExcInternalError()); - - // preset the column numbers by a value indicating it is not in use - std::fill_n(colnums.get(), vec_len, invalid_entry); - - compressed = false; - }; - - void - copy_from(const DynamicSparsityPattern &dsp) - { - std::vector row_lengths(dsp.n_rows()); - for (size_type i = 0; i < dsp.n_rows(); ++i) - row_lengths[i] = dsp.row_length(i); - - reinit(dsp.n_rows(), dsp.n_cols(), row_lengths); - - if (n_rows() != 0 && n_cols() != 0) - for (size_type row = 0; row < dsp.n_rows(); ++row) - { - size_type * cols = &colnums[rowstart[row]]; - const unsigned int row_length = dsp.row_length(row); - for (unsigned int index = 0; index < row_length; ++index) - { - const size_type col = dsp.column_number(row, index); - *cols++ = col; - } - } - - compressed = true; - }; -}; - - - -int -main() -{ - initlog(); - deallog << std::setprecision(3) << std::fixed; - - DynamicSparsityPattern dsp(2); - dsp.add(0, 1); - dsp.add(1, 0); - - SparsityPattern sp_usual; - sp_usual.copy_from(dsp); - deallog << "SparsityPattern:" << std::endl; - sp_usual.print(deallog.get_file_stream()); - - SparsityPatternStandard sp; - sp.copy_from(dsp); - - deallog << "SparsityPatternStandard:" << std::endl; - sp.print(deallog.get_file_stream()); -} diff --git a/tests/sparsity/sparsity_pattern_13.output b/tests/sparsity/sparsity_pattern_13.output deleted file mode 100644 index 0da0dd68e1..0000000000 --- a/tests/sparsity/sparsity_pattern_13.output +++ /dev/null @@ -1,7 +0,0 @@ - -DEAL::SparsityPattern: -[0,0,1] -[1,1,0] -DEAL::SparsityPatternStandard: -[0,1] -[1,0]