Reference documentation for deal.II version GIT c00406db16 2023-09-28 18:00:02+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\}}\)
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 
24 #include <deal.II/base/numbers.h>
26 
28 #include <deal.II/lac/exceptions.h>
30 #include <deal.II/lac/vector.h>
32 
33 #include <cmath>
34 #include <cstddef>
35 #include <iterator>
36 #include <type_traits>
37 #include <vector>
38 
40 
41 
47 namespace internal
48 {
49  template <typename T>
50  using has_block_t = decltype(std::declval<const T>().block(0));
51 
52  template <typename T>
53  constexpr bool has_block = internal::is_supported_operation<has_block_t, T>;
54 
55  template <typename T>
56  using has_n_blocks_t = decltype(std::declval<const T>().n_blocks());
57 
58  template <typename T>
59  constexpr bool has_n_blocks =
60  internal::is_supported_operation<has_n_blocks_t, T>;
61 
62  template <typename T>
63  constexpr bool is_block_vector = has_block<T> && has_n_blocks<T>;
64 } // namespace internal
65 
80 template <typename VectorType>
82 {
83 public:
89  static const bool value = internal::is_block_vector<VectorType>;
90 };
91 
92 
93 // instantiation of the static member
94 template <typename VectorType>
96 
97 
98 
99 namespace internal
100 {
104  namespace BlockVectorIterators
105  {
121  template <typename BlockVectorType, bool Constness>
122  class Iterator
123  {
124  public:
129 
135  using value_type =
136  typename std::conditional<Constness,
137  const typename BlockVectorType::value_type,
138  typename BlockVectorType::value_type>::type;
139 
153  using iterator_category = std::random_access_iterator_tag;
154  using difference_type = std::ptrdiff_t;
155  using reference = typename BlockVectorType::reference;
156  using pointer = value_type *;
157 
158  using dereference_type = typename std::conditional<
159  Constness,
160  value_type,
161  typename BlockVectorType::BlockType::reference>::type;
162 
167  using BlockVector = typename std::
168  conditional<Constness, const BlockVectorType, BlockVectorType>::type;
169 
179 
188 
189 
193  Iterator(const Iterator &c);
194 
195  private:
201  const size_type global_index,
202  const size_type current_block,
206 
207  public:
211  Iterator &
212  operator=(const Iterator &c);
213 
220  operator*() const;
221 
228 
233  Iterator &
235 
241  Iterator
243 
248  Iterator &
250 
256  Iterator
258 
263  template <bool OtherConstness>
264  bool
266 
271  template <bool OtherConstness>
272  bool
274 
280  template <bool OtherConstness>
281  bool
283 
287  template <bool OtherConstness>
288  bool
290 
294  template <bool OtherConstness>
295  bool
297 
301  template <bool OtherConstness>
302  bool
304 
308  template <bool OtherConstness>
311 
316  Iterator
317  operator+(const difference_type &d) const;
318 
323  Iterator
324  operator-(const difference_type &d) const;
325 
330  Iterator &
332 
337  Iterator &
339 
350  "Your program tried to compare iterators pointing to "
351  "different block vectors. There is no reasonable way "
352  "to do this.");
353 
355  private:
362 
367 
372  unsigned int current_block;
374 
383 
387  void
389 
393  void
395 
396  // Mark all other instances of this template as friends.
397  template <typename, bool>
398  friend class Iterator;
399  };
400  } // namespace BlockVectorIterators
401 } // namespace internal
402 
403 
441 template <typename VectorType>
443  public ReadVector<typename VectorType::value_type>
444 {
445 public:
449  using BlockType = VectorType;
450 
451  /*
452  * Declare standard types used in
453  * all containers. These types
454  * parallel those in the
455  * <tt>C++</tt> standard
456  * libraries
457  * <tt>std::vector<...></tt>
458  * class. This includes iterator
459  * types.
460  */
461  using value_type = typename BlockType::value_type;
462  using pointer = value_type *;
463  using const_pointer = const value_type *;
464  using iterator =
468  using reference = typename BlockType::reference;
469  using const_reference = typename BlockType::const_reference;
471 
481  using real_type = typename BlockType::real_type;
482 
486  BlockVectorBase() = default;
487 
491  BlockVectorBase(const BlockVectorBase & /*V*/) = default;
492 
498  BlockVectorBase(BlockVectorBase && /*V*/) noexcept = default;
499 
506  void
508 
519  void
521 
525  BlockType &
526  block(const unsigned int i);
527 
531  const BlockType &
532  block(const unsigned int i) const;
533 
539  const BlockIndices &
541 
545  unsigned int
546  n_blocks() const;
547 
552  virtual size_type
553  size() const override;
554 
560  std::size_t
562 
579  IndexSet
581 
585  iterator
586  begin();
587 
593  begin() const;
594 
598  iterator
599  end();
600 
606  end() const;
607 
611  value_type
612  operator()(const size_type i) const;
613 
617  reference
618  operator()(const size_type i);
619 
625  value_type
626  operator[](const size_type i) const;
627 
633  reference
634  operator[](const size_type i);
635 
651  template <typename OtherNumber>
652  void
653  extract_subvector_to(const std::vector<size_type> &indices,
654  std::vector<OtherNumber> &values) const;
655 
656  virtual void
658  ArrayView<value_type> &entries) const override;
659 
687  template <typename ForwardIterator, typename OutputIterator>
688  void
689  extract_subvector_to(ForwardIterator indices_begin,
690  const ForwardIterator indices_end,
691  OutputIterator values_begin) const;
692 
698  operator=(const value_type s);
699 
704  operator=(const BlockVectorBase &V);
705 
712  operator=(BlockVectorBase && /*V*/) = default; // NOLINT
713 
717  template <typename VectorType2>
719  operator=(const BlockVectorBase<VectorType2> &V);
720 
725  operator=(const VectorType &v);
726 
731  template <typename VectorType2>
732  bool
733  operator==(const BlockVectorBase<VectorType2> &v) const;
734 
738  value_type
739  operator*(const BlockVectorBase &V) const;
740 
744  real_type
745  norm_sqr() const;
746 
750  value_type
751  mean_value() const;
752 
756  real_type
757  l1_norm() const;
758 
763  real_type
764  l2_norm() const;
765 
770  real_type
771  linfty_norm() const;
772 
795  value_type
797  const BlockVectorBase &V,
798  const BlockVectorBase &W);
799 
804  bool
806 
812  bool
813  all_zero() const;
814 
820  bool
822 
827  operator+=(const BlockVectorBase &V);
828 
833  operator-=(const BlockVectorBase &V);
834 
835 
840  template <typename Number>
841  void
842  add(const std::vector<size_type> &indices, const std::vector<Number> &values);
843 
848  template <typename Number>
849  void
850  add(const std::vector<size_type> &indices, const Vector<Number> &values);
851 
857  template <typename Number>
858  void
859  add(const size_type n_elements,
860  const size_type *indices,
861  const Number *values);
862 
867  void
868  add(const value_type s);
869 
873  void
874  add(const value_type a, const BlockVectorBase &V);
875 
879  void
880  add(const value_type a,
881  const BlockVectorBase &V,
882  const value_type b,
883  const BlockVectorBase &W);
884 
888  void
889  sadd(const value_type s, const BlockVectorBase &V);
890 
894  void
895  sadd(const value_type s, const value_type a, const BlockVectorBase &V);
896 
900  void
901  sadd(const value_type s,
902  const value_type a,
903  const BlockVectorBase &V,
904  const value_type b,
905  const BlockVectorBase &W);
906 
910  void
911  sadd(const value_type s,
912  const value_type a,
913  const BlockVectorBase &V,
914  const value_type b,
915  const BlockVectorBase &W,
916  const value_type c,
917  const BlockVectorBase &X);
918 
923  operator*=(const value_type factor);
924 
929  operator/=(const value_type factor);
930 
935  template <class BlockVector2>
936  void
937  scale(const BlockVector2 &v);
938 
942  template <class BlockVector2>
943  void
944  equ(const value_type a, const BlockVector2 &V);
945 
950  void
952 
957  std::size_t
959 
960 protected:
964  std::vector<VectorType> components;
965 
971 
972  // Make the iterator class a friend.
973  template <typename N, bool C>
974  friend class ::internal::BlockVectorIterators::Iterator;
975 
976  template <typename>
977  friend class BlockVectorBase;
978 };
979 
980 
983 /*----------------------- Inline functions ----------------------------------*/
984 
985 
986 #ifndef DOXYGEN
987 namespace internal
988 {
989  namespace BlockVectorIterators
990  {
991  template <typename BlockVectorType, bool Constness>
992  inline Iterator<BlockVectorType, Constness>::Iterator(
993  const Iterator<BlockVectorType, Constness> &c)
994  : parent(c.parent)
995  , global_index(c.global_index)
996  , current_block(c.current_block)
997  , index_within_block(c.index_within_block)
998  , next_break_forward(c.next_break_forward)
999  , next_break_backward(c.next_break_backward)
1000  {}
1001 
1002 
1003 
1004  template <typename BlockVectorType, bool Constness>
1005  inline Iterator<BlockVectorType, Constness>::Iterator(
1006  const Iterator<BlockVectorType, !Constness> &c)
1007  : parent(c.parent)
1008  , global_index(c.global_index)
1009  , current_block(c.current_block)
1010  , index_within_block(c.index_within_block)
1011  , next_break_forward(c.next_break_forward)
1012  , next_break_backward(c.next_break_backward)
1013  {
1014  // Only permit copy-constructing const iterators from non-const
1015  // iterators, and not vice versa (i.e., Constness must always be
1016  // true).
1017  static_assert(Constness == true,
1018  "Constructing a non-const iterator from a const iterator "
1019  "does not make sense.");
1020  }
1021 
1022 
1023 
1024  template <typename BlockVectorType, bool Constness>
1025  inline Iterator<BlockVectorType, Constness>::Iterator(
1026  BlockVector &parent,
1027  const size_type global_index,
1028  const size_type current_block,
1029  const size_type index_within_block,
1030  const size_type next_break_forward,
1031  const size_type next_break_backward)
1032  : parent(&parent)
1034  , current_block(current_block)
1035  , index_within_block(index_within_block)
1036  , next_break_forward(next_break_forward)
1037  , next_break_backward(next_break_backward)
1038  {}
1039 
1040 
1041 
1042  template <typename BlockVectorType, bool Constness>
1043  inline Iterator<BlockVectorType, Constness> &
1044  Iterator<BlockVectorType, Constness>::operator=(const Iterator &c)
1045  {
1046  parent = c.parent;
1047  global_index = c.global_index;
1048  index_within_block = c.index_within_block;
1049  current_block = c.current_block;
1050  next_break_forward = c.next_break_forward;
1051  next_break_backward = c.next_break_backward;
1052 
1053  return *this;
1054  }
1055 
1056 
1057 
1058  template <typename BlockVectorType, bool Constness>
1059  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1061  {
1062  return parent->block(current_block)(index_within_block);
1063  }
1064 
1065 
1066 
1067  template <typename BlockVectorType, bool Constness>
1068  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1069  Iterator<BlockVectorType, Constness>::operator[](
1070  const difference_type d) const
1071  {
1072  // if the index pointed to is
1073  // still within the block we
1074  // currently point into, then we
1075  // can save the computation of
1076  // the block
1077  if ((global_index + d >= next_break_backward) &&
1078  (global_index + d <= next_break_forward))
1079  return parent->block(current_block)(index_within_block + d);
1080 
1081  // if the index is not within the
1082  // block of the block vector into
1083  // which we presently point, then
1084  // there is no way: we have to
1085  // search for the block. this can
1086  // be done through the parent
1087  // class as well.
1088  return (*parent)(global_index + d);
1089  }
1090 
1091 
1092 
1093  template <typename BlockVectorType, bool Constness>
1094  inline Iterator<BlockVectorType, Constness> &
1096  {
1097  move_forward();
1098  return *this;
1099  }
1100 
1101 
1102 
1103  template <typename BlockVectorType, bool Constness>
1104  inline Iterator<BlockVectorType, Constness>
1106  {
1107  const Iterator old_value = *this;
1108  move_forward();
1109  return old_value;
1110  }
1111 
1112 
1113 
1114  template <typename BlockVectorType, bool Constness>
1115  inline Iterator<BlockVectorType, Constness> &
1116  Iterator<BlockVectorType, Constness>::operator--()
1117  {
1118  move_backward();
1119  return *this;
1120  }
1121 
1122 
1123 
1124  template <typename BlockVectorType, bool Constness>
1125  inline Iterator<BlockVectorType, Constness>
1126  Iterator<BlockVectorType, Constness>::operator--(int)
1127  {
1128  const Iterator old_value = *this;
1129  move_backward();
1130  return old_value;
1131  }
1132 
1133 
1134 
1135  template <typename BlockVectorType, bool Constness>
1136  template <bool OtherConstness>
1137  inline bool
1139  const Iterator<BlockVectorType, OtherConstness> &i) const
1140  {
1141  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1142 
1143  return (global_index == i.global_index);
1144  }
1145 
1146 
1147 
1148  template <typename BlockVectorType, bool Constness>
1149  template <bool OtherConstness>
1150  inline bool
1152  const Iterator<BlockVectorType, OtherConstness> &i) const
1153  {
1154  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1155 
1156  return (global_index != i.global_index);
1157  }
1158 
1159 
1160 
1161  template <typename BlockVectorType, bool Constness>
1162  template <bool OtherConstness>
1163  inline bool
1165  const Iterator<BlockVectorType, OtherConstness> &i) const
1166  {
1167  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1168 
1169  return (global_index < i.global_index);
1170  }
1171 
1172 
1173 
1174  template <typename BlockVectorType, bool Constness>
1175  template <bool OtherConstness>
1176  inline bool
1178  const Iterator<BlockVectorType, OtherConstness> &i) const
1179  {
1180  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1181 
1182  return (global_index <= i.global_index);
1183  }
1184 
1185 
1186 
1187  template <typename BlockVectorType, bool Constness>
1188  template <bool OtherConstness>
1189  inline bool
1191  const Iterator<BlockVectorType, OtherConstness> &i) const
1192  {
1193  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1194 
1195  return (global_index > i.global_index);
1196  }
1197 
1198 
1199 
1200  template <typename BlockVectorType, bool Constness>
1201  template <bool OtherConstness>
1202  inline bool
1204  const Iterator<BlockVectorType, OtherConstness> &i) const
1205  {
1206  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1207 
1208  return (global_index >= i.global_index);
1209  }
1210 
1211 
1212 
1213  template <typename BlockVectorType, bool Constness>
1214  template <bool OtherConstness>
1215  inline typename Iterator<BlockVectorType, Constness>::difference_type
1217  const Iterator<BlockVectorType, OtherConstness> &i) const
1218  {
1219  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1220 
1221  return (static_cast<signed int>(global_index) -
1222  static_cast<signed int>(i.global_index));
1223  }
1224 
1225 
1226 
1227  template <typename BlockVectorType, bool Constness>
1228  inline Iterator<BlockVectorType, Constness>
1230  const difference_type &d) const
1231  {
1232  // if the index pointed to is
1233  // still within the block we
1234  // currently point into, then we
1235  // can save the computation of
1236  // the block
1237  if ((global_index + d >= next_break_backward) &&
1238  (global_index + d <= next_break_forward))
1239  return Iterator(*parent,
1240  global_index + d,
1241  current_block,
1242  index_within_block + d,
1243  next_break_forward,
1244  next_break_backward);
1245  else
1246  // outside present block, so
1247  // have to seek new block
1248  // anyway
1249  return Iterator(*parent, global_index + d);
1250  }
1251 
1252 
1253 
1254  template <typename BlockVectorType, bool Constness>
1255  inline Iterator<BlockVectorType, Constness>
1257  const difference_type &d) const
1258  {
1259  // if the index pointed to is
1260  // still within the block we
1261  // currently point into, then we
1262  // can save the computation of
1263  // the block
1264  if ((global_index - d >= next_break_backward) &&
1265  (global_index - d <= next_break_forward))
1266  return Iterator(*parent,
1267  global_index - d,
1268  current_block,
1269  index_within_block - d,
1270  next_break_forward,
1271  next_break_backward);
1272  else
1273  // outside present block, so
1274  // have to seek new block
1275  // anyway
1276  return Iterator(*parent, global_index - d);
1277  }
1278 
1279 
1280 
1281  template <typename BlockVectorType, bool Constness>
1282  inline Iterator<BlockVectorType, Constness> &
1283  Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1284  {
1285  // if the index pointed to is
1286  // still within the block we
1287  // currently point into, then we
1288  // can save the computation of
1289  // the block
1290  if ((global_index + d >= next_break_backward) &&
1291  (global_index + d <= next_break_forward))
1292  {
1293  global_index += d;
1294  index_within_block += d;
1295  }
1296  else
1297  // outside present block, so
1298  // have to seek new block
1299  // anyway
1300  *this = Iterator(*parent, global_index + d);
1301 
1302  return *this;
1303  }
1304 
1305 
1306 
1307  template <typename BlockVectorType, bool Constness>
1308  inline Iterator<BlockVectorType, Constness> &
1309  Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1310  {
1311  // if the index pointed to is
1312  // still within the block we
1313  // currently point into, then we
1314  // can save the computation of
1315  // the block
1316  if ((global_index - d >= next_break_backward) &&
1317  (global_index - d <= next_break_forward))
1318  {
1319  global_index -= d;
1320  index_within_block -= d;
1321  }
1322  else
1323  // outside present block, so
1324  // have to seek new block
1325  // anyway
1326  *this = Iterator(*parent, global_index - d);
1327 
1328  return *this;
1329  }
1330 
1331 
1332  template <typename BlockVectorType, bool Constness>
1333  Iterator<BlockVectorType, Constness>::Iterator(BlockVector &parent,
1334  const size_type global_index)
1335  : parent(&parent)
1337  {
1338  // find which block we are
1339  // in. for this, take into
1340  // account that it happens at
1341  // times that people want to
1342  // initialize iterators
1343  // past-the-end
1344  if (global_index < parent.size())
1345  {
1346  const std::pair<size_type, size_type> indices =
1348  current_block = indices.first;
1349  index_within_block = indices.second;
1350 
1351  next_break_backward =
1352  parent.block_indices.local_to_global(current_block, 0);
1353  next_break_forward =
1354  (parent.block_indices.local_to_global(current_block, 0) +
1355  parent.block_indices.block_size(current_block) - 1);
1356  }
1357  else
1358  // past the end. only have one
1359  // value for this
1360  {
1361  this->global_index = parent.size();
1362  current_block = parent.n_blocks();
1363  index_within_block = 0;
1364  next_break_backward = global_index;
1365  next_break_forward = numbers::invalid_size_type;
1366  };
1367  }
1368 
1369 
1370 
1371  template <typename BlockVectorType, bool Constness>
1372  void
1373  Iterator<BlockVectorType, Constness>::move_forward()
1374  {
1375  if (global_index != next_break_forward)
1376  ++index_within_block;
1377  else
1378  {
1379  // ok, we traverse a boundary
1380  // between blocks:
1381  index_within_block = 0;
1382  ++current_block;
1383 
1384  // break backwards is now old
1385  // break forward
1386  next_break_backward = next_break_forward + 1;
1387 
1388  // compute new break forward
1389  if (current_block < parent->block_indices.size())
1390  next_break_forward +=
1391  parent->block_indices.block_size(current_block);
1392  else
1393  // if we are beyond the end,
1394  // then move the next
1395  // boundary arbitrarily far
1396  // away
1397  next_break_forward = numbers::invalid_size_type;
1398  };
1399 
1400  ++global_index;
1401  }
1402 
1403 
1404 
1405  template <typename BlockVectorType, bool Constness>
1406  void
1407  Iterator<BlockVectorType, Constness>::move_backward()
1408  {
1409  if (global_index != next_break_backward)
1410  --index_within_block;
1411  else if (current_block != 0)
1412  {
1413  // ok, we traverse a boundary
1414  // between blocks:
1415  --current_block;
1416  index_within_block =
1417  parent->block_indices.block_size(current_block) - 1;
1418 
1419  // break forwards is now old
1420  // break backward
1421  next_break_forward = next_break_backward - 1;
1422 
1423  // compute new break forward
1424  next_break_backward -=
1425  parent->block_indices.block_size(current_block);
1426  }
1427  else
1428  // current block was 0, we now
1429  // get into unspecified terrain
1430  {
1431  --current_block;
1432  index_within_block = numbers::invalid_size_type;
1433  next_break_forward = 0;
1434  next_break_backward = 0;
1435  };
1436 
1437  --global_index;
1438  }
1439 
1440 
1441  } // namespace BlockVectorIterators
1442 
1443 } // namespace internal
1444 
1445 
1446 
1447 template <typename VectorType>
1450 {
1451  return block_indices.total_size();
1452 }
1453 
1454 
1455 
1456 template <typename VectorType>
1457 inline std::size_t
1459 {
1460  std::size_t local_size = 0;
1461  for (unsigned int b = 0; b < n_blocks(); ++b)
1462  local_size += block(b).locally_owned_size();
1463  return local_size;
1464 }
1465 
1466 
1467 
1468 template <typename VectorType>
1469 inline IndexSet
1471 {
1472  IndexSet is(size());
1473 
1474  // copy index sets from blocks into the global one, shifted
1475  // by the appropriate amount for each block
1476  for (unsigned int b = 0; b < n_blocks(); ++b)
1477  {
1478  IndexSet x = block(b).locally_owned_elements();
1479  is.add_indices(x, block_indices.block_start(b));
1480  }
1481 
1482  is.compress();
1483 
1484  return is;
1485 }
1486 
1487 
1488 
1489 template <typename VectorType>
1490 inline unsigned int
1492 {
1493  return block_indices.size();
1494 }
1495 
1496 
1497 template <typename VectorType>
1499 BlockVectorBase<VectorType>::block(const unsigned int i)
1500 {
1501  AssertIndexRange(i, n_blocks());
1502 
1503  return components[i];
1504 }
1505 
1506 
1507 
1508 template <typename VectorType>
1509 inline const typename BlockVectorBase<VectorType>::BlockType &
1510 BlockVectorBase<VectorType>::block(const unsigned int i) const
1511 {
1512  AssertIndexRange(i, n_blocks());
1513 
1514  return components[i];
1515 }
1516 
1517 
1518 
1519 template <typename VectorType>
1520 inline const BlockIndices &
1522 {
1523  return block_indices;
1524 }
1525 
1526 
1527 template <typename VectorType>
1528 inline void
1530 {
1531  std::vector<size_type> sizes(n_blocks());
1532 
1533  for (size_type i = 0; i < n_blocks(); ++i)
1534  sizes[i] = block(i).size();
1535 
1536  block_indices.reinit(sizes);
1537 }
1538 
1539 
1540 
1541 template <typename VectorType>
1542 inline void
1544 {
1545  for (unsigned int i = 0; i < n_blocks(); ++i)
1546  block(i).compress(operation);
1547 }
1548 
1549 
1550 
1551 template <typename VectorType>
1554 {
1555  return iterator(*this, 0U);
1556 }
1557 
1558 
1559 
1560 template <typename VectorType>
1563 {
1564  return const_iterator(*this, 0U);
1565 }
1566 
1567 
1568 template <typename VectorType>
1571 {
1572  return iterator(*this, size());
1573 }
1574 
1575 
1576 
1577 template <typename VectorType>
1580 {
1581  return const_iterator(*this, size());
1582 }
1583 
1584 
1585 template <typename VectorType>
1586 inline bool
1588 {
1589  const std::pair<size_type, size_type> local_index =
1591 
1592  return components[local_index.first].in_local_range(global_index);
1593 }
1594 
1595 
1596 template <typename VectorType>
1597 bool
1599 {
1600  for (size_type i = 0; i < n_blocks(); ++i)
1601  if (components[i].all_zero() == false)
1602  return false;
1603 
1604  return true;
1605 }
1606 
1607 
1608 
1609 template <typename VectorType>
1610 bool
1612 {
1613  for (size_type i = 0; i < n_blocks(); ++i)
1614  if (components[i].is_non_negative() == false)
1615  return false;
1616 
1617  return true;
1618 }
1619 
1620 
1621 
1622 template <typename VectorType>
1625  const BlockVectorBase<VectorType> &v) const
1626 {
1627  Assert(n_blocks() == v.n_blocks(),
1629 
1630  value_type sum = 0.;
1631  for (size_type i = 0; i < n_blocks(); ++i)
1632  sum += components[i] * v.components[i];
1633 
1634  return sum;
1635 }
1636 
1637 
1638 template <typename VectorType>
1641 {
1642  real_type sum = 0.;
1643  for (size_type i = 0; i < n_blocks(); ++i)
1644  sum += components[i].norm_sqr();
1645 
1646  return sum;
1647 }
1648 
1649 
1650 
1651 template <typename VectorType>
1654 {
1655  value_type sum = 0.;
1656  // need to do static_cast as otherwise it won't work with
1657  // value_type=complex<T>
1658  for (size_type i = 0; i < n_blocks(); ++i)
1659  sum += components[i].mean_value() *
1661  components[i].size()));
1662 
1663  return sum / (typename numbers::NumberTraits<value_type>::real_type(size()));
1664 }
1665 
1666 
1667 
1668 template <typename VectorType>
1671 {
1672  real_type sum = 0.;
1673  for (size_type i = 0; i < n_blocks(); ++i)
1674  sum += components[i].l1_norm();
1675 
1676  return sum;
1677 }
1678 
1679 
1680 
1681 template <typename VectorType>
1684 {
1685  return std::sqrt(norm_sqr());
1686 }
1687 
1688 
1689 
1690 template <typename VectorType>
1693 {
1694  real_type sum = 0.;
1695  for (size_type i = 0; i < n_blocks(); ++i)
1696  {
1697  value_type newval = components[i].linfty_norm();
1698  if (sum < newval)
1699  sum = newval;
1700  }
1701  return sum;
1702 }
1703 
1704 
1705 
1706 template <typename VectorType>
1709  const typename BlockVectorBase<VectorType>::value_type a,
1711  const BlockVectorBase<VectorType> &W)
1712 {
1713  AssertDimension(n_blocks(), V.n_blocks());
1715 
1716  value_type sum = 0.;
1717  for (size_type i = 0; i < n_blocks(); ++i)
1718  sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1719 
1720  return sum;
1721 }
1722 
1723 
1724 
1725 template <typename VectorType>
1728 {
1729  Assert(n_blocks() == v.n_blocks(),
1731 
1732  for (size_type i = 0; i < n_blocks(); ++i)
1733  {
1734  components[i] += v.components[i];
1735  }
1736 
1737  return *this;
1738 }
1739 
1740 
1741 
1742 template <typename VectorType>
1745 {
1746  Assert(n_blocks() == v.n_blocks(),
1748 
1749  for (size_type i = 0; i < n_blocks(); ++i)
1750  {
1751  components[i] -= v.components[i];
1752  }
1753  return *this;
1754 }
1755 
1756 
1757 
1758 template <typename VectorType>
1759 template <typename Number>
1760 inline void
1761 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1762  const std::vector<Number> &values)
1763 {
1764  Assert(indices.size() == values.size(),
1765  ExcDimensionMismatch(indices.size(), values.size()));
1766  add(indices.size(), indices.data(), values.data());
1767 }
1768 
1769 
1770 
1771 template <typename VectorType>
1772 template <typename Number>
1773 inline void
1774 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1775  const Vector<Number> &values)
1776 {
1777  Assert(indices.size() == values.size(),
1778  ExcDimensionMismatch(indices.size(), values.size()));
1779  const size_type n_indices = indices.size();
1780  for (size_type i = 0; i < n_indices; ++i)
1781  (*this)(indices[i]) += values(i);
1782 }
1783 
1784 
1785 
1786 template <typename VectorType>
1787 template <typename Number>
1788 inline void
1790  const size_type *indices,
1791  const Number *values)
1792 {
1793  for (size_type i = 0; i < n_indices; ++i)
1794  (*this)(indices[i]) += values[i];
1795 }
1796 
1797 
1798 
1799 template <typename VectorType>
1800 void
1802 {
1803  AssertIsFinite(a);
1804 
1805  for (size_type i = 0; i < n_blocks(); ++i)
1806  {
1807  components[i].add(a);
1808  }
1809 }
1810 
1811 
1812 
1813 template <typename VectorType>
1814 void
1816  const BlockVectorBase<VectorType> &v)
1817 {
1818  AssertIsFinite(a);
1819 
1820  Assert(n_blocks() == v.n_blocks(),
1822 
1823  for (size_type i = 0; i < n_blocks(); ++i)
1824  {
1825  components[i].add(a, v.components[i]);
1826  }
1827 }
1828 
1829 
1830 
1831 template <typename VectorType>
1832 void
1834  const BlockVectorBase<VectorType> &v,
1835  const value_type b,
1837 {
1838  AssertIsFinite(a);
1839  AssertIsFinite(b);
1840 
1841  Assert(n_blocks() == v.n_blocks(),
1843  Assert(n_blocks() == w.n_blocks(),
1844  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1845 
1846 
1847  for (size_type i = 0; i < n_blocks(); ++i)
1848  {
1849  components[i].add(a, v.components[i], b, w.components[i]);
1850  }
1851 }
1852 
1853 
1854 
1855 template <typename VectorType>
1856 void
1858  const BlockVectorBase<VectorType> &v)
1859 {
1860  AssertIsFinite(x);
1861 
1862  Assert(n_blocks() == v.n_blocks(),
1864 
1865  for (size_type i = 0; i < n_blocks(); ++i)
1866  {
1867  components[i].sadd(x, v.components[i]);
1868  }
1869 }
1870 
1871 
1872 
1873 template <typename VectorType>
1874 void
1876  const value_type a,
1877  const BlockVectorBase<VectorType> &v)
1878 {
1879  AssertIsFinite(x);
1880  AssertIsFinite(a);
1881 
1882  Assert(n_blocks() == v.n_blocks(),
1884 
1885  for (size_type i = 0; i < n_blocks(); ++i)
1886  {
1887  components[i].sadd(x, a, v.components[i]);
1888  }
1889 }
1890 
1891 
1892 
1893 template <typename VectorType>
1894 void
1896  const value_type a,
1897  const BlockVectorBase<VectorType> &v,
1898  const value_type b,
1900 {
1901  AssertIsFinite(x);
1902  AssertIsFinite(a);
1903  AssertIsFinite(b);
1904 
1905  Assert(n_blocks() == v.n_blocks(),
1907  Assert(n_blocks() == w.n_blocks(),
1908  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1909 
1910  for (size_type i = 0; i < n_blocks(); ++i)
1911  {
1912  components[i].sadd(x, a, v.components[i], b, w.components[i]);
1913  }
1914 }
1915 
1916 
1917 
1918 template <typename VectorType>
1919 void
1921  const value_type a,
1922  const BlockVectorBase<VectorType> &v,
1923  const value_type b,
1925  const value_type c,
1926  const BlockVectorBase<VectorType> &y)
1927 {
1928  AssertIsFinite(x);
1929  AssertIsFinite(a);
1930  AssertIsFinite(b);
1931  AssertIsFinite(c);
1932 
1933  Assert(n_blocks() == v.n_blocks(),
1935  Assert(n_blocks() == w.n_blocks(),
1936  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1937  Assert(n_blocks() == y.n_blocks(),
1939 
1940  for (size_type i = 0; i < n_blocks(); ++i)
1941  {
1942  components[i].sadd(
1943  x, a, v.components[i], b, w.components[i], c, y.components[i]);
1944  }
1945 }
1946 
1947 
1948 
1949 template <typename VectorType>
1950 template <class BlockVector2>
1951 void
1952 BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
1953 {
1954  Assert(n_blocks() == v.n_blocks(),
1955  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1956  for (size_type i = 0; i < n_blocks(); ++i)
1957  components[i].scale(v.block(i));
1958 }
1959 
1960 
1961 
1962 template <typename VectorType>
1963 std::size_t
1965 {
1966  return (MemoryConsumption::memory_consumption(this->block_indices) +
1967  MemoryConsumption::memory_consumption(this->components));
1968 }
1969 
1970 
1971 
1972 template <typename VectorType>
1973 template <class BlockVector2>
1974 void
1975 BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
1976 {
1977  AssertIsFinite(a);
1978 
1979  Assert(n_blocks() == v.n_blocks(),
1980  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1981 
1982  for (size_type i = 0; i < n_blocks(); ++i)
1983  components[i].equ(a, v.components[i]);
1984 }
1985 
1986 
1987 
1988 template <typename VectorType>
1989 void
1991 {
1992  for (size_type i = 0; i < n_blocks(); ++i)
1993  block(i).update_ghost_values();
1994 }
1995 
1996 
1997 
1998 template <typename VectorType>
2001 {
2002  AssertIsFinite(s);
2003 
2004  for (size_type i = 0; i < n_blocks(); ++i)
2005  components[i] = s;
2006 
2007  return *this;
2008 }
2009 
2010 
2011 template <typename VectorType>
2014 {
2016 
2017  for (size_type i = 0; i < n_blocks(); ++i)
2018  components[i] = v.components[i];
2019 
2020  return *this;
2021 }
2022 
2023 
2024 template <typename VectorType>
2025 template <typename VectorType2>
2028 {
2030 
2031  for (size_type i = 0; i < n_blocks(); ++i)
2032  components[i] = v.components[i];
2033 
2034  return *this;
2035 }
2036 
2037 
2038 
2039 template <typename VectorType>
2041 BlockVectorBase<VectorType>::operator=(const VectorType &v)
2042 {
2043  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2044 
2045  size_type index_v = 0;
2046  for (size_type b = 0; b < n_blocks(); ++b)
2047  for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2048  block(b)(i) = v(index_v);
2049 
2050  return *this;
2051 }
2052 
2053 
2054 
2055 template <typename VectorType>
2056 template <typename VectorType2>
2057 inline bool
2059  const BlockVectorBase<VectorType2> &v) const
2060 {
2062 
2063  for (size_type i = 0; i < n_blocks(); ++i)
2064  if (!(components[i] == v.components[i]))
2065  return false;
2066 
2067  return true;
2068 }
2069 
2070 
2071 
2072 template <typename VectorType>
2075 {
2076  AssertIsFinite(factor);
2077 
2078  for (size_type i = 0; i < n_blocks(); ++i)
2079  components[i] *= factor;
2080 
2081  return *this;
2082 }
2083 
2084 
2085 
2086 template <typename VectorType>
2089 {
2090  AssertIsFinite(factor);
2091  Assert(factor != 0., ExcDivideByZero());
2092 
2093  const value_type inverse_factor = value_type(1.) / factor;
2094 
2095  for (size_type i = 0; i < n_blocks(); ++i)
2096  components[i] *= inverse_factor;
2097 
2098  return *this;
2099 }
2100 
2101 
2102 template <typename VectorType>
2105 {
2106  const std::pair<unsigned int, size_type> local_index =
2108  return components[local_index.first](local_index.second);
2109 }
2110 
2111 
2112 
2113 template <typename VectorType>
2116 {
2117  const std::pair<unsigned int, size_type> local_index =
2119  return components[local_index.first](local_index.second);
2120 }
2121 
2122 
2123 
2124 template <typename VectorType>
2127 {
2128  return operator()(i);
2129 }
2130 
2131 
2132 
2133 template <typename VectorType>
2136 {
2137  return operator()(i);
2138 }
2139 
2140 
2141 
2142 template <typename VectorType>
2143 template <typename OtherNumber>
2144 inline void
2146  const std::vector<size_type> &indices,
2147  std::vector<OtherNumber> &values) const
2148 {
2149  for (size_type i = 0; i < indices.size(); ++i)
2150  values[i] = operator()(indices[i]);
2151 }
2152 
2153 
2154 
2155 template <typename VectorType>
2156 inline void
2159  ArrayView<value_type> &entries) const
2160 {
2161  AssertDimension(indices.size(), entries.size());
2162  for (unsigned int i = 0; i < indices.size(); ++i)
2163  {
2164  AssertIndexRange(indices[i], size());
2165  entries[i] = (*this)[indices[i]];
2166  }
2167 }
2168 
2169 
2170 
2171 template <typename VectorType>
2172 template <typename ForwardIterator, typename OutputIterator>
2173 inline void
2175  ForwardIterator indices_begin,
2176  const ForwardIterator indices_end,
2177  OutputIterator values_begin) const
2178 {
2179  while (indices_begin != indices_end)
2180  {
2181  *values_begin = operator()(*indices_begin);
2182  indices_begin++;
2183  values_begin++;
2184  }
2185 }
2186 
2187 #endif // DOXYGEN
2188 
2190 
2191 #endif
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
std::size_t size() const
Definition: array_view.h:573
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()=default
BlockVectorBase & operator/=(const value_type factor)
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< OtherNumber > &values) const
const BlockIndices & get_block_indices() 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)
BlockVectorBase & operator+=(const BlockVectorBase &V)
::internal::BlockVectorIterators::Iterator< BlockVectorBase, false > iterator
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
BlockVectorBase & operator=(const value_type s)
value_type operator*(const BlockVectorBase &V) const
void compress(VectorOperation::values operation)
bool all_zero() const
void collect_sizes()
void sadd(const value_type s, const BlockVectorBase &V)
BlockVectorBase & operator*=(const value_type factor)
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)
BlockVectorBase & operator-=(const BlockVectorBase &V)
real_type l2_norm() const
value_type operator()(const size_type i) const
IndexSet locally_owned_elements() const
bool is_non_negative() const
BlockType & block(const unsigned int i)
bool operator==(const BlockVectorBase< VectorType2 > &v) const
std::size_t locally_owned_size() const
value_type * pointer
void equ(const value_type a, const BlockVector2 &V)
Definition: vector.h:110
bool operator==(const Iterator< BlockVectorType, OtherConstness > &i) const
Iterator & operator-=(const difference_type &d)
dereference_type operator[](const difference_type d) const
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
Iterator & operator+=(const difference_type &d)
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
Iterator & operator=(const Iterator &c)
typename std::conditional< Constness, const typename BlockVectorType::value_type, typename BlockVectorType::value_type >::type value_type
bool operator<=(const Iterator< BlockVectorType, OtherConstness > &i) const
typename BlockVectorType::reference reference
Iterator operator+(const difference_type &d) const
bool operator!=(const Iterator< BlockVectorType, OtherConstness > &i) const
bool operator>(const Iterator< BlockVectorType, OtherConstness > &i) const
std::enable_if_t< std::is_floating_point_v< T > &&std::is_floating_point_v< U >, typename ProductType< std::complex< T >, std::complex< U > >::type > operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcPointerToDifferentVectors()
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define AssertIsFinite(number)
Definition: exceptions.h:1884
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:490
static ::ExceptionBase & ExcDifferentBlockIndices()
static const bool value
Expression operator>=(const Expression &lhs, const Expression &rhs)
Expression operator>(const Expression &lhs, const Expression &rhs)
Expression operator<=(const Expression &lhs, const Expression &rhs)
static const char U
static const char N
static const char V
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)
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)
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:222
Definition: types.h:33
unsigned int global_dof_index
Definition: types.h:82
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
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)