Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2968-g5f01c80b02 2025-03-29 13:10:00+00:00
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
block_matrix_base.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2004 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_block_matrix_base_h
16#define dealii_block_matrix_base_h
17
18
19#include <deal.II/base/config.h>
20
22#include <deal.II/base/mutex.h>
24#include <deal.II/base/table.h>
26
31#include <deal.II/lac/vector.h>
33
34#include <cmath>
35#include <mutex>
36
38
39
40// Forward declaration
41#ifndef DOXYGEN
42template <typename>
43class MatrixIterator;
44#endif
45
46
56{
61 template <typename BlockMatrixType>
63 {
64 public:
69
73 using value_type = typename BlockMatrixType::value_type;
74
79
83 unsigned int
84 block_row() const;
85
89 unsigned int
90 block_column() const;
91
92 protected:
96 unsigned int row_block;
97
101 unsigned int col_block;
102 };
103
104
105
109 template <typename BlockMatrixType, bool Constness>
110 class Accessor;
111
112
116 template <typename BlockMatrixType>
117 class Accessor<BlockMatrixType, false> : public AccessorBase<BlockMatrixType>
118 {
119 public:
124
128 using MatrixType = BlockMatrixType;
129
133 using value_type = typename BlockMatrixType::value_type;
134
143 Accessor(BlockMatrixType *m, const size_type row, const size_type col);
144
149 row() const;
150
155 column() const;
156
161 value() const;
162
166 void
167 set_value(value_type newval) const;
168
169 protected:
173 BlockMatrixType *matrix;
174
178 typename BlockMatrixType::BlockType::iterator base_iterator;
179
183 void
185
189 bool
190 operator==(const Accessor &a) const;
191
192 template <typename>
193 friend class ::MatrixIterator;
194
195 friend class Accessor<BlockMatrixType, true>;
196 };
197
198
199
204 template <typename 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
244 row() const;
245
250 column() const;
251
256 value() const;
257
258 protected:
262 const BlockMatrixType *matrix;
263
267 typename BlockMatrixType::BlockType::const_iterator base_iterator;
268
272 void
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
348template <typename MatrixType>
350{
351public:
355 using BlockType = MatrixType;
356
361 using value_type = typename BlockType::value_type;
364 using const_pointer = const value_type *;
368
369 using iterator =
371
374
375
379 BlockMatrixBase() = default;
380
385
402 template <typename 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
425 m() const;
426
432 n() const;
433
434
439 unsigned int
441
446 unsigned int
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
619 operator()(const size_type i, const size_type j) const;
620
630 el(const size_type i, const size_type j) const;
631
643 diag_element(const size_type i) const;
644
653 void
655
660 operator*=(const value_type factor);
661
666 operator/=(const value_type factor);
667
672 template <typename BlockVectorType>
673 void
674 vmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
675
681 template <typename BlockVectorType>
682 void
683 Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
684
697 template <typename BlockVectorType>
699 matrix_norm_square(const BlockVectorType &v) const;
700
707
711 template <typename BlockVectorType>
713 matrix_scalar_product(const BlockVectorType &u,
714 const BlockVectorType &v) const;
715
719 template <typename BlockVectorType>
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
739
745
750 begin(const size_type r);
751
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 &
786
790 const BlockIndices &
792
798 std::size_t
800
810 int,
811 int,
812 int,
813 int,
814 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
815 << ',' << arg4 << "] have differing row numbers.");
820 int,
821 int,
822 int,
823 int,
824 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
825 << ',' << arg4 << "] have differing column numbers.");
827protected:
840 void
842
848
853
872 void
874
885 template <typename BlockVectorType>
886 void
887 vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const;
888
899 template <typename BlockVectorType, typename VectorType>
900 void
901 vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const;
902
913 template <typename BlockVectorType, typename VectorType>
914 void
915 vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const;
916
927 template <typename VectorType>
928 void
929 vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const;
930
942 template <typename BlockVectorType>
943 void
944 Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const;
945
956 template <typename BlockVectorType, typename VectorType>
957 void
958 Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const;
959
970 template <typename BlockVectorType, typename VectorType>
971 void
972 Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const;
973
984 template <typename VectorType>
985 void
986 Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const;
987
988
989protected:
996 void
998
1003 void
1005
1006
1007private:
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
1041
1053 {
1054 return *this;
1055 }
1056 };
1057
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
1082namespace BlockMatrixIterators
1083{
1084 template <typename BlockMatrixType>
1086 : row_block(0)
1087 , col_block(0)
1088 {}
1089
1090
1091 template <typename BlockMatrixType>
1092 inline unsigned int
1093 AccessorBase<BlockMatrixType>::block_row() const
1094 {
1096
1097 return row_block;
1098 }
1099
1100
1101 template <typename BlockMatrixType>
1102 inline unsigned int
1103 AccessorBase<BlockMatrixType>::block_column() const
1104 {
1106
1107 return col_block;
1108 }
1109
1110
1111 template <typename 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 <typename 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 <typename 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 <typename BlockMatrixType>
1187 inline typename Accessor<BlockMatrixType, true>::size_type
1188 Accessor<BlockMatrixType, true>::row() const
1189 {
1190 Assert(this->row_block != numbers::invalid_unsigned_int,
1192
1193 return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1194 base_iterator->row());
1195 }
1196
1197
1198 template <typename BlockMatrixType>
1199 inline typename Accessor<BlockMatrixType, true>::size_type
1200 Accessor<BlockMatrixType, true>::column() const
1201 {
1202 Assert(this->col_block != numbers::invalid_unsigned_int,
1204
1205 return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
1206 base_iterator->column());
1207 }
1208
1209
1210 template <typename 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,
1216 Assert(this->col_block != numbers::invalid_unsigned_int,
1218
1219 return base_iterator->value();
1220 }
1221
1222
1223
1224 template <typename BlockMatrixType>
1225 inline void
1226 Accessor<BlockMatrixType, true>::advance()
1227 {
1228 Assert(this->row_block != numbers::invalid_unsigned_int,
1230 Assert(this->col_block != numbers::invalid_unsigned_int,
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 <typename 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 <typename 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 <typename BlockMatrixType>
1362 inline typename Accessor<BlockMatrixType, false>::size_type
1363 Accessor<BlockMatrixType, false>::row() const
1364 {
1366
1367 return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1368 base_iterator->row());
1369 }
1370
1371
1372 template <typename BlockMatrixType>
1373 inline typename Accessor<BlockMatrixType, false>::size_type
1374 Accessor<BlockMatrixType, false>::column() const
1375 {
1377
1378 return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
1379 base_iterator->column());
1380 }
1381
1382
1383 template <typename BlockMatrixType>
1384 inline typename Accessor<BlockMatrixType, false>::value_type
1385 Accessor<BlockMatrixType, false>::value() const
1386 {
1389
1390 return base_iterator->value();
1391 }
1392
1393
1394
1395 template <typename BlockMatrixType>
1396 inline void
1397 Accessor<BlockMatrixType, false>::set_value(
1398 typename Accessor<BlockMatrixType, false>::value_type newval) const
1399 {
1402
1403 base_iterator->value() = newval;
1404 }
1405
1406
1407
1408 template <typename BlockMatrixType>
1409 inline void
1410 Accessor<BlockMatrixType, false>::advance()
1411 {
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 <typename 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
1496template <typename MatrixType>
1498{
1499 try
1500 {
1501 clear();
1502 }
1503 catch (...)
1504 {}
1505}
1506
1507
1508template <typename MatrixType>
1509template <typename BlockMatrixType>
1511BlockMatrixBase<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
1521template <typename MatrixType>
1522std::size_t
1524{
1525 std::size_t mem =
1526 MemoryConsumption::memory_consumption(row_block_indices) +
1527 MemoryConsumption::memory_consumption(column_block_indices) +
1529 MemoryConsumption::memory_consumption(temporary_data.counter_within_block) +
1530 MemoryConsumption::memory_consumption(temporary_data.column_indices) +
1531 MemoryConsumption::memory_consumption(temporary_data.column_values) +
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
1546template <typename MatrixType>
1547inline 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
1560 row_block_indices = column_block_indices = BlockIndices();
1561}
1562
1563
1564
1565template <typename MatrixType>
1567BlockMatrixBase<MatrixType>::block(const unsigned int row,
1568 const unsigned int column)
1569{
1570 AssertIndexRange(row, n_block_rows());
1571 AssertIndexRange(column, n_block_cols());
1572
1573 return *sub_objects[row][column];
1574}
1575
1576
1577
1578template <typename MatrixType>
1579inline const typename BlockMatrixBase<MatrixType>::BlockType &
1580BlockMatrixBase<MatrixType>::block(const unsigned int row,
1581 const unsigned int column) const
1582{
1583 AssertIndexRange(row, n_block_rows());
1584 AssertIndexRange(column, n_block_cols());
1585
1586 return *sub_objects[row][column];
1587}
1588
1589
1590template <typename MatrixType>
1593{
1594 return row_block_indices.total_size();
1595}
1596
1597
1598
1599template <typename MatrixType>
1602{
1603 return column_block_indices.total_size();
1604}
1605
1606
1607
1608template <typename MatrixType>
1609inline unsigned int
1611{
1612 return column_block_indices.size();
1613}
1614
1615
1616
1617template <typename MatrixType>
1618inline 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.
1629template <typename MatrixType>
1630inline void
1631BlockMatrixBase<MatrixType>::set(const size_type i,
1632 const size_type j,
1633 const value_type value)
1634{
1635 prepare_set_operation();
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
1648template <typename MatrixType>
1649template <typename number>
1650inline void
1651BlockMatrixBase<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
1671template <typename MatrixType>
1672template <typename number>
1673inline void
1674BlockMatrixBase<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
1692template <typename MatrixType>
1693template <typename number>
1694inline void
1695BlockMatrixBase<MatrixType>::set(const size_type row,
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.
1715template <typename MatrixType>
1716template <typename number>
1717inline void
1718BlockMatrixBase<MatrixType>::set(const size_type row,
1719 const size_type n_cols,
1720 const size_type *col_indices,
1721 const number *values,
1722 const bool elide_zero_values)
1723{
1724 prepare_set_operation();
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());
1735 temporary_data.counter_within_block.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)
1760 temporary_data.counter_within_block[i] = 0;
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 if constexpr (running_in_debug_mode())
1790 {
1791 // If in debug mode, do a check whether
1792 // the right length has been obtained.
1793 size_type length = 0;
1794 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1795 length += temporary_data.counter_within_block[i];
1796 Assert(length <= n_cols, ExcInternalError());
1797 }
1798
1799 // Now we found out about where the
1800 // individual columns should start and
1801 // where we should start reading out
1802 // data. Now let's write the data into
1803 // the individual blocks!
1804 const std::pair<unsigned int, size_type> row_index =
1805 this->row_block_indices.global_to_local(row);
1806 for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
1807 {
1808 if (temporary_data.counter_within_block[block_col] == 0)
1809 continue;
1810
1811 block(row_index.first, block_col)
1812 .set(row_index.second,
1813 temporary_data.counter_within_block[block_col],
1814 temporary_data.column_indices[block_col].data(),
1815 temporary_data.column_values[block_col].data(),
1816 false);
1817 }
1818}
1819
1820
1821
1822template <typename MatrixType>
1823inline void
1824BlockMatrixBase<MatrixType>::add(const size_type i,
1825 const size_type j,
1826 const value_type value)
1827{
1828 AssertIsFinite(value);
1829
1830 prepare_add_operation();
1831
1832 // save some cycles for zero additions, but
1833 // only if it is safe for the matrix we are
1834 // working with
1835 using MatrixTraits = typename MatrixType::Traits;
1836 if ((MatrixTraits::zero_addition_can_be_elided == true) &&
1837 (value == value_type()))
1838 return;
1839
1840 const std::pair<unsigned int, size_type>
1841 row_index = row_block_indices.global_to_local(i),
1842 col_index = column_block_indices.global_to_local(j);
1843 block(row_index.first, col_index.first)
1844 .add(row_index.second, col_index.second, value);
1845}
1846
1847
1848
1849template <typename MatrixType>
1850template <typename number>
1851inline void
1852BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &row_indices,
1853 const std::vector<size_type> &col_indices,
1854 const FullMatrix<number> &values,
1855 const bool elide_zero_values)
1856{
1857 Assert(row_indices.size() == values.m(),
1858 ExcDimensionMismatch(row_indices.size(), values.m()));
1859 Assert(col_indices.size() == values.n(),
1860 ExcDimensionMismatch(col_indices.size(), values.n()));
1861
1862 for (size_type i = 0; i < row_indices.size(); ++i)
1863 add(row_indices[i],
1864 col_indices.size(),
1865 col_indices.data(),
1866 &values(i, 0),
1867 elide_zero_values);
1868}
1869
1870
1871
1872template <typename MatrixType>
1873template <typename number>
1874inline void
1875BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &indices,
1876 const FullMatrix<number> &values,
1877 const bool elide_zero_values)
1878{
1879 Assert(indices.size() == values.m(),
1880 ExcDimensionMismatch(indices.size(), values.m()));
1881 Assert(values.n() == values.m(), ExcNotQuadratic());
1882
1883 for (size_type i = 0; i < indices.size(); ++i)
1884 add(indices[i],
1885 indices.size(),
1886 indices.data(),
1887 &values(i, 0),
1888 elide_zero_values);
1889}
1890
1891
1892
1893template <typename MatrixType>
1894template <typename number>
1895inline void
1896BlockMatrixBase<MatrixType>::add(const size_type row,
1897 const std::vector<size_type> &col_indices,
1898 const std::vector<number> &values,
1899 const bool elide_zero_values)
1900{
1901 Assert(col_indices.size() == values.size(),
1902 ExcDimensionMismatch(col_indices.size(), values.size()));
1903
1904 add(row,
1905 col_indices.size(),
1906 col_indices.data(),
1907 values.data(),
1908 elide_zero_values);
1909}
1910
1911
1912
1913// This is a very messy function, since
1914// we need to calculate to each position
1915// the location in the global array.
1916template <typename MatrixType>
1917template <typename number>
1918inline void
1919BlockMatrixBase<MatrixType>::add(const size_type row,
1920 const size_type n_cols,
1921 const size_type *col_indices,
1922 const number *values,
1923 const bool elide_zero_values,
1924 const bool col_indices_are_sorted)
1925{
1926 prepare_add_operation();
1927
1928 // TODO: Look over this to find out
1929 // whether we can do that more
1930 // efficiently.
1931 if (col_indices_are_sorted == true)
1932 {
1933 if constexpr (running_in_debug_mode())
1934 {
1935 // check whether indices really are
1936 // sorted.
1937 size_type before = col_indices[0];
1938 for (size_type i = 1; i < n_cols; ++i)
1939 if (col_indices[i] <= before)
1940 {
1941 Assert(false,
1942 ExcMessage("Flag col_indices_are_sorted is set, but "
1943 "indices appear to not be sorted."));
1944 }
1945 else
1946 before = col_indices[i];
1947 }
1948 const std::pair<unsigned int, size_type> row_index =
1949 this->row_block_indices.global_to_local(row);
1950
1951 if (this->n_block_cols() > 1)
1952 {
1953 const size_type *first_block =
1954 Utilities::lower_bound(col_indices,
1955 col_indices + n_cols,
1956 this->column_block_indices.block_start(1));
1957
1958 const size_type n_zero_block_indices = first_block - col_indices;
1959 block(row_index.first, 0)
1960 .add(row_index.second,
1961 n_zero_block_indices,
1962 col_indices,
1963 values,
1964 elide_zero_values,
1965 col_indices_are_sorted);
1966
1967 if (n_zero_block_indices < n_cols)
1968 this->add(row,
1969 n_cols - n_zero_block_indices,
1970 first_block,
1971 values + n_zero_block_indices,
1972 elide_zero_values,
1973 false);
1974 }
1975 else
1976 {
1977 block(row_index.first, 0)
1978 .add(row_index.second,
1979 n_cols,
1980 col_indices,
1981 values,
1982 elide_zero_values,
1983 col_indices_are_sorted);
1984 }
1985
1986 return;
1987 }
1988
1989 // Lock scratch arrays, then resize them
1990 std::lock_guard<std::mutex> lock(temporary_data.mutex);
1991
1992 if (temporary_data.column_indices.size() < this->n_block_cols())
1993 {
1994 temporary_data.column_indices.resize(this->n_block_cols());
1995 temporary_data.column_values.resize(this->n_block_cols());
1996 temporary_data.counter_within_block.resize(this->n_block_cols());
1997 }
1998
1999 // Resize sub-arrays to n_cols. This
2000 // is a bit wasteful, but we resize
2001 // only a few times (then the maximum
2002 // row length won't increase that
2003 // much any more). At least we know
2004 // that all arrays are going to be of
2005 // the same size, so we can check
2006 // whether the size of one is large
2007 // enough before actually going
2008 // through all of them.
2009 if (temporary_data.column_indices[0].size() < n_cols)
2010 {
2011 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2012 {
2013 temporary_data.column_indices[i].resize(n_cols);
2014 temporary_data.column_values[i].resize(n_cols);
2015 }
2016 }
2017
2018 // Reset the number of added elements
2019 // in each block to zero.
2020 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2021 temporary_data.counter_within_block[i] = 0;
2022
2023 // Go through the column indices to
2024 // find out which portions of the
2025 // values should be written into
2026 // which block of the matrix. We need
2027 // to touch all the data, since we
2028 // can't be sure that the data of one
2029 // block is stored contiguously (in
2030 // fact, data will be intermixed when
2031 // it comes from an element matrix).
2032 for (size_type j = 0; j < n_cols; ++j)
2033 {
2034 number value = values[j];
2035
2036 if (value == number() && elide_zero_values == true)
2037 continue;
2038
2039 const std::pair<unsigned int, size_type> col_index =
2040 this->column_block_indices.global_to_local(col_indices[j]);
2041
2042 const size_type local_index =
2043 temporary_data.counter_within_block[col_index.first]++;
2044
2045 temporary_data.column_indices[col_index.first][local_index] =
2046 col_index.second;
2047 temporary_data.column_values[col_index.first][local_index] = value;
2048 }
2049
2050 if constexpr (running_in_debug_mode())
2051 {
2052 // If in debug mode, do a check whether
2053 // the right length has been obtained.
2054 size_type length = 0;
2055 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2056 length += temporary_data.counter_within_block[i];
2057 Assert(length <= n_cols, ExcInternalError());
2058 }
2059
2060 // Now we found out about where the
2061 // individual columns should start and
2062 // where we should start reading out
2063 // data. Now let's write the data into
2064 // the individual blocks!
2065 const std::pair<unsigned int, size_type> row_index =
2066 this->row_block_indices.global_to_local(row);
2067 for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
2068 {
2069 if (temporary_data.counter_within_block[block_col] == 0)
2070 continue;
2071
2072 block(row_index.first, block_col)
2073 .add(row_index.second,
2074 temporary_data.counter_within_block[block_col],
2075 temporary_data.column_indices[block_col].data(),
2076 temporary_data.column_values[block_col].data(),
2077 false,
2078 col_indices_are_sorted);
2079 }
2080}
2081
2082
2083
2084template <typename MatrixType>
2085inline void
2086BlockMatrixBase<MatrixType>::add(const value_type factor,
2087 const BlockMatrixBase<MatrixType> &matrix)
2088{
2089 AssertIsFinite(factor);
2090
2091 prepare_add_operation();
2092
2093 // save some cycles for zero additions, but
2094 // only if it is safe for the matrix we are
2095 // working with
2096 using MatrixTraits = typename MatrixType::Traits;
2097 if ((MatrixTraits::zero_addition_can_be_elided == true) && (factor == 0))
2098 return;
2099
2100 for (unsigned int row = 0; row < n_block_rows(); ++row)
2101 for (unsigned int col = 0; col < n_block_cols(); ++col)
2102 // This function should throw if the sparsity
2103 // patterns of the two blocks differ
2104 block(row, col).add(factor, matrix.block(row, col));
2105}
2106
2107
2108
2109template <typename MatrixType>
2112 const size_type j) const
2113{
2114 const std::pair<unsigned int, size_type>
2115 row_index = row_block_indices.global_to_local(i),
2116 col_index = column_block_indices.global_to_local(j);
2117 return block(row_index.first, col_index.first)(row_index.second,
2118 col_index.second);
2119}
2120
2121
2122
2123template <typename MatrixType>
2125BlockMatrixBase<MatrixType>::el(const size_type i, const size_type j) const
2126{
2127 const std::pair<unsigned int, size_type>
2128 row_index = row_block_indices.global_to_local(i),
2129 col_index = column_block_indices.global_to_local(j);
2130 return block(row_index.first, col_index.first)
2131 .el(row_index.second, col_index.second);
2132}
2133
2134
2135
2136template <typename MatrixType>
2138BlockMatrixBase<MatrixType>::diag_element(const size_type i) const
2139{
2140 Assert(n_block_rows() == n_block_cols(), ExcNotQuadratic());
2141
2142 const std::pair<unsigned int, size_type> index =
2143 row_block_indices.global_to_local(i);
2144 return block(index.first, index.first).diag_element(index.second);
2145}
2146
2147
2148
2149template <typename MatrixType>
2150inline void
2152{
2153 for (unsigned int r = 0; r < n_block_rows(); ++r)
2154 for (unsigned int c = 0; c < n_block_cols(); ++c)
2155 block(r, c).compress(operation);
2156}
2157
2158
2159
2160template <typename MatrixType>
2162BlockMatrixBase<MatrixType>::operator*=(const value_type factor)
2163{
2164 Assert(n_block_cols() != 0, ExcNotInitialized());
2165 Assert(n_block_rows() != 0, ExcNotInitialized());
2166
2167 for (unsigned int r = 0; r < n_block_rows(); ++r)
2168 for (unsigned int c = 0; c < n_block_cols(); ++c)
2169 block(r, c) *= factor;
2170
2171 return *this;
2172}
2173
2174
2175
2176template <typename MatrixType>
2178BlockMatrixBase<MatrixType>::operator/=(const value_type factor)
2179{
2180 Assert(n_block_cols() != 0, ExcNotInitialized());
2181 Assert(n_block_rows() != 0, ExcNotInitialized());
2182 Assert(factor != 0, ExcDivideByZero());
2183
2184 const value_type factor_inv = 1. / factor;
2185
2186 for (unsigned int r = 0; r < n_block_rows(); ++r)
2187 for (unsigned int c = 0; c < n_block_cols(); ++c)
2188 block(r, c) *= factor_inv;
2189
2190 return *this;
2191}
2192
2193
2194
2195template <typename MatrixType>
2196const BlockIndices &
2198{
2199 return this->row_block_indices;
2200}
2201
2202
2203
2204template <typename MatrixType>
2205const BlockIndices &
2207{
2208 return this->column_block_indices;
2209}
2210
2211
2212
2213template <typename MatrixType>
2214template <typename BlockVectorType>
2215void
2217 const BlockVectorType &src) const
2218{
2219 Assert(dst.n_blocks() == n_block_rows(),
2220 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2221 Assert(src.n_blocks() == n_block_cols(),
2222 ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2223
2224 for (size_type row = 0; row < n_block_rows(); ++row)
2225 {
2226 block(row, 0).vmult(dst.block(row), src.block(0));
2227 for (size_type col = 1; col < n_block_cols(); ++col)
2228 block(row, col).vmult_add(dst.block(row), src.block(col));
2229 };
2230}
2231
2232
2233
2234template <typename MatrixType>
2235template <typename BlockVectorType, typename VectorType>
2236void
2238 VectorType &dst,
2239 const BlockVectorType &src) const
2240{
2241 Assert(n_block_rows() == 1, ExcDimensionMismatch(1, n_block_rows()));
2242 Assert(src.n_blocks() == n_block_cols(),
2243 ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2244
2245 block(0, 0).vmult(dst, src.block(0));
2246 for (size_type col = 1; col < n_block_cols(); ++col)
2247 block(0, col).vmult_add(dst, src.block(col));
2248}
2249
2250
2251
2252template <typename MatrixType>
2253template <typename BlockVectorType, typename VectorType>
2254void
2256 const VectorType &src) const
2257{
2258 Assert(dst.n_blocks() == n_block_rows(),
2259 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2260 Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2261
2262 for (size_type row = 0; row < n_block_rows(); ++row)
2263 block(row, 0).vmult(dst.block(row), src);
2264}
2265
2266
2267
2268template <typename MatrixType>
2269template <typename VectorType>
2270void
2272 VectorType &dst,
2273 const VectorType &src) const
2274{
2275 Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));
2276 Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2277
2278 block(0, 0).vmult(dst, src);
2279}
2280
2281
2282
2283template <typename MatrixType>
2284template <typename BlockVectorType>
2285void
2286BlockMatrixBase<MatrixType>::vmult_add(BlockVectorType &dst,
2287 const BlockVectorType &src) const
2288{
2289 Assert(dst.n_blocks() == n_block_rows(),
2290 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2291 Assert(src.n_blocks() == n_block_cols(),
2292 ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2293
2294 for (unsigned int row = 0; row < n_block_rows(); ++row)
2295 for (unsigned int col = 0; col < n_block_cols(); ++col)
2296 block(row, col).vmult_add(dst.block(row), src.block(col));
2297}
2298
2299
2300
2301template <typename MatrixType>
2302template <typename BlockVectorType>
2303void
2305 BlockVectorType &dst,
2306 const BlockVectorType &src) const
2307{
2308 Assert(dst.n_blocks() == n_block_cols(),
2309 ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2310 Assert(src.n_blocks() == n_block_rows(),
2311 ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2312
2313 dst = 0.;
2314
2315 for (unsigned int row = 0; row < n_block_rows(); ++row)
2316 {
2317 for (unsigned int col = 0; col < n_block_cols(); ++col)
2318 block(row, col).Tvmult_add(dst.block(col), src.block(row));
2319 };
2320}
2321
2322
2323
2324template <typename MatrixType>
2325template <typename BlockVectorType, typename VectorType>
2326void
2328 const VectorType &src) const
2329{
2330 Assert(dst.n_blocks() == n_block_cols(),
2331 ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2332 Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));
2333
2334 dst = 0.;
2335
2336 for (unsigned int col = 0; col < n_block_cols(); ++col)
2337 block(0, col).Tvmult_add(dst.block(col), src);
2338}
2339
2340
2341
2342template <typename MatrixType>
2343template <typename BlockVectorType, typename VectorType>
2344void
2346 VectorType &dst,
2347 const BlockVectorType &src) const
2348{
2349 Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2350 Assert(src.n_blocks() == n_block_rows(),
2351 ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2352
2353 block(0, 0).Tvmult(dst, src.block(0));
2354
2355 for (size_type row = 1; row < n_block_rows(); ++row)
2356 block(row, 0).Tvmult_add(dst, src.block(row));
2357}
2358
2359
2360
2361template <typename MatrixType>
2362template <typename VectorType>
2363void
2365 VectorType &dst,
2366 const VectorType &src) const
2367{
2368 Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2369 Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));
2370
2371 block(0, 0).Tvmult(dst, src);
2372}
2373
2374
2375
2376template <typename MatrixType>
2377template <typename BlockVectorType>
2378void
2379BlockMatrixBase<MatrixType>::Tvmult_add(BlockVectorType &dst,
2380 const BlockVectorType &src) const
2381{
2382 Assert(dst.n_blocks() == n_block_cols(),
2383 ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2384 Assert(src.n_blocks() == n_block_rows(),
2385 ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2386
2387 for (unsigned int row = 0; row < n_block_rows(); ++row)
2388 for (unsigned int col = 0; col < n_block_cols(); ++col)
2389 block(row, col).Tvmult_add(dst.block(col), src.block(row));
2390}
2391
2392
2393
2394template <typename MatrixType>
2395template <typename BlockVectorType>
2397BlockMatrixBase<MatrixType>::matrix_norm_square(const BlockVectorType &v) const
2398{
2399 Assert(n_block_rows() == n_block_cols(), ExcNotQuadratic());
2400 Assert(v.n_blocks() == n_block_rows(),
2401 ExcDimensionMismatch(v.n_blocks(), n_block_rows()));
2402
2403 value_type norm_sqr = 0;
2404 for (unsigned int row = 0; row < n_block_rows(); ++row)
2405 for (unsigned int col = 0; col < n_block_cols(); ++col)
2406 if (row == col)
2407 norm_sqr += block(row, col).matrix_norm_square(v.block(row));
2408 else
2409 norm_sqr +=
2410 block(row, col).matrix_scalar_product(v.block(row), v.block(col));
2411 return norm_sqr;
2412}
2413
2414
2415
2416template <typename MatrixType>
2419{
2420 value_type norm_sqr = 0;
2421
2422 // For each block, get the Frobenius norm, and add the square to the
2423 // accumulator for the full matrix
2424 for (unsigned int row = 0; row < n_block_rows(); ++row)
2425 {
2426 for (unsigned int col = 0; col < n_block_cols(); ++col)
2427 {
2428 const value_type block_norm = block(row, col).frobenius_norm();
2429 norm_sqr += block_norm * block_norm;
2430 }
2431 }
2432
2433 return std::sqrt(norm_sqr);
2434}
2435
2436
2437
2438template <typename MatrixType>
2439template <typename BlockVectorType>
2442 const BlockVectorType &u,
2443 const BlockVectorType &v) const
2444{
2445 Assert(u.n_blocks() == n_block_rows(),
2446 ExcDimensionMismatch(u.n_blocks(), n_block_rows()));
2447 Assert(v.n_blocks() == n_block_cols(),
2448 ExcDimensionMismatch(v.n_blocks(), n_block_cols()));
2449
2450 value_type result = 0;
2451 for (unsigned int row = 0; row < n_block_rows(); ++row)
2452 for (unsigned int col = 0; col < n_block_cols(); ++col)
2453 result +=
2454 block(row, col).matrix_scalar_product(u.block(row), v.block(col));
2455 return result;
2456}
2457
2458
2459
2460template <typename MatrixType>
2461template <typename BlockVectorType>
2463BlockMatrixBase<MatrixType>::residual(BlockVectorType &dst,
2464 const BlockVectorType &x,
2465 const BlockVectorType &b) const
2466{
2467 Assert(dst.n_blocks() == n_block_rows(),
2468 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2469 Assert(b.n_blocks() == n_block_rows(),
2470 ExcDimensionMismatch(b.n_blocks(), n_block_rows()));
2471 Assert(x.n_blocks() == n_block_cols(),
2472 ExcDimensionMismatch(x.n_blocks(), n_block_cols()));
2473 // in block notation, the residual is
2474 // r_i = b_i - \sum_j A_ij x_j.
2475 // this can be written as
2476 // r_i = b_i - A_i0 x_0 - \sum_{j>0} A_ij x_j.
2477 //
2478 // for the first two terms, we can
2479 // call the residual function of
2480 // A_i0. for the other terms, we
2481 // use vmult_add. however, we want
2482 // to subtract, so in order to
2483 // avoid a temporary vector, we
2484 // perform a sign change of the
2485 // first two term before, and after
2486 // adding up
2487 for (unsigned int row = 0; row < n_block_rows(); ++row)
2488 {
2489 block(row, 0).residual(dst.block(row), x.block(0), b.block(row));
2490
2491 for (size_type i = 0; i < dst.block(row).size(); ++i)
2492 dst.block(row)(i) = -dst.block(row)(i);
2493
2494 for (unsigned int col = 1; col < n_block_cols(); ++col)
2495 block(row, col).vmult_add(dst.block(row), x.block(col));
2496
2497 for (size_type i = 0; i < dst.block(row).size(); ++i)
2498 dst.block(row)(i) = -dst.block(row)(i);
2499 };
2500
2501 value_type res = 0;
2502 for (size_type row = 0; row < n_block_rows(); ++row)
2503 res += dst.block(row).norm_sqr();
2504 return std::sqrt(res);
2505}
2506
2507
2508
2509template <typename MatrixType>
2510inline void
2511BlockMatrixBase<MatrixType>::print(std::ostream &out,
2512 const bool alternative_output) const
2513{
2514 for (unsigned int row = 0; row < n_block_rows(); ++row)
2515 for (unsigned int col = 0; col < n_block_cols(); ++col)
2516 {
2517 if (!alternative_output)
2518 out << "Block (" << row << ", " << col << ')' << std::endl;
2519
2520 block(row, col).print(out, alternative_output);
2521 }
2522}
2523
2524
2525
2526template <typename MatrixType>
2529{
2530 return const_iterator(this, 0);
2531}
2532
2533
2534
2535template <typename MatrixType>
2538{
2539 return const_iterator(this, m());
2540}
2541
2542
2543
2544template <typename MatrixType>
2546BlockMatrixBase<MatrixType>::begin(const size_type r) const
2547{
2548 AssertIndexRange(r, m());
2549 return const_iterator(this, r);
2550}
2551
2552
2553
2554template <typename MatrixType>
2556BlockMatrixBase<MatrixType>::end(const size_type r) const
2557{
2558 AssertIndexRange(r, m());
2559 return const_iterator(this, r + 1);
2560}
2561
2562
2563
2564template <typename MatrixType>
2567{
2568 return iterator(this, 0);
2569}
2570
2571
2572
2573template <typename MatrixType>
2576{
2577 return iterator(this, m());
2578}
2579
2580
2581
2582template <typename MatrixType>
2584BlockMatrixBase<MatrixType>::begin(const size_type r)
2585{
2586 AssertIndexRange(r, m());
2587 return iterator(this, r);
2588}
2589
2590
2591
2592template <typename MatrixType>
2594BlockMatrixBase<MatrixType>::end(const size_type r)
2595{
2596 AssertIndexRange(r, m());
2597 return iterator(this, r + 1);
2598}
2599
2600
2601
2602template <typename MatrixType>
2603void
2605{
2606 std::vector<size_type> row_sizes(this->n_block_rows());
2607 std::vector<size_type> col_sizes(this->n_block_cols());
2608
2609 // first find out the row sizes
2610 // from the first block column
2611 for (unsigned int r = 0; r < this->n_block_rows(); ++r)
2612 row_sizes[r] = sub_objects[r][0]->m();
2613 // then check that the following
2614 // block columns have the same
2615 // sizes
2616 for (unsigned int c = 1; c < this->n_block_cols(); ++c)
2617 for (unsigned int r = 0; r < this->n_block_rows(); ++r)
2618 Assert(row_sizes[r] == sub_objects[r][c]->m(),
2619 ExcIncompatibleRowNumbers(r, 0, r, c));
2620
2621 // finally initialize the row
2622 // indices with this array
2623 this->row_block_indices.reinit(row_sizes);
2624
2625
2626 // then do the same with the columns
2627 for (unsigned int c = 0; c < this->n_block_cols(); ++c)
2628 col_sizes[c] = sub_objects[0][c]->n();
2629 for (unsigned int r = 1; r < this->n_block_rows(); ++r)
2630 for (unsigned int c = 0; c < this->n_block_cols(); ++c)
2631 Assert(col_sizes[c] == sub_objects[r][c]->n(),
2632 ExcIncompatibleRowNumbers(0, c, r, c));
2633
2634 // finally initialize the row
2635 // indices with this array
2636 this->column_block_indices.reinit(col_sizes);
2637}
2638
2639
2640
2641template <typename MatrixType>
2642void
2644{
2645 for (unsigned int row = 0; row < n_block_rows(); ++row)
2646 for (unsigned int col = 0; col < n_block_cols(); ++col)
2647 block(row, col).prepare_add();
2648}
2649
2650
2651
2652template <typename MatrixType>
2653void
2655{
2656 for (unsigned int row = 0; row < n_block_rows(); ++row)
2657 for (unsigned int col = 0; col < n_block_cols(); ++col)
2658 block(row, col).prepare_set();
2659}
2660
2661#endif // DOXYGEN
2662
2663
2665
2666#endif // dealii_block_matrix_base_h
void Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
void vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
void add(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number > &values, const bool elide_zero_values=true)
void add(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=true)
value_type & reference
void add(const size_type row, const size_type n_cols, const size_type *col_indices, const number *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
void print(std::ostream &out, const bool alternative_output=false) const
friend class BlockMatrixIterators::Accessor
const_iterator end() const
BlockIndices column_block_indices
BlockMatrixBase & operator*=(const value_type factor)
value_type matrix_norm_square(const BlockVectorType &v) const
void Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
const value_type & const_reference
void Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
unsigned int n_block_rows() const
value_type operator()(const size_type i, const size_type j) const
value_type el(const size_type i, const size_type j) const
const_iterator begin() const
real_type frobenius_norm() const
void vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
void compress(VectorOperation::values operation)
void vmult_add(BlockVectorType &dst, const BlockVectorType &src) const
types::global_dof_index size_type
void Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const
typename BlockType::value_type value_type
void collect_sizes()
void prepare_add_operation()
void vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
value_type residual(BlockVectorType &dst, const BlockVectorType &x, const BlockVectorType &b) const
BlockMatrixBase()=default
void set(const size_type row, const size_type n_cols, const size_type *col_indices, const number *values, const bool elide_zero_values=false)
size_type n() const
void prepare_set_operation()
unsigned int n_block_cols() const
const value_type * const_pointer
iterator begin()
TemporaryData temporary_data
const BlockIndices & get_row_indices() const
iterator end()
void set(const size_type i, const size_type j, const value_type value)
BlockMatrixBase & copy_from(const BlockMatrixType &source)
std::size_t memory_consumption() const
void add(const value_type factor, const BlockMatrixBase< MatrixType > &matrix)
const BlockType & block(const unsigned int row, const unsigned int column) const
typename numbers::NumberTraits< value_type >::real_type real_type
BlockIndices row_block_indices
~BlockMatrixBase() override
BlockType & block(const unsigned int row, const unsigned int column)
size_type m() const
value_type matrix_scalar_product(const BlockVectorType &u, const BlockVectorType &v) const
const BlockIndices & get_column_indices() const
void set(const std::vector< size_type > &indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=false)
void add(const std::vector< size_type > &indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=true)
void add(const size_type i, const size_type j, const value_type value)
void set(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number > &values, const bool elide_zero_values=false)
value_type diag_element(const size_type i) const
void vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
Table< 2, ObserverPointer< BlockType, BlockMatrixBase< MatrixType > > > sub_objects
iterator begin(const size_type r)
const_iterator begin(const size_type r) const
BlockMatrixBase & operator/=(const value_type factor)
void Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
iterator end(const size_type r)
const_iterator end(const size_type r) const
void set(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=false)
unsigned int block_row() const
unsigned int block_column() const
typename BlockMatrixType::value_type value_type
Accessor(BlockMatrixType *m, const size_type row, const size_type col)
Accessor(const Accessor< BlockMatrixType, false > &)
BlockMatrixType::BlockType::const_iterator base_iterator
Accessor(const BlockMatrixType *m, const size_type row, const size_type col)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
constexpr bool running_in_debug_mode()
Definition config.h:78
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
static ::ExceptionBase & ExcIncompatibleColNumbers(int arg1, int arg2, int arg3, int arg4)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcIncompatibleRowNumbers(int arg1, int arg2, int arg3, int arg4)
static ::ExceptionBase & ExcIteratorPastEnd()
#define AssertIsFinite(number)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcMessage(std::string arg1)
std::size_t size
Definition mpi.cc:745
@ matrix
Contents is actually a matrix.
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition utilities.h:1033
constexpr types::global_dof_index invalid_size_type
Definition types.h:244
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition types.h:90
std::vector< std::vector< size_type > > column_indices
std::vector< std::vector< value_type > > column_values
TemporaryData & operator=(const TemporaryData &)
std::vector< size_type > counter_within_block