deal.II version GIT relicensing-2017-g7ae82fad8b 2024-10-22 19:20: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_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 {
1051 parent = c.parent;
1052 global_index = c.global_index;
1053 index_within_block = c.index_within_block;
1054 current_block = c.current_block;
1055 next_break_forward = c.next_break_forward;
1056 next_break_backward = c.next_break_backward;
1057
1058 return *this;
1059 }
1060
1061
1062
1063 template <typename BlockVectorType, bool Constness>
1064 inline typename Iterator<BlockVectorType, Constness>::dereference_type
1065 Iterator<BlockVectorType, Constness>::operator*() const
1066 {
1067 return parent->block(current_block)(index_within_block);
1068 }
1069
1070
1071
1072 template <typename BlockVectorType, bool Constness>
1073 inline typename Iterator<BlockVectorType, Constness>::dereference_type
1074 Iterator<BlockVectorType, Constness>::operator[](
1075 const difference_type d) const
1076 {
1077 // if the index pointed to is
1078 // still within the block we
1079 // currently point into, then we
1080 // can save the computation of
1081 // the block
1082 if ((global_index + d >= next_break_backward) &&
1083 (global_index + d <= next_break_forward))
1084 return parent->block(current_block)(index_within_block + d);
1085
1086 // if the index is not within the
1087 // block of the block vector into
1088 // which we presently point, then
1089 // there is no way: we have to
1090 // search for the block. this can
1091 // be done through the parent
1092 // class as well.
1093 return (*parent)(global_index + d);
1094 }
1095
1096
1097
1098 template <typename BlockVectorType, bool Constness>
1099 inline Iterator<BlockVectorType, Constness> &
1100 Iterator<BlockVectorType, Constness>::operator++()
1101 {
1102 move_forward();
1103 return *this;
1104 }
1105
1106
1107
1108 template <typename BlockVectorType, bool Constness>
1109 inline Iterator<BlockVectorType, Constness>
1110 Iterator<BlockVectorType, Constness>::operator++(int)
1111 {
1112 const Iterator old_value = *this;
1113 move_forward();
1114 return old_value;
1115 }
1116
1117
1118
1119 template <typename BlockVectorType, bool Constness>
1120 inline Iterator<BlockVectorType, Constness> &
1121 Iterator<BlockVectorType, Constness>::operator--()
1122 {
1123 move_backward();
1124 return *this;
1125 }
1126
1127
1128
1129 template <typename BlockVectorType, bool Constness>
1130 inline Iterator<BlockVectorType, Constness>
1131 Iterator<BlockVectorType, Constness>::operator--(int)
1132 {
1133 const Iterator old_value = *this;
1134 move_backward();
1135 return old_value;
1136 }
1137
1138
1139
1140 template <typename BlockVectorType, bool Constness>
1141 template <bool OtherConstness>
1142 inline bool
1143 Iterator<BlockVectorType, Constness>::operator==(
1144 const Iterator<BlockVectorType, OtherConstness> &i) const
1145 {
1146 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1147
1148 return (global_index == i.global_index);
1149 }
1150
1151
1152
1153 template <typename BlockVectorType, bool Constness>
1154 template <bool OtherConstness>
1155 inline bool
1156 Iterator<BlockVectorType, Constness>::operator!=(
1157 const Iterator<BlockVectorType, OtherConstness> &i) const
1158 {
1159 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1160
1161 return (global_index != i.global_index);
1162 }
1163
1164
1165
1166 template <typename BlockVectorType, bool Constness>
1167 template <bool OtherConstness>
1168 inline bool
1169 Iterator<BlockVectorType, Constness>::operator<(
1170 const Iterator<BlockVectorType, OtherConstness> &i) const
1171 {
1172 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1173
1174 return (global_index < i.global_index);
1175 }
1176
1177
1178
1179 template <typename BlockVectorType, bool Constness>
1180 template <bool OtherConstness>
1181 inline bool
1182 Iterator<BlockVectorType, Constness>::operator<=(
1183 const Iterator<BlockVectorType, OtherConstness> &i) const
1184 {
1185 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1186
1187 return (global_index <= i.global_index);
1188 }
1189
1190
1191
1192 template <typename BlockVectorType, bool Constness>
1193 template <bool OtherConstness>
1194 inline bool
1195 Iterator<BlockVectorType, Constness>::operator>(
1196 const Iterator<BlockVectorType, OtherConstness> &i) const
1197 {
1198 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1199
1200 return (global_index > i.global_index);
1201 }
1202
1203
1204
1205 template <typename BlockVectorType, bool Constness>
1206 template <bool OtherConstness>
1207 inline bool
1208 Iterator<BlockVectorType, Constness>::operator>=(
1209 const Iterator<BlockVectorType, OtherConstness> &i) const
1210 {
1211 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1212
1213 return (global_index >= i.global_index);
1214 }
1215
1216
1217
1218 template <typename BlockVectorType, bool Constness>
1219 template <bool OtherConstness>
1220 inline typename Iterator<BlockVectorType, Constness>::difference_type
1221 Iterator<BlockVectorType, Constness>::operator-(
1222 const Iterator<BlockVectorType, OtherConstness> &i) const
1223 {
1224 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1225
1226 return (static_cast<signed int>(global_index) -
1227 static_cast<signed int>(i.global_index));
1228 }
1229
1230
1231
1232 template <typename BlockVectorType, bool Constness>
1233 inline Iterator<BlockVectorType, Constness>
1234 Iterator<BlockVectorType, Constness>::operator+(
1235 const difference_type &d) const
1236 {
1237 // if the index pointed to is
1238 // still within the block we
1239 // currently point into, then we
1240 // can save the computation of
1241 // the block
1242 if ((global_index + d >= next_break_backward) &&
1243 (global_index + d <= next_break_forward))
1244 return Iterator(*parent,
1245 global_index + d,
1246 current_block,
1247 index_within_block + d,
1248 next_break_forward,
1249 next_break_backward);
1250 else
1251 // outside present block, so
1252 // have to seek new block
1253 // anyway
1254 return Iterator(*parent, global_index + d);
1255 }
1256
1257
1258
1259 template <typename BlockVectorType, bool Constness>
1260 inline Iterator<BlockVectorType, Constness>
1261 Iterator<BlockVectorType, Constness>::operator-(
1262 const difference_type &d) const
1263 {
1264 // if the index pointed to is
1265 // still within the block we
1266 // currently point into, then we
1267 // can save the computation of
1268 // the block
1269 if ((global_index - d >= next_break_backward) &&
1270 (global_index - d <= next_break_forward))
1271 return Iterator(*parent,
1272 global_index - d,
1273 current_block,
1274 index_within_block - d,
1275 next_break_forward,
1276 next_break_backward);
1277 else
1278 // outside present block, so
1279 // have to seek new block
1280 // anyway
1281 return Iterator(*parent, global_index - d);
1282 }
1283
1284
1285
1286 template <typename BlockVectorType, bool Constness>
1287 inline Iterator<BlockVectorType, Constness> &
1288 Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1289 {
1290 // if the index pointed to is
1291 // still within the block we
1292 // currently point into, then we
1293 // can save the computation of
1294 // the block
1295 if ((global_index + d >= next_break_backward) &&
1296 (global_index + d <= next_break_forward))
1297 {
1298 global_index += d;
1299 index_within_block += d;
1300 }
1301 else
1302 // outside present block, so
1303 // have to seek new block
1304 // anyway
1305 *this = Iterator(*parent, global_index + d);
1306
1307 return *this;
1308 }
1309
1310
1311
1312 template <typename BlockVectorType, bool Constness>
1313 inline Iterator<BlockVectorType, Constness> &
1314 Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1315 {
1316 // if the index pointed to is
1317 // still within the block we
1318 // currently point into, then we
1319 // can save the computation of
1320 // the block
1321 if ((global_index - d >= next_break_backward) &&
1322 (global_index - d <= next_break_forward))
1323 {
1324 global_index -= d;
1325 index_within_block -= d;
1326 }
1327 else
1328 // outside present block, so
1329 // have to seek new block
1330 // anyway
1331 *this = Iterator(*parent, global_index - d);
1332
1333 return *this;
1334 }
1335
1336
1337 template <typename BlockVectorType, bool Constness>
1338 Iterator<BlockVectorType, Constness>::Iterator(BlockVector &parent,
1339 const size_type global_index)
1340 : parent(&parent)
1341 , global_index(global_index)
1342 {
1343 // find which block we are
1344 // in. for this, take into
1345 // account that it happens at
1346 // times that people want to
1347 // initialize iterators
1348 // past-the-end
1349 if (global_index < parent.size())
1350 {
1351 const std::pair<size_type, size_type> indices =
1352 parent.block_indices.global_to_local(global_index);
1353 current_block = indices.first;
1354 index_within_block = indices.second;
1355
1356 next_break_backward =
1357 parent.block_indices.local_to_global(current_block, 0);
1358 next_break_forward =
1359 (parent.block_indices.local_to_global(current_block, 0) +
1360 parent.block_indices.block_size(current_block) - 1);
1361 }
1362 else
1363 // past the end. only have one
1364 // value for this
1365 {
1366 this->global_index = parent.size();
1367 current_block = parent.n_blocks();
1368 index_within_block = 0;
1369 next_break_backward = global_index;
1370 next_break_forward = numbers::invalid_size_type;
1371 };
1372 }
1373
1374
1375
1376 template <typename BlockVectorType, bool Constness>
1377 void
1378 Iterator<BlockVectorType, Constness>::move_forward()
1379 {
1380 if (global_index != next_break_forward)
1381 ++index_within_block;
1382 else
1383 {
1384 // ok, we traverse a boundary
1385 // between blocks:
1386 index_within_block = 0;
1387 ++current_block;
1388
1389 // break backwards is now old
1390 // break forward
1391 next_break_backward = next_break_forward + 1;
1392
1393 // compute new break forward
1394 if (current_block < parent->block_indices.size())
1395 next_break_forward +=
1396 parent->block_indices.block_size(current_block);
1397 else
1398 // if we are beyond the end,
1399 // then move the next
1400 // boundary arbitrarily far
1401 // away
1402 next_break_forward = numbers::invalid_size_type;
1403 };
1404
1405 ++global_index;
1406 }
1407
1408
1409
1410 template <typename BlockVectorType, bool Constness>
1411 void
1412 Iterator<BlockVectorType, Constness>::move_backward()
1413 {
1414 if (global_index != next_break_backward)
1415 --index_within_block;
1416 else if (current_block != 0)
1417 {
1418 // ok, we traverse a boundary
1419 // between blocks:
1420 --current_block;
1421 index_within_block =
1422 parent->block_indices.block_size(current_block) - 1;
1423
1424 // break forwards is now old
1425 // break backward
1426 next_break_forward = next_break_backward - 1;
1427
1428 // compute new break forward
1429 next_break_backward -=
1430 parent->block_indices.block_size(current_block);
1431 }
1432 else
1433 // current block was 0, we now
1434 // get into unspecified terrain
1435 {
1436 --current_block;
1437 index_within_block = numbers::invalid_size_type;
1438 next_break_forward = 0;
1439 next_break_backward = 0;
1440 };
1441
1442 --global_index;
1443 }
1444
1445
1446 } // namespace BlockVectorIterators
1447
1448} // namespace internal
1449
1450
1451
1452template <typename VectorType>
1455{
1456 return block_indices.total_size();
1457}
1458
1459
1460
1461template <typename VectorType>
1462inline std::size_t
1464{
1465 std::size_t local_size = 0;
1466 for (unsigned int b = 0; b < n_blocks(); ++b)
1467 local_size += block(b).locally_owned_size();
1468 return local_size;
1469}
1470
1471
1472
1473template <typename VectorType>
1474inline IndexSet
1476{
1477 IndexSet is(size());
1478
1479 // copy index sets from blocks into the global one, shifted
1480 // by the appropriate amount for each block
1481 for (unsigned int b = 0; b < n_blocks(); ++b)
1482 {
1483 IndexSet x = block(b).locally_owned_elements();
1485 }
1486
1487 is.compress();
1488
1489 return is;
1490}
1491
1492
1493
1494template <typename VectorType>
1495inline unsigned int
1497{
1498 return block_indices.size();
1499}
1500
1501
1502template <typename VectorType>
1504BlockVectorBase<VectorType>::block(const unsigned int i)
1505{
1507
1508 return components[i];
1509}
1510
1511
1512
1513template <typename VectorType>
1514inline const typename BlockVectorBase<VectorType>::BlockType &
1515BlockVectorBase<VectorType>::block(const unsigned int i) const
1516{
1518
1519 return components[i];
1520}
1521
1522
1523
1524template <typename VectorType>
1525inline const BlockIndices &
1527{
1528 return block_indices;
1529}
1530
1531
1532template <typename VectorType>
1533inline void
1535{
1536 std::vector<size_type> sizes(n_blocks());
1537
1538 for (size_type i = 0; i < n_blocks(); ++i)
1539 sizes[i] = block(i).size();
1540
1541 block_indices.reinit(sizes);
1542}
1543
1544
1545
1546template <typename VectorType>
1547inline void
1549{
1550 for (unsigned int i = 0; i < n_blocks(); ++i)
1551 block(i).compress(operation);
1552}
1553
1554
1555
1556template <typename VectorType>
1559{
1560 return iterator(*this, 0U);
1561}
1562
1563
1564
1565template <typename VectorType>
1568{
1569 return const_iterator(*this, 0U);
1570}
1571
1572
1573template <typename VectorType>
1576{
1577 return iterator(*this, size());
1578}
1579
1580
1581
1582template <typename VectorType>
1585{
1586 return const_iterator(*this, size());
1587}
1588
1589
1590template <typename VectorType>
1591inline bool
1593{
1594 const std::pair<size_type, size_type> local_index =
1595 block_indices.global_to_local(global_index);
1596
1597 return components[local_index.first].in_local_range(global_index);
1598}
1599
1600
1601template <typename VectorType>
1602bool
1604{
1605 for (size_type i = 0; i < n_blocks(); ++i)
1606 if (components[i].all_zero() == false)
1607 return false;
1608
1609 return true;
1610}
1611
1612
1613
1614template <typename VectorType>
1615bool
1617{
1618 for (size_type i = 0; i < n_blocks(); ++i)
1619 if (components[i].is_non_negative() == false)
1620 return false;
1621
1622 return true;
1623}
1624
1625
1626
1627template <typename VectorType>
1630 const BlockVectorBase<VectorType> &v) const
1631{
1632 Assert(n_blocks() == v.n_blocks(),
1634
1635 value_type sum = 0.;
1636 for (size_type i = 0; i < n_blocks(); ++i)
1637 sum += components[i] * v.components[i];
1638
1639 return sum;
1640}
1641
1642
1643template <typename VectorType>
1646{
1647 real_type sum = 0.;
1648 for (size_type i = 0; i < n_blocks(); ++i)
1649 sum += components[i].norm_sqr();
1650
1651 return sum;
1652}
1653
1654
1655
1656template <typename VectorType>
1659{
1660 value_type sum = 0.;
1661 // need to do static_cast as otherwise it won't work with
1662 // value_type=complex<T>
1663 for (size_type i = 0; i < n_blocks(); ++i)
1664 sum += components[i].mean_value() *
1666 components[i].size()));
1667
1669}
1670
1671
1672
1673template <typename VectorType>
1676{
1677 real_type sum = 0.;
1678 for (size_type i = 0; i < n_blocks(); ++i)
1679 sum += components[i].l1_norm();
1680
1681 return sum;
1682}
1683
1684
1685
1686template <typename VectorType>
1689{
1690 return std::sqrt(norm_sqr());
1691}
1692
1693
1694
1695template <typename VectorType>
1698{
1699 real_type sum = 0.;
1700 for (size_type i = 0; i < n_blocks(); ++i)
1701 {
1702 value_type newval = components[i].linfty_norm();
1703 if (sum < newval)
1704 sum = newval;
1705 }
1706 return sum;
1707}
1708
1709
1710
1711template <typename VectorType>
1717{
1720
1721 value_type sum = 0.;
1722 for (size_type i = 0; i < n_blocks(); ++i)
1723 sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1724
1725 return sum;
1726}
1727
1728
1729
1730template <typename VectorType>
1733{
1734 Assert(n_blocks() == v.n_blocks(),
1736
1737 for (size_type i = 0; i < n_blocks(); ++i)
1738 {
1739 components[i] += v.components[i];
1740 }
1741
1742 return *this;
1743}
1744
1745
1746
1747template <typename VectorType>
1750{
1751 Assert(n_blocks() == v.n_blocks(),
1753
1754 for (size_type i = 0; i < n_blocks(); ++i)
1755 {
1756 components[i] -= v.components[i];
1757 }
1758 return *this;
1759}
1760
1761
1762
1763template <typename VectorType>
1764template <typename Number>
1765inline void
1766BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1767 const std::vector<Number> &values)
1768{
1769 Assert(indices.size() == values.size(),
1770 ExcDimensionMismatch(indices.size(), values.size()));
1771 add(indices.size(), indices.data(), values.data());
1772}
1773
1774
1775
1776template <typename VectorType>
1777template <typename Number>
1778inline void
1779BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1780 const Vector<Number> &values)
1781{
1782 Assert(indices.size() == values.size(),
1783 ExcDimensionMismatch(indices.size(), values.size()));
1784 const size_type n_indices = indices.size();
1785 for (size_type i = 0; i < n_indices; ++i)
1786 (*this)(indices[i]) += values(i);
1787}
1788
1789
1790
1791template <typename VectorType>
1792template <typename Number>
1793inline void
1795 const size_type *indices,
1796 const Number *values)
1797{
1798 for (size_type i = 0; i < n_indices; ++i)
1799 (*this)(indices[i]) += values[i];
1800}
1801
1802
1803
1804template <typename VectorType>
1805void
1807{
1808 AssertIsFinite(a);
1809
1810 for (size_type i = 0; i < n_blocks(); ++i)
1811 {
1812 components[i].add(a);
1813 }
1814}
1815
1816
1817
1818template <typename VectorType>
1819void
1822{
1823 AssertIsFinite(a);
1824
1825 Assert(n_blocks() == v.n_blocks(),
1827
1828 for (size_type i = 0; i < n_blocks(); ++i)
1829 {
1830 components[i].add(a, v.components[i]);
1831 }
1832}
1833
1834
1835
1836template <typename VectorType>
1837void
1840 const value_type b,
1842{
1843 AssertIsFinite(a);
1844 AssertIsFinite(b);
1845
1846 Assert(n_blocks() == v.n_blocks(),
1848 Assert(n_blocks() == w.n_blocks(),
1849 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1850
1851
1852 for (size_type i = 0; i < n_blocks(); ++i)
1853 {
1854 components[i].add(a, v.components[i], b, w.components[i]);
1855 }
1856}
1857
1858
1859
1860template <typename VectorType>
1861void
1864{
1865 AssertIsFinite(x);
1866
1867 Assert(n_blocks() == v.n_blocks(),
1869
1870 for (size_type i = 0; i < n_blocks(); ++i)
1871 {
1872 components[i].sadd(x, v.components[i]);
1873 }
1874}
1875
1876
1877
1878template <typename VectorType>
1879void
1881 const value_type a,
1883{
1884 AssertIsFinite(x);
1885 AssertIsFinite(a);
1886
1887 Assert(n_blocks() == v.n_blocks(),
1889
1890 for (size_type i = 0; i < n_blocks(); ++i)
1891 {
1892 components[i].sadd(x, a, v.components[i]);
1893 }
1894}
1895
1896
1897
1898template <typename VectorType>
1899void
1901 const value_type a,
1903 const value_type b,
1905{
1906 AssertIsFinite(x);
1907 AssertIsFinite(a);
1908 AssertIsFinite(b);
1909
1910 Assert(n_blocks() == v.n_blocks(),
1912 Assert(n_blocks() == w.n_blocks(),
1913 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1914
1915 for (size_type i = 0; i < n_blocks(); ++i)
1916 {
1917 components[i].sadd(x, a, v.components[i], b, w.components[i]);
1918 }
1919}
1920
1921
1922
1923template <typename VectorType>
1924void
1926 const value_type a,
1928 const value_type b,
1930 const value_type c,
1932{
1933 AssertIsFinite(x);
1934 AssertIsFinite(a);
1935 AssertIsFinite(b);
1936 AssertIsFinite(c);
1937
1938 Assert(n_blocks() == v.n_blocks(),
1940 Assert(n_blocks() == w.n_blocks(),
1941 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1942 Assert(n_blocks() == y.n_blocks(),
1944
1945 for (size_type i = 0; i < n_blocks(); ++i)
1946 {
1947 components[i].sadd(
1948 x, a, v.components[i], b, w.components[i], c, y.components[i]);
1949 }
1950}
1951
1952
1953
1954template <typename VectorType>
1955template <class BlockVector2>
1956void
1957BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
1958{
1959 Assert(n_blocks() == v.n_blocks(),
1960 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1961 for (size_type i = 0; i < n_blocks(); ++i)
1962 components[i].scale(v.block(i));
1963}
1964
1965
1966
1967template <typename VectorType>
1968std::size_t
1970{
1971 return (MemoryConsumption::memory_consumption(this->block_indices) +
1972 MemoryConsumption::memory_consumption(this->components));
1973}
1974
1975
1976
1977template <typename VectorType>
1978template <class BlockVector2>
1979void
1980BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
1981{
1982 AssertIsFinite(a);
1983
1984 Assert(n_blocks() == v.n_blocks(),
1985 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1986
1987 for (size_type i = 0; i < n_blocks(); ++i)
1988 components[i].equ(a, v.components[i]);
1989}
1990
1991
1992
1993template <typename VectorType>
1994void
1996{
1997 for (size_type i = 0; i < n_blocks(); ++i)
1998 block(i).update_ghost_values();
1999}
2000
2001
2002
2003template <typename VectorType>
2006{
2007 if (n_blocks() > 0)
2008 return block(0).get_mpi_communicator();
2009 else
2010 return MPI_COMM_SELF;
2011}
2012
2013
2014
2015template <typename VectorType>
2018{
2019 AssertIsFinite(s);
2020
2021 for (size_type i = 0; i < n_blocks(); ++i)
2022 components[i] = s;
2023
2024 return *this;
2025}
2026
2027
2028template <typename VectorType>
2031{
2033
2034 for (size_type i = 0; i < n_blocks(); ++i)
2035 components[i] = v.components[i];
2036
2037 return *this;
2038}
2039
2040
2041template <typename VectorType>
2042template <typename VectorType2>
2045{
2047
2048 for (size_type i = 0; i < n_blocks(); ++i)
2049 components[i] = v.components[i];
2050
2051 return *this;
2052}
2053
2054
2055
2056template <typename VectorType>
2058BlockVectorBase<VectorType>::operator=(const VectorType &v)
2059{
2060 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2061
2062 size_type index_v = 0;
2063 for (size_type b = 0; b < n_blocks(); ++b)
2064 for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2065 block(b)(i) = v(index_v);
2066
2067 return *this;
2068}
2069
2070
2071
2072template <typename VectorType>
2073template <typename VectorType2>
2074inline bool
2076 const BlockVectorBase<VectorType2> &v) const
2077{
2079
2080 for (size_type i = 0; i < n_blocks(); ++i)
2081 if (!(components[i] == v.components[i]))
2082 return false;
2083
2084 return true;
2085}
2086
2087
2088
2089template <typename VectorType>
2092{
2093 AssertIsFinite(factor);
2094
2095 for (size_type i = 0; i < n_blocks(); ++i)
2096 components[i] *= factor;
2097
2098 return *this;
2099}
2100
2101
2102
2103template <typename VectorType>
2106{
2107 AssertIsFinite(factor);
2108 Assert(factor != 0., ExcDivideByZero());
2109
2110 const value_type inverse_factor = value_type(1.) / factor;
2111
2112 for (size_type i = 0; i < n_blocks(); ++i)
2113 components[i] *= inverse_factor;
2114
2115 return *this;
2116}
2117
2118
2119template <typename VectorType>
2122{
2123 const std::pair<unsigned int, size_type> local_index =
2125 return components[local_index.first](local_index.second);
2126}
2127
2128
2129
2130template <typename VectorType>
2133{
2134 const std::pair<unsigned int, size_type> local_index =
2136 return components[local_index.first](local_index.second);
2137}
2138
2139
2140
2141template <typename VectorType>
2144{
2145 return operator()(i);
2146}
2147
2148
2149
2150template <typename VectorType>
2153{
2154 return operator()(i);
2155}
2156
2157
2158
2159template <typename VectorType>
2160template <typename OtherNumber>
2161inline void
2163 const std::vector<size_type> &indices,
2164 std::vector<OtherNumber> &values) const
2165{
2166 for (size_type i = 0; i < indices.size(); ++i)
2167 values[i] = operator()(indices[i]);
2168}
2169
2170
2171
2172template <typename VectorType>
2173inline void
2176 const ArrayView<value_type> &entries) const
2177{
2178 AssertDimension(indices.size(), entries.size());
2179 for (unsigned int i = 0; i < indices.size(); ++i)
2180 {
2181 AssertIndexRange(indices[i], size());
2182 entries[i] = (*this)[indices[i]];
2183 }
2184}
2185
2186
2187
2188template <typename VectorType>
2189template <typename ForwardIterator, typename OutputIterator>
2190inline void
2192 ForwardIterator indices_begin,
2193 const ForwardIterator indices_end,
2194 OutputIterator values_begin) const
2195{
2196 while (indices_begin != indices_end)
2197 {
2198 *values_begin = operator()(*indices_begin);
2199 ++indices_begin;
2200 ++values_begin;
2201 }
2202}
2203
2204#endif // DOXYGEN
2205
2207
2208#endif
std::size_t size() const
Definition array_view.h:684
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:1831
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#define Assert(cond, exc)
static ::ExceptionBase & ExcPointerToDifferentVectors()
#define AssertIsFinite(number)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:489
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcDifferentBlockIndices()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static const bool value
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
const types::global_dof_index invalid_size_type
Definition types.h:233
STL namespace.
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition types.h:32
unsigned int global_dof_index
Definition types.h:81