Reference documentation for deal.II version Git 8ba89c874e 2021-08-04 22:06:33 -0600
\(\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 - 2021 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 &
229  operator*() const;
230 
234  const Accessor *
235  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 
321 {
322 public:
327 
336 
342 
349 
360 
368  const size_type n,
369  const IndexSet &rowset = IndexSet());
370 
376  DynamicSparsityPattern(const IndexSet &indexset);
377 
382 
389  operator=(const DynamicSparsityPattern &);
390 
397  void
398  reinit(const size_type m,
399  const size_type n,
400  const IndexSet &rowset = IndexSet());
401 
407  void
408  compress();
409 
414  bool
415  empty() const;
416 
421  size_type
422  max_entries_per_row() const;
423 
427  void
428  add(const size_type i, const size_type j);
429 
434  template <typename ForwardIterator>
435  void
436  add_entries(const size_type row,
437  ForwardIterator begin,
438  ForwardIterator end,
439  const bool indices_are_unique_and_sorted = false);
440 
444  bool
445  exists(const size_type i, const size_type j) const;
446 
454  get_view(const IndexSet &rows) const;
455 
463  void
464  symmetrize();
465 
470  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
471  void
472  compute_mmult_pattern(const SparsityPatternTypeLeft & left,
473  const SparsityPatternTypeRight &right);
474 
479  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
480  void
481  compute_Tmmult_pattern(const SparsityPatternTypeLeft & left,
482  const SparsityPatternTypeRight &right);
483 
489  void
490  print(std::ostream &out) const;
491 
505  void
506  print_gnuplot(std::ostream &out) const;
507 
511  size_type
512  n_rows() const;
513 
518  size_type
519  n_cols() const;
520 
525  size_type
526  row_length(const size_type row) const;
527 
531  void
532  clear_row(const size_type row);
533 
538  size_type
539  column_number(const size_type row, const size_type index) const;
540 
546  size_type
547  column_index(const size_type row, const size_type col) const;
548 
552  // @{
553 
567  iterator
568  begin() const;
569 
573  iterator
574  end() const;
575 
592  iterator
593  begin(const size_type r) const;
594 
603  iterator
604  end(const size_type r) const;
605 
606  // @}
607 
613  size_type
614  bandwidth() const;
615 
620  size_type
621  n_nonzero_elements() const;
622 
628  const IndexSet &
629  row_index_set() const;
630 
637  IndexSet
638  nonempty_cols() const;
639 
646  IndexSet
647  nonempty_rows() const;
648 
659  static bool
660  stores_only_added_elements();
661 
666  size_type
667  memory_consumption() const;
668 
669 private:
674 
679 
684 
690 
691 
698  struct Line
699  {
700  public:
705  std::vector<size_type> entries;
706 
710  void
711  add(const size_type col_num);
712 
716  template <typename ForwardIterator>
717  void
718  add_entries(ForwardIterator begin,
719  ForwardIterator end,
720  const bool indices_are_sorted);
721 
725  size_type
726  memory_consumption() const;
727  };
728 
729 
733  std::vector<Line> lines;
734 
735  // make the accessor class a friend
737 };
738 
740 /*---------------------- Inline functions -----------------------------------*/
741 
742 
744 {
746  const size_type row,
747  const unsigned int index_within_row)
748  : sparsity_pattern(sparsity_pattern)
749  , current_row(row)
750  , current_entry(
751  ((sparsity_pattern->rowset.size() == 0) ?
752  sparsity_pattern->lines[current_row].entries.begin() :
753  sparsity_pattern
754  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
755  .entries.begin()) +
756  index_within_row)
757  , end_of_row(
758  (sparsity_pattern->rowset.size() == 0) ?
759  sparsity_pattern->lines[current_row].entries.end() :
760  sparsity_pattern
761  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
762  .entries.end())
763  {
764  AssertIndexRange(current_row, sparsity_pattern->n_rows());
765  Assert((sparsity_pattern->rowset.size() == 0) ||
766  sparsity_pattern->rowset.is_element(current_row),
767  ExcMessage("You can't create an iterator into a "
768  "DynamicSparsityPattern's row that is not "
769  "actually stored by that sparsity pattern "
770  "based on the IndexSet argument to it."));
772  index_within_row,
773  ((sparsity_pattern->rowset.size() == 0) ?
774  sparsity_pattern->lines[current_row].entries.size() :
775  sparsity_pattern
776  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
777  .entries.size()));
778  }
779 
780 
782  : sparsity_pattern(sparsity_pattern)
784  , current_entry()
785  , end_of_row()
786  {}
787 
788 
789 
791  : sparsity_pattern(nullptr)
793  , current_entry()
794  , end_of_row()
795  {}
796 
797 
798  inline size_type
800  {
801  Assert(sparsity_pattern != nullptr, DummyAccessor());
802  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
803 
804  return current_row;
805  }
806 
807 
808  inline size_type
810  {
811  Assert(sparsity_pattern != nullptr, DummyAccessor());
812  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
813 
814  return *current_entry;
815  }
816 
817 
818  inline size_type
820  {
821  Assert(sparsity_pattern != nullptr, DummyAccessor());
822  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
823 
824  return (current_entry -
825  ((sparsity_pattern->rowset.size() == 0) ?
826  sparsity_pattern->lines[current_row].entries.begin() :
829  .entries.begin()));
830  }
831 
832 
833 
834  inline bool
835  Accessor::operator==(const Accessor &other) const
836  {
837  Assert(sparsity_pattern != nullptr, DummyAccessor());
838  Assert(other.sparsity_pattern != nullptr, DummyAccessor());
839  // compare the sparsity pattern the iterator points into, the
840  // current row, and the location within this row. ignore the
841  // latter if the row is past-the-end because in that case the
842  // current_entry field may not point to a deterministic location
843  return (sparsity_pattern == other.sparsity_pattern &&
844  current_row == other.current_row &&
846  (current_entry == other.current_entry)));
847  }
848 
849 
850 
851  inline bool
852  Accessor::operator<(const Accessor &other) const
853  {
854  Assert(sparsity_pattern != nullptr, DummyAccessor());
855  Assert(other.sparsity_pattern != nullptr, DummyAccessor());
857 
858  // if *this is past-the-end, then it is less than no one
860  return (false);
861  // now *this should be an valid value
862  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
863 
864  // if other is past-the-end
866  return (true);
867  // now other should be an valid value
869 
870  // both iterators are not one-past-the-end
871  return ((current_row < other.current_row) ||
872  ((current_row == other.current_row) &&
873  (current_entry < other.current_entry)));
874  }
875 
876 
877  inline void
879  {
880  Assert(sparsity_pattern != nullptr, DummyAccessor());
881  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
882 
883  // move to the next element in this row
884  ++current_entry;
885 
886  // if this moves us beyond the end of the row, go to the next row
887  // if possible, or set the iterator to an invalid state if not.
888  //
889  // going to the next row is a bit complicated because we may have
890  // to skip over empty rows, and because we also have to avoid rows
891  // that aren't listed in a possibly passed IndexSet argument of
892  // the sparsity pattern. consequently, rather than trying to
893  // duplicate code here, just call the begin() function of the
894  // sparsity pattern itself
895  if (current_entry == end_of_row)
896  {
897  if (current_row + 1 < sparsity_pattern->n_rows())
898  *this = *sparsity_pattern->begin(current_row + 1);
899  else
900  *this = Accessor(sparsity_pattern); // invalid object
901  }
902  }
903 
904 
905 
907  const size_type row,
908  const unsigned int index_within_row)
909  : accessor(sparsity_pattern, row, index_within_row)
910  {}
911 
912 
913 
914  inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern)
915  : accessor(sparsity_pattern)
916  {}
917 
918 
919 
920  inline Iterator &
922  {
923  accessor.advance();
924  return *this;
925  }
926 
927 
928 
929  inline Iterator
931  {
932  const Iterator iter = *this;
933  accessor.advance();
934  return iter;
935  }
936 
937 
938 
939  inline const Accessor &
941  {
942  return accessor;
943  }
944 
945 
946 
947  inline const Accessor *
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  AssertIndexRange(i, rows);
1034  AssertIndexRange(j, 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  AssertIndexRange(row, 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  AssertIndexRange(row, 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  AssertIndexRange(row, 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  AssertIndexRange(index, lines[local_row].entries.size());
1098  return lines[local_row].entries[index];
1099 }
1100 
1101 
1102 
1105 {
1106  if (n_rows() > 0)
1107  return begin(0);
1108  else
1109  return end();
1110 }
1111 
1112 
1115 {
1116  return {this};
1117 }
1118 
1119 
1120 
1123 {
1124  AssertIndexRange(r, n_rows());
1125 
1126  if (!have_entries)
1127  return {this};
1128 
1129  if (rowset.size() > 0)
1130  {
1131  // We have an IndexSet that describes the locally owned set. For
1132  // performance reasons we need to make sure that we don't do a
1133  // linear search over 0..n_rows(). Instead, find the first entry
1134  // >= row r in the locally owned set (this is done in log
1135  // n_ranges time inside at()). From there, we move forward until
1136  // we find a non-empty row. By iterating over the IndexSet instead
1137  // of incrementing the row index, we potentially skip over entries
1138  // not in the rowset.
1139  IndexSet::ElementIterator it = rowset.at(r);
1140  if (it == rowset.end())
1141  return end(); // we don't own any row between r and the end
1142 
1143  // Instead of using row_length(*it)==0 in the while loop below,
1144  // which involves an expensive index_within_set() call, we
1145  // look at the lines vector directly. This works, because we are
1146  // walking over this vector entry by entry anyways.
1147  size_type rowindex = rowset.index_within_set(*it);
1148 
1149  while (it != rowset.end() && lines[rowindex].entries.size() == 0)
1150  {
1151  ++it;
1152  ++rowindex;
1153  }
1154 
1155  if (it == rowset.end())
1156  return end();
1157  else
1158  return {this, *it, 0};
1159  }
1160 
1161  // Without an index set we have to do a linear search starting at
1162  // row r until we find a non-empty one. We will check the lines vector
1163  // directly instead of going through the slower row_length() function
1164  size_type row = r;
1165 
1166  while (row < n_rows() && lines[row].entries.size() == 0)
1167  {
1168  ++row;
1169  }
1170 
1171  if (row == n_rows())
1172  return {this};
1173  else
1174  return {this, row, 0};
1175 }
1176 
1177 
1178 
1181 {
1182  AssertIndexRange(r, n_rows());
1183 
1184  const size_type row = r + 1;
1185  if (row == n_rows())
1186  return {this};
1187  else
1188  return begin(row);
1189 }
1190 
1191 
1192 
1193 inline const IndexSet &
1195 {
1196  return rowset;
1197 }
1198 
1199 
1200 
1201 inline bool
1203 {
1204  return true;
1205 }
1206 
1207 
1209 
1210 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1111
const types::global_dof_index invalid_size_type
Definition: types.h:205
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
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
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:1724
size_type size() const
Definition: index_set.h:1639
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:1467
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1926
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
const IndexSet & row_index_set() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
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:400
VectorType::value_type * begin(VectorType &V)
static ::ExceptionBase & ExcNotImplemented()
bool is_element(const size_type index) const
Definition: index_set.h:1770
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()