Reference documentation for deal.II version GIT relicensing-426-g7976cfd195 2024-04-18 21:10:01+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
25
29#include <deal.II/lac/vector.h>
31
32#include <cmath>
33#include <cstddef>
34#include <iterator>
35#include <type_traits>
36#include <vector>
37
39
40
46namespace internal
47{
48 template <typename T>
49 using has_block_t = decltype(std::declval<const T>().block(0));
50
51 template <typename T>
52 constexpr bool has_block = internal::is_supported_operation<has_block_t, T>;
53
54 template <typename T>
55 using has_n_blocks_t = decltype(std::declval<const T>().n_blocks());
56
57 template <typename T>
58 constexpr bool has_n_blocks =
59 internal::is_supported_operation<has_n_blocks_t, T>;
60
61 template <typename T>
62 constexpr bool is_block_vector = has_block<T> && has_n_blocks<T>;
63} // namespace internal
64
79template <typename VectorType>
81{
82public:
88 static const bool value = internal::is_block_vector<VectorType>;
89};
90
91
92// instantiation of the static member
93template <typename VectorType>
95
96
97
98namespace internal
99{
103 namespace BlockVectorIterators
104 {
120 template <typename BlockVectorType, bool Constness>
122 {
123 public:
128
135 std::conditional_t<Constness,
136 const typename BlockVectorType::value_type,
137 typename BlockVectorType::value_type>;
138
152 using iterator_category = std::random_access_iterator_tag;
153 using difference_type = std::ptrdiff_t;
154 using reference = typename BlockVectorType::reference;
156
158 std::conditional_t<Constness,
160 typename BlockVectorType::BlockType::reference>;
161
167 std::conditional_t<Constness, const BlockVectorType, BlockVectorType>;
168
178
187
188
193
194 private:
205
206 public:
210 Iterator &
212
219 operator*() const;
220
227
232 Iterator &
234
242
247 Iterator &
249
257
262 template <bool OtherConstness>
263 bool
265
270 template <bool OtherConstness>
271 bool
273
279 template <bool OtherConstness>
280 bool
282
286 template <bool OtherConstness>
287 bool
289
293 template <bool OtherConstness>
294 bool
296
300 template <bool OtherConstness>
301 bool
303
307 template <bool OtherConstness>
310
316 operator+(const difference_type &d) const;
317
323 operator-(const difference_type &d) const;
324
329 Iterator &
331
336 Iterator &
338
349 "Your program tried to compare iterators pointing to "
350 "different block vectors. There is no reasonable way "
351 "to do this.");
352
354 private:
361
366
371 unsigned int current_block;
373
382
386 void
388
392 void
394
395 // Mark all other instances of this template as friends.
396 template <typename, bool>
397 friend class Iterator;
398 };
399 } // namespace BlockVectorIterators
400} // namespace internal
401
402
440template <typename VectorType>
442 public ReadVector<typename VectorType::value_type>
443{
444public:
448 using BlockType = VectorType;
449
450 /*
451 * Declare standard types used in
452 * all containers. These types
453 * parallel those in the
454 * <tt>C++</tt> standard
455 * libraries
456 * <tt>std::vector<...></tt>
457 * class. This includes iterator
458 * types.
459 */
460 using value_type = typename BlockType::value_type;
462 using const_pointer = const value_type *;
463 using iterator =
467 using reference = typename BlockType::reference;
468 using const_reference = typename BlockType::const_reference;
470
480 using real_type = typename BlockType::real_type;
481
485 BlockVectorBase() = default;
486
490 BlockVectorBase(const BlockVectorBase & /*V*/) = default;
491
497 BlockVectorBase(BlockVectorBase && /*V*/) noexcept = default;
498
505 void
507
518 void
519 compress(VectorOperation::values operation);
520
524 BlockType &
525 block(const unsigned int i);
526
530 const BlockType &
531 block(const unsigned int i) const;
532
538 const BlockIndices &
540
544 unsigned int
545 n_blocks() const;
546
551 virtual size_type
552 size() const override;
553
559 std::size_t
561
580
586
592 begin() const;
593
599
605 end() const;
606
611 operator()(const size_type i) const;
612
617 operator()(const size_type i);
618
625 operator[](const size_type i) const;
626
633 operator[](const size_type i);
634
650 template <typename OtherNumber>
651 void
652 extract_subvector_to(const std::vector<size_type> &indices,
653 std::vector<OtherNumber> &values) const;
654
655 virtual void
656 extract_subvector_to(const ArrayView<const types::global_dof_index> &indices,
657 ArrayView<value_type> &entries) const override;
658
686 template <typename ForwardIterator, typename OutputIterator>
687 void
688 extract_subvector_to(ForwardIterator indices_begin,
689 const ForwardIterator indices_end,
690 OutputIterator values_begin) const;
691
697 operator=(const value_type s);
698
703 operator=(const BlockVectorBase &V);
704
711 operator=(BlockVectorBase && /*V*/) = default; // NOLINT
712
716 template <typename VectorType2>
718 operator=(const BlockVectorBase<VectorType2> &V);
719
724 operator=(const VectorType &v);
725
730 template <typename VectorType2>
731 bool
732 operator==(const BlockVectorBase<VectorType2> &v) const;
733
738 operator*(const BlockVectorBase &V) const;
739
744 norm_sqr() const;
745
750 mean_value() const;
751
756 l1_norm() const;
757
763 l2_norm() const;
764
770 linfty_norm() const;
771
796 const BlockVectorBase &V,
797 const BlockVectorBase &W);
798
803 bool
804 in_local_range(const size_type global_index) const;
805
811 bool
812 all_zero() const;
813
819 bool
821
826 operator+=(const BlockVectorBase &V);
827
832 operator-=(const BlockVectorBase &V);
833
834
839 template <typename Number>
840 void
841 add(const std::vector<size_type> &indices, const std::vector<Number> &values);
842
847 template <typename Number>
848 void
849 add(const std::vector<size_type> &indices, const Vector<Number> &values);
850
856 template <typename Number>
857 void
858 add(const size_type n_elements,
859 const size_type *indices,
860 const Number *values);
861
866 void
867 add(const value_type s);
868
872 void
873 add(const value_type a, const BlockVectorBase &V);
874
878 void
880 const BlockVectorBase &V,
881 const value_type b,
882 const BlockVectorBase &W);
883
887 void
888 sadd(const value_type s, const BlockVectorBase &V);
889
893 void
894 sadd(const value_type s, const value_type a, const BlockVectorBase &V);
895
899 void
901 const value_type a,
902 const BlockVectorBase &V,
903 const value_type b,
904 const BlockVectorBase &W);
905
909 void
911 const value_type a,
912 const BlockVectorBase &V,
913 const value_type b,
914 const BlockVectorBase &W,
915 const value_type c,
916 const BlockVectorBase &X);
917
922 operator*=(const value_type factor);
923
928 operator/=(const value_type factor);
929
934 template <class BlockVector2>
935 void
936 scale(const BlockVector2 &v);
937
941 template <class BlockVector2>
942 void
943 equ(const value_type a, const BlockVector2 &V);
944
949 void
951
956 std::size_t
958
959protected:
963 std::vector<VectorType> components;
964
970
971 // Make the iterator class a friend.
972 template <typename N, bool C>
973 friend class ::internal::BlockVectorIterators::Iterator;
974
975 template <typename>
976 friend class BlockVectorBase;
977};
978
979
982/*----------------------- Inline functions ----------------------------------*/
983
984
985#ifndef DOXYGEN
986namespace internal
987{
988 namespace BlockVectorIterators
989 {
990 template <typename BlockVectorType, bool Constness>
991 inline Iterator<BlockVectorType, Constness>::Iterator(
992 const Iterator<BlockVectorType, Constness> &c)
993 : parent(c.parent)
994 , global_index(c.global_index)
995 , current_block(c.current_block)
996 , index_within_block(c.index_within_block)
997 , next_break_forward(c.next_break_forward)
998 , next_break_backward(c.next_break_backward)
999 {}
1000
1001
1002
1003 template <typename BlockVectorType, bool Constness>
1004 inline Iterator<BlockVectorType, Constness>::Iterator(
1005 const Iterator<BlockVectorType, !Constness> &c)
1006 : parent(c.parent)
1007 , global_index(c.global_index)
1008 , current_block(c.current_block)
1009 , index_within_block(c.index_within_block)
1010 , next_break_forward(c.next_break_forward)
1011 , next_break_backward(c.next_break_backward)
1012 {
1013 // Only permit copy-constructing const iterators from non-const
1014 // iterators, and not vice versa (i.e., Constness must always be
1015 // true).
1016 static_assert(Constness == true,
1017 "Constructing a non-const iterator from a const iterator "
1018 "does not make sense.");
1019 }
1020
1021
1022
1023 template <typename BlockVectorType, bool Constness>
1024 inline Iterator<BlockVectorType, Constness>::Iterator(
1025 BlockVector &parent,
1026 const size_type global_index,
1027 const size_type current_block,
1028 const size_type index_within_block,
1029 const size_type next_break_forward,
1030 const size_type next_break_backward)
1031 : parent(&parent)
1032 , global_index(global_index)
1033 , current_block(current_block)
1034 , index_within_block(index_within_block)
1035 , next_break_forward(next_break_forward)
1036 , next_break_backward(next_break_backward)
1037 {}
1038
1039
1040
1041 template <typename BlockVectorType, bool Constness>
1042 inline Iterator<BlockVectorType, Constness> &
1043 Iterator<BlockVectorType, Constness>::operator=(const Iterator &c)
1044 {
1045 parent = c.parent;
1046 global_index = c.global_index;
1047 index_within_block = c.index_within_block;
1048 current_block = c.current_block;
1049 next_break_forward = c.next_break_forward;
1050 next_break_backward = c.next_break_backward;
1051
1052 return *this;
1053 }
1054
1055
1056
1057 template <typename BlockVectorType, bool Constness>
1058 inline typename Iterator<BlockVectorType, Constness>::dereference_type
1059 Iterator<BlockVectorType, Constness>::operator*() const
1060 {
1061 return parent->block(current_block)(index_within_block);
1062 }
1063
1064
1065
1066 template <typename BlockVectorType, bool Constness>
1067 inline typename Iterator<BlockVectorType, Constness>::dereference_type
1068 Iterator<BlockVectorType, Constness>::operator[](
1069 const difference_type d) const
1070 {
1071 // if the index pointed to is
1072 // still within the block we
1073 // currently point into, then we
1074 // can save the computation of
1075 // the block
1076 if ((global_index + d >= next_break_backward) &&
1077 (global_index + d <= next_break_forward))
1078 return parent->block(current_block)(index_within_block + d);
1079
1080 // if the index is not within the
1081 // block of the block vector into
1082 // which we presently point, then
1083 // there is no way: we have to
1084 // search for the block. this can
1085 // be done through the parent
1086 // class as well.
1087 return (*parent)(global_index + d);
1088 }
1089
1090
1091
1092 template <typename BlockVectorType, bool Constness>
1093 inline Iterator<BlockVectorType, Constness> &
1094 Iterator<BlockVectorType, Constness>::operator++()
1095 {
1096 move_forward();
1097 return *this;
1098 }
1099
1100
1101
1102 template <typename BlockVectorType, bool Constness>
1103 inline Iterator<BlockVectorType, Constness>
1104 Iterator<BlockVectorType, Constness>::operator++(int)
1105 {
1106 const Iterator old_value = *this;
1107 move_forward();
1108 return old_value;
1109 }
1110
1111
1112
1113 template <typename BlockVectorType, bool Constness>
1114 inline Iterator<BlockVectorType, Constness> &
1115 Iterator<BlockVectorType, Constness>::operator--()
1116 {
1117 move_backward();
1118 return *this;
1119 }
1120
1121
1122
1123 template <typename BlockVectorType, bool Constness>
1124 inline Iterator<BlockVectorType, Constness>
1125 Iterator<BlockVectorType, Constness>::operator--(int)
1126 {
1127 const Iterator old_value = *this;
1128 move_backward();
1129 return old_value;
1130 }
1131
1132
1133
1134 template <typename BlockVectorType, bool Constness>
1135 template <bool OtherConstness>
1136 inline bool
1137 Iterator<BlockVectorType, Constness>::operator==(
1138 const Iterator<BlockVectorType, OtherConstness> &i) const
1139 {
1140 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1141
1142 return (global_index == i.global_index);
1143 }
1144
1145
1146
1147 template <typename BlockVectorType, bool Constness>
1148 template <bool OtherConstness>
1149 inline bool
1150 Iterator<BlockVectorType, Constness>::operator!=(
1151 const Iterator<BlockVectorType, OtherConstness> &i) const
1152 {
1153 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1154
1155 return (global_index != i.global_index);
1156 }
1157
1158
1159
1160 template <typename BlockVectorType, bool Constness>
1161 template <bool OtherConstness>
1162 inline bool
1163 Iterator<BlockVectorType, Constness>::operator<(
1164 const Iterator<BlockVectorType, OtherConstness> &i) const
1165 {
1166 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1167
1168 return (global_index < i.global_index);
1169 }
1170
1171
1172
1173 template <typename BlockVectorType, bool Constness>
1174 template <bool OtherConstness>
1175 inline bool
1176 Iterator<BlockVectorType, Constness>::operator<=(
1177 const Iterator<BlockVectorType, OtherConstness> &i) const
1178 {
1179 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1180
1181 return (global_index <= i.global_index);
1182 }
1183
1184
1185
1186 template <typename BlockVectorType, bool Constness>
1187 template <bool OtherConstness>
1188 inline bool
1189 Iterator<BlockVectorType, Constness>::operator>(
1190 const Iterator<BlockVectorType, OtherConstness> &i) const
1191 {
1192 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1193
1194 return (global_index > i.global_index);
1195 }
1196
1197
1198
1199 template <typename BlockVectorType, bool Constness>
1200 template <bool OtherConstness>
1201 inline bool
1202 Iterator<BlockVectorType, Constness>::operator>=(
1203 const Iterator<BlockVectorType, OtherConstness> &i) const
1204 {
1205 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1206
1207 return (global_index >= i.global_index);
1208 }
1209
1210
1211
1212 template <typename BlockVectorType, bool Constness>
1213 template <bool OtherConstness>
1214 inline typename Iterator<BlockVectorType, Constness>::difference_type
1215 Iterator<BlockVectorType, Constness>::operator-(
1216 const Iterator<BlockVectorType, OtherConstness> &i) const
1217 {
1218 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1219
1220 return (static_cast<signed int>(global_index) -
1221 static_cast<signed int>(i.global_index));
1222 }
1223
1224
1225
1226 template <typename BlockVectorType, bool Constness>
1227 inline Iterator<BlockVectorType, Constness>
1228 Iterator<BlockVectorType, Constness>::operator+(
1229 const difference_type &d) const
1230 {
1231 // if the index pointed to is
1232 // still within the block we
1233 // currently point into, then we
1234 // can save the computation of
1235 // the block
1236 if ((global_index + d >= next_break_backward) &&
1237 (global_index + d <= next_break_forward))
1238 return Iterator(*parent,
1239 global_index + d,
1240 current_block,
1241 index_within_block + d,
1242 next_break_forward,
1243 next_break_backward);
1244 else
1245 // outside present block, so
1246 // have to seek new block
1247 // anyway
1248 return Iterator(*parent, global_index + d);
1249 }
1250
1251
1252
1253 template <typename BlockVectorType, bool Constness>
1254 inline Iterator<BlockVectorType, Constness>
1255 Iterator<BlockVectorType, Constness>::operator-(
1256 const difference_type &d) const
1257 {
1258 // if the index pointed to is
1259 // still within the block we
1260 // currently point into, then we
1261 // can save the computation of
1262 // the block
1263 if ((global_index - d >= next_break_backward) &&
1264 (global_index - d <= next_break_forward))
1265 return Iterator(*parent,
1266 global_index - d,
1267 current_block,
1268 index_within_block - d,
1269 next_break_forward,
1270 next_break_backward);
1271 else
1272 // outside present block, so
1273 // have to seek new block
1274 // anyway
1275 return Iterator(*parent, global_index - d);
1276 }
1277
1278
1279
1280 template <typename BlockVectorType, bool Constness>
1281 inline Iterator<BlockVectorType, Constness> &
1282 Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1283 {
1284 // if the index pointed to is
1285 // still within the block we
1286 // currently point into, then we
1287 // can save the computation of
1288 // the block
1289 if ((global_index + d >= next_break_backward) &&
1290 (global_index + d <= next_break_forward))
1291 {
1292 global_index += d;
1293 index_within_block += d;
1294 }
1295 else
1296 // outside present block, so
1297 // have to seek new block
1298 // anyway
1299 *this = Iterator(*parent, global_index + d);
1300
1301 return *this;
1302 }
1303
1304
1305
1306 template <typename BlockVectorType, bool Constness>
1307 inline Iterator<BlockVectorType, Constness> &
1308 Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1309 {
1310 // if the index pointed to is
1311 // still within the block we
1312 // currently point into, then we
1313 // can save the computation of
1314 // the block
1315 if ((global_index - d >= next_break_backward) &&
1316 (global_index - d <= next_break_forward))
1317 {
1318 global_index -= d;
1319 index_within_block -= d;
1320 }
1321 else
1322 // outside present block, so
1323 // have to seek new block
1324 // anyway
1325 *this = Iterator(*parent, global_index - d);
1326
1327 return *this;
1328 }
1329
1330
1331 template <typename BlockVectorType, bool Constness>
1332 Iterator<BlockVectorType, Constness>::Iterator(BlockVector &parent,
1333 const size_type global_index)
1334 : parent(&parent)
1335 , global_index(global_index)
1336 {
1337 // find which block we are
1338 // in. for this, take into
1339 // account that it happens at
1340 // times that people want to
1341 // initialize iterators
1342 // past-the-end
1343 if (global_index < parent.size())
1344 {
1345 const std::pair<size_type, size_type> indices =
1346 parent.block_indices.global_to_local(global_index);
1347 current_block = indices.first;
1348 index_within_block = indices.second;
1349
1350 next_break_backward =
1351 parent.block_indices.local_to_global(current_block, 0);
1352 next_break_forward =
1353 (parent.block_indices.local_to_global(current_block, 0) +
1354 parent.block_indices.block_size(current_block) - 1);
1355 }
1356 else
1357 // past the end. only have one
1358 // value for this
1359 {
1360 this->global_index = parent.size();
1361 current_block = parent.n_blocks();
1362 index_within_block = 0;
1363 next_break_backward = global_index;
1364 next_break_forward = numbers::invalid_size_type;
1365 };
1366 }
1367
1368
1369
1370 template <typename BlockVectorType, bool Constness>
1371 void
1372 Iterator<BlockVectorType, Constness>::move_forward()
1373 {
1374 if (global_index != next_break_forward)
1375 ++index_within_block;
1376 else
1377 {
1378 // ok, we traverse a boundary
1379 // between blocks:
1380 index_within_block = 0;
1381 ++current_block;
1382
1383 // break backwards is now old
1384 // break forward
1385 next_break_backward = next_break_forward + 1;
1386
1387 // compute new break forward
1388 if (current_block < parent->block_indices.size())
1389 next_break_forward +=
1390 parent->block_indices.block_size(current_block);
1391 else
1392 // if we are beyond the end,
1393 // then move the next
1394 // boundary arbitrarily far
1395 // away
1396 next_break_forward = numbers::invalid_size_type;
1397 };
1398
1399 ++global_index;
1400 }
1401
1402
1403
1404 template <typename BlockVectorType, bool Constness>
1405 void
1406 Iterator<BlockVectorType, Constness>::move_backward()
1407 {
1408 if (global_index != next_break_backward)
1409 --index_within_block;
1410 else if (current_block != 0)
1411 {
1412 // ok, we traverse a boundary
1413 // between blocks:
1414 --current_block;
1415 index_within_block =
1416 parent->block_indices.block_size(current_block) - 1;
1417
1418 // break forwards is now old
1419 // break backward
1420 next_break_forward = next_break_backward - 1;
1421
1422 // compute new break forward
1423 next_break_backward -=
1424 parent->block_indices.block_size(current_block);
1425 }
1426 else
1427 // current block was 0, we now
1428 // get into unspecified terrain
1429 {
1430 --current_block;
1431 index_within_block = numbers::invalid_size_type;
1432 next_break_forward = 0;
1433 next_break_backward = 0;
1434 };
1435
1436 --global_index;
1437 }
1438
1439
1440 } // namespace BlockVectorIterators
1441
1442} // namespace internal
1443
1444
1445
1446template <typename VectorType>
1449{
1450 return block_indices.total_size();
1451}
1452
1453
1454
1455template <typename VectorType>
1456inline std::size_t
1458{
1459 std::size_t local_size = 0;
1460 for (unsigned int b = 0; b < n_blocks(); ++b)
1461 local_size += block(b).locally_owned_size();
1462 return local_size;
1463}
1464
1465
1466
1467template <typename VectorType>
1468inline IndexSet
1470{
1471 IndexSet is(size());
1472
1473 // copy index sets from blocks into the global one, shifted
1474 // by the appropriate amount for each block
1475 for (unsigned int b = 0; b < n_blocks(); ++b)
1476 {
1477 IndexSet x = block(b).locally_owned_elements();
1479 }
1480
1481 is.compress();
1482
1483 return is;
1484}
1485
1486
1487
1488template <typename VectorType>
1489inline unsigned int
1491{
1492 return block_indices.size();
1493}
1494
1495
1496template <typename VectorType>
1498BlockVectorBase<VectorType>::block(const unsigned int i)
1499{
1501
1502 return components[i];
1503}
1504
1505
1506
1507template <typename VectorType>
1508inline const typename BlockVectorBase<VectorType>::BlockType &
1509BlockVectorBase<VectorType>::block(const unsigned int i) const
1510{
1512
1513 return components[i];
1514}
1515
1516
1517
1518template <typename VectorType>
1519inline const BlockIndices &
1521{
1522 return block_indices;
1523}
1524
1525
1526template <typename VectorType>
1527inline void
1529{
1530 std::vector<size_type> sizes(n_blocks());
1531
1532 for (size_type i = 0; i < n_blocks(); ++i)
1533 sizes[i] = block(i).size();
1534
1535 block_indices.reinit(sizes);
1536}
1537
1538
1539
1540template <typename VectorType>
1541inline void
1543{
1544 for (unsigned int i = 0; i < n_blocks(); ++i)
1545 block(i).compress(operation);
1546}
1547
1548
1549
1550template <typename VectorType>
1553{
1554 return iterator(*this, 0U);
1555}
1556
1557
1558
1559template <typename VectorType>
1562{
1563 return const_iterator(*this, 0U);
1564}
1565
1566
1567template <typename VectorType>
1570{
1571 return iterator(*this, size());
1572}
1573
1574
1575
1576template <typename VectorType>
1579{
1580 return const_iterator(*this, size());
1581}
1582
1583
1584template <typename VectorType>
1585inline bool
1587{
1588 const std::pair<size_type, size_type> local_index =
1589 block_indices.global_to_local(global_index);
1590
1591 return components[local_index.first].in_local_range(global_index);
1592}
1593
1594
1595template <typename VectorType>
1596bool
1598{
1599 for (size_type i = 0; i < n_blocks(); ++i)
1600 if (components[i].all_zero() == false)
1601 return false;
1602
1603 return true;
1604}
1605
1606
1607
1608template <typename VectorType>
1609bool
1611{
1612 for (size_type i = 0; i < n_blocks(); ++i)
1613 if (components[i].is_non_negative() == false)
1614 return false;
1615
1616 return true;
1617}
1618
1619
1620
1621template <typename VectorType>
1624 const BlockVectorBase<VectorType> &v) const
1625{
1626 Assert(n_blocks() == v.n_blocks(),
1628
1629 value_type sum = 0.;
1630 for (size_type i = 0; i < n_blocks(); ++i)
1631 sum += components[i] * v.components[i];
1632
1633 return sum;
1634}
1635
1636
1637template <typename VectorType>
1640{
1641 real_type sum = 0.;
1642 for (size_type i = 0; i < n_blocks(); ++i)
1643 sum += components[i].norm_sqr();
1644
1645 return sum;
1646}
1647
1648
1649
1650template <typename VectorType>
1653{
1654 value_type sum = 0.;
1655 // need to do static_cast as otherwise it won't work with
1656 // value_type=complex<T>
1657 for (size_type i = 0; i < n_blocks(); ++i)
1658 sum += components[i].mean_value() *
1660 components[i].size()));
1661
1663}
1664
1665
1666
1667template <typename VectorType>
1670{
1671 real_type sum = 0.;
1672 for (size_type i = 0; i < n_blocks(); ++i)
1673 sum += components[i].l1_norm();
1674
1675 return sum;
1676}
1677
1678
1679
1680template <typename VectorType>
1683{
1684 return std::sqrt(norm_sqr());
1685}
1686
1687
1688
1689template <typename VectorType>
1692{
1693 real_type sum = 0.;
1694 for (size_type i = 0; i < n_blocks(); ++i)
1695 {
1696 value_type newval = components[i].linfty_norm();
1697 if (sum < newval)
1698 sum = newval;
1699 }
1700 return sum;
1701}
1702
1703
1704
1705template <typename VectorType>
1711{
1714
1715 value_type sum = 0.;
1716 for (size_type i = 0; i < n_blocks(); ++i)
1717 sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1718
1719 return sum;
1720}
1721
1722
1723
1724template <typename VectorType>
1727{
1728 Assert(n_blocks() == v.n_blocks(),
1730
1731 for (size_type i = 0; i < n_blocks(); ++i)
1732 {
1733 components[i] += v.components[i];
1734 }
1735
1736 return *this;
1737}
1738
1739
1740
1741template <typename VectorType>
1744{
1745 Assert(n_blocks() == v.n_blocks(),
1747
1748 for (size_type i = 0; i < n_blocks(); ++i)
1749 {
1750 components[i] -= v.components[i];
1751 }
1752 return *this;
1753}
1754
1755
1756
1757template <typename VectorType>
1758template <typename Number>
1759inline void
1760BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1761 const std::vector<Number> &values)
1762{
1763 Assert(indices.size() == values.size(),
1764 ExcDimensionMismatch(indices.size(), values.size()));
1765 add(indices.size(), indices.data(), values.data());
1766}
1767
1768
1769
1770template <typename VectorType>
1771template <typename Number>
1772inline void
1773BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1774 const Vector<Number> &values)
1775{
1776 Assert(indices.size() == values.size(),
1777 ExcDimensionMismatch(indices.size(), values.size()));
1778 const size_type n_indices = indices.size();
1779 for (size_type i = 0; i < n_indices; ++i)
1780 (*this)(indices[i]) += values(i);
1781}
1782
1783
1784
1785template <typename VectorType>
1786template <typename Number>
1787inline void
1789 const size_type *indices,
1790 const Number *values)
1791{
1792 for (size_type i = 0; i < n_indices; ++i)
1793 (*this)(indices[i]) += values[i];
1794}
1795
1796
1797
1798template <typename VectorType>
1799void
1801{
1802 AssertIsFinite(a);
1803
1804 for (size_type i = 0; i < n_blocks(); ++i)
1805 {
1806 components[i].add(a);
1807 }
1808}
1809
1810
1811
1812template <typename VectorType>
1813void
1816{
1817 AssertIsFinite(a);
1818
1819 Assert(n_blocks() == v.n_blocks(),
1821
1822 for (size_type i = 0; i < n_blocks(); ++i)
1823 {
1824 components[i].add(a, v.components[i]);
1825 }
1826}
1827
1828
1829
1830template <typename VectorType>
1831void
1834 const value_type b,
1836{
1837 AssertIsFinite(a);
1838 AssertIsFinite(b);
1839
1840 Assert(n_blocks() == v.n_blocks(),
1842 Assert(n_blocks() == w.n_blocks(),
1843 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1844
1845
1846 for (size_type i = 0; i < n_blocks(); ++i)
1847 {
1848 components[i].add(a, v.components[i], b, w.components[i]);
1849 }
1850}
1851
1852
1853
1854template <typename VectorType>
1855void
1858{
1859 AssertIsFinite(x);
1860
1861 Assert(n_blocks() == v.n_blocks(),
1863
1864 for (size_type i = 0; i < n_blocks(); ++i)
1865 {
1866 components[i].sadd(x, v.components[i]);
1867 }
1868}
1869
1870
1871
1872template <typename VectorType>
1873void
1875 const value_type a,
1877{
1878 AssertIsFinite(x);
1879 AssertIsFinite(a);
1880
1881 Assert(n_blocks() == v.n_blocks(),
1883
1884 for (size_type i = 0; i < n_blocks(); ++i)
1885 {
1886 components[i].sadd(x, a, v.components[i]);
1887 }
1888}
1889
1890
1891
1892template <typename VectorType>
1893void
1895 const value_type a,
1897 const value_type b,
1899{
1900 AssertIsFinite(x);
1901 AssertIsFinite(a);
1902 AssertIsFinite(b);
1903
1904 Assert(n_blocks() == v.n_blocks(),
1906 Assert(n_blocks() == w.n_blocks(),
1907 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1908
1909 for (size_type i = 0; i < n_blocks(); ++i)
1910 {
1911 components[i].sadd(x, a, v.components[i], b, w.components[i]);
1912 }
1913}
1914
1915
1916
1917template <typename VectorType>
1918void
1920 const value_type a,
1922 const value_type b,
1924 const value_type c,
1926{
1927 AssertIsFinite(x);
1928 AssertIsFinite(a);
1929 AssertIsFinite(b);
1930 AssertIsFinite(c);
1931
1932 Assert(n_blocks() == v.n_blocks(),
1934 Assert(n_blocks() == w.n_blocks(),
1935 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1936 Assert(n_blocks() == y.n_blocks(),
1938
1939 for (size_type i = 0; i < n_blocks(); ++i)
1940 {
1941 components[i].sadd(
1942 x, a, v.components[i], b, w.components[i], c, y.components[i]);
1943 }
1944}
1945
1946
1947
1948template <typename VectorType>
1949template <class BlockVector2>
1950void
1951BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
1952{
1953 Assert(n_blocks() == v.n_blocks(),
1954 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1955 for (size_type i = 0; i < n_blocks(); ++i)
1956 components[i].scale(v.block(i));
1957}
1958
1959
1960
1961template <typename VectorType>
1962std::size_t
1964{
1965 return (MemoryConsumption::memory_consumption(this->block_indices) +
1966 MemoryConsumption::memory_consumption(this->components));
1967}
1968
1969
1970
1971template <typename VectorType>
1972template <class BlockVector2>
1973void
1974BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
1975{
1976 AssertIsFinite(a);
1977
1978 Assert(n_blocks() == v.n_blocks(),
1979 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1980
1981 for (size_type i = 0; i < n_blocks(); ++i)
1982 components[i].equ(a, v.components[i]);
1983}
1984
1985
1986
1987template <typename VectorType>
1988void
1990{
1991 for (size_type i = 0; i < n_blocks(); ++i)
1992 block(i).update_ghost_values();
1993}
1994
1995
1996
1997template <typename VectorType>
2000{
2001 AssertIsFinite(s);
2002
2003 for (size_type i = 0; i < n_blocks(); ++i)
2004 components[i] = s;
2005
2006 return *this;
2007}
2008
2009
2010template <typename VectorType>
2013{
2015
2016 for (size_type i = 0; i < n_blocks(); ++i)
2017 components[i] = v.components[i];
2018
2019 return *this;
2020}
2021
2022
2023template <typename VectorType>
2024template <typename VectorType2>
2027{
2029
2030 for (size_type i = 0; i < n_blocks(); ++i)
2031 components[i] = v.components[i];
2032
2033 return *this;
2034}
2035
2036
2037
2038template <typename VectorType>
2040BlockVectorBase<VectorType>::operator=(const VectorType &v)
2041{
2042 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2043
2044 size_type index_v = 0;
2045 for (size_type b = 0; b < n_blocks(); ++b)
2046 for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2047 block(b)(i) = v(index_v);
2048
2049 return *this;
2050}
2051
2052
2053
2054template <typename VectorType>
2055template <typename VectorType2>
2056inline bool
2058 const BlockVectorBase<VectorType2> &v) const
2059{
2061
2062 for (size_type i = 0; i < n_blocks(); ++i)
2063 if (!(components[i] == v.components[i]))
2064 return false;
2065
2066 return true;
2067}
2068
2069
2070
2071template <typename VectorType>
2074{
2075 AssertIsFinite(factor);
2076
2077 for (size_type i = 0; i < n_blocks(); ++i)
2078 components[i] *= factor;
2079
2080 return *this;
2081}
2082
2083
2084
2085template <typename VectorType>
2088{
2089 AssertIsFinite(factor);
2090 Assert(factor != 0., ExcDivideByZero());
2091
2092 const value_type inverse_factor = value_type(1.) / factor;
2093
2094 for (size_type i = 0; i < n_blocks(); ++i)
2095 components[i] *= inverse_factor;
2096
2097 return *this;
2098}
2099
2100
2101template <typename VectorType>
2104{
2105 const std::pair<unsigned int, size_type> local_index =
2107 return components[local_index.first](local_index.second);
2108}
2109
2110
2111
2112template <typename VectorType>
2115{
2116 const std::pair<unsigned int, size_type> local_index =
2118 return components[local_index.first](local_index.second);
2119}
2120
2121
2122
2123template <typename VectorType>
2126{
2127 return operator()(i);
2128}
2129
2130
2131
2132template <typename VectorType>
2135{
2136 return operator()(i);
2137}
2138
2139
2140
2141template <typename VectorType>
2142template <typename OtherNumber>
2143inline void
2145 const std::vector<size_type> &indices,
2146 std::vector<OtherNumber> &values) const
2147{
2148 for (size_type i = 0; i < indices.size(); ++i)
2149 values[i] = operator()(indices[i]);
2150}
2151
2152
2153
2154template <typename VectorType>
2155inline void
2158 ArrayView<value_type> &entries) const
2159{
2160 AssertDimension(indices.size(), entries.size());
2161 for (unsigned int i = 0; i < indices.size(); ++i)
2162 {
2163 AssertIndexRange(indices[i], size());
2164 entries[i] = (*this)[indices[i]];
2165 }
2166}
2167
2168
2169
2170template <typename VectorType>
2171template <typename ForwardIterator, typename OutputIterator>
2172inline void
2174 ForwardIterator indices_begin,
2175 const ForwardIterator indices_end,
2176 OutputIterator values_begin) const
2177{
2178 while (indices_begin != indices_end)
2179 {
2180 *values_begin = operator()(*indices_begin);
2181 ++indices_begin;
2182 ++values_begin;
2183 }
2184}
2185
2186#endif // DOXYGEN
2187
2189
2190#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
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:1820
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#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:494
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