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