Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
chunk_sparsity_pattern.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2019 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 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/subscriptor.h>
24 
25 #include <deal.II/lac/sparsity_pattern.h>
26 
27 #include <iostream>
28 #include <vector>
29 
30 DEAL_II_NAMESPACE_OPEN
31 
32 
33 template <typename>
34 class ChunkSparseMatrix;
35 
36 
47 {
48  // forward declaration
49  class Iterator;
50 
63  class Accessor
64  {
65  public:
70 
74  Accessor(const ChunkSparsityPattern *matrix, const size_type row);
75 
79  Accessor(const ChunkSparsityPattern *matrix);
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
156  advance();
157 
161  friend class Iterator;
162  };
163 
164 
165 
169  class Iterator
170  {
171  public:
176 
181  Iterator(const ChunkSparsityPattern *sp, const size_type row);
182 
186  Iterator &
187  operator++();
188 
192  Iterator
193  operator++(int);
194 
198  const Accessor &operator*() const;
199 
203  const Accessor *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 
248 {
249 public:
259 
268 
284 
291 
308 
316  const size_type n,
317  const size_type max_chunks_per_row,
318  const size_type chunk_size);
319 
328  const size_type n,
329  const std::vector<size_type> &row_lengths,
330  const size_type chunk_size);
331 
341  const size_type max_per_row,
342  const size_type chunk_size);
343 
352  const std::vector<size_type> &row_lengths,
353  const size_type chunk_size);
354 
358  ~ChunkSparsityPattern() override = default;
359 
367 
376  void
377  reinit(const size_type m,
378  const size_type n,
379  const size_type max_per_row,
380  const size_type chunk_size);
381 
396  void
397  reinit(const size_type m,
398  const size_type n,
399  const std::vector<size_type> &row_lengths,
400  const size_type chunk_size);
401 
405  void
406  reinit(const size_type m,
407  const size_type n,
408  const ArrayView<const size_type> &row_lengths,
409  const size_type chunk_size);
410 
423  void
424  compress();
425 
501  template <typename ForwardIterator>
502  void
503  copy_from(const size_type n_rows,
504  const size_type n_cols,
505  const ForwardIterator begin,
506  const ForwardIterator end,
507  const size_type chunk_size);
508 
514  template <typename SparsityPatternType>
515  void
516  copy_from(const SparsityPatternType &dsp, const size_type chunk_size);
517 
525  template <typename number>
526  void
527  copy_from(const FullMatrix<number> &matrix, const size_type chunk_size);
528 
546  template <typename Sparsity>
547  void
548  create_from(const size_type m,
549  const size_type n,
550  const Sparsity &sparsity_pattern_for_chunks,
551  const size_type chunk_size,
552  const bool optimize_diagonal = true);
553 
558  bool
559  empty() const;
560 
564  size_type
565  get_chunk_size() const;
566 
572  size_type
573  max_entries_per_row() const;
574 
581  void
582  add(const size_type i, const size_type j);
583 
591  void
592  symmetrize();
593 
598  inline size_type
599  n_rows() const;
600 
605  inline size_type
606  n_cols() const;
607 
611  bool
612  exists(const size_type i, const size_type j) const;
613 
617  size_type
618  row_length(const size_type row) const;
619 
626  size_type
627  bandwidth() const;
628 
637  size_type
638  n_nonzero_elements() const;
639 
643  bool
644  is_compressed() const;
645 
658  bool
660 
666  iterator
667  begin() const;
668 
672  iterator
673  end() const;
674 
683  iterator
684  begin(const size_type r) const;
685 
694  iterator
695  end(const size_type r) const;
696 
707  void
708  block_write(std::ostream &out) const;
709 
723  void
724  block_read(std::istream &in);
725 
731  void
732  print(std::ostream &out) const;
733 
747  void
748  print_gnuplot(std::ostream &out) const;
749 
754  std::size_t
755  memory_consumption() const;
756 
765  size_type,
766  << "The provided number is invalid here: " << arg1);
771  size_type,
772  size_type,
773  << "The given index " << arg1 << " should be less than "
774  << arg2 << ".");
779  size_type,
780  size_type,
781  << "Upon entering a new entry to row " << arg1
782  << ": there was no free entry any more. " << std::endl
783  << "(Maximum number of entries for this row: " << arg2
784  << "; maybe the matrix is already compressed?)");
791  "The operation you attempted is only allowed after the SparsityPattern "
792  "has been set up and compress() was called.");
799  "The operation you attempted changes the structure of the SparsityPattern "
800  "and is not possible after compress() has been called.");
809  size_type,
810  size_type,
811  << "The iterators denote a range of " << arg1
812  << " elements, but the given number of rows was " << arg2);
821  size_type,
822  << "The number of partitions you gave is " << arg1
823  << ", but must be greater than zero.");
828  size_type,
829  size_type,
830  << "The array has size " << arg1 << " but should have size "
831  << arg2);
833 private:
838 
843 
848 
854 
858  template <typename>
859  friend class ChunkSparseMatrix;
860 
865 };
866 
867 
869 /*---------------------- Inline functions -----------------------------------*/
870 
871 #ifndef DOXYGEN
872 
874 {
875  inline Accessor::Accessor(const ChunkSparsityPattern *sparsity_pattern,
876  const size_type row)
877  : sparsity_pattern(sparsity_pattern)
878  , reduced_accessor(row == sparsity_pattern->n_rows() ?
879  *sparsity_pattern->sparsity_pattern.end() :
880  *sparsity_pattern->sparsity_pattern.begin(
881  row / sparsity_pattern->get_chunk_size()))
882  , chunk_row(row == sparsity_pattern->n_rows() ?
883  0 :
884  row % sparsity_pattern->get_chunk_size())
885  , chunk_col(0)
886  {}
887 
888 
889 
890  inline Accessor::Accessor(const ChunkSparsityPattern *sparsity_pattern)
891  : sparsity_pattern(sparsity_pattern)
892  , reduced_accessor(*sparsity_pattern->sparsity_pattern.end())
893  , chunk_row(0)
894  , chunk_col(0)
895  {}
896 
897 
898 
899  inline bool
900  Accessor::is_valid_entry() const
901  {
902  return reduced_accessor.is_valid_entry() &&
903  sparsity_pattern->get_chunk_size() * reduced_accessor.row() +
904  chunk_row <
905  sparsity_pattern->n_rows() &&
906  sparsity_pattern->get_chunk_size() * reduced_accessor.column() +
907  chunk_col <
908  sparsity_pattern->n_cols();
909  }
910 
911 
912 
913  inline Accessor::size_type
914  Accessor::row() const
915  {
916  Assert(is_valid_entry() == true, ExcInvalidIterator());
917 
918  return sparsity_pattern->get_chunk_size() * reduced_accessor.row() +
919  chunk_row;
920  }
921 
922 
923 
924  inline Accessor::size_type
925  Accessor::column() const
926  {
927  Assert(is_valid_entry() == true, ExcInvalidIterator());
928 
929  return sparsity_pattern->get_chunk_size() * reduced_accessor.column() +
930  chunk_col;
931  }
932 
933 
934 
935  inline std::size_t
936  Accessor::reduced_index() const
937  {
938  Assert(is_valid_entry() == true, ExcInvalidIterator());
939 
940  return reduced_accessor.linear_index;
941  }
942 
943 
944 
945  inline bool
946  Accessor::operator==(const Accessor &other) const
947  {
948  // no need to check for equality of sparsity patterns as this is done in
949  // the reduced case already and every ChunkSparsityPattern has its own
950  // reduced sparsity pattern
951  return (reduced_accessor == other.reduced_accessor &&
952  chunk_row == other.chunk_row && chunk_col == other.chunk_col);
953  }
954 
955 
956 
957  inline bool
958  Accessor::operator<(const Accessor &other) const
959  {
960  Assert(sparsity_pattern == other.sparsity_pattern, ExcInternalError());
961 
962  if (chunk_row != other.chunk_row)
963  {
964  if (reduced_accessor.linear_index ==
965  reduced_accessor.container->n_nonzero_elements())
966  return false;
967  if (other.reduced_accessor.linear_index ==
968  reduced_accessor.container->n_nonzero_elements())
969  return true;
970 
971  const auto global_row = sparsity_pattern->get_chunk_size() *
972  reduced_accessor.row() +
973  chunk_row,
974  other_global_row = sparsity_pattern->get_chunk_size() *
975  other.reduced_accessor.row() +
976  other.chunk_row;
977  if (global_row < other_global_row)
978  return true;
979  else if (global_row > other_global_row)
980  return false;
981  }
982 
983  return (
984  reduced_accessor.linear_index < other.reduced_accessor.linear_index ||
985  (reduced_accessor.linear_index == other.reduced_accessor.linear_index &&
986  chunk_col < other.chunk_col));
987  }
988 
989 
990  inline void
991  Accessor::advance()
992  {
993  const auto chunk_size = sparsity_pattern->get_chunk_size();
994  Assert(chunk_row < chunk_size && chunk_col < chunk_size,
996  Assert(reduced_accessor.row() * chunk_size + chunk_row <
997  sparsity_pattern->n_rows() &&
998  reduced_accessor.column() * chunk_size + chunk_col <
999  sparsity_pattern->n_cols(),
1000  ExcIteratorPastEnd());
1001  if (chunk_size == 1)
1002  {
1003  reduced_accessor.advance();
1004  return;
1005  }
1006 
1007  ++chunk_col;
1008 
1009  // end of chunk
1010  if (chunk_col == chunk_size ||
1011  reduced_accessor.column() * chunk_size + chunk_col ==
1012  sparsity_pattern->n_cols())
1013  {
1014  const auto reduced_row = reduced_accessor.row();
1015  // end of row
1016  if (reduced_accessor.linear_index + 1 ==
1017  reduced_accessor.container->rowstart[reduced_row + 1])
1018  {
1019  ++chunk_row;
1020 
1021  chunk_col = 0;
1022 
1023  // end of chunk rows or end of matrix
1024  if (chunk_row == chunk_size ||
1025  (reduced_row * chunk_size + chunk_row ==
1026  sparsity_pattern->n_rows()))
1027  {
1028  chunk_row = 0;
1029  reduced_accessor.advance();
1030  }
1031  // go back to the beginning of the same reduced row but with
1032  // chunk_row increased by one
1033  else
1034  reduced_accessor.linear_index =
1035  reduced_accessor.container->rowstart[reduced_row];
1036  }
1037  // advance within chunk
1038  else
1039  {
1040  reduced_accessor.advance();
1041  chunk_col = 0;
1042  }
1043  }
1044  }
1045 
1046 
1047 
1048  inline Iterator::Iterator(const ChunkSparsityPattern *sparsity_pattern,
1049  const size_type row)
1050  : accessor(sparsity_pattern, row)
1051  {}
1052 
1053 
1054 
1055  inline Iterator &
1056  Iterator::operator++()
1057  {
1058  accessor.advance();
1059  return *this;
1060  }
1061 
1062 
1063 
1064  inline Iterator
1065  Iterator::operator++(int)
1066  {
1067  const Iterator iter = *this;
1068  accessor.advance();
1069  return iter;
1070  }
1071 
1072 
1073 
1074  inline const Accessor &Iterator::operator*() const
1075  {
1076  return accessor;
1077  }
1078 
1079 
1080 
1081  inline const Accessor *Iterator::operator->() const
1082  {
1083  return &accessor;
1084  }
1085 
1086 
1087  inline bool
1088  Iterator::operator==(const Iterator &other) const
1089  {
1090  return (accessor == other.accessor);
1091  }
1092 
1093 
1094 
1095  inline bool
1096  Iterator::operator!=(const Iterator &other) const
1097  {
1098  return !(accessor == other.accessor);
1099  }
1100 
1101 
1102  inline bool
1103  Iterator::operator<(const Iterator &other) const
1104  {
1105  return accessor < other.accessor;
1106  }
1107 
1108 } // namespace ChunkSparsityPatternIterators
1109 
1110 
1111 
1114 {
1115  return {this, 0};
1116 }
1117 
1118 
1121 {
1122  return {this, n_rows()};
1123 }
1124 
1125 
1126 
1128 ChunkSparsityPattern::begin(const size_type r) const
1129 {
1130  Assert(r < n_rows(), ExcIndexRange(r, 0, n_rows()));
1131  return {this, r};
1132 }
1133 
1134 
1135 
1137 ChunkSparsityPattern::end(const size_type r) const
1138 {
1139  Assert(r < n_rows(), ExcIndexRange(r, 0, n_rows()));
1140  return {this, r + 1};
1141 }
1142 
1143 
1144 
1147 {
1148  return rows;
1149 }
1150 
1151 
1154 {
1155  return cols;
1156 }
1157 
1158 
1159 
1162 {
1163  return chunk_size;
1164 }
1165 
1166 
1167 
1168 inline bool
1170 {
1172 }
1173 
1174 
1175 
1176 template <typename ForwardIterator>
1177 void
1178 ChunkSparsityPattern::copy_from(const size_type n_rows,
1179  const size_type n_cols,
1180  const ForwardIterator begin,
1181  const ForwardIterator end,
1182  const size_type chunk_size)
1183 {
1184  Assert(static_cast<size_type>(std::distance(begin, end)) == n_rows,
1185  ExcIteratorRange(std::distance(begin, end), n_rows));
1186 
1187  // first determine row lengths for each row. if the matrix is quadratic,
1188  // then we might have to add an additional entry for the diagonal, if that
1189  // is not yet present. as we have to call compress anyway later on, don't
1190  // bother to check whether that diagonal entry is in a certain row or not
1191  const bool is_square = (n_rows == n_cols);
1192  std::vector<size_type> row_lengths;
1193  row_lengths.reserve(n_rows);
1194  for (ForwardIterator i = begin; i != end; ++i)
1195  row_lengths.push_back(std::distance(i->begin(), i->end()) +
1196  (is_square ? 1 : 0));
1197  reinit(n_rows, n_cols, row_lengths, chunk_size);
1198 
1199  // now enter all the elements into the matrix
1200  size_type row = 0;
1201  using inner_iterator =
1202  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1203  for (ForwardIterator i = begin; i != end; ++i, ++row)
1204  {
1205  const inner_iterator end_of_row = i->end();
1206  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1207  {
1208  const size_type col =
1209  internal::SparsityPatternTools::get_column_index_from_iterator(*j);
1210  Assert(col < n_cols, ExcInvalidIndex(col, n_cols));
1211 
1212  add(row, col);
1213  }
1214  }
1215 
1216  // finally compress everything. this also sorts the entries within each row
1217  compress();
1218 }
1219 
1220 
1221 #endif // DOXYGEN
1222 
1223 DEAL_II_NAMESPACE_CLOSE
1224 
1225 #endif
bool operator==(const Iterator &) const
size_type get_chunk_size() const
void block_write(std::ostream &out) 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)
std::size_t memory_consumption() const
static ::ExceptionBase & ExcEmptyObject()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
static ::ExceptionBase & ExcNotEnoughSpace(size_type arg1, size_type arg2)
types::global_dof_index size_type
iterator end() const
void block_read(std::istream &in)
bool operator!=(const Iterator &) const
Iterator(const ChunkSparsityPattern *sp, const size_type row)
bool exists(const size_type i, const size_type j) 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)
const ChunkSparsityPattern * sparsity_pattern
void add(const size_type i, const size_type j)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
static ::ExceptionBase & ExcInvalidNumber(size_type arg1)
SparsityPattern sparsity_pattern
bool operator<(const Iterator &) const
static ::ExceptionBase & ExcInvalidIterator()
size_type n_cols() const
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
Accessor(const ChunkSparsityPattern *matrix, const size_type row)
bool is_compressed() const
#define Assert(cond, exc)
Definition: exceptions.h:1407
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:496
#define DeclException0(Exception0)
Definition: exceptions.h:473
size_type row_length(const size_type row) const
static ::ExceptionBase & ExcInvalidNumberOfPartitions(size_type arg1)
size_type n_rows() const
size_type n_nonzero_elements() const
const Accessor & operator*() const
static ::ExceptionBase & ExcInvalidArraySize(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
ChunkSparsityPattern & operator=(const ChunkSparsityPattern &)
static ::ExceptionBase & ExcIteratorPastEnd()
unsigned int global_dof_index
Definition: types.h:89
static ::ExceptionBase & ExcMatrixIsCompressed()
MatrixTableIterators::Iterator< TransposeTable< T >, Constness, MatrixTableIterators::Storage::column_major > Iterator
Definition: table.h:1912
static ::ExceptionBase & ExcIteratorRange(size_type arg1, size_type arg2)
constexpr int chunk_size
Definition: cuda_size.h:35
static ::ExceptionBase & ExcMETISNotInstalled()
iterator begin() const
~ChunkSparsityPattern() override=default
static ::ExceptionBase & ExcNotCompressed()
static const size_type invalid_entry
static const size_type invalid_entry
bool operator<(const Accessor &) const
void print(std::ostream &out) const
SparsityPatternIterators::Accessor reduced_accessor
void print_gnuplot(std::ostream &out) const
void reinit(const size_type m, const size_type n, const size_type max_per_row, const size_type chunk_size)
bool operator==(const Accessor &) const
static ::ExceptionBase & ExcInternalError()
const Accessor * operator->() const
size_type max_entries_per_row() const