Reference documentation for deal.II version GIT d2cc530c04 2023-03-22 15:10: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 - 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 
24 #include <deal.II/base/numbers.h>
26 
28 #include <deal.II/lac/exceptions.h>
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 
46 namespace 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 
79 template <typename VectorType>
81 {
82 public:
88  static const bool value = internal::is_block_vector<VectorType>;
89 };
90 
91 
92 // instantiation of the static member
93 template <typename VectorType>
95 
96 
97 
98 namespace internal
99 {
103  namespace BlockVectorIterators
104  {
120  template <class BlockVectorType, bool Constness>
121  class Iterator
122  {
123  public:
128 
134  using value_type =
135  typename std::conditional<Constness,
136  const typename BlockVectorType::value_type,
137  typename BlockVectorType::value_type>::type;
138 
145  using iterator_category = std::random_access_iterator_tag;
146  using difference_type = std::ptrdiff_t;
147  using reference = typename BlockVectorType::reference;
148  using pointer = value_type *;
149 
150  using dereference_type = typename std::conditional<
151  Constness,
152  value_type,
153  typename BlockVectorType::BlockType::reference>::type;
154 
159  using BlockVector = typename std::
160  conditional<Constness, const BlockVectorType, BlockVectorType>::type;
161 
171 
180 
181 
185  Iterator(const Iterator &c);
186 
187  private:
193  const size_type global_index,
194  const size_type current_block,
198 
199  public:
203  Iterator &
204  operator=(const Iterator &c);
205 
212  operator*() const;
213 
220 
225  Iterator &
227 
233  Iterator
235 
240  Iterator &
242 
248  Iterator
250 
255  template <bool OtherConstness>
256  bool
258 
263  template <bool OtherConstness>
264  bool
266 
272  template <bool OtherConstness>
273  bool
275 
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>
303 
308  Iterator
309  operator+(const difference_type &d) const;
310 
315  Iterator
316  operator-(const difference_type &d) const;
317 
322  Iterator &
324 
329  Iterator &
331 
342  "Your program tried to compare iterators pointing to "
343  "different block vectors. There is no reasonable way "
344  "to do this.");
345 
347  private:
354 
359 
364  unsigned int current_block;
366 
375 
379  void
381 
385  void
387 
388  // Mark all other instances of this template as friends.
389  template <typename, bool>
390  friend class Iterator;
391  };
392  } // namespace BlockVectorIterators
393 } // namespace internal
394 
395 
433 template <class VectorType>
435 {
436 public:
440  using BlockType = VectorType;
441 
442  /*
443  * Declare standard types used in
444  * all containers. These types
445  * parallel those in the
446  * <tt>C++</tt> standard
447  * libraries
448  * <tt>std::vector<...></tt>
449  * class. This includes iterator
450  * types.
451  */
452  using value_type = typename BlockType::value_type;
453  using pointer = value_type *;
454  using const_pointer = const value_type *;
455  using iterator =
459  using reference = typename BlockType::reference;
460  using const_reference = typename BlockType::const_reference;
462 
472  using real_type = typename BlockType::real_type;
473 
477  BlockVectorBase() = default;
478 
482  BlockVectorBase(const BlockVectorBase & /*V*/) = default;
483 
489  BlockVectorBase(BlockVectorBase && /*V*/) noexcept = default;
490 
497  void
499 
510  void
512 
516  BlockType &
517  block(const unsigned int i);
518 
522  const BlockType &
523  block(const unsigned int i) const;
524 
530  const BlockIndices &
532 
536  unsigned int
537  n_blocks() const;
538 
543  std::size_t
544  size() const;
545 
551  std::size_t
553 
570  IndexSet
572 
576  iterator
577  begin();
578 
584  begin() const;
585 
589  iterator
590  end();
591 
597  end() const;
598 
602  value_type
603  operator()(const size_type i) const;
604 
608  reference
609  operator()(const size_type i);
610 
616  value_type
617  operator[](const size_type i) const;
618 
624  reference
625  operator[](const size_type i);
626 
642  template <typename OtherNumber>
643  void
644  extract_subvector_to(const std::vector<size_type> &indices,
645  std::vector<OtherNumber> & values) const;
646 
674  template <typename ForwardIterator, typename OutputIterator>
675  void
676  extract_subvector_to(ForwardIterator indices_begin,
677  const ForwardIterator indices_end,
678  OutputIterator values_begin) const;
679 
685  operator=(const value_type s);
686 
691  operator=(const BlockVectorBase &V);
692 
699  operator=(BlockVectorBase && /*V*/) = default; // NOLINT
700 
704  template <class VectorType2>
706  operator=(const BlockVectorBase<VectorType2> &V);
707 
712  operator=(const VectorType &v);
713 
718  template <class VectorType2>
719  bool
720  operator==(const BlockVectorBase<VectorType2> &v) const;
721 
725  value_type
726  operator*(const BlockVectorBase &V) const;
727 
731  real_type
732  norm_sqr() const;
733 
737  value_type
738  mean_value() const;
739 
743  real_type
744  l1_norm() const;
745 
750  real_type
751  l2_norm() const;
752 
757  real_type
758  linfty_norm() const;
759 
782  value_type
784  const BlockVectorBase &V,
785  const BlockVectorBase &W);
786 
791  bool
793 
799  bool
800  all_zero() const;
801 
807  bool
809 
814  operator+=(const BlockVectorBase &V);
815 
820  operator-=(const BlockVectorBase &V);
821 
822 
827  template <typename Number>
828  void
829  add(const std::vector<size_type> &indices, const std::vector<Number> &values);
830 
835  template <typename Number>
836  void
837  add(const std::vector<size_type> &indices, const Vector<Number> &values);
838 
844  template <typename Number>
845  void
846  add(const size_type n_elements,
847  const size_type *indices,
848  const Number * values);
849 
854  void
855  add(const value_type s);
856 
860  void
861  add(const value_type a, const BlockVectorBase &V);
862 
866  void
867  add(const value_type a,
868  const BlockVectorBase &V,
869  const value_type b,
870  const BlockVectorBase &W);
871 
875  void
876  sadd(const value_type s, const BlockVectorBase &V);
877 
881  void
882  sadd(const value_type s, const value_type a, const BlockVectorBase &V);
883 
887  void
888  sadd(const value_type s,
889  const value_type a,
890  const BlockVectorBase &V,
891  const value_type b,
892  const BlockVectorBase &W);
893 
897  void
898  sadd(const value_type s,
899  const value_type a,
900  const BlockVectorBase &V,
901  const value_type b,
902  const BlockVectorBase &W,
903  const value_type c,
904  const BlockVectorBase &X);
905 
910  operator*=(const value_type factor);
911 
916  operator/=(const value_type factor);
917 
922  template <class BlockVector2>
923  void
924  scale(const BlockVector2 &v);
925 
929  template <class BlockVector2>
930  void
931  equ(const value_type a, const BlockVector2 &V);
932 
937  void
939 
944  std::size_t
946 
947 protected:
951  std::vector<VectorType> components;
952 
958 
959  // Make the iterator class a friend.
960  template <typename N, bool C>
961  friend class ::internal::BlockVectorIterators::Iterator;
962 
963  template <typename>
964  friend class BlockVectorBase;
965 };
966 
967 
970 /*----------------------- Inline functions ----------------------------------*/
971 
972 
973 #ifndef DOXYGEN
974 namespace internal
975 {
976  namespace BlockVectorIterators
977  {
978  template <class BlockVectorType, bool Constness>
979  inline Iterator<BlockVectorType, Constness>::Iterator(
980  const Iterator<BlockVectorType, Constness> &c)
981  : parent(c.parent)
982  , global_index(c.global_index)
983  , current_block(c.current_block)
984  , index_within_block(c.index_within_block)
985  , next_break_forward(c.next_break_forward)
986  , next_break_backward(c.next_break_backward)
987  {}
988 
989 
990 
991  template <class 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  // Only permit copy-constructing const iterators from non-const
1002  // iterators, and not vice versa (i.e., Constness must always be
1003  // true).
1004  static_assert(Constness == true,
1005  "Constructing a non-const iterator from a const iterator "
1006  "does not make sense.");
1007  }
1008 
1009 
1010 
1011  template <class BlockVectorType, bool Constness>
1012  inline Iterator<BlockVectorType, Constness>::Iterator(
1013  BlockVector & parent,
1014  const size_type global_index,
1015  const size_type current_block,
1016  const size_type index_within_block,
1017  const size_type next_break_forward,
1018  const size_type next_break_backward)
1019  : parent(&parent)
1021  , current_block(current_block)
1022  , index_within_block(index_within_block)
1023  , next_break_forward(next_break_forward)
1024  , next_break_backward(next_break_backward)
1025  {}
1026 
1027 
1028 
1029  template <class BlockVectorType, bool Constness>
1030  inline Iterator<BlockVectorType, Constness> &
1031  Iterator<BlockVectorType, Constness>::operator=(const Iterator &c)
1032  {
1033  parent = c.parent;
1034  global_index = c.global_index;
1035  index_within_block = c.index_within_block;
1036  current_block = c.current_block;
1037  next_break_forward = c.next_break_forward;
1038  next_break_backward = c.next_break_backward;
1039 
1040  return *this;
1041  }
1042 
1043 
1044 
1045  template <class BlockVectorType, bool Constness>
1046  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1048  {
1049  return parent->block(current_block)(index_within_block);
1050  }
1051 
1052 
1053 
1054  template <class BlockVectorType, bool Constness>
1055  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1056  Iterator<BlockVectorType, Constness>::operator[](
1057  const difference_type d) const
1058  {
1059  // if the index pointed to is
1060  // still within the block we
1061  // currently point into, then we
1062  // can save the computation of
1063  // the block
1064  if ((global_index + d >= next_break_backward) &&
1065  (global_index + d <= next_break_forward))
1066  return parent->block(current_block)(index_within_block + d);
1067 
1068  // if the index is not within the
1069  // block of the block vector into
1070  // which we presently point, then
1071  // there is no way: we have to
1072  // search for the block. this can
1073  // be done through the parent
1074  // class as well.
1075  return (*parent)(global_index + d);
1076  }
1077 
1078 
1079 
1080  template <class BlockVectorType, bool Constness>
1081  inline Iterator<BlockVectorType, Constness> &
1083  {
1084  move_forward();
1085  return *this;
1086  }
1087 
1088 
1089 
1090  template <class BlockVectorType, bool Constness>
1091  inline Iterator<BlockVectorType, Constness>
1093  {
1094  const Iterator old_value = *this;
1095  move_forward();
1096  return old_value;
1097  }
1098 
1099 
1100 
1101  template <class BlockVectorType, bool Constness>
1102  inline Iterator<BlockVectorType, Constness> &
1103  Iterator<BlockVectorType, Constness>::operator--()
1104  {
1105  move_backward();
1106  return *this;
1107  }
1108 
1109 
1110 
1111  template <class BlockVectorType, bool Constness>
1112  inline Iterator<BlockVectorType, Constness>
1113  Iterator<BlockVectorType, Constness>::operator--(int)
1114  {
1115  const Iterator old_value = *this;
1116  move_backward();
1117  return old_value;
1118  }
1119 
1120 
1121 
1122  template <class BlockVectorType, bool Constness>
1123  template <bool OtherConstness>
1124  inline bool
1126  const Iterator<BlockVectorType, OtherConstness> &i) const
1127  {
1128  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1129 
1130  return (global_index == i.global_index);
1131  }
1132 
1133 
1134 
1135  template <class 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 <class 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 <class 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 <class 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 <class 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 <class BlockVectorType, bool Constness>
1201  template <bool OtherConstness>
1202  inline typename Iterator<BlockVectorType, Constness>::difference_type
1204  const Iterator<BlockVectorType, OtherConstness> &i) const
1205  {
1206  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1207 
1208  return (static_cast<signed int>(global_index) -
1209  static_cast<signed int>(i.global_index));
1210  }
1211 
1212 
1213 
1214  template <class BlockVectorType, bool Constness>
1215  inline Iterator<BlockVectorType, Constness>
1217  const difference_type &d) const
1218  {
1219  // if the index pointed to is
1220  // still within the block we
1221  // currently point into, then we
1222  // can save the computation of
1223  // the block
1224  if ((global_index + d >= next_break_backward) &&
1225  (global_index + d <= next_break_forward))
1226  return Iterator(*parent,
1227  global_index + d,
1228  current_block,
1229  index_within_block + d,
1230  next_break_forward,
1231  next_break_backward);
1232  else
1233  // outside present block, so
1234  // have to seek new block
1235  // anyway
1236  return Iterator(*parent, global_index + d);
1237  }
1238 
1239 
1240 
1241  template <class BlockVectorType, bool Constness>
1242  inline Iterator<BlockVectorType, Constness>
1244  const difference_type &d) const
1245  {
1246  // if the index pointed to is
1247  // still within the block we
1248  // currently point into, then we
1249  // can save the computation of
1250  // the block
1251  if ((global_index - d >= next_break_backward) &&
1252  (global_index - d <= next_break_forward))
1253  return Iterator(*parent,
1254  global_index - d,
1255  current_block,
1256  index_within_block - d,
1257  next_break_forward,
1258  next_break_backward);
1259  else
1260  // outside present block, so
1261  // have to seek new block
1262  // anyway
1263  return Iterator(*parent, global_index - d);
1264  }
1265 
1266 
1267 
1268  template <class BlockVectorType, bool Constness>
1269  inline Iterator<BlockVectorType, Constness> &
1270  Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1271  {
1272  // if the index pointed to is
1273  // still within the block we
1274  // currently point into, then we
1275  // can save the computation of
1276  // the block
1277  if ((global_index + d >= next_break_backward) &&
1278  (global_index + d <= next_break_forward))
1279  {
1280  global_index += d;
1281  index_within_block += d;
1282  }
1283  else
1284  // outside present block, so
1285  // have to seek new block
1286  // anyway
1287  *this = Iterator(*parent, global_index + d);
1288 
1289  return *this;
1290  }
1291 
1292 
1293 
1294  template <class BlockVectorType, bool Constness>
1295  inline Iterator<BlockVectorType, Constness> &
1296  Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1297  {
1298  // if the index pointed to is
1299  // still within the block we
1300  // currently point into, then we
1301  // can save the computation of
1302  // the block
1303  if ((global_index - d >= next_break_backward) &&
1304  (global_index - d <= next_break_forward))
1305  {
1306  global_index -= d;
1307  index_within_block -= d;
1308  }
1309  else
1310  // outside present block, so
1311  // have to seek new block
1312  // anyway
1313  *this = Iterator(*parent, global_index - d);
1314 
1315  return *this;
1316  }
1317 
1318 
1319  template <class BlockVectorType, bool Constness>
1320  Iterator<BlockVectorType, Constness>::Iterator(BlockVector & parent,
1321  const size_type global_index)
1322  : parent(&parent)
1324  {
1325  // find which block we are
1326  // in. for this, take into
1327  // account that it happens at
1328  // times that people want to
1329  // initialize iterators
1330  // past-the-end
1331  if (global_index < parent.size())
1332  {
1333  const std::pair<size_type, size_type> indices =
1335  current_block = indices.first;
1336  index_within_block = indices.second;
1337 
1338  next_break_backward =
1339  parent.block_indices.local_to_global(current_block, 0);
1340  next_break_forward =
1341  (parent.block_indices.local_to_global(current_block, 0) +
1342  parent.block_indices.block_size(current_block) - 1);
1343  }
1344  else
1345  // past the end. only have one
1346  // value for this
1347  {
1348  this->global_index = parent.size();
1349  current_block = parent.n_blocks();
1350  index_within_block = 0;
1351  next_break_backward = global_index;
1352  next_break_forward = numbers::invalid_size_type;
1353  };
1354  }
1355 
1356 
1357 
1358  template <class BlockVectorType, bool Constness>
1359  void
1360  Iterator<BlockVectorType, Constness>::move_forward()
1361  {
1362  if (global_index != next_break_forward)
1363  ++index_within_block;
1364  else
1365  {
1366  // ok, we traverse a boundary
1367  // between blocks:
1368  index_within_block = 0;
1369  ++current_block;
1370 
1371  // break backwards is now old
1372  // break forward
1373  next_break_backward = next_break_forward + 1;
1374 
1375  // compute new break forward
1376  if (current_block < parent->block_indices.size())
1377  next_break_forward +=
1378  parent->block_indices.block_size(current_block);
1379  else
1380  // if we are beyond the end,
1381  // then move the next
1382  // boundary arbitrarily far
1383  // away
1384  next_break_forward = numbers::invalid_size_type;
1385  };
1386 
1387  ++global_index;
1388  }
1389 
1390 
1391 
1392  template <class BlockVectorType, bool Constness>
1393  void
1394  Iterator<BlockVectorType, Constness>::move_backward()
1395  {
1396  if (global_index != next_break_backward)
1397  --index_within_block;
1398  else if (current_block != 0)
1399  {
1400  // ok, we traverse a boundary
1401  // between blocks:
1402  --current_block;
1403  index_within_block =
1404  parent->block_indices.block_size(current_block) - 1;
1405 
1406  // break forwards is now old
1407  // break backward
1408  next_break_forward = next_break_backward - 1;
1409 
1410  // compute new break forward
1411  next_break_backward -=
1412  parent->block_indices.block_size(current_block);
1413  }
1414  else
1415  // current block was 0, we now
1416  // get into unspecified terrain
1417  {
1418  --current_block;
1419  index_within_block = numbers::invalid_size_type;
1420  next_break_forward = 0;
1421  next_break_backward = 0;
1422  };
1423 
1424  --global_index;
1425  }
1426 
1427 
1428  } // namespace BlockVectorIterators
1429 
1430 } // namespace internal
1431 
1432 
1433 
1434 template <class VectorType>
1435 inline std::size_t
1437 {
1438  return block_indices.total_size();
1439 }
1440 
1441 
1442 
1443 template <class VectorType>
1444 inline std::size_t
1446 {
1447  std::size_t local_size = 0;
1448  for (unsigned int b = 0; b < n_blocks(); ++b)
1449  local_size += block(b).locally_owned_size();
1450  return local_size;
1451 }
1452 
1453 
1454 
1455 template <class VectorType>
1456 inline IndexSet
1458 {
1459  IndexSet is(size());
1460 
1461  // copy index sets from blocks into the global one, shifted
1462  // by the appropriate amount for each block
1463  for (unsigned int b = 0; b < n_blocks(); ++b)
1464  {
1465  IndexSet x = block(b).locally_owned_elements();
1466  is.add_indices(x, block_indices.block_start(b));
1467  }
1468 
1469  is.compress();
1470 
1471  return is;
1472 }
1473 
1474 
1475 
1476 template <class VectorType>
1477 inline unsigned int
1479 {
1480  return block_indices.size();
1481 }
1482 
1483 
1484 template <class VectorType>
1486 BlockVectorBase<VectorType>::block(const unsigned int i)
1487 {
1488  AssertIndexRange(i, n_blocks());
1489 
1490  return components[i];
1491 }
1492 
1493 
1494 
1495 template <class VectorType>
1496 inline const typename BlockVectorBase<VectorType>::BlockType &
1497 BlockVectorBase<VectorType>::block(const unsigned int i) const
1498 {
1499  AssertIndexRange(i, n_blocks());
1500 
1501  return components[i];
1502 }
1503 
1504 
1505 
1506 template <class VectorType>
1507 inline const BlockIndices &
1509 {
1510  return block_indices;
1511 }
1512 
1513 
1514 template <class VectorType>
1515 inline void
1517 {
1518  std::vector<size_type> sizes(n_blocks());
1519 
1520  for (size_type i = 0; i < n_blocks(); ++i)
1521  sizes[i] = block(i).size();
1522 
1523  block_indices.reinit(sizes);
1524 }
1525 
1526 
1527 
1528 template <class VectorType>
1529 inline void
1531  ::VectorOperation::values operation)
1532 {
1533  for (unsigned int i = 0; i < n_blocks(); ++i)
1534  block(i).compress(operation);
1535 }
1536 
1537 
1538 
1539 template <class VectorType>
1542 {
1543  return iterator(*this, 0U);
1544 }
1545 
1546 
1547 
1548 template <class VectorType>
1551 {
1552  return const_iterator(*this, 0U);
1553 }
1554 
1555 
1556 template <class VectorType>
1559 {
1560  return iterator(*this, size());
1561 }
1562 
1563 
1564 
1565 template <class VectorType>
1568 {
1569  return const_iterator(*this, size());
1570 }
1571 
1572 
1573 template <class VectorType>
1574 inline bool
1576 {
1577  const std::pair<size_type, size_type> local_index =
1579 
1580  return components[local_index.first].in_local_range(global_index);
1581 }
1582 
1583 
1584 template <class VectorType>
1585 bool
1587 {
1588  for (size_type i = 0; i < n_blocks(); ++i)
1589  if (components[i].all_zero() == false)
1590  return false;
1591 
1592  return true;
1593 }
1594 
1595 
1596 
1597 template <class VectorType>
1598 bool
1600 {
1601  for (size_type i = 0; i < n_blocks(); ++i)
1602  if (components[i].is_non_negative() == false)
1603  return false;
1604 
1605  return true;
1606 }
1607 
1608 
1609 
1610 template <class VectorType>
1613  const BlockVectorBase<VectorType> &v) const
1614 {
1615  Assert(n_blocks() == v.n_blocks(),
1617 
1618  value_type sum = 0.;
1619  for (size_type i = 0; i < n_blocks(); ++i)
1620  sum += components[i] * v.components[i];
1621 
1622  return sum;
1623 }
1624 
1625 
1626 template <class VectorType>
1629 {
1630  real_type sum = 0.;
1631  for (size_type i = 0; i < n_blocks(); ++i)
1632  sum += components[i].norm_sqr();
1633 
1634  return sum;
1635 }
1636 
1637 
1638 
1639 template <class VectorType>
1642 {
1643  value_type sum = 0.;
1644  // need to do static_cast as otherwise it won't work with
1645  // value_type=complex<T>
1646  for (size_type i = 0; i < n_blocks(); ++i)
1647  sum += components[i].mean_value() *
1649  components[i].size()));
1650 
1651  return sum / (typename numbers::NumberTraits<value_type>::real_type(size()));
1652 }
1653 
1654 
1655 
1656 template <class VectorType>
1659 {
1660  real_type sum = 0.;
1661  for (size_type i = 0; i < n_blocks(); ++i)
1662  sum += components[i].l1_norm();
1663 
1664  return sum;
1665 }
1666 
1667 
1668 
1669 template <class VectorType>
1672 {
1673  return std::sqrt(norm_sqr());
1674 }
1675 
1676 
1677 
1678 template <class VectorType>
1681 {
1682  real_type sum = 0.;
1683  for (size_type i = 0; i < n_blocks(); ++i)
1684  {
1685  value_type newval = components[i].linfty_norm();
1686  if (sum < newval)
1687  sum = newval;
1688  }
1689  return sum;
1690 }
1691 
1692 
1693 
1694 template <class VectorType>
1697  const typename BlockVectorBase<VectorType>::value_type a,
1699  const BlockVectorBase<VectorType> & W)
1700 {
1701  AssertDimension(n_blocks(), V.n_blocks());
1703 
1704  value_type sum = 0.;
1705  for (size_type i = 0; i < n_blocks(); ++i)
1706  sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1707 
1708  return sum;
1709 }
1710 
1711 
1712 
1713 template <class VectorType>
1716 {
1717  Assert(n_blocks() == v.n_blocks(),
1719 
1720  for (size_type i = 0; i < n_blocks(); ++i)
1721  {
1722  components[i] += v.components[i];
1723  }
1724 
1725  return *this;
1726 }
1727 
1728 
1729 
1730 template <class VectorType>
1733 {
1734  Assert(n_blocks() == v.n_blocks(),
1736 
1737  for (size_type i = 0; i < n_blocks(); ++i)
1738  {
1739  components[i] -= v.components[i];
1740  }
1741  return *this;
1742 }
1743 
1744 
1745 
1746 template <class VectorType>
1747 template <typename Number>
1748 inline void
1749 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1750  const std::vector<Number> & values)
1751 {
1752  Assert(indices.size() == values.size(),
1753  ExcDimensionMismatch(indices.size(), values.size()));
1754  add(indices.size(), indices.data(), values.data());
1755 }
1756 
1757 
1758 
1759 template <class VectorType>
1760 template <typename Number>
1761 inline void
1762 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1763  const Vector<Number> & values)
1764 {
1765  Assert(indices.size() == values.size(),
1766  ExcDimensionMismatch(indices.size(), values.size()));
1767  const size_type n_indices = indices.size();
1768  for (size_type i = 0; i < n_indices; ++i)
1769  (*this)(indices[i]) += values(i);
1770 }
1771 
1772 
1773 
1774 template <class VectorType>
1775 template <typename Number>
1776 inline void
1778  const size_type *indices,
1779  const Number * values)
1780 {
1781  for (size_type i = 0; i < n_indices; ++i)
1782  (*this)(indices[i]) += values[i];
1783 }
1784 
1785 
1786 
1787 template <class VectorType>
1788 void
1790 {
1791  AssertIsFinite(a);
1792 
1793  for (size_type i = 0; i < n_blocks(); ++i)
1794  {
1795  components[i].add(a);
1796  }
1797 }
1798 
1799 
1800 
1801 template <class VectorType>
1802 void
1804  const BlockVectorBase<VectorType> &v)
1805 {
1806  AssertIsFinite(a);
1807 
1808  Assert(n_blocks() == v.n_blocks(),
1810 
1811  for (size_type i = 0; i < n_blocks(); ++i)
1812  {
1813  components[i].add(a, v.components[i]);
1814  }
1815 }
1816 
1817 
1818 
1819 template <class VectorType>
1820 void
1822  const BlockVectorBase<VectorType> &v,
1823  const value_type b,
1825 {
1826  AssertIsFinite(a);
1827  AssertIsFinite(b);
1828 
1829  Assert(n_blocks() == v.n_blocks(),
1831  Assert(n_blocks() == w.n_blocks(),
1832  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1833 
1834 
1835  for (size_type i = 0; i < n_blocks(); ++i)
1836  {
1837  components[i].add(a, v.components[i], b, w.components[i]);
1838  }
1839 }
1840 
1841 
1842 
1843 template <class VectorType>
1844 void
1846  const BlockVectorBase<VectorType> &v)
1847 {
1848  AssertIsFinite(x);
1849 
1850  Assert(n_blocks() == v.n_blocks(),
1852 
1853  for (size_type i = 0; i < n_blocks(); ++i)
1854  {
1855  components[i].sadd(x, v.components[i]);
1856  }
1857 }
1858 
1859 
1860 
1861 template <class VectorType>
1862 void
1864  const value_type a,
1865  const BlockVectorBase<VectorType> &v)
1866 {
1867  AssertIsFinite(x);
1868  AssertIsFinite(a);
1869 
1870  Assert(n_blocks() == v.n_blocks(),
1872 
1873  for (size_type i = 0; i < n_blocks(); ++i)
1874  {
1875  components[i].sadd(x, a, v.components[i]);
1876  }
1877 }
1878 
1879 
1880 
1881 template <class VectorType>
1882 void
1884  const value_type a,
1885  const BlockVectorBase<VectorType> &v,
1886  const value_type b,
1888 {
1889  AssertIsFinite(x);
1890  AssertIsFinite(a);
1891  AssertIsFinite(b);
1892 
1893  Assert(n_blocks() == v.n_blocks(),
1895  Assert(n_blocks() == w.n_blocks(),
1896  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1897 
1898  for (size_type i = 0; i < n_blocks(); ++i)
1899  {
1900  components[i].sadd(x, a, v.components[i], b, w.components[i]);
1901  }
1902 }
1903 
1904 
1905 
1906 template <class VectorType>
1907 void
1909  const value_type a,
1910  const BlockVectorBase<VectorType> &v,
1911  const value_type b,
1913  const value_type c,
1914  const BlockVectorBase<VectorType> &y)
1915 {
1916  AssertIsFinite(x);
1917  AssertIsFinite(a);
1918  AssertIsFinite(b);
1919  AssertIsFinite(c);
1920 
1921  Assert(n_blocks() == v.n_blocks(),
1923  Assert(n_blocks() == w.n_blocks(),
1924  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1925  Assert(n_blocks() == y.n_blocks(),
1927 
1928  for (size_type i = 0; i < n_blocks(); ++i)
1929  {
1930  components[i].sadd(
1931  x, a, v.components[i], b, w.components[i], c, y.components[i]);
1932  }
1933 }
1934 
1935 
1936 
1937 template <class VectorType>
1938 template <class BlockVector2>
1939 void
1940 BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
1941 {
1942  Assert(n_blocks() == v.n_blocks(),
1943  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1944  for (size_type i = 0; i < n_blocks(); ++i)
1945  components[i].scale(v.block(i));
1946 }
1947 
1948 
1949 
1950 template <class VectorType>
1951 std::size_t
1953 {
1954  return (MemoryConsumption::memory_consumption(this->block_indices) +
1955  MemoryConsumption::memory_consumption(this->components));
1956 }
1957 
1958 
1959 
1960 template <class VectorType>
1961 template <class BlockVector2>
1962 void
1963 BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
1964 {
1965  AssertIsFinite(a);
1966 
1967  Assert(n_blocks() == v.n_blocks(),
1968  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1969 
1970  for (size_type i = 0; i < n_blocks(); ++i)
1971  components[i].equ(a, v.components[i]);
1972 }
1973 
1974 
1975 
1976 template <class VectorType>
1977 void
1979 {
1980  for (size_type i = 0; i < n_blocks(); ++i)
1981  block(i).update_ghost_values();
1982 }
1983 
1984 
1985 
1986 template <class VectorType>
1989 {
1990  AssertIsFinite(s);
1991 
1992  for (size_type i = 0; i < n_blocks(); ++i)
1993  components[i] = s;
1994 
1995  return *this;
1996 }
1997 
1998 
1999 template <class VectorType>
2002 {
2004 
2005  for (size_type i = 0; i < n_blocks(); ++i)
2006  components[i] = v.components[i];
2007 
2008  return *this;
2009 }
2010 
2011 
2012 template <class VectorType>
2013 template <class VectorType2>
2016 {
2018 
2019  for (size_type i = 0; i < n_blocks(); ++i)
2020  components[i] = v.components[i];
2021 
2022  return *this;
2023 }
2024 
2025 
2026 
2027 template <class VectorType>
2029 BlockVectorBase<VectorType>::operator=(const VectorType &v)
2030 {
2031  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2032 
2033  size_type index_v = 0;
2034  for (size_type b = 0; b < n_blocks(); ++b)
2035  for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2036  block(b)(i) = v(index_v);
2037 
2038  return *this;
2039 }
2040 
2041 
2042 
2043 template <class VectorType>
2044 template <class VectorType2>
2045 inline bool
2047  const BlockVectorBase<VectorType2> &v) const
2048 {
2050 
2051  for (size_type i = 0; i < n_blocks(); ++i)
2052  if (!(components[i] == v.components[i]))
2053  return false;
2054 
2055  return true;
2056 }
2057 
2058 
2059 
2060 template <class VectorType>
2063 {
2064  AssertIsFinite(factor);
2065 
2066  for (size_type i = 0; i < n_blocks(); ++i)
2067  components[i] *= factor;
2068 
2069  return *this;
2070 }
2071 
2072 
2073 
2074 template <class VectorType>
2077 {
2078  AssertIsFinite(factor);
2079  Assert(factor != 0., ExcDivideByZero());
2080 
2081  const value_type inverse_factor = value_type(1.) / factor;
2082 
2083  for (size_type i = 0; i < n_blocks(); ++i)
2084  components[i] *= inverse_factor;
2085 
2086  return *this;
2087 }
2088 
2089 
2090 template <class VectorType>
2093 {
2094  const std::pair<unsigned int, size_type> local_index =
2096  return components[local_index.first](local_index.second);
2097 }
2098 
2099 
2100 
2101 template <class VectorType>
2104 {
2105  const std::pair<unsigned int, size_type> local_index =
2107  return components[local_index.first](local_index.second);
2108 }
2109 
2110 
2111 
2112 template <class VectorType>
2115 {
2116  return operator()(i);
2117 }
2118 
2119 
2120 
2121 template <class VectorType>
2124 {
2125  return operator()(i);
2126 }
2127 
2128 
2129 
2130 template <typename VectorType>
2131 template <typename OtherNumber>
2132 inline void
2134  const std::vector<size_type> &indices,
2135  std::vector<OtherNumber> & values) const
2136 {
2137  for (size_type i = 0; i < indices.size(); ++i)
2138  values[i] = operator()(indices[i]);
2139 }
2140 
2141 
2142 
2143 template <typename VectorType>
2144 template <typename ForwardIterator, typename OutputIterator>
2145 inline void
2147  ForwardIterator indices_begin,
2148  const ForwardIterator indices_end,
2149  OutputIterator values_begin) const
2150 {
2151  while (indices_begin != indices_end)
2152  {
2153  *values_begin = operator()(*indices_begin);
2154  indices_begin++;
2155  values_begin++;
2156  }
2157 }
2158 
2159 #endif // DOXYGEN
2160 
2162 
2163 #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
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
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
bool all_zero() const
std::size_t size() 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)
void compress(::VectorOperation::values operation)
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:109
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< T >::value &&std::is_floating_point< U >::value, 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:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcPointerToDifferentVectors()
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
#define AssertIsFinite(number)
Definition: exceptions.h:1854
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:488
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< T >::value, 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)
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
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)