Reference documentation for deal.II version 8.4.1
sparsity_pattern.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2016 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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii__sparsity_pattern_h
17 #define dealii__sparsity_pattern_h
18 
19 
20 #include <deal.II/base/config.h>
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/subscriptor.h>
23 #include <boost/serialization/array.hpp>
24 #include <boost/serialization/split_member.hpp>
25 
26 #include <vector>
27 #include <iostream>
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 class SparsityPattern;
33 template <typename number> class FullMatrix;
34 template <typename number> class SparseMatrix;
35 template <typename number> class SparseLUDecomposition;
36 template <typename number> class SparseILU;
37 template <typename VectorType> class VectorSlice;
38 
40 {
41  class Accessor;
42 }
43 
44 
49 namespace internals
50 {
51  namespace SparsityPatternTools
52  {
57 
63  size_type
64  get_column_index_from_iterator (const size_type i);
65 
71  template <typename value>
72  size_type
73  get_column_index_from_iterator (const std::pair<size_type, value> &i);
74 
80  template <typename value>
81  size_type
82  get_column_index_from_iterator (const std::pair<const size_type, value> &i);
83 
84  }
85 }
86 
87 
92 {
93  // forward declaration
94  class Iterator;
95 
100 
114  class Accessor
115  {
116  public:
120  Accessor (const SparsityPattern *matrix,
121  const std::size_t index_within_sparsity);
122 
126  Accessor (const SparsityPattern *matrix);
127 
132  size_type row () const;
133 
139  size_type index () const;
140 
145  size_type column () const;
146 
156  bool is_valid_entry () const;
157 
161  bool operator == (const Accessor &) const;
162 
170  bool operator < (const Accessor &) const;
171 
172  protected:
177 
186 
190  void advance ();
191 
195  friend class Iterator;
196 
201  };
202 
203 
204 
229  class Iterator
230  {
231  public:
237  Iterator (const SparsityPattern *sp,
238  const std::size_t index_within_sparsity);
239 
244 
248  Iterator operator++ (int);
249 
253  const Accessor &operator* () const;
254 
258  const Accessor *operator-> () const;
259 
263  bool operator == (const Iterator &) const;
264 
268  bool operator != (const Iterator &) const;
269 
277  bool operator < (const Iterator &) const;
278 
285  int operator - (const Iterator &p) const;
286 
287  private:
292  };
293 }
294 
295 
296 
332 {
333 public:
338 
343  typedef
346 
354  typedef
357 
358 
373  static const size_type invalid_entry = numbers::invalid_size_type;
374 
379 // @{
385  SparsityPattern ();
386 
403 
411  SparsityPattern (const size_type m,
412  const size_type n,
413  const unsigned int max_per_row);
414 
415 
424  SparsityPattern (const size_type m,
425  const size_type n,
426  const std::vector<unsigned int> &row_lengths);
427 
436  SparsityPattern (const size_type m,
437  const unsigned int max_per_row);
438 
446  SparsityPattern (const size_type m,
447  const std::vector<unsigned int> &row_lengths);
448 
471  SparsityPattern (const SparsityPattern &original,
472  const unsigned int max_per_row,
473  const size_type extra_off_diagonals);
474 
478  ~SparsityPattern ();
479 
486 
495  void reinit (const size_type m,
496  const size_type n,
497  const unsigned int max_per_row);
498 
499 
514  void reinit (const size_type m,
515  const size_type n,
516  const std::vector<unsigned int> &row_lengths);
517 
518 
522  void reinit (const size_type m,
523  const size_type n,
524  const VectorSlice<const std::vector<unsigned int> > &row_lengths);
525 
538  void compress ();
539 
540 
617  template <typename ForwardIterator>
618  void copy_from (const size_type n_rows,
619  const size_type n_cols,
620  const ForwardIterator begin,
621  const ForwardIterator end);
622 
630  template <typename SparsityPatternType>
631  void copy_from (const SparsityPatternType &dsp);
632 
633 
641  template <typename number>
642  void copy_from (const FullMatrix<number> &matrix);
643 
651  void symmetrize ();
652 
659  void add (const size_type i,
660  const size_type j);
661 
668  template <typename ForwardIterator>
669  void add_entries (const size_type row,
670  ForwardIterator begin,
671  ForwardIterator end,
672  const bool indices_are_sorted = false);
673 
674 // @}
675 
676 
677 
678 
682 // @{
683 
692  iterator begin () const;
693 
697  iterator end () const;
698 
710  iterator begin (const size_type r) const;
711 
720  iterator end (const size_type r) const;
721 
722 
723 // @}
727 // @{
731  bool operator == (const SparsityPattern &) const;
732 
737  bool empty () const;
738 
744  size_type max_entries_per_row () const;
745 
755  size_type bandwidth () const;
756 
765  size_type n_nonzero_elements () const;
766 
770  bool is_compressed () const;
771 
776  size_type n_rows () const;
777 
782  size_type n_cols () const;
783 
787  unsigned int row_length (const size_type row) const;
788 
801  bool stores_only_added_elements () const;
802 
807  std::size_t memory_consumption () const;
808 
809 // @}
813 // @{
836  size_type operator() (const size_type i,
837  const size_type j) const;
838 
850  std::pair<size_type, size_type>
851  matrix_position (const size_type global_index) const;
852 
856  bool exists (const size_type i,
857  const size_type j) const;
858 
866  size_type row_position(const size_type i,
867  const size_type j) const;
868 
879  size_type column_number (const size_type row,
880  const unsigned int index) const;
881 
882 
883 // @}
887 // @{
898  void block_write (std::ostream &out) const;
899 
913  void block_read (std::istream &in);
914 
920  void print (std::ostream &out) const;
921 
935  void print_gnuplot (std::ostream &out) const;
936 
944  void print_svg (std::ostream &out) const;
945 
946 
951  template <class Archive>
952  void save (Archive &ar, const unsigned int version) const;
953 
958  template <class Archive>
959  void load (Archive &ar, const unsigned int version);
960 
961  BOOST_SERIALIZATION_SPLIT_MEMBER()
962 
963 // @}
964 
965 
972  DeclException2 (ExcNotEnoughSpace,
973  int, int,
974  << "Upon entering a new entry to row " << arg1
975  << ": there was no free entry any more. " << std::endl
976  << "(Maximum number of entries for this row: "
977  << arg2 << "; maybe the matrix is already compressed?)");
982  DeclException0 (ExcNotCompressed);
987  DeclException0 (ExcMatrixIsCompressed);
991  DeclException0 (ExcInvalidConstructorCall);
997  DeclException0 (ExcDiagonalNotOptimized);
1001  DeclException2 (ExcIteratorRange,
1002  int, int,
1003  << "The iterators denote a range of " << arg1
1004  << " elements, but the given number of rows was " << arg2);
1008  DeclException1 (ExcInvalidNumberOfPartitions,
1009  int,
1010  << "The number of partitions you gave is " << arg1
1011  << ", but must be greater than zero.");
1013 private:
1021  size_type max_dim;
1022 
1026  size_type rows;
1027 
1031  size_type cols;
1032 
1038  size_type max_vec_len;
1039 
1047  unsigned int max_row_length;
1048 
1062  std::size_t *rowstart;
1063 
1086  size_type *colnums;
1087 
1091  bool compressed;
1092 
1097 
1101  template <typename number> friend class SparseMatrix;
1102  template <typename number> friend class SparseLUDecomposition;
1103  template <typename number> friend class SparseILU;
1104  template <typename number> friend class ChunkSparseMatrix;
1105 
1106  friend class ChunkSparsityPattern;
1107 
1111  friend class SparsityPatternIterators::Iterator;
1112  friend class SparsityPatternIterators::Accessor;
1113  friend class ChunkSparsityPatternIterators::Accessor;
1114 };
1115 
1116 
1118 /*---------------------- Inline functions -----------------------------------*/
1119 
1120 #ifndef DOXYGEN
1121 
1122 
1123 namespace SparsityPatternIterators
1124 {
1125  inline
1126  Accessor::
1127  Accessor (const SparsityPattern *sparsity_pattern,
1128  const std::size_t i)
1129  :
1130  sparsity_pattern(sparsity_pattern),
1131  index_within_sparsity(i)
1132  {}
1133 
1134 
1135  inline
1136  Accessor::
1137  Accessor (const SparsityPattern *sparsity_pattern)
1138  :
1139  sparsity_pattern(sparsity_pattern),
1140  index_within_sparsity(sparsity_pattern->rowstart[sparsity_pattern->rows])
1141  {}
1142 
1143 
1144  inline
1145  bool
1146  Accessor::is_valid_entry () const
1147  {
1148  return (index_within_sparsity < sparsity_pattern->rowstart[sparsity_pattern->rows]
1149  &&
1150  sparsity_pattern->colnums[index_within_sparsity]
1152  }
1153 
1154 
1155  inline
1156  size_type
1157  Accessor::row() const
1158  {
1159  Assert (is_valid_entry() == true, ExcInvalidIterator());
1160 
1161  const std::size_t *insert_point =
1162  std::upper_bound(sparsity_pattern->rowstart,
1163  sparsity_pattern->rowstart + sparsity_pattern->rows + 1,
1164  index_within_sparsity);
1165  return insert_point - sparsity_pattern->rowstart - 1;
1166  }
1167 
1168 
1169  inline
1170  size_type
1171  Accessor::column() const
1172  {
1173  Assert (is_valid_entry() == true, ExcInvalidIterator());
1174 
1175  return (sparsity_pattern->colnums[index_within_sparsity]);
1176  }
1177 
1178 
1179  inline
1180  size_type
1181  Accessor::index() const
1182  {
1183  Assert (is_valid_entry() == true, ExcInvalidIterator());
1184 
1185  return index_within_sparsity - sparsity_pattern->rowstart[row()];
1186  }
1187 
1188 
1189 
1190 
1191  inline
1192  bool
1193  Accessor::operator == (const Accessor &other) const
1194  {
1195  return (sparsity_pattern == other.sparsity_pattern &&
1196  index_within_sparsity == other.index_within_sparsity);
1197  }
1198 
1199 
1200 
1201  inline
1202  bool
1203  Accessor::operator < (const Accessor &other) const
1204  {
1205  Assert (sparsity_pattern == other.sparsity_pattern,
1206  ExcInternalError());
1207 
1208  return index_within_sparsity < other.index_within_sparsity;
1209  }
1210 
1211 
1212  inline
1213  void
1214  Accessor::advance ()
1215  {
1216  Assert (index_within_sparsity < sparsity_pattern->rowstart[sparsity_pattern->rows],
1217  ExcIteratorPastEnd());
1218  ++index_within_sparsity;
1219  }
1220 
1221 
1222 
1223  inline
1224  Iterator::Iterator (const SparsityPattern *sparsity_pattern,
1225  const std::size_t i)
1226  :
1227  accessor(sparsity_pattern, i)
1228  {}
1229 
1230 
1231 
1232  inline
1233  Iterator &
1234  Iterator::operator++ ()
1235  {
1236  accessor.advance ();
1237  return *this;
1238  }
1239 
1240 
1241 
1242  inline
1243  Iterator
1244  Iterator::operator++ (int)
1245  {
1246  const Iterator iter = *this;
1247  accessor.advance ();
1248  return iter;
1249  }
1250 
1251 
1252 
1253  inline
1254  const Accessor &
1255  Iterator::operator* () const
1256  {
1257  return accessor;
1258  }
1259 
1260 
1261 
1262  inline
1263  const Accessor *
1264  Iterator::operator-> () const
1265  {
1266  return &accessor;
1267  }
1268 
1269 
1270  inline
1271  bool
1272  Iterator::operator == (const Iterator &other) const
1273  {
1274  return (accessor == other.accessor);
1275  }
1276 
1277 
1278 
1279  inline
1280  bool
1281  Iterator::operator != (const Iterator &other) const
1282  {
1283  return ! (*this == other);
1284  }
1285 
1286 
1287  inline
1288  bool
1289  Iterator::operator < (const Iterator &other) const
1290  {
1291  return accessor < other.accessor;
1292  }
1293 
1294 
1295  inline
1296  int
1297  Iterator::operator - (const Iterator &other) const
1298  {
1299  Assert (accessor.sparsity_pattern == other.accessor.sparsity_pattern,
1300  ExcInternalError());
1301 
1302  return (*this)->index_within_sparsity - other->index_within_sparsity;
1303  }
1304 }
1305 
1306 
1307 
1308 inline
1310 SparsityPattern::begin () const
1311 {
1312  return iterator(this, rowstart[0]);
1313 }
1314 
1315 
1316 inline
1318 SparsityPattern::end () const
1319 {
1320  return iterator(this, rowstart[rows]);
1321 }
1322 
1323 
1324 
1325 inline
1327 SparsityPattern::begin (const size_type r) const
1328 {
1329  Assert (r<n_rows(), ExcIndexRangeType<size_type>(r,0,n_rows()));
1330 
1331  return iterator(this, rowstart[r]);
1332 }
1333 
1334 
1335 
1336 inline
1338 SparsityPattern::end (const size_type r) const
1339 {
1340  Assert (r<n_rows(), ExcIndexRangeType<size_type>(r,0,n_rows()));
1341 
1342  return iterator(this, rowstart[r+1]);
1343 }
1344 
1345 
1346 
1347 inline
1349 SparsityPattern::n_rows () const
1350 {
1351  return rows;
1352 }
1353 
1354 
1355 inline
1357 SparsityPattern::n_cols () const
1358 {
1359  return cols;
1360 }
1361 
1362 
1363 inline
1364 bool
1366 {
1367  return compressed;
1368 }
1369 
1370 
1371 inline
1372 bool
1374 {
1375  return (store_diagonal_first_in_row == false);
1376 }
1377 
1378 
1379 
1380 inline
1381 unsigned int
1382 SparsityPattern::row_length (const size_type row) const
1383 {
1384  Assert(row<rows, ExcIndexRangeType<size_type>(row,0,rows));
1385  return rowstart[row+1]-rowstart[row];
1386 }
1387 
1388 
1389 
1390 inline
1392 SparsityPattern::column_number (const size_type row,
1393  const unsigned int index) const
1394 {
1395  Assert(row<rows, ExcIndexRangeType<size_type>(row,0,rows));
1396  Assert(index<row_length(row), ExcIndexRange(index,0,row_length(row)));
1397 
1398  return colnums[rowstart[row]+index];
1399 }
1400 
1401 
1402 inline
1405 {
1406  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
1407  Assert (compressed, ExcNotCompressed());
1408  return rowstart[rows]-rowstart[0];
1409 }
1410 
1411 
1412 
1413 template <class Archive>
1414 inline
1415 void
1416 SparsityPattern::save (Archive &ar, const unsigned int) const
1417 {
1418  // forward to serialization function in the base class.
1419  ar &static_cast<const Subscriptor &>(*this);
1420 
1421  ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed &store_diagonal_first_in_row;
1422 
1423  ar &boost::serialization::make_array(rowstart, max_dim + 1);
1424  ar &boost::serialization::make_array(colnums, max_vec_len);
1425 }
1426 
1427 
1428 
1429 template <class Archive>
1430 inline
1431 void
1432 SparsityPattern::load (Archive &ar, const unsigned int)
1433 {
1434  // forward to serialization function in the base class.
1435  ar &static_cast<Subscriptor &>(*this);
1436 
1437  ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed &store_diagonal_first_in_row;
1438 
1439  rowstart = new std::size_t [max_dim + 1];
1440  colnums = new size_type [max_vec_len];
1441 
1442  ar &boost::serialization::make_array(rowstart, max_dim + 1);
1443  ar &boost::serialization::make_array(colnums, max_vec_len);
1444 }
1445 
1446 
1447 
1448 inline
1449 bool
1451 {
1452  // it isn't quite necessary to compare *all* member variables. by only
1453  // comparing the essential ones, we can say that two sparsity patterns are
1454  // equal even if one is compressed and the other is not (in which case some
1455  // of the member variables are not yet set correctly)
1456  if (rows != sp2.rows ||
1457  cols != sp2.cols ||
1458  compressed != sp2.compressed ||
1459  store_diagonal_first_in_row != sp2.store_diagonal_first_in_row)
1460  return false;
1461 
1462  for (size_type i = 0; i < rows+1; ++i)
1463  if (rowstart[i] != sp2.rowstart[i])
1464  return false;
1465 
1466  for (size_type i = 0; i < rowstart[rows]; ++i)
1467  if (colnums[i] != sp2.colnums[i])
1468  return false;
1469 
1470  return true;
1471 }
1472 
1473 
1474 
1475 namespace internal
1476 {
1477  namespace SparsityPatternTools
1478  {
1483 
1484  inline
1485  size_type
1486  get_column_index_from_iterator (const size_type i)
1487  {
1488  return i;
1489  }
1490 
1491 
1492 
1493  template <typename value>
1494  inline
1495  size_type
1496  get_column_index_from_iterator (const std::pair<size_type, value> &i)
1497  {
1498  return i.first;
1499  }
1500 
1501 
1502 
1503  template <typename value>
1504  inline
1505  size_type
1506  get_column_index_from_iterator (const std::pair<const size_type, value> &i)
1507  {
1508  return i.first;
1509  }
1510  }
1511 }
1512 
1513 
1514 
1515 template <typename ForwardIterator>
1516 void
1517 SparsityPattern::copy_from (const size_type n_rows,
1518  const size_type n_cols,
1519  const ForwardIterator begin,
1520  const ForwardIterator end)
1521 {
1522  Assert (static_cast<size_type>(std::distance (begin, end)) == n_rows,
1523  ExcIteratorRange (std::distance (begin, end), n_rows));
1524 
1525  // first determine row lengths for each row. if the matrix is quadratic,
1526  // then we might have to add an additional entry for the diagonal, if that
1527  // is not yet present. as we have to call compress anyway later on, don't
1528  // bother to check whether that diagonal entry is in a certain row or not
1529  const bool is_square = (n_rows == n_cols);
1530  std::vector<unsigned int> row_lengths;
1531  row_lengths.reserve(n_rows);
1532  for (ForwardIterator i=begin; i!=end; ++i)
1533  row_lengths.push_back (std::distance (i->begin(), i->end())
1534  +
1535  (is_square ? 1 : 0));
1536  reinit (n_rows, n_cols, row_lengths);
1537 
1538  // now enter all the elements into the matrix. note that if the matrix is
1539  // quadratic, then we already have the diagonal element preallocated
1540  //
1541  // for use in the inner loop, we define a typedef to the type of the inner
1542  // iterators
1543  size_type row = 0;
1544  typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
1545  for (ForwardIterator i=begin; i!=end; ++i, ++row)
1546  {
1547  size_type *cols = &colnums[rowstart[row]] + (is_square ? 1 : 0);
1548  const inner_iterator end_of_row = i->end();
1549  for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
1550  {
1551  const size_type col
1552  = internal::SparsityPatternTools::get_column_index_from_iterator(*j);
1553  Assert (col < n_cols, ExcIndexRange(col,0,n_cols));
1554 
1555  if ((col!=row) || !is_square)
1556  *cols++ = col;
1557  }
1558  }
1559 
1560  // finally compress everything. this also sorts the entries within each row
1561  compress ();
1562 }
1563 
1564 
1565 #endif // DOXYGEN
1566 
1567 DEAL_II_NAMESPACE_CLOSE
1568 
1569 #endif
const types::global_dof_index invalid_size_type
Definition: types.h:173
SparsityPatternIterators::Iterator iterator
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
int operator-(const Iterator &p) const
void print_svg(std::ostream &out) const
std::size_t memory_consumption() const
bool operator<(const Accessor &) const
types::global_dof_index size_type
static const size_type invalid_entry
bool operator==(const Accessor &) const
void print_gnuplot(std::ostream &out) const
size_type row_position(const size_type i, const size_type j) const
void save(Archive &ar, const unsigned int version) const
iterator begin() const
size_type operator()(const size_type i, const size_type j) const
DeclException1(ExcInvalidNumberOfPartitions, int,<< "The number of partitions you gave is "<< arg1<< ", but must be greater than zero.")
STL namespace.
size_type n_nonzero_elements() const
DeclException0(ExcNotCompressed)
std::pair< size_type, size_type > matrix_position(const size_type global_index) const
const Accessor & operator*() const
size_type max_entries_per_row() const
size_type bandwidth() const
void add(const size_type i, const size_type j)
unsigned int global_dof_index
Definition: types.h:88
void block_write(std::ostream &out) const
#define Assert(cond, exc)
Definition: exceptions.h:294
size_type n_cols() const
unsigned int row_length(const size_type row) const
size_type column_number(const size_type row, const unsigned int index) const
bool operator==(const Iterator &) const
SparsityPattern & operator=(const SparsityPattern &)
bool stores_only_added_elements() const
iterator end() const
void load(Archive &ar, const unsigned int version)
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
SparsityPatternIterators::Iterator const_iterator
DeclException2(ExcNotEnoughSpace, int, int,<< "Upon entering a new entry to row "<< arg1<< ": there was no free entry any more. "<< std::endl<< "(Maximum number of entries for this row: "<< arg2<< "; maybe the matrix is already compressed?)")
void reinit(const size_type m, const size_type n, const unsigned int max_per_row)
std::size_t * rowstart
unsigned int max_row_length
types::global_dof_index size_type
const SparsityPattern * sparsity_pattern
size_type n_rows() const
void print(std::ostream &out) const
bool operator<(const Iterator &) const
Iterator(const SparsityPattern *sp, const std::size_t index_within_sparsity)
void block_read(std::istream &in)
Accessor(const SparsityPattern *matrix, const std::size_t index_within_sparsity)
bool exists(const size_type i, const size_type j) const
bool empty() const
bool operator==(const SparsityPattern &) const
bool is_compressed() const
bool operator!=(const Iterator &) const
const Accessor * operator->() const