deal.II version GIT relicensing-2289-g1e5549a87a 2024-12-21 21:30: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\}}\)
Loading...
Searching...
No Matches
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
21#include <deal.II/base/memory_consumption.h>
22#include <deal.II/base/mutex.h>
23#include <deal.II/base/observer_pointer.h>
24#include <deal.II/base/table.h>
25#include <deal.II/base/utilities.h>
26
27#include <deal.II/lac/block_indices.h>
28#include <deal.II/lac/exceptions.h>
29#include <deal.II/lac/full_matrix.h>
30#include <deal.II/lac/matrix_iterator.h>
31#include <deal.II/lac/vector.h>
32#include <deal.II/lac/vector_operation.h>
33
34#include <cmath>
35#include <mutex>
36
37DEAL_II_NAMESPACE_OPEN
38
39
40// Forward declaration
41#ifndef DOXYGEN
42template <typename>
43class MatrixIterator;
44#endif
45
46
55namespace BlockMatrixIterators
56{
61 template <typename BlockMatrixType>
62 class AccessorBase
63 {
64 public:
68 using size_type = types::global_dof_index;
69
73 using value_type = typename BlockMatrixType::value_type;
74
78 AccessorBase();
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:
123 using size_type = types::global_dof_index;
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
148 size_type
149 row() const;
150
154 size_type
155 column() const;
156
160 value_type
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
184 advance();
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:
211 using size_type = types::global_dof_index;
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
238 Accessor(const Accessor<BlockMatrixType, false> &);
239
243 size_type
244 row() const;
245
249 size_type
250 column() const;
251
255 value_type
256 value() const;
257
258 protected:
262 const BlockMatrixType *matrix;
263
267 typename BlockMatrixType::BlockType::const_iterator base_iterator;
268
272 void
273 advance();
274
278 bool
279 operator==(const Accessor &a) const;
280
281 // Let the iterator class be a friend.
282 template <typename>
283 friend class ::MatrixIterator;
284 };
285} // namespace BlockMatrixIterators
286
287
288
348template <typename MatrixType>
349class BlockMatrixBase : public EnableObserverPointer
350{
351public:
355 using BlockType = MatrixType;
356
361 using value_type = typename BlockType::value_type;
362 using real_type = typename numbers::NumberTraits<value_type>::real_type;
363 using pointer = value_type *;
364 using const_pointer = const value_type *;
365 using reference = value_type &;
366 using const_reference = const value_type &;
367 using size_type = types::global_dof_index;
368
369 using iterator =
370 MatrixIterator<BlockMatrixIterators::Accessor<BlockMatrixBase, false>>;
371
372 using const_iterator =
373 MatrixIterator<BlockMatrixIterators::Accessor<BlockMatrixBase, true>>;
374
375
379 BlockMatrixBase() = default;
380
384 ~BlockMatrixBase() override;
385
402 template <typename BlockMatrixType>
403 BlockMatrixBase &
404 copy_from(const BlockMatrixType &source);
405
409 BlockType &
410 block(const unsigned int row, const unsigned int column);
411
412
417 const BlockType &
418 block(const unsigned int row, const unsigned int column) const;
419
424 size_type
425 m() const;
426
431 size_type
432 n() const;
433
434
439 unsigned int
440 n_block_rows() const;
441
446 unsigned int
447 n_block_cols() const;
448
454 void
455 set(const size_type i, const size_type j, const value_type value);
456
472 template <typename number>
473 void
474 set(const std::vector<size_type> &indices,
475 const FullMatrix<number> &full_matrix,
476 const bool elide_zero_values = false);
477
483 template <typename number>
484 void
485 set(const std::vector<size_type> &row_indices,
486 const std::vector<size_type> &col_indices,
487 const FullMatrix<number> &full_matrix,
488 const bool elide_zero_values = false);
489
500 template <typename number>
501 void
502 set(const size_type row,
503 const std::vector<size_type> &col_indices,
504 const std::vector<number> &values,
505 const bool elide_zero_values = false);
506
516 template <typename number>
517 void
518 set(const size_type row,
519 const size_type n_cols,
520 const size_type *col_indices,
521 const number *values,
522 const bool elide_zero_values = false);
523
529 void
530 add(const size_type i, const size_type j, const value_type value);
531
546 template <typename number>
547 void
548 add(const std::vector<size_type> &indices,
549 const FullMatrix<number> &full_matrix,
550 const bool elide_zero_values = true);
551
557 template <typename number>
558 void
559 add(const std::vector<size_type> &row_indices,
560 const std::vector<size_type> &col_indices,
561 const FullMatrix<number> &full_matrix,
562 const bool elide_zero_values = true);
563
573 template <typename number>
574 void
575 add(const size_type row,
576 const std::vector<size_type> &col_indices,
577 const std::vector<number> &values,
578 const bool elide_zero_values = true);
579
589 template <typename number>
590 void
591 add(const size_type row,
592 const size_type n_cols,
593 const size_type *col_indices,
594 const number *values,
595 const bool elide_zero_values = true,
596 const bool col_indices_are_sorted = false);
597
609 void
610 add(const value_type factor, const BlockMatrixBase<MatrixType> &matrix);
611
618 value_type
619 operator()(const size_type i, const size_type j) const;
620
629 value_type
630 el(const size_type i, const size_type j) const;
631
642 value_type
643 diag_element(const size_type i) const;
644
653 void
654 compress(VectorOperation::values operation);
655
659 BlockMatrixBase &
660 operator*=(const value_type factor);
661
665 BlockMatrixBase &
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>
698 value_type
699 matrix_norm_square(const BlockVectorType &v) const;
700
705 real_type
706 frobenius_norm() const;
707
711 template <typename BlockVectorType>
712 value_type
713 matrix_scalar_product(const BlockVectorType &u,
714 const BlockVectorType &v) const;
715
719 template <typename BlockVectorType>
720 value_type
721 residual(BlockVectorType &dst,
722 const BlockVectorType &x,
723 const BlockVectorType &b) const;
724
731 void
732 print(std::ostream &out, const bool alternative_output = false) const;
733
737 iterator
738 begin();
739
743 iterator
744 end();
745
749 iterator
750 begin(const size_type r);
751
755 iterator
756 end(const size_type r);
760 const_iterator
761 begin() const;
762
766 const_iterator
767 end() const;
768
772 const_iterator
773 begin(const size_type r) const;
774
778 const_iterator
779 end(const size_type r) const;
780
784 const BlockIndices &
785 get_row_indices() const;
786
790 const BlockIndices &
791 get_column_indices() const;
792
798 std::size_t
799 memory_consumption() const;
800
809 DeclException4(ExcIncompatibleRowNumbers,
810 int,
811 int,
812 int,
813 int,
814 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
815 << ',' << arg4 << "] have differing row numbers.");
819 DeclException4(ExcIncompatibleColNumbers,
820 int,
821 int,
822 int,
823 int,
824 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
825 << ',' << arg4 << "] have differing column numbers.");
827protected:
840 void
841 clear();
842
846 BlockIndices row_block_indices;
847 BlockIndices column_block_indices;
848
852 Table<2, ObserverPointer<BlockType, BlockMatrixBase<MatrixType>>> sub_objects;
853
872 void
873 collect_sizes();
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
997 prepare_add_operation();
998
1003 void
1004 prepare_set_operation();
1005
1006
1007private:
1016 struct TemporaryData
1017 {
1022 std::vector<size_type> counter_within_block;
1023
1028 std::vector<std::vector<size_type>> column_indices;
1029
1034 std::vector<std::vector<value_type>> column_values;
1035
1040 Threads::Mutex mutex;
1041
1051 TemporaryData &
1052 operator=(const TemporaryData &)
1053 {
1054 return *this;
1055 }
1056 };
1057
1064 TemporaryData temporary_data;
1065
1066 // Make the iterator class a friend. We have to work around a compiler bug
1067 // here again.
1068 template <typename, bool>
1069 friend class BlockMatrixIterators::Accessor;
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>
1085 inline AccessorBase<BlockMatrixType>::AccessorBase()
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 {
1095 Assert(row_block != numbers::invalid_unsigned_int, ExcIteratorPastEnd());
1096
1097 return row_block;
1098 }
1099
1100
1101 template <typename BlockMatrixType>
1102 inline unsigned int
1103 AccessorBase<BlockMatrixType>::block_column() const
1104 {
1105 Assert(col_block != numbers::invalid_unsigned_int, ExcIteratorPastEnd());
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,
1191 ExcIteratorPastEnd());
1192
1193 return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1194 base_iterator->row());
1195 }
1196
1197
1198 template <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,
1203 ExcIteratorPastEnd());
1204
1205 return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
1206 base_iterator->column());
1207 }
1208
1209
1210 template <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,
1215 ExcIteratorPastEnd());
1216 Assert(this->col_block != numbers::invalid_unsigned_int,
1217 ExcIteratorPastEnd());
1218
1219 return base_iterator->value();
1220 }
1221
1222
1223
1224 template <typename BlockMatrixType>
1225 inline void
1226 Accessor<BlockMatrixType, true>::advance()
1227 {
1228 Assert(this->row_block != numbers::invalid_unsigned_int,
1229 ExcIteratorPastEnd());
1230 Assert(this->col_block != numbers::invalid_unsigned_int,
1231 ExcIteratorPastEnd());
1232
1233 // Remember current row inside block
1234 size_type local_row = base_iterator->row();
1235
1236 // Advance one element inside the
1237 // current block
1238 ++base_iterator;
1239
1240 // while we hit the end of the row of a
1241 // block (which may happen multiple
1242 // times if rows inside a block are
1243 // empty), we have to jump to the next
1244 // block and take the
1245 while (base_iterator ==
1246 matrix->block(this->row_block, this->col_block).end(local_row))
1247 {
1248 // jump to next block in this block
1249 // row, if possible, otherwise go
1250 // to next row
1251 if (this->col_block < matrix->n_block_cols() - 1)
1252 {
1253 ++this->col_block;
1254 base_iterator =
1255 matrix->block(this->row_block, this->col_block).begin(local_row);
1256 }
1257 else
1258 {
1259 // jump back to next row in
1260 // first block column
1261 this->col_block = 0;
1262 ++local_row;
1263
1264 // see if this has brought us
1265 // past the number of rows in
1266 // this block. if so see
1267 // whether we've just fallen
1268 // off the end of the whole
1269 // matrix
1270 if (local_row ==
1271 matrix->block(this->row_block, this->col_block).m())
1272 {
1273 local_row = 0;
1274 ++this->row_block;
1275 if (this->row_block == matrix->n_block_rows())
1276 {
1277 this->row_block = numbers::invalid_unsigned_int;
1278 this->col_block = numbers::invalid_unsigned_int;
1279 return;
1280 }
1281 }
1282
1283 base_iterator =
1284 matrix->block(this->row_block, this->col_block).begin(local_row);
1285 }
1286 }
1287 }
1288
1289
1290 template <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 {
1365 Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1366
1367 return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1368 base_iterator->row());
1369 }
1370
1371
1372 template <typename BlockMatrixType>
1373 inline typename Accessor<BlockMatrixType, false>::size_type
1374 Accessor<BlockMatrixType, false>::column() const
1375 {
1376 Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1377
1378 return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
1379 base_iterator->column());
1380 }
1381
1382
1383 template <typename BlockMatrixType>
1384 inline typename Accessor<BlockMatrixType, false>::value_type
1385 Accessor<BlockMatrixType, false>::value() const
1386 {
1387 Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1388 Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1389
1390 return base_iterator->value();
1391 }
1392
1393
1394
1395 template <typename BlockMatrixType>
1396 inline void
1397 Accessor<BlockMatrixType, false>::set_value(
1398 typename Accessor<BlockMatrixType, false>::value_type newval) const
1399 {
1400 Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1401 Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1402
1403 base_iterator->value() = newval;
1404 }
1405
1406
1407
1408 template <typename BlockMatrixType>
1409 inline void
1410 Accessor<BlockMatrixType, false>::advance()
1411 {
1412 Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1413 Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());
1414
1415 // Remember current row inside block
1416 size_type local_row = base_iterator->row();
1417
1418 // Advance one element inside the
1419 // current block
1420 ++base_iterator;
1421
1422 // while we hit the end of the row of a
1423 // block (which may happen multiple
1424 // times if rows inside a block are
1425 // empty), we have to jump to the next
1426 // block and take the
1427 while (base_iterator ==
1428 matrix->block(this->row_block, this->col_block).end(local_row))
1429 {
1430 // jump to next block in this block
1431 // row, if possible, otherwise go
1432 // to next row
1433 if (this->col_block < matrix->n_block_cols() - 1)
1434 {
1435 ++this->col_block;
1436 base_iterator =
1437 matrix->block(this->row_block, this->col_block).begin(local_row);
1438 }
1439 else
1440 {
1441 // jump back to next row in
1442 // first block column
1443 this->col_block = 0;
1444 ++local_row;
1445
1446 // see if this has brought us
1447 // past the number of rows in
1448 // this block. if so see
1449 // whether we've just fallen
1450 // off the end of the whole
1451 // matrix
1452 if (local_row ==
1453 matrix->block(this->row_block, this->col_block).m())
1454 {
1455 local_row = 0;
1456 ++this->row_block;
1457 if (this->row_block == matrix->n_block_rows())
1458 {
1459 this->row_block = numbers::invalid_size_type;
1460 this->col_block = numbers::invalid_size_type;
1461 return;
1462 }
1463 }
1464
1465 base_iterator =
1466 matrix->block(this->row_block, this->col_block).begin(local_row);
1467 }
1468 }
1469 }
1470
1471
1472
1473 template <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>
1497inline BlockMatrixBase<MatrixType>::~BlockMatrixBase()
1498{
1499 try
1500 {
1501 clear();
1502 }
1503 catch (...)
1504 {}
1505}
1506
1507
1508template <typename MatrixType>
1509template <typename BlockMatrixType>
1510inline BlockMatrixBase<MatrixType> &
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
1523BlockMatrixBase<MatrixType>::memory_consumption() const
1524{
1525 std::size_t mem =
1526 MemoryConsumption::memory_consumption(row_block_indices) +
1527 MemoryConsumption::memory_consumption(column_block_indices) +
1528 MemoryConsumption::memory_consumption(sub_objects) +
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];
1538 mem += MemoryConsumption::memory_consumption(*p);
1539 }
1540
1541 return mem;
1542}
1543
1544
1545
1546template <typename MatrixType>
1547inline void
1548BlockMatrixBase<MatrixType>::clear()
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>
1566inline typename BlockMatrixBase<MatrixType>::BlockType &
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>
1591inline typename BlockMatrixBase<MatrixType>::size_type
1592BlockMatrixBase<MatrixType>::m() const
1593{
1594 return row_block_indices.total_size();
1595}
1596
1597
1598
1599template <typename MatrixType>
1600inline typename BlockMatrixBase<MatrixType>::size_type
1601BlockMatrixBase<MatrixType>::n() const
1602{
1603 return column_block_indices.total_size();
1604}
1605
1606
1607
1608template <typename MatrixType>
1609inline unsigned int
1610BlockMatrixBase<MatrixType>::n_block_cols() const
1611{
1612 return column_block_indices.size();
1613}
1614
1615
1616
1617template <typename MatrixType>
1618inline unsigned int
1619BlockMatrixBase<MatrixType>::n_block_rows() const
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# ifdef DEBUG
1790 // If in debug mode, do a check whether
1791 // the right length has been obtained.
1792 size_type length = 0;
1793 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
1794 length += temporary_data.counter_within_block[i];
1795 Assert(length <= n_cols, ExcInternalError());
1796# endif
1797
1798 // Now we found out about where the
1799 // individual columns should start and
1800 // where we should start reading out
1801 // data. Now let's write the data into
1802 // the individual blocks!
1803 const std::pair<unsigned int, size_type> row_index =
1804 this->row_block_indices.global_to_local(row);
1805 for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
1806 {
1807 if (temporary_data.counter_within_block[block_col] == 0)
1808 continue;
1809
1810 block(row_index.first, block_col)
1811 .set(row_index.second,
1812 temporary_data.counter_within_block[block_col],
1813 temporary_data.column_indices[block_col].data(),
1814 temporary_data.column_values[block_col].data(),
1815 false);
1816 }
1817}
1818
1819
1820
1821template <typename MatrixType>
1822inline void
1823BlockMatrixBase<MatrixType>::add(const size_type i,
1824 const size_type j,
1825 const value_type value)
1826{
1827 AssertIsFinite(value);
1828
1829 prepare_add_operation();
1830
1831 // save some cycles for zero additions, but
1832 // only if it is safe for the matrix we are
1833 // working with
1834 using MatrixTraits = typename MatrixType::Traits;
1835 if ((MatrixTraits::zero_addition_can_be_elided == true) &&
1836 (value == value_type()))
1837 return;
1838
1839 const std::pair<unsigned int, size_type>
1840 row_index = row_block_indices.global_to_local(i),
1841 col_index = column_block_indices.global_to_local(j);
1842 block(row_index.first, col_index.first)
1843 .add(row_index.second, col_index.second, value);
1844}
1845
1846
1847
1848template <typename MatrixType>
1849template <typename number>
1850inline void
1851BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &row_indices,
1852 const std::vector<size_type> &col_indices,
1853 const FullMatrix<number> &values,
1854 const bool elide_zero_values)
1855{
1856 Assert(row_indices.size() == values.m(),
1857 ExcDimensionMismatch(row_indices.size(), values.m()));
1858 Assert(col_indices.size() == values.n(),
1859 ExcDimensionMismatch(col_indices.size(), values.n()));
1860
1861 for (size_type i = 0; i < row_indices.size(); ++i)
1862 add(row_indices[i],
1863 col_indices.size(),
1864 col_indices.data(),
1865 &values(i, 0),
1866 elide_zero_values);
1867}
1868
1869
1870
1871template <typename MatrixType>
1872template <typename number>
1873inline void
1874BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &indices,
1875 const FullMatrix<number> &values,
1876 const bool elide_zero_values)
1877{
1878 Assert(indices.size() == values.m(),
1879 ExcDimensionMismatch(indices.size(), values.m()));
1880 Assert(values.n() == values.m(), ExcNotQuadratic());
1881
1882 for (size_type i = 0; i < indices.size(); ++i)
1883 add(indices[i],
1884 indices.size(),
1885 indices.data(),
1886 &values(i, 0),
1887 elide_zero_values);
1888}
1889
1890
1891
1892template <typename MatrixType>
1893template <typename number>
1894inline void
1895BlockMatrixBase<MatrixType>::add(const size_type row,
1896 const std::vector<size_type> &col_indices,
1897 const std::vector<number> &values,
1898 const bool elide_zero_values)
1899{
1900 Assert(col_indices.size() == values.size(),
1901 ExcDimensionMismatch(col_indices.size(), values.size()));
1902
1903 add(row,
1904 col_indices.size(),
1905 col_indices.data(),
1906 values.data(),
1907 elide_zero_values);
1908}
1909
1910
1911
1912// This is a very messy function, since
1913// we need to calculate to each position
1914// the location in the global array.
1915template <typename MatrixType>
1916template <typename number>
1917inline void
1918BlockMatrixBase<MatrixType>::add(const size_type row,
1919 const size_type n_cols,
1920 const size_type *col_indices,
1921 const number *values,
1922 const bool elide_zero_values,
1923 const bool col_indices_are_sorted)
1924{
1925 prepare_add_operation();
1926
1927 // TODO: Look over this to find out
1928 // whether we can do that more
1929 // efficiently.
1930 if (col_indices_are_sorted == true)
1931 {
1932# ifdef DEBUG
1933 // check whether indices really are
1934 // sorted.
1935 size_type before = col_indices[0];
1936 for (size_type i = 1; i < n_cols; ++i)
1937 if (col_indices[i] <= before)
1938 {
1939 Assert(false,
1940 ExcMessage("Flag col_indices_are_sorted is set, but "
1941 "indices appear to not be sorted."));
1942 }
1943 else
1944 before = col_indices[i];
1945# endif
1946 const std::pair<unsigned int, size_type> row_index =
1947 this->row_block_indices.global_to_local(row);
1948
1949 if (this->n_block_cols() > 1)
1950 {
1951 const size_type *first_block =
1952 Utilities::lower_bound(col_indices,
1953 col_indices + n_cols,
1954 this->column_block_indices.block_start(1));
1955
1956 const size_type n_zero_block_indices = first_block - col_indices;
1957 block(row_index.first, 0)
1958 .add(row_index.second,
1959 n_zero_block_indices,
1960 col_indices,
1961 values,
1962 elide_zero_values,
1963 col_indices_are_sorted);
1964
1965 if (n_zero_block_indices < n_cols)
1966 this->add(row,
1967 n_cols - n_zero_block_indices,
1968 first_block,
1969 values + n_zero_block_indices,
1970 elide_zero_values,
1971 false);
1972 }
1973 else
1974 {
1975 block(row_index.first, 0)
1976 .add(row_index.second,
1977 n_cols,
1978 col_indices,
1979 values,
1980 elide_zero_values,
1981 col_indices_are_sorted);
1982 }
1983
1984 return;
1985 }
1986
1987 // Lock scratch arrays, then resize them
1988 std::lock_guard<std::mutex> lock(temporary_data.mutex);
1989
1990 if (temporary_data.column_indices.size() < this->n_block_cols())
1991 {
1992 temporary_data.column_indices.resize(this->n_block_cols());
1993 temporary_data.column_values.resize(this->n_block_cols());
1994 temporary_data.counter_within_block.resize(this->n_block_cols());
1995 }
1996
1997 // Resize sub-arrays to n_cols. This
1998 // is a bit wasteful, but we resize
1999 // only a few times (then the maximum
2000 // row length won't increase that
2001 // much any more). At least we know
2002 // that all arrays are going to be of
2003 // the same size, so we can check
2004 // whether the size of one is large
2005 // enough before actually going
2006 // through all of them.
2007 if (temporary_data.column_indices[0].size() < n_cols)
2008 {
2009 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2010 {
2011 temporary_data.column_indices[i].resize(n_cols);
2012 temporary_data.column_values[i].resize(n_cols);
2013 }
2014 }
2015
2016 // Reset the number of added elements
2017 // in each block to zero.
2018 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2019 temporary_data.counter_within_block[i] = 0;
2020
2021 // Go through the column indices to
2022 // find out which portions of the
2023 // values should be written into
2024 // which block of the matrix. We need
2025 // to touch all the data, since we
2026 // can't be sure that the data of one
2027 // block is stored contiguously (in
2028 // fact, data will be intermixed when
2029 // it comes from an element matrix).
2030 for (size_type j = 0; j < n_cols; ++j)
2031 {
2032 number value = values[j];
2033
2034 if (value == number() && elide_zero_values == true)
2035 continue;
2036
2037 const std::pair<unsigned int, size_type> col_index =
2038 this->column_block_indices.global_to_local(col_indices[j]);
2039
2040 const size_type local_index =
2041 temporary_data.counter_within_block[col_index.first]++;
2042
2043 temporary_data.column_indices[col_index.first][local_index] =
2044 col_index.second;
2045 temporary_data.column_values[col_index.first][local_index] = value;
2046 }
2047
2048# ifdef DEBUG
2049 // If in debug mode, do a check whether
2050 // the right length has been obtained.
2051 size_type length = 0;
2052 for (unsigned int i = 0; i < this->n_block_cols(); ++i)
2053 length += temporary_data.counter_within_block[i];
2054 Assert(length <= n_cols, ExcInternalError());
2055# endif
2056
2057 // Now we found out about where the
2058 // individual columns should start and
2059 // where we should start reading out
2060 // data. Now let's write the data into
2061 // the individual blocks!
2062 const std::pair<unsigned int, size_type> row_index =
2063 this->row_block_indices.global_to_local(row);
2064 for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
2065 {
2066 if (temporary_data.counter_within_block[block_col] == 0)
2067 continue;
2068
2069 block(row_index.first, block_col)
2070 .add(row_index.second,
2071 temporary_data.counter_within_block[block_col],
2072 temporary_data.column_indices[block_col].data(),
2073 temporary_data.column_values[block_col].data(),
2074 false,
2075 col_indices_are_sorted);
2076 }
2077}
2078
2079
2080
2081template <typename MatrixType>
2082inline void
2083BlockMatrixBase<MatrixType>::add(const value_type factor,
2084 const BlockMatrixBase<MatrixType> &matrix)
2085{
2086 AssertIsFinite(factor);
2087
2088 prepare_add_operation();
2089
2090 // save some cycles for zero additions, but
2091 // only if it is safe for the matrix we are
2092 // working with
2093 using MatrixTraits = typename MatrixType::Traits;
2094 if ((MatrixTraits::zero_addition_can_be_elided == true) && (factor == 0))
2095 return;
2096
2097 for (unsigned int row = 0; row < n_block_rows(); ++row)
2098 for (unsigned int col = 0; col < n_block_cols(); ++col)
2099 // This function should throw if the sparsity
2100 // patterns of the two blocks differ
2101 block(row, col).add(factor, matrix.block(row, col));
2102}
2103
2104
2105
2106template <typename MatrixType>
2107inline typename BlockMatrixBase<MatrixType>::value_type
2108BlockMatrixBase<MatrixType>::operator()(const size_type i,
2109 const size_type j) const
2110{
2111 const std::pair<unsigned int, size_type>
2112 row_index = row_block_indices.global_to_local(i),
2113 col_index = column_block_indices.global_to_local(j);
2114 return block(row_index.first, col_index.first)(row_index.second,
2115 col_index.second);
2116}
2117
2118
2119
2120template <typename MatrixType>
2121inline typename BlockMatrixBase<MatrixType>::value_type
2122BlockMatrixBase<MatrixType>::el(const size_type i, const size_type j) const
2123{
2124 const std::pair<unsigned int, size_type>
2125 row_index = row_block_indices.global_to_local(i),
2126 col_index = column_block_indices.global_to_local(j);
2127 return block(row_index.first, col_index.first)
2128 .el(row_index.second, col_index.second);
2129}
2130
2131
2132
2133template <typename MatrixType>
2134inline typename BlockMatrixBase<MatrixType>::value_type
2135BlockMatrixBase<MatrixType>::diag_element(const size_type i) const
2136{
2137 Assert(n_block_rows() == n_block_cols(), ExcNotQuadratic());
2138
2139 const std::pair<unsigned int, size_type> index =
2140 row_block_indices.global_to_local(i);
2141 return block(index.first, index.first).diag_element(index.second);
2142}
2143
2144
2145
2146template <typename MatrixType>
2147inline void
2148BlockMatrixBase<MatrixType>::compress(VectorOperation::values operation)
2149{
2150 for (unsigned int r = 0; r < n_block_rows(); ++r)
2151 for (unsigned int c = 0; c < n_block_cols(); ++c)
2152 block(r, c).compress(operation);
2153}
2154
2155
2156
2157template <typename MatrixType>
2158inline BlockMatrixBase<MatrixType> &
2159BlockMatrixBase<MatrixType>::operator*=(const value_type factor)
2160{
2161 Assert(n_block_cols() != 0, ExcNotInitialized());
2162 Assert(n_block_rows() != 0, ExcNotInitialized());
2163
2164 for (unsigned int r = 0; r < n_block_rows(); ++r)
2165 for (unsigned int c = 0; c < n_block_cols(); ++c)
2166 block(r, c) *= factor;
2167
2168 return *this;
2169}
2170
2171
2172
2173template <typename MatrixType>
2174inline BlockMatrixBase<MatrixType> &
2175BlockMatrixBase<MatrixType>::operator/=(const value_type factor)
2176{
2177 Assert(n_block_cols() != 0, ExcNotInitialized());
2178 Assert(n_block_rows() != 0, ExcNotInitialized());
2179 Assert(factor != 0, ExcDivideByZero());
2180
2181 const value_type factor_inv = 1. / factor;
2182
2183 for (unsigned int r = 0; r < n_block_rows(); ++r)
2184 for (unsigned int c = 0; c < n_block_cols(); ++c)
2185 block(r, c) *= factor_inv;
2186
2187 return *this;
2188}
2189
2190
2191
2192template <typename MatrixType>
2193const BlockIndices &
2194BlockMatrixBase<MatrixType>::get_row_indices() const
2195{
2196 return this->row_block_indices;
2197}
2198
2199
2200
2201template <typename MatrixType>
2202const BlockIndices &
2203BlockMatrixBase<MatrixType>::get_column_indices() const
2204{
2205 return this->column_block_indices;
2206}
2207
2208
2209
2210template <typename MatrixType>
2211template <typename BlockVectorType>
2212void
2213BlockMatrixBase<MatrixType>::vmult_block_block(BlockVectorType &dst,
2214 const BlockVectorType &src) const
2215{
2216 Assert(dst.n_blocks() == n_block_rows(),
2217 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2218 Assert(src.n_blocks() == n_block_cols(),
2219 ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2220
2221 for (size_type row = 0; row < n_block_rows(); ++row)
2222 {
2223 block(row, 0).vmult(dst.block(row), src.block(0));
2224 for (size_type col = 1; col < n_block_cols(); ++col)
2225 block(row, col).vmult_add(dst.block(row), src.block(col));
2226 };
2227}
2228
2229
2230
2231template <typename MatrixType>
2232template <typename BlockVectorType, typename VectorType>
2233void
2234BlockMatrixBase<MatrixType>::vmult_nonblock_block(
2235 VectorType &dst,
2236 const BlockVectorType &src) const
2237{
2238 Assert(n_block_rows() == 1, ExcDimensionMismatch(1, n_block_rows()));
2239 Assert(src.n_blocks() == n_block_cols(),
2240 ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2241
2242 block(0, 0).vmult(dst, src.block(0));
2243 for (size_type col = 1; col < n_block_cols(); ++col)
2244 block(0, col).vmult_add(dst, src.block(col));
2245}
2246
2247
2248
2249template <typename MatrixType>
2250template <typename BlockVectorType, typename VectorType>
2251void
2252BlockMatrixBase<MatrixType>::vmult_block_nonblock(BlockVectorType &dst,
2253 const VectorType &src) const
2254{
2255 Assert(dst.n_blocks() == n_block_rows(),
2256 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2257 Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2258
2259 for (size_type row = 0; row < n_block_rows(); ++row)
2260 block(row, 0).vmult(dst.block(row), src);
2261}
2262
2263
2264
2265template <typename MatrixType>
2266template <typename VectorType>
2267void
2268BlockMatrixBase<MatrixType>::vmult_nonblock_nonblock(
2269 VectorType &dst,
2270 const VectorType &src) const
2271{
2272 Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));
2273 Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2274
2275 block(0, 0).vmult(dst, src);
2276}
2277
2278
2279
2280template <typename MatrixType>
2281template <typename BlockVectorType>
2282void
2283BlockMatrixBase<MatrixType>::vmult_add(BlockVectorType &dst,
2284 const BlockVectorType &src) const
2285{
2286 Assert(dst.n_blocks() == n_block_rows(),
2287 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2288 Assert(src.n_blocks() == n_block_cols(),
2289 ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2290
2291 for (unsigned int row = 0; row < n_block_rows(); ++row)
2292 for (unsigned int col = 0; col < n_block_cols(); ++col)
2293 block(row, col).vmult_add(dst.block(row), src.block(col));
2294}
2295
2296
2297
2298template <typename MatrixType>
2299template <typename BlockVectorType>
2300void
2301BlockMatrixBase<MatrixType>::Tvmult_block_block(
2302 BlockVectorType &dst,
2303 const BlockVectorType &src) const
2304{
2305 Assert(dst.n_blocks() == n_block_cols(),
2306 ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2307 Assert(src.n_blocks() == n_block_rows(),
2308 ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2309
2310 dst = 0.;
2311
2312 for (unsigned int row = 0; row < n_block_rows(); ++row)
2313 {
2314 for (unsigned int col = 0; col < n_block_cols(); ++col)
2315 block(row, col).Tvmult_add(dst.block(col), src.block(row));
2316 };
2317}
2318
2319
2320
2321template <typename MatrixType>
2322template <typename BlockVectorType, typename VectorType>
2323void
2324BlockMatrixBase<MatrixType>::Tvmult_block_nonblock(BlockVectorType &dst,
2325 const VectorType &src) const
2326{
2327 Assert(dst.n_blocks() == n_block_cols(),
2328 ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2329 Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));
2330
2331 dst = 0.;
2332
2333 for (unsigned int col = 0; col < n_block_cols(); ++col)
2334 block(0, col).Tvmult_add(dst.block(col), src);
2335}
2336
2337
2338
2339template <typename MatrixType>
2340template <typename BlockVectorType, typename VectorType>
2341void
2342BlockMatrixBase<MatrixType>::Tvmult_nonblock_block(
2343 VectorType &dst,
2344 const BlockVectorType &src) const
2345{
2346 Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2347 Assert(src.n_blocks() == n_block_rows(),
2348 ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2349
2350 block(0, 0).Tvmult(dst, src.block(0));
2351
2352 for (size_type row = 1; row < n_block_rows(); ++row)
2353 block(row, 0).Tvmult_add(dst, src.block(row));
2354}
2355
2356
2357
2358template <typename MatrixType>
2359template <typename VectorType>
2360void
2361BlockMatrixBase<MatrixType>::Tvmult_nonblock_nonblock(
2362 VectorType &dst,
2363 const VectorType &src) const
2364{
2365 Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
2366 Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));
2367
2368 block(0, 0).Tvmult(dst, src);
2369}
2370
2371
2372
2373template <typename MatrixType>
2374template <typename BlockVectorType>
2375void
2376BlockMatrixBase<MatrixType>::Tvmult_add(BlockVectorType &dst,
2377 const BlockVectorType &src) const
2378{
2379 Assert(dst.n_blocks() == n_block_cols(),
2380 ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2381 Assert(src.n_blocks() == n_block_rows(),
2382 ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2383
2384 for (unsigned int row = 0; row < n_block_rows(); ++row)
2385 for (unsigned int col = 0; col < n_block_cols(); ++col)
2386 block(row, col).Tvmult_add(dst.block(col), src.block(row));
2387}
2388
2389
2390
2391template <typename MatrixType>
2392template <typename BlockVectorType>
2393typename BlockMatrixBase<MatrixType>::value_type
2394BlockMatrixBase<MatrixType>::matrix_norm_square(const BlockVectorType &v) const
2395{
2396 Assert(n_block_rows() == n_block_cols(), ExcNotQuadratic());
2397 Assert(v.n_blocks() == n_block_rows(),
2398 ExcDimensionMismatch(v.n_blocks(), n_block_rows()));
2399
2400 value_type norm_sqr = 0;
2401 for (unsigned int row = 0; row < n_block_rows(); ++row)
2402 for (unsigned int col = 0; col < n_block_cols(); ++col)
2403 if (row == col)
2404 norm_sqr += block(row, col).matrix_norm_square(v.block(row));
2405 else
2406 norm_sqr +=
2407 block(row, col).matrix_scalar_product(v.block(row), v.block(col));
2408 return norm_sqr;
2409}
2410
2411
2412
2413template <typename MatrixType>
2414typename BlockMatrixBase<MatrixType>::real_type
2415BlockMatrixBase<MatrixType>::frobenius_norm() const
2416{
2417 value_type norm_sqr = 0;
2418
2419 // For each block, get the Frobenius norm, and add the square to the
2420 // accumulator for the full matrix
2421 for (unsigned int row = 0; row < n_block_rows(); ++row)
2422 {
2423 for (unsigned int col = 0; col < n_block_cols(); ++col)
2424 {
2425 const value_type block_norm = block(row, col).frobenius_norm();
2426 norm_sqr += block_norm * block_norm;
2427 }
2428 }
2429
2430 return std::sqrt(norm_sqr);
2431}
2432
2433
2434
2435template <typename MatrixType>
2436template <typename BlockVectorType>
2437typename BlockMatrixBase<MatrixType>::value_type
2438BlockMatrixBase<MatrixType>::matrix_scalar_product(
2439 const BlockVectorType &u,
2440 const BlockVectorType &v) const
2441{
2442 Assert(u.n_blocks() == n_block_rows(),
2443 ExcDimensionMismatch(u.n_blocks(), n_block_rows()));
2444 Assert(v.n_blocks() == n_block_cols(),
2445 ExcDimensionMismatch(v.n_blocks(), n_block_cols()));
2446
2447 value_type result = 0;
2448 for (unsigned int row = 0; row < n_block_rows(); ++row)
2449 for (unsigned int col = 0; col < n_block_cols(); ++col)
2450 result +=
2451 block(row, col).matrix_scalar_product(u.block(row), v.block(col));
2452 return result;
2453}
2454
2455
2456
2457template <typename MatrixType>
2458template <typename BlockVectorType>
2459typename BlockMatrixBase<MatrixType>::value_type
2460BlockMatrixBase<MatrixType>::residual(BlockVectorType &dst,
2461 const BlockVectorType &x,
2462 const BlockVectorType &b) const
2463{
2464 Assert(dst.n_blocks() == n_block_rows(),
2465 ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2466 Assert(b.n_blocks() == n_block_rows(),
2467 ExcDimensionMismatch(b.n_blocks(), n_block_rows()));
2468 Assert(x.n_blocks() == n_block_cols(),
2469 ExcDimensionMismatch(x.n_blocks(), n_block_cols()));
2470 // in block notation, the residual is
2471 // r_i = b_i - \sum_j A_ij x_j.
2472 // this can be written as
2473 // r_i = b_i - A_i0 x_0 - \sum_{j>0} A_ij x_j.
2474 //
2475 // for the first two terms, we can
2476 // call the residual function of
2477 // A_i0. for the other terms, we
2478 // use vmult_add. however, we want
2479 // to subtract, so in order to
2480 // avoid a temporary vector, we
2481 // perform a sign change of the
2482 // first two term before, and after
2483 // adding up
2484 for (unsigned int row = 0; row < n_block_rows(); ++row)
2485 {
2486 block(row, 0).residual(dst.block(row), x.block(0), b.block(row));
2487
2488 for (size_type i = 0; i < dst.block(row).size(); ++i)
2489 dst.block(row)(i) = -dst.block(row)(i);
2490
2491 for (unsigned int col = 1; col < n_block_cols(); ++col)
2492 block(row, col).vmult_add(dst.block(row), x.block(col));
2493
2494 for (size_type i = 0; i < dst.block(row).size(); ++i)
2495 dst.block(row)(i) = -dst.block(row)(i);
2496 };
2497
2498 value_type res = 0;
2499 for (size_type row = 0; row < n_block_rows(); ++row)
2500 res += dst.block(row).norm_sqr();
2501 return std::sqrt(res);
2502}
2503
2504
2505
2506template <typename MatrixType>
2507inline void
2508BlockMatrixBase<MatrixType>::print(std::ostream &out,
2509 const bool alternative_output) const
2510{
2511 for (unsigned int row = 0; row < n_block_rows(); ++row)
2512 for (unsigned int col = 0; col < n_block_cols(); ++col)
2513 {
2514 if (!alternative_output)
2515 out << "Block (" << row << ", " << col << ')' << std::endl;
2516
2517 block(row, col).print(out, alternative_output);
2518 }
2519}
2520
2521
2522
2523template <typename MatrixType>
2524inline typename BlockMatrixBase<MatrixType>::const_iterator
2525BlockMatrixBase<MatrixType>::begin() const
2526{
2527 return const_iterator(this, 0);
2528}
2529
2530
2531
2532template <typename MatrixType>
2533inline typename BlockMatrixBase<MatrixType>::const_iterator
2534BlockMatrixBase<MatrixType>::end() const
2535{
2536 return const_iterator(this, m());
2537}
2538
2539
2540
2541template <typename MatrixType>
2542inline typename BlockMatrixBase<MatrixType>::const_iterator
2543BlockMatrixBase<MatrixType>::begin(const size_type r) const
2544{
2545 AssertIndexRange(r, m());
2546 return const_iterator(this, r);
2547}
2548
2549
2550
2551template <typename MatrixType>
2552inline typename BlockMatrixBase<MatrixType>::const_iterator
2553BlockMatrixBase<MatrixType>::end(const size_type r) const
2554{
2555 AssertIndexRange(r, m());
2556 return const_iterator(this, r + 1);
2557}
2558
2559
2560
2561template <typename MatrixType>
2562inline typename BlockMatrixBase<MatrixType>::iterator
2563BlockMatrixBase<MatrixType>::begin()
2564{
2565 return iterator(this, 0);
2566}
2567
2568
2569
2570template <typename MatrixType>
2571inline typename BlockMatrixBase<MatrixType>::iterator
2572BlockMatrixBase<MatrixType>::end()
2573{
2574 return iterator(this, m());
2575}
2576
2577
2578
2579template <typename MatrixType>
2580inline typename BlockMatrixBase<MatrixType>::iterator
2581BlockMatrixBase<MatrixType>::begin(const size_type r)
2582{
2583 AssertIndexRange(r, m());
2584 return iterator(this, r);
2585}
2586
2587
2588
2589template <typename MatrixType>
2590inline typename BlockMatrixBase<MatrixType>::iterator
2591BlockMatrixBase<MatrixType>::end(const size_type r)
2592{
2593 AssertIndexRange(r, m());
2594 return iterator(this, r + 1);
2595}
2596
2597
2598
2599template <typename MatrixType>
2600void
2601BlockMatrixBase<MatrixType>::collect_sizes()
2602{
2603 std::vector<size_type> row_sizes(this->n_block_rows());
2604 std::vector<size_type> col_sizes(this->n_block_cols());
2605
2606 // first find out the row sizes
2607 // from the first block column
2608 for (unsigned int r = 0; r < this->n_block_rows(); ++r)
2609 row_sizes[r] = sub_objects[r][0]->m();
2610 // then check that the following
2611 // block columns have the same
2612 // sizes
2613 for (unsigned int c = 1; c < this->n_block_cols(); ++c)
2614 for (unsigned int r = 0; r < this->n_block_rows(); ++r)
2615 Assert(row_sizes[r] == sub_objects[r][c]->m(),
2616 ExcIncompatibleRowNumbers(r, 0, r, c));
2617
2618 // finally initialize the row
2619 // indices with this array
2620 this->row_block_indices.reinit(row_sizes);
2621
2622
2623 // then do the same with the columns
2624 for (unsigned int c = 0; c < this->n_block_cols(); ++c)
2625 col_sizes[c] = sub_objects[0][c]->n();
2626 for (unsigned int r = 1; r < this->n_block_rows(); ++r)
2627 for (unsigned int c = 0; c < this->n_block_cols(); ++c)
2628 Assert(col_sizes[c] == sub_objects[r][c]->n(),
2629 ExcIncompatibleRowNumbers(0, c, r, c));
2630
2631 // finally initialize the row
2632 // indices with this array
2633 this->column_block_indices.reinit(col_sizes);
2634}
2635
2636
2637
2638template <typename MatrixType>
2639void
2640BlockMatrixBase<MatrixType>::prepare_add_operation()
2641{
2642 for (unsigned int row = 0; row < n_block_rows(); ++row)
2643 for (unsigned int col = 0; col < n_block_cols(); ++col)
2644 block(row, col).prepare_add();
2645}
2646
2647
2648
2649template <typename MatrixType>
2650void
2651BlockMatrixBase<MatrixType>::prepare_set_operation()
2652{
2653 for (unsigned int row = 0; row < n_block_rows(); ++row)
2654 for (unsigned int col = 0; col < n_block_cols(); ++col)
2655 block(row, col).prepare_set();
2656}
2657
2658#endif // DOXYGEN
2659
2660
2661DEAL_II_NAMESPACE_CLOSE
2662
2663#endif // dealii_block_matrix_base_h
block_indices.h
BlockMatrixBase::Tvmult_nonblock_nonblock
void Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
BlockMatrixBase::vmult_nonblock_nonblock
void vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
BlockMatrixBase::add
void add(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number > &values, const bool elide_zero_values=true)
BlockMatrixBase::add
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)
BlockMatrixBase::reference
value_type & reference
Definition block_matrix_base.h:365
BlockMatrixBase::add
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)
BlockMatrixBase::print
void print(std::ostream &out, const bool alternative_output=false) const
BlockMatrixBase::Accessor
friend class BlockMatrixIterators::Accessor
Definition block_matrix_base.h:1069
BlockMatrixBase::end
const_iterator end() const
BlockMatrixBase::column_block_indices
BlockIndices column_block_indices
Definition block_matrix_base.h:847
BlockMatrixBase::operator*=
BlockMatrixBase & operator*=(const value_type factor)
BlockMatrixBase::matrix_norm_square
value_type matrix_norm_square(const BlockVectorType &v) const
BlockMatrixBase::Tvmult_block_block
void Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
BlockMatrixBase::const_reference
const value_type & const_reference
Definition block_matrix_base.h:366
BlockMatrixBase::Tvmult_nonblock_block
void Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
BlockMatrixBase::n_block_rows
unsigned int n_block_rows() const
BlockMatrixBase::operator()
value_type operator()(const size_type i, const size_type j) const
BlockMatrixBase::el
value_type el(const size_type i, const size_type j) const
BlockMatrixBase::begin
const_iterator begin() const
BlockMatrixBase::frobenius_norm
real_type frobenius_norm() const
BlockMatrixBase::vmult_block_block
void vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
BlockMatrixBase::compress
void compress(VectorOperation::values operation)
BlockMatrixBase::vmult_add
void vmult_add(BlockVectorType &dst, const BlockVectorType &src) const
BlockMatrixBase::size_type
types::global_dof_index size_type
Definition block_matrix_base.h:367
BlockMatrixBase::Tvmult_add
void Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const
BlockMatrixBase::value_type
typename BlockType::value_type value_type
Definition block_matrix_base.h:361
BlockMatrixBase::collect_sizes
void collect_sizes()
BlockMatrixBase::prepare_add_operation
void prepare_add_operation()
BlockMatrixBase::vmult_nonblock_block
void vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
BlockMatrixBase::residual
value_type residual(BlockVectorType &dst, const BlockVectorType &x, const BlockVectorType &b) const
BlockMatrixBase::BlockMatrixBase
BlockMatrixBase()=default
BlockMatrixBase::set
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)
BlockMatrixBase::n
size_type n() const
BlockMatrixBase::prepare_set_operation
void prepare_set_operation()
BlockMatrixBase::n_block_cols
unsigned int n_block_cols() const
BlockMatrixBase::const_pointer
const value_type * const_pointer
Definition block_matrix_base.h:364
BlockMatrixBase::begin
iterator begin()
BlockMatrixBase::temporary_data
TemporaryData temporary_data
Definition block_matrix_base.h:1064
BlockMatrixBase::get_row_indices
const BlockIndices & get_row_indices() const
BlockMatrixBase::end
iterator end()
BlockMatrixBase::set
void set(const size_type i, const size_type j, const value_type value)
BlockMatrixBase::copy_from
BlockMatrixBase & copy_from(const BlockMatrixType &source)
BlockMatrixBase::memory_consumption
std::size_t memory_consumption() const
BlockMatrixBase::add
void add(const value_type factor, const BlockMatrixBase< MatrixType > &matrix)
BlockMatrixBase::block
const BlockType & block(const unsigned int row, const unsigned int column) const
BlockMatrixBase::real_type
typename numbers::NumberTraits< value_type >::real_type real_type
Definition block_matrix_base.h:362
BlockMatrixBase::row_block_indices
BlockIndices row_block_indices
Definition block_matrix_base.h:846
BlockMatrixBase::~BlockMatrixBase
~BlockMatrixBase() override
BlockMatrixBase::block
BlockType & block(const unsigned int row, const unsigned int column)
BlockMatrixBase::m
size_type m() const
BlockMatrixBase::matrix_scalar_product
value_type matrix_scalar_product(const BlockVectorType &u, const BlockVectorType &v) const
BlockMatrixBase::get_column_indices
const BlockIndices & get_column_indices() const
BlockMatrixBase::set
void set(const std::vector< size_type > &indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=false)
BlockMatrixBase::add
void add(const std::vector< size_type > &indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=true)
BlockMatrixBase::add
void add(const size_type i, const size_type j, const value_type value)
BlockMatrixBase::set
void set(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number > &values, const bool elide_zero_values=false)
BlockMatrixBase::diag_element
value_type diag_element(const size_type i) const
BlockMatrixBase::vmult_block_nonblock
void vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
BlockMatrixBase::sub_objects
Table< 2, ObserverPointer< BlockType, BlockMatrixBase< MatrixType > > > sub_objects
Definition block_matrix_base.h:852
BlockMatrixBase::begin
iterator begin(const size_type r)
BlockMatrixBase::begin
const_iterator begin(const size_type r) const
BlockMatrixBase::operator/=
BlockMatrixBase & operator/=(const value_type factor)
BlockMatrixBase::Tvmult_block_nonblock
void Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
BlockMatrixBase::end
iterator end(const size_type r)
BlockMatrixBase::end
const_iterator end(const size_type r) const
BlockMatrixBase::set
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)
BlockMatrixIterators::AccessorBase::block_row
unsigned int block_row() const
BlockMatrixIterators::AccessorBase::block_column
unsigned int block_column() const
BlockMatrixIterators::AccessorBase::value_type
typename BlockMatrixType::value_type value_type
Definition block_matrix_base.h:73
BlockMatrixIterators::Accessor< BlockMatrixType, false >::Accessor
Accessor(BlockMatrixType *m, const size_type row, const size_type col)
BlockMatrixIterators::Accessor< BlockMatrixType, true >::Accessor
Accessor(const Accessor< BlockMatrixType, false > &)
BlockMatrixIterators::Accessor< BlockMatrixType, true >::base_iterator
BlockMatrixType::BlockType::const_iterator base_iterator
Definition block_matrix_base.h:267
BlockMatrixIterators::Accessor< BlockMatrixType, true >::Accessor
Accessor(const BlockMatrixType *m, const size_type row, const size_type col)
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
full_matrix.h
DeclException4
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:580
BlockMatrixBase::ExcIncompatibleColNumbers
static ::ExceptionBase & ExcIncompatibleColNumbers(int arg1, int arg2, int arg3, int arg4)
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
BlockMatrixBase::ExcIncompatibleRowNumbers
static ::ExceptionBase & ExcIncompatibleRowNumbers(int arg1, int arg2, int arg3, int arg4)
StandardExceptions::ExcIteratorPastEnd
static ::ExceptionBase & ExcIteratorPastEnd()
AssertIsFinite
#define AssertIsFinite(number)
Definition exceptions.h:2142
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2115
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDivideByZero
static ::ExceptionBase & ExcDivideByZero()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
LACExceptions::ExcNotQuadratic
static ::ExceptionBase & ExcNotQuadratic()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
matrix_iterator.h
memory_consumption.h
size
std::size_t size
Definition mpi.cc:734
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition lapack_support.h:57
LinearAlgebra::TpetraWrappers::TpetraTypes::MatrixType
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
Definition trilinos_tpetra_types.h:138
MemoryConsumption::memory_consumption
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
Definition memory_consumption.h:271
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition trilinos_sparse_matrix.cc:50
Utilities::lower_bound
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition utilities.h:1032
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:220
numbers::invalid_size_type
const types::global_dof_index invalid_size_type
Definition types.h:233
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6869
types::global_dof_index
unsigned int global_dof_index
Definition types.h:81
observer_pointer.h
BlockMatrixBase::TemporaryData::column_indices
std::vector< std::vector< size_type > > column_indices
Definition block_matrix_base.h:1028
BlockMatrixBase::TemporaryData::column_values
std::vector< std::vector< value_type > > column_values
Definition block_matrix_base.h:1034
BlockMatrixBase::TemporaryData::operator=
TemporaryData & operator=(const TemporaryData &)
Definition block_matrix_base.h:1052
BlockMatrixBase::TemporaryData::counter_within_block
std::vector< size_type > counter_within_block
Definition block_matrix_base.h:1022
vector_operation.h
Generated by doxygen 1.9.8