Reference documentation for deal.II version GIT bdf8bf8f35 2023-03-27 16:55:01+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\}}\)
block_matrix_base.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 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_block_matrix_base_h
17 #define dealii_block_matrix_base_h
18 
19 
20 #include <deal.II/base/config.h>
21 
23 #include <deal.II/base/mutex.h>
25 #include <deal.II/base/table.h>
26 #include <deal.II/base/utilities.h>
27 
29 #include <deal.II/lac/exceptions.h>
32 #include <deal.II/lac/vector.h>
34 
35 #include <cmath>
36 #include <mutex>
37 
39 
40 
41 // Forward declaration
42 #ifndef DOXYGEN
43 template <typename>
44 class MatrixIterator;
45 #endif
46 
47 
57 {
62  template <class BlockMatrixType>
64  {
65  public:
70 
74  using value_type = typename BlockMatrixType::value_type;
75 
80 
84  unsigned int
85  block_row() const;
86 
90  unsigned int
91  block_column() const;
92 
93  protected:
97  unsigned int row_block;
98 
102  unsigned int col_block;
103 
104  // Let the iterator class be a friend.
105  template <typename>
106  friend class MatrixIterator;
107  };
108 
109 
110 
114  template <class BlockMatrixType, bool Constness>
115  class Accessor;
116 
117 
121  template <class BlockMatrixType>
122  class Accessor<BlockMatrixType, false> : public AccessorBase<BlockMatrixType>
123  {
124  public:
129 
133  using MatrixType = BlockMatrixType;
134 
138  using value_type = typename BlockMatrixType::value_type;
139 
148  Accessor(BlockMatrixType *m, const size_type row, const size_type col);
149 
153  size_type
154  row() const;
155 
159  size_type
160  column() const;
161 
165  value_type
166  value() const;
167 
171  void
172  set_value(value_type newval) const;
173 
174  protected:
178  BlockMatrixType *matrix;
179 
183  typename BlockMatrixType::BlockType::iterator base_iterator;
184 
188  void
190 
194  bool
195  operator==(const Accessor &a) const;
196 
197  template <typename>
198  friend class MatrixIterator;
199  friend class Accessor<BlockMatrixType, true>;
200  };
201 
206  template <class BlockMatrixType>
207  class Accessor<BlockMatrixType, true> : public AccessorBase<BlockMatrixType>
208  {
209  public:
214 
218  using MatrixType = const BlockMatrixType;
219 
223  using value_type = typename BlockMatrixType::value_type;
224 
233  Accessor(const BlockMatrixType *m,
234  const size_type row,
235  const size_type col);
236 
241 
245  size_type
246  row() const;
247 
251  size_type
252  column() const;
253 
257  value_type
258  value() const;
259 
260  protected:
264  const BlockMatrixType *matrix;
265 
269  typename BlockMatrixType::BlockType::const_iterator base_iterator;
270 
274  void
276 
280  bool
281  operator==(const Accessor &a) const;
282 
283  // Let the iterator class be a friend.
284  template <typename>
285  friend class ::MatrixIterator;
286  };
287 } // namespace BlockMatrixIterators
288 
289 
290 
350 template <typename MatrixType>
352 {
353 public:
357  using BlockType = MatrixType;
358 
363  using value_type = typename BlockType::value_type;
365  using pointer = value_type *;
366  using const_pointer = const value_type *;
368  using const_reference = const value_type &;
370 
371  using iterator =
373 
376 
377 
381  BlockMatrixBase() = default;
382 
386  ~BlockMatrixBase() override;
387 
404  template <class BlockMatrixType>
406  copy_from(const BlockMatrixType &source);
407 
411  BlockType &
412  block(const unsigned int row, const unsigned int column);
413 
414 
419  const BlockType &
420  block(const unsigned int row, const unsigned int column) const;
421 
426  size_type
427  m() const;
428 
433  size_type
434  n() const;
435 
436 
441  unsigned int
442  n_block_rows() const;
443 
448  unsigned int
449  n_block_cols() const;
450 
456  void
457  set(const size_type i, const size_type j, const value_type value);
458 
474  template <typename number>
475  void
476  set(const std::vector<size_type> &indices,
477  const FullMatrix<number> & full_matrix,
478  const bool elide_zero_values = false);
479 
485  template <typename number>
486  void
487  set(const std::vector<size_type> &row_indices,
488  const std::vector<size_type> &col_indices,
489  const FullMatrix<number> & full_matrix,
490  const bool elide_zero_values = false);
491 
502  template <typename number>
503  void
504  set(const size_type row,
505  const std::vector<size_type> &col_indices,
506  const std::vector<number> & values,
507  const bool elide_zero_values = false);
508 
518  template <typename number>
519  void
520  set(const size_type row,
521  const size_type n_cols,
522  const size_type *col_indices,
523  const number * values,
524  const bool elide_zero_values = false);
525 
531  void
532  add(const size_type i, const size_type j, const value_type value);
533 
548  template <typename number>
549  void
550  add(const std::vector<size_type> &indices,
551  const FullMatrix<number> & full_matrix,
552  const bool elide_zero_values = true);
553 
559  template <typename number>
560  void
561  add(const std::vector<size_type> &row_indices,
562  const std::vector<size_type> &col_indices,
563  const FullMatrix<number> & full_matrix,
564  const bool elide_zero_values = true);
565 
575  template <typename number>
576  void
577  add(const size_type row,
578  const std::vector<size_type> &col_indices,
579  const std::vector<number> & values,
580  const bool elide_zero_values = true);
581 
591  template <typename number>
592  void
593  add(const size_type row,
594  const size_type n_cols,
595  const size_type *col_indices,
596  const number * values,
597  const bool elide_zero_values = true,
598  const bool col_indices_are_sorted = false);
599 
611  void
613 
620  value_type
621  operator()(const size_type i, const size_type j) const;
622 
631  value_type
632  el(const size_type i, const size_type j) const;
633 
644  value_type
645  diag_element(const size_type i) const;
646 
655  void
657 
662  operator*=(const value_type factor);
663 
668  operator/=(const value_type factor);
669 
674  template <class BlockVectorType>
675  void
676  vmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
677 
683  template <class BlockVectorType>
684  void
685  Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
686 
699  template <class BlockVectorType>
700  value_type
701  matrix_norm_square(const BlockVectorType &v) const;
702 
707  real_type
708  frobenius_norm() const;
709 
713  template <class BlockVectorType>
714  value_type
715  matrix_scalar_product(const BlockVectorType &u,
716  const BlockVectorType &v) const;
717 
721  template <class BlockVectorType>
722  value_type
723  residual(BlockVectorType & dst,
724  const BlockVectorType &x,
725  const BlockVectorType &b) const;
726 
733  void
734  print(std::ostream &out, const bool alternative_output = false) const;
735 
739  iterator
740  begin();
741 
745  iterator
746  end();
747 
751  iterator
752  begin(const size_type r);
753 
757  iterator
758  end(const size_type r);
763  begin() const;
764 
769  end() const;
770 
775  begin(const size_type r) const;
776 
781  end(const size_type r) const;
782 
786  const BlockIndices &
788 
792  const BlockIndices &
794 
800  std::size_t
802 
812  int,
813  int,
814  int,
815  int,
816  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
817  << ',' << arg4 << "] have differing row numbers.");
822  int,
823  int,
824  int,
825  int,
826  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
827  << ',' << arg4 << "] have differing column numbers.");
829 protected:
842  void
843  clear();
844 
850 
855 
874  void
876 
887  template <class BlockVectorType>
888  void
889  vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const;
890 
901  template <class BlockVectorType, class VectorType>
902  void
903  vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const;
904 
915  template <class BlockVectorType, class VectorType>
916  void
917  vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const;
918 
929  template <class VectorType>
930  void
931  vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const;
932 
944  template <class BlockVectorType>
945  void
946  Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const;
947 
958  template <class BlockVectorType, class VectorType>
959  void
960  Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const;
961 
972  template <class BlockVectorType, class VectorType>
973  void
974  Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const;
975 
986  template <class VectorType>
987  void
988  Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const;
989 
990 
991 protected:
998  void
1000 
1005  void
1007 
1008 
1009 private:
1019  {
1024  std::vector<size_type> counter_within_block;
1025 
1030  std::vector<std::vector<size_type>> column_indices;
1031 
1036  std::vector<std::vector<value_type>> column_values;
1037 
1043 
1053  TemporaryData &
1055  {
1056  return *this;
1057  }
1058  };
1059 
1067 
1068  // Make the iterator class a friend. We have to work around a compiler bug
1069  // here again.
1070  template <typename, bool>
1072 
1073  template <typename>
1074  friend class MatrixIterator;
1075 };
1076 
1077 
1080 #ifndef DOXYGEN
1081 /* ------------------------- Template functions ---------------------- */
1082 
1083 
1084 namespace BlockMatrixIterators
1085 {
1086  template <class BlockMatrixType>
1088  : row_block(0)
1089  , col_block(0)
1090  {}
1091 
1092 
1093  template <class BlockMatrixType>
1094  inline unsigned int
1095  AccessorBase<BlockMatrixType>::block_row() const
1096  {
1098 
1099  return row_block;
1100  }
1101 
1102 
1103  template <class BlockMatrixType>
1104  inline unsigned int
1105  AccessorBase<BlockMatrixType>::block_column() const
1106  {
1108 
1109  return col_block;
1110  }
1111 
1112 
1113  template <class BlockMatrixType>
1114  inline Accessor<BlockMatrixType, true>::Accessor(
1115  const BlockMatrixType *matrix,
1116  const size_type row,
1117  const size_type col)
1118  : matrix(matrix)
1119  , base_iterator(matrix->block(0, 0).begin())
1120  {
1121  (void)col;
1122  Assert(col == 0, ExcNotImplemented());
1123 
1124  // check if this is a regular row or
1125  // the end of the matrix
1126  if (row < matrix->m())
1127  {
1128  const std::pair<unsigned int, size_type> indices =
1129  matrix->row_block_indices.global_to_local(row);
1130 
1131  // find the first block that does
1132  // have an entry in this row
1133  for (unsigned int bc = 0; bc < matrix->n_block_cols(); ++bc)
1134  {
1135  base_iterator =
1136  matrix->block(indices.first, bc).begin(indices.second);
1137  if (base_iterator !=
1138  matrix->block(indices.first, bc).end(indices.second))
1139  {
1140  this->row_block = indices.first;
1141  this->col_block = bc;
1142  return;
1143  }
1144  }
1145 
1146  // hm, there is no block that has
1147  // an entry in this column. we need
1148  // to take the next entry then,
1149  // which may be the first entry of
1150  // the next row, or recursively the
1151  // next row, or so on
1152  *this = Accessor(matrix, row + 1, 0);
1153  }
1154  else
1155  {
1156  // we were asked to create the end
1157  // iterator for this matrix
1158  this->row_block = numbers::invalid_unsigned_int;
1159  this->col_block = numbers::invalid_unsigned_int;
1160  }
1161  }
1162 
1163 
1164  // template <class BlockMatrixType>
1165  // inline
1166  // Accessor<BlockMatrixType, true>::Accessor (const
1167  // Accessor<BlockMatrixType, true>& other)
1168  // :
1169  // matrix(other.matrix),
1170  // base_iterator(other.base_iterator)
1171  // {
1172  // this->row_block = other.row_block;
1173  // this->col_block = other.col_block;
1174  // }
1175 
1176 
1177  template <class BlockMatrixType>
1178  inline Accessor<BlockMatrixType, true>::Accessor(
1179  const Accessor<BlockMatrixType, false> &other)
1180  : matrix(other.matrix)
1181  , base_iterator(other.base_iterator)
1182  {
1183  this->row_block = other.row_block;
1184  this->col_block = other.col_block;
1185  }
1186 
1187 
1188  template <class BlockMatrixType>
1190  Accessor<BlockMatrixType, true>::row() const
1191  {
1192  Assert(this->row_block != numbers::invalid_unsigned_int,
1193  ExcIteratorPastEnd());
1194 
1195  return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1196  base_iterator->row());
1197  }
1198 
1199 
1200  template <class BlockMatrixType>
1202  Accessor<BlockMatrixType, true>::column() const
1203  {
1204  Assert(this->col_block != numbers::invalid_unsigned_int,
1205  ExcIteratorPastEnd());
1206 
1207  return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
1208  base_iterator->column());
1209  }
1210 
1211 
1212  template <class BlockMatrixType>
1213  inline typename Accessor<BlockMatrixType, true>::value_type
1214  Accessor<BlockMatrixType, true>::value() const
1215  {
1216  Assert(this->row_block != numbers::invalid_unsigned_int,
1217  ExcIteratorPastEnd());
1218  Assert(this->col_block != numbers::invalid_unsigned_int,
1219  ExcIteratorPastEnd());
1220 
1221  return base_iterator->value();
1222  }
1223 
1224 
1225 
1226  template <class BlockMatrixType>
1227  inline void
1229  {
1230  Assert(this->row_block != numbers::invalid_unsigned_int,
1231  ExcIteratorPastEnd());
1232  Assert(this->col_block != numbers::invalid_unsigned_int,
1233  ExcIteratorPastEnd());
1234 
1235  // Remember current row inside block
1236  size_type local_row = base_iterator->row();
1237 
1238  // Advance one element inside the
1239  // current block
1240  ++base_iterator;
1241 
1242  // while we hit the end of the row of a
1243  // block (which may happen multiple
1244  // times if rows inside a block are
1245  // empty), we have to jump to the next
1246  // block and take the
1247  while (base_iterator ==
1248  matrix->block(this->row_block, this->col_block).end(local_row))
1249  {
1250  // jump to next block in this block
1251  // row, if possible, otherwise go
1252  // to next row
1253  if (this->col_block < matrix->n_block_cols() - 1)
1254  {
1255  ++this->col_block;
1256  base_iterator =
1257  matrix->block(this->row_block, this->col_block).begin(local_row);
1258  }
1259  else
1260  {
1261  // jump back to next row in
1262  // first block column
1263  this->col_block = 0;
1264  ++local_row;
1265 
1266  // see if this has brought us
1267  // past the number of rows in
1268  // this block. if so see
1269  // whether we've just fallen
1270  // off the end of the whole
1271  // matrix
1272  if (local_row ==
1273  matrix->block(this->row_block, this->col_block).m())
1274  {
1275  local_row = 0;
1276  ++this->row_block;
1277  if (this->row_block == matrix->n_block_rows())
1278  {
1279  this->row_block = numbers::invalid_unsigned_int;
1280  this->col_block = numbers::invalid_unsigned_int;
1281  return;
1282  }
1283  }
1284 
1285  base_iterator =
1286  matrix->block(this->row_block, this->col_block).begin(local_row);
1287  }
1288  }
1289  }
1290 
1291 
1292  template <class BlockMatrixType>
1293  inline bool
1294  Accessor<BlockMatrixType, true>::operator==(const Accessor &a) const
1295  {
1296  if (matrix != a.matrix)
1297  return false;
1298 
1299  if (this->row_block == a.row_block && this->col_block == a.col_block)
1300  // end iterators do not necessarily
1301  // have to have the same
1302  // base_iterator representation, but
1303  // valid iterators have to
1304  return (((this->row_block == numbers::invalid_unsigned_int) &&
1305  (this->col_block == numbers::invalid_unsigned_int)) ||
1306  (base_iterator == a.base_iterator));
1307 
1308  return false;
1309  }
1310 
1311  //----------------------------------------------------------------------//
1312 
1313 
1314  template <class BlockMatrixType>
1315  inline Accessor<BlockMatrixType, false>::Accessor(BlockMatrixType *matrix,
1316  const size_type row,
1317  const size_type col)
1318  : matrix(matrix)
1319  , base_iterator(matrix->block(0, 0).begin())
1320  {
1321  (void)col;
1322  Assert(col == 0, ExcNotImplemented());
1323  // check if this is a regular row or
1324  // the end of the matrix
1325  if (row < matrix->m())
1326  {
1327  const std::pair<unsigned int, size_type> indices =
1328  matrix->row_block_indices.global_to_local(row);
1329 
1330  // find the first block that does
1331  // have an entry in this row
1332  for (size_type bc = 0; bc < matrix->n_block_cols(); ++bc)
1333  {
1334  base_iterator =
1335  matrix->block(indices.first, bc).begin(indices.second);
1336  if (base_iterator !=
1337  matrix->block(indices.first, bc).end(indices.second))
1338  {
1339  this->row_block = indices.first;
1340  this->col_block = bc;
1341  return;
1342  }
1343  }
1344 
1345  // hm, there is no block that has
1346  // an entry in this column. we need
1347  // to take the next entry then,
1348  // which may be the first entry of
1349  // the next row, or recursively the
1350  // next row, or so on
1351  *this = Accessor(matrix, row + 1, 0);
1352  }
1353  else
1354  {
1355  // we were asked to create the end
1356  // iterator for this matrix
1357  this->row_block = numbers::invalid_size_type;
1358  this->col_block = numbers::invalid_size_type;
1359  }
1360  }
1361 
1362 
1363  template <class BlockMatrixType>
1365  Accessor<BlockMatrixType, false>::row() const
1366  {
1367  Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1368 
1369  return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1370  base_iterator->row());
1371  }
1372 
1373 
1374  template <class BlockMatrixType>
1376  Accessor<BlockMatrixType, false>::column() const
1377  {
1378  Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1379 
1380  return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
1381  base_iterator->column());
1382  }
1383 
1384 
1385  template <class BlockMatrixType>
1386  inline typename Accessor<BlockMatrixType, false>::value_type
1387  Accessor<BlockMatrixType, false>::value() const
1388  {
1389  Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1390  Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1391 
1392  return base_iterator->value();
1393  }
1394 
1395 
1396 
1397  template <class BlockMatrixType>
1398  inline void
1399  Accessor<BlockMatrixType, false>::set_value(
1400  typename Accessor<BlockMatrixType, false>::value_type newval) const
1401  {
1402  Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1403  Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1404 
1405  base_iterator->value() = newval;
1406  }
1407 
1408 
1409 
1410  template <class BlockMatrixType>
1411  inline void
1413  {
1414  Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1415  Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1416 
1417  // Remember current row inside block
1418  size_type local_row = base_iterator->row();
1419 
1420  // Advance one element inside the
1421  // current block
1422  ++base_iterator;
1423 
1424  // while we hit the end of the row of a
1425  // block (which may happen multiple
1426  // times if rows inside a block are
1427  // empty), we have to jump to the next
1428  // block and take the
1429  while (base_iterator ==
1430  matrix->block(this->row_block, this->col_block).end(local_row))
1431  {
1432  // jump to next block in this block
1433  // row, if possible, otherwise go
1434  // to next row
1435  if (this->col_block < matrix->n_block_cols() - 1)
1436  {
1437  ++this->col_block;
1438  base_iterator =
1439  matrix->block(this->row_block, this->col_block).begin(local_row);
1440  }
1441  else
1442  {
1443  // jump back to next row in
1444  // first block column
1445  this->col_block = 0;
1446  ++local_row;
1447 
1448  // see if this has brought us
1449  // past the number of rows in
1450  // this block. if so see
1451  // whether we've just fallen
1452  // off the end of the whole
1453  // matrix
1454  if (local_row ==
1455  matrix->block(this->row_block, this->col_block).m())
1456  {
1457  local_row = 0;
1458  ++this->row_block;
1459  if (this->row_block == matrix->n_block_rows())
1460  {
1461  this->row_block = numbers::invalid_size_type;
1462  this->col_block = numbers::invalid_size_type;
1463  return;
1464  }
1465  }
1466 
1467  base_iterator =
1468  matrix->block(this->row_block, this->col_block).begin(local_row);
1469  }
1470  }
1471  }
1472 
1473 
1474 
1475  template <class BlockMatrixType>
1476  inline bool
1477  Accessor<BlockMatrixType, false>::operator==(const Accessor &a) const
1478  {
1479  if (matrix != a.matrix)
1480  return false;
1481 
1482  if (this->row_block == a.row_block && this->col_block == a.col_block)
1483  // end iterators do not necessarily
1484  // have to have the same
1485  // base_iterator representation, but
1486  // valid iterators have to
1487  return (((this->row_block == numbers::invalid_size_type) &&
1488  (this->col_block == numbers::invalid_size_type)) ||
1489  (base_iterator == a.base_iterator));
1490 
1491  return false;
1492  }
1493 } // namespace BlockMatrixIterators
1494 
1495 
1496 //---------------------------------------------------------------------------
1497 
1498 template <typename MatrixType>
1500 {
1501  try
1502  {
1503  clear();
1504  }
1505  catch (...)
1506  {}
1507 }
1508 
1509 
1510 template <class MatrixType>
1511 template <class BlockMatrixType>
1513 BlockMatrixBase<MatrixType>::copy_from(const BlockMatrixType &source)
1514 {
1515  for (unsigned int r = 0; r < n_block_rows(); ++r)
1516  for (unsigned int c = 0; c < n_block_cols(); ++c)
1517  block(r, c).copy_from(source.block(r, c));
1518 
1519  return *this;
1520 }
1521 
1522 
1523 template <class MatrixType>
1524 std::size_t
1526 {
1527  std::size_t mem =
1528  MemoryConsumption::memory_consumption(row_block_indices) +
1529  MemoryConsumption::memory_consumption(column_block_indices) +
1531  MemoryConsumption::memory_consumption(temporary_data.counter_within_block) +
1532  MemoryConsumption::memory_consumption(temporary_data.column_indices) +
1533  MemoryConsumption::memory_consumption(temporary_data.column_values) +
1534  sizeof(temporary_data.mutex);
1535 
1536  for (unsigned int r = 0; r < n_block_rows(); ++r)
1537  for (unsigned int c = 0; c < n_block_cols(); ++c)
1538  {
1539  MatrixType *p = this->sub_objects[r][c];
1541  }
1542 
1543  return mem;
1544 }
1545 
1546 
1547 
1548 template <class MatrixType>
1549 inline void
1551 {
1552  for (unsigned int r = 0; r < n_block_rows(); ++r)
1553  for (unsigned int c = 0; c < n_block_cols(); ++c)
1554  {
1555  MatrixType *p = this->sub_objects[r][c];
1556  this->sub_objects[r][c] = nullptr;
1557  delete p;
1558  }
1559  sub_objects.reinit(0, 0);
1560 
1561  // reset block indices to empty
1562  row_block_indices = column_block_indices = BlockIndices();
1563 }
1564 
1565 
1566 
1567 template <class MatrixType>
1569 BlockMatrixBase<MatrixType>::block(const unsigned int row,
1570  const unsigned int column)
1571 {
1572  AssertIndexRange(row, n_block_rows());
1573  AssertIndexRange(column, n_block_cols());
1574 
1575  return *sub_objects[row][column];
1576 }
1577 
1578 
1579 
1580 template <class MatrixType>
1581 inline const typename BlockMatrixBase<MatrixType>::BlockType &
1582 BlockMatrixBase<MatrixType>::block(const unsigned int row,
1583  const unsigned int column) const
1584 {
1585  AssertIndexRange(row, n_block_rows());
1586  AssertIndexRange(column, n_block_cols());
1587 
1588  return *sub_objects[row][column];
1589 }
1590 
1591 
1592 template <class MatrixType>
1595 {
1596  return row_block_indices.total_size();
1597 }
1598 
1599 
1600 
1601 template <class MatrixType>
1604 {
1605  return column_block_indices.total_size();
1606 }
1607 
1608 
1609 
1610 template <class MatrixType>
1611 inline unsigned int
1613 {
1614  return column_block_indices.size();
1615 }
1616 
1617 
1618 
1619 template <class MatrixType>
1620 inline unsigned int
1622 {
1623  return row_block_indices.size();
1624 }
1625 
1626 
1627 
1628 // Write the single set manually,
1629 // since the other function has a lot
1630 // of overhead in that case.
1631 template <class MatrixType>
1632 inline void
1634  const size_type j,
1635  const value_type value)
1636 {
1637  prepare_set_operation();
1638 
1639  AssertIsFinite(value);
1640 
1641  const std::pair<unsigned int, size_type>
1642  row_index = row_block_indices.global_to_local(i),
1643  col_index = column_block_indices.global_to_local(j);
1644  block(row_index.first, col_index.first)
1645  .set(row_index.second, col_index.second, value);
1646 }
1647 
1648 
1649 
1650 template <class MatrixType>
1651 template <typename number>
1652 inline void
1653 BlockMatrixBase<MatrixType>::set(const std::vector<size_type> &row_indices,
1654  const std::vector<size_type> &col_indices,
1655  const FullMatrix<number> & values,
1656  const bool elide_zero_values)
1657 {
1658  Assert(row_indices.size() == values.m(),
1659  ExcDimensionMismatch(row_indices.size(), values.m()));
1660  Assert(col_indices.size() == values.n(),
1661  ExcDimensionMismatch(col_indices.size(), values.n()));
1662 
1663  for (size_type i = 0; i < row_indices.size(); ++i)
1664  set(row_indices[i],
1665  col_indices.size(),
1666  col_indices.data(),
1667  &values(i, 0),
1668  elide_zero_values);
1669 }
1670 
1671 
1672 
1673 template <class MatrixType>
1674 template <typename number>
1675 inline void
1676 BlockMatrixBase<MatrixType>::set(const std::vector<size_type> &indices,
1677  const FullMatrix<number> & values,
1678  const bool elide_zero_values)
1679 {
1680  Assert(indices.size() == values.m(),
1681  ExcDimensionMismatch(indices.size(), values.m()));
1682  Assert(values.n() == values.m(), ExcNotQuadratic());
1683 
1684  for (size_type i = 0; i < indices.size(); ++i)
1685  set(indices[i],
1686  indices.size(),
1687  indices.data(),
1688  &values(i, 0),
1689  elide_zero_values);
1690 }
1691 
1692 
1693 
1694 template <class MatrixType>
1695 template <typename number>
1696 inline void
1698  const std::vector<size_type> &col_indices,
1699  const std::vector<number> & values,
1700  const bool elide_zero_values)
1701 {
1702  Assert(col_indices.size() == values.size(),
1703  ExcDimensionMismatch(col_indices.size(), values.size()));
1704 
1705  set(row,
1706  col_indices.size(),
1707  col_indices.data(),
1708  values.data(),
1709  elide_zero_values);
1710 }
1711 
1712 
1713 
1714 // This is a very messy function, since
1715 // we need to calculate to each position
1716 // the location in the global array.
1717 template <class MatrixType>
1718 template <typename number>
1719 inline void
1721  const size_type n_cols,
1722  const size_type *col_indices,
1723  const number * values,
1724  const bool elide_zero_values)
1725 {
1726  prepare_set_operation();
1727 
1728  // lock access to the temporary data structure to
1729  // allow multiple threads to call this function concurrently
1730  std::lock_guard<std::mutex> lock(temporary_data.mutex);
1731 
1732  // Resize scratch arrays
1733  if (temporary_data.column_indices.size() < this->n_block_cols())
1734  {
1735  temporary_data.column_indices.resize(this->n_block_cols());
1736  temporary_data.column_values.resize(this->n_block_cols());
1737  temporary_data.counter_within_block.resize(this->n_block_cols());
1738  }
1739 
1740  // Resize sub-arrays to n_cols. This
1741  // is a bit wasteful, but we resize
1742  // only a few times (then the maximum
1743  // row length won't increase that
1744  // much any more). At least we know
1745  // that all arrays are going to be of
1746  // the same size, so we can check
1747  // whether the size of one is large
1748  // enough before actually going
1749  // through all of them.
1750  if (temporary_data.column_indices[0].size() < n_cols)
1751  {
1752  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1753  {
1754  temporary_data.column_indices[i].resize(n_cols);
1755  temporary_data.column_values[i].resize(n_cols);
1756  }
1757  }
1758 
1759  // Reset the number of added elements
1760  // in each block to zero.
1761  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1762  temporary_data.counter_within_block[i] = 0;
1763 
1764  // Go through the column indices to
1765  // find out which portions of the
1766  // values should be set in which
1767  // block of the matrix. We need to
1768  // touch all the data, since we can't
1769  // be sure that the data of one block
1770  // is stored contiguously (in fact,
1771  // indices will be intermixed when it
1772  // comes from an element matrix).
1773  for (size_type j = 0; j < n_cols; ++j)
1774  {
1775  number value = values[j];
1776 
1777  if (value == number() && elide_zero_values == true)
1778  continue;
1779 
1780  const std::pair<unsigned int, size_type> col_index =
1781  this->column_block_indices.global_to_local(col_indices[j]);
1782 
1783  const size_type local_index =
1784  temporary_data.counter_within_block[col_index.first]++;
1785 
1786  temporary_data.column_indices[col_index.first][local_index] =
1787  col_index.second;
1788  temporary_data.column_values[col_index.first][local_index] = value;
1789  }
1790 
1791 # ifdef DEBUG
1792  // If in debug mode, do a check whether
1793  // the right length has been obtained.
1794  size_type length = 0;
1795  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1796  length += temporary_data.counter_within_block[i];
1797  Assert(length <= n_cols, ExcInternalError());
1798 # endif
1799 
1800  // Now we found out about where the
1801  // individual columns should start and
1802  // where we should start reading out
1803  // data. Now let's write the data into
1804  // the individual blocks!
1805  const std::pair<unsigned int, size_type> row_index =
1806  this->row_block_indices.global_to_local(row);
1807  for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
1808  {
1809  if (temporary_data.counter_within_block[block_col] == 0)
1810  continue;
1811 
1812  block(row_index.first, block_col)
1813  .set(row_index.second,
1814  temporary_data.counter_within_block[block_col],
1815  temporary_data.column_indices[block_col].data(),
1816  temporary_data.column_values[block_col].data(),
1817  false);
1818  }
1819 }
1820 
1821 
1822 
1823 template <class MatrixType>
1824 inline void
1826  const size_type j,
1827  const value_type value)
1828 {
1829  AssertIsFinite(value);
1830 
1831  prepare_add_operation();
1832 
1833  // save some cycles for zero additions, but
1834  // only if it is safe for the matrix we are
1835  // working with
1836  using MatrixTraits = typename MatrixType::Traits;
1837  if ((MatrixTraits::zero_addition_can_be_elided == true) &&
1838  (value == value_type()))
1839  return;
1840 
1841  const std::pair<unsigned int, size_type>
1842  row_index = row_block_indices.global_to_local(i),
1843  col_index = column_block_indices.global_to_local(j);
1844  block(row_index.first, col_index.first)
1845  .add(row_index.second, col_index.second, value);
1846 }
1847 
1848 
1849 
1850 template <class MatrixType>
1851 template <typename number>
1852 inline void
1853 BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &row_indices,
1854  const std::vector<size_type> &col_indices,
1855  const FullMatrix<number> & values,
1856  const bool elide_zero_values)
1857 {
1858  Assert(row_indices.size() == values.m(),
1859  ExcDimensionMismatch(row_indices.size(), values.m()));
1860  Assert(col_indices.size() == values.n(),
1861  ExcDimensionMismatch(col_indices.size(), values.n()));
1862 
1863  for (size_type i = 0; i < row_indices.size(); ++i)
1864  add(row_indices[i],
1865  col_indices.size(),
1866  col_indices.data(),
1867  &values(i, 0),
1868  elide_zero_values);
1869 }
1870 
1871 
1872 
1873 template <class MatrixType>
1874 template <typename number>
1875 inline void
1876 BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &indices,
1877  const FullMatrix<number> & values,
1878  const bool elide_zero_values)
1879 {
1880  Assert(indices.size() == values.m(),
1881  ExcDimensionMismatch(indices.size(), values.m()));
1882  Assert(values.n() == values.m(), ExcNotQuadratic());
1883 
1884  for (size_type i = 0; i < indices.size(); ++i)
1885  add(indices[i],
1886  indices.size(),
1887  indices.data(),
1888  &values(i, 0),
1889  elide_zero_values);
1890 }
1891 
1892 
1893 
1894 template <class MatrixType>
1895 template <typename number>
1896 inline void
1898  const std::vector<size_type> &col_indices,
1899  const std::vector<number> & values,
1900  const bool elide_zero_values)
1901 {
1902  Assert(col_indices.size() == values.size(),
1903  ExcDimensionMismatch(col_indices.size(), values.size()));
1904 
1905  add(row,
1906  col_indices.size(),
1907  col_indices.data(),
1908  values.data(),
1909  elide_zero_values);
1910 }
1911 
1912 
1913 
1914 // This is a very messy function, since
1915 // we need to calculate to each position
1916 // the location in the global array.
1917 template <class MatrixType>
1918 template <typename number>
1919 inline void
1921  const size_type n_cols,
1922  const size_type *col_indices,
1923  const number * values,
1924  const bool elide_zero_values,
1925  const bool col_indices_are_sorted)
1926 {
1927  prepare_add_operation();
1928 
1929  // TODO: Look over this to find out
1930  // whether we can do that more
1931  // efficiently.
1932  if (col_indices_are_sorted == true)
1933  {
1934 # ifdef DEBUG
1935  // check whether indices really are
1936  // sorted.
1937  size_type before = col_indices[0];
1938  for (size_type i = 1; i < n_cols; ++i)
1939  if (col_indices[i] <= before)
1940  Assert(false,
1941  ExcMessage("Flag col_indices_are_sorted is set, but "
1942  "indices appear to not be sorted.")) else before =
1943  col_indices[i];
1944 # endif
1945  const std::pair<unsigned int, size_type> row_index =
1946  this->row_block_indices.global_to_local(row);
1947 
1948  if (this->n_block_cols() > 1)
1949  {
1950  const size_type *first_block =
1951  Utilities::lower_bound(col_indices,
1952  col_indices + n_cols,
1953  this->column_block_indices.block_start(1));
1954 
1955  const size_type n_zero_block_indices = first_block - col_indices;
1956  block(row_index.first, 0)
1957  .add(row_index.second,
1958  n_zero_block_indices,
1959  col_indices,
1960  values,
1961  elide_zero_values,
1962  col_indices_are_sorted);
1963 
1964  if (n_zero_block_indices < n_cols)
1965  this->add(row,
1966  n_cols - n_zero_block_indices,
1967  first_block,
1968  values + n_zero_block_indices,
1969  elide_zero_values,
1970  false);
1971  }
1972  else
1973  {
1974  block(row_index.first, 0)
1975  .add(row_index.second,
1976  n_cols,
1977  col_indices,
1978  values,
1979  elide_zero_values,
1980  col_indices_are_sorted);
1981  }
1982 
1983  return;
1984  }
1985 
1986  // Lock scratch arrays, then resize them
1987  std::lock_guard<std::mutex> lock(temporary_data.mutex);
1988 
1989  if (temporary_data.column_indices.size() < this->n_block_cols())
1990  {
1991  temporary_data.column_indices.resize(this->n_block_cols());
1992  temporary_data.column_values.resize(this->n_block_cols());
1993  temporary_data.counter_within_block.resize(this->n_block_cols());
1994  }
1995 
1996  // Resize sub-arrays to n_cols. This
1997  // is a bit wasteful, but we resize
1998  // only a few times (then the maximum
1999  // row length won't increase that
2000  // much any more). At least we know
2001  // that all arrays are going to be of
2002  // the same size, so we can check
2003  // whether the size of one is large
2004  // enough before actually going
2005  // through all of them.
2006  if (temporary_data.column_indices[0].size() < n_cols)
2007  {
2008  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2009  {
2010  temporary_data.column_indices[i].resize(n_cols);
2011  temporary_data.column_values[i].resize(n_cols);
2012  }
2013  }
2014 
2015  // Reset the number of added elements
2016  // in each block to zero.
2017  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2018  temporary_data.counter_within_block[i] = 0;
2019 
2020  // Go through the column indices to
2021  // find out which portions of the
2022  // values should be written into
2023  // which block of the matrix. We need
2024  // to touch all the data, since we
2025  // can't be sure that the data of one
2026  // block is stored contiguously (in
2027  // fact, data will be intermixed when
2028  // it comes from an element matrix).
2029  for (size_type j = 0; j < n_cols; ++j)
2030  {
2031  number value = values[j];
2032 
2033  if (value == number() && elide_zero_values == true)
2034  continue;
2035 
2036  const std::pair<unsigned int, size_type> col_index =
2037  this->column_block_indices.global_to_local(col_indices[j]);
2038 
2039  const size_type local_index =
2040  temporary_data.counter_within_block[col_index.first]++;
2041 
2042  temporary_data.column_indices[col_index.first][local_index] =
2043  col_index.second;
2044  temporary_data.column_values[col_index.first][local_index] = value;
2045  }
2046 
2047 # ifdef DEBUG
2048  // If in debug mode, do a check whether
2049  // the right length has been obtained.
2050  size_type length = 0;
2051  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2052  length += temporary_data.counter_within_block[i];
2053  Assert(length <= n_cols, ExcInternalError());
2054 # endif
2055 
2056  // Now we found out about where the
2057  // individual columns should start and
2058  // where we should start reading out
2059  // data. Now let's write the data into
2060  // the individual blocks!
2061  const std::pair<unsigned int, size_type> row_index =
2062  this->row_block_indices.global_to_local(row);
2063  for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
2064  {
2065  if (temporary_data.counter_within_block[block_col] == 0)
2066  continue;
2067 
2068  block(row_index.first, block_col)
2069  .add(row_index.second,
2070  temporary_data.counter_within_block[block_col],
2071  temporary_data.column_indices[block_col].data(),
2072  temporary_data.column_values[block_col].data(),
2073  false,
2074  col_indices_are_sorted);
2075  }
2076 }
2077 
2078 
2079 
2080 template <class MatrixType>
2081 inline void
2082 BlockMatrixBase<MatrixType>::add(const value_type factor,
2084 {
2085  AssertIsFinite(factor);
2086 
2087  prepare_add_operation();
2088 
2089  // save some cycles for zero additions, but
2090  // only if it is safe for the matrix we are
2091  // working with
2092  using MatrixTraits = typename MatrixType::Traits;
2093  if ((MatrixTraits::zero_addition_can_be_elided == true) && (factor == 0))
2094  return;
2095 
2096  for (unsigned int row = 0; row < n_block_rows(); ++row)
2097  for (unsigned int col = 0; col < n_block_cols(); ++col)
2098  // This function should throw if the sparsity
2099  // patterns of the two blocks differ
2100  block(row, col).add(factor, matrix.block(row, col));
2101 }
2102 
2103 
2104 
2105 template <class MatrixType>
2108  const size_type j) const
2109 {
2110  const std::pair<unsigned int, size_type>
2111  row_index = row_block_indices.global_to_local(i),
2112  col_index = column_block_indices.global_to_local(j);
2113  return block(row_index.first, col_index.first)(row_index.second,
2114  col_index.second);
2115 }
2116 
2117 
2118 
2119 template <class MatrixType>
2121 BlockMatrixBase<MatrixType>::el(const size_type i, const size_type j) const
2122 {
2123  const std::pair<unsigned int, size_type>
2124  row_index = row_block_indices.global_to_local(i),
2125  col_index = column_block_indices.global_to_local(j);
2126  return block(row_index.first, col_index.first)
2127  .el(row_index.second, col_index.second);
2128 }
2129 
2130 
2131 
2132 template <class MatrixType>
2135 {
2136  Assert(n_block_rows() == n_block_cols(), ExcNotQuadratic());
2137 
2138  const std::pair<unsigned int, size_type> index =
2139  row_block_indices.global_to_local(i);
2140  return block(index.first, index.first).diag_element(index.second);
2141 }
2142 
2143 
2144 
2145 template <class MatrixType>
2146 inline void
2148 {
2149  for (unsigned int r = 0; r < n_block_rows(); ++r)
2150  for (unsigned int c = 0; c < n_block_cols(); ++c)
2151  block(r, c).compress(operation);
2152 }
2153 
2154 
2155 
2156 template <class MatrixType>
2158 BlockMatrixBase<MatrixType>::operator*=(const value_type factor)
2159 {
2160  Assert(n_block_cols() != 0, ExcNotInitialized());
2161  Assert(n_block_rows() != 0, ExcNotInitialized());
2162 
2163  for (unsigned int r = 0; r < n_block_rows(); ++r)
2164  for (unsigned int c = 0; c < n_block_cols(); ++c)
2165  block(r, c) *= factor;
2166 
2167  return *this;
2168 }
2169 
2170 
2171 
2172 template <class MatrixType>
2174 BlockMatrixBase<MatrixType>::operator/=(const value_type factor)
2175 {
2176  Assert(n_block_cols() != 0, ExcNotInitialized());
2177  Assert(n_block_rows() != 0, ExcNotInitialized());
2178  Assert(factor != 0, ExcDivideByZero());
2179 
2180  const value_type factor_inv = 1. / factor;
2181 
2182  for (unsigned int r = 0; r < n_block_rows(); ++r)
2183  for (unsigned int c = 0; c < n_block_cols(); ++c)
2184  block(r, c) *= factor_inv;
2185 
2186  return *this;
2187 }
2188 
2189 
2190 
2191 template <class MatrixType>
2192 const BlockIndices &
2194 {
2195  return this->row_block_indices;
2196 }
2197 
2198 
2199 
2200 template <class MatrixType>
2201 const BlockIndices &
2203 {
2204  return this->column_block_indices;
2205 }
2206 
2207 
2208 
2209 template <class MatrixType>
2210 template <class BlockVectorType>
2211 void
2213  const BlockVectorType &src) const
2214 {
2215  Assert(dst.n_blocks() == n_block_rows(),
2216  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2217  Assert(src.n_blocks() == n_block_cols(),
2218  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2219 
2220  for (size_type row = 0; row < n_block_rows(); ++row)
2221  {
2222  block(row, 0).vmult(dst.block(row), src.block(0));
2223  for (size_type col = 1; col < n_block_cols(); ++col)
2224  block(row, col).vmult_add(dst.block(row), src.block(col));
2225  };
2226 }
2227 
2228 
2229 
2230 template <class MatrixType>
2231 template <class BlockVectorType, class VectorType>
2232 void
2234  VectorType & dst,
2235  const BlockVectorType &src) const
2236 {
2237  Assert(n_block_rows() == 1, ExcDimensionMismatch(1, n_block_rows()));
2238  Assert(src.n_blocks() == n_block_cols(),
2239  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2240 
2241  block(0, 0).vmult(dst, src.block(0));
2242  for (size_type col = 1; col < n_block_cols(); ++col)
2243  block(0, col).vmult_add(dst, src.block(col));
2244 }
2245 
2246 
2247 
2248 template <class MatrixType>
2249 template <class BlockVectorType, class VectorType>
2250 void
2252  const VectorType &src) const
2253 {
2254  Assert(dst.n_blocks() == n_block_rows(),
2255  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2256  Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2257 
2258  for (size_type row = 0; row < n_block_rows(); ++row)
2259  block(row, 0).vmult(dst.block(row), src);
2260 }
2261 
2262 
2263 
2264 template <class MatrixType>
2265 template <class VectorType>
2266 void
2268  VectorType & dst,
2269  const VectorType &src) const
2270 {
2271  Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));
2272  Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2273 
2274  block(0, 0).vmult(dst, src);
2275 }
2276 
2277 
2278 
2279 template <class MatrixType>
2280 template <class BlockVectorType>
2281 void
2282 BlockMatrixBase<MatrixType>::vmult_add(BlockVectorType & dst,
2283  const BlockVectorType &src) const
2284 {
2285  Assert(dst.n_blocks() == n_block_rows(),
2286  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2287  Assert(src.n_blocks() == n_block_cols(),
2288  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2289 
2290  for (unsigned int row = 0; row < n_block_rows(); ++row)
2291  for (unsigned int col = 0; col < n_block_cols(); ++col)
2292  block(row, col).vmult_add(dst.block(row), src.block(col));
2293 }
2294 
2295 
2296 
2297 template <class MatrixType>
2298 template <class BlockVectorType>
2299 void
2301  BlockVectorType & dst,
2302  const BlockVectorType &src) const
2303 {
2304  Assert(dst.n_blocks() == n_block_cols(),
2305  ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2306  Assert(src.n_blocks() == n_block_rows(),
2307  ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2308 
2309  dst = 0.;
2310 
2311  for (unsigned int row = 0; row < n_block_rows(); ++row)
2312  {
2313  for (unsigned int col = 0; col < n_block_cols(); ++col)
2314  block(row, col).Tvmult_add(dst.block(col), src.block(row));
2315  };
2316 }
2317 
2318 
2319 
2320 template <class MatrixType>
2321 template <class BlockVectorType, class VectorType>
2322 void
2324  const VectorType &src) const
2325 {
2326  Assert(dst.n_blocks() == n_block_cols(),
2327  ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2328  Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));
2329 
2330  dst = 0.;
2331 
2332  for (unsigned int col = 0; col < n_block_cols(); ++col)
2333  block(0, col).Tvmult_add(dst.block(col), src);
2334 }
2335 
2336 
2337 
2338 template <class MatrixType>
2339 template <class BlockVectorType, class VectorType>
2340 void
2342  VectorType & dst,
2343  const BlockVectorType &src) const
2344 {
2345  Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2346  Assert(src.n_blocks() == n_block_rows(),
2347  ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2348 
2349  block(0, 0).Tvmult(dst, src.block(0));
2350 
2351  for (size_type row = 1; row < n_block_rows(); ++row)
2352  block(row, 0).Tvmult_add(dst, src.block(row));
2353 }
2354 
2355 
2356 
2357 template <class MatrixType>
2358 template <class VectorType>
2359 void
2361  VectorType & dst,
2362  const VectorType &src) const
2363 {
2364  Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2365  Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));
2366 
2367  block(0, 0).Tvmult(dst, src);
2368 }
2369 
2370 
2371 
2372 template <class MatrixType>
2373 template <class BlockVectorType>
2374 void
2375 BlockMatrixBase<MatrixType>::Tvmult_add(BlockVectorType & dst,
2376  const BlockVectorType &src) const
2377 {
2378  Assert(dst.n_blocks() == n_block_cols(),
2379  ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2380  Assert(src.n_blocks() == n_block_rows(),
2381  ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2382 
2383  for (unsigned int row = 0; row < n_block_rows(); ++row)
2384  for (unsigned int col = 0; col < n_block_cols(); ++col)
2385  block(row, col).Tvmult_add(dst.block(col), src.block(row));
2386 }
2387 
2388 
2389 
2390 template <class MatrixType>
2391 template <class BlockVectorType>
2393 BlockMatrixBase<MatrixType>::matrix_norm_square(const BlockVectorType &v) const
2394 {
2395  Assert(n_block_rows() == n_block_cols(), ExcNotQuadratic());
2396  Assert(v.n_blocks() == n_block_rows(),
2397  ExcDimensionMismatch(v.n_blocks(), n_block_rows()));
2398 
2399  value_type norm_sqr = 0;
2400  for (unsigned int row = 0; row < n_block_rows(); ++row)
2401  for (unsigned int col = 0; col < n_block_cols(); ++col)
2402  if (row == col)
2403  norm_sqr += block(row, col).matrix_norm_square(v.block(row));
2404  else
2405  norm_sqr +=
2406  block(row, col).matrix_scalar_product(v.block(row), v.block(col));
2407  return norm_sqr;
2408 }
2409 
2410 
2411 
2412 template <class MatrixType>
2415 {
2416  value_type norm_sqr = 0;
2417 
2418  // For each block, get the Frobenius norm, and add the square to the
2419  // accumulator for the full matrix
2420  for (unsigned int row = 0; row < n_block_rows(); ++row)
2421  {
2422  for (unsigned int col = 0; col < n_block_cols(); ++col)
2423  {
2424  const value_type block_norm = block(row, col).frobenius_norm();
2425  norm_sqr += block_norm * block_norm;
2426  }
2427  }
2428 
2429  return std::sqrt(norm_sqr);
2430 }
2431 
2432 
2433 
2434 template <class MatrixType>
2435 template <class BlockVectorType>
2438  const BlockVectorType &u,
2439  const BlockVectorType &v) const
2440 {
2441  Assert(u.n_blocks() == n_block_rows(),
2442  ExcDimensionMismatch(u.n_blocks(), n_block_rows()));
2443  Assert(v.n_blocks() == n_block_cols(),
2444  ExcDimensionMismatch(v.n_blocks(), n_block_cols()));
2445 
2446  value_type result = 0;
2447  for (unsigned int row = 0; row < n_block_rows(); ++row)
2448  for (unsigned int col = 0; col < n_block_cols(); ++col)
2449  result +=
2450  block(row, col).matrix_scalar_product(u.block(row), v.block(col));
2451  return result;
2452 }
2453 
2454 
2455 
2456 template <class MatrixType>
2457 template <class BlockVectorType>
2459 BlockMatrixBase<MatrixType>::residual(BlockVectorType & dst,
2460  const BlockVectorType &x,
2461  const BlockVectorType &b) const
2462 {
2463  Assert(dst.n_blocks() == n_block_rows(),
2464  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2465  Assert(b.n_blocks() == n_block_rows(),
2466  ExcDimensionMismatch(b.n_blocks(), n_block_rows()));
2467  Assert(x.n_blocks() == n_block_cols(),
2468  ExcDimensionMismatch(x.n_blocks(), n_block_cols()));
2469  // in block notation, the residual is
2470  // r_i = b_i - \sum_j A_ij x_j.
2471  // this can be written as
2472  // r_i = b_i - A_i0 x_0 - \sum_{j>0} A_ij x_j.
2473  //
2474  // for the first two terms, we can
2475  // call the residual function of
2476  // A_i0. for the other terms, we
2477  // use vmult_add. however, we want
2478  // to subtract, so in order to
2479  // avoid a temporary vector, we
2480  // perform a sign change of the
2481  // first two term before, and after
2482  // adding up
2483  for (unsigned int row = 0; row < n_block_rows(); ++row)
2484  {
2485  block(row, 0).residual(dst.block(row), x.block(0), b.block(row));
2486 
2487  for (size_type i = 0; i < dst.block(row).size(); ++i)
2488  dst.block(row)(i) = -dst.block(row)(i);
2489 
2490  for (unsigned int col = 1; col < n_block_cols(); ++col)
2491  block(row, col).vmult_add(dst.block(row), x.block(col));
2492 
2493  for (size_type i = 0; i < dst.block(row).size(); ++i)
2494  dst.block(row)(i) = -dst.block(row)(i);
2495  };
2496 
2497  value_type res = 0;
2498  for (size_type row = 0; row < n_block_rows(); ++row)
2499  res += dst.block(row).norm_sqr();
2500  return std::sqrt(res);
2501 }
2502 
2503 
2504 
2505 template <class MatrixType>
2506 inline void
2507 BlockMatrixBase<MatrixType>::print(std::ostream &out,
2508  const bool alternative_output) const
2509 {
2510  for (unsigned int row = 0; row < n_block_rows(); ++row)
2511  for (unsigned int col = 0; col < n_block_cols(); ++col)
2512  {
2513  if (!alternative_output)
2514  out << "Block (" << row << ", " << col << ')' << std::endl;
2515 
2516  block(row, col).print(out, alternative_output);
2517  }
2518 }
2519 
2520 
2521 
2522 template <class MatrixType>
2525 {
2526  return const_iterator(this, 0);
2527 }
2528 
2529 
2530 
2531 template <class MatrixType>
2534 {
2535  return const_iterator(this, m());
2536 }
2537 
2538 
2539 
2540 template <class MatrixType>
2543 {
2544  AssertIndexRange(r, m());
2545  return const_iterator(this, r);
2546 }
2547 
2548 
2549 
2550 template <class MatrixType>
2553 {
2554  AssertIndexRange(r, m());
2555  return const_iterator(this, r + 1);
2556 }
2557 
2558 
2559 
2560 template <class MatrixType>
2563 {
2564  return iterator(this, 0);
2565 }
2566 
2567 
2568 
2569 template <class MatrixType>
2572 {
2573  return iterator(this, m());
2574 }
2575 
2576 
2577 
2578 template <class MatrixType>
2581 {
2582  AssertIndexRange(r, m());
2583  return iterator(this, r);
2584 }
2585 
2586 
2587 
2588 template <class MatrixType>
2591 {
2592  AssertIndexRange(r, m());
2593  return iterator(this, r + 1);
2594 }
2595 
2596 
2597 
2598 template <class MatrixType>
2599 void
2601 {
2602  std::vector<size_type> row_sizes(this->n_block_rows());
2603  std::vector<size_type> col_sizes(this->n_block_cols());
2604 
2605  // first find out the row sizes
2606  // from the first block column
2607  for (unsigned int r = 0; r < this->n_block_rows(); ++r)
2608  row_sizes[r] = sub_objects[r][0]->m();
2609  // then check that the following
2610  // block columns have the same
2611  // sizes
2612  for (unsigned int c = 1; c < this->n_block_cols(); ++c)
2613  for (unsigned int r = 0; r < this->n_block_rows(); ++r)
2614  Assert(row_sizes[r] == sub_objects[r][c]->m(),
2615  ExcIncompatibleRowNumbers(r, 0, r, c));
2616 
2617  // finally initialize the row
2618  // indices with this array
2619  this->row_block_indices.reinit(row_sizes);
2620 
2621 
2622  // then do the same with the columns
2623  for (unsigned int c = 0; c < this->n_block_cols(); ++c)
2624  col_sizes[c] = sub_objects[0][c]->n();
2625  for (unsigned int r = 1; r < this->n_block_rows(); ++r)
2626  for (unsigned int c = 0; c < this->n_block_cols(); ++c)
2627  Assert(col_sizes[c] == sub_objects[r][c]->n(),
2628  ExcIncompatibleRowNumbers(0, c, r, c));
2629 
2630  // finally initialize the row
2631  // indices with this array
2632  this->column_block_indices.reinit(col_sizes);
2633 }
2634 
2635 
2636 
2637 template <class MatrixType>
2638 void
2640 {
2641  for (unsigned int row = 0; row < n_block_rows(); ++row)
2642  for (unsigned int col = 0; col < n_block_cols(); ++col)
2643  block(row, col).prepare_add();
2644 }
2645 
2646 
2647 
2648 template <class MatrixType>
2649 void
2651 {
2652  for (unsigned int row = 0; row < n_block_rows(); ++row)
2653  for (unsigned int col = 0; col < n_block_cols(); ++col)
2654  block(row, col).prepare_set();
2655 }
2656 
2657 #endif // DOXYGEN
2658 
2659 
2661 
2662 #endif // dealii_block_matrix_base_h
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
void Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
void vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
void add(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number > &values, const bool elide_zero_values=true)
void add(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=true)
value_type & reference
void add(const size_type row, const size_type n_cols, const size_type *col_indices, const number *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
void print(std::ostream &out, const bool alternative_output=false) const
const_iterator end() const
BlockIndices column_block_indices
value_type matrix_norm_square(const BlockVectorType &v) const
void Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
BlockType & block(const unsigned int row, const unsigned int column)
const value_type & const_reference
void Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
BlockMatrixBase & operator*=(const value_type factor)
unsigned int n_block_rows() const
value_type operator()(const size_type i, const size_type j) const
value_type el(const size_type i, const size_type j) const
const_iterator begin() const
real_type frobenius_norm() const
void vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
void compress(VectorOperation::values operation)
void vmult_add(BlockVectorType &dst, const BlockVectorType &src) const
BlockMatrixBase & operator/=(const value_type factor)
types::global_dof_index size_type
Table< 2, SmartPointer< BlockType, BlockMatrixBase< MatrixType > > > sub_objects
void Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const
typename BlockType::value_type value_type
void collect_sizes()
void prepare_add_operation()
void vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
value_type residual(BlockVectorType &dst, const BlockVectorType &x, const BlockVectorType &b) const
BlockMatrixBase()=default
void set(const size_type row, const size_type n_cols, const size_type *col_indices, const number *values, const bool elide_zero_values=false)
size_type n() const
void prepare_set_operation()
unsigned int n_block_cols() const
const value_type * const_pointer
iterator begin()
TemporaryData temporary_data
iterator end()
void set(const size_type i, const size_type j, const value_type value)
std::size_t memory_consumption() const
void add(const value_type factor, const BlockMatrixBase< MatrixType > &matrix)
typename numbers::NumberTraits< value_type >::real_type real_type
BlockIndices row_block_indices
const BlockIndices & get_row_indices() const
~BlockMatrixBase() override
size_type m() const
value_type matrix_scalar_product(const BlockVectorType &u, const BlockVectorType &v) const
BlockMatrixBase & copy_from(const BlockMatrixType &source)
void set(const std::vector< size_type > &indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=false)
void add(const std::vector< size_type > &indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=true)
void add(const size_type i, const size_type j, const value_type value)
void set(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number > &values, const bool elide_zero_values=false)
value_type diag_element(const size_type i) const
void vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
value_type * pointer
iterator begin(const size_type r)
const_iterator begin(const size_type r) const
const BlockType & block(const unsigned int row, const unsigned int column) const
void Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
iterator end(const size_type r)
const_iterator end(const size_type r) const
void set(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=false)
const BlockIndices & get_column_indices() const
unsigned int block_row() const
types::global_dof_index size_type
unsigned int block_column() const
typename BlockMatrixType::value_type value_type
Accessor(BlockMatrixType *m, const size_type row, const size_type col)
Accessor(const Accessor< BlockMatrixType, false > &)
BlockMatrixType::BlockType::const_iterator base_iterator
Accessor(const BlockMatrixType *m, const size_type row, const size_type col)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcInternalError()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:579
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcIteratorPastEnd()
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
Definition: exceptions.h:1854
static ::ExceptionBase & ExcIncompatibleColNumbers(int arg1, int arg2, int arg3, int arg4)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
static ::ExceptionBase & ExcIncompatibleRowNumbers(int arg1, int arg2, int arg3, int arg4)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNotQuadratic()
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
Definition: cuda_kernels.h:45
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:999
static const unsigned int invalid_unsigned_int
Definition: types.h:213
const types::global_dof_index invalid_size_type
Definition: types.h:222
unsigned int global_dof_index
Definition: types.h:82
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
std::vector< std::vector< size_type > > column_indices
std::vector< std::vector< value_type > > column_values
TemporaryData & operator=(const TemporaryData &)
std::vector< size_type > counter_within_block
void advance(std::tuple< I1, I2 > &t, const unsigned int n)