Reference documentation for deal.II version 8.1.0
chunk_sparsity_pattern.h
1 // ---------------------------------------------------------------------
2 // @f$Id: chunk_sparsity_pattern.h 31414 2013-10-24 22:15:50Z bangerth @f$
3 //
4 // Copyright (C) 2008 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__chunk_sparsity_pattern_h
18 #define __deal2__chunk_sparsity_pattern_h
19 
20 
21 #include <deal.II/base/config.h>
23 #include <deal.II/base/subscriptor.h>
24 #include <deal.II/base/vector_slice.h>
25 
26 #include <deal.II/lac/sparsity_pattern.h>
27 
28 #include <vector>
29 #include <iostream>
30 
32 
33 
34 template <typename> class ChunkSparseMatrix;
35 
36 
47 {
48  // forward declaration
49  class Iterator;
50 
63  class Accessor
64  {
65  public:
69  Accessor (const ChunkSparsityPattern *matrix,
70  const unsigned int row);
71 
75  Accessor (const ChunkSparsityPattern *matrix);
76 
81  unsigned int row () const;
82 
86  std::size_t reduced_index() const;
87 
92  unsigned int column () const;
93 
103  bool is_valid_entry () const;
104 
105 
109  bool operator == (const Accessor &) const;
110 
111 
119  bool operator < (const Accessor &) const;
120 
121  protected:
126 
131 
135  unsigned int chunk_row;
136 
140  unsigned int chunk_col;
141 
145  void advance ();
146 
150  friend class Iterator;
151  };
152 
153 
154 
158  class Iterator
159  {
160  public:
165  Iterator (const ChunkSparsityPattern *sp,
166  const unsigned int row);
167 
172 
176  Iterator operator++ (int);
177 
181  const Accessor &operator* () const;
182 
186  const Accessor *operator-> () const;
187 
191  bool operator == (const Iterator &) const;
192 
196  bool operator != (const Iterator &) const;
197 
205  bool operator < (const Iterator &) const;
206 
207  private:
212  };
213 }
214 
215 
216 
227 {
228 public:
238 
247 
263 
270 
287 
295  ChunkSparsityPattern (const size_type m,
296  const size_type n,
297  const size_type max_chunks_per_row,
298  const size_type chunk_size);
299 
304  ChunkSparsityPattern (const size_type m,
305  const size_type n,
306  const size_type max_chunks_per_row,
307  const size_type chunk_size,
309 
318  ChunkSparsityPattern (const size_type m,
319  const size_type n,
320  const std::vector<size_type> &row_lengths,
321  const size_type chunk_size);
322 
327  ChunkSparsityPattern (const size_type m,
328  const size_type n,
329  const std::vector<size_type> &row_lengths,
330  const size_type chunk_size,
332 
341  ChunkSparsityPattern (const size_type n,
342  const size_type max_per_row,
343  const size_type chunk_size);
344 
352  ChunkSparsityPattern (const size_type m,
353  const std::vector<size_type> &row_lengths,
354  const size_type chunk_size);
355 
360  ChunkSparsityPattern (const size_type m,
361  const std::vector<size_type> &row_lengths,
362  const size_type chunk_size,
364 
369 
376 
385  void reinit (const size_type m,
386  const size_type n,
387  const size_type max_per_row,
388  const size_type chunk_size);
389 
394  void reinit (const size_type m,
395  const size_type n,
396  const size_type max_per_row,
397  const size_type chunk_size,
399 
417  void reinit (const size_type m,
418  const size_type n,
419  const std::vector<size_type> &row_lengths,
420  const size_type chunk_size);
421 
426  void reinit (const size_type m,
427  const size_type n,
428  const std::vector<size_type > &row_lengths,
429  const size_type chunk_size,
431 
435  void reinit (const size_type m,
436  const size_type n,
437  const VectorSlice<const std::vector<size_type> > &row_lengths,
438  const size_type chunk_size);
439 
444  void reinit (const size_type m,
445  const size_type n,
446  const VectorSlice<const std::vector<size_type> > &row_lengths,
447  const size_type chunk_size,
449 
462  void compress ();
463 
540  template <typename ForwardIterator>
541  void copy_from (const size_type n_rows,
542  const size_type n_cols,
543  const ForwardIterator begin,
544  const ForwardIterator end,
545  const size_type chunk_size);
546 
551  template <typename ForwardIterator>
552  void copy_from (const size_type n_rows,
553  const size_type n_cols,
554  const ForwardIterator begin,
555  const ForwardIterator end,
556  const size_type chunk_size,
558 
565  template <typename SparsityType>
566  void copy_from (const SparsityType &csp,
567  const size_type chunk_size);
568 
573  template <typename SparsityType>
574  void copy_from (const SparsityType &csp,
575  const size_type chunk_size,
577 
585  template <typename number>
586  void copy_from (const FullMatrix<number> &matrix,
587  const size_type chunk_size);
588 
593  template <typename number>
594  void copy_from (const FullMatrix<number> &matrix,
595  const size_type chunk_size,
597 
615  template <typename Sparsity>
616  void create_from (const unsigned int m,
617  const unsigned int n,
618  const Sparsity &sparsity_pattern_for_chunks,
619  const unsigned int chunk_size,
620  const bool optimize_diagonal = true);
621 
626  bool empty () const;
627 
631  size_type get_chunk_size () const;
632 
638  size_type max_entries_per_row () const;
639 
646  void add (const size_type i,
647  const size_type j);
648 
656  void symmetrize ();
657 
662  inline size_type n_rows () const;
663 
668  inline size_type n_cols () const;
669 
673  bool exists (const size_type i,
674  const size_type j) const;
675 
679  size_type row_length (const size_type row) const;
680 
687  size_type bandwidth () const;
688 
697  size_type n_nonzero_elements () const;
698 
702  bool is_compressed () const;
703 
711  bool optimize_diagonal () const DEAL_II_DEPRECATED;
712 
724  bool stores_only_added_elements () const;
725 
731  iterator begin () const;
732 
736  iterator end () const;
737 
746  iterator begin (const unsigned int r) const;
747 
756  iterator end (const unsigned int r) const;
757 
768  void block_write (std::ostream &out) const;
769 
783  void block_read (std::istream &in);
784 
790  void print (std::ostream &out) const;
791 
805  void print_gnuplot (std::ostream &out) const;
806 
811  std::size_t memory_consumption () const;
812 
818  DeclException1 (ExcInvalidNumber,
819  int,
820  << "The provided number is invalid here: " << arg1);
824  DeclException2 (ExcInvalidIndex,
825  int, int,
826  << "The given index " << arg1
827  << " should be less than " << arg2 << ".");
831  DeclException2 (ExcNotEnoughSpace,
832  int, int,
833  << "Upon entering a new entry to row " << arg1
834  << ": there was no free entry any more. " << std::endl
835  << "(Maximum number of entries for this row: "
836  << arg2 << "; maybe the matrix is already compressed?)");
840  DeclException0 (ExcNotCompressed);
844  DeclException0 (ExcMatrixIsCompressed);
856  DeclException2 (ExcIteratorRange,
857  int, int,
858  << "The iterators denote a range of " << arg1
859  << " elements, but the given number of rows was " << arg2);
863  DeclException0 (ExcMETISNotInstalled);
867  DeclException1 (ExcInvalidNumberOfPartitions,
868  int,
869  << "The number of partitions you gave is " << arg1
870  << ", but must be greater than zero.");
874  DeclException2 (ExcInvalidArraySize,
875  int, int,
876  << "The array has size " << arg1 << " but should have size "
877  << arg2);
879 private:
883  size_type rows;
884 
888  size_type cols;
889 
893  size_type chunk_size;
894 
900 
904  template <typename> friend class ChunkSparseMatrix;
905 
909  friend class ChunkSparsityPatternIterators::Accessor;
910 };
911 
912 
914 /*---------------------- Inline functions -----------------------------------*/
915 
916 #ifndef DOXYGEN
917 
919 {
920  inline
921  Accessor::
922  Accessor (const ChunkSparsityPattern *sparsity_pattern,
923  const unsigned int row)
924  :
925  sparsity_pattern(sparsity_pattern),
926  reduced_accessor(row==sparsity_pattern->n_rows() ?
927  *sparsity_pattern->sparsity_pattern.end() :
928  *sparsity_pattern->sparsity_pattern.
929  begin(row/sparsity_pattern->get_chunk_size())),
930  chunk_row (row==sparsity_pattern->n_rows() ? 0 :
931  row%sparsity_pattern->get_chunk_size()),
932  chunk_col (0)
933  {}
934 
935 
936 
937  inline
938  Accessor::
939  Accessor (const ChunkSparsityPattern *sparsity_pattern)
940  :
941  sparsity_pattern(sparsity_pattern),
942  reduced_accessor(*sparsity_pattern->sparsity_pattern.end()),
943  chunk_row (0),
944  chunk_col (0)
945  {}
946 
947 
948 
949  inline
950  bool
951  Accessor::is_valid_entry () const
952  {
953  return reduced_accessor.is_valid_entry()
954  &&
955  sparsity_pattern->get_chunk_size()*reduced_accessor.row()+chunk_row <
956  sparsity_pattern->n_rows()
957  &&
958  sparsity_pattern->get_chunk_size()*reduced_accessor.column()+chunk_col <
959  sparsity_pattern->n_cols();
960  }
961 
962 
963 
964  inline
965  unsigned int
966  Accessor::row() const
967  {
968  Assert (is_valid_entry() == true, ExcInvalidIterator());
969 
970  return sparsity_pattern->get_chunk_size()*reduced_accessor.row() +
971  chunk_row;
972  }
973 
974 
975 
976  inline
977  unsigned int
978  Accessor::column() const
979  {
980  Assert (is_valid_entry() == true, ExcInvalidIterator());
981 
982  return sparsity_pattern->get_chunk_size()*reduced_accessor.column() +
983  chunk_col;
984  }
985 
986 
987 
988  inline
989  std::size_t
990  Accessor::reduced_index() const
991  {
992  Assert (is_valid_entry() == true, ExcInvalidIterator());
993 
994  return reduced_accessor.index_within_sparsity;
995  }
996 
997 
998 
999 
1000  inline
1001  bool
1002  Accessor::operator == (const Accessor &other) const
1003  {
1004  // no need to check for equality of sparsity patterns as this is done in
1005  // the reduced case already and every ChunkSparsityPattern has its own
1006  // reduced sparsity pattern
1007  return (reduced_accessor == other.reduced_accessor &&
1008  chunk_row == other.chunk_row &&
1009  chunk_col == other.chunk_col);
1010  }
1011 
1012 
1013 
1014  inline
1015  bool
1016  Accessor::operator < (const Accessor &other) const
1017  {
1018  Assert (sparsity_pattern == other.sparsity_pattern,
1019  ExcInternalError());
1020 
1021  if (chunk_row != other.chunk_row)
1022  {
1023  if (reduced_accessor.index_within_sparsity ==
1024  reduced_accessor.sparsity_pattern->n_nonzero_elements())
1025  return false;
1026  if (other.reduced_accessor.index_within_sparsity ==
1027  reduced_accessor.sparsity_pattern->n_nonzero_elements())
1028  return true;
1029 
1030  const unsigned int
1031  global_row = sparsity_pattern->get_chunk_size()*reduced_accessor.row()
1032  +chunk_row,
1033  other_global_row = sparsity_pattern->get_chunk_size()*
1034  other.reduced_accessor.row()+other.chunk_row;
1035  if (global_row < other_global_row)
1036  return true;
1037  else if (global_row > other_global_row)
1038  return false;
1039  }
1040 
1041  return (reduced_accessor.index_within_sparsity <
1042  other.reduced_accessor.index_within_sparsity ||
1043  (reduced_accessor.index_within_sparsity ==
1044  other.reduced_accessor.index_within_sparsity &&
1045  chunk_col < other.chunk_col));
1046  }
1047 
1048 
1049  inline
1050  void
1051  Accessor::advance ()
1052  {
1053  const unsigned int chunk_size = sparsity_pattern->get_chunk_size();
1054  Assert (chunk_row < chunk_size && chunk_col < chunk_size,
1055  ExcIteratorPastEnd());
1056  Assert (reduced_accessor.row() * chunk_size + chunk_row <
1057  sparsity_pattern->n_rows()
1058  &&
1059  reduced_accessor.column() * chunk_size + chunk_col <
1060  sparsity_pattern->n_cols(),
1061  ExcIteratorPastEnd());
1062  if (chunk_size == 1)
1063  {
1064  reduced_accessor.advance();
1065  return;
1066  }
1067 
1068  ++chunk_col;
1069 
1070  // end of chunk
1071  if (chunk_col == chunk_size
1072  ||
1073  reduced_accessor.column() * chunk_size + chunk_col ==
1074  sparsity_pattern->n_cols())
1075  {
1076  const unsigned int reduced_row = reduced_accessor.row();
1077  // end of row
1078  if (reduced_accessor.index_within_sparsity + 1 ==
1079  reduced_accessor.sparsity_pattern->rowstart[reduced_row+1])
1080  {
1081  ++chunk_row;
1082 
1083  chunk_col = 0;
1084 
1085  // end of chunk rows or end of matrix
1086  if (chunk_row == chunk_size ||
1087  (reduced_row * chunk_size + chunk_row ==
1088  sparsity_pattern->n_rows()))
1089  {
1090  chunk_row = 0;
1091  reduced_accessor.advance();
1092  }
1093  // go back to the beginning of the same reduced row but with
1094  // chunk_row increased by one
1095  else
1096  reduced_accessor.index_within_sparsity =
1097  reduced_accessor.sparsity_pattern->rowstart[reduced_row];
1098  }
1099  // advance within chunk
1100  else
1101  {
1102  reduced_accessor.advance();
1103  chunk_col = 0;
1104  }
1105  }
1106  }
1107 
1108 
1109 
1110  inline
1111  Iterator::Iterator (const ChunkSparsityPattern *sparsity_pattern,
1112  const unsigned int row)
1113  :
1114  accessor(sparsity_pattern, row)
1115  {}
1116 
1117 
1118 
1119  inline
1120  Iterator &
1121  Iterator::operator++ ()
1122  {
1123  accessor.advance ();
1124  return *this;
1125  }
1126 
1127 
1128 
1129  inline
1130  Iterator
1131  Iterator::operator++ (int)
1132  {
1133  const Iterator iter = *this;
1134  accessor.advance ();
1135  return iter;
1136  }
1137 
1138 
1139 
1140  inline
1141  const Accessor &
1142  Iterator::operator* () const
1143  {
1144  return accessor;
1145  }
1146 
1147 
1148 
1149  inline
1150  const Accessor *
1151  Iterator::operator-> () const
1152  {
1153  return &accessor;
1154  }
1155 
1156 
1157  inline
1158  bool
1159  Iterator::operator == (const Iterator &other) const
1160  {
1161  return (accessor == other.accessor);
1162  }
1163 
1164 
1165 
1166  inline
1167  bool
1168  Iterator::operator != (const Iterator &other) const
1169  {
1170  return ! (accessor == other.accessor);
1171  }
1172 
1173 
1174  inline
1175  bool
1176  Iterator::operator < (const Iterator &other) const
1177  {
1178  return accessor < other.accessor;
1179  }
1180 
1181 }
1182 
1183 
1184 
1185 inline
1188 {
1189  return iterator(this, 0);
1190 }
1191 
1192 
1193 inline
1196 {
1197  return iterator(this, n_rows());
1198 }
1199 
1200 
1201 
1202 inline
1204 ChunkSparsityPattern::begin (const unsigned int r) const
1205 {
1206  Assert (r<n_rows(), ExcIndexRange(r,0,n_rows()));
1207  return iterator(this, r);
1208 }
1209 
1210 
1211 
1212 inline
1214 ChunkSparsityPattern::end (const unsigned int r) const
1215 {
1216  Assert (r<n_rows(), ExcIndexRange(r,0,n_rows()))
1217  return iterator(this, r+1);
1218 }
1219 
1220 
1221 
1222 inline
1223 ChunkSparsityPattern::size_type
1224 ChunkSparsityPattern::n_rows () const
1225 {
1226  return rows;
1227 }
1228 
1229 
1230 inline
1233 {
1234  return cols;
1235 }
1236 
1237 
1238 
1239 inline
1242 {
1243  return chunk_size;
1244 }
1245 
1246 
1247 
1248 inline
1249 bool
1251 {
1252  return sparsity_pattern.compressed;
1253 }
1254 
1255 
1256 
1257 inline
1258 bool
1260 {
1261  return sparsity_pattern.store_diagonal_first_in_row;
1262 }
1263 
1264 
1265 
1266 template <typename ForwardIterator>
1267 inline
1268 void
1269 ChunkSparsityPattern::copy_from (const size_type n_rows,
1270  const size_type n_cols,
1271  const ForwardIterator begin,
1272  const ForwardIterator end,
1273  const size_type chunk_size,
1274  const bool)
1275 {
1276  copy_from (n_rows, n_cols, begin, end, chunk_size);
1277 }
1278 
1279 
1280 
1281 template <typename ForwardIterator>
1282 void
1283 ChunkSparsityPattern::copy_from (const size_type n_rows,
1284  const size_type n_cols,
1285  const ForwardIterator begin,
1286  const ForwardIterator end,
1287  const size_type chunk_size)
1288 {
1289  Assert (static_cast<size_type>(std::distance (begin, end)) == n_rows,
1290  ExcIteratorRange (std::distance (begin, end), n_rows));
1291 
1292  // first determine row lengths for each row. if the matrix is quadratic,
1293  // then we might have to add an additional entry for the diagonal, if that
1294  // is not yet present. as we have to call compress anyway later on, don't
1295  // bother to check whether that diagonal entry is in a certain row or not
1296  const bool is_square = (n_rows == n_cols);
1297  std::vector<size_type> row_lengths;
1298  row_lengths.reserve(n_rows);
1299  for (ForwardIterator i=begin; i!=end; ++i)
1300  row_lengths.push_back (std::distance (i->begin(), i->end())
1301  +
1302  (is_square ? 1 : 0));
1303  reinit (n_rows, n_cols, row_lengths, chunk_size);
1304 
1305  // now enter all the elements into the matrix
1306  size_type row = 0;
1307  typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
1308  for (ForwardIterator i=begin; i!=end; ++i, ++row)
1309  {
1310  const inner_iterator end_of_row = i->end();
1311  for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
1312  {
1313  const size_type col
1314  = internal::SparsityPatternTools::get_column_index_from_iterator(*j);
1315  Assert (col < n_cols, ExcInvalidIndex(col,n_cols));
1316 
1317  add (row, col);
1318  }
1319  }
1320 
1321  // finally compress everything. this also sorts the entries within each row
1322  compress ();
1323 }
1324 
1325 
1326 #endif // DOXYGEN
1327 
1328 DEAL_II_NAMESPACE_CLOSE
1329 
1330 #endif
size_type n_rows() const
bool optimize_diagonal() const DEAL_II_DEPRECATED
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 block_write(std::ostream &out) const
void create_from(const unsigned int m, const unsigned int n, const Sparsity &sparsity_pattern_for_chunks, const unsigned int chunk_size, const bool optimize_diagonal=true)
std::size_t memory_consumption() const
Iterator(const ChunkSparsityPattern *sp, const unsigned int row)
void print_gnuplot(std::ostream &out) const
ChunkSparsityPattern & operator=(const ChunkSparsityPattern &)
bool is_compressed() const
bool exists(const size_type i, const size_type j) const
ChunkSparsityPatternIterators::Iterator const_iterator
void block_read(std::istream &in)
bool operator==(const Accessor &) const
STL namespace.
const Accessor * operator->() const
const ChunkSparsityPattern * sparsity_pattern
void add(const size_type i, const size_type j)
size_type n_nonzero_elements() const
SparsityPattern sparsity_pattern
iterator end() const
DeclException2(ExcInvalidIndex, int, int,<< "The given index "<< arg1<< " should be less than "<< arg2<< ".")
iterator begin() const
bool stores_only_added_elements() const
size_type get_chunk_size() const
Iterator is invalid, probably due to an error.
unsigned int global_dof_index
Definition: types.h:100
size_type row_length(const size_type row) const
ChunkSparsityPatternIterators::Iterator iterator
#define Assert(cond, exc)
Definition: exceptions.h:299
::ExceptionBase & ExcEmptyObject()
bool operator<(const Accessor &) const
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
::ExceptionBase & ExcInvalidConstructorCall()
void print(std::ostream &out) const
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
iterator end() const
bool operator<(const Iterator &) const
types::global_dof_index size_type
Accessor(const ChunkSparsityPattern *matrix, const unsigned int row)
::ExceptionBase & ExcIteratorPastEnd()
size_type n_cols() const
bool operator==(const Iterator &) const
DeclException0(ExcNotCompressed)
size_type bandwidth() const
static const size_type invalid_entry
size_type max_entries_per_row() const
::ExceptionBase & ExcInvalidIterator()
bool operator!=(const Iterator &) const
::ExceptionBase & ExcInternalError()
DeclException1(ExcInvalidNumber, int,<< "The provided number is invalid here: "<< arg1)
SparsityPatternIterators::Accessor reduced_accessor
bool empty() const
void reinit(const size_type m, const size_type n, const size_type max_per_row, const size_type chunk_size)
const Accessor & operator*() const
static const size_type invalid_entry