Reference documentation for deal.II version GIT relicensing-446-ge11b936273 2024-04-20 12:30:02+00:00
\(\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\}}\)
Loading...
Searching...
No Matches
dynamic_sparsity_pattern.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2011 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_dynamic_sparsity_pattern_h
16#define dealii_dynamic_sparsity_pattern_h
17
18
19#include <deal.II/base/config.h>
20
24
27
28#include <algorithm>
29#include <iostream>
30#include <vector>
31
33
34// Forward declaration
35#ifndef DOXYGEN
37#endif
38
49{
50 // forward declaration
51 class Iterator;
52
57
67 {
68 public:
73 const size_type row,
74 const unsigned int index_within_row);
75
80
86 Accessor();
87
92 row() const;
93
98 index() const;
99
104 column() const;
105
109 bool
110 operator==(const Accessor &) const;
111
119 bool
120 operator<(const Accessor &) const;
121
122 protected:
124 "The instance of this class was initialized"
125 " without DynamicSparsityPattern object, which"
126 " means that it is a dummy accessor that can"
127 " not do any operations.");
128
133
138
143 std::vector<size_type>::const_iterator current_entry;
144
151 std::vector<size_type>::const_iterator end_of_row;
152
156 void
157 advance();
158
159 // Grant access to iterator class.
160 friend class Iterator;
161 };
162
163
164
190 {
191 public:
198 const size_type row,
199 const unsigned int index_within_row);
200
206
212 Iterator() = default;
213
217 Iterator &
218 operator++();
219
224 operator++(int);
225
229 const Accessor &
230 operator*() const;
231
235 const Accessor *
236 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
322{
323public:
328
337
343
350
361
369 const size_type n,
370 const IndexSet &rowset = IndexSet());
371
377 DynamicSparsityPattern(const IndexSet &indexset);
378
383
391
398 void
399 reinit(const size_type m,
400 const size_type n,
401 const IndexSet &rowset = IndexSet());
402
408 void
409 compress();
410
415 bool
416 empty() const;
417
423 max_entries_per_row() const;
424
428 void
429 add(const size_type i, const size_type j);
430
435 template <typename ForwardIterator>
436 void
437 add_entries(const size_type row,
438 ForwardIterator begin,
439 ForwardIterator end,
440 const bool indices_are_unique_and_sorted = false);
441
442 virtual void
443 add_row_entries(const size_type &row,
444 const ArrayView<const size_type> &columns,
445 const bool indices_are_sorted = false) override;
446
448
452 bool
453 exists(const size_type i, const size_type j) const;
454
462 get_view(const IndexSet &rows) const;
463
471 void
472 symmetrize();
473
478 template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
479 void
480 compute_mmult_pattern(const SparsityPatternTypeLeft &left,
481 const SparsityPatternTypeRight &right);
482
487 template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
488 void
489 compute_Tmmult_pattern(const SparsityPatternTypeLeft &left,
490 const SparsityPatternTypeRight &right);
491
497 void
498 print(std::ostream &out) const;
499
513 void
514 print_gnuplot(std::ostream &out) const;
515
521 row_length(const size_type row) const;
522
526 void
527 clear_row(const size_type row);
528
534 column_number(const size_type row, const size_type index) const;
535
542 column_index(const size_type row, const size_type col) const;
543
563 begin() const;
564
569 end() const;
570
588 begin(const size_type r) const;
589
599 end(const size_type r) const;
600
611 bandwidth() const;
612
618 n_nonzero_elements() const;
619
625 const IndexSet &
626 row_index_set() const;
627
635 nonempty_cols() const;
636
644 nonempty_rows() const;
645
656 static bool
658
664 memory_consumption() const;
665
666private:
671
677
678
685 struct Line
686 {
687 public:
692 std::vector<size_type> entries;
693
697 void
698 add(const size_type col_num);
699
703 template <typename ForwardIterator>
704 void
705 add_entries(ForwardIterator begin,
706 ForwardIterator end,
707 const bool indices_are_sorted);
708
713 memory_consumption() const;
714 };
715
716
720 std::vector<Line> lines;
721
722 // make the accessor class a friend
724};
725
727/*---------------------- Inline functions -----------------------------------*/
728
729
731{
732 inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern,
733 const size_type row,
734 const unsigned int index_within_row)
735 : sparsity_pattern(sparsity_pattern)
736 , current_row(row)
737 , current_entry(
738 ((sparsity_pattern->rowset.size() == 0) ?
739 sparsity_pattern->lines[current_row].entries.begin() :
740 sparsity_pattern
741 ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
742 .entries.begin()) +
743 index_within_row)
744 , end_of_row(
745 (sparsity_pattern->rowset.size() == 0) ?
746 sparsity_pattern->lines[current_row].entries.end() :
747 sparsity_pattern
748 ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
749 .entries.end())
750 {
754 ExcMessage("You can't create an iterator into a "
755 "DynamicSparsityPattern's row that is not "
756 "actually stored by that sparsity pattern "
757 "based on the IndexSet argument to it."));
759 index_within_row,
760 ((sparsity_pattern->rowset.size() == 0) ?
761 sparsity_pattern->lines[current_row].entries.size() :
764 .entries.size()));
765 }
766
767
768 inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern)
769 : sparsity_pattern(sparsity_pattern)
770 , current_row(numbers::invalid_size_type)
771 , current_entry()
772 , end_of_row()
773 {}
774
775
776
778 : sparsity_pattern(nullptr)
779 , current_row(numbers::invalid_size_type)
780 , current_entry()
781 , end_of_row()
782 {}
783
784
785 inline size_type
787 {
789 Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
790
791 return current_row;
792 }
793
794
795 inline size_type
797 {
799 Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
800
801 return *current_entry;
802 }
803
804
805 inline size_type
807 {
809 Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
810
811 return (current_entry -
812 ((sparsity_pattern->rowset.size() == 0) ?
813 sparsity_pattern->lines[current_row].entries.begin() :
816 .entries.begin()));
817 }
818
819
820
821 inline bool
822 Accessor::operator==(const Accessor &other) const
823 {
825 Assert(other.sparsity_pattern != nullptr, DummyAccessor());
826 // compare the sparsity pattern the iterator points into, the
827 // current row, and the location within this row. ignore the
828 // latter if the row is past-the-end because in that case the
829 // current_entry field may not point to a deterministic location
830 return (sparsity_pattern == other.sparsity_pattern &&
831 current_row == other.current_row &&
833 (current_entry == other.current_entry)));
834 }
835
836
837
838 inline bool
839 Accessor::operator<(const Accessor &other) const
840 {
842 Assert(other.sparsity_pattern != nullptr, DummyAccessor());
844
845 // if *this is past-the-end, then it is less than no one
847 return (false);
848 // now *this should be an valid value
849 Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
850
851 // if other is past-the-end
853 return (true);
854 // now other should be an valid value
856
857 // both iterators are not one-past-the-end
858 return ((current_row < other.current_row) ||
859 ((current_row == other.current_row) &&
860 (current_entry < other.current_entry)));
861 }
862
863
864 inline void
866 {
868 Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
869
870 // move to the next element in this row
872
873 // if this moves us beyond the end of the row, go to the next row
874 // if possible, or set the iterator to an invalid state if not.
875 //
876 // going to the next row is a bit complicated because we may have
877 // to skip over empty rows, and because we also have to avoid rows
878 // that aren't listed in a possibly passed IndexSet argument of
879 // the sparsity pattern. consequently, rather than trying to
880 // duplicate code here, just call the begin() function of the
881 // sparsity pattern itself
883 {
885 *this = *sparsity_pattern->begin(current_row + 1);
886 else
887 *this = Accessor(sparsity_pattern); // invalid object
888 }
889 }
890
891
892
893 inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern,
894 const size_type row,
895 const unsigned int index_within_row)
896 : accessor(sparsity_pattern, row, index_within_row)
897 {}
898
899
900
901 inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern)
902 : accessor(sparsity_pattern)
903 {}
904
905
906
907 inline Iterator &
909 {
911 return *this;
912 }
913
914
915
916 inline Iterator
918 {
919 const Iterator iter = *this;
921 return iter;
922 }
923
924
925
926 inline const Accessor &
928 {
929 return accessor;
930 }
931
932
933
934 inline const Accessor *
936 {
937 return &accessor;
938 }
939
940
941 inline bool
942 Iterator::operator==(const Iterator &other) const
943 {
944 return (accessor == other.accessor);
945 }
946
947
948
949 inline bool
950 Iterator::operator!=(const Iterator &other) const
951 {
952 return !(*this == other);
953 }
954
955
956 inline bool
957 Iterator::operator<(const Iterator &other) const
958 {
959 return accessor < other.accessor;
960 }
961
962
963 inline int
964 Iterator::operator-(const Iterator &other) const
965 {
966 (void)other;
970
971 return 0;
972 }
973} // namespace DynamicSparsityPatternIterators
974
975
976inline void
978{
979 // first check the last element (or if line is still empty)
980 if ((entries.empty()) || (entries.back() < j))
981 {
982 entries.push_back(j);
983 return;
984 }
985
986 // do a binary search to find the place where to insert:
987 std::vector<size_type>::iterator it =
988 Utilities::lower_bound(entries.begin(), entries.end(), j);
989
990 // If this entry is a duplicate, exit immediately
991 if (*it == j)
992 return;
993
994 // Insert at the right place in the vector. Vector grows automatically to
995 // fit elements. Always doubles its size.
996 entries.insert(it, j);
997}
998
999
1000
1001inline void
1003{
1006
1007 if (rowset.size() > 0 && !rowset.is_element(i))
1008 return;
1009
1010 have_entries = true;
1011
1012 const size_type rowindex =
1013 rowset.size() == 0 ? i : rowset.index_within_set(i);
1014 lines[rowindex].add(j);
1015}
1016
1017
1018
1019template <typename ForwardIterator>
1020inline void
1022 ForwardIterator begin,
1023 ForwardIterator end,
1024 const bool indices_are_sorted)
1025{
1026 AssertIndexRange(row, rows);
1027
1028 if (rowset.size() > 0 && !rowset.is_element(row))
1029 return;
1030
1031 if (!have_entries && begin < end)
1032 have_entries = true;
1033
1034 const size_type rowindex =
1035 rowset.size() == 0 ? row : rowset.index_within_set(row);
1036 lines[rowindex].add_entries(begin, end, indices_are_sorted);
1037}
1038
1039
1040
1043{
1044 AssertIndexRange(row, n_rows());
1045
1046 if (!have_entries)
1047 return 0;
1048
1049 if (rowset.size() > 0 && !rowset.is_element(row))
1050 return 0;
1051
1052 const size_type rowindex =
1053 rowset.size() == 0 ? row : rowset.index_within_set(row);
1054 return lines[rowindex].entries.size();
1055}
1056
1057
1058
1061 const size_type index) const
1062{
1063 AssertIndexRange(row, n_rows());
1065
1066 const size_type local_row =
1067 rowset.size() != 0u ? rowset.index_within_set(row) : row;
1068 AssertIndexRange(index, lines[local_row].entries.size());
1069 return lines[local_row].entries[index];
1070}
1071
1072
1073
1076{
1077 if (n_rows() > 0)
1078 return begin(0);
1079 else
1080 return end();
1081}
1082
1083
1086{
1087 return {this};
1088}
1089
1090
1091
1094{
1096
1097 if (!have_entries)
1098 return {this};
1099
1100 if (rowset.size() > 0)
1101 {
1102 // We have an IndexSet that describes the locally owned set. For
1103 // performance reasons we need to make sure that we don't do a
1104 // linear search over 0..n_rows(). Instead, find the first entry
1105 // >= row r in the locally owned set (this is done in log
1106 // n_ranges time inside at()). From there, we move forward until
1107 // we find a non-empty row. By iterating over the IndexSet instead
1108 // of incrementing the row index, we potentially skip over entries
1109 // not in the rowset.
1111 if (it == rowset.end())
1112 return end(); // we don't own any row between r and the end
1113
1114 // Instead of using row_length(*it)==0 in the while loop below,
1115 // which involves an expensive index_within_set() call, we
1116 // look at the lines vector directly. This works, because we are
1117 // walking over this vector entry by entry anyways.
1118 size_type rowindex = rowset.index_within_set(*it);
1119
1120 while (it != rowset.end() && lines[rowindex].entries.empty())
1121 {
1122 ++it;
1123 ++rowindex;
1124 }
1125
1126 if (it == rowset.end())
1127 return end();
1128 else
1129 return {this, *it, 0};
1130 }
1131
1132 // Without an index set we have to do a linear search starting at
1133 // row r until we find a non-empty one. We will check the lines vector
1134 // directly instead of going through the slower row_length() function
1135 size_type row = r;
1136
1137 while (row < n_rows() && lines[row].entries.empty())
1138 {
1139 ++row;
1140 }
1141
1142 if (row == n_rows())
1143 return {this};
1144 else
1145 return {this, row, 0};
1146}
1147
1148
1149
1152{
1154
1155 const size_type row = r + 1;
1156 if (row == n_rows())
1157 return {this};
1158 else
1159 return begin(row);
1160}
1161
1162
1163
1164inline const IndexSet &
1166{
1167 return rowset;
1168}
1169
1170
1171
1172inline bool
1177
1178
1180
1181#endif
std::vector< size_type >::const_iterator end_of_row
std::vector< size_type >::const_iterator current_entry
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
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
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false) override
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 reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
void clear_row(const size_type 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)
ElementIterator at(const size_type global_index) const
Definition index_set.cc:873
size_type size() const
Definition index_set.h:1765
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1990
bool is_element(const size_type index) const
Definition index_set.h:1883
ElementIterator end() const
Definition index_set.h:1711
size_type n_rows() const
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type > > &entries)
size_type n_cols() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:494
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & DummyAccessor()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_NOT_IMPLEMENTED()
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition utilities.h:1032
const types::global_dof_index invalid_size_type
Definition types.h:233
unsigned int global_dof_index
Definition types.h:81
void add_entries(ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted)
void add(const size_type col_num)