Reference documentation for deal.II version GIT 6bdf9a4b6c 2022-08-12 19:20: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\}}\)
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
35 template <typename>
36 class ChunkSparseMatrix;
37 #endif
38 
50 {
51  // forward declaration
52  class Iterator;
53 
63  class Accessor
64  {
65  public:
70 
75 
80 
85  size_type
86  row() const;
87 
91  std::size_t
92  reduced_index() const;
93 
98  size_type
99  column() const;
100 
110  bool
111  is_valid_entry() const;
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 
167  class Iterator
168  {
169  public:
174 
179  Iterator(const ChunkSparsityPattern *sp, const size_type row);
180 
184  Iterator &
186 
190  Iterator
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 {
247 public:
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
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 
562  size_type
563  get_chunk_size() const;
564 
570  size_type
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 
615  size_type
616  row_length(const size_type row) const;
617 
624  size_type
625  bandwidth() const;
626 
635  size_type
636  n_nonzero_elements() const;
637 
641  bool
642  is_compressed() const;
643 
656  bool
658 
664  iterator
665  begin() const;
666 
670  iterator
671  end() const;
672 
681  iterator
682  begin(const size_type r) const;
683 
692  iterator
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);
831 private:
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
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 &
1051  {
1052  accessor.advance();
1053  return *this;
1054  }
1055 
1056 
1057 
1058  inline Iterator
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 
1125 {
1126  AssertIndexRange(r, n_rows());
1127  return {this, r};
1128 }
1129 
1130 
1131 
1133 ChunkSparsityPattern::end(const size_type r) const
1134 {
1135  AssertIndexRange(r, n_rows());
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 
1164 inline bool
1166 {
1168 }
1169 
1170 
1171 
1172 template <typename ForwardIterator>
1173 void
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 =
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
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
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
const ChunkSparsityPattern * sparsity_pattern
bool operator!=(const Iterator &) const
bool operator==(const Iterator &) const
const Accessor * operator->() const
const Accessor & operator*() const
Iterator(const ChunkSparsityPattern *sp, const size_type row)
bool operator<(const Iterator &) 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
SparsityPattern sparsity_pattern
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 const size_type invalid_entry
std::enable_if_t< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type > operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcMETISNotInstalled()
#define DeclException0(Exception0)
Definition: exceptions.h:464
static ::ExceptionBase & ExcInvalidArraySize(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcNotCompressed()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotEnoughSpace(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcIteratorRange(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcIteratorPastEnd()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcInvalidIterator()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
static ::ExceptionBase & ExcInvalidNumberOfPartitions(size_type arg1)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcEmptyObject()
static ::ExceptionBase & ExcMatrixIsCompressed()
static ::ExceptionBase & ExcInvalidNumber(size_type arg1)
constexpr int chunk_size
Definition: cuda_size.h:35
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
Definition: cuda_kernels.h:45
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
size_type get_column_index_from_iterator(const size_type i)
unsigned int global_dof_index
Definition: types.h:76
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
void advance(std::tuple< I1, I2 > &t, const unsigned int n)