]> https://gitweb.dealii.org/ - dealii.git/commitdiff
split SparsityPattern into a base class and derived class with a special treatment... 7766/head
authorDenis Davydov <davydden@gmail.com>
Wed, 27 Feb 2019 21:04:58 +0000 (22:04 +0100)
committerDenis Davydov <davydden@gmail.com>
Wed, 20 Mar 2019 21:43:37 +0000 (22:43 +0100)
include/deal.II/lac/sparsity_pattern.h
source/lac/sparsity_pattern.cc
tests/serialization/sparsity_pattern.output
tests/sparsity/sparsity_pattern_13.cc [new file with mode: 0644]
tests/sparsity/sparsity_pattern_13.output [new file with mode: 0644]

index 573c6c29475bbe9c6fa533bada2c941c0d76f29e..e87ca818e207142a00b922efbfb2abd6afab609a 100644 (file)
@@ -42,6 +42,7 @@
 DEAL_II_NAMESPACE_OPEN
 
 class SparsityPattern;
+class SparsityPatternBase;
 class DynamicSparsityPattern;
 class ChunkSparsityPattern;
 template <typename number>
@@ -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 <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 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 <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
+   */
+  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
+   */
+  template <class Archive>
+  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
+   * <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
+ * 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
  * 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 <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.
-   *
-   * 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<unsigned int> &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<const unsigned int> &row_lengths);
+         const ArrayView<const unsigned int> &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 <typename number>
-  void
-  copy_from(const FullMatrix<number> &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 <typename ForwardIterator>
-  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 <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.
-   *
-   * 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 <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 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 <typename number>
+  void
+  copy_from(const FullMatrix<number> &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 <typename ForwardIterator>
+  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 <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;
-
-  /**
-   * 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
-   * <tt>j</tt>.
-   */
-  size_type
-  row_position(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;
-
-
   // @}
   /**
    * @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 <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
@@ -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
-   * <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;
-
   /**
    * 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<Iterator, Accessor>(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<size_type>(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<size_type>(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<size_type>(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<size_type>(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 <class Archive>
 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<const Subscriptor &>(*this);
+  ar &boost::serialization::base_object<const Subscriptor>(*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 <class Archive>
 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<Subscriptor &>(*this);
+  ar &boost::serialization::base_object<Subscriptor>(*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<std::size_t[]>(max_dim + 1);
   colnums  = std_cxx14::make_unique<size_type[]>(max_vec_len);
@@ -1433,15 +1602,36 @@ SparsityPattern::load(Archive &ar, const unsigned int)
 
 
 
+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 SparsityPatternBase>(*this);
+  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<SparsityPatternBase>(*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<const SparsityPatternBase &>(*this) == sp2) &&
+         (store_diagonal_first_in_row == sp2.store_diagonal_first_in_row);
+}
+
+
+
 namespace internal
 {
   namespace SparsityPatternTools
index 2eb565a4542dcf09635dee41ee2ed21bfcf91fb9..ee331f9ead2766877db70928f0d92e993cf21791 100644 (file)
@@ -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<unsigned int> &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<unsigned int> &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<unsigned int> row_lengths(m, max_per_row);
@@ -601,9 +586,9 @@ SparsityPattern::copy_from(const FullMatrix<number> &matrix)
 
 
 void
-SparsityPattern::reinit(const size_type                  m,
-                        const size_type                  n,
-                        const std::vector<unsigned int> &row_lengths)
+SparsityPatternBase::reinit(const size_type                  m,
+                            const size_type                  n,
+                            const std::vector<unsigned int> &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::size_type, SparsityPattern::size_type>
-SparsityPattern::matrix_position(const std::size_t global_index) const
+std::pair<SparsityPatternBase::size_type, SparsityPatternBase::size_type>
+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
index cf500822409cdda14ecce3ab78d5a8ccce245cf3..abc3bde154b59c23037c62f3d8ee7605e9c9f271 100644 (file)
@@ -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 (file)
index 0000000..c1b0a8f
--- /dev/null
@@ -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<const unsigned int> &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<size_type>(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<size_type[]>(max_vec_len);
+      }
+
+    max_row_length =
+      (row_lengths.size() == 0 ?
+         0 :
+         std::min(static_cast<size_type>(
+                    *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<std::size_t[]>(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<size_type[]>(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<size_type>(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<unsigned int> 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 (file)
index 0000000..0da0dd6
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL::SparsityPattern:
+[0,0,1]
+[1,1,0]
+DEAL::SparsityPatternStandard:
+[0,1]
+[1,0]

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.