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