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