Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3067-g24f3489fdc 2025-04-14 21:10:00+00:00
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
block_vector_base.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2004 - 2024 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_vector_base_h
16#define dealii_block_vector_base_h
17
18
19#include <deal.II/base/config.h>
20
24
28#include <deal.II/lac/vector.h>
30
31#include <cmath>
32#include <cstddef>
33#include <iterator>
34#include <type_traits>
35#include <vector>
36
38
39
45namespace internal
46{
47 template <typename T>
48 using has_block_t = decltype(std::declval<const T>().block(0));
49
50 template <typename T>
51 constexpr bool has_block = internal::is_supported_operation<has_block_t, T>;
52
53 template <typename T>
54 using has_n_blocks_t = decltype(std::declval<const T>().n_blocks());
55
56 template <typename T>
57 constexpr bool has_n_blocks =
58 internal::is_supported_operation<has_n_blocks_t, T>;
59
60 template <typename T>
61 constexpr bool is_block_vector = has_block<T> && has_n_blocks<T>;
62} // namespace internal
63
78template <typename VectorType>
80{
81public:
87 static const bool value = internal::is_block_vector<VectorType>;
88};
89
90
91// instantiation of the static member
92template <typename VectorType>
94
95
96
97namespace internal
98{
102 namespace BlockVectorIterators
103 {
119 template <typename BlockVectorType, bool Constness>
121 {
122 public:
127
134 std::conditional_t<Constness,
135 const typename BlockVectorType::value_type,
136 typename BlockVectorType::value_type>;
137
151 using iterator_category = std::random_access_iterator_tag;
152 using difference_type = std::ptrdiff_t;
153 using reference = typename BlockVectorType::reference;
155
157 std::conditional_t<Constness,
159 typename BlockVectorType::BlockType::reference>;
160
166 std::conditional_t<Constness, const BlockVectorType, BlockVectorType>;
167
177
186
187
192
193 private:
204
205 public:
209 Iterator &
211
218 operator*() const;
219
226
231 Iterator &
233
241
246 Iterator &
248
256
261 template <bool OtherConstness>
262 bool
264
269 template <bool OtherConstness>
270 bool
272
278 template <bool OtherConstness>
279 bool
281
285 template <bool OtherConstness>
286 bool
288
292 template <bool OtherConstness>
293 bool
295
299 template <bool OtherConstness>
300 bool
302
306 template <bool OtherConstness>
309
315 operator+(const difference_type &d) const;
316
322 operator-(const difference_type &d) const;
323
328 Iterator &
330
335 Iterator &
337
348 "Your program tried to compare iterators pointing to "
349 "different block vectors. There is no reasonable way "
350 "to do this.");
351
353 private:
360
365
370 unsigned int current_block;
372
381
385 void
387
391 void
393
394 // Mark all other instances of this template as friends.
395 template <typename, bool>
396 friend class Iterator;
397 };
398 } // namespace BlockVectorIterators
399} // namespace internal
400
401
439template <typename VectorType>
440class BlockVectorBase : public ReadVector<typename VectorType::value_type>
441{
442public:
446 using BlockType = VectorType;
447
448 /*
449 * Declare standard types used in
450 * all containers. These types
451 * parallel those in the
452 * <tt>C++</tt> standard
453 * libraries
454 * <tt>std::vector<...></tt>
455 * class. This includes iterator
456 * types.
457 */
458 using value_type = typename BlockType::value_type;
460 using const_pointer = const value_type *;
461 using iterator =
465 using reference = typename BlockType::reference;
466 using const_reference = typename BlockType::const_reference;
468
478 using real_type = typename BlockType::real_type;
479
483 BlockVectorBase() = default;
484
488 BlockVectorBase(const BlockVectorBase & /*V*/) = default;
489
495 BlockVectorBase(BlockVectorBase && /*V*/) noexcept = default;
496
503 void
505
516 void
517 compress(VectorOperation::values operation);
518
522 BlockType &
523 block(const unsigned int i);
524
528 const BlockType &
529 block(const unsigned int i) const;
530
536 const BlockIndices &
538
542 unsigned int
543 n_blocks() const;
544
549 virtual size_type
550 size() const override;
551
557 std::size_t
559
578
584
590 begin() const;
591
597
603 end() const;
604
609 operator()(const size_type i) const;
610
615 operator()(const size_type i);
616
623 operator[](const size_type i) const;
624
631 operator[](const size_type i);
632
648 template <typename OtherNumber>
649 void
650 extract_subvector_to(const std::vector<size_type> &indices,
651 std::vector<OtherNumber> &values) const;
652
653 virtual void
654 extract_subvector_to(const ArrayView<const types::global_dof_index> &indices,
655 const ArrayView<value_type> &entries) const override;
656
684 template <typename ForwardIterator, typename OutputIterator>
685 void
686 extract_subvector_to(ForwardIterator indices_begin,
687 const ForwardIterator indices_end,
688 OutputIterator values_begin) const;
689
695 operator=(const value_type s);
696
701 operator=(const BlockVectorBase &V);
702
709 operator=(BlockVectorBase && /*V*/) = default; // NOLINT
710
714 template <typename VectorType2>
716 operator=(const BlockVectorBase<VectorType2> &V);
717
722 operator=(const VectorType &v);
723
728 template <typename VectorType2>
729 bool
730 operator==(const BlockVectorBase<VectorType2> &v) const;
731
736 operator*(const BlockVectorBase &V) const;
737
742 norm_sqr() const;
743
748 mean_value() const;
749
754 l1_norm() const;
755
761 l2_norm() const;
762
768 linfty_norm() const;
769
794 const BlockVectorBase &V,
795 const BlockVectorBase &W);
796
801 bool
802 in_local_range(const size_type global_index) const;
803
809 bool
810 all_zero() const;
811
817 bool
819
824 operator+=(const BlockVectorBase &V);
825
830 operator-=(const BlockVectorBase &V);
831
832
837 template <typename Number>
838 void
839 add(const std::vector<size_type> &indices, const std::vector<Number> &values);
840
845 template <typename Number>
846 void
847 add(const std::vector<size_type> &indices, const Vector<Number> &values);
848
854 template <typename Number>
855 void
856 add(const size_type n_elements,
857 const size_type *indices,
858 const Number *values);
859
864 void
865 add(const value_type s);
866
870 void
871 add(const value_type a, const BlockVectorBase &V);
872
876 void
878 const BlockVectorBase &V,
879 const value_type b,
880 const BlockVectorBase &W);
881
885 void
886 sadd(const value_type s, const BlockVectorBase &V);
887
891 void
892 sadd(const value_type s, const value_type a, const BlockVectorBase &V);
893
897 void
899 const value_type a,
900 const BlockVectorBase &V,
901 const value_type b,
902 const BlockVectorBase &W);
903
907 void
909 const value_type a,
910 const BlockVectorBase &V,
911 const value_type b,
912 const BlockVectorBase &W,
913 const value_type c,
914 const BlockVectorBase &X);
915
920 operator*=(const value_type factor);
921
926 operator/=(const value_type factor);
927
932 template <class BlockVector2>
933 void
934 scale(const BlockVector2 &v);
935
939 template <class BlockVector2>
940 void
941 equ(const value_type a, const BlockVector2 &V);
942
947 void
949
957
962 std::size_t
964
965protected:
969 std::vector<VectorType> components;
970
976
977 // Make the iterator class a friend.
978 template <typename N, bool C>
979 friend class ::internal::BlockVectorIterators::Iterator;
980
981 template <typename>
982 friend class BlockVectorBase;
983};
984
985
988/*----------------------- Inline functions ----------------------------------*/
989
990
991#ifndef DOXYGEN
992namespace internal
993{
994 namespace BlockVectorIterators
995 {
996 template <typename BlockVectorType, bool Constness>
997 inline Iterator<BlockVectorType, Constness>::Iterator(
998 const Iterator<BlockVectorType, Constness> &c)
999 : parent(c.parent)
1000 , global_index(c.global_index)
1001 , current_block(c.current_block)
1002 , index_within_block(c.index_within_block)
1003 , next_break_forward(c.next_break_forward)
1004 , next_break_backward(c.next_break_backward)
1005 {}
1006
1007
1008
1009 template <typename BlockVectorType, bool Constness>
1010 inline Iterator<BlockVectorType, Constness>::Iterator(
1011 const Iterator<BlockVectorType, !Constness> &c)
1012 : parent(c.parent)
1013 , global_index(c.global_index)
1014 , current_block(c.current_block)
1015 , index_within_block(c.index_within_block)
1016 , next_break_forward(c.next_break_forward)
1017 , next_break_backward(c.next_break_backward)
1018 {
1019 // Only permit copy-constructing const iterators from non-const
1020 // iterators, and not vice versa (i.e., Constness must always be
1021 // true).
1022 static_assert(Constness == true,
1023 "Constructing a non-const iterator from a const iterator "
1024 "does not make sense.");
1025 }
1026
1027
1028
1029 template <typename BlockVectorType, bool Constness>
1030 inline Iterator<BlockVectorType, Constness>::Iterator(
1031 BlockVector &parent,
1032 const size_type global_index,
1033 const size_type current_block,
1034 const size_type index_within_block,
1035 const size_type next_break_forward,
1036 const size_type next_break_backward)
1037 : parent(&parent)
1038 , global_index(global_index)
1039 , current_block(current_block)
1040 , index_within_block(index_within_block)
1041 , next_break_forward(next_break_forward)
1042 , next_break_backward(next_break_backward)
1043 {}
1044
1045
1046
1047 template <typename BlockVectorType, bool Constness>
1048 inline Iterator<BlockVectorType, Constness> &
1049 Iterator<BlockVectorType, Constness>::operator=(const Iterator &c) =
1050 default;
1051
1052
1053
1054 template <typename BlockVectorType, bool Constness>
1055 inline typename Iterator<BlockVectorType, Constness>::dereference_type
1056 Iterator<BlockVectorType, Constness>::operator*() const
1057 {
1058 return parent->block(current_block)(index_within_block);
1059 }
1060
1061
1062
1063 template <typename BlockVectorType, bool Constness>
1064 inline typename Iterator<BlockVectorType, Constness>::dereference_type
1065 Iterator<BlockVectorType, Constness>::operator[](
1066 const difference_type d) const
1067 {
1068 // if the index pointed to is
1069 // still within the block we
1070 // currently point into, then we
1071 // can save the computation of
1072 // the block
1073 if ((global_index + d >= next_break_backward) &&
1074 (global_index + d <= next_break_forward))
1075 return parent->block(current_block)(index_within_block + d);
1076
1077 // if the index is not within the
1078 // block of the block vector into
1079 // which we presently point, then
1080 // there is no way: we have to
1081 // search for the block. this can
1082 // be done through the parent
1083 // class as well.
1084 return (*parent)(global_index + d);
1085 }
1086
1087
1088
1089 template <typename BlockVectorType, bool Constness>
1090 inline Iterator<BlockVectorType, Constness> &
1091 Iterator<BlockVectorType, Constness>::operator++()
1092 {
1093 move_forward();
1094 return *this;
1095 }
1096
1097
1098
1099 template <typename BlockVectorType, bool Constness>
1100 inline Iterator<BlockVectorType, Constness>
1101 Iterator<BlockVectorType, Constness>::operator++(int)
1102 {
1103 const Iterator old_value = *this;
1104 move_forward();
1105 return old_value;
1106 }
1107
1108
1109
1110 template <typename BlockVectorType, bool Constness>
1111 inline Iterator<BlockVectorType, Constness> &
1112 Iterator<BlockVectorType, Constness>::operator--()
1113 {
1114 move_backward();
1115 return *this;
1116 }
1117
1118
1119
1120 template <typename BlockVectorType, bool Constness>
1121 inline Iterator<BlockVectorType, Constness>
1122 Iterator<BlockVectorType, Constness>::operator--(int)
1123 {
1124 const Iterator old_value = *this;
1125 move_backward();
1126 return old_value;
1127 }
1128
1129
1130
1131 template <typename BlockVectorType, bool Constness>
1132 template <bool OtherConstness>
1133 inline bool
1134 Iterator<BlockVectorType, Constness>::operator==(
1135 const Iterator<BlockVectorType, OtherConstness> &i) const
1136 {
1137 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1138
1139 return (global_index == i.global_index);
1140 }
1141
1142
1143
1144 template <typename BlockVectorType, bool Constness>
1145 template <bool OtherConstness>
1146 inline bool
1147 Iterator<BlockVectorType, Constness>::operator!=(
1148 const Iterator<BlockVectorType, OtherConstness> &i) const
1149 {
1150 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1151
1152 return (global_index != i.global_index);
1153 }
1154
1155
1156
1157 template <typename BlockVectorType, bool Constness>
1158 template <bool OtherConstness>
1159 inline bool
1160 Iterator<BlockVectorType, Constness>::operator<(
1161 const Iterator<BlockVectorType, OtherConstness> &i) const
1162 {
1163 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1164
1165 return (global_index < i.global_index);
1166 }
1167
1168
1169
1170 template <typename BlockVectorType, bool Constness>
1171 template <bool OtherConstness>
1172 inline bool
1173 Iterator<BlockVectorType, Constness>::operator<=(
1174 const Iterator<BlockVectorType, OtherConstness> &i) const
1175 {
1176 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1177
1178 return (global_index <= i.global_index);
1179 }
1180
1181
1182
1183 template <typename BlockVectorType, bool Constness>
1184 template <bool OtherConstness>
1185 inline bool
1186 Iterator<BlockVectorType, Constness>::operator>(
1187 const Iterator<BlockVectorType, OtherConstness> &i) const
1188 {
1189 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1190
1191 return (global_index > i.global_index);
1192 }
1193
1194
1195
1196 template <typename BlockVectorType, bool Constness>
1197 template <bool OtherConstness>
1198 inline bool
1199 Iterator<BlockVectorType, Constness>::operator>=(
1200 const Iterator<BlockVectorType, OtherConstness> &i) const
1201 {
1202 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1203
1204 return (global_index >= i.global_index);
1205 }
1206
1207
1208
1209 template <typename BlockVectorType, bool Constness>
1210 template <bool OtherConstness>
1211 inline typename Iterator<BlockVectorType, Constness>::difference_type
1212 Iterator<BlockVectorType, Constness>::operator-(
1213 const Iterator<BlockVectorType, OtherConstness> &i) const
1214 {
1215 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1216
1217 return (static_cast<signed int>(global_index) -
1218 static_cast<signed int>(i.global_index));
1219 }
1220
1221
1222
1223 template <typename BlockVectorType, bool Constness>
1224 inline Iterator<BlockVectorType, Constness>
1225 Iterator<BlockVectorType, Constness>::operator+(
1226 const difference_type &d) const
1227 {
1228 // if the index pointed to is
1229 // still within the block we
1230 // currently point into, then we
1231 // can save the computation of
1232 // the block
1233 if ((global_index + d >= next_break_backward) &&
1234 (global_index + d <= next_break_forward))
1235 return Iterator(*parent,
1236 global_index + d,
1237 current_block,
1238 index_within_block + d,
1239 next_break_forward,
1240 next_break_backward);
1241 else
1242 // outside present block, so
1243 // have to seek new block
1244 // anyway
1245 return Iterator(*parent, global_index + d);
1246 }
1247
1248
1249
1250 template <typename BlockVectorType, bool Constness>
1251 inline Iterator<BlockVectorType, Constness>
1252 Iterator<BlockVectorType, Constness>::operator-(
1253 const difference_type &d) const
1254 {
1255 // if the index pointed to is
1256 // still within the block we
1257 // currently point into, then we
1258 // can save the computation of
1259 // the block
1260 if ((global_index - d >= next_break_backward) &&
1261 (global_index - d <= next_break_forward))
1262 return Iterator(*parent,
1263 global_index - d,
1264 current_block,
1265 index_within_block - d,
1266 next_break_forward,
1267 next_break_backward);
1268 else
1269 // outside present block, so
1270 // have to seek new block
1271 // anyway
1272 return Iterator(*parent, global_index - d);
1273 }
1274
1275
1276
1277 template <typename BlockVectorType, bool Constness>
1278 inline Iterator<BlockVectorType, Constness> &
1279 Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1280 {
1281 // if the index pointed to is
1282 // still within the block we
1283 // currently point into, then we
1284 // can save the computation of
1285 // the block
1286 if ((global_index + d >= next_break_backward) &&
1287 (global_index + d <= next_break_forward))
1288 {
1289 global_index += d;
1290 index_within_block += d;
1291 }
1292 else
1293 // outside present block, so
1294 // have to seek new block
1295 // anyway
1296 *this = Iterator(*parent, global_index + d);
1297
1298 return *this;
1299 }
1300
1301
1302
1303 template <typename BlockVectorType, bool Constness>
1304 inline Iterator<BlockVectorType, Constness> &
1305 Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1306 {
1307 // if the index pointed to is
1308 // still within the block we
1309 // currently point into, then we
1310 // can save the computation of
1311 // the block
1312 if ((global_index - d >= next_break_backward) &&
1313 (global_index - d <= next_break_forward))
1314 {
1315 global_index -= d;
1316 index_within_block -= d;
1317 }
1318 else
1319 // outside present block, so
1320 // have to seek new block
1321 // anyway
1322 *this = Iterator(*parent, global_index - d);
1323
1324 return *this;
1325 }
1326
1327
1328 template <typename BlockVectorType, bool Constness>
1329 Iterator<BlockVectorType, Constness>::Iterator(BlockVector &parent,
1330 const size_type global_index)
1331 : parent(&parent)
1332 , global_index(global_index)
1333 {
1334 // find which block we are
1335 // in. for this, take into
1336 // account that it happens at
1337 // times that people want to
1338 // initialize iterators
1339 // past-the-end
1340 if (global_index < parent.size())
1341 {
1342 const std::pair<size_type, size_type> indices =
1343 parent.block_indices.global_to_local(global_index);
1344 current_block = indices.first;
1345 index_within_block = indices.second;
1346
1347 next_break_backward =
1348 parent.block_indices.local_to_global(current_block, 0);
1349 next_break_forward =
1350 (parent.block_indices.local_to_global(current_block, 0) +
1351 parent.block_indices.block_size(current_block) - 1);
1352 }
1353 else
1354 // past the end. only have one
1355 // value for this
1356 {
1357 this->global_index = parent.size();
1358 current_block = parent.n_blocks();
1359 index_within_block = 0;
1360 next_break_backward = global_index;
1361 next_break_forward = numbers::invalid_size_type;
1362 };
1363 }
1364
1365
1366
1367 template <typename BlockVectorType, bool Constness>
1368 void
1369 Iterator<BlockVectorType, Constness>::move_forward()
1370 {
1371 if (global_index != next_break_forward)
1372 ++index_within_block;
1373 else
1374 {
1375 // ok, we traverse a boundary
1376 // between blocks:
1377 index_within_block = 0;
1378 ++current_block;
1379
1380 // break backwards is now old
1381 // break forward
1382 next_break_backward = next_break_forward + 1;
1383
1384 // compute new break forward
1385 if (current_block < parent->block_indices.size())
1386 next_break_forward +=
1387 parent->block_indices.block_size(current_block);
1388 else
1389 // if we are beyond the end,
1390 // then move the next
1391 // boundary arbitrarily far
1392 // away
1393 next_break_forward = numbers::invalid_size_type;
1394 };
1395
1396 ++global_index;
1397 }
1398
1399
1400
1401 template <typename BlockVectorType, bool Constness>
1402 void
1403 Iterator<BlockVectorType, Constness>::move_backward()
1404 {
1405 if (global_index != next_break_backward)
1406 --index_within_block;
1407 else if (current_block != 0)
1408 {
1409 // ok, we traverse a boundary
1410 // between blocks:
1411 --current_block;
1412 index_within_block =
1413 parent->block_indices.block_size(current_block) - 1;
1414
1415 // break forwards is now old
1416 // break backward
1417 next_break_forward = next_break_backward - 1;
1418
1419 // compute new break forward
1420 next_break_backward -=
1421 parent->block_indices.block_size(current_block);
1422 }
1423 else
1424 // current block was 0, we now
1425 // get into unspecified terrain
1426 {
1427 --current_block;
1428 index_within_block = numbers::invalid_size_type;
1429 next_break_forward = 0;
1430 next_break_backward = 0;
1431 };
1432
1433 --global_index;
1434 }
1435
1436
1437 } // namespace BlockVectorIterators
1438
1439} // namespace internal
1440
1441
1442
1443template <typename VectorType>
1446{
1447 return block_indices.total_size();
1448}
1449
1450
1451
1452template <typename VectorType>
1453inline std::size_t
1455{
1456 std::size_t local_size = 0;
1457 for (unsigned int b = 0; b < n_blocks(); ++b)
1458 local_size += block(b).locally_owned_size();
1459 return local_size;
1460}
1461
1462
1463
1464template <typename VectorType>
1465inline IndexSet
1467{
1468 IndexSet is(size());
1469
1470 // copy index sets from blocks into the global one, shifted
1471 // by the appropriate amount for each block
1472 for (unsigned int b = 0; b < n_blocks(); ++b)
1473 {
1474 IndexSet x = block(b).locally_owned_elements();
1476 }
1477
1478 is.compress();
1479
1480 return is;
1481}
1482
1483
1484
1485template <typename VectorType>
1486inline unsigned int
1488{
1489 return block_indices.size();
1490}
1491
1492
1493template <typename VectorType>
1495BlockVectorBase<VectorType>::block(const unsigned int i)
1496{
1498
1499 return components[i];
1500}
1501
1502
1503
1504template <typename VectorType>
1505inline const typename BlockVectorBase<VectorType>::BlockType &
1506BlockVectorBase<VectorType>::block(const unsigned int i) const
1507{
1509
1510 return components[i];
1511}
1512
1513
1514
1515template <typename VectorType>
1516inline const BlockIndices &
1518{
1519 return block_indices;
1520}
1521
1522
1523template <typename VectorType>
1524inline void
1526{
1527 std::vector<size_type> sizes(n_blocks());
1528
1529 for (size_type i = 0; i < n_blocks(); ++i)
1530 sizes[i] = block(i).size();
1531
1532 block_indices.reinit(sizes);
1533}
1534
1535
1536
1537template <typename VectorType>
1538inline void
1540{
1541 for (unsigned int i = 0; i < n_blocks(); ++i)
1542 block(i).compress(operation);
1543}
1544
1545
1546
1547template <typename VectorType>
1550{
1551 return iterator(*this, 0U);
1552}
1553
1554
1555
1556template <typename VectorType>
1559{
1560 return const_iterator(*this, 0U);
1561}
1562
1563
1564template <typename VectorType>
1567{
1568 return iterator(*this, size());
1569}
1570
1571
1572
1573template <typename VectorType>
1576{
1577 return const_iterator(*this, size());
1578}
1579
1580
1581template <typename VectorType>
1582inline bool
1584{
1585 const std::pair<size_type, size_type> local_index =
1586 block_indices.global_to_local(global_index);
1587
1588 return components[local_index.first].in_local_range(global_index);
1589}
1590
1591
1592template <typename VectorType>
1593bool
1595{
1596 for (size_type i = 0; i < n_blocks(); ++i)
1597 if (components[i].all_zero() == false)
1598 return false;
1599
1600 return true;
1601}
1602
1603
1604
1605template <typename VectorType>
1606bool
1608{
1609 for (size_type i = 0; i < n_blocks(); ++i)
1610 if (components[i].is_non_negative() == false)
1611 return false;
1612
1613 return true;
1614}
1615
1616
1617
1618template <typename VectorType>
1621 const BlockVectorBase<VectorType> &v) const
1622{
1623 Assert(n_blocks() == v.n_blocks(),
1625
1626 value_type sum = 0.;
1627 for (size_type i = 0; i < n_blocks(); ++i)
1628 sum += components[i] * v.components[i];
1629
1630 return sum;
1631}
1632
1633
1634template <typename VectorType>
1637{
1638 real_type sum = 0.;
1639 for (size_type i = 0; i < n_blocks(); ++i)
1640 sum += components[i].norm_sqr();
1641
1642 return sum;
1643}
1644
1645
1646
1647template <typename VectorType>
1650{
1651 value_type sum = 0.;
1652 // need to do static_cast as otherwise it won't work with
1653 // value_type=complex<T>
1654 for (size_type i = 0; i < n_blocks(); ++i)
1655 sum += components[i].mean_value() *
1657 components[i].size()));
1658
1660}
1661
1662
1663
1664template <typename VectorType>
1667{
1668 real_type sum = 0.;
1669 for (size_type i = 0; i < n_blocks(); ++i)
1670 sum += components[i].l1_norm();
1671
1672 return sum;
1673}
1674
1675
1676
1677template <typename VectorType>
1680{
1681 return std::sqrt(norm_sqr());
1682}
1683
1684
1685
1686template <typename VectorType>
1689{
1690 real_type sum = 0.;
1691 for (size_type i = 0; i < n_blocks(); ++i)
1692 {
1693 value_type newval = components[i].linfty_norm();
1694 if (sum < newval)
1695 sum = newval;
1696 }
1697 return sum;
1698}
1699
1700
1701
1702template <typename VectorType>
1708{
1709 AssertDimension(n_blocks(), V.n_blocks());
1711
1712 value_type sum = 0.;
1713 for (size_type i = 0; i < n_blocks(); ++i)
1714 sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1715
1716 return sum;
1717}
1718
1719
1720
1721template <typename VectorType>
1724{
1725 Assert(n_blocks() == v.n_blocks(),
1727
1728 for (size_type i = 0; i < n_blocks(); ++i)
1729 {
1730 components[i] += v.components[i];
1731 }
1732
1733 return *this;
1734}
1735
1736
1737
1738template <typename VectorType>
1741{
1742 Assert(n_blocks() == v.n_blocks(),
1744
1745 for (size_type i = 0; i < n_blocks(); ++i)
1746 {
1747 components[i] -= v.components[i];
1748 }
1749 return *this;
1750}
1751
1752
1753
1754template <typename VectorType>
1755template <typename Number>
1756inline void
1757BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1758 const std::vector<Number> &values)
1759{
1760 Assert(indices.size() == values.size(),
1761 ExcDimensionMismatch(indices.size(), values.size()));
1762 add(indices.size(), indices.data(), values.data());
1763}
1764
1765
1766
1767template <typename VectorType>
1768template <typename Number>
1769inline void
1770BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1771 const Vector<Number> &values)
1772{
1773 Assert(indices.size() == values.size(),
1774 ExcDimensionMismatch(indices.size(), values.size()));
1775 const size_type n_indices = indices.size();
1776 for (size_type i = 0; i < n_indices; ++i)
1777 (*this)(indices[i]) += values(i);
1778}
1779
1780
1781
1782template <typename VectorType>
1783template <typename Number>
1784inline void
1786 const size_type *indices,
1787 const Number *values)
1788{
1789 for (size_type i = 0; i < n_indices; ++i)
1790 (*this)(indices[i]) += values[i];
1791}
1792
1793
1794
1795template <typename VectorType>
1796void
1798{
1799 AssertIsFinite(a);
1800
1801 for (size_type i = 0; i < n_blocks(); ++i)
1802 {
1803 components[i].add(a);
1804 }
1805}
1806
1807
1808
1809template <typename VectorType>
1810void
1813{
1814 AssertIsFinite(a);
1815
1816 Assert(n_blocks() == v.n_blocks(),
1818
1819 for (size_type i = 0; i < n_blocks(); ++i)
1820 {
1821 components[i].add(a, v.components[i]);
1822 }
1823}
1824
1825
1826
1827template <typename VectorType>
1828void
1831 const value_type b,
1833{
1834 AssertIsFinite(a);
1835 AssertIsFinite(b);
1836
1837 Assert(n_blocks() == v.n_blocks(),
1839 Assert(n_blocks() == w.n_blocks(),
1840 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1841
1842
1843 for (size_type i = 0; i < n_blocks(); ++i)
1844 {
1845 components[i].add(a, v.components[i], b, w.components[i]);
1846 }
1847}
1848
1849
1850
1851template <typename VectorType>
1852void
1855{
1856 AssertIsFinite(x);
1857
1858 Assert(n_blocks() == v.n_blocks(),
1860
1861 for (size_type i = 0; i < n_blocks(); ++i)
1862 {
1863 components[i].sadd(x, v.components[i]);
1864 }
1865}
1866
1867
1868
1869template <typename VectorType>
1870void
1872 const value_type a,
1874{
1875 AssertIsFinite(x);
1876 AssertIsFinite(a);
1877
1878 Assert(n_blocks() == v.n_blocks(),
1880
1881 for (size_type i = 0; i < n_blocks(); ++i)
1882 {
1883 components[i].sadd(x, a, v.components[i]);
1884 }
1885}
1886
1887
1888
1889template <typename VectorType>
1890void
1892 const value_type a,
1894 const value_type b,
1896{
1897 AssertIsFinite(x);
1898 AssertIsFinite(a);
1899 AssertIsFinite(b);
1900
1901 Assert(n_blocks() == v.n_blocks(),
1903 Assert(n_blocks() == w.n_blocks(),
1904 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1905
1906 for (size_type i = 0; i < n_blocks(); ++i)
1907 {
1908 components[i].sadd(x, a, v.components[i], b, w.components[i]);
1909 }
1910}
1911
1912
1913
1914template <typename VectorType>
1915void
1917 const value_type a,
1919 const value_type b,
1921 const value_type c,
1923{
1924 AssertIsFinite(x);
1925 AssertIsFinite(a);
1926 AssertIsFinite(b);
1927 AssertIsFinite(c);
1928
1929 Assert(n_blocks() == v.n_blocks(),
1931 Assert(n_blocks() == w.n_blocks(),
1932 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1933 Assert(n_blocks() == y.n_blocks(),
1935
1936 for (size_type i = 0; i < n_blocks(); ++i)
1937 {
1938 components[i].sadd(
1939 x, a, v.components[i], b, w.components[i], c, y.components[i]);
1940 }
1941}
1942
1943
1944
1945template <typename VectorType>
1946template <class BlockVector2>
1947void
1948BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
1949{
1950 Assert(n_blocks() == v.n_blocks(),
1951 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1952 for (size_type i = 0; i < n_blocks(); ++i)
1953 components[i].scale(v.block(i));
1954}
1955
1956
1957
1958template <typename VectorType>
1959std::size_t
1961{
1962 return (MemoryConsumption::memory_consumption(this->block_indices) +
1963 MemoryConsumption::memory_consumption(this->components));
1964}
1965
1966
1967
1968template <typename VectorType>
1969template <class BlockVector2>
1970void
1971BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
1972{
1973 AssertIsFinite(a);
1974
1975 Assert(n_blocks() == v.n_blocks(),
1976 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1977
1978 for (size_type i = 0; i < n_blocks(); ++i)
1979 components[i].equ(a, v.components[i]);
1980}
1981
1982
1983
1984template <typename VectorType>
1985void
1987{
1988 for (size_type i = 0; i < n_blocks(); ++i)
1989 block(i).update_ghost_values();
1990}
1991
1992
1993
1994template <typename VectorType>
1997{
1998 if (n_blocks() > 0)
1999 return block(0).get_mpi_communicator();
2000 else
2001 return MPI_COMM_SELF;
2002}
2003
2004
2005
2006template <typename VectorType>
2009{
2010 AssertIsFinite(s);
2011
2012 for (size_type i = 0; i < n_blocks(); ++i)
2013 components[i] = s;
2014
2015 return *this;
2016}
2017
2018
2019template <typename VectorType>
2022{
2024
2025 for (size_type i = 0; i < n_blocks(); ++i)
2026 components[i] = v.components[i];
2027
2028 return *this;
2029}
2030
2031
2032template <typename VectorType>
2033template <typename VectorType2>
2036{
2038
2039 for (size_type i = 0; i < n_blocks(); ++i)
2040 components[i] = v.components[i];
2041
2042 return *this;
2043}
2044
2045
2046
2047template <typename VectorType>
2049BlockVectorBase<VectorType>::operator=(const VectorType &v)
2050{
2051 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2052
2053 size_type index_v = 0;
2054 for (size_type b = 0; b < n_blocks(); ++b)
2055 for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2056 block(b)(i) = v(index_v);
2057
2058 return *this;
2059}
2060
2061
2062
2063template <typename VectorType>
2064template <typename VectorType2>
2065inline bool
2067 const BlockVectorBase<VectorType2> &v) const
2068{
2070
2071 for (size_type i = 0; i < n_blocks(); ++i)
2072 if (!(components[i] == v.components[i]))
2073 return false;
2074
2075 return true;
2076}
2077
2078
2079
2080template <typename VectorType>
2083{
2084 AssertIsFinite(factor);
2085
2086 for (size_type i = 0; i < n_blocks(); ++i)
2087 components[i] *= factor;
2088
2089 return *this;
2090}
2091
2092
2093
2094template <typename VectorType>
2097{
2098 AssertIsFinite(factor);
2099 Assert(factor != 0., ExcDivideByZero());
2100
2101 const value_type inverse_factor = value_type(1.) / factor;
2102
2103 for (size_type i = 0; i < n_blocks(); ++i)
2104 components[i] *= inverse_factor;
2105
2106 return *this;
2107}
2108
2109
2110template <typename VectorType>
2113{
2114 const std::pair<unsigned int, size_type> local_index =
2116 return components[local_index.first](local_index.second);
2117}
2118
2119
2120
2121template <typename VectorType>
2124{
2125 const std::pair<unsigned int, size_type> local_index =
2127 return components[local_index.first](local_index.second);
2128}
2129
2130
2131
2132template <typename VectorType>
2135{
2136 return operator()(i);
2137}
2138
2139
2140
2141template <typename VectorType>
2144{
2145 return operator()(i);
2146}
2147
2148
2149
2150template <typename VectorType>
2151template <typename OtherNumber>
2152inline void
2154 const std::vector<size_type> &indices,
2155 std::vector<OtherNumber> &values) const
2156{
2157 for (size_type i = 0; i < indices.size(); ++i)
2158 values[i] = operator()(indices[i]);
2159}
2160
2161
2162
2163template <typename VectorType>
2164inline void
2167 const ArrayView<value_type> &entries) const
2168{
2169 AssertDimension(indices.size(), entries.size());
2170 for (unsigned int i = 0; i < indices.size(); ++i)
2171 {
2172 AssertIndexRange(indices[i], size());
2173 entries[i] = (*this)[indices[i]];
2174 }
2175}
2176
2177
2178
2179template <typename VectorType>
2180template <typename ForwardIterator, typename OutputIterator>
2181inline void
2183 ForwardIterator indices_begin,
2184 const ForwardIterator indices_end,
2185 OutputIterator values_begin) const
2186{
2187 while (indices_begin != indices_end)
2188 {
2189 *values_begin = operator()(*indices_begin);
2190 ++indices_begin;
2191 ++values_begin;
2192 }
2193}
2194
2195#endif // DOXYGEN
2196
2198
2199#endif
std::size_t size() const
Definition array_view.h:689
size_type block_size(const unsigned int i) const
size_type total_size() const
void reinit(const unsigned int n_blocks, const size_type n_elements_per_block)
unsigned int size() const
size_type local_to_global(const unsigned int block, const size_type index) const
size_type block_start(const unsigned int i) const
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
BlockVectorBase & operator+=(const BlockVectorBase &V)
BlockVectorBase()=default
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< OtherNumber > &values) const
value_type mean_value() const
virtual size_type size() const override
void add(const std::vector< size_type > &indices, const std::vector< Number > &values)
::internal::BlockVectorIterators::Iterator< BlockVectorBase, false > iterator
BlockVectorBase & operator/=(const value_type factor)
unsigned int n_blocks() const
real_type norm_sqr() const
const value_type * const_pointer
typename BlockType::const_reference const_reference
void update_ghost_values() const
std::size_t memory_consumption() const
real_type l1_norm() const
types::global_dof_index size_type
real_type linfty_norm() const
value_type add_and_dot(const value_type a, const BlockVectorBase &V, const BlockVectorBase &W)
typename BlockType::value_type value_type
value_type operator*(const BlockVectorBase &V) const
void compress(VectorOperation::values operation)
bool all_zero() const
BlockVectorBase & operator=(const value_type s)
void collect_sizes()
BlockVectorBase & operator*=(const value_type factor)
void sadd(const value_type s, const BlockVectorBase &V)
std::vector< VectorType > components
iterator begin()
bool in_local_range(const size_type global_index) const
::internal::BlockVectorIterators::Iterator< BlockVectorBase, true > const_iterator
MPI_Comm get_mpi_communicator() const
BlockVectorBase(BlockVectorBase &&) noexcept=default
BlockVectorBase(const BlockVectorBase &)=default
typename BlockType::reference reference
iterator end()
value_type operator[](const size_type i) const
BlockIndices block_indices
typename BlockType::real_type real_type
void scale(const BlockVector2 &v)
BlockType & block(const unsigned int i)
real_type l2_norm() const
value_type operator()(const size_type i) const
IndexSet locally_owned_elements() const
bool is_non_negative() const
BlockVectorBase & operator-=(const BlockVectorBase &V)
friend class ::internal::BlockVectorIterators::Iterator
bool operator==(const BlockVectorBase< VectorType2 > &v) const
std::size_t locally_owned_size() const
const BlockIndices & get_block_indices() const
void equ(const value_type a, const BlockVector2 &V)
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1842
bool operator==(const Iterator< BlockVectorType, OtherConstness > &i) const
dereference_type operator[](const difference_type d) const
Iterator & operator=(const Iterator &c)
bool operator<(const Iterator< BlockVectorType, OtherConstness > &i) const
Iterator(BlockVector &parent, const size_type global_index, const size_type current_block, const size_type index_within_block, const size_type next_break_forward, const size_type next_break_backward)
std::conditional_t< Constness, const typename BlockVectorType::value_type, typename BlockVectorType::value_type > value_type
Iterator(const Iterator< BlockVectorType, !Constness > &c)
std::conditional_t< Constness, const BlockVectorType, BlockVectorType > BlockVector
dereference_type operator*() const
Iterator operator-(const difference_type &d) const
difference_type operator-(const Iterator< BlockVectorType, OtherConstness > &i) const
std::random_access_iterator_tag iterator_category
Iterator(BlockVector &parent, const size_type global_index)
bool operator>=(const Iterator< BlockVectorType, OtherConstness > &i) const
Iterator & operator-=(const difference_type &d)
std::conditional_t< Constness, value_type, typename BlockVectorType::BlockType::reference > dereference_type
bool operator<=(const Iterator< BlockVectorType, OtherConstness > &i) const
typename BlockVectorType::reference reference
Iterator operator+(const difference_type &d) const
Iterator & operator+=(const difference_type &d)
bool operator!=(const Iterator< BlockVectorType, OtherConstness > &i) const
bool operator>(const Iterator< BlockVectorType, OtherConstness > &i) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define Assert(cond, exc)
static ::ExceptionBase & ExcPointerToDifferentVectors()
#define AssertIsFinite(number)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcDifferentBlockIndices()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static const bool value
constexpr char V
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
T sum(const T &t, const MPI_Comm mpi_communicator)
constexpr bool has_n_blocks
decltype(std::declval< const T >().n_blocks()) has_n_blocks_t
decltype(std::declval< const T >().block(0)) has_block_t
constexpr bool is_block_vector
constexpr bool has_block
constexpr types::global_dof_index invalid_size_type
Definition types.h:250
STL namespace.
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition types.h:32
unsigned int global_dof_index
Definition types.h:94