Reference documentation for deal.II version Git e8f99f7e31 2020-07-08 12:50:43 +0200
\(\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 - 2020 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 
48 {
49  // forward declaration
50  class Iterator;
51 
56 
65  class Accessor
66  {
67  public:
72  const size_type row,
73  const unsigned int index_within_row);
74 
78  Accessor(const DynamicSparsityPattern *sparsity_pattern);
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 
158  // Grant access to iterator class.
159  friend class Iterator;
160  };
161 
162 
163 
188  class Iterator
189  {
190  public:
197  const size_type row,
198  const unsigned int index_within_row);
199 
204  Iterator(const DynamicSparsityPattern *sp);
205 
211  Iterator() = default;
212 
216  Iterator &
217  operator++();
218 
222  Iterator
223  operator++(int);
224 
228  const Accessor &operator*() const;
229 
233  const Accessor *operator->() const;
234 
238  bool
239  operator==(const Iterator &) const;
240 
244  bool
245  operator!=(const Iterator &) const;
246 
254  bool
255  operator<(const Iterator &) const;
256 
263  int
264  operator-(const Iterator &p) const;
265 
266  private:
271  };
272 } // namespace DynamicSparsityPatternIterators
273 
274 
319 {
320 public:
325 
334 
340 
347 
358 
366  const size_type n,
367  const IndexSet &rowset = IndexSet());
368 
374  DynamicSparsityPattern(const IndexSet &indexset);
375 
380 
387  operator=(const DynamicSparsityPattern &);
388 
395  void
396  reinit(const size_type m,
397  const size_type n,
398  const IndexSet &rowset = IndexSet());
399 
405  void
406  compress();
407 
412  bool
413  empty() const;
414 
419  size_type
420  max_entries_per_row() const;
421 
425  void
426  add(const size_type i, const size_type j);
427 
432  template <typename ForwardIterator>
433  void
434  add_entries(const size_type row,
435  ForwardIterator begin,
436  ForwardIterator end,
437  const bool indices_are_unique_and_sorted = false);
438 
442  bool
443  exists(const size_type i, const size_type j) const;
444 
452  get_view(const IndexSet &rows) const;
453 
461  void
462  symmetrize();
463 
468  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
469  void
470  compute_mmult_pattern(const SparsityPatternTypeLeft & left,
471  const SparsityPatternTypeRight &right);
472 
477  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
478  void
479  compute_Tmmult_pattern(const SparsityPatternTypeLeft & left,
480  const SparsityPatternTypeRight &right);
481 
487  void
488  print(std::ostream &out) const;
489 
503  void
504  print_gnuplot(std::ostream &out) const;
505 
509  size_type
510  n_rows() const;
511 
516  size_type
517  n_cols() const;
518 
523  size_type
524  row_length(const size_type row) const;
525 
529  void
530  clear_row(const size_type row);
531 
536  size_type
537  column_number(const size_type row, const size_type index) const;
538 
544  size_type
545  column_index(const size_type row, const size_type col) const;
546 
550  // @{
551 
565  iterator
566  begin() const;
567 
571  iterator
572  end() const;
573 
590  iterator
591  begin(const size_type r) const;
592 
601  iterator
602  end(const size_type r) const;
603 
604  // @}
605 
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
658  stores_only_added_elements();
659 
664  size_type
665  memory_consumption() const;
666 
667 private:
672 
677 
682 
688 
689 
696  struct Line
697  {
698  public:
703  std::vector<size_type> entries;
704 
708  void
709  add(const size_type col_num);
710 
714  template <typename ForwardIterator>
715  void
716  add_entries(ForwardIterator begin,
717  ForwardIterator end,
718  const bool indices_are_sorted);
719 
723  size_type
724  memory_consumption() const;
725  };
726 
727 
731  std::vector<Line> lines;
732 
733  // make the accessor class a friend
735 };
736 
738 /*---------------------- Inline functions -----------------------------------*/
739 
740 
742 {
744  const size_type row,
745  const unsigned int index_within_row)
746  : sparsity_pattern(sparsity_pattern)
747  , current_row(row)
748  , current_entry(
749  ((sparsity_pattern->rowset.size() == 0) ?
750  sparsity_pattern->lines[current_row].entries.begin() :
751  sparsity_pattern
752  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
753  .entries.begin()) +
754  index_within_row)
755  , end_of_row(
756  (sparsity_pattern->rowset.size() == 0) ?
757  sparsity_pattern->lines[current_row].entries.end() :
758  sparsity_pattern
759  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
760  .entries.end())
761  {
762  AssertIndexRange(current_row, sparsity_pattern->n_rows());
763  Assert((sparsity_pattern->rowset.size() == 0) ||
764  sparsity_pattern->rowset.is_element(current_row),
765  ExcMessage("You can't create an iterator into a "
766  "DynamicSparsityPattern's row that is not "
767  "actually stored by that sparsity pattern "
768  "based on the IndexSet argument to it."));
770  index_within_row,
771  ((sparsity_pattern->rowset.size() == 0) ?
772  sparsity_pattern->lines[current_row].entries.size() :
773  sparsity_pattern
774  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
775  .entries.size()));
776  }
777 
778 
780  : sparsity_pattern(sparsity_pattern)
782  , current_entry()
783  , end_of_row()
784  {}
785 
786 
787 
789  : sparsity_pattern(nullptr)
791  , current_entry()
792  , end_of_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_row;
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  }
814 
815 
816  inline size_type
818  {
819  Assert(sparsity_pattern != nullptr, DummyAccessor());
820  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
821 
822  return (current_entry -
823  ((sparsity_pattern->rowset.size() == 0) ?
824  sparsity_pattern->lines[current_row].entries.begin() :
827  .entries.begin()));
828  }
829 
830 
831 
832  inline bool
833  Accessor::operator==(const Accessor &other) const
834  {
835  Assert(sparsity_pattern != nullptr, DummyAccessor());
836  Assert(other.sparsity_pattern != nullptr, DummyAccessor());
837  // compare the sparsity pattern the iterator points into, the
838  // current row, and the location within this row. ignore the
839  // latter if the row is past-the-end because in that case the
840  // current_entry field may not point to a deterministic location
841  return (sparsity_pattern == other.sparsity_pattern &&
842  current_row == other.current_row &&
844  (current_entry == other.current_entry)));
845  }
846 
847 
848 
849  inline bool
850  Accessor::operator<(const Accessor &other) const
851  {
852  Assert(sparsity_pattern != nullptr, DummyAccessor());
853  Assert(other.sparsity_pattern != nullptr, DummyAccessor());
855 
856  // if *this is past-the-end, then it is less than no one
858  return (false);
859  // now *this should be an valid value
860  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
861 
862  // if other is past-the-end
864  return (true);
865  // now other should be an valid value
867 
868  // both iterators are not one-past-the-end
869  return ((current_row < other.current_row) ||
870  ((current_row == other.current_row) &&
871  (current_entry < other.current_entry)));
872  }
873 
874 
875  inline void
877  {
878  Assert(sparsity_pattern != nullptr, DummyAccessor());
879  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
880 
881  // move to the next element in this row
882  ++current_entry;
883 
884  // if this moves us beyond the end of the row, go to the next row
885  // if possible, or set the iterator to an invalid state if not.
886  //
887  // going to the next row is a bit complicated because we may have
888  // to skip over empty rows, and because we also have to avoid rows
889  // that aren't listed in a possibly passed IndexSet argument of
890  // the sparsity pattern. consequently, rather than trying to
891  // duplicate code here, just call the begin() function of the
892  // sparsity pattern itself
893  if (current_entry == end_of_row)
894  {
895  if (current_row + 1 < sparsity_pattern->n_rows())
896  *this = *sparsity_pattern->begin(current_row + 1);
897  else
898  *this = Accessor(sparsity_pattern); // invalid object
899  }
900  }
901 
902 
903 
905  const size_type row,
906  const unsigned int index_within_row)
907  : accessor(sparsity_pattern, row, index_within_row)
908  {}
909 
910 
911 
912  inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern)
913  : accessor(sparsity_pattern)
914  {}
915 
916 
917 
918  inline Iterator &
920  {
921  accessor.advance();
922  return *this;
923  }
924 
925 
926 
927  inline Iterator
929  {
930  const Iterator iter = *this;
931  accessor.advance();
932  return iter;
933  }
934 
935 
936 
937  inline const Accessor &Iterator::operator*() const
938  {
939  return accessor;
940  }
941 
942 
943 
944  inline const Accessor *Iterator::operator->() const
945  {
946  return &accessor;
947  }
948 
949 
950  inline bool
951  Iterator::operator==(const Iterator &other) const
952  {
953  return (accessor == other.accessor);
954  }
955 
956 
957 
958  inline bool
959  Iterator::operator!=(const Iterator &other) const
960  {
961  return !(*this == other);
962  }
963 
964 
965  inline bool
966  Iterator::operator<(const Iterator &other) const
967  {
968  return accessor < other.accessor;
969  }
970 
971 
972  inline int
973  Iterator::operator-(const Iterator &other) const
974  {
975  (void)other;
977  ExcInternalError());
978  Assert(false, ExcNotImplemented());
979 
980  return 0;
981  }
982 } // namespace DynamicSparsityPatternIterators
983 
984 
985 inline void
987 {
988  // first check the last element (or if line is still empty)
989  if ((entries.size() == 0) || (entries.back() < j))
990  {
991  entries.push_back(j);
992  return;
993  }
994 
995  // do a binary search to find the place where to insert:
996  std::vector<size_type>::iterator it =
997  Utilities::lower_bound(entries.begin(), entries.end(), j);
998 
999  // If this entry is a duplicate, exit immediately
1000  if (*it == j)
1001  return;
1002 
1003  // Insert at the right place in the vector. Vector grows automatically to
1004  // fit elements. Always doubles its size.
1005  entries.insert(it, j);
1006 }
1007 
1008 
1009 
1012 {
1013  return rows;
1014 }
1015 
1016 
1017 
1020 {
1021  return cols;
1022 }
1023 
1024 
1025 
1026 inline void
1028 {
1029  AssertIndexRange(i, rows);
1030  AssertIndexRange(j, cols);
1031 
1032  if (rowset.size() > 0 && !rowset.is_element(i))
1033  return;
1034 
1035  have_entries = true;
1036 
1037  const size_type rowindex =
1038  rowset.size() == 0 ? i : rowset.index_within_set(i);
1039  lines[rowindex].add(j);
1040 }
1041 
1042 
1043 
1044 template <typename ForwardIterator>
1045 inline void
1047  ForwardIterator begin,
1048  ForwardIterator end,
1049  const bool indices_are_sorted)
1050 {
1051  AssertIndexRange(row, rows);
1052 
1053  if (rowset.size() > 0 && !rowset.is_element(row))
1054  return;
1055 
1056  if (!have_entries && begin < end)
1057  have_entries = true;
1058 
1059  const size_type rowindex =
1060  rowset.size() == 0 ? row : rowset.index_within_set(row);
1061  lines[rowindex].add_entries(begin, end, indices_are_sorted);
1062 }
1063 
1064 
1065 
1068 {
1069  AssertIndexRange(row, n_rows());
1070 
1071  if (!have_entries)
1072  return 0;
1073 
1074  if (rowset.size() > 0 && !rowset.is_element(row))
1075  return 0;
1076 
1077  const size_type rowindex =
1078  rowset.size() == 0 ? row : rowset.index_within_set(row);
1079  return lines[rowindex].entries.size();
1080 }
1081 
1082 
1083 
1086  const size_type index) const
1087 {
1088  AssertIndexRange(row, n_rows());
1089  Assert(rowset.size() == 0 || rowset.is_element(row), ExcInternalError());
1090 
1091  const size_type local_row =
1092  rowset.size() ? rowset.index_within_set(row) : row;
1093  AssertIndexRange(index, lines[local_row].entries.size());
1094  return lines[local_row].entries[index];
1095 }
1096 
1097 
1098 
1101 {
1102  if (n_rows() > 0)
1103  return begin(0);
1104  else
1105  return end();
1106 }
1107 
1108 
1111 {
1112  return {this};
1113 }
1114 
1115 
1116 
1119 {
1120  AssertIndexRange(r, n_rows());
1121 
1122  if (!have_entries)
1123  return {this};
1124 
1125  if (rowset.size() > 0)
1126  {
1127  // We have an IndexSet that describes the locally owned set. For
1128  // performance reasons we need to make sure that we don't do a
1129  // linear search over 0..n_rows(). Instead, find the first entry
1130  // >= row r in the locally owned set (this is done in log
1131  // n_ranges time inside at()). From there, we move forward until
1132  // we find a non-empty row. By iterating over the IndexSet instead
1133  // of incrementing the row index, we potentially skip over entries
1134  // not in the rowset.
1135  IndexSet::ElementIterator it = rowset.at(r);
1136  if (it == rowset.end())
1137  return end(); // we don't own any row between r and the end
1138 
1139  // Instead of using row_length(*it)==0 in the while loop below,
1140  // which involves an expensive index_within_set() call, we
1141  // look at the lines vector directly. This works, because we are
1142  // walking over this vector entry by entry anyways.
1143  size_type rowindex = rowset.index_within_set(*it);
1144 
1145  while (it != rowset.end() && lines[rowindex].entries.size() == 0)
1146  {
1147  ++it;
1148  ++rowindex;
1149  }
1150 
1151  if (it == rowset.end())
1152  return end();
1153  else
1154  return {this, *it, 0};
1155  }
1156 
1157  // Without an index set we have to do a linear search starting at
1158  // row r until we find a non-empty one. We will check the lines vector
1159  // directly instead of going through the slower row_length() function
1160  size_type row = r;
1161 
1162  while (row < n_rows() && lines[row].entries.size() == 0)
1163  {
1164  ++row;
1165  }
1166 
1167  if (row == n_rows())
1168  return {this};
1169  else
1170  return {this, row, 0};
1171 }
1172 
1173 
1174 
1177 {
1178  AssertIndexRange(r, n_rows());
1179 
1180  unsigned int row = r + 1;
1181  if (row == n_rows())
1182  return {this};
1183  else
1184  return begin(row);
1185 }
1186 
1187 
1188 
1189 inline const IndexSet &
1191 {
1192  return rowset;
1193 }
1194 
1195 
1196 
1197 inline bool
1199 {
1200  return true;
1201 }
1202 
1203 
1205 
1206 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1075
const types::global_dof_index invalid_size_type
Definition: types.h:200
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
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
size_type column_number(const size_type row, const size_type index) const
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
void add(const size_type i, const size_type j)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1628
size_type size() const
Definition: index_set.h:1632
static ::ExceptionBase & ExcMessage(std::string arg1)
std::string compress(const std::string &input)
Definition: utilities.cc:392
static ::ExceptionBase & DummyAccessor()
#define Assert(cond, exc)
Definition: exceptions.h:1403
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1919
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
const IndexSet & row_index_set() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
VectorType::value_type * end(VectorType &V)
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
std::vector< size_type >::const_iterator current_entry
unsigned int global_dof_index
Definition: types.h:76
size_type row_length(const size_type row) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
VectorType::value_type * begin(VectorType &V)
static ::ExceptionBase & ExcNotImplemented()
bool is_element(const size_type index) const
Definition: index_set.h:1763
void add(const size_type col_num)
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()