Reference documentation for deal.II version 9.4.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
block_vector_base.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2004 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_block_vector_base_h
17#define dealii_block_vector_base_h
18
19
20#include <deal.II/base/config.h>
21
26
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
45namespace internal
46{
47 template <typename T>
48 using has_block_t = decltype(std::declval<T const>().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<T const>().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 <class BlockVectorType, bool Constness>
121 {
122 public:
127
134 typename std::conditional<Constness,
135 const typename BlockVectorType::value_type,
136 typename BlockVectorType::value_type>::type;
137
144 using iterator_category = std::random_access_iterator_tag;
145 using difference_type = std::ptrdiff_t;
146 using reference = typename BlockVectorType::reference;
148
149 using dereference_type = typename std::conditional<
150 Constness,
152 typename BlockVectorType::BlockType::reference>::type;
153
158 using BlockVector = typename std::
159 conditional<Constness, const BlockVectorType, BlockVectorType>::type;
160
170
179
180
185
186 private:
197
198 public:
202 Iterator &
204
211 operator*() const;
212
219
224 Iterator &
226
234
239 Iterator &
241
249
254 template <bool OtherConstness>
255 bool
257
262 template <bool OtherConstness>
263 bool
265
271 template <bool OtherConstness>
272 bool
274
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>
302
309
316
321 Iterator &
323
328 Iterator &
330
341 "Your program tried to compare iterators pointing to "
342 "different block vectors. There is no reasonable way "
343 "to do this.");
344
346 private:
353
358
363 unsigned int current_block;
365
374
378 void
380
384 void
386
387 // Mark all other instances of this template as friends.
388 template <typename, bool>
389 friend class Iterator;
390 };
391 } // namespace BlockVectorIterators
392} // namespace internal
393
394
432template <class VectorType>
434{
435public:
439 using BlockType = VectorType;
440
441 /*
442 * Declare standard types used in
443 * all containers. These types
444 * parallel those in the
445 * <tt>C++</tt> standard
446 * libraries
447 * <tt>std::vector<...></tt>
448 * class. This includes iterator
449 * types.
450 */
451 using value_type = typename BlockType::value_type;
453 using const_pointer = const value_type *;
454 using iterator =
458 using reference = typename BlockType::reference;
459 using const_reference = typename BlockType::const_reference;
461
471 using real_type = typename BlockType::real_type;
472
476 BlockVectorBase() = default;
477
481 BlockVectorBase(const BlockVectorBase & /*V*/) = default;
482
488 BlockVectorBase(BlockVectorBase && /*V*/) noexcept = default;
489
496 void
498
509 void
511
515 BlockType &
516 block(const unsigned int i);
517
521 const BlockType &
522 block(const unsigned int i) const;
523
529 const BlockIndices &
531
535 unsigned int
536 n_blocks() const;
537
542 std::size_t
543 size() const;
544
550 std::size_t
552
571
577
583 begin() const;
584
590
596 end() const;
597
602 operator()(const size_type i) const;
603
608 operator()(const size_type i);
609
616 operator[](const size_type i) const;
617
624 operator[](const size_type i);
625
641 template <typename OtherNumber>
642 void
643 extract_subvector_to(const std::vector<size_type> &indices,
644 std::vector<OtherNumber> & values) const;
645
673 template <typename ForwardIterator, typename OutputIterator>
674 void
675 extract_subvector_to(ForwardIterator indices_begin,
676 const ForwardIterator indices_end,
677 OutputIterator values_begin) const;
678
684 operator=(const value_type s);
685
690 operator=(const BlockVectorBase &V);
691
698 operator=(BlockVectorBase && /*V*/) = default; // NOLINT
699
703 template <class VectorType2>
705 operator=(const BlockVectorBase<VectorType2> &V);
706
711 operator=(const VectorType &v);
712
717 template <class VectorType2>
718 bool
719 operator==(const BlockVectorBase<VectorType2> &v) const;
720
725 operator*(const BlockVectorBase &V) const;
726
731 norm_sqr() const;
732
737 mean_value() const;
738
743 l1_norm() const;
744
750 l2_norm() const;
751
757 linfty_norm() const;
758
783 const BlockVectorBase &V,
784 const BlockVectorBase &W);
785
790 bool
792
798 bool
799 all_zero() const;
800
806 bool
808
813 operator+=(const BlockVectorBase &V);
814
819 operator-=(const BlockVectorBase &V);
820
821
826 template <typename Number>
827 void
828 add(const std::vector<size_type> &indices, const std::vector<Number> &values);
829
834 template <typename Number>
835 void
836 add(const std::vector<size_type> &indices, const Vector<Number> &values);
837
843 template <typename Number>
844 void
845 add(const size_type n_elements,
846 const size_type *indices,
847 const Number * values);
848
853 void
854 add(const value_type s);
855
859 void
860 add(const value_type a, const BlockVectorBase &V);
861
865 void
867 const BlockVectorBase &V,
868 const value_type b,
869 const BlockVectorBase &W);
870
874 void
876
880 void
881 sadd(const value_type s, const value_type a, const BlockVectorBase &V);
882
886 void
888 const value_type a,
889 const BlockVectorBase &V,
890 const value_type b,
891 const BlockVectorBase &W);
892
896 void
898 const value_type a,
899 const BlockVectorBase &V,
900 const value_type b,
901 const BlockVectorBase &W,
902 const value_type c,
903 const BlockVectorBase &X);
904
909 operator*=(const value_type factor);
910
915 operator/=(const value_type factor);
916
921 template <class BlockVector2>
922 void
923 scale(const BlockVector2 &v);
924
928 template <class BlockVector2>
929 void
930 equ(const value_type a, const BlockVector2 &V);
931
936 void
938
943 std::size_t
945
946protected:
950 std::vector<VectorType> components;
951
957
958 // Make the iterator class a friend.
959 template <typename N, bool C>
960 friend class ::internal::BlockVectorIterators::Iterator;
961
962 template <typename>
963 friend class BlockVectorBase;
964};
965
966
969/*----------------------- Inline functions ----------------------------------*/
970
971
972#ifndef DOXYGEN
973namespace internal
974{
975 namespace BlockVectorIterators
976 {
977 template <class BlockVectorType, bool Constness>
978 inline Iterator<BlockVectorType, Constness>::Iterator(
979 const Iterator<BlockVectorType, Constness> &c)
980 : parent(c.parent)
981 , global_index(c.global_index)
982 , current_block(c.current_block)
983 , index_within_block(c.index_within_block)
984 , next_break_forward(c.next_break_forward)
985 , next_break_backward(c.next_break_backward)
986 {}
987
988
989
990 template <class 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 // Only permit copy-constructing const iterators from non-const
1001 // iterators, and not vice versa (i.e., Constness must always be
1002 // true).
1003 static_assert(Constness == true,
1004 "Constructing a non-const iterator from a const iterator "
1005 "does not make sense.");
1006 }
1007
1008
1009
1010 template <class BlockVectorType, bool Constness>
1011 inline Iterator<BlockVectorType, Constness>::Iterator(
1012 BlockVector & parent,
1013 const size_type global_index,
1014 const size_type current_block,
1015 const size_type index_within_block,
1016 const size_type next_break_forward,
1017 const size_type next_break_backward)
1018 : parent(&parent)
1020 , current_block(current_block)
1021 , index_within_block(index_within_block)
1022 , next_break_forward(next_break_forward)
1023 , next_break_backward(next_break_backward)
1024 {}
1025
1026
1027
1028 template <class BlockVectorType, bool Constness>
1029 inline Iterator<BlockVectorType, Constness> &
1030 Iterator<BlockVectorType, Constness>::operator=(const Iterator &c)
1031 {
1032 parent = c.parent;
1033 global_index = c.global_index;
1034 index_within_block = c.index_within_block;
1035 current_block = c.current_block;
1036 next_break_forward = c.next_break_forward;
1037 next_break_backward = c.next_break_backward;
1038
1039 return *this;
1040 }
1041
1042
1043
1044 template <class BlockVectorType, bool Constness>
1045 inline typename Iterator<BlockVectorType, Constness>::dereference_type
1047 {
1048 return parent->block(current_block)(index_within_block);
1049 }
1050
1051
1052
1053 template <class BlockVectorType, bool Constness>
1054 inline typename Iterator<BlockVectorType, Constness>::dereference_type
1055 Iterator<BlockVectorType, Constness>::operator[](
1056 const difference_type d) const
1057 {
1058 // if the index pointed to is
1059 // still within the block we
1060 // currently point into, then we
1061 // can save the computation of
1062 // the block
1063 if ((global_index + d >= next_break_backward) &&
1064 (global_index + d <= next_break_forward))
1065 return parent->block(current_block)(index_within_block + d);
1066
1067 // if the index is not within the
1068 // block of the block vector into
1069 // which we presently point, then
1070 // there is no way: we have to
1071 // search for the block. this can
1072 // be done through the parent
1073 // class as well.
1074 return (*parent)(global_index + d);
1075 }
1076
1077
1078
1079 template <class BlockVectorType, bool Constness>
1080 inline Iterator<BlockVectorType, Constness> &
1082 {
1083 move_forward();
1084 return *this;
1085 }
1086
1087
1088
1089 template <class BlockVectorType, bool Constness>
1090 inline Iterator<BlockVectorType, Constness>
1092 {
1093 const Iterator old_value = *this;
1094 move_forward();
1095 return old_value;
1096 }
1097
1098
1099
1100 template <class BlockVectorType, bool Constness>
1101 inline Iterator<BlockVectorType, Constness> &
1102 Iterator<BlockVectorType, Constness>::operator--()
1103 {
1104 move_backward();
1105 return *this;
1106 }
1107
1108
1109
1110 template <class BlockVectorType, bool Constness>
1111 inline Iterator<BlockVectorType, Constness>
1112 Iterator<BlockVectorType, Constness>::operator--(int)
1113 {
1114 const Iterator old_value = *this;
1115 move_backward();
1116 return old_value;
1117 }
1118
1119
1120
1121 template <class BlockVectorType, bool Constness>
1122 template <bool OtherConstness>
1123 inline bool
1125 const Iterator<BlockVectorType, OtherConstness> &i) const
1126 {
1127 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1128
1129 return (global_index == i.global_index);
1130 }
1131
1132
1133
1134 template <class BlockVectorType, bool Constness>
1135 template <bool OtherConstness>
1136 inline bool
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 <class BlockVectorType, bool Constness>
1148 template <bool OtherConstness>
1149 inline bool
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 <class BlockVectorType, bool Constness>
1161 template <bool OtherConstness>
1162 inline bool
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 <class BlockVectorType, bool Constness>
1174 template <bool OtherConstness>
1175 inline bool
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 <class BlockVectorType, bool Constness>
1187 template <bool OtherConstness>
1188 inline bool
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 <class BlockVectorType, bool Constness>
1200 template <bool OtherConstness>
1201 inline typename Iterator<BlockVectorType, Constness>::difference_type
1203 const Iterator<BlockVectorType, OtherConstness> &i) const
1204 {
1205 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1206
1207 return (static_cast<signed int>(global_index) -
1208 static_cast<signed int>(i.global_index));
1209 }
1210
1211
1212
1213 template <class BlockVectorType, bool Constness>
1214 inline Iterator<BlockVectorType, Constness>
1216 const difference_type &d) const
1217 {
1218 // if the index pointed to is
1219 // still within the block we
1220 // currently point into, then we
1221 // can save the computation of
1222 // the block
1223 if ((global_index + d >= next_break_backward) &&
1224 (global_index + d <= next_break_forward))
1225 return Iterator(*parent,
1226 global_index + d,
1227 current_block,
1228 index_within_block + d,
1229 next_break_forward,
1230 next_break_backward);
1231 else
1232 // outside present block, so
1233 // have to seek new block
1234 // anyway
1235 return Iterator(*parent, global_index + d);
1236 }
1237
1238
1239
1240 template <class BlockVectorType, bool Constness>
1241 inline Iterator<BlockVectorType, Constness>
1243 const difference_type &d) const
1244 {
1245 // if the index pointed to is
1246 // still within the block we
1247 // currently point into, then we
1248 // can save the computation of
1249 // the block
1250 if ((global_index - d >= next_break_backward) &&
1251 (global_index - d <= next_break_forward))
1252 return Iterator(*parent,
1253 global_index - d,
1254 current_block,
1255 index_within_block - d,
1256 next_break_forward,
1257 next_break_backward);
1258 else
1259 // outside present block, so
1260 // have to seek new block
1261 // anyway
1262 return Iterator(*parent, global_index - d);
1263 }
1264
1265
1266
1267 template <class BlockVectorType, bool Constness>
1268 inline Iterator<BlockVectorType, Constness> &
1269 Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1270 {
1271 // if the index pointed to is
1272 // still within the block we
1273 // currently point into, then we
1274 // can save the computation of
1275 // the block
1276 if ((global_index + d >= next_break_backward) &&
1277 (global_index + d <= next_break_forward))
1278 {
1279 global_index += d;
1280 index_within_block += d;
1281 }
1282 else
1283 // outside present block, so
1284 // have to seek new block
1285 // anyway
1286 *this = Iterator(*parent, global_index + d);
1287
1288 return *this;
1289 }
1290
1291
1292
1293 template <class BlockVectorType, bool Constness>
1294 inline Iterator<BlockVectorType, Constness> &
1295 Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1296 {
1297 // if the index pointed to is
1298 // still within the block we
1299 // currently point into, then we
1300 // can save the computation of
1301 // the block
1302 if ((global_index - d >= next_break_backward) &&
1303 (global_index - d <= next_break_forward))
1304 {
1305 global_index -= d;
1306 index_within_block -= d;
1307 }
1308 else
1309 // outside present block, so
1310 // have to seek new block
1311 // anyway
1312 *this = Iterator(*parent, global_index - d);
1313
1314 return *this;
1315 }
1316
1317
1318 template <class BlockVectorType, bool Constness>
1319 Iterator<BlockVectorType, Constness>::Iterator(BlockVector & parent,
1320 const size_type global_index)
1321 : parent(&parent)
1323 {
1324 // find which block we are
1325 // in. for this, take into
1326 // account that it happens at
1327 // times that people want to
1328 // initialize iterators
1329 // past-the-end
1330 if (global_index < parent.size())
1331 {
1332 const std::pair<size_type, size_type> indices =
1334 current_block = indices.first;
1335 index_within_block = indices.second;
1336
1337 next_break_backward =
1338 parent.block_indices.local_to_global(current_block, 0);
1339 next_break_forward =
1340 (parent.block_indices.local_to_global(current_block, 0) +
1341 parent.block_indices.block_size(current_block) - 1);
1342 }
1343 else
1344 // past the end. only have one
1345 // value for this
1346 {
1347 this->global_index = parent.size();
1348 current_block = parent.n_blocks();
1349 index_within_block = 0;
1350 next_break_backward = global_index;
1351 next_break_forward = numbers::invalid_size_type;
1352 };
1353 }
1354
1355
1356
1357 template <class BlockVectorType, bool Constness>
1358 void
1359 Iterator<BlockVectorType, Constness>::move_forward()
1360 {
1361 if (global_index != next_break_forward)
1362 ++index_within_block;
1363 else
1364 {
1365 // ok, we traverse a boundary
1366 // between blocks:
1367 index_within_block = 0;
1368 ++current_block;
1369
1370 // break backwards is now old
1371 // break forward
1372 next_break_backward = next_break_forward + 1;
1373
1374 // compute new break forward
1375 if (current_block < parent->block_indices.size())
1376 next_break_forward +=
1377 parent->block_indices.block_size(current_block);
1378 else
1379 // if we are beyond the end,
1380 // then move the next
1381 // boundary arbitrarily far
1382 // away
1383 next_break_forward = numbers::invalid_size_type;
1384 };
1385
1386 ++global_index;
1387 }
1388
1389
1390
1391 template <class BlockVectorType, bool Constness>
1392 void
1393 Iterator<BlockVectorType, Constness>::move_backward()
1394 {
1395 if (global_index != next_break_backward)
1396 --index_within_block;
1397 else if (current_block != 0)
1398 {
1399 // ok, we traverse a boundary
1400 // between blocks:
1401 --current_block;
1402 index_within_block =
1403 parent->block_indices.block_size(current_block) - 1;
1404
1405 // break forwards is now old
1406 // break backward
1407 next_break_forward = next_break_backward - 1;
1408
1409 // compute new break forward
1410 next_break_backward -=
1411 parent->block_indices.block_size(current_block);
1412 }
1413 else
1414 // current block was 0, we now
1415 // get into unspecified terrain
1416 {
1417 --current_block;
1418 index_within_block = numbers::invalid_size_type;
1419 next_break_forward = 0;
1420 next_break_backward = 0;
1421 };
1422
1423 --global_index;
1424 }
1425
1426
1427 } // namespace BlockVectorIterators
1428
1429} // namespace internal
1430
1431
1432
1433template <class VectorType>
1434inline std::size_t
1436{
1437 return block_indices.total_size();
1438}
1439
1440
1441
1442template <class VectorType>
1443inline std::size_t
1445{
1446 std::size_t local_size = 0;
1447 for (unsigned int b = 0; b < n_blocks(); ++b)
1448 local_size += block(b).locally_owned_size();
1449 return local_size;
1450}
1451
1452
1453
1454template <class VectorType>
1455inline IndexSet
1457{
1458 IndexSet is(size());
1459
1460 // copy index sets from blocks into the global one, shifted
1461 // by the appropriate amount for each block
1462 for (unsigned int b = 0; b < n_blocks(); ++b)
1463 {
1464 IndexSet x = block(b).locally_owned_elements();
1465 is.add_indices(x, block_indices.block_start(b));
1466 }
1467
1468 is.compress();
1469
1470 return is;
1471}
1472
1473
1474
1475template <class VectorType>
1476inline unsigned int
1478{
1479 return block_indices.size();
1480}
1481
1482
1483template <class VectorType>
1485BlockVectorBase<VectorType>::block(const unsigned int i)
1486{
1488
1489 return components[i];
1490}
1491
1492
1493
1494template <class VectorType>
1495inline const typename BlockVectorBase<VectorType>::BlockType &
1496BlockVectorBase<VectorType>::block(const unsigned int i) const
1497{
1499
1500 return components[i];
1501}
1502
1503
1504
1505template <class VectorType>
1506inline const BlockIndices &
1508{
1509 return block_indices;
1510}
1511
1512
1513template <class VectorType>
1514inline void
1516{
1517 std::vector<size_type> sizes(n_blocks());
1518
1519 for (size_type i = 0; i < n_blocks(); ++i)
1520 sizes[i] = block(i).size();
1521
1522 block_indices.reinit(sizes);
1523}
1524
1525
1526
1527template <class VectorType>
1528inline void
1530 ::VectorOperation::values operation)
1531{
1532 for (unsigned int i = 0; i < n_blocks(); ++i)
1533 block(i).compress(operation);
1534}
1535
1536
1537
1538template <class VectorType>
1541{
1542 return iterator(*this, 0U);
1543}
1544
1545
1546
1547template <class VectorType>
1550{
1551 return const_iterator(*this, 0U);
1552}
1553
1554
1555template <class VectorType>
1558{
1559 return iterator(*this, size());
1560}
1561
1562
1563
1564template <class VectorType>
1567{
1568 return const_iterator(*this, size());
1569}
1570
1571
1572template <class VectorType>
1573inline bool
1575{
1576 const std::pair<size_type, size_type> local_index =
1578
1579 return components[local_index.first].in_local_range(global_index);
1580}
1581
1582
1583template <class VectorType>
1584bool
1586{
1587 for (size_type i = 0; i < n_blocks(); ++i)
1588 if (components[i].all_zero() == false)
1589 return false;
1590
1591 return true;
1592}
1593
1594
1595
1596template <class VectorType>
1597bool
1599{
1600 for (size_type i = 0; i < n_blocks(); ++i)
1601 if (components[i].is_non_negative() == false)
1602 return false;
1603
1604 return true;
1605}
1606
1607
1608
1609template <class VectorType>
1612 const BlockVectorBase<VectorType> &v) const
1613{
1614 Assert(n_blocks() == v.n_blocks(),
1616
1617 value_type sum = 0.;
1618 for (size_type i = 0; i < n_blocks(); ++i)
1619 sum += components[i] * v.components[i];
1620
1621 return sum;
1622}
1623
1624
1625template <class VectorType>
1628{
1629 real_type sum = 0.;
1630 for (size_type i = 0; i < n_blocks(); ++i)
1631 sum += components[i].norm_sqr();
1632
1633 return sum;
1634}
1635
1636
1637
1638template <class VectorType>
1641{
1642 value_type sum = 0.;
1643 // need to do static_cast as otherwise it won't work with
1644 // value_type=complex<T>
1645 for (size_type i = 0; i < n_blocks(); ++i)
1646 sum += components[i].mean_value() *
1648 components[i].size()));
1649
1651}
1652
1653
1654
1655template <class VectorType>
1658{
1659 real_type sum = 0.;
1660 for (size_type i = 0; i < n_blocks(); ++i)
1661 sum += components[i].l1_norm();
1662
1663 return sum;
1664}
1665
1666
1667
1668template <class VectorType>
1671{
1672 return std::sqrt(norm_sqr());
1673}
1674
1675
1676
1677template <class VectorType>
1680{
1681 real_type sum = 0.;
1682 for (size_type i = 0; i < n_blocks(); ++i)
1683 {
1684 value_type newval = components[i].linfty_norm();
1685 if (sum < newval)
1686 sum = newval;
1687 }
1688 return sum;
1689}
1690
1691
1692
1693template <class VectorType>
1699{
1700 AssertDimension(n_blocks(), V.n_blocks());
1702
1703 value_type sum = 0.;
1704 for (size_type i = 0; i < n_blocks(); ++i)
1705 sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1706
1707 return sum;
1708}
1709
1710
1711
1712template <class VectorType>
1715{
1716 Assert(n_blocks() == v.n_blocks(),
1718
1719 for (size_type i = 0; i < n_blocks(); ++i)
1720 {
1721 components[i] += v.components[i];
1722 }
1723
1724 return *this;
1725}
1726
1727
1728
1729template <class VectorType>
1732{
1733 Assert(n_blocks() == v.n_blocks(),
1735
1736 for (size_type i = 0; i < n_blocks(); ++i)
1737 {
1738 components[i] -= v.components[i];
1739 }
1740 return *this;
1741}
1742
1743
1744
1745template <class VectorType>
1746template <typename Number>
1747inline void
1748BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1749 const std::vector<Number> & values)
1750{
1751 Assert(indices.size() == values.size(),
1752 ExcDimensionMismatch(indices.size(), values.size()));
1753 add(indices.size(), indices.data(), values.data());
1754}
1755
1756
1757
1758template <class VectorType>
1759template <typename Number>
1760inline void
1761BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1762 const Vector<Number> & values)
1763{
1764 Assert(indices.size() == values.size(),
1765 ExcDimensionMismatch(indices.size(), values.size()));
1766 const size_type n_indices = indices.size();
1767 for (size_type i = 0; i < n_indices; ++i)
1768 (*this)(indices[i]) += values(i);
1769}
1770
1771
1772
1773template <class VectorType>
1774template <typename Number>
1775inline void
1777 const size_type *indices,
1778 const Number * values)
1779{
1780 for (size_type i = 0; i < n_indices; ++i)
1781 (*this)(indices[i]) += values[i];
1782}
1783
1784
1785
1786template <class VectorType>
1787void
1789{
1790 AssertIsFinite(a);
1791
1792 for (size_type i = 0; i < n_blocks(); ++i)
1793 {
1794 components[i].add(a);
1795 }
1796}
1797
1798
1799
1800template <class VectorType>
1801void
1804{
1805 AssertIsFinite(a);
1806
1807 Assert(n_blocks() == v.n_blocks(),
1809
1810 for (size_type i = 0; i < n_blocks(); ++i)
1811 {
1812 components[i].add(a, v.components[i]);
1813 }
1814}
1815
1816
1817
1818template <class VectorType>
1819void
1822 const value_type b,
1824{
1825 AssertIsFinite(a);
1827
1828 Assert(n_blocks() == v.n_blocks(),
1830 Assert(n_blocks() == w.n_blocks(),
1831 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1832
1833
1834 for (size_type i = 0; i < n_blocks(); ++i)
1835 {
1836 components[i].add(a, v.components[i], b, w.components[i]);
1837 }
1838}
1839
1840
1841
1842template <class VectorType>
1843void
1846{
1847 AssertIsFinite(x);
1848
1849 Assert(n_blocks() == v.n_blocks(),
1851
1852 for (size_type i = 0; i < n_blocks(); ++i)
1853 {
1854 components[i].sadd(x, v.components[i]);
1855 }
1856}
1857
1858
1859
1860template <class VectorType>
1861void
1863 const value_type a,
1865{
1866 AssertIsFinite(x);
1867 AssertIsFinite(a);
1868
1869 Assert(n_blocks() == v.n_blocks(),
1871
1872 for (size_type i = 0; i < n_blocks(); ++i)
1873 {
1874 components[i].sadd(x, a, v.components[i]);
1875 }
1876}
1877
1878
1879
1880template <class VectorType>
1881void
1883 const value_type a,
1885 const value_type b,
1887{
1888 AssertIsFinite(x);
1889 AssertIsFinite(a);
1891
1892 Assert(n_blocks() == v.n_blocks(),
1894 Assert(n_blocks() == w.n_blocks(),
1895 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1896
1897 for (size_type i = 0; i < n_blocks(); ++i)
1898 {
1899 components[i].sadd(x, a, v.components[i], b, w.components[i]);
1900 }
1901}
1902
1903
1904
1905template <class VectorType>
1906void
1908 const value_type a,
1910 const value_type b,
1912 const value_type c,
1914{
1915 AssertIsFinite(x);
1916 AssertIsFinite(a);
1918 AssertIsFinite(c);
1919
1920 Assert(n_blocks() == v.n_blocks(),
1922 Assert(n_blocks() == w.n_blocks(),
1923 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1924 Assert(n_blocks() == y.n_blocks(),
1926
1927 for (size_type i = 0; i < n_blocks(); ++i)
1928 {
1929 components[i].sadd(
1930 x, a, v.components[i], b, w.components[i], c, y.components[i]);
1931 }
1932}
1933
1934
1935
1936template <class VectorType>
1937template <class BlockVector2>
1938void
1939BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
1940{
1941 Assert(n_blocks() == v.n_blocks(),
1942 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1943 for (size_type i = 0; i < n_blocks(); ++i)
1944 components[i].scale(v.block(i));
1945}
1946
1947
1948
1949template <class VectorType>
1950std::size_t
1952{
1953 return (MemoryConsumption::memory_consumption(this->block_indices) +
1954 MemoryConsumption::memory_consumption(this->components));
1955}
1956
1957
1958
1959template <class VectorType>
1960template <class BlockVector2>
1961void
1962BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
1963{
1964 AssertIsFinite(a);
1965
1966 Assert(n_blocks() == v.n_blocks(),
1967 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1968
1969 for (size_type i = 0; i < n_blocks(); ++i)
1970 components[i].equ(a, v.components[i]);
1971}
1972
1973
1974
1975template <class VectorType>
1976void
1978{
1979 for (size_type i = 0; i < n_blocks(); ++i)
1980 block(i).update_ghost_values();
1981}
1982
1983
1984
1985template <class VectorType>
1988{
1989 AssertIsFinite(s);
1990
1991 for (size_type i = 0; i < n_blocks(); ++i)
1992 components[i] = s;
1993
1994 return *this;
1995}
1996
1997
1998template <class VectorType>
2001{
2003
2004 for (size_type i = 0; i < n_blocks(); ++i)
2005 components[i] = v.components[i];
2006
2007 return *this;
2008}
2009
2010
2011template <class VectorType>
2012template <class VectorType2>
2015{
2017
2018 for (size_type i = 0; i < n_blocks(); ++i)
2019 components[i] = v.components[i];
2020
2021 return *this;
2022}
2023
2024
2025
2026template <class VectorType>
2028BlockVectorBase<VectorType>::operator=(const VectorType &v)
2029{
2030 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2031
2032 size_type index_v = 0;
2033 for (size_type b = 0; b < n_blocks(); ++b)
2034 for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2035 block(b)(i) = v(index_v);
2036
2037 return *this;
2038}
2039
2040
2041
2042template <class VectorType>
2043template <class VectorType2>
2044inline bool
2046 const BlockVectorBase<VectorType2> &v) const
2047{
2049
2050 for (size_type i = 0; i < n_blocks(); ++i)
2051 if (!(components[i] == v.components[i]))
2052 return false;
2053
2054 return true;
2055}
2056
2057
2058
2059template <class VectorType>
2062{
2063 AssertIsFinite(factor);
2064
2065 for (size_type i = 0; i < n_blocks(); ++i)
2066 components[i] *= factor;
2067
2068 return *this;
2069}
2070
2071
2072
2073template <class VectorType>
2076{
2077 AssertIsFinite(factor);
2078 Assert(factor != 0., ExcDivideByZero());
2079
2080 const value_type inverse_factor = value_type(1.) / factor;
2081
2082 for (size_type i = 0; i < n_blocks(); ++i)
2083 components[i] *= inverse_factor;
2084
2085 return *this;
2086}
2087
2088
2089template <class VectorType>
2092{
2093 const std::pair<unsigned int, size_type> local_index =
2095 return components[local_index.first](local_index.second);
2096}
2097
2098
2099
2100template <class VectorType>
2103{
2104 const std::pair<unsigned int, size_type> local_index =
2106 return components[local_index.first](local_index.second);
2107}
2108
2109
2110
2111template <class VectorType>
2114{
2115 return operator()(i);
2116}
2117
2118
2119
2120template <class VectorType>
2123{
2124 return operator()(i);
2125}
2126
2127
2128
2129template <typename VectorType>
2130template <typename OtherNumber>
2131inline void
2133 const std::vector<size_type> &indices,
2134 std::vector<OtherNumber> & values) const
2135{
2136 for (size_type i = 0; i < indices.size(); ++i)
2137 values[i] = operator()(indices[i]);
2138}
2139
2140
2141
2142template <typename VectorType>
2143template <typename ForwardIterator, typename OutputIterator>
2144inline void
2146 ForwardIterator indices_begin,
2147 const ForwardIterator indices_end,
2148 OutputIterator values_begin) const
2149{
2150 while (indices_begin != indices_end)
2151 {
2152 *values_begin = operator()(*indices_begin);
2153 indices_begin++;
2154 values_begin++;
2155 }
2156}
2157
2158#endif // DOXYGEN
2159
2161
2162#endif
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
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
Definition: vector.h:109
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typenameProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
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
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
bool all_zero() const
BlockVectorBase & operator=(const value_type s)
std::size_t size() const
void collect_sizes()
#define Assert(cond, exc)
Definition: exceptions.h:1473
BlockVectorBase & operator*=(const value_type factor)
void sadd(const value_type s, const BlockVectorBase &V)
std::vector< VectorType > components
static ::ExceptionBase & ExcPointerToDifferentVectors()
iterator begin()
Expression operator>=(const Expression &lhs, const Expression &rhs)
bool in_local_range(const size_type global_index) const
#define AssertIsFinite(number)
Definition: exceptions.h:1758
::internal::BlockVectorIterators::Iterator< BlockVectorBase, true > const_iterator
BlockVectorBase(BlockVectorBase &&) noexcept=default
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
BlockVectorBase(const BlockVectorBase &)=default
typename BlockType::reference reference
iterator end()
value_type operator[](const size_type i) const
BlockIndices block_indices
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
typename BlockType::real_type real_type
Expression operator>(const Expression &lhs, const Expression &rhs)
void scale(const BlockVector2 &v)
void compress(::VectorOperation::values operation)
BlockType & block(const unsigned int i)
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcDifferentBlockIndices()
Expression operator<=(const Expression &lhs, const Expression &rhs)
real_type l2_norm() const
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
value_type operator()(const size_type i) const
IndexSet locally_owned_elements() const
bool is_non_negative() const
BlockVectorBase & operator-=(const BlockVectorBase &V)
bool operator==(const BlockVectorBase< VectorType2 > &v) const
std::size_t locally_owned_size() const
const BlockIndices & get_block_indices() const
value_type * pointer
void equ(const value_type a, const BlockVector2 &V)
decltype(std::declval< T const >().block(0)) has_block_t
constexpr bool has_n_blocks
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)
Iterator(const Iterator< BlockVectorType, !Constness > &c)
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
typename std::conditional< Constness, value_type, typename BlockVectorType::BlockType::reference >::type dereference_type
Iterator(BlockVector &parent, const size_type global_index)
bool operator>=(const Iterator< BlockVectorType, OtherConstness > &i) const
decltype(std::declval< T const >().n_blocks()) has_n_blocks_t
typename std::conditional< Constness, const typename BlockVectorType::value_type, typename BlockVectorType::value_type >::type value_type
Iterator & operator-=(const difference_type &d)
constexpr bool is_block_vector
bool operator<=(const Iterator< BlockVectorType, OtherConstness > &i) const
typename BlockVectorType::reference reference
static const bool value
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
constexpr bool has_block
static const char U
static const char N
static const char V
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
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)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
T sum(const T &t, const MPI_Comm &mpi_communicator)
const types::global_dof_index invalid_size_type
Definition: types.h:210
STL namespace.
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition: types.h:76
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
BarycentricPolynomial< dim, Number1 > operator+(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)