Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
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
28
29#include <algorithm>
30#include <iostream>
31#include <vector>
32
34
35// Forward declaration
36#ifndef DOXYGEN
38#endif
39
50{
51 // forward declaration
52 class Iterator;
53
58
68 {
69 public:
74 const size_type row,
75 const unsigned int index_within_row);
76
81
87 Accessor();
88
93 row() const;
94
99 index() const;
100
105 column() const;
106
110 bool
111 operator==(const Accessor &) const;
112
120 bool
121 operator<(const Accessor &) const;
122
123 protected:
125 "The instance of this class was initialized"
126 " without DynamicSparsityPattern object, which"
127 " means that it is a dummy accessor that can"
128 " not do any operations.");
129
134
139
144 std::vector<size_type>::const_iterator current_entry;
145
152 std::vector<size_type>::const_iterator end_of_row;
153
157 void
158 advance();
159
160 // Grant access to iterator class.
161 friend class Iterator;
162 };
163
164
165
191 {
192 public:
199 const size_type row,
200 const unsigned int index_within_row);
201
207
213 Iterator() = default;
214
218 Iterator &
219 operator++();
220
225 operator++(int);
226
230 const Accessor &
231 operator*() const;
232
236 const Accessor *
237 operator->() const;
238
242 bool
243 operator==(const Iterator &) const;
244
248 bool
249 operator!=(const Iterator &) const;
250
258 bool
259 operator<(const Iterator &) const;
260
267 int
268 operator-(const Iterator &p) const;
269
270 private:
275 };
276} // namespace DynamicSparsityPatternIterators
277
278
323{
324public:
329
338
344
351
362
370 const size_type n,
371 const IndexSet &rowset = IndexSet());
372
378 DynamicSparsityPattern(const IndexSet &indexset);
379
384
392
399 void
400 reinit(const size_type m,
401 const size_type n,
402 const IndexSet &rowset = IndexSet());
403
409 void
410 compress();
411
416 bool
417 empty() const;
418
424 max_entries_per_row() const;
425
429 void
430 add(const size_type i, const size_type j);
431
436 template <typename ForwardIterator>
437 void
438 add_entries(const size_type row,
439 ForwardIterator begin,
440 ForwardIterator end,
441 const bool indices_are_unique_and_sorted = false);
442
443 virtual void
444 add_row_entries(const size_type & row,
445 const ArrayView<const size_type> &columns,
446 const bool indices_are_sorted = false) override;
447
449
453 bool
454 exists(const size_type i, const size_type j) const;
455
463 get_view(const IndexSet &rows) const;
464
472 void
473 symmetrize();
474
479 template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
480 void
481 compute_mmult_pattern(const SparsityPatternTypeLeft & left,
482 const SparsityPatternTypeRight &right);
483
488 template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
489 void
490 compute_Tmmult_pattern(const SparsityPatternTypeLeft & left,
491 const SparsityPatternTypeRight &right);
492
498 void
499 print(std::ostream &out) const;
500
514 void
515 print_gnuplot(std::ostream &out) const;
516
522 row_length(const size_type row) const;
523
527 void
528 clear_row(const size_type row);
529
535 column_number(const size_type row, const size_type index) const;
536
543 column_index(const size_type row, const size_type col) const;
544
564 begin() const;
565
570 end() const;
571
589 begin(const size_type r) const;
590
600 end(const size_type r) const;
601
612 bandwidth() const;
613
619 n_nonzero_elements() const;
620
626 const IndexSet &
627 row_index_set() const;
628
636 nonempty_cols() const;
637
645 nonempty_rows() const;
646
657 static bool
659
665 memory_consumption() const;
666
667private:
672
678
679
686 struct Line
687 {
688 public:
693 std::vector<size_type> entries;
694
698 void
699 add(const size_type col_num);
700
704 template <typename ForwardIterator>
705 void
706 add_entries(ForwardIterator begin,
707 ForwardIterator end,
708 const bool indices_are_sorted);
709
714 memory_consumption() const;
715 };
716
717
721 std::vector<Line> lines;
722
723 // make the accessor class a friend
725};
726
728/*---------------------- Inline functions -----------------------------------*/
729
730
732{
733 inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern,
734 const size_type row,
735 const unsigned int index_within_row)
736 : sparsity_pattern(sparsity_pattern)
737 , current_row(row)
738 , current_entry(
739 ((sparsity_pattern->rowset.size() == 0) ?
740 sparsity_pattern->lines[current_row].entries.begin() :
741 sparsity_pattern
742 ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
743 .entries.begin()) +
744 index_within_row)
745 , end_of_row(
746 (sparsity_pattern->rowset.size() == 0) ?
747 sparsity_pattern->lines[current_row].entries.end() :
748 sparsity_pattern
749 ->lines[sparsity_pattern->rowset.index_within_set(current_row)]
750 .entries.end())
751 {
755 ExcMessage("You can't create an iterator into a "
756 "DynamicSparsityPattern's row that is not "
757 "actually stored by that sparsity pattern "
758 "based on the IndexSet argument to it."));
760 index_within_row,
761 ((sparsity_pattern->rowset.size() == 0) ?
762 sparsity_pattern->lines[current_row].entries.size() :
765 .entries.size()));
766 }
767
768
769 inline Accessor::Accessor(const DynamicSparsityPattern *sparsity_pattern)
770 : sparsity_pattern(sparsity_pattern)
771 , current_row(numbers::invalid_size_type)
772 , current_entry()
773 , end_of_row()
774 {}
775
776
777
779 : sparsity_pattern(nullptr)
780 , current_row(numbers::invalid_size_type)
781 , current_entry()
782 , end_of_row()
783 {}
784
785
786 inline size_type
788 {
790 Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
791
792 return current_row;
793 }
794
795
796 inline size_type
798 {
800 Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
801
802 return *current_entry;
803 }
804
805
806 inline size_type
808 {
810 Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
811
812 return (current_entry -
813 ((sparsity_pattern->rowset.size() == 0) ?
814 sparsity_pattern->lines[current_row].entries.begin() :
817 .entries.begin()));
818 }
819
820
821
822 inline bool
823 Accessor::operator==(const Accessor &other) const
824 {
826 Assert(other.sparsity_pattern != nullptr, DummyAccessor());
827 // compare the sparsity pattern the iterator points into, the
828 // current row, and the location within this row. ignore the
829 // latter if the row is past-the-end because in that case the
830 // current_entry field may not point to a deterministic location
831 return (sparsity_pattern == other.sparsity_pattern &&
832 current_row == other.current_row &&
834 (current_entry == other.current_entry)));
835 }
836
837
838
839 inline bool
840 Accessor::operator<(const Accessor &other) const
841 {
843 Assert(other.sparsity_pattern != nullptr, DummyAccessor());
845
846 // if *this is past-the-end, then it is less than no one
848 return (false);
849 // now *this should be an valid value
850 Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
851
852 // if other is past-the-end
854 return (true);
855 // now other should be an valid value
857
858 // both iterators are not one-past-the-end
859 return ((current_row < other.current_row) ||
860 ((current_row == other.current_row) &&
861 (current_entry < other.current_entry)));
862 }
863
864
865 inline void
867 {
869 Assert(current_row < sparsity_pattern->n_rows(), ExcInternalError());
870
871 // move to the next element in this row
873
874 // if this moves us beyond the end of the row, go to the next row
875 // if possible, or set the iterator to an invalid state if not.
876 //
877 // going to the next row is a bit complicated because we may have
878 // to skip over empty rows, and because we also have to avoid rows
879 // that aren't listed in a possibly passed IndexSet argument of
880 // the sparsity pattern. consequently, rather than trying to
881 // duplicate code here, just call the begin() function of the
882 // sparsity pattern itself
884 {
886 *this = *sparsity_pattern->begin(current_row + 1);
887 else
888 *this = Accessor(sparsity_pattern); // invalid object
889 }
890 }
891
892
893
894 inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern,
895 const size_type row,
896 const unsigned int index_within_row)
897 : accessor(sparsity_pattern, row, index_within_row)
898 {}
899
900
901
902 inline Iterator::Iterator(const DynamicSparsityPattern *sparsity_pattern)
903 : accessor(sparsity_pattern)
904 {}
905
906
907
908 inline Iterator &
910 {
912 return *this;
913 }
914
915
916
917 inline Iterator
919 {
920 const Iterator iter = *this;
922 return iter;
923 }
924
925
926
927 inline const Accessor &
929 {
930 return accessor;
931 }
932
933
934
935 inline const Accessor *
937 {
938 return &accessor;
939 }
940
941
942 inline bool
943 Iterator::operator==(const Iterator &other) const
944 {
945 return (accessor == other.accessor);
946 }
947
948
949
950 inline bool
951 Iterator::operator!=(const Iterator &other) const
952 {
953 return !(*this == other);
954 }
955
956
957 inline bool
958 Iterator::operator<(const Iterator &other) const
959 {
960 return accessor < other.accessor;
961 }
962
963
964 inline int
965 Iterator::operator-(const Iterator &other) const
966 {
967 (void)other;
970 Assert(false, ExcNotImplemented());
971
972 return 0;
973 }
974} // namespace DynamicSparsityPatternIterators
975
976
977inline void
979{
980 // first check the last element (or if line is still empty)
981 if ((entries.size() == 0) || (entries.back() < j))
982 {
983 entries.push_back(j);
984 return;
985 }
986
987 // do a binary search to find the place where to insert:
988 std::vector<size_type>::iterator it =
989 Utilities::lower_bound(entries.begin(), entries.end(), j);
990
991 // If this entry is a duplicate, exit immediately
992 if (*it == j)
993 return;
994
995 // Insert at the right place in the vector. Vector grows automatically to
996 // fit elements. Always doubles its size.
997 entries.insert(it, j);
998}
999
1000
1001
1002inline void
1004{
1007
1008 if (rowset.size() > 0 && !rowset.is_element(i))
1009 return;
1010
1011 have_entries = true;
1012
1013 const size_type rowindex =
1014 rowset.size() == 0 ? i : rowset.index_within_set(i);
1015 lines[rowindex].add(j);
1016}
1017
1018
1019
1020template <typename ForwardIterator>
1021inline void
1023 ForwardIterator begin,
1024 ForwardIterator end,
1025 const bool indices_are_sorted)
1026{
1027 AssertIndexRange(row, rows);
1028
1029 if (rowset.size() > 0 && !rowset.is_element(row))
1030 return;
1031
1032 if (!have_entries && begin < end)
1033 have_entries = true;
1034
1035 const size_type rowindex =
1036 rowset.size() == 0 ? row : rowset.index_within_set(row);
1037 lines[rowindex].add_entries(begin, end, indices_are_sorted);
1038}
1039
1040
1041
1044{
1045 AssertIndexRange(row, n_rows());
1046
1047 if (!have_entries)
1048 return 0;
1049
1050 if (rowset.size() > 0 && !rowset.is_element(row))
1051 return 0;
1052
1053 const size_type rowindex =
1054 rowset.size() == 0 ? row : rowset.index_within_set(row);
1055 return lines[rowindex].entries.size();
1056}
1057
1058
1059
1062 const size_type index) const
1063{
1064 AssertIndexRange(row, n_rows());
1066
1067 const size_type local_row =
1068 rowset.size() != 0u ? rowset.index_within_set(row) : row;
1069 AssertIndexRange(index, lines[local_row].entries.size());
1070 return lines[local_row].entries[index];
1071}
1072
1073
1074
1077{
1078 if (n_rows() > 0)
1079 return begin(0);
1080 else
1081 return end();
1082}
1083
1084
1087{
1088 return {this};
1089}
1090
1091
1092
1095{
1097
1098 if (!have_entries)
1099 return {this};
1100
1101 if (rowset.size() > 0)
1102 {
1103 // We have an IndexSet that describes the locally owned set. For
1104 // performance reasons we need to make sure that we don't do a
1105 // linear search over 0..n_rows(). Instead, find the first entry
1106 // >= row r in the locally owned set (this is done in log
1107 // n_ranges time inside at()). From there, we move forward until
1108 // we find a non-empty row. By iterating over the IndexSet instead
1109 // of incrementing the row index, we potentially skip over entries
1110 // not in the rowset.
1112 if (it == rowset.end())
1113 return end(); // we don't own any row between r and the end
1114
1115 // Instead of using row_length(*it)==0 in the while loop below,
1116 // which involves an expensive index_within_set() call, we
1117 // look at the lines vector directly. This works, because we are
1118 // walking over this vector entry by entry anyways.
1119 size_type rowindex = rowset.index_within_set(*it);
1120
1121 while (it != rowset.end() && lines[rowindex].entries.size() == 0)
1122 {
1123 ++it;
1124 ++rowindex;
1125 }
1126
1127 if (it == rowset.end())
1128 return end();
1129 else
1130 return {this, *it, 0};
1131 }
1132
1133 // Without an index set we have to do a linear search starting at
1134 // row r until we find a non-empty one. We will check the lines vector
1135 // directly instead of going through the slower row_length() function
1136 size_type row = r;
1137
1138 while (row < n_rows() && lines[row].entries.size() == 0)
1139 {
1140 ++row;
1141 }
1142
1143 if (row == n_rows())
1144 return {this};
1145 else
1146 return {this, row, 0};
1147}
1148
1149
1150
1153{
1155
1156 const size_type row = r + 1;
1157 if (row == n_rows())
1158 return {this};
1159 else
1160 return begin(row);
1161}
1162
1163
1164
1165inline const IndexSet &
1167{
1168 return rowset;
1169}
1170
1171
1172
1173inline bool
1175{
1176 return true;
1177}
1178
1179
1181
1182#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:640
size_type size() const
Definition index_set.h:1661
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1883
bool is_element(const size_type index) const
Definition index_set.h:1776
ElementIterator end() const
Definition index_set.h:1607
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:488
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & DummyAccessor()
static ::ExceptionBase & ExcMessage(std::string arg1)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition utilities.h:1016
const types::global_dof_index invalid_size_type
Definition types.h:222
unsigned int global_dof_index
Definition types.h:82
void add_entries(ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted)
void add(const size_type col_num)