Reference documentation for deal.II version 9.1.1
\(\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 template <typename>
42 
43 
44 
55 {
60  template <class BlockMatrixType>
62  {
63  public:
68 
72  using value_type = typename BlockMatrixType::value_type;
73 
77  AccessorBase();
78 
82  unsigned int
83  block_row() const;
84 
88  unsigned int
89  block_column() const;
90 
91  protected:
95  unsigned int row_block;
96 
100  unsigned int col_block;
101 
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 
286  template <typename>
287  friend class ::MatrixIterator;
288  };
289 } // namespace BlockMatrixIterators
290 
291 
292 
353 template <typename MatrixType>
354 class BlockMatrixBase : public Subscriptor
355 {
356 public:
360  using BlockType = MatrixType;
361 
367  using pointer = value_type *;
368  using const_pointer = const value_type *;
369  using reference = value_type &;
370  using const_reference = const value_type &;
371  using size_type = types::global_dof_index;
372 
373  using iterator =
375 
376  using const_iterator =
378 
379 
383  BlockMatrixBase() = default;
384 
388  ~BlockMatrixBase() override;
389 
406  template <class BlockMatrixType>
408  copy_from(const BlockMatrixType &source);
409 
413  BlockType &
414  block(const unsigned int row, const unsigned int column);
415 
416 
421  const BlockType &
422  block(const unsigned int row, const unsigned int column) const;
423 
428  size_type
429  m() const;
430 
435  size_type
436  n() const;
437 
438 
443  unsigned int
444  n_block_rows() const;
445 
450  unsigned int
451  n_block_cols() const;
452 
458  void
459  set(const size_type i, const size_type j, const value_type value);
460 
476  template <typename number>
477  void
478  set(const std::vector<size_type> &indices,
479  const FullMatrix<number> & full_matrix,
480  const bool elide_zero_values = false);
481 
487  template <typename number>
488  void
489  set(const std::vector<size_type> &row_indices,
490  const std::vector<size_type> &col_indices,
491  const FullMatrix<number> & full_matrix,
492  const bool elide_zero_values = false);
493 
504  template <typename number>
505  void
506  set(const size_type row,
507  const std::vector<size_type> &col_indices,
508  const std::vector<number> & values,
509  const bool elide_zero_values = false);
510 
520  template <typename number>
521  void
522  set(const size_type row,
523  const size_type n_cols,
524  const size_type *col_indices,
525  const number * values,
526  const bool elide_zero_values = false);
527 
533  void
534  add(const size_type i, const size_type j, const value_type value);
535 
550  template <typename number>
551  void
552  add(const std::vector<size_type> &indices,
553  const FullMatrix<number> & full_matrix,
554  const bool elide_zero_values = true);
555 
561  template <typename number>
562  void
563  add(const std::vector<size_type> &row_indices,
564  const std::vector<size_type> &col_indices,
565  const FullMatrix<number> & full_matrix,
566  const bool elide_zero_values = true);
567 
577  template <typename number>
578  void
579  add(const size_type row,
580  const std::vector<size_type> &col_indices,
581  const std::vector<number> & values,
582  const bool elide_zero_values = true);
583 
593  template <typename number>
594  void
595  add(const size_type row,
596  const size_type n_cols,
597  const size_type *col_indices,
598  const number * values,
599  const bool elide_zero_values = true,
600  const bool col_indices_are_sorted = false);
601 
613  void
614  add(const value_type factor, const BlockMatrixBase<MatrixType> &matrix);
615 
622  value_type
623  operator()(const size_type i, const size_type j) const;
624 
633  value_type
634  el(const size_type i, const size_type j) const;
635 
646  value_type
647  diag_element(const size_type i) const;
648 
657  void
658  compress(::VectorOperation::values operation);
659 
664  operator*=(const value_type factor);
665 
670  operator/=(const value_type factor);
671 
676  template <class BlockVectorType>
677  void
678  vmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
679 
685  template <class BlockVectorType>
686  void
687  Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
688 
701  template <class BlockVectorType>
702  value_type
703  matrix_norm_square(const BlockVectorType &v) const;
704 
708  template <class BlockVectorType>
709  value_type
710  matrix_scalar_product(const BlockVectorType &u,
711  const BlockVectorType &v) const;
712 
716  template <class BlockVectorType>
717  value_type
718  residual(BlockVectorType & dst,
719  const BlockVectorType &x,
720  const BlockVectorType &b) const;
721 
728  void
729  print(std::ostream &out, const bool alternative_output = false) const;
730 
734  iterator
735  begin();
736 
740  iterator
741  end();
742 
746  iterator
747  begin(const size_type r);
748 
752  iterator
753  end(const size_type r);
758  begin() const;
759 
764  end() const;
765 
770  begin(const size_type r) const;
771 
776  end(const size_type r) const;
777 
781  const BlockIndices &
782  get_row_indices() const;
783 
787  const BlockIndices &
788  get_column_indices() const;
789 
795  std::size_t
796  memory_consumption() const;
797 
806  DeclException4(ExcIncompatibleRowNumbers,
807  int,
808  int,
809  int,
810  int,
811  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
812  << ',' << arg4 << "] have differing row numbers.");
816  DeclException4(ExcIncompatibleColNumbers,
817  int,
818  int,
819  int,
820  int,
821  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
822  << ',' << arg4 << "] have differing column numbers.");
824 protected:
837  void
838  clear();
839 
844  BlockIndices column_block_indices;
845 
850 
869  void
870  collect_sizes();
871 
882  template <class BlockVectorType>
883  void
884  vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const;
885 
896  template <class BlockVectorType, class VectorType>
897  void
898  vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const;
899 
910  template <class BlockVectorType, class VectorType>
911  void
912  vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const;
913 
924  template <class VectorType>
925  void
926  vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const;
927 
939  template <class BlockVectorType>
940  void
941  Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const;
942 
953  template <class BlockVectorType, class VectorType>
954  void
955  Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const;
956 
967  template <class BlockVectorType, class VectorType>
968  void
969  Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const;
970 
981  template <class VectorType>
982  void
983  Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const;
984 
985 
986 protected:
993  void
994  prepare_add_operation();
995 
1000  void
1001  prepare_set_operation();
1002 
1003 
1004 private:
1014  {
1019  std::vector<size_type> counter_within_block;
1020 
1025  std::vector<std::vector<size_type>> column_indices;
1026 
1031  std::vector<std::vector<value_type>> column_values;
1032 
1038 
1048  TemporaryData &
1050  {
1051  return *this;
1052  }
1053  };
1054 
1061  TemporaryData temporary_data;
1062 
1067  template <typename, bool>
1069 
1070  template <typename>
1071  friend class MatrixIterator;
1072 };
1073 
1074 
1077 #ifndef DOXYGEN
1078 /* ------------------------- Template functions ---------------------- */
1079 
1080 
1081 namespace BlockMatrixIterators
1082 {
1083  template <class BlockMatrixType>
1085  : row_block(0)
1086  , col_block(0)
1087  {}
1088 
1089 
1090  template <class BlockMatrixType>
1091  inline unsigned int
1092  AccessorBase<BlockMatrixType>::block_row() const
1093  {
1095 
1096  return row_block;
1097  }
1098 
1099 
1100  template <class BlockMatrixType>
1101  inline unsigned int
1102  AccessorBase<BlockMatrixType>::block_column() const
1103  {
1105 
1106  return col_block;
1107  }
1108 
1109 
1110  template <class BlockMatrixType>
1111  inline Accessor<BlockMatrixType, true>::Accessor(
1112  const BlockMatrixType *matrix,
1113  const size_type row,
1114  const size_type col)
1115  : matrix(matrix)
1116  , base_iterator(matrix->block(0, 0).begin())
1117  {
1118  (void)col;
1119  Assert(col == 0, ExcNotImplemented());
1120 
1121  // check if this is a regular row or
1122  // the end of the matrix
1123  if (row < matrix->m())
1124  {
1125  const std::pair<unsigned int, size_type> indices =
1126  matrix->row_block_indices.global_to_local(row);
1127 
1128  // find the first block that does
1129  // have an entry in this row
1130  for (unsigned int bc = 0; bc < matrix->n_block_cols(); ++bc)
1131  {
1132  base_iterator =
1133  matrix->block(indices.first, bc).begin(indices.second);
1134  if (base_iterator !=
1135  matrix->block(indices.first, bc).end(indices.second))
1136  {
1137  this->row_block = indices.first;
1138  this->col_block = bc;
1139  return;
1140  }
1141  }
1142 
1143  // hm, there is no block that has
1144  // an entry in this column. we need
1145  // to take the next entry then,
1146  // which may be the first entry of
1147  // the next row, or recursively the
1148  // next row, or so on
1149  *this = Accessor(matrix, row + 1, 0);
1150  }
1151  else
1152  {
1153  // we were asked to create the end
1154  // iterator for this matrix
1155  this->row_block = numbers::invalid_unsigned_int;
1156  this->col_block = numbers::invalid_unsigned_int;
1157  }
1158  }
1159 
1160 
1161  // template <class BlockMatrixType>
1162  // inline
1163  // Accessor<BlockMatrixType, true>::Accessor (const
1164  // Accessor<BlockMatrixType, true>& other)
1165  // :
1166  // matrix(other.matrix),
1167  // base_iterator(other.base_iterator)
1168  // {
1169  // this->row_block = other.row_block;
1170  // this->col_block = other.col_block;
1171  // }
1172 
1173 
1174  template <class BlockMatrixType>
1175  inline Accessor<BlockMatrixType, true>::Accessor(
1176  const Accessor<BlockMatrixType, false> &other)
1177  : matrix(other.matrix)
1178  , base_iterator(other.base_iterator)
1179  {
1180  this->row_block = other.row_block;
1181  this->col_block = other.col_block;
1182  }
1183 
1184 
1185  template <class BlockMatrixType>
1186  inline typename Accessor<BlockMatrixType, true>::size_type
1187  Accessor<BlockMatrixType, true>::row() const
1188  {
1189  Assert(this->row_block != numbers::invalid_unsigned_int,
1190  ExcIteratorPastEnd());
1191 
1192  return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1193  base_iterator->row());
1194  }
1195 
1196 
1197  template <class BlockMatrixType>
1198  inline typename Accessor<BlockMatrixType, true>::size_type
1199  Accessor<BlockMatrixType, true>::column() const
1200  {
1201  Assert(this->col_block != numbers::invalid_unsigned_int,
1202  ExcIteratorPastEnd());
1203 
1204  return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
1205  base_iterator->column());
1206  }
1207 
1208 
1209  template <class BlockMatrixType>
1210  inline typename Accessor<BlockMatrixType, true>::value_type
1211  Accessor<BlockMatrixType, true>::value() const
1212  {
1213  Assert(this->row_block != numbers::invalid_unsigned_int,
1214  ExcIteratorPastEnd());
1215  Assert(this->col_block != numbers::invalid_unsigned_int,
1216  ExcIteratorPastEnd());
1217 
1218  return base_iterator->value();
1219  }
1220 
1221 
1222 
1223  template <class BlockMatrixType>
1224  inline void
1225  Accessor<BlockMatrixType, true>::advance()
1226  {
1227  Assert(this->row_block != numbers::invalid_unsigned_int,
1228  ExcIteratorPastEnd());
1229  Assert(this->col_block != numbers::invalid_unsigned_int,
1230  ExcIteratorPastEnd());
1231 
1232  // Remember current row inside block
1233  size_type local_row = base_iterator->row();
1234 
1235  // Advance one element inside the
1236  // current block
1237  ++base_iterator;
1238 
1239  // while we hit the end of the row of a
1240  // block (which may happen multiple
1241  // times if rows inside a block are
1242  // empty), we have to jump to the next
1243  // block and take the
1244  while (base_iterator ==
1245  matrix->block(this->row_block, this->col_block).end(local_row))
1246  {
1247  // jump to next block in this block
1248  // row, if possible, otherwise go
1249  // to next row
1250  if (this->col_block < matrix->n_block_cols() - 1)
1251  {
1252  ++this->col_block;
1253  base_iterator =
1254  matrix->block(this->row_block, this->col_block).begin(local_row);
1255  }
1256  else
1257  {
1258  // jump back to next row in
1259  // first block column
1260  this->col_block = 0;
1261  ++local_row;
1262 
1263  // see if this has brought us
1264  // past the number of rows in
1265  // this block. if so see
1266  // whether we've just fallen
1267  // off the end of the whole
1268  // matrix
1269  if (local_row ==
1270  matrix->block(this->row_block, this->col_block).m())
1271  {
1272  local_row = 0;
1273  ++this->row_block;
1274  if (this->row_block == matrix->n_block_rows())
1275  {
1276  this->row_block = numbers::invalid_unsigned_int;
1277  this->col_block = numbers::invalid_unsigned_int;
1278  return;
1279  }
1280  }
1281 
1282  base_iterator =
1283  matrix->block(this->row_block, this->col_block).begin(local_row);
1284  }
1285  }
1286  }
1287 
1288 
1289  template <class BlockMatrixType>
1290  inline bool
1291  Accessor<BlockMatrixType, true>::operator==(const Accessor &a) const
1292  {
1293  if (matrix != a.matrix)
1294  return false;
1295 
1296  if (this->row_block == a.row_block && this->col_block == a.col_block)
1297  // end iterators do not necessarily
1298  // have to have the same
1299  // base_iterator representation, but
1300  // valid iterators have to
1301  return (((this->row_block == numbers::invalid_unsigned_int) &&
1302  (this->col_block == numbers::invalid_unsigned_int)) ||
1303  (base_iterator == a.base_iterator));
1304 
1305  return false;
1306  }
1307 
1308  //----------------------------------------------------------------------//
1309 
1310 
1311  template <class BlockMatrixType>
1312  inline Accessor<BlockMatrixType, false>::Accessor(BlockMatrixType *matrix,
1313  const size_type row,
1314  const size_type col)
1315  : matrix(matrix)
1316  , base_iterator(matrix->block(0, 0).begin())
1317  {
1318  (void)col;
1319  Assert(col == 0, ExcNotImplemented());
1320  // check if this is a regular row or
1321  // the end of the matrix
1322  if (row < matrix->m())
1323  {
1324  const std::pair<unsigned int, size_type> indices =
1325  matrix->row_block_indices.global_to_local(row);
1326 
1327  // find the first block that does
1328  // have an entry in this row
1329  for (size_type bc = 0; bc < matrix->n_block_cols(); ++bc)
1330  {
1331  base_iterator =
1332  matrix->block(indices.first, bc).begin(indices.second);
1333  if (base_iterator !=
1334  matrix->block(indices.first, bc).end(indices.second))
1335  {
1336  this->row_block = indices.first;
1337  this->col_block = bc;
1338  return;
1339  }
1340  }
1341 
1342  // hm, there is no block that has
1343  // an entry in this column. we need
1344  // to take the next entry then,
1345  // which may be the first entry of
1346  // the next row, or recursively the
1347  // next row, or so on
1348  *this = Accessor(matrix, row + 1, 0);
1349  }
1350  else
1351  {
1352  // we were asked to create the end
1353  // iterator for this matrix
1354  this->row_block = numbers::invalid_size_type;
1355  this->col_block = numbers::invalid_size_type;
1356  }
1357  }
1358 
1359 
1360  template <class BlockMatrixType>
1361  inline typename Accessor<BlockMatrixType, false>::size_type
1362  Accessor<BlockMatrixType, false>::row() const
1363  {
1364  Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1365 
1366  return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1367  base_iterator->row());
1368  }
1369 
1370 
1371  template <class BlockMatrixType>
1372  inline typename Accessor<BlockMatrixType, false>::size_type
1373  Accessor<BlockMatrixType, false>::column() const
1374  {
1375  Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1376 
1377  return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
1378  base_iterator->column());
1379  }
1380 
1381 
1382  template <class BlockMatrixType>
1383  inline typename Accessor<BlockMatrixType, false>::value_type
1384  Accessor<BlockMatrixType, false>::value() const
1385  {
1386  Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1387  Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1388 
1389  return base_iterator->value();
1390  }
1391 
1392 
1393 
1394  template <class BlockMatrixType>
1395  inline void
1396  Accessor<BlockMatrixType, false>::set_value(
1397  typename Accessor<BlockMatrixType, false>::value_type newval) const
1398  {
1399  Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1400  Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1401 
1402  base_iterator->value() = newval;
1403  }
1404 
1405 
1406 
1407  template <class BlockMatrixType>
1408  inline void
1409  Accessor<BlockMatrixType, false>::advance()
1410  {
1411  Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1412  Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1413 
1414  // Remember current row inside block
1415  size_type local_row = base_iterator->row();
1416 
1417  // Advance one element inside the
1418  // current block
1419  ++base_iterator;
1420 
1421  // while we hit the end of the row of a
1422  // block (which may happen multiple
1423  // times if rows inside a block are
1424  // empty), we have to jump to the next
1425  // block and take the
1426  while (base_iterator ==
1427  matrix->block(this->row_block, this->col_block).end(local_row))
1428  {
1429  // jump to next block in this block
1430  // row, if possible, otherwise go
1431  // to next row
1432  if (this->col_block < matrix->n_block_cols() - 1)
1433  {
1434  ++this->col_block;
1435  base_iterator =
1436  matrix->block(this->row_block, this->col_block).begin(local_row);
1437  }
1438  else
1439  {
1440  // jump back to next row in
1441  // first block column
1442  this->col_block = 0;
1443  ++local_row;
1444 
1445  // see if this has brought us
1446  // past the number of rows in
1447  // this block. if so see
1448  // whether we've just fallen
1449  // off the end of the whole
1450  // matrix
1451  if (local_row ==
1452  matrix->block(this->row_block, this->col_block).m())
1453  {
1454  local_row = 0;
1455  ++this->row_block;
1456  if (this->row_block == matrix->n_block_rows())
1457  {
1458  this->row_block = numbers::invalid_size_type;
1459  this->col_block = numbers::invalid_size_type;
1460  return;
1461  }
1462  }
1463 
1464  base_iterator =
1465  matrix->block(this->row_block, this->col_block).begin(local_row);
1466  }
1467  }
1468  }
1469 
1470 
1471 
1472  template <class BlockMatrixType>
1473  inline bool
1474  Accessor<BlockMatrixType, false>::operator==(const Accessor &a) const
1475  {
1476  if (matrix != a.matrix)
1477  return false;
1478 
1479  if (this->row_block == a.row_block && this->col_block == a.col_block)
1480  // end iterators do not necessarily
1481  // have to have the same
1482  // base_iterator representation, but
1483  // valid iterators have to
1484  return (((this->row_block == numbers::invalid_size_type) &&
1485  (this->col_block == numbers::invalid_size_type)) ||
1486  (base_iterator == a.base_iterator));
1487 
1488  return false;
1489  }
1490 } // namespace BlockMatrixIterators
1491 
1492 
1493 //---------------------------------------------------------------------------
1494 
1495 template <typename MatrixType>
1497 {
1498  try
1499  {
1500  clear();
1501  }
1502  catch (...)
1503  {}
1504 }
1505 
1506 
1507 template <class MatrixType>
1508 template <class BlockMatrixType>
1510 BlockMatrixBase<MatrixType>::copy_from(const BlockMatrixType &source)
1511 {
1512  for (unsigned int r = 0; r < n_block_rows(); ++r)
1513  for (unsigned int c = 0; c < n_block_cols(); ++c)
1514  block(r, c).copy_from(source.block(r, c));
1515 
1516  return *this;
1517 }
1518 
1519 
1520 template <class MatrixType>
1521 std::size_t
1523 {
1524  std::size_t mem =
1526  MemoryConsumption::memory_consumption(column_block_indices) +
1531  sizeof(temporary_data.mutex);
1532 
1533  for (unsigned int r = 0; r < n_block_rows(); ++r)
1534  for (unsigned int c = 0; c < n_block_cols(); ++c)
1535  {
1536  MatrixType *p = this->sub_objects[r][c];
1538  }
1539 
1540  return mem;
1541 }
1542 
1543 
1544 
1545 template <class MatrixType>
1546 inline void
1548 {
1549  for (unsigned int r = 0; r < n_block_rows(); ++r)
1550  for (unsigned int c = 0; c < n_block_cols(); ++c)
1551  {
1552  MatrixType *p = this->sub_objects[r][c];
1553  this->sub_objects[r][c] = nullptr;
1554  delete p;
1555  }
1556  sub_objects.reinit(0, 0);
1557 
1558  // reset block indices to empty
1559  row_block_indices = column_block_indices = BlockIndices();
1560 }
1561 
1562 
1563 
1564 template <class MatrixType>
1566 BlockMatrixBase<MatrixType>::block(const unsigned int row,
1567  const unsigned int column)
1568 {
1569  Assert(row < n_block_rows(), ExcIndexRange(row, 0, n_block_rows()));
1570  Assert(column < n_block_cols(), ExcIndexRange(column, 0, n_block_cols()));
1571 
1572  return *sub_objects[row][column];
1573 }
1574 
1575 
1576 
1577 template <class MatrixType>
1578 inline const typename BlockMatrixBase<MatrixType>::BlockType &
1579 BlockMatrixBase<MatrixType>::block(const unsigned int row,
1580  const unsigned int column) const
1581 {
1582  Assert(row < n_block_rows(), ExcIndexRange(row, 0, n_block_rows()));
1583  Assert(column < n_block_cols(), ExcIndexRange(column, 0, n_block_cols()));
1584 
1585  return *sub_objects[row][column];
1586 }
1587 
1588 
1589 template <class MatrixType>
1590 inline typename BlockMatrixBase<MatrixType>::size_type
1592 {
1593  return row_block_indices.total_size();
1594 }
1595 
1596 
1597 
1598 template <class MatrixType>
1599 inline typename BlockMatrixBase<MatrixType>::size_type
1601 {
1602  return column_block_indices.total_size();
1603 }
1604 
1605 
1606 
1607 template <class MatrixType>
1608 inline unsigned int
1610 {
1611  return column_block_indices.size();
1612 }
1613 
1614 
1615 
1616 template <class MatrixType>
1617 inline unsigned int
1619 {
1620  return row_block_indices.size();
1621 }
1622 
1623 
1624 
1625 // Write the single set manually,
1626 // since the other function has a lot
1627 // of overhead in that case.
1628 template <class MatrixType>
1629 inline void
1630 BlockMatrixBase<MatrixType>::set(const size_type i,
1631  const size_type j,
1632  const value_type value)
1633 {
1635 
1636  AssertIsFinite(value);
1637 
1638  const std::pair<unsigned int, size_type>
1639  row_index = row_block_indices.global_to_local(i),
1640  col_index = column_block_indices.global_to_local(j);
1641  block(row_index.first, col_index.first)
1642  .set(row_index.second, col_index.second, value);
1643 }
1644 
1645 
1646 
1647 template <class MatrixType>
1648 template <typename number>
1649 inline void
1650 BlockMatrixBase<MatrixType>::set(const std::vector<size_type> &row_indices,
1651  const std::vector<size_type> &col_indices,
1652  const FullMatrix<number> & values,
1653  const bool elide_zero_values)
1654 {
1655  Assert(row_indices.size() == values.m(),
1656  ExcDimensionMismatch(row_indices.size(), values.m()));
1657  Assert(col_indices.size() == values.n(),
1658  ExcDimensionMismatch(col_indices.size(), values.n()));
1659 
1660  for (size_type i = 0; i < row_indices.size(); ++i)
1661  set(row_indices[i],
1662  col_indices.size(),
1663  col_indices.data(),
1664  &values(i, 0),
1665  elide_zero_values);
1666 }
1667 
1668 
1669 
1670 template <class MatrixType>
1671 template <typename number>
1672 inline void
1673 BlockMatrixBase<MatrixType>::set(const std::vector<size_type> &indices,
1674  const FullMatrix<number> & values,
1675  const bool elide_zero_values)
1676 {
1677  Assert(indices.size() == values.m(),
1678  ExcDimensionMismatch(indices.size(), values.m()));
1679  Assert(values.n() == values.m(), ExcNotQuadratic());
1680 
1681  for (size_type i = 0; i < indices.size(); ++i)
1682  set(indices[i],
1683  indices.size(),
1684  indices.data(),
1685  &values(i, 0),
1686  elide_zero_values);
1687 }
1688 
1689 
1690 
1691 template <class MatrixType>
1692 template <typename number>
1693 inline void
1694 BlockMatrixBase<MatrixType>::set(const size_type row,
1695  const std::vector<size_type> &col_indices,
1696  const std::vector<number> & values,
1697  const bool elide_zero_values)
1698 {
1699  Assert(col_indices.size() == values.size(),
1700  ExcDimensionMismatch(col_indices.size(), values.size()));
1701 
1702  set(row,
1703  col_indices.size(),
1704  col_indices.data(),
1705  values.data(),
1706  elide_zero_values);
1707 }
1708 
1709 
1710 
1711 // This is a very messy function, since
1712 // we need to calculate to each position
1713 // the location in the global array.
1714 template <class MatrixType>
1715 template <typename number>
1716 inline void
1717 BlockMatrixBase<MatrixType>::set(const size_type row,
1718  const size_type n_cols,
1719  const size_type *col_indices,
1720  const number * values,
1721  const bool elide_zero_values)
1722 {
1724 
1725  // lock access to the temporary data structure to
1726  // allow multiple threads to call this function concurrently
1727  std::lock_guard<std::mutex> lock(temporary_data.mutex);
1728 
1729  // Resize scratch arrays
1730  if (temporary_data.column_indices.size() < this->n_block_cols())
1731  {
1732  temporary_data.column_indices.resize(this->n_block_cols());
1733  temporary_data.column_values.resize(this->n_block_cols());
1735  }
1736 
1737  // Resize sub-arrays to n_cols. This
1738  // is a bit wasteful, but we resize
1739  // only a few times (then the maximum
1740  // row length won't increase that
1741  // much any more). At least we know
1742  // that all arrays are going to be of
1743  // the same size, so we can check
1744  // whether the size of one is large
1745  // enough before actually going
1746  // through all of them.
1747  if (temporary_data.column_indices[0].size() < n_cols)
1748  {
1749  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1750  {
1751  temporary_data.column_indices[i].resize(n_cols);
1752  temporary_data.column_values[i].resize(n_cols);
1753  }
1754  }
1755 
1756  // Reset the number of added elements
1757  // in each block to zero.
1758  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1760 
1761  // Go through the column indices to
1762  // find out which portions of the
1763  // values should be set in which
1764  // block of the matrix. We need to
1765  // touch all the data, since we can't
1766  // be sure that the data of one block
1767  // is stored contiguously (in fact,
1768  // indices will be intermixed when it
1769  // comes from an element matrix).
1770  for (size_type j = 0; j < n_cols; ++j)
1771  {
1772  number value = values[j];
1773 
1774  if (value == number() && elide_zero_values == true)
1775  continue;
1776 
1777  const std::pair<unsigned int, size_type> col_index =
1778  this->column_block_indices.global_to_local(col_indices[j]);
1779 
1780  const size_type local_index =
1781  temporary_data.counter_within_block[col_index.first]++;
1782 
1783  temporary_data.column_indices[col_index.first][local_index] =
1784  col_index.second;
1785  temporary_data.column_values[col_index.first][local_index] = value;
1786  }
1787 
1788 # ifdef DEBUG
1789  // If in debug mode, do a check whether
1790  // the right length has been obtained.
1791  size_type length = 0;
1792  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1793  length += temporary_data.counter_within_block[i];
1794  Assert(length <= n_cols, ExcInternalError());
1795 # endif
1796 
1797  // Now we found out about where the
1798  // individual columns should start and
1799  // where we should start reading out
1800  // data. Now let's write the data into
1801  // the individual blocks!
1802  const std::pair<unsigned int, size_type> row_index =
1804  for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
1805  {
1806  if (temporary_data.counter_within_block[block_col] == 0)
1807  continue;
1808 
1809  block(row_index.first, block_col)
1810  .set(row_index.second,
1812  temporary_data.column_indices[block_col].data(),
1813  temporary_data.column_values[block_col].data(),
1814  false);
1815  }
1816 }
1817 
1818 
1819 
1820 template <class MatrixType>
1821 inline void
1822 BlockMatrixBase<MatrixType>::add(const size_type i,
1823  const size_type j,
1824  const value_type value)
1825 {
1826  AssertIsFinite(value);
1827 
1829 
1830  // save some cycles for zero additions, but
1831  // only if it is safe for the matrix we are
1832  // working with
1833  using MatrixTraits = typename MatrixType::Traits;
1834  if ((MatrixTraits::zero_addition_can_be_elided == true) &&
1835  (value == value_type()))
1836  return;
1837 
1838  const std::pair<unsigned int, size_type>
1839  row_index = row_block_indices.global_to_local(i),
1840  col_index = column_block_indices.global_to_local(j);
1841  block(row_index.first, col_index.first)
1842  .add(row_index.second, col_index.second, value);
1843 }
1844 
1845 
1846 
1847 template <class MatrixType>
1848 template <typename number>
1849 inline void
1850 BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &row_indices,
1851  const std::vector<size_type> &col_indices,
1852  const FullMatrix<number> & values,
1853  const bool elide_zero_values)
1854 {
1855  Assert(row_indices.size() == values.m(),
1856  ExcDimensionMismatch(row_indices.size(), values.m()));
1857  Assert(col_indices.size() == values.n(),
1858  ExcDimensionMismatch(col_indices.size(), values.n()));
1859 
1860  for (size_type i = 0; i < row_indices.size(); ++i)
1861  add(row_indices[i],
1862  col_indices.size(),
1863  col_indices.data(),
1864  &values(i, 0),
1865  elide_zero_values);
1866 }
1867 
1868 
1869 
1870 template <class MatrixType>
1871 template <typename number>
1872 inline void
1873 BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &indices,
1874  const FullMatrix<number> & values,
1875  const bool elide_zero_values)
1876 {
1877  Assert(indices.size() == values.m(),
1878  ExcDimensionMismatch(indices.size(), values.m()));
1879  Assert(values.n() == values.m(), ExcNotQuadratic());
1880 
1881  for (size_type i = 0; i < indices.size(); ++i)
1882  add(indices[i],
1883  indices.size(),
1884  indices.data(),
1885  &values(i, 0),
1886  elide_zero_values);
1887 }
1888 
1889 
1890 
1891 template <class MatrixType>
1892 template <typename number>
1893 inline void
1894 BlockMatrixBase<MatrixType>::add(const size_type row,
1895  const std::vector<size_type> &col_indices,
1896  const std::vector<number> & values,
1897  const bool elide_zero_values)
1898 {
1899  Assert(col_indices.size() == values.size(),
1900  ExcDimensionMismatch(col_indices.size(), values.size()));
1901 
1902  add(row,
1903  col_indices.size(),
1904  col_indices.data(),
1905  values.data(),
1906  elide_zero_values);
1907 }
1908 
1909 
1910 
1911 // This is a very messy function, since
1912 // we need to calculate to each position
1913 // the location in the global array.
1914 template <class MatrixType>
1915 template <typename number>
1916 inline void
1917 BlockMatrixBase<MatrixType>::add(const size_type row,
1918  const size_type n_cols,
1919  const size_type *col_indices,
1920  const number * values,
1921  const bool elide_zero_values,
1922  const bool col_indices_are_sorted)
1923 {
1925 
1926  // TODO: Look over this to find out
1927  // whether we can do that more
1928  // efficiently.
1929  if (col_indices_are_sorted == true)
1930  {
1931 # ifdef DEBUG
1932  // check whether indices really are
1933  // sorted.
1934  size_type before = col_indices[0];
1935  for (size_type i = 1; i < n_cols; ++i)
1936  if (col_indices[i] <= before)
1937  Assert(false,
1938  ExcMessage("Flag col_indices_are_sorted is set, but "
1939  "indices appear to not be sorted.")) else before =
1940  col_indices[i];
1941 # endif
1942  const std::pair<unsigned int, size_type> row_index =
1944 
1945  if (this->n_block_cols() > 1)
1946  {
1947  const size_type *first_block =
1948  Utilities::lower_bound(col_indices,
1949  col_indices + n_cols,
1950  this->column_block_indices.block_start(1));
1951 
1952  const size_type n_zero_block_indices = first_block - col_indices;
1953  block(row_index.first, 0)
1954  .add(row_index.second,
1955  n_zero_block_indices,
1956  col_indices,
1957  values,
1958  elide_zero_values,
1959  col_indices_are_sorted);
1960 
1961  if (n_zero_block_indices < n_cols)
1962  this->add(row,
1963  n_cols - n_zero_block_indices,
1964  first_block,
1965  values + n_zero_block_indices,
1966  elide_zero_values,
1967  false);
1968  }
1969  else
1970  {
1971  block(row_index.first, 0)
1972  .add(row_index.second,
1973  n_cols,
1974  col_indices,
1975  values,
1976  elide_zero_values,
1977  col_indices_are_sorted);
1978  }
1979 
1980  return;
1981  }
1982 
1983  // Lock scratch arrays, then resize them
1984  std::lock_guard<std::mutex> lock(temporary_data.mutex);
1985 
1986  if (temporary_data.column_indices.size() < this->n_block_cols())
1987  {
1988  temporary_data.column_indices.resize(this->n_block_cols());
1989  temporary_data.column_values.resize(this->n_block_cols());
1991  }
1992 
1993  // Resize sub-arrays to n_cols. This
1994  // is a bit wasteful, but we resize
1995  // only a few times (then the maximum
1996  // row length won't increase that
1997  // much any more). At least we know
1998  // that all arrays are going to be of
1999  // the same size, so we can check
2000  // whether the size of one is large
2001  // enough before actually going
2002  // through all of them.
2003  if (temporary_data.column_indices[0].size() < n_cols)
2004  {
2005  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2006  {
2007  temporary_data.column_indices[i].resize(n_cols);
2008  temporary_data.column_values[i].resize(n_cols);
2009  }
2010  }
2011 
2012  // Reset the number of added elements
2013  // in each block to zero.
2014  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2016 
2017  // Go through the column indices to
2018  // find out which portions of the
2019  // values should be written into
2020  // which block of the matrix. We need
2021  // to touch all the data, since we
2022  // can't be sure that the data of one
2023  // block is stored contiguously (in
2024  // fact, data will be intermixed when
2025  // it comes from an element matrix).
2026  for (size_type j = 0; j < n_cols; ++j)
2027  {
2028  number value = values[j];
2029 
2030  if (value == number() && elide_zero_values == true)
2031  continue;
2032 
2033  const std::pair<unsigned int, size_type> col_index =
2034  this->column_block_indices.global_to_local(col_indices[j]);
2035 
2036  const size_type local_index =
2037  temporary_data.counter_within_block[col_index.first]++;
2038 
2039  temporary_data.column_indices[col_index.first][local_index] =
2040  col_index.second;
2041  temporary_data.column_values[col_index.first][local_index] = value;
2042  }
2043 
2044 # ifdef DEBUG
2045  // If in debug mode, do a check whether
2046  // the right length has been obtained.
2047  size_type length = 0;
2048  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2049  length += temporary_data.counter_within_block[i];
2050  Assert(length <= n_cols, ExcInternalError());
2051 # endif
2052 
2053  // Now we found out about where the
2054  // individual columns should start and
2055  // where we should start reading out
2056  // data. Now let's write the data into
2057  // the individual blocks!
2058  const std::pair<unsigned int, size_type> row_index =
2060  for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
2061  {
2062  if (temporary_data.counter_within_block[block_col] == 0)
2063  continue;
2064 
2065  block(row_index.first, block_col)
2066  .add(row_index.second,
2068  temporary_data.column_indices[block_col].data(),
2069  temporary_data.column_values[block_col].data(),
2070  false,
2071  col_indices_are_sorted);
2072  }
2073 }
2074 
2075 
2076 
2077 template <class MatrixType>
2078 inline void
2080  const BlockMatrixBase<MatrixType> &matrix)
2081 {
2082  AssertIsFinite(factor);
2083 
2085 
2086  // save some cycles for zero additions, but
2087  // only if it is safe for the matrix we are
2088  // working with
2089  using MatrixTraits = typename MatrixType::Traits;
2090  if ((MatrixTraits::zero_addition_can_be_elided == true) && (factor == 0))
2091  return;
2092 
2093  for (unsigned int row = 0; row < n_block_rows(); ++row)
2094  for (unsigned int col = 0; col < n_block_cols(); ++col)
2095  // This function should throw if the sparsity
2096  // patterns of the two blocks differ
2097  block(row, col).add(factor, matrix.block(row, col));
2098 }
2099 
2100 
2101 
2102 template <class MatrixType>
2104 BlockMatrixBase<MatrixType>::operator()(const size_type i,
2105  const size_type j) const
2106 {
2107  const std::pair<unsigned int, size_type>
2108  row_index = row_block_indices.global_to_local(i),
2109  col_index = column_block_indices.global_to_local(j);
2110  return block(row_index.first, col_index.first)(row_index.second,
2111  col_index.second);
2112 }
2113 
2114 
2115 
2116 template <class MatrixType>
2118 BlockMatrixBase<MatrixType>::el(const size_type i, const size_type j) const
2119 {
2120  const std::pair<unsigned int, size_type>
2121  row_index = row_block_indices.global_to_local(i),
2122  col_index = column_block_indices.global_to_local(j);
2123  return block(row_index.first, col_index.first)
2124  .el(row_index.second, col_index.second);
2125 }
2126 
2127 
2128 
2129 template <class MatrixType>
2131 BlockMatrixBase<MatrixType>::diag_element(const size_type i) const
2132 {
2134 
2135  const std::pair<unsigned int, size_type> index =
2137  return block(index.first, index.first).diag_element(index.second);
2138 }
2139 
2140 
2141 
2142 template <class MatrixType>
2143 inline void
2145  ::VectorOperation::values operation)
2146 {
2147  for (unsigned int r = 0; r < n_block_rows(); ++r)
2148  for (unsigned int c = 0; c < n_block_cols(); ++c)
2149  block(r, c).compress(operation);
2150 }
2151 
2152 
2153 
2154 template <class MatrixType>
2157 {
2160 
2161  for (unsigned int r = 0; r < n_block_rows(); ++r)
2162  for (unsigned int c = 0; c < n_block_cols(); ++c)
2163  block(r, c) *= factor;
2164 
2165  return *this;
2166 }
2167 
2168 
2169 
2170 template <class MatrixType>
2173 {
2176  Assert(factor != 0, ExcDivideByZero());
2177 
2178  const value_type factor_inv = 1. / factor;
2179 
2180  for (unsigned int r = 0; r < n_block_rows(); ++r)
2181  for (unsigned int c = 0; c < n_block_cols(); ++c)
2182  block(r, c) *= factor_inv;
2183 
2184  return *this;
2185 }
2186 
2187 
2188 
2189 template <class MatrixType>
2190 const BlockIndices &
2192 {
2193  return this->row_block_indices;
2194 }
2195 
2196 
2197 
2198 template <class MatrixType>
2199 const BlockIndices &
2201 {
2202  return this->column_block_indices;
2203 }
2204 
2205 
2206 
2207 template <class MatrixType>
2208 template <class BlockVectorType>
2209 void
2211  const BlockVectorType &src) const
2212 {
2213  Assert(dst.n_blocks() == n_block_rows(),
2214  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2215  Assert(src.n_blocks() == n_block_cols(),
2216  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2217 
2218  for (size_type row = 0; row < n_block_rows(); ++row)
2219  {
2220  block(row, 0).vmult(dst.block(row), src.block(0));
2221  for (size_type col = 1; col < n_block_cols(); ++col)
2222  block(row, col).vmult_add(dst.block(row), src.block(col));
2223  };
2224 }
2225 
2226 
2227 
2228 template <class MatrixType>
2229 template <class BlockVectorType, class VectorType>
2230 void
2232  VectorType & dst,
2233  const BlockVectorType &src) const
2234 {
2236  Assert(src.n_blocks() == n_block_cols(),
2237  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2238 
2239  block(0, 0).vmult(dst, src.block(0));
2240  for (size_type col = 1; col < n_block_cols(); ++col)
2241  block(0, col).vmult_add(dst, src.block(col));
2242 }
2243 
2244 
2245 
2246 template <class MatrixType>
2247 template <class BlockVectorType, class VectorType>
2248 void
2250  const VectorType &src) const
2251 {
2252  Assert(dst.n_blocks() == n_block_rows(),
2253  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2255 
2256  for (size_type row = 0; row < n_block_rows(); ++row)
2257  block(row, 0).vmult(dst.block(row), src);
2258 }
2259 
2260 
2261 
2262 template <class MatrixType>
2263 template <class VectorType>
2264 void
2266  VectorType & dst,
2267  const VectorType &src) const
2268 {
2271 
2272  block(0, 0).vmult(dst, src);
2273 }
2274 
2275 
2276 
2277 template <class MatrixType>
2278 template <class BlockVectorType>
2279 void
2280 BlockMatrixBase<MatrixType>::vmult_add(BlockVectorType & dst,
2281  const BlockVectorType &src) const
2282 {
2283  Assert(dst.n_blocks() == n_block_rows(),
2284  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2285  Assert(src.n_blocks() == n_block_cols(),
2286  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2287 
2288  for (unsigned int row = 0; row < n_block_rows(); ++row)
2289  for (unsigned int col = 0; col < n_block_cols(); ++col)
2290  block(row, col).vmult_add(dst.block(row), src.block(col));
2291 }
2292 
2293 
2294 
2295 template <class MatrixType>
2296 template <class BlockVectorType>
2297 void
2299  BlockVectorType & dst,
2300  const BlockVectorType &src) const
2301 {
2302  Assert(dst.n_blocks() == n_block_cols(),
2303  ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2304  Assert(src.n_blocks() == n_block_rows(),
2305  ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2306 
2307  dst = 0.;
2308 
2309  for (unsigned int row = 0; row < n_block_rows(); ++row)
2310  {
2311  for (unsigned int col = 0; col < n_block_cols(); ++col)
2312  block(row, col).Tvmult_add(dst.block(col), src.block(row));
2313  };
2314 }
2315 
2316 
2317 
2318 template <class MatrixType>
2319 template <class BlockVectorType, class VectorType>
2320 void
2322  const VectorType &src) const
2323 {
2324  Assert(dst.n_blocks() == n_block_cols(),
2325  ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2327 
2328  dst = 0.;
2329 
2330  for (unsigned int col = 0; col < n_block_cols(); ++col)
2331  block(0, col).Tvmult_add(dst.block(col), src);
2332 }
2333 
2334 
2335 
2336 template <class MatrixType>
2337 template <class BlockVectorType, class VectorType>
2338 void
2340  VectorType & dst,
2341  const BlockVectorType &src) const
2342 {
2344  Assert(src.n_blocks() == n_block_rows(),
2345  ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2346 
2347  block(0, 0).Tvmult(dst, src.block(0));
2348 
2349  for (size_type row = 1; row < n_block_rows(); ++row)
2350  block(row, 0).Tvmult_add(dst, src.block(row));
2351 }
2352 
2353 
2354 
2355 template <class MatrixType>
2356 template <class VectorType>
2357 void
2359  VectorType & dst,
2360  const VectorType &src) const
2361 {
2364 
2365  block(0, 0).Tvmult(dst, src);
2366 }
2367 
2368 
2369 
2370 template <class MatrixType>
2371 template <class BlockVectorType>
2372 void
2373 BlockMatrixBase<MatrixType>::Tvmult_add(BlockVectorType & dst,
2374  const BlockVectorType &src) const
2375 {
2376  Assert(dst.n_blocks() == n_block_cols(),
2377  ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2378  Assert(src.n_blocks() == n_block_rows(),
2379  ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2380 
2381  for (unsigned int row = 0; row < n_block_rows(); ++row)
2382  for (unsigned int col = 0; col < n_block_cols(); ++col)
2383  block(row, col).Tvmult_add(dst.block(col), src.block(row));
2384 }
2385 
2386 
2387 
2388 template <class MatrixType>
2389 template <class BlockVectorType>
2391 BlockMatrixBase<MatrixType>::matrix_norm_square(const BlockVectorType &v) const
2392 {
2394  Assert(v.n_blocks() == n_block_rows(),
2395  ExcDimensionMismatch(v.n_blocks(), n_block_rows()));
2396 
2397  value_type norm_sqr = 0;
2398  for (unsigned int row = 0; row < n_block_rows(); ++row)
2399  for (unsigned int col = 0; col < n_block_cols(); ++col)
2400  if (row == col)
2401  norm_sqr += block(row, col).matrix_norm_square(v.block(row));
2402  else
2403  norm_sqr +=
2404  block(row, col).matrix_scalar_product(v.block(row), v.block(col));
2405  return norm_sqr;
2406 }
2407 
2408 
2409 
2410 template <class MatrixType>
2411 template <class BlockVectorType>
2414  const BlockVectorType &u,
2415  const BlockVectorType &v) const
2416 {
2417  Assert(u.n_blocks() == n_block_rows(),
2418  ExcDimensionMismatch(u.n_blocks(), n_block_rows()));
2419  Assert(v.n_blocks() == n_block_cols(),
2420  ExcDimensionMismatch(v.n_blocks(), n_block_cols()));
2421 
2422  value_type result = 0;
2423  for (unsigned int row = 0; row < n_block_rows(); ++row)
2424  for (unsigned int col = 0; col < n_block_cols(); ++col)
2425  result +=
2426  block(row, col).matrix_scalar_product(u.block(row), v.block(col));
2427  return result;
2428 }
2429 
2430 
2431 
2432 template <class MatrixType>
2433 template <class BlockVectorType>
2435 BlockMatrixBase<MatrixType>::residual(BlockVectorType & dst,
2436  const BlockVectorType &x,
2437  const BlockVectorType &b) const
2438 {
2439  Assert(dst.n_blocks() == n_block_rows(),
2440  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2441  Assert(b.n_blocks() == n_block_rows(),
2442  ExcDimensionMismatch(b.n_blocks(), n_block_rows()));
2443  Assert(x.n_blocks() == n_block_cols(),
2444  ExcDimensionMismatch(x.n_blocks(), n_block_cols()));
2445  // in block notation, the residual is
2446  // r_i = b_i - \sum_j A_ij x_j.
2447  // this can be written as
2448  // r_i = b_i - A_i0 x_0 - \sum_{j>0} A_ij x_j.
2449  //
2450  // for the first two terms, we can
2451  // call the residual function of
2452  // A_i0. for the other terms, we
2453  // use vmult_add. however, we want
2454  // to subtract, so in order to
2455  // avoid a temporary vector, we
2456  // perform a sign change of the
2457  // first two term before, and after
2458  // adding up
2459  for (unsigned int row = 0; row < n_block_rows(); ++row)
2460  {
2461  block(row, 0).residual(dst.block(row), x.block(0), b.block(row));
2462 
2463  for (size_type i = 0; i < dst.block(row).size(); ++i)
2464  dst.block(row)(i) = -dst.block(row)(i);
2465 
2466  for (unsigned int col = 1; col < n_block_cols(); ++col)
2467  block(row, col).vmult_add(dst.block(row), x.block(col));
2468 
2469  for (size_type i = 0; i < dst.block(row).size(); ++i)
2470  dst.block(row)(i) = -dst.block(row)(i);
2471  };
2472 
2473  value_type res = 0;
2474  for (size_type row = 0; row < n_block_rows(); ++row)
2475  res += dst.block(row).norm_sqr();
2476  return std::sqrt(res);
2477 }
2478 
2479 
2480 
2481 template <class MatrixType>
2482 inline void
2483 BlockMatrixBase<MatrixType>::print(std::ostream &out,
2484  const bool alternative_output) const
2485 {
2486  for (unsigned int row = 0; row < n_block_rows(); ++row)
2487  for (unsigned int col = 0; col < n_block_cols(); ++col)
2488  {
2489  if (!alternative_output)
2490  out << "Block (" << row << ", " << col << ")" << std::endl;
2491 
2492  block(row, col).print(out, alternative_output);
2493  }
2494 }
2495 
2496 
2497 
2498 template <class MatrixType>
2501 {
2502  return const_iterator(this, 0);
2503 }
2504 
2505 
2506 
2507 template <class MatrixType>
2510 {
2511  return const_iterator(this, m());
2512 }
2513 
2514 
2515 
2516 template <class MatrixType>
2518 BlockMatrixBase<MatrixType>::begin(const size_type r) const
2519 {
2520  Assert(r < m(), ExcIndexRange(r, 0, m()));
2521  return const_iterator(this, r);
2522 }
2523 
2524 
2525 
2526 template <class MatrixType>
2528 BlockMatrixBase<MatrixType>::end(const size_type r) const
2529 {
2530  Assert(r < m(), ExcIndexRange(r, 0, m()));
2531  return const_iterator(this, r + 1);
2532 }
2533 
2534 
2535 
2536 template <class MatrixType>
2539 {
2540  return iterator(this, 0);
2541 }
2542 
2543 
2544 
2545 template <class MatrixType>
2548 {
2549  return iterator(this, m());
2550 }
2551 
2552 
2553 
2554 template <class MatrixType>
2556 BlockMatrixBase<MatrixType>::begin(const size_type r)
2557 {
2558  Assert(r < m(), ExcIndexRange(r, 0, m()));
2559  return iterator(this, r);
2560 }
2561 
2562 
2563 
2564 template <class MatrixType>
2566 BlockMatrixBase<MatrixType>::end(const size_type r)
2567 {
2568  Assert(r < m(), ExcIndexRange(r, 0, m()));
2569  return iterator(this, r + 1);
2570 }
2571 
2572 
2573 
2574 template <class MatrixType>
2575 void
2577 {
2578  std::vector<size_type> row_sizes(this->n_block_rows());
2579  std::vector<size_type> col_sizes(this->n_block_cols());
2580 
2581  // first find out the row sizes
2582  // from the first block column
2583  for (unsigned int r = 0; r < this->n_block_rows(); ++r)
2584  row_sizes[r] = sub_objects[r][0]->m();
2585  // then check that the following
2586  // block columns have the same
2587  // sizes
2588  for (unsigned int c = 1; c < this->n_block_cols(); ++c)
2589  for (unsigned int r = 0; r < this->n_block_rows(); ++r)
2590  Assert(row_sizes[r] == sub_objects[r][c]->m(),
2591  ExcIncompatibleRowNumbers(r, 0, r, c));
2592 
2593  // finally initialize the row
2594  // indices with this array
2595  this->row_block_indices.reinit(row_sizes);
2596 
2597 
2598  // then do the same with the columns
2599  for (unsigned int c = 0; c < this->n_block_cols(); ++c)
2600  col_sizes[c] = sub_objects[0][c]->n();
2601  for (unsigned int r = 1; r < this->n_block_rows(); ++r)
2602  for (unsigned int c = 0; c < this->n_block_cols(); ++c)
2603  Assert(col_sizes[c] == sub_objects[r][c]->n(),
2604  ExcIncompatibleRowNumbers(0, c, r, c));
2605 
2606  // finally initialize the row
2607  // indices with this array
2608  this->column_block_indices.reinit(col_sizes);
2609 }
2610 
2611 
2612 
2613 template <class MatrixType>
2614 void
2616 {
2617  for (unsigned int row = 0; row < n_block_rows(); ++row)
2618  for (unsigned int col = 0; col < n_block_cols(); ++col)
2619  block(row, col).prepare_add();
2620 }
2621 
2622 
2623 
2624 template <class MatrixType>
2625 void
2627 {
2628  for (unsigned int row = 0; row < n_block_rows(); ++row)
2629  for (unsigned int col = 0; col < n_block_cols(); ++col)
2630  block(row, col).prepare_set();
2631 }
2632 
2633 #endif // DOXYGEN
2634 
2635 
2636 DEAL_II_NAMESPACE_CLOSE
2637 
2638 #endif // dealii_block_matrix_base_h
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1067
size_type m() const
const types::global_dof_index invalid_size_type
Definition: types.h:182
static const unsigned int invalid_unsigned_int
Definition: types.h:173
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
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:37
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()