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
chunk_sparsity_pattern.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2008 - 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_chunk_sparsity_pattern_h
17#define dealii_chunk_sparsity_pattern_h
18
19
20#include <deal.II/base/config.h>
21
24
26
27#include <iostream>
28#include <vector>
29
31
32
33// Forward declaration
34#ifndef DOXYGEN
35template <typename>
37#endif
38
50{
51 // forward declaration
52 class Iterator;
53
64 {
65 public:
70
75
80
86 row() const;
87
91 std::size_t
93
99 column() const;
100
110 bool
112
113
117 bool
118 operator==(const Accessor &) const;
119
120
128 bool
129 operator<(const Accessor &) const;
130
131 protected:
136
141
146
151
155 void
157
158 // Grant access to iterator class.
159 friend class Iterator;
160 };
161
162
163
168 {
169 public:
174
180
184 Iterator &
186
192
196 const Accessor &
197 operator*() const;
198
202 const Accessor *
203 operator->() const;
204
208 bool
209 operator==(const Iterator &) const;
210
214 bool
215 operator!=(const Iterator &) const;
216
224 bool
225 operator<(const Iterator &) const;
226
227 private:
232 };
233} // namespace ChunkSparsityPatternIterators
234
235
236
246{
247public:
257
266
282
289
306
314 const size_type n,
315 const size_type max_chunks_per_row,
316 const size_type chunk_size);
317
326 const size_type n,
327 const std::vector<size_type> &row_lengths,
328 const size_type chunk_size);
329
339 const size_type max_per_row,
340 const size_type chunk_size);
341
350 const std::vector<size_type> &row_lengths,
351 const size_type chunk_size);
352
356 ~ChunkSparsityPattern() override = default;
357
365
374 void
375 reinit(const size_type m,
376 const size_type n,
377 const size_type max_per_row,
378 const size_type chunk_size);
379
394 void
395 reinit(const size_type m,
396 const size_type n,
397 const std::vector<size_type> &row_lengths,
398 const size_type chunk_size);
399
403 void
404 reinit(const size_type m,
405 const size_type n,
406 const ArrayView<const size_type> &row_lengths,
407 const size_type chunk_size);
408
421 void
422 compress();
423
499 template <typename ForwardIterator>
500 void
502 const size_type n_cols,
503 const ForwardIterator begin,
504 const ForwardIterator end,
505 const size_type chunk_size);
506
512 template <typename SparsityPatternType>
513 void
514 copy_from(const SparsityPatternType &dsp, const size_type chunk_size);
515
523 template <typename number>
524 void
525 copy_from(const FullMatrix<number> &matrix, const size_type chunk_size);
526
544 template <typename Sparsity>
545 void
546 create_from(const size_type m,
547 const size_type n,
548 const Sparsity &sparsity_pattern_for_chunks,
549 const size_type chunk_size,
550 const bool optimize_diagonal = true);
551
556 bool
557 empty() const;
558
564
571 max_entries_per_row() const;
572
579 void
580 add(const size_type i, const size_type j);
581
589 void
590 symmetrize();
591
596 inline size_type
597 n_rows() const;
598
603 inline size_type
604 n_cols() const;
605
609 bool
610 exists(const size_type i, const size_type j) const;
611
616 row_length(const size_type row) const;
617
625 bandwidth() const;
626
636 n_nonzero_elements() const;
637
641 bool
643
656 bool
658
665 begin() const;
666
671 end() const;
672
682 begin(const size_type r) const;
683
693 end(const size_type r) const;
694
705 void
706 block_write(std::ostream &out) const;
707
721 void
722 block_read(std::istream &in);
723
729 void
730 print(std::ostream &out) const;
731
745 void
746 print_gnuplot(std::ostream &out) const;
747
752 std::size_t
753 memory_consumption() const;
754
763 size_type,
764 << "The provided number is invalid here: " << arg1);
769 size_type,
770 size_type,
771 << "The given index " << arg1 << " should be less than "
772 << arg2 << '.');
777 size_type,
778 size_type,
779 << "Upon entering a new entry to row " << arg1
780 << ": there was no free entry any more. " << std::endl
781 << "(Maximum number of entries for this row: " << arg2
782 << "; maybe the matrix is already compressed?)");
789 "The operation you attempted is only allowed after the SparsityPattern "
790 "has been set up and compress() was called.");
797 "The operation you attempted changes the structure of the SparsityPattern "
798 "and is not possible after compress() has been called.");
807 size_type,
808 size_type,
809 << "The iterators denote a range of " << arg1
810 << " elements, but the given number of rows was " << arg2);
819 size_type,
820 << "The number of partitions you gave is " << arg1
821 << ", but must be greater than zero.");
826 size_type,
827 size_type,
828 << "The array has size " << arg1 << " but should have size "
829 << arg2);
831private:
836
841
846
852
853 // Make all the chunk sparse matrix kinds friends.
854 template <typename>
855 friend class ChunkSparseMatrix;
856
857 // Make the accessor class a friend.
859};
860
861
863/*---------------------- Inline functions -----------------------------------*/
864
865#ifndef DOXYGEN
866
868{
869 inline Accessor::Accessor(const ChunkSparsityPattern *sparsity_pattern,
870 const size_type row)
871 : sparsity_pattern(sparsity_pattern)
872 , reduced_accessor(row == sparsity_pattern->n_rows() ?
873 *sparsity_pattern->sparsity_pattern.end() :
874 *sparsity_pattern->sparsity_pattern.begin(
875 row / sparsity_pattern->get_chunk_size()))
876 , chunk_row(row == sparsity_pattern->n_rows() ?
877 0 :
878 row % sparsity_pattern->get_chunk_size())
879 , chunk_col(0)
880 {}
881
882
883
884 inline Accessor::Accessor(const ChunkSparsityPattern *sparsity_pattern)
885 : sparsity_pattern(sparsity_pattern)
886 , reduced_accessor(*sparsity_pattern->sparsity_pattern.end())
887 , chunk_row(0)
888 , chunk_col(0)
889 {}
890
891
892
893 inline bool
894 Accessor::is_valid_entry() const
895 {
896 return reduced_accessor.is_valid_entry() &&
897 sparsity_pattern->get_chunk_size() * reduced_accessor.row() +
898 chunk_row <
899 sparsity_pattern->n_rows() &&
900 sparsity_pattern->get_chunk_size() * reduced_accessor.column() +
901 chunk_col <
902 sparsity_pattern->n_cols();
903 }
904
905
906
907 inline Accessor::size_type
908 Accessor::row() const
909 {
910 Assert(is_valid_entry() == true, ExcInvalidIterator());
911
912 return sparsity_pattern->get_chunk_size() * reduced_accessor.row() +
913 chunk_row;
914 }
915
916
917
918 inline Accessor::size_type
919 Accessor::column() const
920 {
921 Assert(is_valid_entry() == true, ExcInvalidIterator());
922
923 return sparsity_pattern->get_chunk_size() * reduced_accessor.column() +
924 chunk_col;
925 }
926
927
928
929 inline std::size_t
930 Accessor::reduced_index() const
931 {
932 Assert(is_valid_entry() == true, ExcInvalidIterator());
933
934 return reduced_accessor.linear_index;
935 }
936
937
938
939 inline bool
940 Accessor::operator==(const Accessor &other) const
941 {
942 // no need to check for equality of sparsity patterns as this is done in
943 // the reduced case already and every ChunkSparsityPattern has its own
944 // reduced sparsity pattern
945 return (reduced_accessor == other.reduced_accessor &&
946 chunk_row == other.chunk_row && chunk_col == other.chunk_col);
947 }
948
949
950
951 inline bool
952 Accessor::operator<(const Accessor &other) const
953 {
954 Assert(sparsity_pattern == other.sparsity_pattern, ExcInternalError());
955
956 if (chunk_row != other.chunk_row)
957 {
958 if (reduced_accessor.linear_index ==
959 reduced_accessor.container->n_nonzero_elements())
960 return false;
961 if (other.reduced_accessor.linear_index ==
962 reduced_accessor.container->n_nonzero_elements())
963 return true;
964
965 const auto global_row = sparsity_pattern->get_chunk_size() *
966 reduced_accessor.row() +
967 chunk_row,
968 other_global_row = sparsity_pattern->get_chunk_size() *
969 other.reduced_accessor.row() +
970 other.chunk_row;
971 if (global_row < other_global_row)
972 return true;
973 else if (global_row > other_global_row)
974 return false;
975 }
976
977 return (
978 reduced_accessor.linear_index < other.reduced_accessor.linear_index ||
979 (reduced_accessor.linear_index == other.reduced_accessor.linear_index &&
980 chunk_col < other.chunk_col));
981 }
982
983
984 inline void
985 Accessor::advance()
986 {
987 const auto chunk_size = sparsity_pattern->get_chunk_size();
988 Assert(chunk_row < chunk_size && chunk_col < chunk_size,
990 Assert(reduced_accessor.row() * chunk_size + chunk_row <
991 sparsity_pattern->n_rows() &&
992 reduced_accessor.column() * chunk_size + chunk_col <
993 sparsity_pattern->n_cols(),
995 if (chunk_size == 1)
996 {
997 reduced_accessor.advance();
998 return;
999 }
1000
1001 ++chunk_col;
1002
1003 // end of chunk
1004 if (chunk_col == chunk_size ||
1005 reduced_accessor.column() * chunk_size + chunk_col ==
1006 sparsity_pattern->n_cols())
1007 {
1008 const auto reduced_row = reduced_accessor.row();
1009 // end of row
1010 if (reduced_accessor.linear_index + 1 ==
1011 reduced_accessor.container->rowstart[reduced_row + 1])
1012 {
1013 ++chunk_row;
1014
1015 chunk_col = 0;
1016
1017 // end of chunk rows or end of matrix
1018 if (chunk_row == chunk_size ||
1019 (reduced_row * chunk_size + chunk_row ==
1020 sparsity_pattern->n_rows()))
1021 {
1022 chunk_row = 0;
1023 reduced_accessor.advance();
1024 }
1025 // go back to the beginning of the same reduced row but with
1026 // chunk_row increased by one
1027 else
1028 reduced_accessor.linear_index =
1029 reduced_accessor.container->rowstart[reduced_row];
1030 }
1031 // advance within chunk
1032 else
1033 {
1034 reduced_accessor.advance();
1035 chunk_col = 0;
1036 }
1037 }
1038 }
1039
1040
1041
1042 inline Iterator::Iterator(const ChunkSparsityPattern *sparsity_pattern,
1043 const size_type row)
1044 : accessor(sparsity_pattern, row)
1045 {}
1046
1047
1048
1049 inline Iterator &
1050 Iterator::operator++()
1051 {
1052 accessor.advance();
1053 return *this;
1054 }
1055
1056
1057
1058 inline Iterator
1059 Iterator::operator++(int)
1060 {
1061 const Iterator iter = *this;
1062 accessor.advance();
1063 return iter;
1064 }
1065
1066
1067
1068 inline const Accessor &
1069 Iterator::operator*() const
1070 {
1071 return accessor;
1072 }
1073
1074
1075
1076 inline const Accessor *
1077 Iterator::operator->() const
1078 {
1079 return &accessor;
1080 }
1081
1082
1083 inline bool
1084 Iterator::operator==(const Iterator &other) const
1085 {
1086 return (accessor == other.accessor);
1087 }
1088
1089
1090
1091 inline bool
1092 Iterator::operator!=(const Iterator &other) const
1093 {
1094 return !(accessor == other.accessor);
1095 }
1096
1097
1098 inline bool
1099 Iterator::operator<(const Iterator &other) const
1100 {
1101 return accessor < other.accessor;
1102 }
1103
1104} // namespace ChunkSparsityPatternIterators
1105
1106
1107
1110{
1111 return {this, 0};
1112}
1113
1114
1117{
1118 return {this, n_rows()};
1119}
1120
1121
1122
1124ChunkSparsityPattern::begin(const size_type r) const
1125{
1127 return {this, r};
1128}
1129
1130
1131
1133ChunkSparsityPattern::end(const size_type r) const
1134{
1136 return {this, r + 1};
1137}
1138
1139
1140
1143{
1144 return rows;
1145}
1146
1147
1150{
1151 return cols;
1152}
1153
1154
1155
1158{
1159 return chunk_size;
1160}
1161
1162
1163
1164inline bool
1166{
1168}
1169
1170
1171
1172template <typename ForwardIterator>
1173void
1174ChunkSparsityPattern::copy_from(const size_type n_rows,
1175 const size_type n_cols,
1176 const ForwardIterator begin,
1177 const ForwardIterator end,
1178 const size_type chunk_size)
1179{
1180 Assert(static_cast<size_type>(std::distance(begin, end)) == n_rows,
1181 ExcIteratorRange(std::distance(begin, end), n_rows));
1182
1183 // first determine row lengths for each row. if the matrix is quadratic,
1184 // then we might have to add an additional entry for the diagonal, if that
1185 // is not yet present. as we have to call compress anyway later on, don't
1186 // bother to check whether that diagonal entry is in a certain row or not
1187 const bool is_square = (n_rows == n_cols);
1188 std::vector<size_type> row_lengths;
1189 row_lengths.reserve(n_rows);
1190 for (ForwardIterator i = begin; i != end; ++i)
1191 row_lengths.push_back(std::distance(i->begin(), i->end()) +
1192 (is_square ? 1 : 0));
1193 reinit(n_rows, n_cols, row_lengths, chunk_size);
1194
1195 // now enter all the elements into the matrix
1196 size_type row = 0;
1197 using inner_iterator =
1198 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1199 for (ForwardIterator i = begin; i != end; ++i, ++row)
1200 {
1201 const inner_iterator end_of_row = i->end();
1202 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1203 {
1204 const size_type col =
1205 internal::SparsityPatternTools::get_column_index_from_iterator(*j);
1206 Assert(col < n_cols, ExcInvalidIndex(col, n_cols));
1207
1208 add(row, col);
1209 }
1210 }
1211
1212 // finally compress everything. this also sorts the entries within each row
1213 compress();
1214}
1215
1216
1217#endif // DOXYGEN
1218
1220
1221#endif
SparsityPatternIterators::Accessor reduced_accessor
Accessor(const ChunkSparsityPattern *matrix, const size_type row)
bool operator<(const Accessor &) const
Accessor(const ChunkSparsityPattern *matrix)
bool operator==(const Accessor &) const
bool operator!=(const Iterator &) const
bool operator==(const Iterator &) const
Iterator(const ChunkSparsityPattern *sp, const size_type row)
bool operator<(const Iterator &) const
const Accessor * operator->() const
const Accessor & operator*() const
void create_from(const size_type m, const size_type n, const Sparsity &sparsity_pattern_for_chunks, const size_type chunk_size, const bool optimize_diagonal=true)
void add(const size_type i, const size_type j)
void block_write(std::ostream &out) const
~ChunkSparsityPattern() override=default
static const size_type invalid_entry
std::size_t memory_consumption() const
types::global_dof_index size_type
bool exists(const size_type i, const size_type j) const
iterator end() const
iterator end(const size_type r) const
void reinit(const size_type m, const size_type n, const size_type max_per_row, const size_type chunk_size)
void print_gnuplot(std::ostream &out) const
void block_read(std::istream &in)
bool is_compressed() const
ChunkSparsityPattern & operator=(const ChunkSparsityPattern &)
size_type max_entries_per_row() const
size_type n_cols() const
size_type n_nonzero_elements() const
size_type row_length(const size_type row) const
iterator begin() const
void print(std::ostream &out) const
size_type get_chunk_size() const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end, const size_type chunk_size)
size_type n_rows() const
iterator begin(const size_type r) const
static constexpr size_type invalid_entry
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define DeclException0(Exception0)
Definition exceptions.h:465
static ::ExceptionBase & ExcInvalidIterator()
static ::ExceptionBase & ExcInvalidNumberOfPartitions(size_type arg1)
static ::ExceptionBase & ExcNotEnoughSpace(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcMETISNotInstalled()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotCompressed()
static ::ExceptionBase & ExcIteratorPastEnd()
static ::ExceptionBase & ExcMatrixIsCompressed()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:533
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:488
static ::ExceptionBase & ExcInvalidArraySize(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcEmptyObject()
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:510
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcIteratorRange(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcInvalidNumber(size_type arg1)
constexpr int chunk_size
Definition cuda_size.h:35
VectorType::value_type * end(VectorType &V)
unsigned int global_dof_index
Definition types.h:82