From 5213023ca165bfbd420423551c184a3e31886223 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Wed, 27 Feb 2019 22:04:58 +0100 Subject: [PATCH] split SparsityPattern into a base class and derived class with a special treatment for diagonal matrices --- include/deal.II/lac/sparsity_pattern.h | 1041 +++++++++++-------- source/lac/sparsity_pattern.cc | 94 +- tests/serialization/sparsity_pattern.output | 4 +- tests/sparsity/sparsity_pattern_13.cc | 171 +++ tests/sparsity/sparsity_pattern_13.output | 7 + 5 files changed, 843 insertions(+), 474 deletions(-) create mode 100644 tests/sparsity/sparsity_pattern_13.cc create mode 100644 tests/sparsity/sparsity_pattern_13.output diff --git a/include/deal.II/lac/sparsity_pattern.h b/include/deal.II/lac/sparsity_pattern.h index 573c6c2947..e87ca818e2 100644 --- a/include/deal.II/lac/sparsity_pattern.h +++ b/include/deal.II/lac/sparsity_pattern.h @@ -42,6 +42,7 @@ DEAL_II_NAMESPACE_OPEN class SparsityPattern; +class SparsityPatternBase; class DynamicSparsityPattern; class ChunkSparsityPattern; template @@ -139,12 +140,12 @@ namespace SparsityPatternIterators /** * Constructor. */ - Accessor(const SparsityPattern *matrix, const std::size_t linear_index); + Accessor(const SparsityPatternBase *matrix, const std::size_t linear_index); /** * Constructor. Construct the end accessor for the given sparsity pattern. */ - Accessor(const SparsityPattern *matrix); + Accessor(const SparsityPatternBase *matrix); /** * Default constructor creating a dummy accessor. This constructor is here @@ -225,7 +226,7 @@ namespace SparsityPatternIterators /** * The sparsity pattern we operate on accessed. */ - const SparsityPattern *container; + const SparsityPatternBase *container; /** * Index in global sparsity pattern. This index represents the location @@ -267,7 +268,7 @@ namespace SparsityPatternIterators * SparsityPattern class for more information. * * @note This class operates directly on the internal data structures of the - * SparsityPattern class. As a consequence, some operations are cheap and + * SparsityPatternBase 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 @@ -290,14 +291,14 @@ namespace SparsityPatternIterators /** * Type of the stored pointer. */ - using container_pointer_type = SparsityPattern *; + using container_pointer_type = SparsityPatternBase *; /** * 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 SparsityPattern *sp, const std::size_t linear_index); + Iterator(const SparsityPatternBase *sp, const std::size_t linear_index); /** * Constructor. Create an iterator into the sparsity pattern @p sp for @@ -318,6 +319,516 @@ namespace SparsityPatternIterators * It uses the compressed row storage * (CSR) format to store data, and is used as the basis for the + * derived SparsityPattern class and 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. + * + * @author Wolfgang Bangerth, Guido Kanschat and others + */ +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 setup 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 + */ + template + void + save(Archive &ar, const unsigned int version) const; + + /** + * Read the data of this object from a stream for the purpose of + * serialization + */ + template + void + load(Archive &ar, const unsigned int version); + + BOOST_SERIALIZATION_SPLIT_MEMBER() + + // @} + + /** + * @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 compressed row storage + * (CSR) format to store data, and is used as the basis for the * SparseMatrix class. * * The elements of a SparsityPattern, corresponding to the places where @@ -343,19 +854,19 @@ namespace SparsityPatternIterators * * @author Wolfgang Bangerth, Guido Kanschat and others */ -class SparsityPattern : public Subscriptor +class SparsityPattern : public SparsityPatternBase { public: /** * Declare type for container size. */ - using size_type = types::global_dof_index; + using size_type = SparsityPatternBase::size_type; /** * Typedef an iterator class that allows to walk over all nonzero elements * of a sparsity pattern. */ - using const_iterator = SparsityPatternIterators::Iterator; + using const_iterator = SparsityPatternBase::const_iterator; /** * Typedef an iterator class that allows to walk over all nonzero elements @@ -364,8 +875,14 @@ 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 = SparsityPatternIterators::Iterator; + using iterator = SparsityPatternBase::iterator; + /** + * 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; /** * Define a value which is used to indicate that a certain value in the @@ -381,7 +898,7 @@ public: * perform some optimizations, but the actual value of the variable may * change over time. */ - static const size_type invalid_entry = numbers::invalid_size_type; + static const size_type invalid_entry = SparsityPatternBase::invalid_entry; /** * @name Construction and setup Constructors, destructor; functions @@ -495,45 +1012,15 @@ public: SparsityPattern & operator=(const SparsityPattern &); - /** - * 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. - * - * If the number of rows equals the number of columns and the last parameter - * is true, diagonal elements are stored first in each row to allow - * optimized access in relaxation methods of SparseMatrix. - */ - void - reinit(const size_type m, - const size_type n, - const std::vector &row_lengths); - - - /** - * Same as above, but with an ArrayView argument instead. + * 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$. */ - void + virtual void reinit(const size_type m, const size_type n, - const ArrayView &row_lengths); + const ArrayView &row_lengths) override; /** * This function compresses the sparsity structure that this object @@ -649,174 +1136,41 @@ public: copy_from(const SparsityPattern &sp); /** - * Take a full matrix and use its nonzero entries to generate a sparse - * matrix entry pattern for this object. - * - * Previous content of this object is lost, and the sparsity pattern is in - * compressed mode afterwards. - */ - template - void - copy_from(const FullMatrix &matrix); - - /** - * 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); - - /** - * Add several nonzero entries to the specified matrix row. This function - * may only be called for non-compressed sparsity patterns. - * - * If some of the entries already exist, nothing bad happens. - */ - template - void - add_entries(const size_type row, - ForwardIterator begin, - ForwardIterator end, - const bool indices_are_sorted = false); - - // @} - - - - /** - * @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. - * - * Note the discussion in the general documentation of this class about the - * order in which elements are accessed. - */ - 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. - * - * Note also the discussion in the general documentation of this class about - * the order in which elements are accessed. - */ - 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 SparsityPattern &) 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; - - /** - * 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. + * Take a full matrix and use its nonzero entries to generate a sparse + * matrix entry pattern for this object. + * + * Previous content of this object is lost, and the sparsity pattern is in + * compressed mode afterwards. */ - bool - is_compressed() const; + template + void + copy_from(const FullMatrix &matrix); /** - * Return number of rows of this matrix, which equals the dimension of the - * image space. + * Add several nonzero entries to the specified matrix row. This function + * may only be called for non-compressed sparsity patterns. + * + * If some of the entries already exist, nothing bad happens. */ - size_type - n_rows() const; + template + void + add_entries(const size_type row, + ForwardIterator begin, + ForwardIterator end, + const bool indices_are_sorted = false); + + // @} + /** - * Return number of columns of this matrix, which equals the dimension of - * the range space. + * @name Querying information */ - size_type - n_cols() const; - + // @{ /** - * Number of entries in a specific row. + * Test for equality of two SparsityPatterns. */ - unsigned int - row_length(const size_type row) const; + bool + operator==(const SparsityPattern &) const; /** * Return whether this object stores only those entries that have been added @@ -870,55 +1224,12 @@ public: size_type operator()(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; - - /** - * Check if a value at a certain position may be non-zero. - */ - bool - exists(const size_type i, const size_type j) 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; - - /** - * 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; - - // @} /** * @name Input/Output */ // @{ + /** * 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 @@ -948,41 +1259,6 @@ public: void block_read(std::istream &in); - /** - * 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 @@ -1007,32 +1283,6 @@ public: * @addtogroup Exceptions * @{ */ - /** - * 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?)"); - /** - * 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."); - /** - * 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."); /** * Exception */ @@ -1050,85 +1300,6 @@ public: << ", but must be greater than zero."); //@} private: - /** - * 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; - /** * Is special treatment of diagonals enabled? */ @@ -1166,15 +1337,15 @@ private: namespace SparsityPatternIterators { - inline Accessor::Accessor(const SparsityPattern *sparsity_pattern, - const std::size_t i) + inline Accessor::Accessor(const SparsityPatternBase *sparsity_pattern, + const std::size_t i) : container(sparsity_pattern) , linear_index(i) {} - inline Accessor::Accessor(const SparsityPattern *sparsity_pattern) + inline Accessor::Accessor(const SparsityPatternBase *sparsity_pattern) : container(sparsity_pattern) , linear_index(container->rowstart[container->rows]) {} @@ -1274,8 +1445,8 @@ namespace SparsityPatternIterators } - inline Iterator::Iterator(const SparsityPattern *sp, - const std::size_t linear_index) + inline Iterator::Iterator(const SparsityPatternBase *sp, + const std::size_t linear_index) : LinearIndexIterator(Accessor(sp, linear_index)) {} @@ -1289,8 +1460,8 @@ namespace SparsityPatternIterators -inline SparsityPattern::iterator -SparsityPattern::begin() const +inline SparsityPatternBase::iterator +SparsityPatternBase::begin() const { if (n_rows() > 0) return {this, rowstart[0]}; @@ -1300,8 +1471,8 @@ SparsityPattern::begin() const -inline SparsityPattern::iterator -SparsityPattern::end() const +inline SparsityPatternBase::iterator +SparsityPatternBase::end() const { if (n_rows() > 0) return {this, rowstart[rows]}; @@ -1311,8 +1482,8 @@ SparsityPattern::end() const -inline SparsityPattern::iterator -SparsityPattern::begin(const size_type r) const +inline SparsityPatternBase::iterator +SparsityPatternBase::begin(const size_type r) const { Assert(r < n_rows(), ExcIndexRangeType(r, 0, n_rows())); @@ -1321,8 +1492,8 @@ SparsityPattern::begin(const size_type r) const -inline SparsityPattern::iterator -SparsityPattern::end(const size_type r) const +inline SparsityPatternBase::iterator +SparsityPatternBase::end(const size_type r) const { Assert(r < n_rows(), ExcIndexRangeType(r, 0, n_rows())); @@ -1331,16 +1502,16 @@ SparsityPattern::end(const size_type r) const -inline SparsityPattern::size_type -SparsityPattern::n_rows() const +inline SparsityPatternBase::size_type +SparsityPatternBase::n_rows() const { return rows; } -inline SparsityPattern::size_type -SparsityPattern::n_cols() const +inline SparsityPatternBase::size_type +SparsityPatternBase::n_cols() const { return cols; } @@ -1348,7 +1519,7 @@ SparsityPattern::n_cols() const inline bool -SparsityPattern::is_compressed() const +SparsityPatternBase::is_compressed() const { return compressed; } @@ -1364,7 +1535,7 @@ SparsityPattern::stores_only_added_elements() const inline unsigned int -SparsityPattern::row_length(const size_type row) const +SparsityPatternBase::row_length(const size_type row) const { Assert(row < rows, ExcIndexRangeType(row, 0, rows)); return rowstart[row + 1] - rowstart[row]; @@ -1373,8 +1544,8 @@ SparsityPattern::row_length(const size_type row) const inline SparsityPattern::size_type -SparsityPattern::column_number(const size_type row, - const unsigned int index) const +SparsityPatternBase::column_number(const size_type row, + const unsigned int index) const { Assert(row < rows, ExcIndexRangeType(row, 0, rows)); Assert(index < row_length(row), ExcIndexRange(index, 0, row_length(row))); @@ -1385,7 +1556,7 @@ SparsityPattern::column_number(const size_type row, inline std::size_t -SparsityPattern::n_nonzero_elements() const +SparsityPatternBase::n_nonzero_elements() const { Assert(compressed, ExcNotCompressed()); @@ -1400,13 +1571,12 @@ SparsityPattern::n_nonzero_elements() const template inline void -SparsityPattern::save(Archive &ar, const unsigned int) const +SparsityPatternBase::save(Archive &ar, const unsigned int) const { // forward to serialization function in the base class. - ar &static_cast(*this); + ar &boost::serialization::base_object(*this); - ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed - &store_diagonal_first_in_row; + ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed; ar &boost::serialization::make_array(rowstart.get(), max_dim + 1); ar &boost::serialization::make_array(colnums.get(), max_vec_len); @@ -1416,13 +1586,12 @@ SparsityPattern::save(Archive &ar, const unsigned int) const template inline void -SparsityPattern::load(Archive &ar, const unsigned int) +SparsityPatternBase::load(Archive &ar, const unsigned int) { // forward to serialization function in the base class. - ar &static_cast(*this); + ar &boost::serialization::base_object(*this); - ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed - &store_diagonal_first_in_row; + ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed; rowstart = std_cxx14::make_unique(max_dim + 1); colnums = std_cxx14::make_unique(max_vec_len); @@ -1433,15 +1602,36 @@ SparsityPattern::load(Archive &ar, const unsigned int) +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 &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 &store_diagonal_first_in_row; +} + + + inline bool -SparsityPattern::operator==(const SparsityPattern &sp2) const +SparsityPatternBase::operator==(const SparsityPatternBase &sp2) const { // 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 // of the member variables are not yet set correctly) - if (rows != sp2.rows || cols != sp2.cols || compressed != sp2.compressed || - store_diagonal_first_in_row != sp2.store_diagonal_first_in_row) + if (rows != sp2.rows || cols != sp2.cols || compressed != sp2.compressed) return false; for (size_type i = 0; i < rows + 1; ++i) @@ -1457,6 +1647,15 @@ SparsityPattern::operator==(const SparsityPattern &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 diff --git a/source/lac/sparsity_pattern.cc b/source/lac/sparsity_pattern.cc index 2eb565a454..ee331f9ead 100644 --- a/source/lac/sparsity_pattern.cc +++ b/source/lac/sparsity_pattern.cc @@ -36,15 +36,21 @@ DEAL_II_NAMESPACE_OPEN __declspec(selectany) // Weak extern binding due to multiple link error #endif const SparsityPattern::size_type SparsityPattern::invalid_entry; +const SparsityPatternBase::size_type SparsityPatternBase::invalid_entry; - -SparsityPattern::SparsityPattern() +SparsityPatternBase::SparsityPatternBase() : max_dim(0) , max_vec_len(0) , rowstart(nullptr) , colnums(nullptr) , compressed(false) +{} + + + +SparsityPattern::SparsityPattern() + : SparsityPatternBase() , store_diagonal_first_in_row(false) { reinit(0, 0, 0); @@ -53,12 +59,7 @@ SparsityPattern::SparsityPattern() SparsityPattern::SparsityPattern(const SparsityPattern &s) - : Subscriptor() - , max_dim(0) - , max_vec_len(0) - , rowstart(nullptr) - , colnums(nullptr) - , compressed(false) + : SparsityPatternBase() , store_diagonal_first_in_row(false) { (void)s; @@ -76,11 +77,7 @@ SparsityPattern::SparsityPattern(const SparsityPattern &s) SparsityPattern::SparsityPattern(const size_type m, const size_type n, const unsigned int max_per_row) - : max_dim(0) - , max_vec_len(0) - , rowstart(nullptr) - , colnums(nullptr) - , compressed(false) + : SparsityPatternBase() , store_diagonal_first_in_row(m == n) { reinit(m, n, max_per_row); @@ -91,10 +88,7 @@ SparsityPattern::SparsityPattern(const size_type m, SparsityPattern::SparsityPattern(const size_type m, const size_type n, const std::vector &row_lengths) - : max_dim(0) - , max_vec_len(0) - , rowstart(nullptr) - , colnums(nullptr) + : SparsityPatternBase() , store_diagonal_first_in_row(m == n) { reinit(m, n, row_lengths); @@ -104,10 +98,7 @@ SparsityPattern::SparsityPattern(const size_type m, SparsityPattern::SparsityPattern(const size_type m, const unsigned int max_per_row) - : max_dim(0) - , max_vec_len(0) - , rowstart(nullptr) - , colnums(nullptr) + : SparsityPatternBase() { reinit(m, m, max_per_row); } @@ -116,10 +107,7 @@ SparsityPattern::SparsityPattern(const size_type m, SparsityPattern::SparsityPattern(const size_type m, const std::vector &row_lengths) - : max_dim(0) - , max_vec_len(0) - , rowstart(nullptr) - , colnums(nullptr) + : SparsityPatternBase() { reinit(m, m, row_lengths); } @@ -129,10 +117,7 @@ SparsityPattern::SparsityPattern(const size_type m, SparsityPattern::SparsityPattern(const SparsityPattern &original, const unsigned int max_per_row, const size_type extra_off_diagonals) - : max_dim(0) - , max_vec_len(0) - , rowstart(nullptr) - , colnums(nullptr) + : SparsityPattern() { Assert(original.rows == original.cols, ExcNotQuadratic()); Assert(original.is_compressed(), ExcNotCompressed()); @@ -228,9 +213,9 @@ SparsityPattern::operator=(const SparsityPattern &s) void -SparsityPattern::reinit(const size_type m, - const size_type n, - const unsigned int max_per_row) +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); @@ -601,9 +586,9 @@ SparsityPattern::copy_from(const FullMatrix &matrix) void -SparsityPattern::reinit(const size_type m, - const size_type n, - const std::vector &row_lengths) +SparsityPatternBase::reinit(const size_type m, + const size_type n, + const std::vector &row_lengths) { reinit(m, n, make_slice(row_lengths)); } @@ -611,7 +596,7 @@ SparsityPattern::reinit(const size_type m, bool -SparsityPattern::empty() const +SparsityPatternBase::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 @@ -633,8 +618,8 @@ SparsityPattern::empty() const -SparsityPattern::size_type -SparsityPattern::max_entries_per_row() const +SparsityPatternBase::size_type +SparsityPatternBase::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 @@ -691,7 +676,7 @@ SparsityPattern::operator()(const size_type i, const size_type j) const void -SparsityPattern::add(const size_type i, const size_type j) +SparsityPatternBase::add(const size_type i, const size_type j) { Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); Assert(i < rows, ExcIndexRange(i, 0, rows)); @@ -768,7 +753,7 @@ SparsityPattern::add_entries(const size_type row, bool -SparsityPattern::exists(const size_type i, const size_type j) const +SparsityPatternBase::exists(const size_type i, const size_type j) const { Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); Assert(i < rows, ExcIndexRange(i, 0, rows)); @@ -785,8 +770,8 @@ SparsityPattern::exists(const size_type i, const size_type j) const -SparsityPattern::size_type -SparsityPattern::row_position(const size_type i, const size_type j) const +SparsityPatternBase::size_type +SparsityPatternBase::row_position(const size_type i, const size_type j) const { Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); Assert(i < rows, ExcIndexRange(i, 0, rows)); @@ -803,8 +788,8 @@ SparsityPattern::row_position(const size_type i, const size_type j) const -std::pair -SparsityPattern::matrix_position(const std::size_t global_index) const +std::pair +SparsityPatternBase::matrix_position(const std::size_t global_index) const { Assert(compressed == true, ExcNotCompressed()); Assert(global_index < n_nonzero_elements(), @@ -829,7 +814,7 @@ SparsityPattern::matrix_position(const std::size_t global_index) const void -SparsityPattern::symmetrize() +SparsityPatternBase::symmetrize() { Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); Assert(compressed == false, ExcMatrixIsCompressed()); @@ -863,7 +848,7 @@ SparsityPattern::symmetrize() void -SparsityPattern::print(std::ostream &out) const +SparsityPatternBase::print(std::ostream &out) const { Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); @@ -884,7 +869,7 @@ SparsityPattern::print(std::ostream &out) const void -SparsityPattern::print_gnuplot(std::ostream &out) const +SparsityPatternBase::print_gnuplot(std::ostream &out) const { Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); @@ -902,7 +887,7 @@ SparsityPattern::print_gnuplot(std::ostream &out) const } void -SparsityPattern::print_svg(std::ostream &out) const +SparsityPatternBase::print_svg(std::ostream &out) const { unsigned int m = this->n_rows(); unsigned int n = this->n_cols(); @@ -935,8 +920,8 @@ SparsityPattern::print_svg(std::ostream &out) const -SparsityPattern::size_type -SparsityPattern::bandwidth() const +SparsityPatternBase::size_type +SparsityPatternBase::bandwidth() const { Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject()); size_type b = 0; @@ -1019,13 +1004,20 @@ SparsityPattern::block_read(std::istream &in) std::size_t -SparsityPattern::memory_consumption() const +SparsityPatternBase::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(); +} + + // explicit instantiations template void diff --git a/tests/serialization/sparsity_pattern.output b/tests/serialization/sparsity_pattern.output index cf50082240..abc3bde154 100644 --- a/tests/serialization/sparsity_pattern.output +++ b/tests/serialization/sparsity_pattern.output @@ -1,6 +1,6 @@ -DEAL::0 0 0 0 16 16 16 64 5 1 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 +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 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 +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::OK diff --git a/tests/sparsity/sparsity_pattern_13.cc b/tests/sparsity/sparsity_pattern_13.cc new file mode 100644 index 0000000000..c1b0a8f1ff --- /dev/null +++ b/tests/sparsity/sparsity_pattern_13.cc @@ -0,0 +1,171 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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_cxx14::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_cxx14::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_cxx14::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() +{ + std::ofstream logfile("output"); + logfile.setf(std::ios::fixed); + deallog << std::setprecision(3); + deallog.attach(logfile); + + 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 new file mode 100644 index 0000000000..0da0dd68e1 --- /dev/null +++ b/tests/sparsity/sparsity_pattern_13.output @@ -0,0 +1,7 @@ + +DEAL::SparsityPattern: +[0,0,1] +[1,1,0] +DEAL::SparsityPatternStandard: +[0,1] +[1,0] -- 2.39.5