Reference documentation for deal.II version GIT e22cb6a53e 2023-09-21 14:00:02+00:00
\(\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\}}\)
sparsity_pattern.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2022 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_sparsity_pattern_h
17 #define dealii_sparsity_pattern_h
18 
19 
20 #include <deal.II/base/config.h>
21 
26 
28 
29 // boost::serialization::make_array used to be in array.hpp, but was
30 // moved to a different file in BOOST 1.64
31 #include <boost/version.hpp>
32 #if BOOST_VERSION >= 106400
33 # include <boost/serialization/array_wrapper.hpp>
34 #else
35 # include <boost/serialization/array.hpp>
36 #endif
37 #include <boost/serialization/split_member.hpp>
38 
39 #include <algorithm>
40 #include <iostream>
41 #include <memory>
42 #include <vector>
43 
45 
46 // Forward declarations
47 #ifndef DOXYGEN
48 class SparsityPattern;
51 template <typename number>
52 class FullMatrix;
53 template <typename number>
54 class SparseMatrix;
55 template <typename number>
57 template <typename number>
58 class SparseILU;
59 
61 {
62  class Accessor;
63 }
64 #endif
65 
71 namespace internals
72 {
73  namespace SparsityPatternTools
74  {
79 
85  size_type
87 
93  template <typename value>
94  size_type
95  get_column_index_from_iterator(const std::pair<size_type, value> &i);
96 
102  template <typename value>
103  size_type
104  get_column_index_from_iterator(const std::pair<const size_type, value> &i);
105 
106  } // namespace SparsityPatternTools
107 } // namespace internals
108 
109 
114 {
115  // forward declaration
116  class Iterator;
117 
122 
133  class Accessor
134  {
135  public:
140 
144  Accessor(const SparsityPattern *matrix, const std::size_t linear_index);
145 
150 
157 
162  size_type
163  row() const;
164 
170  size_type
171  index() const;
172 
182  size_type
183  global_index() const;
184 
189  size_type
190  column() const;
191 
201  bool
202  is_valid_entry() const;
203 
207  bool
208  operator==(const Accessor &) const;
209 
217  bool
218  operator<(const Accessor &) const;
219 
220  protected:
222  "The instance of this class was initialized"
223  " without SparsityPattern object, which"
224  " means that it is a dummy accessor that can"
225  " not do any operations.");
226 
231 
239  std::size_t linear_index;
240 
244  void
246 
247  // Grant access to iterator class.
248  friend class LinearIndexIterator<Iterator, Accessor>;
249 
250  // Grant access to accessor class of ChunkSparsityPattern.
252  };
253 
254 
255 
280  class Iterator : public LinearIndexIterator<Iterator, Accessor>
281  {
282  public:
287 
292 
298  Iterator(const SparsityPattern *sp, const std::size_t linear_index);
299 
305  };
306 } // namespace SparsityPatternIterators
307 
344 {
345 public:
350 
366 
372 
381 
396  iterator
397  begin() const;
398 
402  iterator
403  end() const;
404 
416  iterator
417  begin(const size_type r) const;
418 
427  iterator
428  end(const size_type r) const;
429 
449  SparsityPattern();
450 
467 
475  SparsityPattern(const size_type m,
476  const size_type n,
477  const unsigned int max_per_row);
478 
479 
488  SparsityPattern(const size_type m,
489  const size_type n,
490  const std::vector<unsigned int> &row_lengths);
491 
500  SparsityPattern(const size_type m, const unsigned int max_per_row);
501 
509  SparsityPattern(const size_type m,
510  const std::vector<unsigned int> &row_lengths);
511 
534  SparsityPattern(const SparsityPattern &original,
535  const unsigned int max_per_row,
536  const size_type extra_off_diagonals);
537 
544  operator=(const SparsityPattern &);
545 
556  void
557  reinit(const size_type m,
558  const size_type n,
559  const ArrayView<const unsigned int> &row_lengths);
560 
564  void
565  reinit(const size_type m,
566  const size_type n,
567  const std::vector<unsigned int> &row_lengths);
568 
577  void
578  reinit(const size_type m, const size_type n, const unsigned int max_per_row);
579 
592  void
593  compress();
594 
595 
672  template <typename ForwardIterator>
673  void
675  const size_type n_cols,
676  const ForwardIterator begin,
677  const ForwardIterator end);
678 
683  void
684  copy_from(const DynamicSparsityPattern &dsp);
685 
690  void
691  copy_from(const SparsityPattern &sp);
692 
707  template <typename number>
708  void
710 
717  void
718  add(const size_type i, const size_type j);
719 
726  template <typename ForwardIterator>
727  void
728  add_entries(const size_type row,
729  ForwardIterator begin,
730  ForwardIterator end,
731  const bool indices_are_sorted = false);
732 
733  virtual void
734  add_row_entries(const size_type &row,
735  const ArrayView<const size_type> &columns,
736  const bool indices_are_sorted = false) override;
737 
739 
747  void
748  symmetrize();
749 
770  std::size_t
772 
777  bool
778  empty() const;
779 
783  bool
784  exists(const size_type i, const size_type j) const;
785 
797  std::pair<size_type, size_type>
798  matrix_position(const std::size_t global_index) const;
799 
809  size_type
810  bandwidth() const;
811 
815  bool
816  is_compressed() const;
817 
823  size_type
824  max_entries_per_row() const;
825 
829  bool
830  operator==(const SparsityPattern &) const;
831 
844  bool
846 
850  unsigned int
851  row_length(const size_type row) const;
852 
857  std::size_t
858  memory_consumption() const;
859 
892  size_type
893  operator()(const size_type i, const size_type j) const;
894 
905  size_type
906  column_number(const size_type row, const unsigned int index) const;
907 
915  size_type
916  row_position(const size_type i, const size_type j) const;
917 
934  void
935  print(std::ostream &out) const;
936 
950  void
951  print_gnuplot(std::ostream &out) const;
952 
960  void
961  print_svg(std::ostream &out) const;
962 
973  void
974  block_write(std::ostream &out) const;
975 
989  void
990  block_read(std::istream &in);
991 
996  template <class Archive>
997  void
998  save(Archive &ar, const unsigned int version) const;
999 
1004  template <class Archive>
1005  void
1006  load(Archive &ar, const unsigned int version);
1007 
1008 #ifdef DOXYGEN
1013  template <class Archive>
1014  void
1015  serialize(Archive &archive, const unsigned int version);
1016 #else
1017  // This macro defines the serialize() method that is compatible with
1018  // the templated save() and load() method that have been implemented.
1019  BOOST_SERIALIZATION_SPLIT_MEMBER()
1020 #endif
1021 
1035  int,
1036  int,
1037  << "The iterators denote a range of " << arg1
1038  << " elements, but the given number of rows was " << arg2);
1043  int,
1044  << "The number of partitions you gave is " << arg1
1045  << ", but must be greater than zero.");
1046 
1051  int,
1052  int,
1053  << "Upon entering a new entry to row " << arg1
1054  << ": there was no free entry any more. " << std::endl
1055  << "(Maximum number of entries for this row: " << arg2
1056  << "; maybe the matrix is already compressed?)");
1057 
1064  "The operation you attempted changes the structure of the SparsityPattern "
1065  "and is not possible after compress() has been called.");
1066 
1073  "The operation you attempted is only allowed after the SparsityPattern "
1074  "has been set up and compress() was called.");
1075 
1079 private:
1084 
1093 
1099  std::size_t max_vec_len;
1100 
1108  unsigned int max_row_length;
1109 
1123  std::unique_ptr<std::size_t[]> rowstart;
1124 
1147  std::unique_ptr<size_type[]> colnums;
1148 
1153 
1154  // Make all sparse matrices friends of this class.
1155  template <typename number>
1156  friend class SparseMatrix;
1157  template <typename number>
1159  template <typename number>
1160  friend class SparseILU;
1161  template <typename number>
1162  friend class ChunkSparseMatrix;
1163 
1164  friend class ChunkSparsityPattern;
1166 
1167  // Also give access to internal details to the iterator/accessor classes.
1171 };
1172 
1173 
1177 /*---------------------- Inline functions -----------------------------------*/
1178 
1179 #ifndef DOXYGEN
1180 
1181 
1182 namespace SparsityPatternIterators
1183 {
1184  inline Accessor::Accessor(const SparsityPattern *sparsity_pattern,
1185  const std::size_t i)
1186  : container(sparsity_pattern)
1187  , linear_index(i)
1188  {}
1189 
1190 
1191 
1192  inline Accessor::Accessor(const SparsityPattern *sparsity_pattern)
1193  : container(sparsity_pattern)
1194  , linear_index(container->rowstart[container->rows])
1195  {}
1196 
1197 
1198 
1199  inline Accessor::Accessor()
1200  : container(nullptr)
1201  , linear_index(numbers::invalid_size_type)
1202  {}
1203 
1204 
1205 
1206  inline bool
1207  Accessor::is_valid_entry() const
1208  {
1209  Assert(container != nullptr, DummyAccessor());
1210  return (linear_index < container->rowstart[container->rows] &&
1211  container->colnums[linear_index] != SparsityPattern::invalid_entry);
1212  }
1213 
1214 
1215 
1216  inline size_type
1217  Accessor::row() const
1218  {
1219  Assert(is_valid_entry() == true, ExcInvalidIterator());
1220 
1221  const std::size_t *insert_point =
1222  std::upper_bound(container->rowstart.get(),
1223  container->rowstart.get() + container->rows + 1,
1224  linear_index);
1225  return insert_point - container->rowstart.get() - 1;
1226  }
1227 
1228 
1229 
1230  inline size_type
1231  Accessor::column() const
1232  {
1233  Assert(is_valid_entry() == true, ExcInvalidIterator());
1234 
1235  return (container->colnums[linear_index]);
1236  }
1237 
1238 
1239 
1240  inline size_type
1241  Accessor::index() const
1242  {
1243  Assert(is_valid_entry() == true, ExcInvalidIterator());
1244 
1245  return linear_index - container->rowstart[row()];
1246  }
1247 
1248 
1249 
1250  inline size_type
1251  Accessor::global_index() const
1252  {
1253  Assert(is_valid_entry() == true, ExcInvalidIterator());
1254 
1255  return linear_index;
1256  }
1257 
1258 
1259 
1260  inline bool
1261  Accessor::operator==(const Accessor &other) const
1262  {
1263  return (container == other.container && linear_index == other.linear_index);
1264  }
1265 
1266 
1267 
1268  inline bool
1269  Accessor::operator<(const Accessor &other) const
1270  {
1271  Assert(container != nullptr, DummyAccessor());
1272  Assert(other.container != nullptr, DummyAccessor());
1273  Assert(container == other.container, ExcInternalError());
1274 
1275  return linear_index < other.linear_index;
1276  }
1277 
1278 
1279 
1280  inline void
1282  {
1283  Assert(container != nullptr, DummyAccessor());
1284  Assert(linear_index < container->rowstart[container->rows],
1285  ExcIteratorPastEnd());
1286  ++linear_index;
1287  }
1288 
1289 
1290  inline Iterator::Iterator(const SparsityPattern *sp,
1291  const std::size_t linear_index)
1292  : LinearIndexIterator<Iterator, Accessor>(Accessor(sp, linear_index))
1293  {}
1294 
1295 
1296  inline Iterator::Iterator(const Accessor &accessor)
1297  : LinearIndexIterator<Iterator, Accessor>(accessor)
1298  {}
1299 
1300 
1301 } // namespace SparsityPatternIterators
1302 
1303 
1304 
1305 inline std::size_t
1307 {
1309 
1310  if ((rowstart != nullptr) && (colnums != nullptr))
1311  return rowstart[rows] - rowstart[0];
1312  else
1313  // the object is empty or has zero size
1314  return 0;
1315 }
1316 
1317 
1318 
1319 inline bool
1321 {
1322  return compressed;
1323 }
1324 
1325 
1326 
1327 inline bool
1329 {
1330  return (store_diagonal_first_in_row == false);
1331 }
1332 
1333 
1334 
1335 inline unsigned int
1336 SparsityPattern::row_length(const size_type row) const
1337 {
1338  AssertIndexRange(row, rows);
1339  return rowstart[row + 1] - rowstart[row];
1340 }
1341 
1342 
1343 
1346  const unsigned int index) const
1347 {
1348  AssertIndexRange(row, rows);
1349  AssertIndexRange(index, row_length(row));
1350 
1351  return colnums[rowstart[row] + index];
1352 }
1353 
1354 
1355 
1357 SparsityPattern::begin() const
1358 {
1359  if (n_rows() > 0)
1360  return {this, rowstart[0]};
1361  else
1362  return end();
1363 }
1364 
1365 
1366 
1368 SparsityPattern::end() const
1369 {
1370  if (n_rows() > 0)
1371  return {this, rowstart[rows]};
1372  else
1373  return {nullptr, 0};
1374 }
1375 
1376 
1377 
1379 SparsityPattern::begin(const size_type r) const
1380 {
1381  AssertIndexRange(r, n_rows());
1382 
1383  return {this, rowstart[r]};
1384 }
1385 
1386 
1387 
1389 SparsityPattern::end(const size_type r) const
1390 {
1391  AssertIndexRange(r, n_rows());
1392 
1393  return {this, rowstart[r + 1]};
1394 }
1395 
1396 
1397 
1398 inline bool
1400 {
1402  return false;
1403 
1404  // it isn't quite necessary to compare *all* member variables. by only
1405  // comparing the essential ones, we can say that two sparsity patterns are
1406  // equal even if one is compressed and the other is not (in which case some
1407  // of the member variables are not yet set correctly)
1408  if (rows != sp2.rows || cols != sp2.cols || compressed != sp2.compressed)
1409  return false;
1410 
1411  if (rows > 0)
1412  for (size_type i = 0; i < rows + 1; ++i)
1413  if (rowstart[i] != sp2.rowstart[i])
1414  return false;
1415 
1416  if (rows > 0)
1417  for (size_type i = 0; i < rowstart[rows]; ++i)
1418  if (colnums[i] != sp2.colnums[i])
1419  return false;
1420 
1421  return true;
1422 }
1423 
1424 
1425 
1426 namespace internal
1427 {
1428  namespace SparsityPatternTools
1429  {
1434 
1435  inline size_type
1437  {
1438  return i;
1439  }
1440 
1441 
1442 
1443  template <typename value>
1444  inline size_type
1445  get_column_index_from_iterator(const std::pair<size_type, value> &i)
1446  {
1447  return i.first;
1448  }
1449 
1450 
1451 
1452  template <typename value>
1453  inline size_type
1454  get_column_index_from_iterator(const std::pair<const size_type, value> &i)
1455  {
1456  return i.first;
1457  }
1458  } // namespace SparsityPatternTools
1459 } // namespace internal
1460 
1461 
1462 
1463 template <typename ForwardIterator>
1464 void
1466  const size_type n_cols,
1467  const ForwardIterator begin,
1468  const ForwardIterator end)
1469 {
1470  Assert(static_cast<size_type>(std::distance(begin, end)) == n_rows,
1471  ExcIteratorRange(std::distance(begin, end), n_rows));
1472 
1473  // first determine row lengths for each row. if the matrix is quadratic,
1474  // then we might have to add an additional entry for the diagonal, if that
1475  // is not yet present. as we have to call compress anyway later on, don't
1476  // bother to check whether that diagonal entry is in a certain row or not
1477  const bool is_square = (n_rows == n_cols);
1478  std::vector<unsigned int> row_lengths;
1479  row_lengths.reserve(n_rows);
1480  for (ForwardIterator i = begin; i != end; ++i)
1481  row_lengths.push_back(std::distance(i->begin(), i->end()) +
1482  (is_square ? 1 : 0));
1483  reinit(n_rows, n_cols, row_lengths);
1484 
1485  // now enter all the elements into the matrix. note that if the matrix is
1486  // quadratic, then we already have the diagonal element preallocated
1487  //
1488  // for use in the inner loop, we define an alias to the type of the inner
1489  // iterators
1490  size_type row = 0;
1491  using inner_iterator =
1492  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1493  for (ForwardIterator i = begin; i != end; ++i, ++row)
1494  {
1495  size_type *cols = &colnums[rowstart[row]] + (is_square ? 1 : 0);
1496  const inner_iterator end_of_row = i->end();
1497  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1498  {
1499  const size_type col =
1501  AssertIndexRange(col, n_cols);
1502 
1503  if ((col != row) || !is_square)
1504  *cols++ = col;
1505  }
1506  }
1507 
1508  // finally compress everything. this also sorts the entries within each row
1509  compress();
1510 }
1511 
1512 
1513 
1514 template <class Archive>
1515 inline void
1516 SparsityPattern::save(Archive &ar, const unsigned int) const
1517 {
1518  // forward to serialization function in the base class.
1519  ar &boost::serialization::base_object<const Subscriptor>(*this);
1520 
1522 
1523  if (max_dim != 0)
1524  ar &boost::serialization::make_array(rowstart.get(), max_dim + 1);
1525  else
1526  Assert(rowstart.get() == nullptr, ExcInternalError());
1527 
1528  if (max_vec_len != 0)
1529  ar &boost::serialization::make_array(colnums.get(), max_vec_len);
1530  else
1531  Assert(colnums.get() == nullptr, ExcInternalError());
1533 }
1534 
1535 
1536 
1537 template <class Archive>
1538 inline void
1539 SparsityPattern::load(Archive &ar, const unsigned int)
1540 {
1541  // forward to serialization function in the base class.
1542  ar &boost::serialization::base_object<Subscriptor>(*this);
1543 
1545 
1546  if (max_dim != 0)
1547  {
1548  rowstart = std::make_unique<std::size_t[]>(max_dim + 1);
1549  ar &boost::serialization::make_array(rowstart.get(), max_dim + 1);
1550  }
1551  else
1552  rowstart.reset();
1553 
1554  if (max_vec_len != 0)
1555  {
1556  colnums = std::make_unique<size_type[]>(max_vec_len);
1557  ar &boost::serialization::make_array(colnums.get(), max_vec_len);
1558  }
1559  else
1560  colnums.reset();
1562 }
1563 
1564 
1565 #endif // DOXYGEN
1566 
1568 
1569 #endif
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
size_type n_rows() const
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type >> &entries)
size_type n_cols() const
bool operator==(const Accessor &) const
SparsityPatternIterators::size_type size_type
const SparsityPattern * container
Accessor(const SparsityPattern *matrix)
Accessor(const SparsityPattern *matrix, const std::size_t linear_index)
bool operator<(const Accessor &) const
Iterator(const SparsityPattern *sp, const std::size_t linear_index)
types::global_dof_index size_type
Iterator(const Accessor &accessor)
void load(Archive &ar, const unsigned int version)
std::pair< size_type, size_type > matrix_position(const std::size_t global_index) const
iterator begin(const size_type r) const
void block_write(std::ostream &out) const
void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths)
void print_svg(std::ostream &out) const
bool stores_only_added_elements() const
void serialize(Archive &archive, const unsigned int version)
size_type bandwidth() const
void print_gnuplot(std::ostream &out) const
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false) override
bool is_compressed() const
SparsityPattern & operator=(const SparsityPattern &)
std::size_t n_nonzero_elements() const
bool exists(const size_type i, const size_type j) const
iterator begin() const
std::unique_ptr< size_type[]> colnums
std::size_t max_vec_len
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type >> &entries)
size_type column_number(const size_type row, const unsigned int index) const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
static constexpr size_type invalid_entry
size_type row_position(const size_type i, const size_type j) const
std::unique_ptr< std::size_t[]> rowstart
void print(std::ostream &out) const
iterator end(const size_type r) const
void add(const size_type i, const size_type j)
size_type max_entries_per_row() const
bool empty() const
unsigned int max_row_length
iterator end() const
size_type operator()(const size_type i, const size_type j) const
void save(Archive &ar, const unsigned int version) const
types::global_dof_index size_type
unsigned int row_length(const size_type row) const
std::size_t memory_consumption() const
bool operator==(const SparsityPattern &) const
void block_read(std::istream &in)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotCompressed()
static ::ExceptionBase & ExcNotEnoughSpace(int arg1, int arg2)
static ::ExceptionBase & ExcIteratorPastEnd()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcInvalidIterator()
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:535
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:490
static ::ExceptionBase & ExcMatrixIsCompressed()
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:512
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static ::ExceptionBase & DummyAccessor()
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
types::global_dof_index size_type
size_type get_column_index_from_iterator(const size_type i)
size_type get_column_index_from_iterator(const std::pair< const size_type, value > &i)
const types::global_dof_index invalid_size_type
Definition: types.h:222
unsigned int global_dof_index
Definition: types.h:82
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
void advance(std::tuple< I1, I2 > &t, const unsigned int n)