Reference documentation for deal.II version 9.3.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\}}\)
chunk_sparsity_pattern.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2020 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 
49 {
50  // forward declaration
51  class Iterator;
52 
62  class Accessor
63  {
64  public:
69 
74 
79 
84  size_type
85  row() const;
86 
90  std::size_t
91  reduced_index() const;
92 
97  size_type
98  column() const;
99 
109  bool
110  is_valid_entry() const;
111 
112 
116  bool
117  operator==(const Accessor &) const;
118 
119 
127  bool
128  operator<(const Accessor &) const;
129 
130  protected:
135 
140 
145 
150 
154  void
155  advance();
156 
157  // Grant access to iterator class.
158  friend class Iterator;
159  };
160 
161 
162 
166  class Iterator
167  {
168  public:
173 
178  Iterator(const ChunkSparsityPattern *sp, const size_type row);
179 
183  Iterator &
184  operator++();
185 
189  Iterator
190  operator++(int);
191 
195  const Accessor &operator*() const;
196 
200  const Accessor *operator->() const;
201 
205  bool
206  operator==(const Iterator &) const;
207 
211  bool
212  operator!=(const Iterator &) const;
213 
221  bool
222  operator<(const Iterator &) const;
223 
224  private:
229  };
230 } // namespace ChunkSparsityPatternIterators
231 
232 
233 
243 {
244 public:
254 
263 
278  static const size_type invalid_entry = SparsityPattern::invalid_entry;
279 
286 
303 
311  const size_type n,
312  const size_type max_chunks_per_row,
313  const size_type chunk_size);
314 
323  const size_type n,
324  const std::vector<size_type> &row_lengths,
325  const size_type chunk_size);
326 
336  const size_type max_per_row,
337  const size_type chunk_size);
338 
347  const std::vector<size_type> &row_lengths,
348  const size_type chunk_size);
349 
353  ~ChunkSparsityPattern() override = default;
354 
361  operator=(const ChunkSparsityPattern &);
362 
371  void
372  reinit(const size_type m,
373  const size_type n,
374  const size_type max_per_row,
375  const size_type chunk_size);
376 
391  void
392  reinit(const size_type m,
393  const size_type n,
394  const std::vector<size_type> &row_lengths,
395  const size_type chunk_size);
396 
400  void
401  reinit(const size_type m,
402  const size_type n,
403  const ArrayView<const size_type> &row_lengths,
404  const size_type chunk_size);
405 
418  void
419  compress();
420 
496  template <typename ForwardIterator>
497  void
498  copy_from(const size_type n_rows,
499  const size_type n_cols,
500  const ForwardIterator begin,
501  const ForwardIterator end,
502  const size_type chunk_size);
503 
509  template <typename SparsityPatternType>
510  void
511  copy_from(const SparsityPatternType &dsp, const size_type chunk_size);
512 
520  template <typename number>
521  void
522  copy_from(const FullMatrix<number> &matrix, const size_type chunk_size);
523 
541  template <typename Sparsity>
542  void
543  create_from(const size_type m,
544  const size_type n,
545  const Sparsity &sparsity_pattern_for_chunks,
546  const size_type chunk_size,
547  const bool optimize_diagonal = true);
548 
553  bool
554  empty() const;
555 
559  size_type
560  get_chunk_size() const;
561 
567  size_type
568  max_entries_per_row() const;
569 
576  void
577  add(const size_type i, const size_type j);
578 
586  void
587  symmetrize();
588 
593  inline size_type
594  n_rows() const;
595 
600  inline size_type
601  n_cols() const;
602 
606  bool
607  exists(const size_type i, const size_type j) const;
608 
612  size_type
613  row_length(const size_type row) const;
614 
621  size_type
622  bandwidth() const;
623 
632  size_type
633  n_nonzero_elements() const;
634 
638  bool
639  is_compressed() const;
640 
653  bool
654  stores_only_added_elements() const;
655 
661  iterator
662  begin() const;
663 
667  iterator
668  end() const;
669 
678  iterator
679  begin(const size_type r) const;
680 
689  iterator
690  end(const size_type r) const;
691 
702  void
703  block_write(std::ostream &out) const;
704 
718  void
719  block_read(std::istream &in);
720 
726  void
727  print(std::ostream &out) const;
728 
742  void
743  print_gnuplot(std::ostream &out) const;
744 
749  std::size_t
750  memory_consumption() const;
751 
760  size_type,
761  << "The provided number is invalid here: " << arg1);
765  DeclException2(ExcInvalidIndex,
766  size_type,
767  size_type,
768  << "The given index " << arg1 << " should be less than "
769  << arg2 << ".");
773  DeclException2(ExcNotEnoughSpace,
774  size_type,
775  size_type,
776  << "Upon entering a new entry to row " << arg1
777  << ": there was no free entry any more. " << std::endl
778  << "(Maximum number of entries for this row: " << arg2
779  << "; maybe the matrix is already compressed?)");
785  ExcNotCompressed,
786  "The operation you attempted is only allowed after the SparsityPattern "
787  "has been set up and compress() was called.");
793  ExcMatrixIsCompressed,
794  "The operation you attempted changes the structure of the SparsityPattern "
795  "and is not possible after compress() has been called.");
803  DeclException2(ExcIteratorRange,
804  size_type,
805  size_type,
806  << "The iterators denote a range of " << arg1
807  << " elements, but the given number of rows was " << arg2);
816  size_type,
817  << "The number of partitions you gave is " << arg1
818  << ", but must be greater than zero.");
823  size_type,
824  size_type,
825  << "The array has size " << arg1 << " but should have size "
826  << arg2);
828 private:
833 
838 
843 
849 
850  // Make all the chunk sparse matrix kinds friends.
851  template <typename>
852  friend class ChunkSparseMatrix;
853 
854  // Make the accessor class a friend.
856 };
857 
858 
860 /*---------------------- Inline functions -----------------------------------*/
861 
862 #ifndef DOXYGEN
863 
865 {
867  const size_type row)
868  : sparsity_pattern(sparsity_pattern)
869  , reduced_accessor(row == sparsity_pattern->n_rows() ?
870  *sparsity_pattern->sparsity_pattern.end() :
871  *sparsity_pattern->sparsity_pattern.begin(
872  row / sparsity_pattern->get_chunk_size()))
873  , chunk_row(row == sparsity_pattern->n_rows() ?
874  0 :
875  row % sparsity_pattern->get_chunk_size())
876  , chunk_col(0)
877  {}
878 
879 
880 
881  inline Accessor::Accessor(const ChunkSparsityPattern *sparsity_pattern)
882  : sparsity_pattern(sparsity_pattern)
883  , reduced_accessor(*sparsity_pattern->sparsity_pattern.end())
884  , chunk_row(0)
885  , chunk_col(0)
886  {}
887 
888 
889 
890  inline bool
891  Accessor::is_valid_entry() const
892  {
893  return reduced_accessor.is_valid_entry() &&
894  sparsity_pattern->get_chunk_size() * reduced_accessor.row() +
895  chunk_row <
897  sparsity_pattern->get_chunk_size() * reduced_accessor.column() +
898  chunk_col <
900  }
901 
902 
903 
904  inline Accessor::size_type
905  Accessor::row() const
906  {
907  Assert(is_valid_entry() == true, ExcInvalidIterator());
908 
909  return sparsity_pattern->get_chunk_size() * reduced_accessor.row() +
910  chunk_row;
911  }
912 
913 
914 
915  inline Accessor::size_type
916  Accessor::column() const
917  {
918  Assert(is_valid_entry() == true, ExcInvalidIterator());
919 
920  return sparsity_pattern->get_chunk_size() * reduced_accessor.column() +
921  chunk_col;
922  }
923 
924 
925 
926  inline std::size_t
927  Accessor::reduced_index() const
928  {
929  Assert(is_valid_entry() == true, ExcInvalidIterator());
930 
931  return reduced_accessor.linear_index;
932  }
933 
934 
935 
936  inline bool
937  Accessor::operator==(const Accessor &other) const
938  {
939  // no need to check for equality of sparsity patterns as this is done in
940  // the reduced case already and every ChunkSparsityPattern has its own
941  // reduced sparsity pattern
942  return (reduced_accessor == other.reduced_accessor &&
943  chunk_row == other.chunk_row && chunk_col == other.chunk_col);
944  }
945 
946 
947 
948  inline bool
949  Accessor::operator<(const Accessor &other) const
950  {
951  Assert(sparsity_pattern == other.sparsity_pattern, ExcInternalError());
952 
953  if (chunk_row != other.chunk_row)
954  {
955  if (reduced_accessor.linear_index ==
956  reduced_accessor.container->n_nonzero_elements())
957  return false;
958  if (other.reduced_accessor.linear_index ==
959  reduced_accessor.container->n_nonzero_elements())
960  return true;
961 
962  const auto global_row = sparsity_pattern->get_chunk_size() *
963  reduced_accessor.row() +
964  chunk_row,
965  other_global_row = sparsity_pattern->get_chunk_size() *
966  other.reduced_accessor.row() +
967  other.chunk_row;
968  if (global_row < other_global_row)
969  return true;
970  else if (global_row > other_global_row)
971  return false;
972  }
973 
974  return (
975  reduced_accessor.linear_index < other.reduced_accessor.linear_index ||
976  (reduced_accessor.linear_index == other.reduced_accessor.linear_index &&
977  chunk_col < other.chunk_col));
978  }
979 
980 
981  inline void
983  {
984  const auto chunk_size = sparsity_pattern->get_chunk_size();
985  Assert(chunk_row < chunk_size && chunk_col < chunk_size,
987  Assert(reduced_accessor.row() * chunk_size + chunk_row <
989  reduced_accessor.column() * chunk_size + chunk_col <
992  if (chunk_size == 1)
993  {
994  reduced_accessor.advance();
995  return;
996  }
997 
998  ++chunk_col;
999 
1000  // end of chunk
1001  if (chunk_col == chunk_size ||
1002  reduced_accessor.column() * chunk_size + chunk_col ==
1004  {
1005  const auto reduced_row = reduced_accessor.row();
1006  // end of row
1007  if (reduced_accessor.linear_index + 1 ==
1008  reduced_accessor.container->rowstart[reduced_row + 1])
1009  {
1010  ++chunk_row;
1011 
1012  chunk_col = 0;
1013 
1014  // end of chunk rows or end of matrix
1015  if (chunk_row == chunk_size ||
1016  (reduced_row * chunk_size + chunk_row ==
1018  {
1019  chunk_row = 0;
1020  reduced_accessor.advance();
1021  }
1022  // go back to the beginning of the same reduced row but with
1023  // chunk_row increased by one
1024  else
1025  reduced_accessor.linear_index =
1026  reduced_accessor.container->rowstart[reduced_row];
1027  }
1028  // advance within chunk
1029  else
1030  {
1031  reduced_accessor.advance();
1032  chunk_col = 0;
1033  }
1034  }
1035  }
1036 
1037 
1038 
1039  inline Iterator::Iterator(const ChunkSparsityPattern *sparsity_pattern,
1040  const size_type row)
1041  : accessor(sparsity_pattern, row)
1042  {}
1043 
1044 
1045 
1046  inline Iterator &
1048  {
1049  accessor.advance();
1050  return *this;
1051  }
1052 
1053 
1054 
1055  inline Iterator
1057  {
1058  const Iterator iter = *this;
1059  accessor.advance();
1060  return iter;
1061  }
1062 
1063 
1064 
1065  inline const Accessor &Iterator::operator*() const
1066  {
1067  return accessor;
1068  }
1069 
1070 
1071 
1072  inline const Accessor *Iterator::operator->() const
1073  {
1074  return &accessor;
1075  }
1076 
1077 
1078  inline bool
1079  Iterator::operator==(const Iterator &other) const
1080  {
1081  return (accessor == other.accessor);
1082  }
1083 
1084 
1085 
1086  inline bool
1087  Iterator::operator!=(const Iterator &other) const
1088  {
1089  return !(accessor == other.accessor);
1090  }
1091 
1092 
1093  inline bool
1094  Iterator::operator<(const Iterator &other) const
1095  {
1096  return accessor < other.accessor;
1097  }
1098 
1099 } // namespace ChunkSparsityPatternIterators
1100 
1101 
1102 
1105 {
1106  return {this, 0};
1107 }
1108 
1109 
1112 {
1113  return {this, n_rows()};
1114 }
1115 
1116 
1117 
1120 {
1121  AssertIndexRange(r, n_rows());
1122  return {this, r};
1123 }
1124 
1125 
1126 
1128 ChunkSparsityPattern::end(const size_type r) const
1129 {
1130  AssertIndexRange(r, n_rows());
1131  return {this, r + 1};
1132 }
1133 
1134 
1135 
1138 {
1139  return rows;
1140 }
1141 
1142 
1145 {
1146  return cols;
1147 }
1148 
1149 
1150 
1153 {
1154  return chunk_size;
1155 }
1156 
1157 
1158 
1159 inline bool
1161 {
1163 }
1164 
1165 
1166 
1167 template <typename ForwardIterator>
1168 void
1170  const size_type n_cols,
1171  const ForwardIterator begin,
1172  const ForwardIterator end,
1173  const size_type chunk_size)
1174 {
1175  Assert(static_cast<size_type>(std::distance(begin, end)) == n_rows,
1176  ExcIteratorRange(std::distance(begin, end), n_rows));
1177 
1178  // first determine row lengths for each row. if the matrix is quadratic,
1179  // then we might have to add an additional entry for the diagonal, if that
1180  // is not yet present. as we have to call compress anyway later on, don't
1181  // bother to check whether that diagonal entry is in a certain row or not
1182  const bool is_square = (n_rows == n_cols);
1183  std::vector<size_type> row_lengths;
1184  row_lengths.reserve(n_rows);
1185  for (ForwardIterator i = begin; i != end; ++i)
1186  row_lengths.push_back(std::distance(i->begin(), i->end()) +
1187  (is_square ? 1 : 0));
1188  reinit(n_rows, n_cols, row_lengths, chunk_size);
1189 
1190  // now enter all the elements into the matrix
1191  size_type row = 0;
1192  using inner_iterator =
1193  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1194  for (ForwardIterator i = begin; i != end; ++i, ++row)
1195  {
1196  const inner_iterator end_of_row = i->end();
1197  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1198  {
1199  const size_type col =
1201  Assert(col < n_cols, ExcInvalidIndex(col, n_cols));
1202 
1203  add(row, col);
1204  }
1205  }
1206 
1207  // finally compress everything. this also sorts the entries within each row
1208  compress();
1209 }
1210 
1211 
1212 #endif // DOXYGEN
1213 
1215 
1216 #endif
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)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
Contents is actually a matrix.
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
types::global_dof_index size_type
Definition: cuda_kernels.h:45
types::global_dof_index size_type
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
static ::ExceptionBase & ExcMETISNotInstalled()
iterator end() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
const ChunkSparsityPattern * sparsity_pattern
void add(const size_type i, const size_type j)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
SparsityPattern sparsity_pattern
static ::ExceptionBase & ExcInvalidIterator()
size_type get_column_index_from_iterator(const size_type i)
size_type n_cols() const
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
std::string compress(const std::string &input)
Definition: utilities.cc:392
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
Accessor(const ChunkSparsityPattern *matrix, const size_type row)
bool is_compressed() const
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcInvalidArraySize(int arg1, int arg2)
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
#define DeclException0(Exception0)
Definition: exceptions.h:470
static ::ExceptionBase & ExcInvalidNumber(unsigned int arg1)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
VectorType::value_type * end(VectorType &V)
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
size_type n_rows() const
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcIteratorPastEnd()
unsigned int global_dof_index
Definition: types.h:76
static ::ExceptionBase & ExcIteratorRange(size_type arg1, size_type arg2)
void advance(std::tuple< I1, I2 > &t, const unsigned int n)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
VectorType::value_type * begin(VectorType &V)
constexpr int chunk_size
Definition: cuda_size.h:35
static ::ExceptionBase & ExcEmptyObject()
iterator begin() const
static const size_type invalid_entry
bool operator<(const Accessor &) const
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
SparsityPatternIterators::Accessor reduced_accessor
void reinit(const size_type m, const size_type n, const size_type max_per_row, const size_type chunk_size)
size_type n_rows() const
bool operator==(const Accessor &) const
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()
size_type n_cols() const