Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
block_matrix_base.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2004 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_block_matrix_base_h
17#define dealii_block_matrix_base_h
18
19
20#include <deal.II/base/config.h>
21
23#include <deal.II/base/mutex.h>
25#include <deal.II/base/table.h>
27
32#include <deal.II/lac/vector.h>
34
35#include <cmath>
36#include <mutex>
37
39
40
41// Forward declaration
42#ifndef DOXYGEN
43template <typename>
44class MatrixIterator;
45#endif
46
47
57{
62 template <class BlockMatrixType>
64 {
65 public:
70
74 using value_type = typename BlockMatrixType::value_type;
75
80
84 unsigned int
85 block_row() const;
86
90 unsigned int
91 block_column() const;
92
93 protected:
97 unsigned int row_block;
98
102 unsigned int col_block;
103
104 // Let the iterator class be a friend.
105 template <typename>
106 friend class MatrixIterator;
107 };
108
109
110
114 template <class BlockMatrixType, bool Constness>
115 class Accessor;
116
117
121 template <class BlockMatrixType>
122 class Accessor<BlockMatrixType, false> : public AccessorBase<BlockMatrixType>
123 {
124 public:
129
133 using MatrixType = BlockMatrixType;
134
138 using value_type = typename BlockMatrixType::value_type;
139
148 Accessor(BlockMatrixType *m, const size_type row, const size_type col);
149
154 row() const;
155
160 column() const;
161
166 value() const;
167
171 void
172 set_value(value_type newval) const;
173
174 protected:
178 BlockMatrixType *matrix;
179
183 typename BlockMatrixType::BlockType::iterator base_iterator;
184
188 void
190
194 bool
195 operator==(const Accessor &a) const;
196
197 template <typename>
198 friend class MatrixIterator;
199 friend class Accessor<BlockMatrixType, true>;
200 };
201
206 template <class BlockMatrixType>
207 class Accessor<BlockMatrixType, true> : public AccessorBase<BlockMatrixType>
208 {
209 public:
214
218 using MatrixType = const BlockMatrixType;
219
223 using value_type = typename BlockMatrixType::value_type;
224
233 Accessor(const BlockMatrixType *m,
234 const size_type row,
235 const size_type col);
236
241
246 row() const;
247
252 column() const;
253
258 value() const;
259
260 protected:
264 const BlockMatrixType *matrix;
265
269 typename BlockMatrixType::BlockType::const_iterator base_iterator;
270
274 void
276
280 bool
281 operator==(const Accessor &a) const;
282
283 // Let the iterator class be a friend.
284 template <typename>
285 friend class ::MatrixIterator;
286 };
287} // namespace BlockMatrixIterators
288
289
290
350template <typename MatrixType>
352{
353public:
357 using BlockType = MatrixType;
358
363 using value_type = typename BlockType::value_type;
366 using const_pointer = const value_type *;
370
371 using iterator =
373
376
377
381 BlockMatrixBase() = default;
382
387
404 template <class BlockMatrixType>
406 copy_from(const BlockMatrixType &source);
407
411 BlockType &
412 block(const unsigned int row, const unsigned int column);
413
414
419 const BlockType &
420 block(const unsigned int row, const unsigned int column) const;
421
427 m() const;
428
434 n() const;
435
436
441 unsigned int
443
448 unsigned int
450
456 void
457 set(const size_type i, const size_type j, const value_type value);
458
474 template <typename number>
475 void
476 set(const std::vector<size_type> &indices,
477 const FullMatrix<number> & full_matrix,
478 const bool elide_zero_values = false);
479
485 template <typename number>
486 void
487 set(const std::vector<size_type> &row_indices,
488 const std::vector<size_type> &col_indices,
489 const FullMatrix<number> & full_matrix,
490 const bool elide_zero_values = false);
491
502 template <typename number>
503 void
504 set(const size_type row,
505 const std::vector<size_type> &col_indices,
506 const std::vector<number> & values,
507 const bool elide_zero_values = false);
508
518 template <typename number>
519 void
520 set(const size_type row,
521 const size_type n_cols,
522 const size_type *col_indices,
523 const number * values,
524 const bool elide_zero_values = false);
525
531 void
532 add(const size_type i, const size_type j, const value_type value);
533
548 template <typename number>
549 void
550 add(const std::vector<size_type> &indices,
551 const FullMatrix<number> & full_matrix,
552 const bool elide_zero_values = true);
553
559 template <typename number>
560 void
561 add(const std::vector<size_type> &row_indices,
562 const std::vector<size_type> &col_indices,
563 const FullMatrix<number> & full_matrix,
564 const bool elide_zero_values = true);
565
575 template <typename number>
576 void
577 add(const size_type row,
578 const std::vector<size_type> &col_indices,
579 const std::vector<number> & values,
580 const bool elide_zero_values = true);
581
591 template <typename number>
592 void
593 add(const size_type row,
594 const size_type n_cols,
595 const size_type *col_indices,
596 const number * values,
597 const bool elide_zero_values = true,
598 const bool col_indices_are_sorted = false);
599
611 void
612 add(const value_type factor, const BlockMatrixBase<MatrixType> &matrix);
613
621 operator()(const size_type i, const size_type j) const;
622
632 el(const size_type i, const size_type j) const;
633
645 diag_element(const size_type i) const;
646
655 void
657
662 operator*=(const value_type factor);
663
668 operator/=(const value_type factor);
669
674 template <class BlockVectorType>
675 void
676 vmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
677
683 template <class BlockVectorType>
684 void
685 Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
686
699 template <class BlockVectorType>
701 matrix_norm_square(const BlockVectorType &v) const;
702
709
713 template <class BlockVectorType>
715 matrix_scalar_product(const BlockVectorType &u,
716 const BlockVectorType &v) const;
717
721 template <class BlockVectorType>
723 residual(BlockVectorType & dst,
724 const BlockVectorType &x,
725 const BlockVectorType &b) const;
726
733 void
734 print(std::ostream &out, const bool alternative_output = false) const;
735
741
747
752 begin(const size_type r);
753
758 end(const size_type r);
763 begin() const;
764
769 end() const;
770
775 begin(const size_type r) const;
776
781 end(const size_type r) const;
782
786 const BlockIndices &
788
792 const BlockIndices &
794
800 std::size_t
802
812 int,
813 int,
814 int,
815 int,
816 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
817 << ',' << arg4 << "] have differing row numbers.");
822 int,
823 int,
824 int,
825 int,
826 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
827 << ',' << arg4 << "] have differing column numbers.");
829protected:
842 void
844
850
855
874 void
876
887 template <class BlockVectorType>
888 void
889 vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const;
890
901 template <class BlockVectorType, class VectorType>
902 void
903 vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const;
904
915 template <class BlockVectorType, class VectorType>
916 void
917 vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const;
918
929 template <class VectorType>
930 void
931 vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const;
932
944 template <class BlockVectorType>
945 void
946 Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const;
947
958 template <class BlockVectorType, class VectorType>
959 void
960 Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const;
961
972 template <class BlockVectorType, class VectorType>
973 void
974 Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const;
975
986 template <class VectorType>
987 void
988 Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const;
989
990
991protected:
998 void
1000
1005 void
1007
1008
1009private:
1019 {
1024 std::vector<size_type> counter_within_block;
1025
1030 std::vector<std::vector<size_type>> column_indices;
1031
1036 std::vector<std::vector<value_type>> column_values;
1037
1043
1055 {
1056 return *this;
1057 }
1058 };
1059
1067
1068 // Make the iterator class a friend. We have to work around a compiler bug
1069 // here again.
1070 template <typename, bool>
1072
1073 template <typename>
1074 friend class MatrixIterator;
1075};
1076
1077
1080#ifndef DOXYGEN
1081/* ------------------------- Template functions ---------------------- */
1082
1083
1084namespace BlockMatrixIterators
1085{
1086 template <class BlockMatrixType>
1088 : row_block(0)
1089 , col_block(0)
1090 {}
1091
1092
1093 template <class BlockMatrixType>
1094 inline unsigned int
1095 AccessorBase<BlockMatrixType>::block_row() const
1096 {
1098
1099 return row_block;
1100 }
1101
1102
1103 template <class BlockMatrixType>
1104 inline unsigned int
1105 AccessorBase<BlockMatrixType>::block_column() const
1106 {
1108
1109 return col_block;
1110 }
1111
1112
1113 template <class BlockMatrixType>
1114 inline Accessor<BlockMatrixType, true>::Accessor(
1115 const BlockMatrixType *matrix,
1116 const size_type row,
1117 const size_type col)
1118 : matrix(matrix)
1119 , base_iterator(matrix->block(0, 0).begin())
1120 {
1121 (void)col;
1122 Assert(col == 0, ExcNotImplemented());
1123
1124 // check if this is a regular row or
1125 // the end of the matrix
1126 if (row < matrix->m())
1127 {
1128 const std::pair<unsigned int, size_type> indices =
1129 matrix->row_block_indices.global_to_local(row);
1130
1131 // find the first block that does
1132 // have an entry in this row
1133 for (unsigned int bc = 0; bc < matrix->n_block_cols(); ++bc)
1134 {
1135 base_iterator =
1136 matrix->block(indices.first, bc).begin(indices.second);
1137 if (base_iterator !=
1138 matrix->block(indices.first, bc).end(indices.second))
1139 {
1140 this->row_block = indices.first;
1141 this->col_block = bc;
1142 return;
1143 }
1144 }
1145
1146 // hm, there is no block that has
1147 // an entry in this column. we need
1148 // to take the next entry then,
1149 // which may be the first entry of
1150 // the next row, or recursively the
1151 // next row, or so on
1152 *this = Accessor(matrix, row + 1, 0);
1153 }
1154 else
1155 {
1156 // we were asked to create the end
1157 // iterator for this matrix
1158 this->row_block = numbers::invalid_unsigned_int;
1159 this->col_block = numbers::invalid_unsigned_int;
1160 }
1161 }
1162
1163
1164 // template <class BlockMatrixType>
1165 // inline
1166 // Accessor<BlockMatrixType, true>::Accessor (const
1167 // Accessor<BlockMatrixType, true>& other)
1168 // :
1169 // matrix(other.matrix),
1170 // base_iterator(other.base_iterator)
1171 // {
1172 // this->row_block = other.row_block;
1173 // this->col_block = other.col_block;
1174 // }
1175
1176
1177 template <class BlockMatrixType>
1178 inline Accessor<BlockMatrixType, true>::Accessor(
1179 const Accessor<BlockMatrixType, false> &other)
1180 : matrix(other.matrix)
1181 , base_iterator(other.base_iterator)
1182 {
1183 this->row_block = other.row_block;
1184 this->col_block = other.col_block;
1185 }
1186
1187
1188 template <class BlockMatrixType>
1189 inline typename Accessor<BlockMatrixType, true>::size_type
1190 Accessor<BlockMatrixType, true>::row() const
1191 {
1192 Assert(this->row_block != numbers::invalid_unsigned_int,
1194
1195 return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1196 base_iterator->row());
1197 }
1198
1199
1200 template <class BlockMatrixType>
1201 inline typename Accessor<BlockMatrixType, true>::size_type
1202 Accessor<BlockMatrixType, true>::column() const
1203 {
1204 Assert(this->col_block != numbers::invalid_unsigned_int,
1206
1207 return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
1208 base_iterator->column());
1209 }
1210
1211
1212 template <class BlockMatrixType>
1213 inline typename Accessor<BlockMatrixType, true>::value_type
1214 Accessor<BlockMatrixType, true>::value() const
1215 {
1216 Assert(this->row_block != numbers::invalid_unsigned_int,
1218 Assert(this->col_block != numbers::invalid_unsigned_int,
1220
1221 return base_iterator->value();
1222 }
1223
1224
1225
1226 template <class BlockMatrixType>
1227 inline void
1228 Accessor<BlockMatrixType, true>::advance()
1229 {
1230 Assert(this->row_block != numbers::invalid_unsigned_int,
1232 Assert(this->col_block != numbers::invalid_unsigned_int,
1234
1235 // Remember current row inside block
1236 size_type local_row = base_iterator->row();
1237
1238 // Advance one element inside the
1239 // current block
1240 ++base_iterator;
1241
1242 // while we hit the end of the row of a
1243 // block (which may happen multiple
1244 // times if rows inside a block are
1245 // empty), we have to jump to the next
1246 // block and take the
1247 while (base_iterator ==
1248 matrix->block(this->row_block, this->col_block).end(local_row))
1249 {
1250 // jump to next block in this block
1251 // row, if possible, otherwise go
1252 // to next row
1253 if (this->col_block < matrix->n_block_cols() - 1)
1254 {
1255 ++this->col_block;
1256 base_iterator =
1257 matrix->block(this->row_block, this->col_block).begin(local_row);
1258 }
1259 else
1260 {
1261 // jump back to next row in
1262 // first block column
1263 this->col_block = 0;
1264 ++local_row;
1265
1266 // see if this has brought us
1267 // past the number of rows in
1268 // this block. if so see
1269 // whether we've just fallen
1270 // off the end of the whole
1271 // matrix
1272 if (local_row ==
1273 matrix->block(this->row_block, this->col_block).m())
1274 {
1275 local_row = 0;
1276 ++this->row_block;
1277 if (this->row_block == matrix->n_block_rows())
1278 {
1279 this->row_block = numbers::invalid_unsigned_int;
1280 this->col_block = numbers::invalid_unsigned_int;
1281 return;
1282 }
1283 }
1284
1285 base_iterator =
1286 matrix->block(this->row_block, this->col_block).begin(local_row);
1287 }
1288 }
1289 }
1290
1291
1292 template <class BlockMatrixType>
1293 inline bool
1294 Accessor<BlockMatrixType, true>::operator==(const Accessor &a) const
1295 {
1296 if (matrix != a.matrix)
1297 return false;
1298
1299 if (this->row_block == a.row_block && this->col_block == a.col_block)
1300 // end iterators do not necessarily
1301 // have to have the same
1302 // base_iterator representation, but
1303 // valid iterators have to
1304 return (((this->row_block == numbers::invalid_unsigned_int) &&
1305 (this->col_block == numbers::invalid_unsigned_int)) ||
1306 (base_iterator == a.base_iterator));
1307
1308 return false;
1309 }
1310
1311 //----------------------------------------------------------------------//
1312
1313
1314 template <class BlockMatrixType>
1315 inline Accessor<BlockMatrixType, false>::Accessor(BlockMatrixType *matrix,
1316 const size_type row,
1317 const size_type col)
1318 : matrix(matrix)
1319 , base_iterator(matrix->block(0, 0).begin())
1320 {
1321 (void)col;
1322 Assert(col == 0, ExcNotImplemented());
1323 // check if this is a regular row or
1324 // the end of the matrix
1325 if (row < matrix->m())
1326 {
1327 const std::pair<unsigned int, size_type> indices =
1328 matrix->row_block_indices.global_to_local(row);
1329
1330 // find the first block that does
1331 // have an entry in this row
1332 for (size_type bc = 0; bc < matrix->n_block_cols(); ++bc)
1333 {
1334 base_iterator =
1335 matrix->block(indices.first, bc).begin(indices.second);
1336 if (base_iterator !=
1337 matrix->block(indices.first, bc).end(indices.second))
1338 {
1339 this->row_block = indices.first;
1340 this->col_block = bc;
1341 return;
1342 }
1343 }
1344
1345 // hm, there is no block that has
1346 // an entry in this column. we need
1347 // to take the next entry then,
1348 // which may be the first entry of
1349 // the next row, or recursively the
1350 // next row, or so on
1351 *this = Accessor(matrix, row + 1, 0);
1352 }
1353 else
1354 {
1355 // we were asked to create the end
1356 // iterator for this matrix
1357 this->row_block = numbers::invalid_size_type;
1358 this->col_block = numbers::invalid_size_type;
1359 }
1360 }
1361
1362
1363 template <class BlockMatrixType>
1364 inline typename Accessor<BlockMatrixType, false>::size_type
1365 Accessor<BlockMatrixType, false>::row() const
1366 {
1368
1369 return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1370 base_iterator->row());
1371 }
1372
1373
1374 template <class BlockMatrixType>
1375 inline typename Accessor<BlockMatrixType, false>::size_type
1376 Accessor<BlockMatrixType, false>::column() const
1377 {
1379
1380 return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
1381 base_iterator->column());
1382 }
1383
1384
1385 template <class BlockMatrixType>
1386 inline typename Accessor<BlockMatrixType, false>::value_type
1387 Accessor<BlockMatrixType, false>::value() const
1388 {
1391
1392 return base_iterator->value();
1393 }
1394
1395
1396
1397 template <class BlockMatrixType>
1398 inline void
1399 Accessor<BlockMatrixType, false>::set_value(
1400 typename Accessor<BlockMatrixType, false>::value_type newval) const
1401 {
1404
1405 base_iterator->value() = newval;
1406 }
1407
1408
1409
1410 template <class BlockMatrixType>
1411 inline void
1412 Accessor<BlockMatrixType, false>::advance()
1413 {
1416
1417 // Remember current row inside block
1418 size_type local_row = base_iterator->row();
1419
1420 // Advance one element inside the
1421 // current block
1422 ++base_iterator;
1423
1424 // while we hit the end of the row of a
1425 // block (which may happen multiple
1426 // times if rows inside a block are
1427 // empty), we have to jump to the next
1428 // block and take the
1429 while (base_iterator ==
1430 matrix->block(this->row_block, this->col_block).end(local_row))
1431 {
1432 // jump to next block in this block
1433 // row, if possible, otherwise go
1434 // to next row
1435 if (this->col_block < matrix->n_block_cols() - 1)
1436 {
1437 ++this->col_block;
1438 base_iterator =
1439 matrix->block(this->row_block, this->col_block).begin(local_row);
1440 }
1441 else
1442 {
1443 // jump back to next row in
1444 // first block column
1445 this->col_block = 0;
1446 ++local_row;
1447
1448 // see if this has brought us
1449 // past the number of rows in
1450 // this block. if so see
1451 // whether we've just fallen
1452 // off the end of the whole
1453 // matrix
1454 if (local_row ==
1455 matrix->block(this->row_block, this->col_block).m())
1456 {
1457 local_row = 0;
1458 ++this->row_block;
1459 if (this->row_block == matrix->n_block_rows())
1460 {
1461 this->row_block = numbers::invalid_size_type;
1462 this->col_block = numbers::invalid_size_type;
1463 return;
1464 }
1465 }
1466
1467 base_iterator =
1468 matrix->block(this->row_block, this->col_block).begin(local_row);
1469 }
1470 }
1471 }
1472
1473
1474
1475 template <class BlockMatrixType>
1476 inline bool
1477 Accessor<BlockMatrixType, false>::operator==(const Accessor &a) const
1478 {
1479 if (matrix != a.matrix)
1480 return false;
1481
1482 if (this->row_block == a.row_block && this->col_block == a.col_block)
1483 // end iterators do not necessarily
1484 // have to have the same
1485 // base_iterator representation, but
1486 // valid iterators have to
1487 return (((this->row_block == numbers::invalid_size_type) &&
1488 (this->col_block == numbers::invalid_size_type)) ||
1489 (base_iterator == a.base_iterator));
1490
1491 return false;
1492 }
1493} // namespace BlockMatrixIterators
1494
1495
1496//---------------------------------------------------------------------------
1497
1498template <typename MatrixType>
1500{
1501 try
1502 {
1503 clear();
1504 }
1505 catch (...)
1506 {}
1507}
1508
1509
1510template <class MatrixType>
1511template <class BlockMatrixType>
1513BlockMatrixBase<MatrixType>::copy_from(const BlockMatrixType &source)
1514{
1515 for (unsigned int r = 0; r < n_block_rows(); ++r)
1516 for (unsigned int c = 0; c < n_block_cols(); ++c)
1517 block(r, c).copy_from(source.block(r, c));
1518
1519 return *this;
1520}
1521
1522
1523template <class MatrixType>
1524std::size_t
1526{
1527 std::size_t mem =
1528 MemoryConsumption::memory_consumption(row_block_indices) +
1529 MemoryConsumption::memory_consumption(column_block_indices) +
1531 MemoryConsumption::memory_consumption(temporary_data.counter_within_block) +
1532 MemoryConsumption::memory_consumption(temporary_data.column_indices) +
1533 MemoryConsumption::memory_consumption(temporary_data.column_values) +
1534 sizeof(temporary_data.mutex);
1535
1536 for (unsigned int r = 0; r < n_block_rows(); ++r)
1537 for (unsigned int c = 0; c < n_block_cols(); ++c)
1538 {
1539 MatrixType *p = this->sub_objects[r][c];
1541 }
1542
1543 return mem;
1544}
1545
1546
1547
1548template <class MatrixType>
1549inline void
1551{
1552 for (unsigned int r = 0; r < n_block_rows(); ++r)
1553 for (unsigned int c = 0; c < n_block_cols(); ++c)
1554 {
1555 MatrixType *p = this->sub_objects[r][c];
1556 this->sub_objects[r][c] = nullptr;
1557 delete p;
1558 }
1559 sub_objects.reinit(0, 0);
1560
1561 // reset block indices to empty
1562 row_block_indices = column_block_indices = BlockIndices();
1563}
1564
1565
1566
1567template <class MatrixType>
1569BlockMatrixBase<MatrixType>::block(const unsigned int row,
1570 const unsigned int column)
1571{
1572 AssertIndexRange(row, n_block_rows());
1573 AssertIndexRange(column, n_block_cols());
1574
1575 return *sub_objects[row][column];
1576}
1577
1578
1579
1580template <class MatrixType>
1581inline const typename BlockMatrixBase<MatrixType>::BlockType &
1582BlockMatrixBase<MatrixType>::block(const unsigned int row,
1583 const unsigned int column) const
1584{
1585 AssertIndexRange(row, n_block_rows());
1586 AssertIndexRange(column, n_block_cols());
1587
1588 return *sub_objects[row][column];
1589}
1590
1591
1592template <class MatrixType>
1595{
1596 return row_block_indices.total_size();
1597}
1598
1599
1600
1601template <class MatrixType>
1604{
1605 return column_block_indices.total_size();
1606}
1607
1608
1609
1610template <class MatrixType>
1611inline unsigned int
1613{
1614 return column_block_indices.size();
1615}
1616
1617
1618
1619template <class MatrixType>
1620inline unsigned int
1622{
1623 return row_block_indices.size();
1624}
1625
1626
1627
1628// Write the single set manually,
1629// since the other function has a lot
1630// of overhead in that case.
1631template <class MatrixType>
1632inline void
1633BlockMatrixBase<MatrixType>::set(const size_type i,
1634 const size_type j,
1635 const value_type value)
1636{
1637 prepare_set_operation();
1638
1639 AssertIsFinite(value);
1640
1641 const std::pair<unsigned int, size_type>
1642 row_index = row_block_indices.global_to_local(i),
1643 col_index = column_block_indices.global_to_local(j);
1644 block(row_index.first, col_index.first)
1645 .set(row_index.second, col_index.second, value);
1646}
1647
1648
1649
1650template <class MatrixType>
1651template <typename number>
1652inline void
1653BlockMatrixBase<MatrixType>::set(const std::vector<size_type> &row_indices,
1654 const std::vector<size_type> &col_indices,
1655 const FullMatrix<number> & values,
1656 const bool elide_zero_values)
1657{
1658 Assert(row_indices.size() == values.m(),
1659 ExcDimensionMismatch(row_indices.size(), values.m()));
1660 Assert(col_indices.size() == values.n(),
1661 ExcDimensionMismatch(col_indices.size(), values.n()));
1662
1663 for (size_type i = 0; i < row_indices.size(); ++i)
1664 set(row_indices[i],
1665 col_indices.size(),
1666 col_indices.data(),
1667 &values(i, 0),
1668 elide_zero_values);
1669}
1670
1671
1672
1673template <class MatrixType>
1674template <typename number>
1675inline void
1676BlockMatrixBase<MatrixType>::set(const std::vector<size_type> &indices,
1677 const FullMatrix<number> & values,
1678 const bool elide_zero_values)
1679{
1680 Assert(indices.size() == values.m(),
1681 ExcDimensionMismatch(indices.size(), values.m()));
1682 Assert(values.n() == values.m(), ExcNotQuadratic());
1683
1684 for (size_type i = 0; i < indices.size(); ++i)
1685 set(indices[i],
1686 indices.size(),
1687 indices.data(),
1688 &values(i, 0),
1689 elide_zero_values);
1690}
1691
1692
1693
1694template <class MatrixType>
1695template <typename number>
1696inline void
1697BlockMatrixBase<MatrixType>::set(const size_type row,
1698 const std::vector<size_type> &col_indices,
1699 const std::vector<number> & values,
1700 const bool elide_zero_values)
1701{
1702 Assert(col_indices.size() == values.size(),
1703 ExcDimensionMismatch(col_indices.size(), values.size()));
1704
1705 set(row,
1706 col_indices.size(),
1707 col_indices.data(),
1708 values.data(),
1709 elide_zero_values);
1710}
1711
1712
1713
1714// This is a very messy function, since
1715// we need to calculate to each position
1716// the location in the global array.
1717template <class MatrixType>
1718template <typename number>
1719inline void
1720BlockMatrixBase<MatrixType>::set(const size_type row,
1721 const size_type n_cols,
1722 const size_type *col_indices,
1723 const number * values,
1724 const bool elide_zero_values)
1725{
1726 prepare_set_operation();
1727
1728 // lock access to the temporary data structure to
1729 // allow multiple threads to call this function concurrently
1730 std::lock_guard<std::mutex> lock(temporary_data.mutex);
1731
1732 // Resize scratch arrays
1733 if (temporary_data.column_indices.size() < this->n_block_cols())
1734 {
1735 temporary_data.column_indices.resize(this->n_block_cols());
1736 temporary_data.column_values.resize(this->n_block_cols());
1737 temporary_data.counter_within_block.resize(this->n_block_cols());
1738 }
1739
1740 // Resize sub-arrays to n_cols. This
1741 // is a bit wasteful, but we resize
1742 // only a few times (then the maximum
1743 // row length won't increase that
1744 // much any more). At least we know
1745 // that all arrays are going to be of
1746 // the same size, so we can check
1747 // whether the size of one is large
1748 // enough before actually going
1749 // through all of them.
1750 if (temporary_data.column_indices[0].size() < n_cols)
1751 {
1752 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1753 {
1754 temporary_data.column_indices[i].resize(n_cols);
1755 temporary_data.column_values[i].resize(n_cols);
1756 }
1757 }
1758
1759 // Reset the number of added elements
1760 // in each block to zero.
1761 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1762 temporary_data.counter_within_block[i] = 0;
1763
1764 // Go through the column indices to
1765 // find out which portions of the
1766 // values should be set in which
1767 // block of the matrix. We need to
1768 // touch all the data, since we can't
1769 // be sure that the data of one block
1770 // is stored contiguously (in fact,
1771 // indices will be intermixed when it
1772 // comes from an element matrix).
1773 for (size_type j = 0; j < n_cols; ++j)
1774 {
1775 number value = values[j];
1776
1777 if (value == number() && elide_zero_values == true)
1778 continue;
1779
1780 const std::pair<unsigned int, size_type> col_index =
1781 this->column_block_indices.global_to_local(col_indices[j]);
1782
1783 const size_type local_index =
1784 temporary_data.counter_within_block[col_index.first]++;
1785
1786 temporary_data.column_indices[col_index.first][local_index] =
1787 col_index.second;
1788 temporary_data.column_values[col_index.first][local_index] = value;
1789 }
1790
1791# ifdef DEBUG
1792 // If in debug mode, do a check whether
1793 // the right length has been obtained.
1794 size_type length = 0;
1795 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1796 length += temporary_data.counter_within_block[i];
1797 Assert(length <= n_cols, ExcInternalError());
1798# endif
1799
1800 // Now we found out about where the
1801 // individual columns should start and
1802 // where we should start reading out
1803 // data. Now let's write the data into
1804 // the individual blocks!
1805 const std::pair<unsigned int, size_type> row_index =
1806 this->row_block_indices.global_to_local(row);
1807 for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
1808 {
1809 if (temporary_data.counter_within_block[block_col] == 0)
1810 continue;
1811
1812 block(row_index.first, block_col)
1813 .set(row_index.second,
1814 temporary_data.counter_within_block[block_col],
1815 temporary_data.column_indices[block_col].data(),
1816 temporary_data.column_values[block_col].data(),
1817 false);
1818 }
1819}
1820
1821
1822
1823template <class MatrixType>
1824inline void
1825BlockMatrixBase<MatrixType>::add(const size_type i,
1826 const size_type j,
1827 const value_type value)
1828{
1829 AssertIsFinite(value);
1830
1831 prepare_add_operation();
1832
1833 // save some cycles for zero additions, but
1834 // only if it is safe for the matrix we are
1835 // working with
1836 using MatrixTraits = typename MatrixType::Traits;
1837 if ((MatrixTraits::zero_addition_can_be_elided == true) &&
1838 (value == value_type()))
1839 return;
1840
1841 const std::pair<unsigned int, size_type>
1842 row_index = row_block_indices.global_to_local(i),
1843 col_index = column_block_indices.global_to_local(j);
1844 block(row_index.first, col_index.first)
1845 .add(row_index.second, col_index.second, value);
1846}
1847
1848
1849
1850template <class MatrixType>
1851template <typename number>
1852inline void
1853BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &row_indices,
1854 const std::vector<size_type> &col_indices,
1855 const FullMatrix<number> & values,
1856 const bool elide_zero_values)
1857{
1858 Assert(row_indices.size() == values.m(),
1859 ExcDimensionMismatch(row_indices.size(), values.m()));
1860 Assert(col_indices.size() == values.n(),
1861 ExcDimensionMismatch(col_indices.size(), values.n()));
1862
1863 for (size_type i = 0; i < row_indices.size(); ++i)
1864 add(row_indices[i],
1865 col_indices.size(),
1866 col_indices.data(),
1867 &values(i, 0),
1868 elide_zero_values);
1869}
1870
1871
1872
1873template <class MatrixType>
1874template <typename number>
1875inline void
1876BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &indices,
1877 const FullMatrix<number> & values,
1878 const bool elide_zero_values)
1879{
1880 Assert(indices.size() == values.m(),
1881 ExcDimensionMismatch(indices.size(), values.m()));
1882 Assert(values.n() == values.m(), ExcNotQuadratic());
1883
1884 for (size_type i = 0; i < indices.size(); ++i)
1885 add(indices[i],
1886 indices.size(),
1887 indices.data(),
1888 &values(i, 0),
1889 elide_zero_values);
1890}
1891
1892
1893
1894template <class MatrixType>
1895template <typename number>
1896inline void
1897BlockMatrixBase<MatrixType>::add(const size_type row,
1898 const std::vector<size_type> &col_indices,
1899 const std::vector<number> & values,
1900 const bool elide_zero_values)
1901{
1902 Assert(col_indices.size() == values.size(),
1903 ExcDimensionMismatch(col_indices.size(), values.size()));
1904
1905 add(row,
1906 col_indices.size(),
1907 col_indices.data(),
1908 values.data(),
1909 elide_zero_values);
1910}
1911
1912
1913
1914// This is a very messy function, since
1915// we need to calculate to each position
1916// the location in the global array.
1917template <class MatrixType>
1918template <typename number>
1919inline void
1920BlockMatrixBase<MatrixType>::add(const size_type row,
1921 const size_type n_cols,
1922 const size_type *col_indices,
1923 const number * values,
1924 const bool elide_zero_values,
1925 const bool col_indices_are_sorted)
1926{
1927 prepare_add_operation();
1928
1929 // TODO: Look over this to find out
1930 // whether we can do that more
1931 // efficiently.
1932 if (col_indices_are_sorted == true)
1933 {
1934# ifdef DEBUG
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 Assert(false,
1941 ExcMessage("Flag col_indices_are_sorted is set, but "
1942 "indices appear to not be sorted.")) else before =
1943 col_indices[i];
1944# endif
1945 const std::pair<unsigned int, size_type> row_index =
1946 this->row_block_indices.global_to_local(row);
1947
1948 if (this->n_block_cols() > 1)
1949 {
1950 const size_type *first_block =
1951 Utilities::lower_bound(col_indices,
1952 col_indices + n_cols,
1953 this->column_block_indices.block_start(1));
1954
1955 const size_type n_zero_block_indices = first_block - col_indices;
1956 block(row_index.first, 0)
1957 .add(row_index.second,
1958 n_zero_block_indices,
1959 col_indices,
1960 values,
1961 elide_zero_values,
1962 col_indices_are_sorted);
1963
1964 if (n_zero_block_indices < n_cols)
1965 this->add(row,
1966 n_cols - n_zero_block_indices,
1967 first_block,
1968 values + n_zero_block_indices,
1969 elide_zero_values,
1970 false);
1971 }
1972 else
1973 {
1974 block(row_index.first, 0)
1975 .add(row_index.second,
1976 n_cols,
1977 col_indices,
1978 values,
1979 elide_zero_values,
1980 col_indices_are_sorted);
1981 }
1982
1983 return;
1984 }
1985
1986 // Lock scratch arrays, then resize them
1987 std::lock_guard<std::mutex> lock(temporary_data.mutex);
1988
1989 if (temporary_data.column_indices.size() < this->n_block_cols())
1990 {
1991 temporary_data.column_indices.resize(this->n_block_cols());
1992 temporary_data.column_values.resize(this->n_block_cols());
1993 temporary_data.counter_within_block.resize(this->n_block_cols());
1994 }
1995
1996 // Resize sub-arrays to n_cols. This
1997 // is a bit wasteful, but we resize
1998 // only a few times (then the maximum
1999 // row length won't increase that
2000 // much any more). At least we know
2001 // that all arrays are going to be of
2002 // the same size, so we can check
2003 // whether the size of one is large
2004 // enough before actually going
2005 // through all of them.
2006 if (temporary_data.column_indices[0].size() < n_cols)
2007 {
2008 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2009 {
2010 temporary_data.column_indices[i].resize(n_cols);
2011 temporary_data.column_values[i].resize(n_cols);
2012 }
2013 }
2014
2015 // Reset the number of added elements
2016 // in each block to zero.
2017 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2018 temporary_data.counter_within_block[i] = 0;
2019
2020 // Go through the column indices to
2021 // find out which portions of the
2022 // values should be written into
2023 // which block of the matrix. We need
2024 // to touch all the data, since we
2025 // can't be sure that the data of one
2026 // block is stored contiguously (in
2027 // fact, data will be intermixed when
2028 // it comes from an element matrix).
2029 for (size_type j = 0; j < n_cols; ++j)
2030 {
2031 number value = values[j];
2032
2033 if (value == number() && elide_zero_values == true)
2034 continue;
2035
2036 const std::pair<unsigned int, size_type> col_index =
2037 this->column_block_indices.global_to_local(col_indices[j]);
2038
2039 const size_type local_index =
2040 temporary_data.counter_within_block[col_index.first]++;
2041
2042 temporary_data.column_indices[col_index.first][local_index] =
2043 col_index.second;
2044 temporary_data.column_values[col_index.first][local_index] = value;
2045 }
2046
2047# ifdef DEBUG
2048 // If in debug mode, do a check whether
2049 // the right length has been obtained.
2050 size_type length = 0;
2051 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2052 length += temporary_data.counter_within_block[i];
2053 Assert(length <= n_cols, ExcInternalError());
2054# endif
2055
2056 // Now we found out about where the
2057 // individual columns should start and
2058 // where we should start reading out
2059 // data. Now let's write the data into
2060 // the individual blocks!
2061 const std::pair<unsigned int, size_type> row_index =
2062 this->row_block_indices.global_to_local(row);
2063 for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
2064 {
2065 if (temporary_data.counter_within_block[block_col] == 0)
2066 continue;
2067
2068 block(row_index.first, block_col)
2069 .add(row_index.second,
2070 temporary_data.counter_within_block[block_col],
2071 temporary_data.column_indices[block_col].data(),
2072 temporary_data.column_values[block_col].data(),
2073 false,
2074 col_indices_are_sorted);
2075 }
2076}
2077
2078
2079
2080template <class MatrixType>
2081inline void
2082BlockMatrixBase<MatrixType>::add(const value_type factor,
2083 const BlockMatrixBase<MatrixType> &matrix)
2084{
2085 AssertIsFinite(factor);
2086
2087 prepare_add_operation();
2088
2089 // save some cycles for zero additions, but
2090 // only if it is safe for the matrix we are
2091 // working with
2092 using MatrixTraits = typename MatrixType::Traits;
2093 if ((MatrixTraits::zero_addition_can_be_elided == true) && (factor == 0))
2094 return;
2095
2096 for (unsigned int row = 0; row < n_block_rows(); ++row)
2097 for (unsigned int col = 0; col < n_block_cols(); ++col)
2098 // This function should throw if the sparsity
2099 // patterns of the two blocks differ
2100 block(row, col).add(factor, matrix.block(row, col));
2101}
2102
2103
2104
2105template <class MatrixType>
2108 const size_type j) const
2109{
2110 const std::pair<unsigned int, size_type>
2111 row_index = row_block_indices.global_to_local(i),
2112 col_index = column_block_indices.global_to_local(j);
2113 return block(row_index.first, col_index.first)(row_index.second,
2114 col_index.second);
2115}
2116
2117
2118
2119template <class MatrixType>
2121BlockMatrixBase<MatrixType>::el(const size_type i, const size_type j) const
2122{
2123 const std::pair<unsigned int, size_type>
2124 row_index = row_block_indices.global_to_local(i),
2125 col_index = column_block_indices.global_to_local(j);
2126 return block(row_index.first, col_index.first)
2127 .el(row_index.second, col_index.second);
2128}
2129
2130
2131
2132template <class MatrixType>
2134BlockMatrixBase<MatrixType>::diag_element(const size_type i) const
2135{
2136 Assert(n_block_rows() == n_block_cols(), ExcNotQuadratic());
2137
2138 const std::pair<unsigned int, size_type> index =
2139 row_block_indices.global_to_local(i);
2140 return block(index.first, index.first).diag_element(index.second);
2141}
2142
2143
2144
2145template <class MatrixType>
2146inline void
2148{
2149 for (unsigned int r = 0; r < n_block_rows(); ++r)
2150 for (unsigned int c = 0; c < n_block_cols(); ++c)
2151 block(r, c).compress(operation);
2152}
2153
2154
2155
2156template <class MatrixType>
2158BlockMatrixBase<MatrixType>::operator*=(const value_type factor)
2159{
2160 Assert(n_block_cols() != 0, ExcNotInitialized());
2161 Assert(n_block_rows() != 0, ExcNotInitialized());
2162
2163 for (unsigned int r = 0; r < n_block_rows(); ++r)
2164 for (unsigned int c = 0; c < n_block_cols(); ++c)
2165 block(r, c) *= factor;
2166
2167 return *this;
2168}
2169
2170
2171
2172template <class MatrixType>
2174BlockMatrixBase<MatrixType>::operator/=(const value_type factor)
2175{
2176 Assert(n_block_cols() != 0, ExcNotInitialized());
2177 Assert(n_block_rows() != 0, ExcNotInitialized());
2178 Assert(factor != 0, ExcDivideByZero());
2179
2180 const value_type factor_inv = 1. / factor;
2181
2182 for (unsigned int r = 0; r < n_block_rows(); ++r)
2183 for (unsigned int c = 0; c < n_block_cols(); ++c)
2184 block(r, c) *= factor_inv;
2185
2186 return *this;
2187}
2188
2189
2190
2191template <class MatrixType>
2192const BlockIndices &
2194{
2195 return this->row_block_indices;
2196}
2197
2198
2199
2200template <class MatrixType>
2201const BlockIndices &
2203{
2204 return this->column_block_indices;
2205}
2206
2207
2208
2209template <class MatrixType>
2210template <class BlockVectorType>
2211void
2213 const BlockVectorType &src) const
2214{
2215 Assert(dst.n_blocks() == n_block_rows(),
2216 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2217 Assert(src.n_blocks() == n_block_cols(),
2218 ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2219
2220 for (size_type row = 0; row < n_block_rows(); ++row)
2221 {
2222 block(row, 0).vmult(dst.block(row), src.block(0));
2223 for (size_type col = 1; col < n_block_cols(); ++col)
2224 block(row, col).vmult_add(dst.block(row), src.block(col));
2225 };
2226}
2227
2228
2229
2230template <class MatrixType>
2231template <class BlockVectorType, class VectorType>
2232void
2234 VectorType & dst,
2235 const BlockVectorType &src) const
2236{
2237 Assert(n_block_rows() == 1, ExcDimensionMismatch(1, n_block_rows()));
2238 Assert(src.n_blocks() == n_block_cols(),
2239 ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2240
2241 block(0, 0).vmult(dst, src.block(0));
2242 for (size_type col = 1; col < n_block_cols(); ++col)
2243 block(0, col).vmult_add(dst, src.block(col));
2244}
2245
2246
2247
2248template <class MatrixType>
2249template <class BlockVectorType, class VectorType>
2250void
2252 const VectorType &src) const
2253{
2254 Assert(dst.n_blocks() == n_block_rows(),
2255 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2256 Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2257
2258 for (size_type row = 0; row < n_block_rows(); ++row)
2259 block(row, 0).vmult(dst.block(row), src);
2260}
2261
2262
2263
2264template <class MatrixType>
2265template <class VectorType>
2266void
2268 VectorType & dst,
2269 const VectorType &src) const
2270{
2271 Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));
2272 Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2273
2274 block(0, 0).vmult(dst, src);
2275}
2276
2277
2278
2279template <class MatrixType>
2280template <class BlockVectorType>
2281void
2282BlockMatrixBase<MatrixType>::vmult_add(BlockVectorType & dst,
2283 const BlockVectorType &src) const
2284{
2285 Assert(dst.n_blocks() == n_block_rows(),
2286 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2287 Assert(src.n_blocks() == n_block_cols(),
2288 ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2289
2290 for (unsigned int row = 0; row < n_block_rows(); ++row)
2291 for (unsigned int col = 0; col < n_block_cols(); ++col)
2292 block(row, col).vmult_add(dst.block(row), src.block(col));
2293}
2294
2295
2296
2297template <class MatrixType>
2298template <class BlockVectorType>
2299void
2301 BlockVectorType & dst,
2302 const BlockVectorType &src) const
2303{
2304 Assert(dst.n_blocks() == n_block_cols(),
2305 ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2306 Assert(src.n_blocks() == n_block_rows(),
2307 ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2308
2309 dst = 0.;
2310
2311 for (unsigned int row = 0; row < n_block_rows(); ++row)
2312 {
2313 for (unsigned int col = 0; col < n_block_cols(); ++col)
2314 block(row, col).Tvmult_add(dst.block(col), src.block(row));
2315 };
2316}
2317
2318
2319
2320template <class MatrixType>
2321template <class BlockVectorType, class VectorType>
2322void
2324 const VectorType &src) const
2325{
2326 Assert(dst.n_blocks() == n_block_cols(),
2327 ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2328 Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));
2329
2330 dst = 0.;
2331
2332 for (unsigned int col = 0; col < n_block_cols(); ++col)
2333 block(0, col).Tvmult_add(dst.block(col), src);
2334}
2335
2336
2337
2338template <class MatrixType>
2339template <class BlockVectorType, class VectorType>
2340void
2342 VectorType & dst,
2343 const BlockVectorType &src) const
2344{
2345 Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2346 Assert(src.n_blocks() == n_block_rows(),
2347 ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2348
2349 block(0, 0).Tvmult(dst, src.block(0));
2350
2351 for (size_type row = 1; row < n_block_rows(); ++row)
2352 block(row, 0).Tvmult_add(dst, src.block(row));
2353}
2354
2355
2356
2357template <class MatrixType>
2358template <class VectorType>
2359void
2361 VectorType & dst,
2362 const VectorType &src) const
2363{
2364 Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2365 Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));
2366
2367 block(0, 0).Tvmult(dst, src);
2368}
2369
2370
2371
2372template <class MatrixType>
2373template <class BlockVectorType>
2374void
2375BlockMatrixBase<MatrixType>::Tvmult_add(BlockVectorType & dst,
2376 const BlockVectorType &src) const
2377{
2378 Assert(dst.n_blocks() == n_block_cols(),
2379 ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2380 Assert(src.n_blocks() == n_block_rows(),
2381 ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2382
2383 for (unsigned int row = 0; row < n_block_rows(); ++row)
2384 for (unsigned int col = 0; col < n_block_cols(); ++col)
2385 block(row, col).Tvmult_add(dst.block(col), src.block(row));
2386}
2387
2388
2389
2390template <class MatrixType>
2391template <class BlockVectorType>
2393BlockMatrixBase<MatrixType>::matrix_norm_square(const BlockVectorType &v) const
2394{
2395 Assert(n_block_rows() == n_block_cols(), ExcNotQuadratic());
2396 Assert(v.n_blocks() == n_block_rows(),
2397 ExcDimensionMismatch(v.n_blocks(), n_block_rows()));
2398
2399 value_type norm_sqr = 0;
2400 for (unsigned int row = 0; row < n_block_rows(); ++row)
2401 for (unsigned int col = 0; col < n_block_cols(); ++col)
2402 if (row == col)
2403 norm_sqr += block(row, col).matrix_norm_square(v.block(row));
2404 else
2405 norm_sqr +=
2406 block(row, col).matrix_scalar_product(v.block(row), v.block(col));
2407 return norm_sqr;
2408}
2409
2410
2411
2412template <class MatrixType>
2415{
2416 value_type norm_sqr = 0;
2417
2418 // For each block, get the Frobenius norm, and add the square to the
2419 // accumulator for the full matrix
2420 for (unsigned int row = 0; row < n_block_rows(); ++row)
2421 {
2422 for (unsigned int col = 0; col < n_block_cols(); ++col)
2423 {
2424 const value_type block_norm = block(row, col).frobenius_norm();
2425 norm_sqr += block_norm * block_norm;
2426 }
2427 }
2428
2429 return std::sqrt(norm_sqr);
2430}
2431
2432
2433
2434template <class MatrixType>
2435template <class BlockVectorType>
2438 const BlockVectorType &u,
2439 const BlockVectorType &v) const
2440{
2441 Assert(u.n_blocks() == n_block_rows(),
2442 ExcDimensionMismatch(u.n_blocks(), n_block_rows()));
2443 Assert(v.n_blocks() == n_block_cols(),
2444 ExcDimensionMismatch(v.n_blocks(), n_block_cols()));
2445
2446 value_type result = 0;
2447 for (unsigned int row = 0; row < n_block_rows(); ++row)
2448 for (unsigned int col = 0; col < n_block_cols(); ++col)
2449 result +=
2450 block(row, col).matrix_scalar_product(u.block(row), v.block(col));
2451 return result;
2452}
2453
2454
2455
2456template <class MatrixType>
2457template <class BlockVectorType>
2459BlockMatrixBase<MatrixType>::residual(BlockVectorType & dst,
2460 const BlockVectorType &x,
2461 const BlockVectorType &b) const
2462{
2463 Assert(dst.n_blocks() == n_block_rows(),
2464 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2465 Assert(b.n_blocks() == n_block_rows(),
2466 ExcDimensionMismatch(b.n_blocks(), n_block_rows()));
2467 Assert(x.n_blocks() == n_block_cols(),
2468 ExcDimensionMismatch(x.n_blocks(), n_block_cols()));
2469 // in block notation, the residual is
2470 // r_i = b_i - \sum_j A_ij x_j.
2471 // this can be written as
2472 // r_i = b_i - A_i0 x_0 - \sum_{j>0} A_ij x_j.
2473 //
2474 // for the first two terms, we can
2475 // call the residual function of
2476 // A_i0. for the other terms, we
2477 // use vmult_add. however, we want
2478 // to subtract, so in order to
2479 // avoid a temporary vector, we
2480 // perform a sign change of the
2481 // first two term before, and after
2482 // adding up
2483 for (unsigned int row = 0; row < n_block_rows(); ++row)
2484 {
2485 block(row, 0).residual(dst.block(row), x.block(0), b.block(row));
2486
2487 for (size_type i = 0; i < dst.block(row).size(); ++i)
2488 dst.block(row)(i) = -dst.block(row)(i);
2489
2490 for (unsigned int col = 1; col < n_block_cols(); ++col)
2491 block(row, col).vmult_add(dst.block(row), x.block(col));
2492
2493 for (size_type i = 0; i < dst.block(row).size(); ++i)
2494 dst.block(row)(i) = -dst.block(row)(i);
2495 };
2496
2497 value_type res = 0;
2498 for (size_type row = 0; row < n_block_rows(); ++row)
2499 res += dst.block(row).norm_sqr();
2500 return std::sqrt(res);
2501}
2502
2503
2504
2505template <class MatrixType>
2506inline void
2507BlockMatrixBase<MatrixType>::print(std::ostream &out,
2508 const bool alternative_output) const
2509{
2510 for (unsigned int row = 0; row < n_block_rows(); ++row)
2511 for (unsigned int col = 0; col < n_block_cols(); ++col)
2512 {
2513 if (!alternative_output)
2514 out << "Block (" << row << ", " << col << ')' << std::endl;
2515
2516 block(row, col).print(out, alternative_output);
2517 }
2518}
2519
2520
2521
2522template <class MatrixType>
2525{
2526 return const_iterator(this, 0);
2527}
2528
2529
2530
2531template <class MatrixType>
2534{
2535 return const_iterator(this, m());
2536}
2537
2538
2539
2540template <class MatrixType>
2542BlockMatrixBase<MatrixType>::begin(const size_type r) const
2543{
2544 AssertIndexRange(r, m());
2545 return const_iterator(this, r);
2546}
2547
2548
2549
2550template <class MatrixType>
2552BlockMatrixBase<MatrixType>::end(const size_type r) const
2553{
2554 AssertIndexRange(r, m());
2555 return const_iterator(this, r + 1);
2556}
2557
2558
2559
2560template <class MatrixType>
2563{
2564 return iterator(this, 0);
2565}
2566
2567
2568
2569template <class MatrixType>
2572{
2573 return iterator(this, m());
2574}
2575
2576
2577
2578template <class MatrixType>
2580BlockMatrixBase<MatrixType>::begin(const size_type r)
2581{
2582 AssertIndexRange(r, m());
2583 return iterator(this, r);
2584}
2585
2586
2587
2588template <class MatrixType>
2590BlockMatrixBase<MatrixType>::end(const size_type r)
2591{
2592 AssertIndexRange(r, m());
2593 return iterator(this, r + 1);
2594}
2595
2596
2597
2598template <class MatrixType>
2599void
2601{
2602 std::vector<size_type> row_sizes(this->n_block_rows());
2603 std::vector<size_type> col_sizes(this->n_block_cols());
2604
2605 // first find out the row sizes
2606 // from the first block column
2607 for (unsigned int r = 0; r < this->n_block_rows(); ++r)
2608 row_sizes[r] = sub_objects[r][0]->m();
2609 // then check that the following
2610 // block columns have the same
2611 // sizes
2612 for (unsigned int c = 1; c < this->n_block_cols(); ++c)
2613 for (unsigned int r = 0; r < this->n_block_rows(); ++r)
2614 Assert(row_sizes[r] == sub_objects[r][c]->m(),
2615 ExcIncompatibleRowNumbers(r, 0, r, c));
2616
2617 // finally initialize the row
2618 // indices with this array
2619 this->row_block_indices.reinit(row_sizes);
2620
2621
2622 // then do the same with the columns
2623 for (unsigned int c = 0; c < this->n_block_cols(); ++c)
2624 col_sizes[c] = sub_objects[0][c]->n();
2625 for (unsigned int r = 1; r < this->n_block_rows(); ++r)
2626 for (unsigned int c = 0; c < this->n_block_cols(); ++c)
2627 Assert(col_sizes[c] == sub_objects[r][c]->n(),
2628 ExcIncompatibleRowNumbers(0, c, r, c));
2629
2630 // finally initialize the row
2631 // indices with this array
2632 this->column_block_indices.reinit(col_sizes);
2633}
2634
2635
2636
2637template <class MatrixType>
2638void
2640{
2641 for (unsigned int row = 0; row < n_block_rows(); ++row)
2642 for (unsigned int col = 0; col < n_block_cols(); ++col)
2643 block(row, col).prepare_add();
2644}
2645
2646
2647
2648template <class MatrixType>
2649void
2651{
2652 for (unsigned int row = 0; row < n_block_rows(); ++row)
2653 for (unsigned int col = 0; col < n_block_cols(); ++col)
2654 block(row, col).prepare_set();
2655}
2656
2657#endif // DOXYGEN
2658
2659
2661
2662#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
Table< 2, SmartPointer< BlockType, BlockMatrixBase< MatrixType > > > sub_objects
void Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const
typename BlockType::value_type value_type
void collect_sizes()
void prepare_add_operation()
void vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
value_type residual(BlockVectorType &dst, const BlockVectorType &x, const BlockVectorType &b) const
BlockMatrixBase()=default
void set(const size_type row, const size_type n_cols, const size_type *col_indices, const number *values, const bool elide_zero_values=false)
size_type n() const
void prepare_set_operation()
unsigned int n_block_cols() const
const value_type * const_pointer
iterator begin()
TemporaryData temporary_data
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
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
__global__ void set(Number *val, const Number s, const size_type N)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:579
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)
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
std::enable_if_t< std::is_fundamental< T >::value, 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:1016
static const unsigned int invalid_unsigned_int
Definition types.h:213
const types::global_dof_index invalid_size_type
Definition types.h:222
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition types.h:82
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