Reference documentation for deal.II version GIT eef19498cf 2022-10-05 14:45:02+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>
27 
28 #include <algorithm>
29 #include <iostream>
30 #include <vector>
31 
33 
34 // Forward declaration
35 #ifndef DOXYGEN
37 #endif
38 
49 {
50  // forward declaration
51  class Iterator;
52 
57 
66  class Accessor
67  {
68  public:
73  const size_type row,
74  const unsigned int index_within_row);
75 
80 
86  Accessor();
87 
91  size_type
92  row() const;
93 
97  size_type
98  index() const;
99 
103  size_type
104  column() const;
105 
109  bool
110  operator==(const Accessor &) const;
111 
119  bool
120  operator<(const Accessor &) const;
121 
122  protected:
124  "The instance of this class was initialized"
125  " without DynamicSparsityPattern object, which"
126  " means that it is a dummy accessor that can"
127  " not do any operations.");
128 
133 
138 
143  std::vector<size_type>::const_iterator current_entry;
144 
151  std::vector<size_type>::const_iterator end_of_row;
152 
156  void
157  advance();
158 
159  // Grant access to iterator class.
160  friend class Iterator;
161  };
162 
163 
164 
189  class Iterator
190  {
191  public:
198  const size_type row,
199  const unsigned int index_within_row);
200 
205  Iterator(const DynamicSparsityPattern *sp);
206 
212  Iterator() = default;
213 
217  Iterator &
218  operator++();
219 
223  Iterator
224  operator++(int);
225 
229  const Accessor &
230  operator*() const;
231 
235  const Accessor *
236  operator->() const;
237 
241  bool
242  operator==(const Iterator &) const;
243 
247  bool
248  operator!=(const Iterator &) const;
249 
257  bool
258  operator<(const Iterator &) const;
259 
266  int
267  operator-(const Iterator &p) const;
268 
269  private:
274  };
275 } // namespace DynamicSparsityPatternIterators
276 
277 
322 {
323 public:
328 
337 
343 
350 
361 
369  const size_type n,
370  const IndexSet &rowset = IndexSet());
371 
377  DynamicSparsityPattern(const IndexSet &indexset);
378 
383 
391 
398  void
399  reinit(const size_type m,
400  const size_type n,
401  const IndexSet &rowset = IndexSet());
402 
408  void
409  compress();
410 
415  bool
416  empty() const;
417 
422  size_type
423  max_entries_per_row() const;
424 
428  void
429  add(const size_type i, const size_type j);
430 
435  template <typename ForwardIterator>
436  void
437  add_entries(const size_type row,
438  ForwardIterator begin,
439  ForwardIterator end,
440  const bool indices_are_unique_and_sorted = false);
441 
445  bool
446  exists(const size_type i, const size_type j) const;
447 
455  get_view(const IndexSet &rows) const;
456 
464  void
465  symmetrize();
466 
471  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
472  void
473  compute_mmult_pattern(const SparsityPatternTypeLeft & left,
474  const SparsityPatternTypeRight &right);
475 
480  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
481  void
482  compute_Tmmult_pattern(const SparsityPatternTypeLeft & left,
483  const SparsityPatternTypeRight &right);
484 
490  void
491  print(std::ostream &out) const;
492 
506  void
507  print_gnuplot(std::ostream &out) const;
508 
512  size_type
513  n_rows() const;
514 
519  size_type
520  n_cols() const;
521 
526  size_type
527  row_length(const size_type row) const;
528 
532  void
533  clear_row(const size_type row);
534 
539  size_type
540  column_number(const size_type row, const size_type index) const;
541 
547  size_type
548  column_index(const size_type row, const size_type col) const;
549 
553  // @{
554 
568  iterator
569  begin() const;
570 
574  iterator
575  end() const;
576 
593  iterator
594  begin(const size_type r) const;
595 
604  iterator
605  end(const size_type r) const;
606 
614  size_type
615  bandwidth() const;
616 
621  size_type
622  n_nonzero_elements() const;
623 
629  const IndexSet &
630  row_index_set() const;
631 
638  IndexSet
639  nonempty_cols() const;
640 
647  IndexSet
648  nonempty_rows() const;
649 
660  static bool
662 
667  size_type
668  memory_consumption() const;
669 
670 private:
675 
680 
685 
691 
692 
699  struct Line
700  {
701  public:
706  std::vector<size_type> entries;
707 
711  void
712  add(const size_type col_num);
713 
717  template <typename ForwardIterator>
718  void
719  add_entries(ForwardIterator begin,
720  ForwardIterator end,
721  const bool indices_are_sorted);
722 
726  size_type
727  memory_consumption() const;
728  };
729 
730 
734  std::vector<Line> lines;
735 
736  // make the accessor class a friend
738 };
739 
741 /*---------------------- Inline functions -----------------------------------*/
742 
743 
745 {
746  inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern,
747  const size_type row,
748  const unsigned int index_within_row)
749  : sparsity_pattern(sparsity_pattern)
750  , current_row(row)
751  , current_entry(
752  ((sparsity_pattern->rowset.size() == 0) ?
753  sparsity_pattern->lines[current_row].entries.begin() :
754  sparsity_pattern
755  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
756  .entries.begin()) +
757  index_within_row)
758  , end_of_row(
759  (sparsity_pattern->rowset.size() == 0) ?
760  sparsity_pattern->lines[current_row].entries.end() :
761  sparsity_pattern
762  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
763  .entries.end())
764  {
766  Assert((sparsity_pattern->rowset.size() == 0) ||
768  ExcMessage("You can't create an iterator into a "
769  "DynamicSparsityPattern's row that is not "
770  "actually stored by that sparsity pattern "
771  "based on the IndexSet argument to it."));
773  index_within_row,
774  ((sparsity_pattern->rowset.size() == 0) ?
775  sparsity_pattern->lines[current_row].entries.size() :
778  .entries.size()));
779  }
780 
781 
782  inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern)
783  : sparsity_pattern(sparsity_pattern)
784  , current_row(numbers::invalid_size_type)
785  , current_entry()
786  , end_of_row()
787  {}
788 
789 
790 
792  : sparsity_pattern(nullptr)
793  , current_row(numbers::invalid_size_type)
794  , current_entry()
795  , end_of_row()
796  {}
797 
798 
799  inline size_type
801  {
802  Assert(sparsity_pattern != nullptr, DummyAccessor());
803  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
804 
805  return current_row;
806  }
807 
808 
809  inline size_type
811  {
812  Assert(sparsity_pattern != nullptr, DummyAccessor());
813  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
814 
815  return *current_entry;
816  }
817 
818 
819  inline size_type
821  {
822  Assert(sparsity_pattern != nullptr, DummyAccessor());
823  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
824 
825  return (current_entry -
826  ((sparsity_pattern->rowset.size() == 0) ?
827  sparsity_pattern->lines[current_row].entries.begin() :
830  .entries.begin()));
831  }
832 
833 
834 
835  inline bool
836  Accessor::operator==(const Accessor &other) const
837  {
838  Assert(sparsity_pattern != nullptr, DummyAccessor());
839  Assert(other.sparsity_pattern != nullptr, DummyAccessor());
840  // compare the sparsity pattern the iterator points into, the
841  // current row, and the location within this row. ignore the
842  // latter if the row is past-the-end because in that case the
843  // current_entry field may not point to a deterministic location
844  return (sparsity_pattern == other.sparsity_pattern &&
845  current_row == other.current_row &&
847  (current_entry == other.current_entry)));
848  }
849 
850 
851 
852  inline bool
853  Accessor::operator<(const Accessor &other) const
854  {
855  Assert(sparsity_pattern != nullptr, DummyAccessor());
856  Assert(other.sparsity_pattern != nullptr, DummyAccessor());
858 
859  // if *this is past-the-end, then it is less than no one
861  return (false);
862  // now *this should be an valid value
863  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
864 
865  // if other is past-the-end
867  return (true);
868  // now other should be an valid value
870 
871  // both iterators are not one-past-the-end
872  return ((current_row < other.current_row) ||
873  ((current_row == other.current_row) &&
874  (current_entry < other.current_entry)));
875  }
876 
877 
878  inline void
880  {
881  Assert(sparsity_pattern != nullptr, DummyAccessor());
882  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
883 
884  // move to the next element in this row
885  ++current_entry;
886 
887  // if this moves us beyond the end of the row, go to the next row
888  // if possible, or set the iterator to an invalid state if not.
889  //
890  // going to the next row is a bit complicated because we may have
891  // to skip over empty rows, and because we also have to avoid rows
892  // that aren't listed in a possibly passed IndexSet argument of
893  // the sparsity pattern. consequently, rather than trying to
894  // duplicate code here, just call the begin() function of the
895  // sparsity pattern itself
896  if (current_entry == end_of_row)
897  {
898  if (current_row + 1 < sparsity_pattern->n_rows())
899  *this = *sparsity_pattern->begin(current_row + 1);
900  else
901  *this = Accessor(sparsity_pattern); // invalid object
902  }
903  }
904 
905 
906 
907  inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern,
908  const size_type row,
909  const unsigned int index_within_row)
910  : accessor(sparsity_pattern, row, index_within_row)
911  {}
912 
913 
914 
915  inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern)
916  : accessor(sparsity_pattern)
917  {}
918 
919 
920 
921  inline Iterator &
923  {
924  accessor.advance();
925  return *this;
926  }
927 
928 
929 
930  inline Iterator
932  {
933  const Iterator iter = *this;
934  accessor.advance();
935  return iter;
936  }
937 
938 
939 
940  inline const Accessor &
942  {
943  return accessor;
944  }
945 
946 
947 
948  inline const Accessor *
950  {
951  return &accessor;
952  }
953 
954 
955  inline bool
956  Iterator::operator==(const Iterator &other) const
957  {
958  return (accessor == other.accessor);
959  }
960 
961 
962 
963  inline bool
964  Iterator::operator!=(const Iterator &other) const
965  {
966  return !(*this == other);
967  }
968 
969 
970  inline bool
971  Iterator::operator<(const Iterator &other) const
972  {
973  return accessor < other.accessor;
974  }
975 
976 
977  inline int
978  Iterator::operator-(const Iterator &other) const
979  {
980  (void)other;
982  ExcInternalError());
983  Assert(false, ExcNotImplemented());
984 
985  return 0;
986  }
987 } // namespace DynamicSparsityPatternIterators
988 
989 
990 inline void
992 {
993  // first check the last element (or if line is still empty)
994  if ((entries.size() == 0) || (entries.back() < j))
995  {
996  entries.push_back(j);
997  return;
998  }
999 
1000  // do a binary search to find the place where to insert:
1001  std::vector<size_type>::iterator it =
1002  Utilities::lower_bound(entries.begin(), entries.end(), j);
1003 
1004  // If this entry is a duplicate, exit immediately
1005  if (*it == j)
1006  return;
1007 
1008  // Insert at the right place in the vector. Vector grows automatically to
1009  // fit elements. Always doubles its size.
1010  entries.insert(it, j);
1011 }
1012 
1013 
1014 
1017 {
1018  return rows;
1019 }
1020 
1021 
1022 
1025 {
1026  return cols;
1027 }
1028 
1029 
1030 
1031 inline void
1033 {
1034  AssertIndexRange(i, rows);
1035  AssertIndexRange(j, cols);
1036 
1037  if (rowset.size() > 0 && !rowset.is_element(i))
1038  return;
1039 
1040  have_entries = true;
1041 
1042  const size_type rowindex =
1043  rowset.size() == 0 ? i : rowset.index_within_set(i);
1044  lines[rowindex].add(j);
1045 }
1046 
1047 
1048 
1049 template <typename ForwardIterator>
1050 inline void
1052  ForwardIterator begin,
1053  ForwardIterator end,
1054  const bool indices_are_sorted)
1055 {
1056  AssertIndexRange(row, rows);
1057 
1058  if (rowset.size() > 0 && !rowset.is_element(row))
1059  return;
1060 
1061  if (!have_entries && begin < end)
1062  have_entries = true;
1063 
1064  const size_type rowindex =
1065  rowset.size() == 0 ? row : rowset.index_within_set(row);
1066  lines[rowindex].add_entries(begin, end, indices_are_sorted);
1067 }
1068 
1069 
1070 
1073 {
1074  AssertIndexRange(row, n_rows());
1075 
1076  if (!have_entries)
1077  return 0;
1078 
1079  if (rowset.size() > 0 && !rowset.is_element(row))
1080  return 0;
1081 
1082  const size_type rowindex =
1083  rowset.size() == 0 ? row : rowset.index_within_set(row);
1084  return lines[rowindex].entries.size();
1085 }
1086 
1087 
1088 
1091  const size_type index) const
1092 {
1093  AssertIndexRange(row, n_rows());
1094  Assert(rowset.size() == 0 || rowset.is_element(row), ExcInternalError());
1095 
1096  const size_type local_row =
1097  rowset.size() != 0u ? rowset.index_within_set(row) : row;
1098  AssertIndexRange(index, lines[local_row].entries.size());
1099  return lines[local_row].entries[index];
1100 }
1101 
1102 
1103 
1106 {
1107  if (n_rows() > 0)
1108  return begin(0);
1109  else
1110  return end();
1111 }
1112 
1113 
1116 {
1117  return {this};
1118 }
1119 
1120 
1121 
1124 {
1125  AssertIndexRange(r, n_rows());
1126 
1127  if (!have_entries)
1128  return {this};
1129 
1130  if (rowset.size() > 0)
1131  {
1132  // We have an IndexSet that describes the locally owned set. For
1133  // performance reasons we need to make sure that we don't do a
1134  // linear search over 0..n_rows(). Instead, find the first entry
1135  // >= row r in the locally owned set (this is done in log
1136  // n_ranges time inside at()). From there, we move forward until
1137  // we find a non-empty row. By iterating over the IndexSet instead
1138  // of incrementing the row index, we potentially skip over entries
1139  // not in the rowset.
1141  if (it == rowset.end())
1142  return end(); // we don't own any row between r and the end
1143 
1144  // Instead of using row_length(*it)==0 in the while loop below,
1145  // which involves an expensive index_within_set() call, we
1146  // look at the lines vector directly. This works, because we are
1147  // walking over this vector entry by entry anyways.
1148  size_type rowindex = rowset.index_within_set(*it);
1149 
1150  while (it != rowset.end() && lines[rowindex].entries.size() == 0)
1151  {
1152  ++it;
1153  ++rowindex;
1154  }
1155 
1156  if (it == rowset.end())
1157  return end();
1158  else
1159  return {this, *it, 0};
1160  }
1161 
1162  // Without an index set we have to do a linear search starting at
1163  // row r until we find a non-empty one. We will check the lines vector
1164  // directly instead of going through the slower row_length() function
1165  size_type row = r;
1166 
1167  while (row < n_rows() && lines[row].entries.size() == 0)
1168  {
1169  ++row;
1170  }
1171 
1172  if (row == n_rows())
1173  return {this};
1174  else
1175  return {this, row, 0};
1176 }
1177 
1178 
1179 
1182 {
1183  AssertIndexRange(r, n_rows());
1184 
1185  const size_type row = r + 1;
1186  if (row == n_rows())
1187  return {this};
1188  else
1189  return begin(row);
1190 }
1191 
1192 
1193 
1194 inline const IndexSet &
1196 {
1197  return rowset;
1198 }
1199 
1200 
1201 
1202 inline bool
1204 {
1205  return true;
1206 }
1207 
1208 
1210 
1211 #endif
std::vector< size_type >::const_iterator end_of_row
std::vector< size_type >::const_iterator current_entry
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
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
size_type column_index(const size_type row, const size_type col) const
DynamicSparsityPattern & operator=(const DynamicSparsityPattern &)
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:638
size_type size() const
Definition: index_set.h:1626
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1841
bool is_element(const size_type index) const
Definition: index_set.h:1734
ElementIterator end() const
Definition: index_set.h:1572
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:457
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:458
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & DummyAccessor()
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
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:999
const types::global_dof_index invalid_size_type
Definition: types.h:215
unsigned int global_dof_index
Definition: types.h:81
void add_entries(ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted)
void add(const size_type col_num)