Reference documentation for deal.II version 9.4.0
\(\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
25
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
66 {
67 public:
72 const size_type row,
73 const unsigned int index_within_row);
74
79
85 Accessor();
86
91 row() const;
92
97 index() const;
98
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
189 {
190 public:
197 const size_type row,
198 const unsigned int index_within_row);
199
205
211 Iterator() = default;
212
216 Iterator &
217 operator++();
218
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{
322public:
327
336
342
349
360
368 const size_type n,
369 const IndexSet &rowset = IndexSet());
370
376 DynamicSparsityPattern(const IndexSet &indexset);
377
382
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
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
512 n_rows() const;
513
519 n_cols() const;
520
526 row_length(const size_type row) const;
527
531 void
532 clear_row(const size_type row);
533
539 column_number(const size_type row, const size_type index) const;
540
547 column_index(const size_type row, const size_type col) const;
548
552 // @{
553
568 begin() const;
569
574 end() const;
575
593 begin(const size_type r) const;
594
604 end(const size_type r) const;
605
606 // @}
607
614 bandwidth() const;
615
621 n_nonzero_elements() const;
622
628 const IndexSet &
629 row_index_set() const;
630
638 nonempty_cols() const;
639
647 nonempty_rows() const;
648
659 static bool
661
667 memory_consumption() const;
668
669private:
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
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{
745 inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern,
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 {
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() :
777 .entries.size()));
778 }
779
780
781 inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern)
782 : sparsity_pattern(sparsity_pattern)
783 , current_row(numbers::invalid_size_type)
784 , current_entry()
785 , end_of_row()
786 {}
787
788
789
791 : sparsity_pattern(nullptr)
792 , current_row(numbers::invalid_size_type)
793 , current_entry()
794 , end_of_row()
795 {}
796
797
798 inline size_type
800 {
802 Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
803
804 return current_row;
805 }
806
807
808 inline size_type
810 {
812 Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
813
814 return *current_entry;
815 }
816
817
818 inline size_type
820 {
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 {
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 {
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 {
881 Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
882
883 // move to the next element in this row
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
896 {
898 *this = *sparsity_pattern->begin(current_row + 1);
899 else
900 *this = Accessor(sparsity_pattern); // invalid object
901 }
902 }
903
904
905
906 inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern,
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 {
924 return *this;
925 }
926
927
928
929 inline Iterator
931 {
932 const Iterator iter = *this;
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;
982 Assert(false, ExcNotImplemented());
983
984 return 0;
985 }
986} // namespace DynamicSparsityPatternIterators
987
988
989inline 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
1030inline void
1032{
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
1048template <typename ForwardIterator>
1049inline 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());
1094
1095 const size_type local_row =
1096 rowset.size() != 0u ? 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{
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.
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{
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
1193inline const IndexSet &
1195{
1196 return rowset;
1197}
1198
1199
1200
1201inline bool
1203{
1204 return true;
1205}
1206
1207
1209
1210#endif
ElementIterator at(const size_type global_index) const
Definition: index_set.h:1532
size_type size() const
Definition: index_set.h:1636
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1923
bool is_element(const size_type index) const
Definition: index_set.h:1767
ElementIterator end() const
Definition: index_set.h:1582
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & DummyAccessor()
static ::ExceptionBase & ExcMessage(std::string arg1)
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
void add_entries(ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted)
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 add(const size_type col_num)
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
void clear_row(const size_type row)
std::vector< size_type >::const_iterator end_of_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)
std::vector< size_type >::const_iterator current_entry
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1153
const types::global_dof_index invalid_size_type
Definition: types.h:210
unsigned int global_dof_index
Definition: types.h:76