Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\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 // Forward declaration
34 #ifndef DOXYGEN
35 template <typename>
36 class ChunkSparseMatrix;
37 #endif
38 
49 {
50  // forward declaration
51  class Iterator;
52 
65  class Accessor
66  {
67  public:
72 
76  Accessor(const ChunkSparsityPattern *matrix, const size_type row);
77 
81  Accessor(const ChunkSparsityPattern *matrix);
82 
87  size_type
88  row() const;
89 
93  std::size_t
94  reduced_index() const;
95 
100  size_type
101  column() const;
102 
112  bool
113  is_valid_entry() const;
114 
115 
119  bool
120  operator==(const Accessor &) const;
121 
122 
130  bool
131  operator<(const Accessor &) const;
132 
133  protected:
138 
143 
148 
153 
157  void
158  advance();
159 
160  // Grant access to iterator class.
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 
283  static const size_type invalid_entry = SparsityPattern::invalid_entry;
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 
366  operator=(const ChunkSparsityPattern &);
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
659  stores_only_added_elements() const;
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);
770  DeclException2(ExcInvalidIndex,
771  size_type,
772  size_type,
773  << "The given index " << arg1 << " should be less than "
774  << arg2 << ".");
778  DeclException2(ExcNotEnoughSpace,
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?)");
790  ExcNotCompressed,
791  "The operation you attempted is only allowed after the SparsityPattern "
792  "has been set up and compress() was called.");
798  ExcMatrixIsCompressed,
799  "The operation you attempted changes the structure of the SparsityPattern "
800  "and is not possible after compress() has been called.");
808  DeclException2(ExcIteratorRange,
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 
855  // Make all the chunk sparse matrix kinds friends.
856  template <typename>
857  friend class ChunkSparseMatrix;
858 
859  // Make the accessor class a friend.
861 };
862 
863 
865 /*---------------------- Inline functions -----------------------------------*/
866 
867 #ifndef DOXYGEN
868 
870 {
872  const size_type row)
873  : sparsity_pattern(sparsity_pattern)
874  , reduced_accessor(row == sparsity_pattern->n_rows() ?
875  *sparsity_pattern->sparsity_pattern.end() :
876  *sparsity_pattern->sparsity_pattern.begin(
877  row / sparsity_pattern->get_chunk_size()))
878  , chunk_row(row == sparsity_pattern->n_rows() ?
879  0 :
880  row % sparsity_pattern->get_chunk_size())
881  , chunk_col(0)
882  {}
883 
884 
885 
886  inline Accessor::Accessor(const ChunkSparsityPattern *sparsity_pattern)
887  : sparsity_pattern(sparsity_pattern)
888  , reduced_accessor(*sparsity_pattern->sparsity_pattern.end())
889  , chunk_row(0)
890  , chunk_col(0)
891  {}
892 
893 
894 
895  inline bool
896  Accessor::is_valid_entry() const
897  {
898  return reduced_accessor.is_valid_entry() &&
899  sparsity_pattern->get_chunk_size() * reduced_accessor.row() +
900  chunk_row <
902  sparsity_pattern->get_chunk_size() * reduced_accessor.column() +
903  chunk_col <
905  }
906 
907 
908 
909  inline Accessor::size_type
910  Accessor::row() const
911  {
912  Assert(is_valid_entry() == true, ExcInvalidIterator());
913 
914  return sparsity_pattern->get_chunk_size() * reduced_accessor.row() +
915  chunk_row;
916  }
917 
918 
919 
920  inline Accessor::size_type
921  Accessor::column() const
922  {
923  Assert(is_valid_entry() == true, ExcInvalidIterator());
924 
925  return sparsity_pattern->get_chunk_size() * reduced_accessor.column() +
926  chunk_col;
927  }
928 
929 
930 
931  inline std::size_t
932  Accessor::reduced_index() const
933  {
934  Assert(is_valid_entry() == true, ExcInvalidIterator());
935 
936  return reduced_accessor.linear_index;
937  }
938 
939 
940 
941  inline bool
942  Accessor::operator==(const Accessor &other) const
943  {
944  // no need to check for equality of sparsity patterns as this is done in
945  // the reduced case already and every ChunkSparsityPattern has its own
946  // reduced sparsity pattern
947  return (reduced_accessor == other.reduced_accessor &&
948  chunk_row == other.chunk_row && chunk_col == other.chunk_col);
949  }
950 
951 
952 
953  inline bool
954  Accessor::operator<(const Accessor &other) const
955  {
956  Assert(sparsity_pattern == other.sparsity_pattern, ExcInternalError());
957 
958  if (chunk_row != other.chunk_row)
959  {
960  if (reduced_accessor.linear_index ==
961  reduced_accessor.container->n_nonzero_elements())
962  return false;
963  if (other.reduced_accessor.linear_index ==
964  reduced_accessor.container->n_nonzero_elements())
965  return true;
966 
967  const auto global_row = sparsity_pattern->get_chunk_size() *
968  reduced_accessor.row() +
969  chunk_row,
970  other_global_row = sparsity_pattern->get_chunk_size() *
971  other.reduced_accessor.row() +
972  other.chunk_row;
973  if (global_row < other_global_row)
974  return true;
975  else if (global_row > other_global_row)
976  return false;
977  }
978 
979  return (
980  reduced_accessor.linear_index < other.reduced_accessor.linear_index ||
981  (reduced_accessor.linear_index == other.reduced_accessor.linear_index &&
982  chunk_col < other.chunk_col));
983  }
984 
985 
986  inline void
987  Accessor::advance()
988  {
989  const auto chunk_size = sparsity_pattern->get_chunk_size();
990  Assert(chunk_row < chunk_size && chunk_col < chunk_size,
992  Assert(reduced_accessor.row() * chunk_size + chunk_row <
994  reduced_accessor.column() * chunk_size + chunk_col <
997  if (chunk_size == 1)
998  {
999  reduced_accessor.advance();
1000  return;
1001  }
1002 
1003  ++chunk_col;
1004 
1005  // end of chunk
1006  if (chunk_col == chunk_size ||
1007  reduced_accessor.column() * chunk_size + chunk_col ==
1009  {
1010  const auto reduced_row = reduced_accessor.row();
1011  // end of row
1012  if (reduced_accessor.linear_index + 1 ==
1013  reduced_accessor.container->rowstart[reduced_row + 1])
1014  {
1015  ++chunk_row;
1016 
1017  chunk_col = 0;
1018 
1019  // end of chunk rows or end of matrix
1020  if (chunk_row == chunk_size ||
1021  (reduced_row * chunk_size + chunk_row ==
1023  {
1024  chunk_row = 0;
1025  reduced_accessor.advance();
1026  }
1027  // go back to the beginning of the same reduced row but with
1028  // chunk_row increased by one
1029  else
1030  reduced_accessor.linear_index =
1031  reduced_accessor.container->rowstart[reduced_row];
1032  }
1033  // advance within chunk
1034  else
1035  {
1036  reduced_accessor.advance();
1037  chunk_col = 0;
1038  }
1039  }
1040  }
1041 
1042 
1043 
1044  inline Iterator::Iterator(const ChunkSparsityPattern *sparsity_pattern,
1045  const size_type row)
1046  : accessor(sparsity_pattern, row)
1047  {}
1048 
1049 
1050 
1051  inline Iterator &
1052  Iterator::operator++()
1053  {
1054  accessor.advance();
1055  return *this;
1056  }
1057 
1058 
1059 
1060  inline Iterator
1061  Iterator::operator++(int)
1062  {
1063  const Iterator iter = *this;
1064  accessor.advance();
1065  return iter;
1066  }
1067 
1068 
1069 
1070  inline const Accessor &Iterator::operator*() const
1071  {
1072  return accessor;
1073  }
1074 
1075 
1076 
1077  inline const Accessor *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  Assert(r < n_rows(), ExcIndexRange(r, 0, n_rows()));
1127  return {this, r};
1128 }
1129 
1130 
1131 
1133 ChunkSparsityPattern::end(const size_type r) const
1134 {
1135  Assert(r < n_rows(), ExcIndexRange(r, 0, 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 =
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 
1219 DEAL_II_NAMESPACE_CLOSE
1220 
1221 #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)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
types::global_dof_index size_type
static ::ExceptionBase & ExcMETISNotInstalled()
iterator end() const
const ChunkSparsityPattern * sparsity_pattern
void add(const size_type i, const size_type j)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
SparsityPattern sparsity_pattern
static ::ExceptionBase & ExcInvalidIterator()
size_type n_cols() const
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
#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
static ::ExceptionBase & ExcInvalidArraySize(int arg1, int arg2)
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:496
#define DeclException0(Exception0)
Definition: exceptions.h:473
static ::ExceptionBase & ExcInvalidNumber(unsigned int arg1)
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:89
static ::ExceptionBase & ExcIteratorRange(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcEmptyObject()
iterator begin() const
static const size_type invalid_entry
bool operator<(const Accessor &) const
DEAL_II_CONSTEXPR SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
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
static ::ExceptionBase & ExcInternalError()
size_type n_cols() const