Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
block_matrix_base.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2019 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 
22 #include <deal.II/base/memory_consumption.h>
23 #include <deal.II/base/smartpointer.h>
24 #include <deal.II/base/table.h>
25 #include <deal.II/base/thread_management.h>
26 #include <deal.II/base/utilities.h>
27 
28 #include <deal.II/lac/block_indices.h>
29 #include <deal.II/lac/exceptions.h>
30 #include <deal.II/lac/full_matrix.h>
31 #include <deal.II/lac/matrix_iterator.h>
32 #include <deal.II/lac/vector.h>
33 #include <deal.II/lac/vector_operation.h>
34 
35 #include <cmath>
36 
37 DEAL_II_NAMESPACE_OPEN
38 
39 
40 // Forward declaration
41 #ifndef DOXYGEN
42 template <typename>
43 class MatrixIterator;
44 #endif
45 
46 
57 {
62  template <class BlockMatrixType>
64  {
65  public:
70 
74  using value_type = typename BlockMatrixType::value_type;
75 
79  AccessorBase();
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
189  advance();
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
275  advance();
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 
351 template <typename MatrixType>
352 class BlockMatrixBase : public Subscriptor
353 {
354 public:
358  using BlockType = MatrixType;
359 
365  using real_type = typename numbers::NumberTraits<value_type>::real_type;
366  using pointer = value_type *;
367  using const_pointer = const value_type *;
368  using reference = value_type &;
369  using const_reference = const value_type &;
370  using size_type = types::global_dof_index;
371 
372  using iterator =
374 
375  using const_iterator =
377 
378 
382  BlockMatrixBase() = default;
383 
387  ~BlockMatrixBase() override;
388 
405  template <class BlockMatrixType>
407  copy_from(const BlockMatrixType &source);
408 
412  BlockType &
413  block(const unsigned int row, const unsigned int column);
414 
415 
420  const BlockType &
421  block(const unsigned int row, const unsigned int column) const;
422 
427  size_type
428  m() const;
429 
434  size_type
435  n() const;
436 
437 
442  unsigned int
443  n_block_rows() const;
444 
449  unsigned int
450  n_block_cols() const;
451 
457  void
458  set(const size_type i, const size_type j, const value_type value);
459 
475  template <typename number>
476  void
477  set(const std::vector<size_type> &indices,
478  const FullMatrix<number> & full_matrix,
479  const bool elide_zero_values = false);
480 
486  template <typename number>
487  void
488  set(const std::vector<size_type> &row_indices,
489  const std::vector<size_type> &col_indices,
490  const FullMatrix<number> & full_matrix,
491  const bool elide_zero_values = false);
492 
503  template <typename number>
504  void
505  set(const size_type row,
506  const std::vector<size_type> &col_indices,
507  const std::vector<number> & values,
508  const bool elide_zero_values = false);
509 
519  template <typename number>
520  void
521  set(const size_type row,
522  const size_type n_cols,
523  const size_type *col_indices,
524  const number * values,
525  const bool elide_zero_values = false);
526 
532  void
533  add(const size_type i, const size_type j, const value_type value);
534 
549  template <typename number>
550  void
551  add(const std::vector<size_type> &indices,
552  const FullMatrix<number> & full_matrix,
553  const bool elide_zero_values = true);
554 
560  template <typename number>
561  void
562  add(const std::vector<size_type> &row_indices,
563  const std::vector<size_type> &col_indices,
564  const FullMatrix<number> & full_matrix,
565  const bool elide_zero_values = true);
566 
576  template <typename number>
577  void
578  add(const size_type row,
579  const std::vector<size_type> &col_indices,
580  const std::vector<number> & values,
581  const bool elide_zero_values = true);
582 
592  template <typename number>
593  void
594  add(const size_type row,
595  const size_type n_cols,
596  const size_type *col_indices,
597  const number * values,
598  const bool elide_zero_values = true,
599  const bool col_indices_are_sorted = false);
600 
612  void
613  add(const value_type factor, const BlockMatrixBase<MatrixType> &matrix);
614 
621  value_type
622  operator()(const size_type i, const size_type j) const;
623 
632  value_type
633  el(const size_type i, const size_type j) const;
634 
645  value_type
646  diag_element(const size_type i) const;
647 
656  void
657  compress(::VectorOperation::values operation);
658 
663  operator*=(const value_type factor);
664 
669  operator/=(const value_type factor);
670 
675  template <class BlockVectorType>
676  void
677  vmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
678 
684  template <class BlockVectorType>
685  void
686  Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
687 
700  template <class BlockVectorType>
701  value_type
702  matrix_norm_square(const BlockVectorType &v) const;
703 
708  real_type
709  frobenius_norm() const;
710 
714  template <class BlockVectorType>
715  value_type
716  matrix_scalar_product(const BlockVectorType &u,
717  const BlockVectorType &v) const;
718 
722  template <class BlockVectorType>
723  value_type
724  residual(BlockVectorType & dst,
725  const BlockVectorType &x,
726  const BlockVectorType &b) const;
727 
734  void
735  print(std::ostream &out, const bool alternative_output = false) const;
736 
740  iterator
741  begin();
742 
746  iterator
747  end();
748 
752  iterator
753  begin(const size_type r);
754 
758  iterator
759  end(const size_type r);
764  begin() const;
765 
770  end() const;
771 
776  begin(const size_type r) const;
777 
782  end(const size_type r) const;
783 
787  const BlockIndices &
788  get_row_indices() const;
789 
793  const BlockIndices &
794  get_column_indices() const;
795 
801  std::size_t
802  memory_consumption() const;
803 
812  DeclException4(ExcIncompatibleRowNumbers,
813  int,
814  int,
815  int,
816  int,
817  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
818  << ',' << arg4 << "] have differing row numbers.");
822  DeclException4(ExcIncompatibleColNumbers,
823  int,
824  int,
825  int,
826  int,
827  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
828  << ',' << arg4 << "] have differing column numbers.");
830 protected:
843  void
844  clear();
845 
850  BlockIndices column_block_indices;
851 
856 
875  void
876  collect_sizes();
877 
888  template <class BlockVectorType>
889  void
890  vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const;
891 
902  template <class BlockVectorType, class VectorType>
903  void
904  vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const;
905 
916  template <class BlockVectorType, class VectorType>
917  void
918  vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const;
919 
930  template <class VectorType>
931  void
932  vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const;
933 
945  template <class BlockVectorType>
946  void
947  Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const;
948 
959  template <class BlockVectorType, class VectorType>
960  void
961  Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const;
962 
973  template <class BlockVectorType, class VectorType>
974  void
975  Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const;
976 
987  template <class VectorType>
988  void
989  Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const;
990 
991 
992 protected:
999  void
1000  prepare_add_operation();
1001 
1006  void
1007  prepare_set_operation();
1008 
1009 
1010 private:
1020  {
1025  std::vector<size_type> counter_within_block;
1026 
1031  std::vector<std::vector<size_type>> column_indices;
1032 
1037  std::vector<std::vector<value_type>> column_values;
1038 
1044 
1054  TemporaryData &
1056  {
1057  return *this;
1058  }
1059  };
1060 
1067  TemporaryData temporary_data;
1068 
1069  // Make the iterator class a friend. We have to work around a compiler bug
1070  // here again.
1071  template <typename, bool>
1072  friend class BlockMatrixIterators::Accessor;
1073 
1074  template <typename>
1075  friend class MatrixIterator;
1076 };
1077 
1078 
1081 #ifndef DOXYGEN
1082 /* ------------------------- Template functions ---------------------- */
1083 
1084 
1085 namespace BlockMatrixIterators
1086 {
1087  template <class BlockMatrixType>
1089  : row_block(0)
1090  , col_block(0)
1091  {}
1092 
1093 
1094  template <class BlockMatrixType>
1095  inline unsigned int
1096  AccessorBase<BlockMatrixType>::block_row() const
1097  {
1099 
1100  return row_block;
1101  }
1102 
1103 
1104  template <class BlockMatrixType>
1105  inline unsigned int
1106  AccessorBase<BlockMatrixType>::block_column() const
1107  {
1109 
1110  return col_block;
1111  }
1112 
1113 
1114  template <class BlockMatrixType>
1115  inline Accessor<BlockMatrixType, true>::Accessor(
1116  const BlockMatrixType *matrix,
1117  const size_type row,
1118  const size_type col)
1119  : matrix(matrix)
1120  , base_iterator(matrix->block(0, 0).begin())
1121  {
1122  (void)col;
1123  Assert(col == 0, ExcNotImplemented());
1124 
1125  // check if this is a regular row or
1126  // the end of the matrix
1127  if (row < matrix->m())
1128  {
1129  const std::pair<unsigned int, size_type> indices =
1130  matrix->row_block_indices.global_to_local(row);
1131 
1132  // find the first block that does
1133  // have an entry in this row
1134  for (unsigned int bc = 0; bc < matrix->n_block_cols(); ++bc)
1135  {
1136  base_iterator =
1137  matrix->block(indices.first, bc).begin(indices.second);
1138  if (base_iterator !=
1139  matrix->block(indices.first, bc).end(indices.second))
1140  {
1141  this->row_block = indices.first;
1142  this->col_block = bc;
1143  return;
1144  }
1145  }
1146 
1147  // hm, there is no block that has
1148  // an entry in this column. we need
1149  // to take the next entry then,
1150  // which may be the first entry of
1151  // the next row, or recursively the
1152  // next row, or so on
1153  *this = Accessor(matrix, row + 1, 0);
1154  }
1155  else
1156  {
1157  // we were asked to create the end
1158  // iterator for this matrix
1159  this->row_block = numbers::invalid_unsigned_int;
1160  this->col_block = numbers::invalid_unsigned_int;
1161  }
1162  }
1163 
1164 
1165  // template <class BlockMatrixType>
1166  // inline
1167  // Accessor<BlockMatrixType, true>::Accessor (const
1168  // Accessor<BlockMatrixType, true>& other)
1169  // :
1170  // matrix(other.matrix),
1171  // base_iterator(other.base_iterator)
1172  // {
1173  // this->row_block = other.row_block;
1174  // this->col_block = other.col_block;
1175  // }
1176 
1177 
1178  template <class BlockMatrixType>
1179  inline Accessor<BlockMatrixType, true>::Accessor(
1180  const Accessor<BlockMatrixType, false> &other)
1181  : matrix(other.matrix)
1182  , base_iterator(other.base_iterator)
1183  {
1184  this->row_block = other.row_block;
1185  this->col_block = other.col_block;
1186  }
1187 
1188 
1189  template <class BlockMatrixType>
1190  inline typename Accessor<BlockMatrixType, true>::size_type
1191  Accessor<BlockMatrixType, true>::row() const
1192  {
1193  Assert(this->row_block != numbers::invalid_unsigned_int,
1194  ExcIteratorPastEnd());
1195 
1196  return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1197  base_iterator->row());
1198  }
1199 
1200 
1201  template <class BlockMatrixType>
1202  inline typename Accessor<BlockMatrixType, true>::size_type
1203  Accessor<BlockMatrixType, true>::column() const
1204  {
1205  Assert(this->col_block != numbers::invalid_unsigned_int,
1206  ExcIteratorPastEnd());
1207 
1208  return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
1209  base_iterator->column());
1210  }
1211 
1212 
1213  template <class BlockMatrixType>
1214  inline typename Accessor<BlockMatrixType, true>::value_type
1215  Accessor<BlockMatrixType, true>::value() const
1216  {
1217  Assert(this->row_block != numbers::invalid_unsigned_int,
1218  ExcIteratorPastEnd());
1219  Assert(this->col_block != numbers::invalid_unsigned_int,
1220  ExcIteratorPastEnd());
1221 
1222  return base_iterator->value();
1223  }
1224 
1225 
1226 
1227  template <class BlockMatrixType>
1228  inline void
1229  Accessor<BlockMatrixType, true>::advance()
1230  {
1231  Assert(this->row_block != numbers::invalid_unsigned_int,
1232  ExcIteratorPastEnd());
1233  Assert(this->col_block != numbers::invalid_unsigned_int,
1234  ExcIteratorPastEnd());
1235 
1236  // Remember current row inside block
1237  size_type local_row = base_iterator->row();
1238 
1239  // Advance one element inside the
1240  // current block
1241  ++base_iterator;
1242 
1243  // while we hit the end of the row of a
1244  // block (which may happen multiple
1245  // times if rows inside a block are
1246  // empty), we have to jump to the next
1247  // block and take the
1248  while (base_iterator ==
1249  matrix->block(this->row_block, this->col_block).end(local_row))
1250  {
1251  // jump to next block in this block
1252  // row, if possible, otherwise go
1253  // to next row
1254  if (this->col_block < matrix->n_block_cols() - 1)
1255  {
1256  ++this->col_block;
1257  base_iterator =
1258  matrix->block(this->row_block, this->col_block).begin(local_row);
1259  }
1260  else
1261  {
1262  // jump back to next row in
1263  // first block column
1264  this->col_block = 0;
1265  ++local_row;
1266 
1267  // see if this has brought us
1268  // past the number of rows in
1269  // this block. if so see
1270  // whether we've just fallen
1271  // off the end of the whole
1272  // matrix
1273  if (local_row ==
1274  matrix->block(this->row_block, this->col_block).m())
1275  {
1276  local_row = 0;
1277  ++this->row_block;
1278  if (this->row_block == matrix->n_block_rows())
1279  {
1280  this->row_block = numbers::invalid_unsigned_int;
1281  this->col_block = numbers::invalid_unsigned_int;
1282  return;
1283  }
1284  }
1285 
1286  base_iterator =
1287  matrix->block(this->row_block, this->col_block).begin(local_row);
1288  }
1289  }
1290  }
1291 
1292 
1293  template <class BlockMatrixType>
1294  inline bool
1295  Accessor<BlockMatrixType, true>::operator==(const Accessor &a) const
1296  {
1297  if (matrix != a.matrix)
1298  return false;
1299 
1300  if (this->row_block == a.row_block && this->col_block == a.col_block)
1301  // end iterators do not necessarily
1302  // have to have the same
1303  // base_iterator representation, but
1304  // valid iterators have to
1305  return (((this->row_block == numbers::invalid_unsigned_int) &&
1306  (this->col_block == numbers::invalid_unsigned_int)) ||
1307  (base_iterator == a.base_iterator));
1308 
1309  return false;
1310  }
1311 
1312  //----------------------------------------------------------------------//
1313 
1314 
1315  template <class BlockMatrixType>
1316  inline Accessor<BlockMatrixType, false>::Accessor(BlockMatrixType *matrix,
1317  const size_type row,
1318  const size_type col)
1319  : matrix(matrix)
1320  , base_iterator(matrix->block(0, 0).begin())
1321  {
1322  (void)col;
1323  Assert(col == 0, ExcNotImplemented());
1324  // check if this is a regular row or
1325  // the end of the matrix
1326  if (row < matrix->m())
1327  {
1328  const std::pair<unsigned int, size_type> indices =
1329  matrix->row_block_indices.global_to_local(row);
1330 
1331  // find the first block that does
1332  // have an entry in this row
1333  for (size_type bc = 0; bc < matrix->n_block_cols(); ++bc)
1334  {
1335  base_iterator =
1336  matrix->block(indices.first, bc).begin(indices.second);
1337  if (base_iterator !=
1338  matrix->block(indices.first, bc).end(indices.second))
1339  {
1340  this->row_block = indices.first;
1341  this->col_block = bc;
1342  return;
1343  }
1344  }
1345 
1346  // hm, there is no block that has
1347  // an entry in this column. we need
1348  // to take the next entry then,
1349  // which may be the first entry of
1350  // the next row, or recursively the
1351  // next row, or so on
1352  *this = Accessor(matrix, row + 1, 0);
1353  }
1354  else
1355  {
1356  // we were asked to create the end
1357  // iterator for this matrix
1358  this->row_block = numbers::invalid_size_type;
1359  this->col_block = numbers::invalid_size_type;
1360  }
1361  }
1362 
1363 
1364  template <class BlockMatrixType>
1365  inline typename Accessor<BlockMatrixType, false>::size_type
1366  Accessor<BlockMatrixType, false>::row() const
1367  {
1368  Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1369 
1370  return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1371  base_iterator->row());
1372  }
1373 
1374 
1375  template <class BlockMatrixType>
1376  inline typename Accessor<BlockMatrixType, false>::size_type
1377  Accessor<BlockMatrixType, false>::column() const
1378  {
1379  Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1380 
1381  return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
1382  base_iterator->column());
1383  }
1384 
1385 
1386  template <class BlockMatrixType>
1387  inline typename Accessor<BlockMatrixType, false>::value_type
1388  Accessor<BlockMatrixType, false>::value() const
1389  {
1390  Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1391  Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1392 
1393  return base_iterator->value();
1394  }
1395 
1396 
1397 
1398  template <class BlockMatrixType>
1399  inline void
1400  Accessor<BlockMatrixType, false>::set_value(
1401  typename Accessor<BlockMatrixType, false>::value_type newval) const
1402  {
1403  Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1404  Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1405 
1406  base_iterator->value() = newval;
1407  }
1408 
1409 
1410 
1411  template <class BlockMatrixType>
1412  inline void
1413  Accessor<BlockMatrixType, false>::advance()
1414  {
1415  Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1416  Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1417 
1418  // Remember current row inside block
1419  size_type local_row = base_iterator->row();
1420 
1421  // Advance one element inside the
1422  // current block
1423  ++base_iterator;
1424 
1425  // while we hit the end of the row of a
1426  // block (which may happen multiple
1427  // times if rows inside a block are
1428  // empty), we have to jump to the next
1429  // block and take the
1430  while (base_iterator ==
1431  matrix->block(this->row_block, this->col_block).end(local_row))
1432  {
1433  // jump to next block in this block
1434  // row, if possible, otherwise go
1435  // to next row
1436  if (this->col_block < matrix->n_block_cols() - 1)
1437  {
1438  ++this->col_block;
1439  base_iterator =
1440  matrix->block(this->row_block, this->col_block).begin(local_row);
1441  }
1442  else
1443  {
1444  // jump back to next row in
1445  // first block column
1446  this->col_block = 0;
1447  ++local_row;
1448 
1449  // see if this has brought us
1450  // past the number of rows in
1451  // this block. if so see
1452  // whether we've just fallen
1453  // off the end of the whole
1454  // matrix
1455  if (local_row ==
1456  matrix->block(this->row_block, this->col_block).m())
1457  {
1458  local_row = 0;
1459  ++this->row_block;
1460  if (this->row_block == matrix->n_block_rows())
1461  {
1462  this->row_block = numbers::invalid_size_type;
1463  this->col_block = numbers::invalid_size_type;
1464  return;
1465  }
1466  }
1467 
1468  base_iterator =
1469  matrix->block(this->row_block, this->col_block).begin(local_row);
1470  }
1471  }
1472  }
1473 
1474 
1475 
1476  template <class BlockMatrixType>
1477  inline bool
1478  Accessor<BlockMatrixType, false>::operator==(const Accessor &a) const
1479  {
1480  if (matrix != a.matrix)
1481  return false;
1482 
1483  if (this->row_block == a.row_block && this->col_block == a.col_block)
1484  // end iterators do not necessarily
1485  // have to have the same
1486  // base_iterator representation, but
1487  // valid iterators have to
1488  return (((this->row_block == numbers::invalid_size_type) &&
1489  (this->col_block == numbers::invalid_size_type)) ||
1490  (base_iterator == a.base_iterator));
1491 
1492  return false;
1493  }
1494 } // namespace BlockMatrixIterators
1495 
1496 
1497 //---------------------------------------------------------------------------
1498 
1499 template <typename MatrixType>
1501 {
1502  try
1503  {
1504  clear();
1505  }
1506  catch (...)
1507  {}
1508 }
1509 
1510 
1511 template <class MatrixType>
1512 template <class BlockMatrixType>
1514 BlockMatrixBase<MatrixType>::copy_from(const BlockMatrixType &source)
1515 {
1516  for (unsigned int r = 0; r < n_block_rows(); ++r)
1517  for (unsigned int c = 0; c < n_block_cols(); ++c)
1518  block(r, c).copy_from(source.block(r, c));
1519 
1520  return *this;
1521 }
1522 
1523 
1524 template <class MatrixType>
1525 std::size_t
1527 {
1528  std::size_t mem =
1530  MemoryConsumption::memory_consumption(column_block_indices) +
1535  sizeof(temporary_data.mutex);
1536 
1537  for (unsigned int r = 0; r < n_block_rows(); ++r)
1538  for (unsigned int c = 0; c < n_block_cols(); ++c)
1539  {
1540  MatrixType *p = this->sub_objects[r][c];
1542  }
1543 
1544  return mem;
1545 }
1546 
1547 
1548 
1549 template <class MatrixType>
1550 inline void
1552 {
1553  for (unsigned int r = 0; r < n_block_rows(); ++r)
1554  for (unsigned int c = 0; c < n_block_cols(); ++c)
1555  {
1556  MatrixType *p = this->sub_objects[r][c];
1557  this->sub_objects[r][c] = nullptr;
1558  delete p;
1559  }
1560  sub_objects.reinit(0, 0);
1561 
1562  // reset block indices to empty
1563  row_block_indices = column_block_indices = BlockIndices();
1564 }
1565 
1566 
1567 
1568 template <class MatrixType>
1570 BlockMatrixBase<MatrixType>::block(const unsigned int row,
1571  const unsigned int column)
1572 {
1573  Assert(row < n_block_rows(), ExcIndexRange(row, 0, n_block_rows()));
1574  Assert(column < n_block_cols(), ExcIndexRange(column, 0, n_block_cols()));
1575 
1576  return *sub_objects[row][column];
1577 }
1578 
1579 
1580 
1581 template <class MatrixType>
1582 inline const typename BlockMatrixBase<MatrixType>::BlockType &
1583 BlockMatrixBase<MatrixType>::block(const unsigned int row,
1584  const unsigned int column) const
1585 {
1586  Assert(row < n_block_rows(), ExcIndexRange(row, 0, n_block_rows()));
1587  Assert(column < n_block_cols(), ExcIndexRange(column, 0, n_block_cols()));
1588 
1589  return *sub_objects[row][column];
1590 }
1591 
1592 
1593 template <class MatrixType>
1594 inline typename BlockMatrixBase<MatrixType>::size_type
1596 {
1597  return row_block_indices.total_size();
1598 }
1599 
1600 
1601 
1602 template <class MatrixType>
1603 inline typename BlockMatrixBase<MatrixType>::size_type
1605 {
1606  return column_block_indices.total_size();
1607 }
1608 
1609 
1610 
1611 template <class MatrixType>
1612 inline unsigned int
1614 {
1615  return column_block_indices.size();
1616 }
1617 
1618 
1619 
1620 template <class MatrixType>
1621 inline unsigned int
1623 {
1624  return row_block_indices.size();
1625 }
1626 
1627 
1628 
1629 // Write the single set manually,
1630 // since the other function has a lot
1631 // of overhead in that case.
1632 template <class MatrixType>
1633 inline void
1634 BlockMatrixBase<MatrixType>::set(const size_type i,
1635  const size_type j,
1636  const value_type value)
1637 {
1639 
1640  AssertIsFinite(value);
1641 
1642  const std::pair<unsigned int, size_type>
1643  row_index = row_block_indices.global_to_local(i),
1644  col_index = column_block_indices.global_to_local(j);
1645  block(row_index.first, col_index.first)
1646  .set(row_index.second, col_index.second, value);
1647 }
1648 
1649 
1650 
1651 template <class MatrixType>
1652 template <typename number>
1653 inline void
1654 BlockMatrixBase<MatrixType>::set(const std::vector<size_type> &row_indices,
1655  const std::vector<size_type> &col_indices,
1656  const FullMatrix<number> & values,
1657  const bool elide_zero_values)
1658 {
1659  Assert(row_indices.size() == values.m(),
1660  ExcDimensionMismatch(row_indices.size(), values.m()));
1661  Assert(col_indices.size() == values.n(),
1662  ExcDimensionMismatch(col_indices.size(), values.n()));
1663 
1664  for (size_type i = 0; i < row_indices.size(); ++i)
1665  set(row_indices[i],
1666  col_indices.size(),
1667  col_indices.data(),
1668  &values(i, 0),
1669  elide_zero_values);
1670 }
1671 
1672 
1673 
1674 template <class MatrixType>
1675 template <typename number>
1676 inline void
1677 BlockMatrixBase<MatrixType>::set(const std::vector<size_type> &indices,
1678  const FullMatrix<number> & values,
1679  const bool elide_zero_values)
1680 {
1681  Assert(indices.size() == values.m(),
1682  ExcDimensionMismatch(indices.size(), values.m()));
1683  Assert(values.n() == values.m(), ExcNotQuadratic());
1684 
1685  for (size_type i = 0; i < indices.size(); ++i)
1686  set(indices[i],
1687  indices.size(),
1688  indices.data(),
1689  &values(i, 0),
1690  elide_zero_values);
1691 }
1692 
1693 
1694 
1695 template <class MatrixType>
1696 template <typename number>
1697 inline void
1698 BlockMatrixBase<MatrixType>::set(const size_type row,
1699  const std::vector<size_type> &col_indices,
1700  const std::vector<number> & values,
1701  const bool elide_zero_values)
1702 {
1703  Assert(col_indices.size() == values.size(),
1704  ExcDimensionMismatch(col_indices.size(), values.size()));
1705 
1706  set(row,
1707  col_indices.size(),
1708  col_indices.data(),
1709  values.data(),
1710  elide_zero_values);
1711 }
1712 
1713 
1714 
1715 // This is a very messy function, since
1716 // we need to calculate to each position
1717 // the location in the global array.
1718 template <class MatrixType>
1719 template <typename number>
1720 inline void
1721 BlockMatrixBase<MatrixType>::set(const size_type row,
1722  const size_type n_cols,
1723  const size_type *col_indices,
1724  const number * values,
1725  const bool elide_zero_values)
1726 {
1728 
1729  // lock access to the temporary data structure to
1730  // allow multiple threads to call this function concurrently
1731  std::lock_guard<std::mutex> lock(temporary_data.mutex);
1732 
1733  // Resize scratch arrays
1734  if (temporary_data.column_indices.size() < this->n_block_cols())
1735  {
1736  temporary_data.column_indices.resize(this->n_block_cols());
1737  temporary_data.column_values.resize(this->n_block_cols());
1739  }
1740 
1741  // Resize sub-arrays to n_cols. This
1742  // is a bit wasteful, but we resize
1743  // only a few times (then the maximum
1744  // row length won't increase that
1745  // much any more). At least we know
1746  // that all arrays are going to be of
1747  // the same size, so we can check
1748  // whether the size of one is large
1749  // enough before actually going
1750  // through all of them.
1751  if (temporary_data.column_indices[0].size() < n_cols)
1752  {
1753  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1754  {
1755  temporary_data.column_indices[i].resize(n_cols);
1756  temporary_data.column_values[i].resize(n_cols);
1757  }
1758  }
1759 
1760  // Reset the number of added elements
1761  // in each block to zero.
1762  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1764 
1765  // Go through the column indices to
1766  // find out which portions of the
1767  // values should be set in which
1768  // block of the matrix. We need to
1769  // touch all the data, since we can't
1770  // be sure that the data of one block
1771  // is stored contiguously (in fact,
1772  // indices will be intermixed when it
1773  // comes from an element matrix).
1774  for (size_type j = 0; j < n_cols; ++j)
1775  {
1776  number value = values[j];
1777 
1778  if (value == number() && elide_zero_values == true)
1779  continue;
1780 
1781  const std::pair<unsigned int, size_type> col_index =
1782  this->column_block_indices.global_to_local(col_indices[j]);
1783 
1784  const size_type local_index =
1785  temporary_data.counter_within_block[col_index.first]++;
1786 
1787  temporary_data.column_indices[col_index.first][local_index] =
1788  col_index.second;
1789  temporary_data.column_values[col_index.first][local_index] = value;
1790  }
1791 
1792 # ifdef DEBUG
1793  // If in debug mode, do a check whether
1794  // the right length has been obtained.
1795  size_type length = 0;
1796  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1797  length += temporary_data.counter_within_block[i];
1798  Assert(length <= n_cols, ExcInternalError());
1799 # endif
1800 
1801  // Now we found out about where the
1802  // individual columns should start and
1803  // where we should start reading out
1804  // data. Now let's write the data into
1805  // the individual blocks!
1806  const std::pair<unsigned int, size_type> row_index =
1808  for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
1809  {
1810  if (temporary_data.counter_within_block[block_col] == 0)
1811  continue;
1812 
1813  block(row_index.first, block_col)
1814  .set(row_index.second,
1816  temporary_data.column_indices[block_col].data(),
1817  temporary_data.column_values[block_col].data(),
1818  false);
1819  }
1820 }
1821 
1822 
1823 
1824 template <class MatrixType>
1825 inline void
1826 BlockMatrixBase<MatrixType>::add(const size_type i,
1827  const size_type j,
1828  const value_type value)
1829 {
1830  AssertIsFinite(value);
1831 
1833 
1834  // save some cycles for zero additions, but
1835  // only if it is safe for the matrix we are
1836  // working with
1837  using MatrixTraits = typename MatrixType::Traits;
1838  if ((MatrixTraits::zero_addition_can_be_elided == true) &&
1839  (value == value_type()))
1840  return;
1841 
1842  const std::pair<unsigned int, size_type>
1843  row_index = row_block_indices.global_to_local(i),
1844  col_index = column_block_indices.global_to_local(j);
1845  block(row_index.first, col_index.first)
1846  .add(row_index.second, col_index.second, value);
1847 }
1848 
1849 
1850 
1851 template <class MatrixType>
1852 template <typename number>
1853 inline void
1854 BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &row_indices,
1855  const std::vector<size_type> &col_indices,
1856  const FullMatrix<number> & values,
1857  const bool elide_zero_values)
1858 {
1859  Assert(row_indices.size() == values.m(),
1860  ExcDimensionMismatch(row_indices.size(), values.m()));
1861  Assert(col_indices.size() == values.n(),
1862  ExcDimensionMismatch(col_indices.size(), values.n()));
1863 
1864  for (size_type i = 0; i < row_indices.size(); ++i)
1865  add(row_indices[i],
1866  col_indices.size(),
1867  col_indices.data(),
1868  &values(i, 0),
1869  elide_zero_values);
1870 }
1871 
1872 
1873 
1874 template <class MatrixType>
1875 template <typename number>
1876 inline void
1877 BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &indices,
1878  const FullMatrix<number> & values,
1879  const bool elide_zero_values)
1880 {
1881  Assert(indices.size() == values.m(),
1882  ExcDimensionMismatch(indices.size(), values.m()));
1883  Assert(values.n() == values.m(), ExcNotQuadratic());
1884 
1885  for (size_type i = 0; i < indices.size(); ++i)
1886  add(indices[i],
1887  indices.size(),
1888  indices.data(),
1889  &values(i, 0),
1890  elide_zero_values);
1891 }
1892 
1893 
1894 
1895 template <class MatrixType>
1896 template <typename number>
1897 inline void
1898 BlockMatrixBase<MatrixType>::add(const size_type row,
1899  const std::vector<size_type> &col_indices,
1900  const std::vector<number> & values,
1901  const bool elide_zero_values)
1902 {
1903  Assert(col_indices.size() == values.size(),
1904  ExcDimensionMismatch(col_indices.size(), values.size()));
1905 
1906  add(row,
1907  col_indices.size(),
1908  col_indices.data(),
1909  values.data(),
1910  elide_zero_values);
1911 }
1912 
1913 
1914 
1915 // This is a very messy function, since
1916 // we need to calculate to each position
1917 // the location in the global array.
1918 template <class MatrixType>
1919 template <typename number>
1920 inline void
1921 BlockMatrixBase<MatrixType>::add(const size_type row,
1922  const size_type n_cols,
1923  const size_type *col_indices,
1924  const number * values,
1925  const bool elide_zero_values,
1926  const bool col_indices_are_sorted)
1927 {
1929 
1930  // TODO: Look over this to find out
1931  // whether we can do that more
1932  // efficiently.
1933  if (col_indices_are_sorted == true)
1934  {
1935 # ifdef DEBUG
1936  // check whether indices really are
1937  // sorted.
1938  size_type before = col_indices[0];
1939  for (size_type i = 1; i < n_cols; ++i)
1940  if (col_indices[i] <= before)
1941  Assert(false,
1942  ExcMessage("Flag col_indices_are_sorted is set, but "
1943  "indices appear to not be sorted.")) else before =
1944  col_indices[i];
1945 # endif
1946  const std::pair<unsigned int, size_type> row_index =
1948 
1949  if (this->n_block_cols() > 1)
1950  {
1951  const size_type *first_block =
1952  Utilities::lower_bound(col_indices,
1953  col_indices + n_cols,
1954  this->column_block_indices.block_start(1));
1955 
1956  const size_type n_zero_block_indices = first_block - col_indices;
1957  block(row_index.first, 0)
1958  .add(row_index.second,
1959  n_zero_block_indices,
1960  col_indices,
1961  values,
1962  elide_zero_values,
1963  col_indices_are_sorted);
1964 
1965  if (n_zero_block_indices < n_cols)
1966  this->add(row,
1967  n_cols - n_zero_block_indices,
1968  first_block,
1969  values + n_zero_block_indices,
1970  elide_zero_values,
1971  false);
1972  }
1973  else
1974  {
1975  block(row_index.first, 0)
1976  .add(row_index.second,
1977  n_cols,
1978  col_indices,
1979  values,
1980  elide_zero_values,
1981  col_indices_are_sorted);
1982  }
1983 
1984  return;
1985  }
1986 
1987  // Lock scratch arrays, then resize them
1988  std::lock_guard<std::mutex> lock(temporary_data.mutex);
1989 
1990  if (temporary_data.column_indices.size() < this->n_block_cols())
1991  {
1992  temporary_data.column_indices.resize(this->n_block_cols());
1993  temporary_data.column_values.resize(this->n_block_cols());
1995  }
1996 
1997  // Resize sub-arrays to n_cols. This
1998  // is a bit wasteful, but we resize
1999  // only a few times (then the maximum
2000  // row length won't increase that
2001  // much any more). At least we know
2002  // that all arrays are going to be of
2003  // the same size, so we can check
2004  // whether the size of one is large
2005  // enough before actually going
2006  // through all of them.
2007  if (temporary_data.column_indices[0].size() < n_cols)
2008  {
2009  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2010  {
2011  temporary_data.column_indices[i].resize(n_cols);
2012  temporary_data.column_values[i].resize(n_cols);
2013  }
2014  }
2015 
2016  // Reset the number of added elements
2017  // in each block to zero.
2018  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2020 
2021  // Go through the column indices to
2022  // find out which portions of the
2023  // values should be written into
2024  // which block of the matrix. We need
2025  // to touch all the data, since we
2026  // can't be sure that the data of one
2027  // block is stored contiguously (in
2028  // fact, data will be intermixed when
2029  // it comes from an element matrix).
2030  for (size_type j = 0; j < n_cols; ++j)
2031  {
2032  number value = values[j];
2033 
2034  if (value == number() && elide_zero_values == true)
2035  continue;
2036 
2037  const std::pair<unsigned int, size_type> col_index =
2038  this->column_block_indices.global_to_local(col_indices[j]);
2039 
2040  const size_type local_index =
2041  temporary_data.counter_within_block[col_index.first]++;
2042 
2043  temporary_data.column_indices[col_index.first][local_index] =
2044  col_index.second;
2045  temporary_data.column_values[col_index.first][local_index] = value;
2046  }
2047 
2048 # ifdef DEBUG
2049  // If in debug mode, do a check whether
2050  // the right length has been obtained.
2051  size_type length = 0;
2052  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2053  length += temporary_data.counter_within_block[i];
2054  Assert(length <= n_cols, ExcInternalError());
2055 # endif
2056 
2057  // Now we found out about where the
2058  // individual columns should start and
2059  // where we should start reading out
2060  // data. Now let's write the data into
2061  // the individual blocks!
2062  const std::pair<unsigned int, size_type> row_index =
2064  for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
2065  {
2066  if (temporary_data.counter_within_block[block_col] == 0)
2067  continue;
2068 
2069  block(row_index.first, block_col)
2070  .add(row_index.second,
2072  temporary_data.column_indices[block_col].data(),
2073  temporary_data.column_values[block_col].data(),
2074  false,
2075  col_indices_are_sorted);
2076  }
2077 }
2078 
2079 
2080 
2081 template <class MatrixType>
2082 inline void
2084  const BlockMatrixBase<MatrixType> &matrix)
2085 {
2086  AssertIsFinite(factor);
2087 
2089 
2090  // save some cycles for zero additions, but
2091  // only if it is safe for the matrix we are
2092  // working with
2093  using MatrixTraits = typename MatrixType::Traits;
2094  if ((MatrixTraits::zero_addition_can_be_elided == true) && (factor == 0))
2095  return;
2096 
2097  for (unsigned int row = 0; row < n_block_rows(); ++row)
2098  for (unsigned int col = 0; col < n_block_cols(); ++col)
2099  // This function should throw if the sparsity
2100  // patterns of the two blocks differ
2101  block(row, col).add(factor, matrix.block(row, col));
2102 }
2103 
2104 
2105 
2106 template <class MatrixType>
2108 BlockMatrixBase<MatrixType>::operator()(const size_type i,
2109  const size_type j) const
2110 {
2111  const std::pair<unsigned int, size_type>
2112  row_index = row_block_indices.global_to_local(i),
2113  col_index = column_block_indices.global_to_local(j);
2114  return block(row_index.first, col_index.first)(row_index.second,
2115  col_index.second);
2116 }
2117 
2118 
2119 
2120 template <class MatrixType>
2122 BlockMatrixBase<MatrixType>::el(const size_type i, const size_type j) const
2123 {
2124  const std::pair<unsigned int, size_type>
2125  row_index = row_block_indices.global_to_local(i),
2126  col_index = column_block_indices.global_to_local(j);
2127  return block(row_index.first, col_index.first)
2128  .el(row_index.second, col_index.second);
2129 }
2130 
2131 
2132 
2133 template <class MatrixType>
2135 BlockMatrixBase<MatrixType>::diag_element(const size_type i) const
2136 {
2138 
2139  const std::pair<unsigned int, size_type> index =
2141  return block(index.first, index.first).diag_element(index.second);
2142 }
2143 
2144 
2145 
2146 template <class MatrixType>
2147 inline void
2149  ::VectorOperation::values operation)
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 <class MatrixType>
2161 {
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 <class MatrixType>
2177 {
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 <class MatrixType>
2194 const BlockIndices &
2196 {
2197  return this->row_block_indices;
2198 }
2199 
2200 
2201 
2202 template <class MatrixType>
2203 const BlockIndices &
2205 {
2206  return this->column_block_indices;
2207 }
2208 
2209 
2210 
2211 template <class MatrixType>
2212 template <class 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 <class MatrixType>
2233 template <class BlockVectorType, class VectorType>
2234 void
2236  VectorType & dst,
2237  const BlockVectorType &src) const
2238 {
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 <class MatrixType>
2251 template <class BlockVectorType, class 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()));
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 <class MatrixType>
2267 template <class VectorType>
2268 void
2270  VectorType & dst,
2271  const VectorType &src) const
2272 {
2275 
2276  block(0, 0).vmult(dst, src);
2277 }
2278 
2279 
2280 
2281 template <class MatrixType>
2282 template <class 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 <class MatrixType>
2300 template <class 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 <class MatrixType>
2323 template <class BlockVectorType, class 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()));
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 <class MatrixType>
2341 template <class BlockVectorType, class VectorType>
2342 void
2344  VectorType & dst,
2345  const BlockVectorType &src) const
2346 {
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 <class MatrixType>
2360 template <class VectorType>
2361 void
2363  VectorType & dst,
2364  const VectorType &src) const
2365 {
2368 
2369  block(0, 0).Tvmult(dst, src);
2370 }
2371 
2372 
2373 
2374 template <class MatrixType>
2375 template <class 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 <class MatrixType>
2393 template <class BlockVectorType>
2395 BlockMatrixBase<MatrixType>::matrix_norm_square(const BlockVectorType &v) const
2396 {
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 <class MatrixType>
2415 typename BlockMatrixBase<MatrixType>::real_type
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 <class MatrixType>
2437 template <class 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 <class MatrixType>
2459 template <class 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 <class 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 <class MatrixType>
2527 {
2528  return const_iterator(this, 0);
2529 }
2530 
2531 
2532 
2533 template <class MatrixType>
2536 {
2537  return const_iterator(this, m());
2538 }
2539 
2540 
2541 
2542 template <class MatrixType>
2544 BlockMatrixBase<MatrixType>::begin(const size_type r) const
2545 {
2546  Assert(r < m(), ExcIndexRange(r, 0, m()));
2547  return const_iterator(this, r);
2548 }
2549 
2550 
2551 
2552 template <class MatrixType>
2554 BlockMatrixBase<MatrixType>::end(const size_type r) const
2555 {
2556  Assert(r < m(), ExcIndexRange(r, 0, m()));
2557  return const_iterator(this, r + 1);
2558 }
2559 
2560 
2561 
2562 template <class MatrixType>
2565 {
2566  return iterator(this, 0);
2567 }
2568 
2569 
2570 
2571 template <class MatrixType>
2574 {
2575  return iterator(this, m());
2576 }
2577 
2578 
2579 
2580 template <class MatrixType>
2582 BlockMatrixBase<MatrixType>::begin(const size_type r)
2583 {
2584  Assert(r < m(), ExcIndexRange(r, 0, m()));
2585  return iterator(this, r);
2586 }
2587 
2588 
2589 
2590 template <class MatrixType>
2592 BlockMatrixBase<MatrixType>::end(const size_type r)
2593 {
2594  Assert(r < m(), ExcIndexRange(r, 0, m()));
2595  return iterator(this, r + 1);
2596 }
2597 
2598 
2599 
2600 template <class 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 <class 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 <class 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 
2662 DEAL_II_NAMESPACE_CLOSE
2663 
2664 #endif // dealii_block_matrix_base_h
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1063
size_type m() const
const types::global_dof_index invalid_size_type
Definition: types.h:196
static const unsigned int invalid_unsigned_int
Definition: types.h:187
void collect_sizes()
void vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
BlockMatrixBase & copy_from(const BlockMatrixType &source)
void Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const
TemporaryData & operator=(const TemporaryData &)
std::vector< size_type > counter_within_block
std::vector< std::vector< size_type > > column_indices
void set(const size_type i, const size_type j, const value_type value)
value_type operator()(const size_type i, const size_type j) const
number value_type
real_type frobenius_norm() const
void vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
value_type diag_element(const size_type i) const
static ::ExceptionBase & ExcNotInitialized()
void Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
void add(const size_type i, const size_type j, const value_type value)
value_type matrix_norm_square(const BlockVectorType &v) const
void Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
TemporaryData temporary_data
static ::ExceptionBase & ExcDivideByZero()
std::size_t memory_consumption() const
void vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
size_type total_size() const
BlockMatrixBase & operator*=(const value_type factor)
typename BlockType::value_type value_type
void print(std::ostream &out, const bool alternative_output=false) const
void vmult_add(BlockVectorType &dst, const BlockVectorType &src) const
void Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
const BlockIndices & get_row_indices() const
std::vector< std::vector< value_type > > column_values
BlockMatrixBase & operator/=(const value_type factor)
Table< 2, SmartPointer< BlockType, BlockMatrixBase< MatrixType > > > sub_objects
unsigned int n_block_cols() const
static ::ExceptionBase & ExcMessage(std::string arg1)
size_type n() const
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
~BlockMatrixBase() override
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void reinit(const unsigned int n_blocks, const size_type n_elements_per_block)
void prepare_add_operation()
value_type el(const size_type i, const size_type j) const
typename BlockMatrixType::value_type value_type
const BlockIndices & get_column_indices() const
iterator end()
BlockMatrixType::BlockType::const_iterator base_iterator
value_type residual(BlockVectorType &dst, const BlockVectorType &x, const BlockVectorType &b) const
unsigned int block_column() const
static ::ExceptionBase & ExcIteratorPastEnd()
static ::ExceptionBase & ExcNotQuadratic()
unsigned int global_dof_index
Definition: types.h:89
unsigned int block_row() const
types::global_dof_index size_type
void vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
BlockType & block(const unsigned int row, const unsigned int column)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:587
size_type block_start(const unsigned int i) const
static ::ExceptionBase & ExcNotImplemented()
BlockIndices row_block_indices
Definition: table.h:699
unsigned int n_block_rows() const
unsigned int size() const
value_type matrix_scalar_product(const BlockVectorType &u, const BlockVectorType &v) const
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
#define AssertIsFinite(number)
Definition: exceptions.h:1669
size_type m() const
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
void compress(::VectorOperation::values operation)
static ::ExceptionBase & ExcIncompatibleRowNumbers(int arg1, int arg2, int arg3, int arg4)
size_type n() const
void prepare_set_operation()
static ::ExceptionBase & ExcInternalError()
iterator begin()