// Forward declarations
#ifndef DOXYGEN
class SparsityPattern;
-class SparsityPatternBase;
class DynamicSparsityPattern;
class ChunkSparsityPattern;
template <typename number>
/**
* 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
/**
* The sparsity pattern we operate on accessed.
*/
- const SparsityPatternBase *container;
+ const SparsityPattern *container;
/**
* Index in global sparsity pattern. This index represents the location
* 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
/**
* 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
};
} // namespace SparsityPatternIterators
-
-
-/**
- * A class that can store which elements of a matrix are nonzero (or, to be
- * precise, <i>may</i> 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 <a
- * href="https://en.wikipedia.org/wiki/Sparse_matrix">compressed row storage
- * (CSR)</a> 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 <tt>m*n==0</tt> 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<unsigned int> &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<const unsigned int> &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 <tt>r</tt>.
- *
- * Note that if the given row is empty, i.e. does not contain any nonzero
- * entries, then the iterator returned by this function equals
- * <tt>end(r)</tt>. 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 <tt>r</tt>. 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
- * <tt>index</tt>th entry in <tt>row</tt>. Note that if diagonal elements
- * are optimized, the first element in each row is the diagonal element,
- * i.e. <tt>column_number(row,0)==row</tt>.
- *
- * If the sparsity pattern is already compressed, then (except for the
- * diagonal element), the entries are sorted by columns, i.e.
- * <tt>column_number(row,i)</tt> <tt><</tt> <tt>column_number(row,i+1)</tt>.
- */
- 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
- * <tt>j</tt>.
- */
- 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 <tt>N</tt> is the number of rows of this matrix, then the complexity
- * of this function is <i>log(N)</i>.
- */
- std::pair<size_type, size_type>
- 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 <tt>[i,j1,j2,j3,...]</tt>. <i>i</i> is the row number and
- * <i>jn</i> are the allocated columns in this row.
- */
- void
- print(std::ostream &out) const;
-
- /**
- * Print the sparsity of the matrix in a format that <tt>gnuplot</tt>
- * understands and which can be used to plot the sparsity pattern in a
- * graphical way. The format consists of pairs <tt>i j</tt> 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 <tt>(0,0)</tt>
- * 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 <tt>plot</tt> 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 <class Archive>
- 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 <class Archive>
- 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 <class Archive>
- 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
- * <tt>row</tt>=#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 <tt>for
- * (i=rowstart[k]; i<rowstart[k+1]; ++i)</tt> 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<std::size_t[]> 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
- * <i>r</i> are stored within the index range
- * #rowstart[<i>r</i>]...#rowstart[<i>r+1</i>]. Therefore to find out
- * whether a given element (<i>r,c</i>) exists, we have to check whether the
- * column number <i>c</i> exists in the above-mentioned range within this
- * array. If it exists, say at position <i>p</i> within this array, the
- * value of the respective element in the sparse matrix will also be at
- * position <i>p</i> 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<size_type[]> colnums;
-
- /**
- * Store whether the compress() function was called for this object.
- */
- bool compressed;
-
- // Make all sparse matrices friends of this class.
- template <typename number>
- friend class SparseMatrix;
- template <typename number>
- friend class SparseLUDecomposition;
- template <typename number>
- friend class SparseILU;
- template <typename number>
- 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 <a
* SparsityPatternIterators::Accessor, like <tt>it->column()</tt> and
* <tt>it->row()</tt>.
*/
-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
* 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 <tt>r</tt>.
+ *
+ * Note that if the given row is empty, i.e. does not contain any nonzero
+ * entries, then the iterator returned by this function equals
+ * <tt>end(r)</tt>. 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 <tt>r</tt>. 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
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
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 <tt>m*n==0</tt> 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<const unsigned int> &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<unsigned int> &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<const unsigned int> &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
void
copy_from(const FullMatrix<number> &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.
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();
+
/**
* @}
*/
* @{
*/
+ /**
+ * 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 <tt>N</tt> is the number of rows of this matrix, then the complexity
+ * of this function is <i>log(N)</i>.
+ */
+ std::pair<size_type, size_type>
+ 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.
*/
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.
size_type
operator()(const size_type i, const size_type j) const;
+ /**
+ * Access to column number field. Return the column number of the
+ * <tt>index</tt>th entry in <tt>row</tt>. Note that if diagonal elements
+ * are optimized, the first element in each row is the diagonal element,
+ * i.e. <tt>column_number(row,0)==row</tt>.
+ *
+ * If the sparsity pattern is already compressed, then (except for the
+ * diagonal element), the entries are sorted by columns, i.e.
+ * <tt>column_number(row,i)</tt> <tt><</tt> <tt>column_number(row,i+1)</tt>.
+ */
+ 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
+ * <tt>j</tt>.
+ */
+ size_type
+ row_position(const size_type i, const size_type j) const;
+
/**
* @}
*/
* @{
*/
+ /**
+ * Print the sparsity of the matrix. The output consists of one line per row
+ * of the format <tt>[i,j1,j2,j3,...]</tt>. <i>i</i> is the row number and
+ * <i>jn</i> are the allocated columns in this row.
+ */
+ void
+ print(std::ostream &out) const;
+
+ /**
+ * Print the sparsity of the matrix in a format that <tt>gnuplot</tt>
+ * understands and which can be used to plot the sparsity pattern in a
+ * graphical way. The format consists of pairs <tt>i j</tt> 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 <tt>(0,0)</tt>
+ * 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 <tt>plot</tt> 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
* Read the data of this object from a stream for the purpose of
* serialization
*/
- template <class Archive>
- void
- load(Archive &ar, const unsigned int version);
+ template <class Archive>
+ 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 <class Archive>
+ 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 <class Archive>
- 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
+ * <tt>row</tt>=#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 <tt>for
+ * (i=rowstart[k]; i<rowstart[k+1]; ++i)</tt> 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<std::size_t[]> 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
+ * <i>r</i> are stored within the index range
+ * #rowstart[<i>r</i>]...#rowstart[<i>r+1</i>]. Therefore to find out
+ * whether a given element (<i>r,c</i>) exists, we have to check whether the
+ * column number <i>c</i> exists in the above-mentioned range within this
+ * array. If it exists, say at position <i>p</i> within this array, the
+ * value of the respective element in the sparse matrix will also be at
+ * position <i>p</i> 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<size_type[]> 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 <typename number>
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])
{}
}
- 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<Iterator, Accessor>(Accessor(sp, linear_index))
{}
-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;
}
inline bool
-SparsityPatternBase::is_compressed() const
+SparsityPattern::is_compressed() const
{
return compressed;
}
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];
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));
-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 <class Archive>
-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<const Subscriptor>(*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 <class Archive>
-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<Subscriptor>(*this);
-
- ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed;
-
- if (max_dim != 0)
- {
- rowstart = std::make_unique<std::size_t[]>(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<size_type[]>(max_vec_len);
- ar &boost::serialization::make_array(colnums.get(), max_vec_len);
- }
- else
- colnums.reset();
+ return {this, rowstart[r]};
}
-template <class Archive>
-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<const SparsityPatternBase>(*this);
- ar &store_diagonal_first_in_row;
-}
-
-
+ AssertIndexRange(r, n_rows());
-template <class Archive>
-inline void
-SparsityPattern::load(Archive &ar, const unsigned int)
-{
- // forward to serialization function in the base class.
- ar &boost::serialization::base_object<SparsityPatternBase>(*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
-inline bool
-SparsityPattern::operator==(const SparsityPattern &sp2) const
-{
- return (static_cast<const SparsityPatternBase &>(*this) == sp2) &&
- (store_diagonal_first_in_row == sp2.store_diagonal_first_in_row);
-}
-
-
-
namespace internal
{
namespace SparsityPatternTools
}
+
+template <class Archive>
+inline void
+SparsityPattern::save(Archive &ar, const unsigned int) const
+{
+ // forward to serialization function in the base class.
+ ar &boost::serialization::base_object<const Subscriptor>(*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 <class Archive>
+inline void
+SparsityPattern::load(Archive &ar, const unsigned int)
+{
+ // forward to serialization function in the base class.
+ ar &boost::serialization::base_object<Subscriptor>(*this);
+
+ ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed;
+
+ if (max_dim != 0)
+ {
+ rowstart = std::make_unique<std::size_t[]>(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<size_type[]>(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
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);
}
SparsityPattern::SparsityPattern(const SparsityPattern &s)
- : SparsityPatternBase()
- , store_diagonal_first_in_row(false)
+ : SparsityPattern()
{
(void)s;
Assert(s.empty(),
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);
}
SparsityPattern::SparsityPattern(const size_type m,
const size_type n,
const std::vector<unsigned int> &row_lengths)
- : SparsityPatternBase()
- , store_diagonal_first_in_row(m == n)
+ : SparsityPattern()
{
reinit(m, n, row_lengths);
}
SparsityPattern::SparsityPattern(const size_type m,
const unsigned int max_per_row)
- : SparsityPatternBase()
+ : SparsityPattern()
{
reinit(m, m, max_per_row);
}
SparsityPattern::SparsityPattern(const size_type m,
const std::vector<unsigned int> &row_lengths)
- : SparsityPatternBase()
+ : SparsityPattern()
{
reinit(m, m, row_lengths);
}
-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<unsigned int> row_lengths(m, max_per_row);
- reinit(m, n, row_lengths);
-}
-
-
-
void
SparsityPattern::reinit(const size_type m,
const size_type n,
+void
+SparsityPattern::reinit(const size_type m,
+ const size_type n,
+ const std::vector<unsigned int> &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<unsigned int> row_lengths(m, max_per_row);
+ reinit(m, n, row_lengths);
+}
+
+
+
void
SparsityPattern::compress()
{
}
-void
-SparsityPatternBase::reinit(const size_type m,
- const size_type n,
- const std::vector<unsigned int> &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
-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::size_type, SparsityPattern::size_type>
+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<size_type>(
+ std::abs(static_cast<int>(i - colnums[j]))) > b)
+ b = std::abs(static_cast<signed int>(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
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);
}
-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::size_type, SparsityPatternBase::size_type>
-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());
+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());
void
-SparsityPatternBase::print_gnuplot(std::ostream &out) const
+SparsityPattern::print_gnuplot(std::ostream &out) const
{
Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
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();
-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<size_type>(
- std::abs(static_cast<int>(i - colnums[j]))) > b)
- b = std::abs(static_cast<signed int>(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
{
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
projects / dealii.git / commitdiff