Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
dynamic_sparsity_pattern.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2019 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>
23 #include <deal.II/base/subscriptor.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 
32 DEAL_II_NAMESPACE_OPEN
33 
35 
45 {
46  // forward declaration
47  class Iterator;
48 
53 
65  class Accessor
66  {
67  public:
72  const size_type row,
73  const unsigned int index_within_row);
74 
79 
85  Accessor();
86 
90  size_type
91  row() const;
92 
96  size_type
97  index() const;
98 
102  size_type
103  column() const;
104 
108  bool
109  operator==(const Accessor &) const;
110 
118  bool
119  operator<(const Accessor &) const;
120 
121  protected:
123  "The instance of this class was initialized"
124  " without DynamicSparsityPattern object, which"
125  " means that it is a dummy accessor that can"
126  " not do any operations.");
127 
132 
137 
142  std::vector<size_type>::const_iterator current_entry;
143 
150  std::vector<size_type>::const_iterator end_of_row;
151 
155  void
156  advance();
157 
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 &operator*() const;
231 
235  const Accessor *operator->() const;
236 
240  bool
241  operator==(const Iterator &) const;
242 
246  bool
247  operator!=(const Iterator &) const;
248 
256  bool
257  operator<(const Iterator &) const;
258 
265  int
266  operator-(const Iterator &p) const;
267 
268  private:
273  };
274 } // namespace DynamicSparsityPatternIterators
275 
276 
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 
446  bool
447  exists(const size_type i, const size_type j) const;
448 
456  get_view(const IndexSet &rows) const;
457 
465  void
466  symmetrize();
467 
472  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
473  void
474  compute_mmult_pattern(const SparsityPatternTypeLeft & left,
475  const SparsityPatternTypeRight &right);
476 
481  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
482  void
483  compute_Tmmult_pattern(const SparsityPatternTypeLeft & left,
484  const SparsityPatternTypeRight &right);
485 
491  void
492  print(std::ostream &out) const;
493 
507  void
508  print_gnuplot(std::ostream &out) const;
509 
513  size_type
514  n_rows() const;
515 
520  size_type
521  n_cols() const;
522 
527  size_type
528  row_length(const size_type row) const;
529 
533  void
534  clear_row(const size_type row);
535 
540  size_type
541  column_number(const size_type row, const size_type index) const;
542 
548  size_type
549  column_index(const size_type row, const size_type col) const;
550 
554  // @{
555 
569  iterator
570  begin() const;
571 
575  iterator
576  end() const;
577 
594  iterator
595  begin(const size_type r) const;
596 
605  iterator
606  end(const size_type r) const;
607 
608  // @}
609 
615  size_type
616  bandwidth() const;
617 
622  size_type
623  n_nonzero_elements() const;
624 
630  const IndexSet &
631  row_index_set() const;
632 
639  IndexSet
640  nonempty_cols() const;
641 
648  IndexSet
649  nonempty_rows() const;
650 
661  static bool
663 
668  size_type
669  memory_consumption() const;
670 
671 private:
676 
681 
686 
692 
693 
700  struct Line
701  {
702  public:
707  std::vector<size_type> entries;
708 
712  void
713  add(const size_type col_num);
714 
718  template <typename ForwardIterator>
719  void
720  add_entries(ForwardIterator begin,
721  ForwardIterator end,
722  const bool indices_are_sorted);
723 
727  size_type
728  memory_consumption() const;
729  };
730 
731 
735  std::vector<Line> lines;
736 
737  // make the accessor class a friend
739 };
740 
742 /*---------------------- Inline functions -----------------------------------*/
743 
744 
746 {
747  inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern,
748  const size_type row,
749  const unsigned int index_within_row)
750  : sparsity_pattern(sparsity_pattern)
751  , current_row(row)
752  , current_entry(
753  ((sparsity_pattern->rowset.size() == 0) ?
754  sparsity_pattern->lines[current_row].entries.begin() :
755  sparsity_pattern
756  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
757  .entries.begin()) +
758  index_within_row)
759  , end_of_row(
760  (sparsity_pattern->rowset.size() == 0) ?
761  sparsity_pattern->lines[current_row].entries.end() :
762  sparsity_pattern
763  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
764  .entries.end())
765  {
767  Assert((sparsity_pattern->rowset.size() == 0) ||
769  ExcMessage("You can't create an iterator into a "
770  "DynamicSparsityPattern's row that is not "
771  "actually stored by that sparsity pattern "
772  "based on the IndexSet argument to it."));
774  index_within_row,
775  ((sparsity_pattern->rowset.size() == 0) ?
776  sparsity_pattern->lines[current_row].entries.size() :
779  .entries.size()));
780  }
781 
782 
783  inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern)
784  : sparsity_pattern(sparsity_pattern)
785  , current_row(numbers::invalid_size_type)
786  , current_entry()
787  , end_of_row()
788  {}
789 
790 
791 
793  : sparsity_pattern(nullptr)
794  , current_row(numbers::invalid_size_type)
795  , current_entry()
796  , end_of_row()
797  {}
798 
799 
800  inline size_type
802  {
803  Assert(sparsity_pattern != nullptr, DummyAccessor());
804  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
805 
806  return current_row;
807  }
808 
809 
810  inline size_type
812  {
813  Assert(sparsity_pattern != nullptr, DummyAccessor());
814  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
815 
816  return *current_entry;
817  }
818 
819 
820  inline size_type
822  {
823  Assert(sparsity_pattern != nullptr, DummyAccessor());
824  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
825 
826  return (current_entry -
827  ((sparsity_pattern->rowset.size() == 0) ?
828  sparsity_pattern->lines[current_row].entries.begin() :
831  .entries.begin()));
832  }
833 
834 
835 
836  inline bool
837  Accessor::operator==(const Accessor &other) const
838  {
839  Assert(sparsity_pattern != nullptr, DummyAccessor());
840  Assert(other.sparsity_pattern != nullptr, DummyAccessor());
841  // compare the sparsity pattern the iterator points into, the
842  // current row, and the location within this row. ignore the
843  // latter if the row is past-the-end because in that case the
844  // current_entry field may not point to a deterministic location
845  return (sparsity_pattern == other.sparsity_pattern &&
846  current_row == other.current_row &&
848  (current_entry == other.current_entry)));
849  }
850 
851 
852 
853  inline bool
854  Accessor::operator<(const Accessor &other) const
855  {
856  Assert(sparsity_pattern != nullptr, DummyAccessor());
857  Assert(other.sparsity_pattern != nullptr, DummyAccessor());
859 
860  // if *this is past-the-end, then it is less than no one
862  return (false);
863  // now *this should be an valid value
864  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
865 
866  // if other is past-the-end
868  return (true);
869  // now other should be an valid value
871 
872  // both iterators are not one-past-the-end
873  return ((current_row < other.current_row) ||
874  ((current_row == other.current_row) &&
875  (current_entry < other.current_entry)));
876  }
877 
878 
879  inline void
881  {
882  Assert(sparsity_pattern != nullptr, DummyAccessor());
883  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
884 
885  // move to the next element in this row
886  ++current_entry;
887 
888  // if this moves us beyond the end of the row, go to the next row
889  // if possible, or set the iterator to an invalid state if not.
890  //
891  // going to the next row is a bit complicated because we may have
892  // to skip over empty rows, and because we also have to avoid rows
893  // that aren't listed in a possibly passed IndexSet argument of
894  // the sparsity pattern. consequently, rather than trying to
895  // duplicate code here, just call the begin() function of the
896  // sparsity pattern itself
897  if (current_entry == end_of_row)
898  {
899  if (current_row + 1 < sparsity_pattern->n_rows())
900  *this = *sparsity_pattern->begin(current_row + 1);
901  else
902  *this = Accessor(sparsity_pattern); // invalid object
903  }
904  }
905 
906 
907 
908  inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern,
909  const size_type row,
910  const unsigned int index_within_row)
911  : accessor(sparsity_pattern, row, index_within_row)
912  {}
913 
914 
915 
916  inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern)
917  : accessor(sparsity_pattern)
918  {}
919 
920 
921 
922  inline Iterator &
924  {
925  accessor.advance();
926  return *this;
927  }
928 
929 
930 
931  inline Iterator
933  {
934  const Iterator iter = *this;
935  accessor.advance();
936  return iter;
937  }
938 
939 
940 
941  inline const Accessor &Iterator::operator*() const
942  {
943  return accessor;
944  }
945 
946 
947 
948  inline const Accessor *Iterator::operator->() const
949  {
950  return &accessor;
951  }
952 
953 
954  inline bool
955  Iterator::operator==(const Iterator &other) const
956  {
957  return (accessor == other.accessor);
958  }
959 
960 
961 
962  inline bool
963  Iterator::operator!=(const Iterator &other) const
964  {
965  return !(*this == other);
966  }
967 
968 
969  inline bool
970  Iterator::operator<(const Iterator &other) const
971  {
972  return accessor < other.accessor;
973  }
974 
975 
976  inline int
977  Iterator::operator-(const Iterator &other) const
978  {
979  (void)other;
981  ExcInternalError());
982  Assert(false, ExcNotImplemented());
983 
984  return 0;
985  }
986 } // namespace DynamicSparsityPatternIterators
987 
988 
989 inline void
991 {
992  // first check the last element (or if line is still empty)
993  if ((entries.size() == 0) || (entries.back() < j))
994  {
995  entries.push_back(j);
996  return;
997  }
998 
999  // do a binary search to find the place where to insert:
1000  std::vector<size_type>::iterator it =
1001  Utilities::lower_bound(entries.begin(), entries.end(), j);
1002 
1003  // If this entry is a duplicate, exit immediately
1004  if (*it == j)
1005  return;
1006 
1007  // Insert at the right place in the vector. Vector grows automatically to
1008  // fit elements. Always doubles its size.
1009  entries.insert(it, j);
1010 }
1011 
1012 
1013 
1016 {
1017  return rows;
1018 }
1019 
1020 
1021 
1024 {
1025  return cols;
1026 }
1027 
1028 
1029 
1030 inline void
1032 {
1033  Assert(i < rows, ExcIndexRangeType<size_type>(i, 0, rows));
1034  Assert(j < cols, ExcIndexRangeType<size_type>(j, 0, cols));
1035 
1036  if (rowset.size() > 0 && !rowset.is_element(i))
1037  return;
1038 
1039  have_entries = true;
1040 
1041  const size_type rowindex =
1042  rowset.size() == 0 ? i : rowset.index_within_set(i);
1043  lines[rowindex].add(j);
1044 }
1045 
1046 
1047 
1048 template <typename ForwardIterator>
1049 inline void
1051  ForwardIterator begin,
1052  ForwardIterator end,
1053  const bool indices_are_sorted)
1054 {
1055  Assert(row < rows, ExcIndexRangeType<size_type>(row, 0, rows));
1056 
1057  if (rowset.size() > 0 && !rowset.is_element(row))
1058  return;
1059 
1060  if (!have_entries && begin < end)
1061  have_entries = true;
1062 
1063  const size_type rowindex =
1064  rowset.size() == 0 ? row : rowset.index_within_set(row);
1065  lines[rowindex].add_entries(begin, end, indices_are_sorted);
1066 }
1067 
1068 
1069 
1072 {
1073  Assert(row < n_rows(), ExcIndexRangeType<size_type>(row, 0, n_rows()));
1074 
1075  if (!have_entries)
1076  return 0;
1077 
1078  if (rowset.size() > 0 && !rowset.is_element(row))
1079  return 0;
1080 
1081  const size_type rowindex =
1082  rowset.size() == 0 ? row : rowset.index_within_set(row);
1083  return lines[rowindex].entries.size();
1084 }
1085 
1086 
1087 
1090  const size_type index) const
1091 {
1092  Assert(row < n_rows(), ExcIndexRangeType<size_type>(row, 0, n_rows()));
1093  Assert(rowset.size() == 0 || rowset.is_element(row), ExcInternalError());
1094 
1095  const size_type local_row =
1096  rowset.size() ? rowset.index_within_set(row) : row;
1097  Assert(index < lines[local_row].entries.size(),
1098  ExcIndexRangeType<size_type>(index,
1099  0,
1100  lines[local_row].entries.size()));
1101  return lines[local_row].entries[index];
1102 }
1103 
1104 
1105 
1108 {
1109  if (n_rows() > 0)
1110  return begin(0);
1111  else
1112  return end();
1113 }
1114 
1115 
1118 {
1119  return {this};
1120 }
1121 
1122 
1123 
1126 {
1127  Assert(r < n_rows(), ExcIndexRangeType<size_type>(r, 0, n_rows()));
1128 
1129  if (!have_entries)
1130  return {this};
1131 
1132  if (rowset.size() > 0)
1133  {
1134  // We have an IndexSet that describes the locally owned set. For
1135  // performance reasons we need to make sure that we don't do a
1136  // linear search over 0..n_rows(). Instead, find the first entry
1137  // >= row r in the locally owned set (this is done in log
1138  // n_ranges time inside at()). From there, we move forward until
1139  // we find a non-empty row. By iterating over the IndexSet instead
1140  // of incrementing the row index, we potentially skip over entries
1141  // not in the rowset.
1143  if (it == rowset.end())
1144  return end(); // we don't own any row between r and the end
1145 
1146  // Instead of using row_length(*it)==0 in the while loop below,
1147  // which involves an expensive index_within_set() call, we
1148  // look at the lines vector directly. This works, because we are
1149  // walking over this vector entry by entry anyways.
1150  size_type rowindex = rowset.index_within_set(*it);
1151 
1152  while (it != rowset.end() && lines[rowindex].entries.size() == 0)
1153  {
1154  ++it;
1155  ++rowindex;
1156  }
1157 
1158  if (it == rowset.end())
1159  return end();
1160  else
1161  return {this, *it, 0};
1162  }
1163 
1164  // Without an index set we have to do a linear search starting at
1165  // row r until we find a non-empty one. We will check the lines vector
1166  // directly instead of going through the slower row_length() function
1167  size_type row = r;
1168 
1169  while (row < n_rows() && lines[row].entries.size() == 0)
1170  {
1171  ++row;
1172  }
1173 
1174  if (row == n_rows())
1175  return {this};
1176  else
1177  return {this, row, 0};
1178 }
1179 
1180 
1181 
1184 {
1185  Assert(r < n_rows(), ExcIndexRangeType<size_type>(r, 0, n_rows()));
1186 
1187  unsigned int row = r + 1;
1188  if (row == n_rows())
1189  return {this};
1190  else
1191  return begin(row);
1192 }
1193 
1194 
1195 
1196 inline const IndexSet &
1198 {
1199  return rowset;
1200 }
1201 
1202 
1203 
1204 inline bool
1206 {
1207  return true;
1208 }
1209 
1210 
1211 DEAL_II_NAMESPACE_CLOSE
1212 
1213 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1067
const types::global_dof_index invalid_size_type
Definition: types.h:182
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
std::vector< size_type >::const_iterator end_of_row
ElementIterator end() const
Definition: index_set.h:1546
void clear_row(const size_type row)
size_type column_number(const size_type row, const size_type index) const
void add(const size_type i, const size_type j)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
size_type size() const
Definition: index_set.h:1600
static ::ExceptionBase & ExcMessage(std::string arg1)
void compute_Tmmult_pattern(const SparsityPatternTypeLeft &left, const SparsityPatternTypeRight &right)
static ::ExceptionBase & DummyAccessor()
DynamicSparsityPattern get_view(const IndexSet &rows) const
#define Assert(cond, exc)
Definition: exceptions.h:1407
void add_entries(ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted)
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1821
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:496
const IndexSet & row_index_set() const
void print_gnuplot(std::ostream &out) const
void compute_mmult_pattern(const SparsityPatternTypeLeft &left, const SparsityPatternTypeRight &right)
std::vector< size_type >::const_iterator current_entry
unsigned int global_dof_index
Definition: types.h:89
types::global_dof_index size_type
size_type row_length(const size_type row) const
DynamicSparsityPattern & operator=(const DynamicSparsityPattern &)
void print(std::ostream &out) const
static ::ExceptionBase & ExcNotImplemented()
bool is_element(const size_type index) const
Definition: index_set.h:1665
void add(const size_type col_num)
ElementIterator at(const size_type global_index) const
Definition: index_set.h:1495
bool exists(const size_type i, const size_type j) const
size_type column_index(const size_type row, const size_type col) const
static ::ExceptionBase & ExcInternalError()