Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
block_vector_base.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2004 - 2023 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
46namespace internal
47{
48 template <typename T>
49 using has_block_t = decltype(std::declval<T const>().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<T const>().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 <class BlockVectorType, bool Constness>
122 {
123 public:
128
135 typename std::conditional<Constness,
136 const typename BlockVectorType::value_type,
137 typename BlockVectorType::value_type>::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
157 using dereference_type = typename std::conditional<
158 Constness,
160 typename BlockVectorType::BlockType::reference>::type;
161
166 using BlockVector = typename std::
167 conditional<Constness, const BlockVectorType, BlockVectorType>::type;
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 <class VectorType>
442{
443public:
447 using BlockType = VectorType;
448
449 /*
450 * Declare standard types used in
451 * all containers. These types
452 * parallel those in the
453 * <tt>C++</tt> standard
454 * libraries
455 * <tt>std::vector<...></tt>
456 * class. This includes iterator
457 * types.
458 */
459 using value_type = typename BlockType::value_type;
461 using const_pointer = const value_type *;
462 using iterator =
466 using reference = typename BlockType::reference;
467 using const_reference = typename BlockType::const_reference;
469
479 using real_type = typename BlockType::real_type;
480
484 BlockVectorBase() = default;
485
489 BlockVectorBase(const BlockVectorBase & /*V*/) = default;
490
496 BlockVectorBase(BlockVectorBase && /*V*/) noexcept = default;
497
504 void
506
517 void
518 compress(VectorOperation::values operation);
519
523 BlockType &
524 block(const unsigned int i);
525
529 const BlockType &
530 block(const unsigned int i) const;
531
537 const BlockIndices &
539
543 unsigned int
544 n_blocks() const;
545
550 std::size_t
551 size() const;
552
558 std::size_t
560
579
585
591 begin() const;
592
598
604 end() const;
605
610 operator()(const size_type i) const;
611
616 operator()(const size_type i);
617
624 operator[](const size_type i) const;
625
632 operator[](const size_type i);
633
649 template <typename OtherNumber>
650 void
651 extract_subvector_to(const std::vector<size_type> &indices,
652 std::vector<OtherNumber> & values) const;
653
681 template <typename ForwardIterator, typename OutputIterator>
682 void
683 extract_subvector_to(ForwardIterator indices_begin,
684 const ForwardIterator indices_end,
685 OutputIterator values_begin) const;
686
692 operator=(const value_type s);
693
698 operator=(const BlockVectorBase &V);
699
706 operator=(BlockVectorBase && /*V*/) = default; // NOLINT
707
711 template <class VectorType2>
713 operator=(const BlockVectorBase<VectorType2> &V);
714
719 operator=(const VectorType &v);
720
725 template <class VectorType2>
726 bool
727 operator==(const BlockVectorBase<VectorType2> &v) const;
728
733 operator*(const BlockVectorBase &V) const;
734
739 norm_sqr() const;
740
745 mean_value() const;
746
751 l1_norm() const;
752
758 l2_norm() const;
759
765 linfty_norm() const;
766
791 const BlockVectorBase &V,
792 const BlockVectorBase &W);
793
798 bool
799 in_local_range(const size_type global_index) const;
800
806 bool
807 all_zero() const;
808
814 bool
816
821 operator+=(const BlockVectorBase &V);
822
827 operator-=(const BlockVectorBase &V);
828
829
834 template <typename Number>
835 void
836 add(const std::vector<size_type> &indices, const std::vector<Number> &values);
837
842 template <typename Number>
843 void
844 add(const std::vector<size_type> &indices, const Vector<Number> &values);
845
851 template <typename Number>
852 void
853 add(const size_type n_elements,
854 const size_type *indices,
855 const Number * values);
856
861 void
862 add(const value_type s);
863
867 void
868 add(const value_type a, const BlockVectorBase &V);
869
873 void
875 const BlockVectorBase &V,
876 const value_type b,
877 const BlockVectorBase &W);
878
882 void
883 sadd(const value_type s, const BlockVectorBase &V);
884
888 void
889 sadd(const value_type s, const value_type a, const BlockVectorBase &V);
890
894 void
896 const value_type a,
897 const BlockVectorBase &V,
898 const value_type b,
899 const BlockVectorBase &W);
900
904 void
906 const value_type a,
907 const BlockVectorBase &V,
908 const value_type b,
909 const BlockVectorBase &W,
910 const value_type c,
911 const BlockVectorBase &X);
912
917 operator*=(const value_type factor);
918
923 operator/=(const value_type factor);
924
929 template <class BlockVector2>
930 void
931 scale(const BlockVector2 &v);
932
936 template <class BlockVector2>
937 void
938 equ(const value_type a, const BlockVector2 &V);
939
944 void
946
951 std::size_t
953
954protected:
958 std::vector<VectorType> components;
959
965
966 // Make the iterator class a friend.
967 template <typename N, bool C>
968 friend class ::internal::BlockVectorIterators::Iterator;
969
970 template <typename>
971 friend class BlockVectorBase;
972};
973
974
977/*----------------------- Inline functions ----------------------------------*/
978
979
980#ifndef DOXYGEN
981namespace internal
982{
983 namespace BlockVectorIterators
984 {
985 template <class BlockVectorType, bool Constness>
986 inline Iterator<BlockVectorType, Constness>::Iterator(
987 const Iterator<BlockVectorType, Constness> &c)
988 : parent(c.parent)
989 , global_index(c.global_index)
990 , current_block(c.current_block)
991 , index_within_block(c.index_within_block)
992 , next_break_forward(c.next_break_forward)
993 , next_break_backward(c.next_break_backward)
994 {}
995
996
997
998 template <class BlockVectorType, bool Constness>
999 inline Iterator<BlockVectorType, Constness>::Iterator(
1000 const Iterator<BlockVectorType, !Constness> &c)
1001 : parent(c.parent)
1002 , global_index(c.global_index)
1003 , current_block(c.current_block)
1004 , index_within_block(c.index_within_block)
1005 , next_break_forward(c.next_break_forward)
1006 , next_break_backward(c.next_break_backward)
1007 {
1008 // Only permit copy-constructing const iterators from non-const
1009 // iterators, and not vice versa (i.e., Constness must always be
1010 // true).
1011 static_assert(Constness == true,
1012 "Constructing a non-const iterator from a const iterator "
1013 "does not make sense.");
1014 }
1015
1016
1017
1018 template <class BlockVectorType, bool Constness>
1019 inline Iterator<BlockVectorType, Constness>::Iterator(
1020 BlockVector & parent,
1021 const size_type global_index,
1022 const size_type current_block,
1023 const size_type index_within_block,
1024 const size_type next_break_forward,
1025 const size_type next_break_backward)
1026 : parent(&parent)
1027 , global_index(global_index)
1028 , current_block(current_block)
1029 , index_within_block(index_within_block)
1030 , next_break_forward(next_break_forward)
1031 , next_break_backward(next_break_backward)
1032 {}
1033
1034
1035
1036 template <class BlockVectorType, bool Constness>
1037 inline Iterator<BlockVectorType, Constness> &
1038 Iterator<BlockVectorType, Constness>::operator=(const Iterator &c)
1039 {
1040 parent = c.parent;
1041 global_index = c.global_index;
1042 index_within_block = c.index_within_block;
1043 current_block = c.current_block;
1044 next_break_forward = c.next_break_forward;
1045 next_break_backward = c.next_break_backward;
1046
1047 return *this;
1048 }
1049
1050
1051
1052 template <class BlockVectorType, bool Constness>
1053 inline typename Iterator<BlockVectorType, Constness>::dereference_type
1054 Iterator<BlockVectorType, Constness>::operator*() const
1055 {
1056 return parent->block(current_block)(index_within_block);
1057 }
1058
1059
1060
1061 template <class BlockVectorType, bool Constness>
1062 inline typename Iterator<BlockVectorType, Constness>::dereference_type
1063 Iterator<BlockVectorType, Constness>::operator[](
1064 const difference_type d) const
1065 {
1066 // if the index pointed to is
1067 // still within the block we
1068 // currently point into, then we
1069 // can save the computation of
1070 // the block
1071 if ((global_index + d >= next_break_backward) &&
1072 (global_index + d <= next_break_forward))
1073 return parent->block(current_block)(index_within_block + d);
1074
1075 // if the index is not within the
1076 // block of the block vector into
1077 // which we presently point, then
1078 // there is no way: we have to
1079 // search for the block. this can
1080 // be done through the parent
1081 // class as well.
1082 return (*parent)(global_index + d);
1083 }
1084
1085
1086
1087 template <class BlockVectorType, bool Constness>
1088 inline Iterator<BlockVectorType, Constness> &
1089 Iterator<BlockVectorType, Constness>::operator++()
1090 {
1091 move_forward();
1092 return *this;
1093 }
1094
1095
1096
1097 template <class BlockVectorType, bool Constness>
1098 inline Iterator<BlockVectorType, Constness>
1099 Iterator<BlockVectorType, Constness>::operator++(int)
1100 {
1101 const Iterator old_value = *this;
1102 move_forward();
1103 return old_value;
1104 }
1105
1106
1107
1108 template <class BlockVectorType, bool Constness>
1109 inline Iterator<BlockVectorType, Constness> &
1110 Iterator<BlockVectorType, Constness>::operator--()
1111 {
1112 move_backward();
1113 return *this;
1114 }
1115
1116
1117
1118 template <class BlockVectorType, bool Constness>
1119 inline Iterator<BlockVectorType, Constness>
1120 Iterator<BlockVectorType, Constness>::operator--(int)
1121 {
1122 const Iterator old_value = *this;
1123 move_backward();
1124 return old_value;
1125 }
1126
1127
1128
1129 template <class BlockVectorType, bool Constness>
1130 template <bool OtherConstness>
1131 inline bool
1132 Iterator<BlockVectorType, Constness>::operator==(
1133 const Iterator<BlockVectorType, OtherConstness> &i) const
1134 {
1135 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1136
1137 return (global_index == i.global_index);
1138 }
1139
1140
1141
1142 template <class BlockVectorType, bool Constness>
1143 template <bool OtherConstness>
1144 inline bool
1145 Iterator<BlockVectorType, Constness>::operator!=(
1146 const Iterator<BlockVectorType, OtherConstness> &i) const
1147 {
1148 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1149
1150 return (global_index != i.global_index);
1151 }
1152
1153
1154
1155 template <class BlockVectorType, bool Constness>
1156 template <bool OtherConstness>
1157 inline bool
1158 Iterator<BlockVectorType, Constness>::operator<(
1159 const Iterator<BlockVectorType, OtherConstness> &i) const
1160 {
1161 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1162
1163 return (global_index < i.global_index);
1164 }
1165
1166
1167
1168 template <class BlockVectorType, bool Constness>
1169 template <bool OtherConstness>
1170 inline bool
1171 Iterator<BlockVectorType, Constness>::operator<=(
1172 const Iterator<BlockVectorType, OtherConstness> &i) const
1173 {
1174 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1175
1176 return (global_index <= i.global_index);
1177 }
1178
1179
1180
1181 template <class BlockVectorType, bool Constness>
1182 template <bool OtherConstness>
1183 inline bool
1184 Iterator<BlockVectorType, Constness>::operator>(
1185 const Iterator<BlockVectorType, OtherConstness> &i) const
1186 {
1187 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1188
1189 return (global_index > i.global_index);
1190 }
1191
1192
1193
1194 template <class BlockVectorType, bool Constness>
1195 template <bool OtherConstness>
1196 inline bool
1197 Iterator<BlockVectorType, Constness>::operator>=(
1198 const Iterator<BlockVectorType, OtherConstness> &i) const
1199 {
1200 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1201
1202 return (global_index >= i.global_index);
1203 }
1204
1205
1206
1207 template <class BlockVectorType, bool Constness>
1208 template <bool OtherConstness>
1209 inline typename Iterator<BlockVectorType, Constness>::difference_type
1210 Iterator<BlockVectorType, Constness>::operator-(
1211 const Iterator<BlockVectorType, OtherConstness> &i) const
1212 {
1213 Assert(parent == i.parent, ExcPointerToDifferentVectors());
1214
1215 return (static_cast<signed int>(global_index) -
1216 static_cast<signed int>(i.global_index));
1217 }
1218
1219
1220
1221 template <class BlockVectorType, bool Constness>
1222 inline Iterator<BlockVectorType, Constness>
1223 Iterator<BlockVectorType, Constness>::operator+(
1224 const difference_type &d) const
1225 {
1226 // if the index pointed to is
1227 // still within the block we
1228 // currently point into, then we
1229 // can save the computation of
1230 // the block
1231 if ((global_index + d >= next_break_backward) &&
1232 (global_index + d <= next_break_forward))
1233 return Iterator(*parent,
1234 global_index + d,
1235 current_block,
1236 index_within_block + d,
1237 next_break_forward,
1238 next_break_backward);
1239 else
1240 // outside present block, so
1241 // have to seek new block
1242 // anyway
1243 return Iterator(*parent, global_index + d);
1244 }
1245
1246
1247
1248 template <class BlockVectorType, bool Constness>
1249 inline Iterator<BlockVectorType, Constness>
1250 Iterator<BlockVectorType, Constness>::operator-(
1251 const difference_type &d) const
1252 {
1253 // if the index pointed to is
1254 // still within the block we
1255 // currently point into, then we
1256 // can save the computation of
1257 // the block
1258 if ((global_index - d >= next_break_backward) &&
1259 (global_index - d <= next_break_forward))
1260 return Iterator(*parent,
1261 global_index - d,
1262 current_block,
1263 index_within_block - d,
1264 next_break_forward,
1265 next_break_backward);
1266 else
1267 // outside present block, so
1268 // have to seek new block
1269 // anyway
1270 return Iterator(*parent, global_index - d);
1271 }
1272
1273
1274
1275 template <class BlockVectorType, bool Constness>
1276 inline Iterator<BlockVectorType, Constness> &
1277 Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1278 {
1279 // if the index pointed to is
1280 // still within the block we
1281 // currently point into, then we
1282 // can save the computation of
1283 // the block
1284 if ((global_index + d >= next_break_backward) &&
1285 (global_index + d <= next_break_forward))
1286 {
1287 global_index += d;
1288 index_within_block += d;
1289 }
1290 else
1291 // outside present block, so
1292 // have to seek new block
1293 // anyway
1294 *this = Iterator(*parent, global_index + d);
1295
1296 return *this;
1297 }
1298
1299
1300
1301 template <class BlockVectorType, bool Constness>
1302 inline Iterator<BlockVectorType, Constness> &
1303 Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1304 {
1305 // if the index pointed to is
1306 // still within the block we
1307 // currently point into, then we
1308 // can save the computation of
1309 // the block
1310 if ((global_index - d >= next_break_backward) &&
1311 (global_index - d <= next_break_forward))
1312 {
1313 global_index -= d;
1314 index_within_block -= d;
1315 }
1316 else
1317 // outside present block, so
1318 // have to seek new block
1319 // anyway
1320 *this = Iterator(*parent, global_index - d);
1321
1322 return *this;
1323 }
1324
1325
1326 template <class BlockVectorType, bool Constness>
1327 Iterator<BlockVectorType, Constness>::Iterator(BlockVector & parent,
1328 const size_type global_index)
1329 : parent(&parent)
1330 , global_index(global_index)
1331 {
1332 // find which block we are
1333 // in. for this, take into
1334 // account that it happens at
1335 // times that people want to
1336 // initialize iterators
1337 // past-the-end
1338 if (global_index < parent.size())
1339 {
1340 const std::pair<size_type, size_type> indices =
1341 parent.block_indices.global_to_local(global_index);
1342 current_block = indices.first;
1343 index_within_block = indices.second;
1344
1345 next_break_backward =
1346 parent.block_indices.local_to_global(current_block, 0);
1347 next_break_forward =
1348 (parent.block_indices.local_to_global(current_block, 0) +
1349 parent.block_indices.block_size(current_block) - 1);
1350 }
1351 else
1352 // past the end. only have one
1353 // value for this
1354 {
1355 this->global_index = parent.size();
1356 current_block = parent.n_blocks();
1357 index_within_block = 0;
1358 next_break_backward = global_index;
1359 next_break_forward = numbers::invalid_size_type;
1360 };
1361 }
1362
1363
1364
1365 template <class BlockVectorType, bool Constness>
1366 void
1367 Iterator<BlockVectorType, Constness>::move_forward()
1368 {
1369 if (global_index != next_break_forward)
1370 ++index_within_block;
1371 else
1372 {
1373 // ok, we traverse a boundary
1374 // between blocks:
1375 index_within_block = 0;
1376 ++current_block;
1377
1378 // break backwards is now old
1379 // break forward
1380 next_break_backward = next_break_forward + 1;
1381
1382 // compute new break forward
1383 if (current_block < parent->block_indices.size())
1384 next_break_forward +=
1385 parent->block_indices.block_size(current_block);
1386 else
1387 // if we are beyond the end,
1388 // then move the next
1389 // boundary arbitrarily far
1390 // away
1391 next_break_forward = numbers::invalid_size_type;
1392 };
1393
1394 ++global_index;
1395 }
1396
1397
1398
1399 template <class BlockVectorType, bool Constness>
1400 void
1401 Iterator<BlockVectorType, Constness>::move_backward()
1402 {
1403 if (global_index != next_break_backward)
1404 --index_within_block;
1405 else if (current_block != 0)
1406 {
1407 // ok, we traverse a boundary
1408 // between blocks:
1409 --current_block;
1410 index_within_block =
1411 parent->block_indices.block_size(current_block) - 1;
1412
1413 // break forwards is now old
1414 // break backward
1415 next_break_forward = next_break_backward - 1;
1416
1417 // compute new break forward
1418 next_break_backward -=
1419 parent->block_indices.block_size(current_block);
1420 }
1421 else
1422 // current block was 0, we now
1423 // get into unspecified terrain
1424 {
1425 --current_block;
1426 index_within_block = numbers::invalid_size_type;
1427 next_break_forward = 0;
1428 next_break_backward = 0;
1429 };
1430
1431 --global_index;
1432 }
1433
1434
1435 } // namespace BlockVectorIterators
1436
1437} // namespace internal
1438
1439
1440
1441template <class VectorType>
1442inline std::size_t
1444{
1445 return block_indices.total_size();
1446}
1447
1448
1449
1450template <class VectorType>
1451inline std::size_t
1453{
1454 std::size_t local_size = 0;
1455 for (unsigned int b = 0; b < n_blocks(); ++b)
1456 local_size += block(b).locally_owned_size();
1457 return local_size;
1458}
1459
1460
1461
1462template <class VectorType>
1463inline IndexSet
1465{
1466 IndexSet is(size());
1467
1468 // copy index sets from blocks into the global one, shifted
1469 // by the appropriate amount for each block
1470 for (unsigned int b = 0; b < n_blocks(); ++b)
1471 {
1472 IndexSet x = block(b).locally_owned_elements();
1473 is.add_indices(x, block_indices.block_start(b));
1474 }
1475
1476 is.compress();
1477
1478 return is;
1479}
1480
1481
1482
1483template <class VectorType>
1484inline unsigned int
1486{
1487 return block_indices.size();
1488}
1489
1490
1491template <class VectorType>
1493BlockVectorBase<VectorType>::block(const unsigned int i)
1494{
1496
1497 return components[i];
1498}
1499
1500
1501
1502template <class VectorType>
1503inline const typename BlockVectorBase<VectorType>::BlockType &
1504BlockVectorBase<VectorType>::block(const unsigned int i) const
1505{
1507
1508 return components[i];
1509}
1510
1511
1512
1513template <class VectorType>
1514inline const BlockIndices &
1516{
1517 return block_indices;
1518}
1519
1520
1521template <class VectorType>
1522inline void
1524{
1525 std::vector<size_type> sizes(n_blocks());
1526
1527 for (size_type i = 0; i < n_blocks(); ++i)
1528 sizes[i] = block(i).size();
1529
1530 block_indices.reinit(sizes);
1531}
1532
1533
1534
1535template <class VectorType>
1536inline void
1538{
1539 for (unsigned int i = 0; i < n_blocks(); ++i)
1540 block(i).compress(operation);
1541}
1542
1543
1544
1545template <class VectorType>
1548{
1549 return iterator(*this, 0U);
1550}
1551
1552
1553
1554template <class VectorType>
1557{
1558 return const_iterator(*this, 0U);
1559}
1560
1561
1562template <class VectorType>
1565{
1566 return iterator(*this, size());
1567}
1568
1569
1570
1571template <class VectorType>
1574{
1575 return const_iterator(*this, size());
1576}
1577
1578
1579template <class VectorType>
1580inline bool
1582{
1583 const std::pair<size_type, size_type> local_index =
1584 block_indices.global_to_local(global_index);
1585
1586 return components[local_index.first].in_local_range(global_index);
1587}
1588
1589
1590template <class VectorType>
1591bool
1593{
1594 for (size_type i = 0; i < n_blocks(); ++i)
1595 if (components[i].all_zero() == false)
1596 return false;
1597
1598 return true;
1599}
1600
1601
1602
1603template <class VectorType>
1604bool
1606{
1607 for (size_type i = 0; i < n_blocks(); ++i)
1608 if (components[i].is_non_negative() == false)
1609 return false;
1610
1611 return true;
1612}
1613
1614
1615
1616template <class VectorType>
1619 const BlockVectorBase<VectorType> &v) const
1620{
1621 Assert(n_blocks() == v.n_blocks(),
1623
1624 value_type sum = 0.;
1625 for (size_type i = 0; i < n_blocks(); ++i)
1626 sum += components[i] * v.components[i];
1627
1628 return sum;
1629}
1630
1631
1632template <class VectorType>
1635{
1636 real_type sum = 0.;
1637 for (size_type i = 0; i < n_blocks(); ++i)
1638 sum += components[i].norm_sqr();
1639
1640 return sum;
1641}
1642
1643
1644
1645template <class VectorType>
1648{
1649 value_type sum = 0.;
1650 // need to do static_cast as otherwise it won't work with
1651 // value_type=complex<T>
1652 for (size_type i = 0; i < n_blocks(); ++i)
1653 sum += components[i].mean_value() *
1655 components[i].size()));
1656
1658}
1659
1660
1661
1662template <class VectorType>
1665{
1666 real_type sum = 0.;
1667 for (size_type i = 0; i < n_blocks(); ++i)
1668 sum += components[i].l1_norm();
1669
1670 return sum;
1671}
1672
1673
1674
1675template <class VectorType>
1678{
1679 return std::sqrt(norm_sqr());
1680}
1681
1682
1683
1684template <class VectorType>
1687{
1688 real_type sum = 0.;
1689 for (size_type i = 0; i < n_blocks(); ++i)
1690 {
1691 value_type newval = components[i].linfty_norm();
1692 if (sum < newval)
1693 sum = newval;
1694 }
1695 return sum;
1696}
1697
1698
1699
1700template <class VectorType>
1706{
1709
1710 value_type sum = 0.;
1711 for (size_type i = 0; i < n_blocks(); ++i)
1712 sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1713
1714 return sum;
1715}
1716
1717
1718
1719template <class VectorType>
1722{
1723 Assert(n_blocks() == v.n_blocks(),
1725
1726 for (size_type i = 0; i < n_blocks(); ++i)
1727 {
1728 components[i] += v.components[i];
1729 }
1730
1731 return *this;
1732}
1733
1734
1735
1736template <class VectorType>
1739{
1740 Assert(n_blocks() == v.n_blocks(),
1742
1743 for (size_type i = 0; i < n_blocks(); ++i)
1744 {
1745 components[i] -= v.components[i];
1746 }
1747 return *this;
1748}
1749
1750
1751
1752template <class VectorType>
1753template <typename Number>
1754inline void
1755BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1756 const std::vector<Number> & values)
1757{
1758 Assert(indices.size() == values.size(),
1759 ExcDimensionMismatch(indices.size(), values.size()));
1760 add(indices.size(), indices.data(), values.data());
1761}
1762
1763
1764
1765template <class VectorType>
1766template <typename Number>
1767inline void
1768BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1769 const Vector<Number> & values)
1770{
1771 Assert(indices.size() == values.size(),
1772 ExcDimensionMismatch(indices.size(), values.size()));
1773 const size_type n_indices = indices.size();
1774 for (size_type i = 0; i < n_indices; ++i)
1775 (*this)(indices[i]) += values(i);
1776}
1777
1778
1779
1780template <class VectorType>
1781template <typename Number>
1782inline void
1784 const size_type *indices,
1785 const Number * values)
1786{
1787 for (size_type i = 0; i < n_indices; ++i)
1788 (*this)(indices[i]) += values[i];
1789}
1790
1791
1792
1793template <class VectorType>
1794void
1796{
1797 AssertIsFinite(a);
1798
1799 for (size_type i = 0; i < n_blocks(); ++i)
1800 {
1801 components[i].add(a);
1802 }
1803}
1804
1805
1806
1807template <class VectorType>
1808void
1811{
1812 AssertIsFinite(a);
1813
1814 Assert(n_blocks() == v.n_blocks(),
1816
1817 for (size_type i = 0; i < n_blocks(); ++i)
1818 {
1819 components[i].add(a, v.components[i]);
1820 }
1821}
1822
1823
1824
1825template <class VectorType>
1826void
1829 const value_type b,
1831{
1832 AssertIsFinite(a);
1833 AssertIsFinite(b);
1834
1835 Assert(n_blocks() == v.n_blocks(),
1837 Assert(n_blocks() == w.n_blocks(),
1838 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1839
1840
1841 for (size_type i = 0; i < n_blocks(); ++i)
1842 {
1843 components[i].add(a, v.components[i], b, w.components[i]);
1844 }
1845}
1846
1847
1848
1849template <class VectorType>
1850void
1853{
1854 AssertIsFinite(x);
1855
1856 Assert(n_blocks() == v.n_blocks(),
1858
1859 for (size_type i = 0; i < n_blocks(); ++i)
1860 {
1861 components[i].sadd(x, v.components[i]);
1862 }
1863}
1864
1865
1866
1867template <class VectorType>
1868void
1870 const value_type a,
1872{
1873 AssertIsFinite(x);
1874 AssertIsFinite(a);
1875
1876 Assert(n_blocks() == v.n_blocks(),
1878
1879 for (size_type i = 0; i < n_blocks(); ++i)
1880 {
1881 components[i].sadd(x, a, v.components[i]);
1882 }
1883}
1884
1885
1886
1887template <class VectorType>
1888void
1890 const value_type a,
1892 const value_type b,
1894{
1895 AssertIsFinite(x);
1896 AssertIsFinite(a);
1897 AssertIsFinite(b);
1898
1899 Assert(n_blocks() == v.n_blocks(),
1901 Assert(n_blocks() == w.n_blocks(),
1902 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1903
1904 for (size_type i = 0; i < n_blocks(); ++i)
1905 {
1906 components[i].sadd(x, a, v.components[i], b, w.components[i]);
1907 }
1908}
1909
1910
1911
1912template <class VectorType>
1913void
1915 const value_type a,
1917 const value_type b,
1919 const value_type c,
1921{
1922 AssertIsFinite(x);
1923 AssertIsFinite(a);
1924 AssertIsFinite(b);
1925 AssertIsFinite(c);
1926
1927 Assert(n_blocks() == v.n_blocks(),
1929 Assert(n_blocks() == w.n_blocks(),
1930 ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1931 Assert(n_blocks() == y.n_blocks(),
1933
1934 for (size_type i = 0; i < n_blocks(); ++i)
1935 {
1936 components[i].sadd(
1937 x, a, v.components[i], b, w.components[i], c, y.components[i]);
1938 }
1939}
1940
1941
1942
1943template <class VectorType>
1944template <class BlockVector2>
1945void
1946BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
1947{
1948 Assert(n_blocks() == v.n_blocks(),
1949 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1950 for (size_type i = 0; i < n_blocks(); ++i)
1951 components[i].scale(v.block(i));
1952}
1953
1954
1955
1956template <class VectorType>
1957std::size_t
1959{
1960 return (MemoryConsumption::memory_consumption(this->block_indices) +
1961 MemoryConsumption::memory_consumption(this->components));
1962}
1963
1964
1965
1966template <class VectorType>
1967template <class BlockVector2>
1968void
1969BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
1970{
1971 AssertIsFinite(a);
1972
1973 Assert(n_blocks() == v.n_blocks(),
1974 ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1975
1976 for (size_type i = 0; i < n_blocks(); ++i)
1977 components[i].equ(a, v.components[i]);
1978}
1979
1980
1981
1982template <class VectorType>
1983void
1985{
1986 for (size_type i = 0; i < n_blocks(); ++i)
1987 block(i).update_ghost_values();
1988}
1989
1990
1991
1992template <class VectorType>
1995{
1996 AssertIsFinite(s);
1997
1998 for (size_type i = 0; i < n_blocks(); ++i)
1999 components[i] = s;
2000
2001 return *this;
2002}
2003
2004
2005template <class VectorType>
2008{
2010
2011 for (size_type i = 0; i < n_blocks(); ++i)
2012 components[i] = v.components[i];
2013
2014 return *this;
2015}
2016
2017
2018template <class VectorType>
2019template <class VectorType2>
2022{
2024
2025 for (size_type i = 0; i < n_blocks(); ++i)
2026 components[i] = v.components[i];
2027
2028 return *this;
2029}
2030
2031
2032
2033template <class VectorType>
2035BlockVectorBase<VectorType>::operator=(const VectorType &v)
2036{
2037 Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2038
2039 size_type index_v = 0;
2040 for (size_type b = 0; b < n_blocks(); ++b)
2041 for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2042 block(b)(i) = v(index_v);
2043
2044 return *this;
2045}
2046
2047
2048
2049template <class VectorType>
2050template <class VectorType2>
2051inline bool
2053 const BlockVectorBase<VectorType2> &v) const
2054{
2056
2057 for (size_type i = 0; i < n_blocks(); ++i)
2058 if (!(components[i] == v.components[i]))
2059 return false;
2060
2061 return true;
2062}
2063
2064
2065
2066template <class VectorType>
2069{
2070 AssertIsFinite(factor);
2071
2072 for (size_type i = 0; i < n_blocks(); ++i)
2073 components[i] *= factor;
2074
2075 return *this;
2076}
2077
2078
2079
2080template <class VectorType>
2083{
2084 AssertIsFinite(factor);
2085 Assert(factor != 0., ExcDivideByZero());
2086
2087 const value_type inverse_factor = value_type(1.) / factor;
2088
2089 for (size_type i = 0; i < n_blocks(); ++i)
2090 components[i] *= inverse_factor;
2091
2092 return *this;
2093}
2094
2095
2096template <class VectorType>
2099{
2100 const std::pair<unsigned int, size_type> local_index =
2102 return components[local_index.first](local_index.second);
2103}
2104
2105
2106
2107template <class VectorType>
2110{
2111 const std::pair<unsigned int, size_type> local_index =
2113 return components[local_index.first](local_index.second);
2114}
2115
2116
2117
2118template <class VectorType>
2121{
2122 return operator()(i);
2123}
2124
2125
2126
2127template <class VectorType>
2130{
2131 return operator()(i);
2132}
2133
2134
2135
2136template <typename VectorType>
2137template <typename OtherNumber>
2138inline void
2140 const std::vector<size_type> &indices,
2141 std::vector<OtherNumber> & values) const
2142{
2143 for (size_type i = 0; i < indices.size(); ++i)
2144 values[i] = operator()(indices[i]);
2145}
2146
2147
2148
2149template <typename VectorType>
2150template <typename ForwardIterator, typename OutputIterator>
2151inline void
2153 ForwardIterator indices_begin,
2154 const ForwardIterator indices_end,
2155 OutputIterator values_begin) const
2156{
2157 while (indices_begin != indices_end)
2158 {
2159 *values_begin = operator()(*indices_begin);
2160 indices_begin++;
2161 values_begin++;
2162 }
2163}
2164
2165#endif // DOXYGEN
2166
2168
2169#endif
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
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)
std::size_t size() const
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)
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
typename std::conditional< Constness, const typename BlockVectorType::value_type, typename BlockVectorType::value_type >::type value_type
Iterator & operator-=(const difference_type &d)
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#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:488
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< T >::value, 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)
decltype(std::declval< T const >().block(0)) has_block_t
constexpr bool has_n_blocks
decltype(std::declval< T const >().n_blocks()) has_n_blocks_t
constexpr bool is_block_vector
constexpr bool has_block
const types::global_dof_index invalid_size_type
Definition types.h:222
STL namespace.
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition types.h:82