Reference documentation for deal.II version GIT 00e6fe71c9 2023-06-04 19:35:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
dynamic_sparsity_pattern.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_dynamic_sparsity_pattern_h
17 #define dealii_dynamic_sparsity_pattern_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/index_set.h>
24 #include <deal.II/base/utilities.h>
25 
26 #include <deal.II/lac/exceptions.h>
28 
29 #include <algorithm>
30 #include <iostream>
31 #include <vector>
32 
34 
35 // Forward declaration
36 #ifndef DOXYGEN
38 #endif
39 
50 {
51  // forward declaration
52  class Iterator;
53 
58 
67  class Accessor
68  {
69  public:
74  const size_type row,
75  const unsigned int index_within_row);
76 
81 
87  Accessor();
88 
92  size_type
93  row() const;
94 
98  size_type
99  index() const;
100 
104  size_type
105  column() const;
106 
110  bool
111  operator==(const Accessor &) const;
112 
120  bool
121  operator<(const Accessor &) const;
122 
123  protected:
125  "The instance of this class was initialized"
126  " without DynamicSparsityPattern object, which"
127  " means that it is a dummy accessor that can"
128  " not do any operations.");
129 
134 
139 
144  std::vector<size_type>::const_iterator current_entry;
145 
152  std::vector<size_type>::const_iterator end_of_row;
153 
157  void
158  advance();
159 
160  // Grant access to iterator class.
161  friend class Iterator;
162  };
163 
164 
165 
190  class Iterator
191  {
192  public:
199  const size_type row,
200  const unsigned int index_within_row);
201 
206  Iterator(const DynamicSparsityPattern *sp);
207 
213  Iterator() = default;
214 
218  Iterator &
219  operator++();
220 
224  Iterator
225  operator++(int);
226 
230  const Accessor &
231  operator*() const;
232 
236  const Accessor *
237  operator->() const;
238 
242  bool
243  operator==(const Iterator &) const;
244 
248  bool
249  operator!=(const Iterator &) const;
250 
258  bool
259  operator<(const Iterator &) const;
260 
267  int
268  operator-(const Iterator &p) const;
269 
270  private:
275  };
276 } // namespace DynamicSparsityPatternIterators
277 
278 
323 {
324 public:
329 
338 
344 
351 
362 
370  const size_type n,
371  const IndexSet &rowset = IndexSet());
372 
378  DynamicSparsityPattern(const IndexSet &indexset);
379 
384 
392 
399  void
400  reinit(const size_type m,
401  const size_type n,
402  const IndexSet &rowset = IndexSet());
403 
409  void
410  compress();
411 
416  bool
417  empty() const;
418 
423  size_type
424  max_entries_per_row() const;
425 
429  void
430  add(const size_type i, const size_type j);
431 
436  template <typename ForwardIterator>
437  void
438  add_entries(const size_type row,
439  ForwardIterator begin,
440  ForwardIterator end,
441  const bool indices_are_unique_and_sorted = false);
442 
443  virtual void
444  add_row_entries(const size_type & row,
445  const ArrayView<const size_type> &columns,
446  const bool indices_are_sorted = false) override;
447 
449 
453  bool
454  exists(const size_type i, const size_type j) const;
455 
463  get_view(const IndexSet &rows) const;
464 
472  void
473  symmetrize();
474 
479  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
480  void
481  compute_mmult_pattern(const SparsityPatternTypeLeft & left,
482  const SparsityPatternTypeRight &right);
483 
488  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
489  void
490  compute_Tmmult_pattern(const SparsityPatternTypeLeft & left,
491  const SparsityPatternTypeRight &right);
492 
498  void
499  print(std::ostream &out) const;
500 
514  void
515  print_gnuplot(std::ostream &out) const;
516 
521  size_type
522  row_length(const size_type row) const;
523 
527  void
528  clear_row(const size_type row);
529 
534  size_type
535  column_number(const size_type row, const size_type index) const;
536 
542  size_type
543  column_index(const size_type row, const size_type col) const;
544 
563  iterator
564  begin() const;
565 
569  iterator
570  end() const;
571 
588  iterator
589  begin(const size_type r) const;
590 
599  iterator
600  end(const size_type r) const;
601 
611  size_type
612  bandwidth() const;
613 
618  size_type
619  n_nonzero_elements() const;
620 
626  const IndexSet &
627  row_index_set() const;
628 
635  IndexSet
636  nonempty_cols() const;
637 
644  IndexSet
645  nonempty_rows() const;
646 
657  static bool
659 
664  size_type
665  memory_consumption() const;
666 
667 private:
672 
678 
679 
686  struct Line
687  {
688  public:
693  std::vector<size_type> entries;
694 
698  void
699  add(const size_type col_num);
700 
704  template <typename ForwardIterator>
705  void
706  add_entries(ForwardIterator begin,
707  ForwardIterator end,
708  const bool indices_are_sorted);
709 
713  size_type
714  memory_consumption() const;
715  };
716 
717 
721  std::vector<Line> lines;
722 
723  // make the accessor class a friend
725 };
726 
728 /*---------------------- Inline functions -----------------------------------*/
729 
730 
732 {
733  inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern,
734  const size_type row,
735  const unsigned int index_within_row)
736  : sparsity_pattern(sparsity_pattern)
737  , current_row(row)
738  , current_entry(
739  ((sparsity_pattern->rowset.size() == 0) ?
740  sparsity_pattern->lines[current_row].entries.begin() :
741  sparsity_pattern
742  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
743  .entries.begin()) +
744  index_within_row)
745  , end_of_row(
746  (sparsity_pattern->rowset.size() == 0) ?
747  sparsity_pattern->lines[current_row].entries.end() :
748  sparsity_pattern
749  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
750  .entries.end())
751  {
753  Assert((sparsity_pattern->rowset.size() == 0) ||
755  ExcMessage("You can't create an iterator into a "
756  "DynamicSparsityPattern's row that is not "
757  "actually stored by that sparsity pattern "
758  "based on the IndexSet argument to it."));
760  index_within_row,
761  ((sparsity_pattern->rowset.size() == 0) ?
762  sparsity_pattern->lines[current_row].entries.size() :
765  .entries.size()));
766  }
767 
768 
769  inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern)
770  : sparsity_pattern(sparsity_pattern)
771  , current_row(numbers::invalid_size_type)
772  , current_entry()
773  , end_of_row()
774  {}
775 
776 
777 
779  : sparsity_pattern(nullptr)
780  , current_row(numbers::invalid_size_type)
781  , current_entry()
782  , end_of_row()
783  {}
784 
785 
786  inline size_type
788  {
789  Assert(sparsity_pattern != nullptr, DummyAccessor());
790  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
791 
792  return current_row;
793  }
794 
795 
796  inline size_type
798  {
799  Assert(sparsity_pattern != nullptr, DummyAccessor());
800  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
801 
802  return *current_entry;
803  }
804 
805 
806  inline size_type
808  {
809  Assert(sparsity_pattern != nullptr, DummyAccessor());
810  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
811 
812  return (current_entry -
813  ((sparsity_pattern->rowset.size() == 0) ?
814  sparsity_pattern->lines[current_row].entries.begin() :
817  .entries.begin()));
818  }
819 
820 
821 
822  inline bool
823  Accessor::operator==(const Accessor &other) const
824  {
825  Assert(sparsity_pattern != nullptr, DummyAccessor());
826  Assert(other.sparsity_pattern != nullptr, DummyAccessor());
827  // compare the sparsity pattern the iterator points into, the
828  // current row, and the location within this row. ignore the
829  // latter if the row is past-the-end because in that case the
830  // current_entry field may not point to a deterministic location
831  return (sparsity_pattern == other.sparsity_pattern &&
832  current_row == other.current_row &&
834  (current_entry == other.current_entry)));
835  }
836 
837 
838 
839  inline bool
840  Accessor::operator<(const Accessor &other) const
841  {
842  Assert(sparsity_pattern != nullptr, DummyAccessor());
843  Assert(other.sparsity_pattern != nullptr, DummyAccessor());
845 
846  // if *this is past-the-end, then it is less than no one
848  return (false);
849  // now *this should be an valid value
850  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
851 
852  // if other is past-the-end
854  return (true);
855  // now other should be an valid value
857 
858  // both iterators are not one-past-the-end
859  return ((current_row < other.current_row) ||
860  ((current_row == other.current_row) &&
861  (current_entry < other.current_entry)));
862  }
863 
864 
865  inline void
867  {
868  Assert(sparsity_pattern != nullptr, DummyAccessor());
869  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
870 
871  // move to the next element in this row
872  ++current_entry;
873 
874  // if this moves us beyond the end of the row, go to the next row
875  // if possible, or set the iterator to an invalid state if not.
876  //
877  // going to the next row is a bit complicated because we may have
878  // to skip over empty rows, and because we also have to avoid rows
879  // that aren't listed in a possibly passed IndexSet argument of
880  // the sparsity pattern. consequently, rather than trying to
881  // duplicate code here, just call the begin() function of the
882  // sparsity pattern itself
883  if (current_entry == end_of_row)
884  {
885  if (current_row + 1 < sparsity_pattern->n_rows())
886  *this = *sparsity_pattern->begin(current_row + 1);
887  else
888  *this = Accessor(sparsity_pattern); // invalid object
889  }
890  }
891 
892 
893 
894  inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern,
895  const size_type row,
896  const unsigned int index_within_row)
897  : accessor(sparsity_pattern, row, index_within_row)
898  {}
899 
900 
901 
902  inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern)
903  : accessor(sparsity_pattern)
904  {}
905 
906 
907 
908  inline Iterator &
910  {
911  accessor.advance();
912  return *this;
913  }
914 
915 
916 
917  inline Iterator
919  {
920  const Iterator iter = *this;
921  accessor.advance();
922  return iter;
923  }
924 
925 
926 
927  inline const Accessor &
929  {
930  return accessor;
931  }
932 
933 
934 
935  inline const Accessor *
937  {
938  return &accessor;
939  }
940 
941 
942  inline bool
943  Iterator::operator==(const Iterator &other) const
944  {
945  return (accessor == other.accessor);
946  }
947 
948 
949 
950  inline bool
951  Iterator::operator!=(const Iterator &other) const
952  {
953  return !(*this == other);
954  }
955 
956 
957  inline bool
958  Iterator::operator<(const Iterator &other) const
959  {
960  return accessor < other.accessor;
961  }
962 
963 
964  inline int
965  Iterator::operator-(const Iterator &other) const
966  {
967  (void)other;
969  ExcInternalError());
970  Assert(false, ExcNotImplemented());
971 
972  return 0;
973  }
974 } // namespace DynamicSparsityPatternIterators
975 
976 
977 inline void
979 {
980  // first check the last element (or if line is still empty)
981  if ((entries.size() == 0) || (entries.back() < j))
982  {
983  entries.push_back(j);
984  return;
985  }
986 
987  // do a binary search to find the place where to insert:
988  std::vector<size_type>::iterator it =
989  Utilities::lower_bound(entries.begin(), entries.end(), j);
990 
991  // If this entry is a duplicate, exit immediately
992  if (*it == j)
993  return;
994 
995  // Insert at the right place in the vector. Vector grows automatically to
996  // fit elements. Always doubles its size.
997  entries.insert(it, j);
998 }
999 
1000 
1001 
1002 inline void
1004 {
1005  AssertIndexRange(i, n_rows());
1006  AssertIndexRange(j, n_cols());
1007 
1008  if (rowset.size() > 0 && !rowset.is_element(i))
1009  return;
1010 
1011  have_entries = true;
1012 
1013  const size_type rowindex =
1014  rowset.size() == 0 ? i : rowset.index_within_set(i);
1015  lines[rowindex].add(j);
1016 }
1017 
1018 
1019 
1020 template <typename ForwardIterator>
1021 inline void
1023  ForwardIterator begin,
1024  ForwardIterator end,
1025  const bool indices_are_sorted)
1026 {
1027  AssertIndexRange(row, rows);
1028 
1029  if (rowset.size() > 0 && !rowset.is_element(row))
1030  return;
1031 
1032  if (!have_entries && begin < end)
1033  have_entries = true;
1034 
1035  const size_type rowindex =
1036  rowset.size() == 0 ? row : rowset.index_within_set(row);
1037  lines[rowindex].add_entries(begin, end, indices_are_sorted);
1038 }
1039 
1040 
1041 
1044 {
1045  AssertIndexRange(row, n_rows());
1046 
1047  if (!have_entries)
1048  return 0;
1049 
1050  if (rowset.size() > 0 && !rowset.is_element(row))
1051  return 0;
1052 
1053  const size_type rowindex =
1054  rowset.size() == 0 ? row : rowset.index_within_set(row);
1055  return lines[rowindex].entries.size();
1056 }
1057 
1058 
1059 
1062  const size_type index) const
1063 {
1064  AssertIndexRange(row, n_rows());
1065  Assert(rowset.size() == 0 || rowset.is_element(row), ExcInternalError());
1066 
1067  const size_type local_row =
1068  rowset.size() != 0u ? rowset.index_within_set(row) : row;
1069  AssertIndexRange(index, lines[local_row].entries.size());
1070  return lines[local_row].entries[index];
1071 }
1072 
1073 
1074 
1077 {
1078  if (n_rows() > 0)
1079  return begin(0);
1080  else
1081  return end();
1082 }
1083 
1084 
1087 {
1088  return {this};
1089 }
1090 
1091 
1092 
1095 {
1096  AssertIndexRange(r, n_rows());
1097 
1098  if (!have_entries)
1099  return {this};
1100 
1101  if (rowset.size() > 0)
1102  {
1103  // We have an IndexSet that describes the locally owned set. For
1104  // performance reasons we need to make sure that we don't do a
1105  // linear search over 0..n_rows(). Instead, find the first entry
1106  // >= row r in the locally owned set (this is done in log
1107  // n_ranges time inside at()). From there, we move forward until
1108  // we find a non-empty row. By iterating over the IndexSet instead
1109  // of incrementing the row index, we potentially skip over entries
1110  // not in the rowset.
1112  if (it == rowset.end())
1113  return end(); // we don't own any row between r and the end
1114 
1115  // Instead of using row_length(*it)==0 in the while loop below,
1116  // which involves an expensive index_within_set() call, we
1117  // look at the lines vector directly. This works, because we are
1118  // walking over this vector entry by entry anyways.
1119  size_type rowindex = rowset.index_within_set(*it);
1120 
1121  while (it != rowset.end() && lines[rowindex].entries.size() == 0)
1122  {
1123  ++it;
1124  ++rowindex;
1125  }
1126 
1127  if (it == rowset.end())
1128  return end();
1129  else
1130  return {this, *it, 0};
1131  }
1132 
1133  // Without an index set we have to do a linear search starting at
1134  // row r until we find a non-empty one. We will check the lines vector
1135  // directly instead of going through the slower row_length() function
1136  size_type row = r;
1137 
1138  while (row < n_rows() && lines[row].entries.size() == 0)
1139  {
1140  ++row;
1141  }
1142 
1143  if (row == n_rows())
1144  return {this};
1145  else
1146  return {this, row, 0};
1147 }
1148 
1149 
1150 
1153 {
1154  AssertIndexRange(r, n_rows());
1155 
1156  const size_type row = r + 1;
1157  if (row == n_rows())
1158  return {this};
1159  else
1160  return begin(row);
1161 }
1162 
1163 
1164 
1165 inline const IndexSet &
1167 {
1168  return rowset;
1169 }
1170 
1171 
1172 
1173 inline bool
1175 {
1176  return true;
1177 }
1178 
1179 
1181 
1182 #endif
std::vector< size_type >::const_iterator end_of_row
std::vector< size_type >::const_iterator current_entry
DynamicSparsityPattern get_view(const IndexSet &rows) const
types::global_dof_index size_type
const IndexSet & row_index_set() const
void compute_mmult_pattern(const SparsityPatternTypeLeft &left, const SparsityPatternTypeRight &right)
size_type row_length(const size_type row) const
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false) override
size_type column_index(const size_type row, const size_type col) const
DynamicSparsityPattern & operator=(const DynamicSparsityPattern &)
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type >> &entries)
size_type column_number(const size_type row, const size_type index) const
void print(std::ostream &out) const
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
void clear_row(const size_type row)
void print_gnuplot(std::ostream &out) const
void compute_Tmmult_pattern(const SparsityPatternTypeLeft &left, const SparsityPatternTypeRight &right)
bool exists(const size_type i, const size_type j) const
void add(const size_type i, const size_type j)
ElementIterator at(const size_type global_index) const
Definition: index_set.cc:640
size_type size() const
Definition: index_set.h:1661
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1883
bool is_element(const size_type index) const
Definition: index_set.h:1776
ElementIterator end() const
Definition: index_set.h:1607
size_type n_rows() const
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type >> &entries)
size_type n_cols() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & DummyAccessor()
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:488
static ::ExceptionBase & ExcMessage(std::string arg1)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1016
const types::global_dof_index invalid_size_type
Definition: types.h:222
unsigned int global_dof_index
Definition: types.h:82
void add_entries(ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted)
void add(const size_type col_num)