Reference documentation for deal.II version Git c93c00f928 2020-04-01 06:59:45 +0200
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
34 // Forward declaration
35 #ifndef DOXYGEN
37 #endif
38 
48 {
49  // forward declaration
50  class Iterator;
51 
56 
68  class Accessor
69  {
70  public:
75  const size_type row,
76  const unsigned int index_within_row);
77 
81  Accessor(const DynamicSparsityPattern *sparsity_pattern);
82 
88  Accessor();
89 
93  size_type
94  row() const;
95 
99  size_type
100  index() const;
101 
105  size_type
106  column() const;
107 
111  bool
112  operator==(const Accessor &) const;
113 
121  bool
122  operator<(const Accessor &) const;
123 
124  protected:
126  "The instance of this class was initialized"
127  " without DynamicSparsityPattern object, which"
128  " means that it is a dummy accessor that can"
129  " not do any operations.");
130 
135 
140 
145  std::vector<size_type>::const_iterator current_entry;
146 
153  std::vector<size_type>::const_iterator end_of_row;
154 
158  void
159  advance();
160 
161  // Grant access to iterator class.
162  friend class Iterator;
163  };
164 
165 
166 
191  class Iterator
192  {
193  public:
200  const size_type row,
201  const unsigned int index_within_row);
202 
207  Iterator(const DynamicSparsityPattern *sp);
208 
214  Iterator() = default;
215 
219  Iterator &
220  operator++();
221 
225  Iterator
226  operator++(int);
227 
231  const Accessor &operator*() const;
232 
236  const Accessor *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 
324 {
325 public:
330 
339 
345 
352 
363 
371  const size_type n,
372  const IndexSet &rowset = IndexSet());
373 
379  DynamicSparsityPattern(const IndexSet &indexset);
380 
385 
392  operator=(const DynamicSparsityPattern &);
393 
400  void
401  reinit(const size_type m,
402  const size_type n,
403  const IndexSet &rowset = IndexSet());
404 
410  void
411  compress();
412 
417  bool
418  empty() const;
419 
424  size_type
425  max_entries_per_row() const;
426 
430  void
431  add(const size_type i, const size_type j);
432 
437  template <typename ForwardIterator>
438  void
439  add_entries(const size_type row,
440  ForwardIterator begin,
441  ForwardIterator end,
442  const bool indices_are_unique_and_sorted = false);
443 
447  bool
448  exists(const size_type i, const size_type j) const;
449 
457  get_view(const IndexSet &rows) const;
458 
466  void
467  symmetrize();
468 
473  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
474  void
475  compute_mmult_pattern(const SparsityPatternTypeLeft & left,
476  const SparsityPatternTypeRight &right);
477 
482  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
483  void
484  compute_Tmmult_pattern(const SparsityPatternTypeLeft & left,
485  const SparsityPatternTypeRight &right);
486 
492  void
493  print(std::ostream &out) const;
494 
508  void
509  print_gnuplot(std::ostream &out) const;
510 
514  size_type
515  n_rows() const;
516 
521  size_type
522  n_cols() const;
523 
528  size_type
529  row_length(const size_type row) const;
530 
534  void
535  clear_row(const size_type row);
536 
541  size_type
542  column_number(const size_type row, const size_type index) const;
543 
549  size_type
550  column_index(const size_type row, const size_type col) const;
551 
555  // @{
556 
570  iterator
571  begin() const;
572 
576  iterator
577  end() const;
578 
595  iterator
596  begin(const size_type r) const;
597 
606  iterator
607  end(const size_type r) const;
608 
609  // @}
610 
616  size_type
617  bandwidth() const;
618 
623  size_type
624  n_nonzero_elements() const;
625 
631  const IndexSet &
632  row_index_set() const;
633 
640  IndexSet
641  nonempty_cols() const;
642 
649  IndexSet
650  nonempty_rows() const;
651 
662  static bool
663  stores_only_added_elements();
664 
669  size_type
670  memory_consumption() const;
671 
672 private:
677 
682 
687 
693 
694 
701  struct Line
702  {
703  public:
708  std::vector<size_type> entries;
709 
713  void
714  add(const size_type col_num);
715 
719  template <typename ForwardIterator>
720  void
721  add_entries(ForwardIterator begin,
722  ForwardIterator end,
723  const bool indices_are_sorted);
724 
728  size_type
729  memory_consumption() const;
730  };
731 
732 
736  std::vector<Line> lines;
737 
738  // make the accessor class a friend
740 };
741 
743 /*---------------------- Inline functions -----------------------------------*/
744 
745 
747 {
749  const size_type row,
750  const unsigned int index_within_row)
751  : sparsity_pattern(sparsity_pattern)
752  , current_row(row)
753  , current_entry(
754  ((sparsity_pattern->rowset.size() == 0) ?
755  sparsity_pattern->lines[current_row].entries.begin() :
756  sparsity_pattern
757  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
758  .entries.begin()) +
759  index_within_row)
760  , end_of_row(
761  (sparsity_pattern->rowset.size() == 0) ?
762  sparsity_pattern->lines[current_row].entries.end() :
763  sparsity_pattern
764  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
765  .entries.end())
766  {
767  AssertIndexRange(current_row, sparsity_pattern->n_rows());
768  Assert((sparsity_pattern->rowset.size() == 0) ||
769  sparsity_pattern->rowset.is_element(current_row),
770  ExcMessage("You can't create an iterator into a "
771  "DynamicSparsityPattern's row that is not "
772  "actually stored by that sparsity pattern "
773  "based on the IndexSet argument to it."));
775  index_within_row,
776  ((sparsity_pattern->rowset.size() == 0) ?
777  sparsity_pattern->lines[current_row].entries.size() :
778  sparsity_pattern
779  ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
780  .entries.size()));
781  }
782 
783 
785  : sparsity_pattern(sparsity_pattern)
786  , current_row(numbers::invalid_size_type)
787  , current_entry()
788  , end_of_row()
789  {}
790 
791 
792 
794  : sparsity_pattern(nullptr)
795  , current_row(numbers::invalid_size_type)
796  , current_entry()
797  , end_of_row()
798  {}
799 
800 
801  inline size_type
803  {
804  Assert(sparsity_pattern != nullptr, DummyAccessor());
805  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
806 
807  return current_row;
808  }
809 
810 
811  inline size_type
813  {
814  Assert(sparsity_pattern != nullptr, DummyAccessor());
815  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
816 
817  return *current_entry;
818  }
819 
820 
821  inline size_type
823  {
824  Assert(sparsity_pattern != nullptr, DummyAccessor());
825  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
826 
827  return (current_entry -
828  ((sparsity_pattern->rowset.size() == 0) ?
829  sparsity_pattern->lines[current_row].entries.begin() :
832  .entries.begin()));
833  }
834 
835 
836 
837  inline bool
838  Accessor::operator==(const Accessor &other) const
839  {
840  Assert(sparsity_pattern != nullptr, DummyAccessor());
841  Assert(other.sparsity_pattern != nullptr, DummyAccessor());
842  // compare the sparsity pattern the iterator points into, the
843  // current row, and the location within this row. ignore the
844  // latter if the row is past-the-end because in that case the
845  // current_entry field may not point to a deterministic location
846  return (sparsity_pattern == other.sparsity_pattern &&
847  current_row == other.current_row &&
849  (current_entry == other.current_entry)));
850  }
851 
852 
853 
854  inline bool
855  Accessor::operator<(const Accessor &other) const
856  {
857  Assert(sparsity_pattern != nullptr, DummyAccessor());
858  Assert(other.sparsity_pattern != nullptr, DummyAccessor());
860 
861  // if *this is past-the-end, then it is less than no one
863  return (false);
864  // now *this should be an valid value
865  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
866 
867  // if other is past-the-end
869  return (true);
870  // now other should be an valid value
872 
873  // both iterators are not one-past-the-end
874  return ((current_row < other.current_row) ||
875  ((current_row == other.current_row) &&
876  (current_entry < other.current_entry)));
877  }
878 
879 
880  inline void
882  {
883  Assert(sparsity_pattern != nullptr, DummyAccessor());
884  Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
885 
886  // move to the next element in this row
887  ++current_entry;
888 
889  // if this moves us beyond the end of the row, go to the next row
890  // if possible, or set the iterator to an invalid state if not.
891  //
892  // going to the next row is a bit complicated because we may have
893  // to skip over empty rows, and because we also have to avoid rows
894  // that aren't listed in a possibly passed IndexSet argument of
895  // the sparsity pattern. consequently, rather than trying to
896  // duplicate code here, just call the begin() function of the
897  // sparsity pattern itself
898  if (current_entry == end_of_row)
899  {
900  if (current_row + 1 < sparsity_pattern->n_rows())
901  *this = *sparsity_pattern->begin(current_row + 1);
902  else
903  *this = Accessor(sparsity_pattern); // invalid object
904  }
905  }
906 
907 
908 
910  const size_type row,
911  const unsigned int index_within_row)
912  : accessor(sparsity_pattern, row, index_within_row)
913  {}
914 
915 
916 
917  inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern)
918  : accessor(sparsity_pattern)
919  {}
920 
921 
922 
923  inline Iterator &
925  {
926  accessor.advance();
927  return *this;
928  }
929 
930 
931 
932  inline Iterator
934  {
935  const Iterator iter = *this;
936  accessor.advance();
937  return iter;
938  }
939 
940 
941 
942  inline const Accessor &Iterator::operator*() const
943  {
944  return accessor;
945  }
946 
947 
948 
949  inline const Accessor *Iterator::operator->() const
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() ? 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.
1140  IndexSet::ElementIterator it = rowset.at(r);
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  unsigned int 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 
1209 DEAL_II_NAMESPACE_CLOSE
1210 
1211 #endif
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1114
const types::global_dof_index invalid_size_type
Definition: types.h:199
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
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
void add(const size_type i, const size_type j)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
size_type size() const
Definition: index_set.h:1635
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & DummyAccessor()
#define Assert(cond, exc)
Definition: exceptions.h:1419
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:1922
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:496
const IndexSet & row_index_set() const
constexpr SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
std::vector< size_type >::const_iterator current_entry
unsigned int global_dof_index
Definition: types.h:77
types::global_dof_index size_type
size_type row_length(const size_type row) const
static ::ExceptionBase & ExcNotImplemented()
bool is_element(const size_type index) const
Definition: index_set.h:1766
void add(const size_type col_num)
static ::ExceptionBase & ExcInternalError()