Reference documentation for deal.II version Git d2841aca35 2021-05-14 14:24:54 -0600
\(\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 - 2021 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 // Forward declaration
46 #ifndef DOXYGEN
47 template <typename>
48 class BlockVectorBase;
49 #endif
50 
65 template <typename VectorType>
67 {
68 private:
73  template <typename T>
74  static std::true_type
76 
81  static std::false_type
83 
84 public:
90  static const bool value =
91  std::is_same<decltype(check_for_block_vector(std::declval<VectorType *>())),
92  std::true_type>::value;
93 };
94 
95 
96 // instantiation of the static member
97 template <typename VectorType>
99 
100 
101 
102 namespace internal
103 {
107  namespace BlockVectorIterators
108  {
124  template <class BlockVectorType, bool Constness>
125  class Iterator
126  {
127  public:
132 
138  using value_type =
139  typename std::conditional<Constness,
140  const typename BlockVectorType::value_type,
141  typename BlockVectorType::value_type>::type;
142 
148  using iterator_category = std::random_access_iterator_tag;
149  using difference_type = std::ptrdiff_t;
150  using reference = typename BlockVectorType::reference;
151  using pointer = value_type *;
152 
153  using dereference_type = typename std::conditional<
154  Constness,
155  value_type,
156  typename BlockVectorType::BlockType::reference>::type;
157 
162  using BlockVector = typename std::
163  conditional<Constness, const BlockVectorType, BlockVectorType>::type;
164 
173  Iterator(BlockVector &parent, const size_type global_index);
174 
183 
184 
188  Iterator(const Iterator &c);
189 
190  private:
195  Iterator(BlockVector & parent,
196  const size_type global_index,
197  const size_type current_block,
198  const size_type index_within_block,
199  const size_type next_break_forward,
200  const size_type next_break_backward);
201 
202  public:
206  Iterator &
207  operator=(const Iterator &c);
208 
214  dereference_type operator*() const;
215 
220  dereference_type operator[](const difference_type d) const;
221 
226  Iterator &
227  operator++();
228 
234  Iterator
235  operator++(int);
236 
241  Iterator &
242  operator--();
243 
249  Iterator
250  operator--(int);
251 
256  template <bool OtherConstness>
257  bool
259 
264  template <bool OtherConstness>
265  bool
267 
273  template <bool OtherConstness>
274  bool
275  operator<(const Iterator<BlockVectorType, OtherConstness> &i) const;
276 
280  template <bool OtherConstness>
281  bool
282  operator<=(const Iterator<BlockVectorType, OtherConstness> &i) const;
283 
287  template <bool OtherConstness>
288  bool
290 
294  template <bool OtherConstness>
295  bool
297 
301  template <bool OtherConstness>
304 
309  Iterator
310  operator+(const difference_type &d) const;
311 
316  Iterator
317  operator-(const difference_type &d) const;
318 
323  Iterator &
324  operator+=(const difference_type &d);
325 
330  Iterator &
331  operator-=(const difference_type &d);
332 
342  DeclExceptionMsg(ExcPointerToDifferentVectors,
343  "Your program tried to compare iterators pointing to "
344  "different block vectors. There is no reasonable way "
345  "to do this.");
346 
348  private:
355 
360 
365  unsigned int current_block;
367 
376 
380  void
381  move_forward();
382 
386  void
387  move_backward();
388 
389  // Mark all other instances of this template as friends.
390  template <typename, bool>
391  friend class Iterator;
392  };
393  } // namespace BlockVectorIterators
394 } // namespace internal
395 
396 
434 template <class VectorType>
436 {
437 public:
442 
443  /*
444  * Declare standard types used in
445  * all containers. These types
446  * parallel those in the
447  * <tt>C++</tt> standard
448  * libraries
449  * <tt>std::vector<...></tt>
450  * class. This includes iterator
451  * types.
452  */
454  using pointer = value_type *;
455  using const_pointer = const value_type *;
456  using iterator =
458  using const_iterator =
460  using reference = typename BlockType::reference;
463 
473  using real_type = typename BlockType::real_type;
474 
478  BlockVectorBase() = default;
479 
483  BlockVectorBase(const BlockVectorBase & /*V*/) = default;
484 
490  BlockVectorBase(BlockVectorBase && /*V*/) noexcept = default;
491 
498  void
499  collect_sizes();
500 
511  void
512  compress(::VectorOperation::values operation);
513 
517  BlockType &
518  block(const unsigned int i);
519 
523  const BlockType &
524  block(const unsigned int i) const;
525 
531  const BlockIndices &
532  get_block_indices() const;
533 
537  unsigned int
538  n_blocks() const;
539 
544  std::size_t
545  size() const;
546 
552  std::size_t
553  locally_owned_size() const;
554 
571  IndexSet
572  locally_owned_elements() const;
573 
577  iterator
578  begin();
579 
585  begin() const;
586 
590  iterator
591  end();
592 
598  end() const;
599 
603  value_type
604  operator()(const size_type i) const;
605 
609  reference
610  operator()(const size_type i);
611 
617  value_type operator[](const size_type i) const;
618 
624  reference 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 operator*(const BlockVectorBase &V) const;
725 
729  real_type
730  norm_sqr() const;
731 
735  value_type
736  mean_value() const;
737 
741  real_type
742  l1_norm() const;
743 
748  real_type
749  l2_norm() const;
750 
755  real_type
756  linfty_norm() const;
757 
780  value_type
781  add_and_dot(const value_type a,
782  const BlockVectorBase &V,
783  const BlockVectorBase &W);
784 
789  bool
790  in_local_range(const size_type global_index) const;
791 
797  bool
798  all_zero() const;
799 
805  bool
806  is_non_negative() const;
807 
812  operator+=(const BlockVectorBase &V);
813 
818  operator-=(const BlockVectorBase &V);
819 
820 
825  template <typename Number>
826  void
827  add(const std::vector<size_type> &indices, const std::vector<Number> &values);
828 
833  template <typename Number>
834  void
835  add(const std::vector<size_type> &indices, const Vector<Number> &values);
836 
842  template <typename Number>
843  void
844  add(const size_type n_elements,
845  const size_type *indices,
846  const Number * values);
847 
852  void
853  add(const value_type s);
854 
858  void
859  add(const value_type a, const BlockVectorBase &V);
860 
864  void
865  add(const value_type a,
866  const BlockVectorBase &V,
867  const value_type b,
868  const BlockVectorBase &W);
869 
873  void
874  sadd(const value_type s, const BlockVectorBase &V);
875 
879  void
880  sadd(const value_type s, const value_type a, const BlockVectorBase &V);
881 
885  void
886  sadd(const value_type s,
887  const value_type a,
888  const BlockVectorBase &V,
889  const value_type b,
890  const BlockVectorBase &W);
891 
895  void
896  sadd(const value_type s,
897  const value_type a,
898  const BlockVectorBase &V,
899  const value_type b,
900  const BlockVectorBase &W,
901  const value_type c,
902  const BlockVectorBase &X);
903 
908  operator*=(const value_type factor);
909 
914  operator/=(const value_type factor);
915 
920  template <class BlockVector2>
921  void
922  scale(const BlockVector2 &v);
923 
927  template <class BlockVector2>
928  void
929  equ(const value_type a, const BlockVector2 &V);
930 
935  void
936  update_ghost_values() const;
937 
942  std::size_t
943  memory_consumption() const;
944 
945 protected:
949  std::vector<VectorType> components;
950 
955  BlockIndices block_indices;
956 
957  // Make the iterator class a friend.
958  template <typename N, bool C>
959  friend class ::internal::BlockVectorIterators::Iterator;
960 
961  template <typename>
962  friend class BlockVectorBase;
963 };
964 
965 
968 /*----------------------- Inline functions ----------------------------------*/
969 
970 
971 #ifndef DOXYGEN
972 namespace internal
973 {
974  namespace BlockVectorIterators
975  {
976  template <class BlockVectorType, bool Constness>
977  inline Iterator<BlockVectorType, Constness>::Iterator(
978  const Iterator<BlockVectorType, Constness> &c)
979  : parent(c.parent)
980  , global_index(c.global_index)
981  , current_block(c.current_block)
982  , index_within_block(c.index_within_block)
983  , next_break_forward(c.next_break_forward)
984  , next_break_backward(c.next_break_backward)
985  {}
986 
987 
988 
989  template <class BlockVectorType, bool Constness>
990  inline Iterator<BlockVectorType, Constness>::Iterator(
991  const Iterator<BlockVectorType, !Constness> &c)
992  : parent(c.parent)
993  , global_index(c.global_index)
994  , current_block(c.current_block)
995  , index_within_block(c.index_within_block)
996  , next_break_forward(c.next_break_forward)
997  , next_break_backward(c.next_break_backward)
998  {
999  // Only permit copy-constructing const iterators from non-const
1000  // iterators, and not vice versa (i.e., Constness must always be
1001  // true).
1002  static_assert(Constness == true,
1003  "Constructing a non-const iterator from a const iterator "
1004  "does not make sense.");
1005  }
1006 
1007 
1008 
1009  template <class BlockVectorType, bool Constness>
1010  inline Iterator<BlockVectorType, Constness>::Iterator(
1011  BlockVector & parent,
1012  const size_type global_index,
1013  const size_type current_block,
1014  const size_type index_within_block,
1015  const size_type next_break_forward,
1016  const size_type next_break_backward)
1017  : parent(&parent)
1018  , global_index(global_index)
1019  , current_block(current_block)
1020  , index_within_block(index_within_block)
1021  , next_break_forward(next_break_forward)
1022  , next_break_backward(next_break_backward)
1023  {}
1024 
1025 
1026 
1027  template <class BlockVectorType, bool Constness>
1028  inline Iterator<BlockVectorType, Constness> &
1029  Iterator<BlockVectorType, Constness>::operator=(const Iterator &c)
1030  {
1031  parent = c.parent;
1032  global_index = c.global_index;
1033  index_within_block = c.index_within_block;
1034  current_block = c.current_block;
1035  next_break_forward = c.next_break_forward;
1036  next_break_backward = c.next_break_backward;
1037 
1038  return *this;
1039  }
1040 
1041 
1042 
1043  template <class BlockVectorType, bool Constness>
1044  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1046  {
1047  return parent->block(current_block)(index_within_block);
1048  }
1049 
1050 
1051 
1052  template <class BlockVectorType, bool Constness>
1053  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1054  Iterator<BlockVectorType, Constness>::
1055  operator[](const difference_type d) const
1056  {
1057  // if the index pointed to is
1058  // still within the block we
1059  // currently point into, then we
1060  // can save the computation of
1061  // the block
1062  if ((global_index + d >= next_break_backward) &&
1063  (global_index + d <= next_break_forward))
1064  return parent->block(current_block)(index_within_block + d);
1065 
1066  // if the index is not within the
1067  // block of the block vector into
1068  // which we presently point, then
1069  // there is no way: we have to
1070  // search for the block. this can
1071  // be done through the parent
1072  // class as well.
1073  return (*parent)(global_index + d);
1074  }
1075 
1076 
1077 
1078  template <class BlockVectorType, bool Constness>
1079  inline Iterator<BlockVectorType, Constness> &
1081  {
1082  move_forward();
1083  return *this;
1084  }
1085 
1086 
1087 
1088  template <class BlockVectorType, bool Constness>
1089  inline Iterator<BlockVectorType, Constness>
1091  {
1092  const Iterator old_value = *this;
1093  move_forward();
1094  return old_value;
1095  }
1096 
1097 
1098 
1099  template <class BlockVectorType, bool Constness>
1100  inline Iterator<BlockVectorType, Constness> &
1101  Iterator<BlockVectorType, Constness>::operator--()
1102  {
1103  move_backward();
1104  return *this;
1105  }
1106 
1107 
1108 
1109  template <class BlockVectorType, bool Constness>
1110  inline Iterator<BlockVectorType, Constness>
1111  Iterator<BlockVectorType, Constness>::operator--(int)
1112  {
1113  const Iterator old_value = *this;
1114  move_backward();
1115  return old_value;
1116  }
1117 
1118 
1119 
1120  template <class BlockVectorType, bool Constness>
1121  template <bool OtherConstness>
1122  inline bool
1124  operator==(const Iterator<BlockVectorType, OtherConstness> &i) const
1125  {
1126  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1127 
1128  return (global_index == i.global_index);
1129  }
1130 
1131 
1132 
1133  template <class BlockVectorType, bool Constness>
1134  template <bool OtherConstness>
1135  inline bool
1137  operator!=(const Iterator<BlockVectorType, OtherConstness> &i) const
1138  {
1139  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1140 
1141  return (global_index != i.global_index);
1142  }
1143 
1144 
1145 
1146  template <class BlockVectorType, bool Constness>
1147  template <bool OtherConstness>
1148  inline bool
1150  operator<(const Iterator<BlockVectorType, OtherConstness> &i) const
1151  {
1152  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1153 
1154  return (global_index < i.global_index);
1155  }
1156 
1157 
1158 
1159  template <class BlockVectorType, bool Constness>
1160  template <bool OtherConstness>
1161  inline bool
1163  operator<=(const Iterator<BlockVectorType, OtherConstness> &i) const
1164  {
1165  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1166 
1167  return (global_index <= i.global_index);
1168  }
1169 
1170 
1171 
1172  template <class BlockVectorType, bool Constness>
1173  template <bool OtherConstness>
1174  inline bool
1176  operator>(const Iterator<BlockVectorType, OtherConstness> &i) const
1177  {
1178  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1179 
1180  return (global_index > i.global_index);
1181  }
1182 
1183 
1184 
1185  template <class BlockVectorType, bool Constness>
1186  template <bool OtherConstness>
1187  inline bool
1189  operator>=(const Iterator<BlockVectorType, OtherConstness> &i) const
1190  {
1191  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1192 
1193  return (global_index >= i.global_index);
1194  }
1195 
1196 
1197 
1198  template <class BlockVectorType, bool Constness>
1199  template <bool OtherConstness>
1200  inline typename Iterator<BlockVectorType, Constness>::difference_type
1202  operator-(const Iterator<BlockVectorType, OtherConstness> &i) const
1203  {
1204  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1205 
1206  return (static_cast<signed int>(global_index) -
1207  static_cast<signed int>(i.global_index));
1208  }
1209 
1210 
1211 
1212  template <class BlockVectorType, bool Constness>
1213  inline Iterator<BlockVectorType, Constness>
1215  operator+(const difference_type &d) const
1216  {
1217  // if the index pointed to is
1218  // still within the block we
1219  // currently point into, then we
1220  // can save the computation of
1221  // the block
1222  if ((global_index + d >= next_break_backward) &&
1223  (global_index + d <= next_break_forward))
1224  return Iterator(*parent,
1225  global_index + d,
1226  current_block,
1227  index_within_block + d,
1228  next_break_forward,
1229  next_break_backward);
1230  else
1231  // outside present block, so
1232  // have to seek new block
1233  // anyway
1234  return Iterator(*parent, global_index + d);
1235  }
1236 
1237 
1238 
1239  template <class BlockVectorType, bool Constness>
1240  inline Iterator<BlockVectorType, Constness>
1242  operator-(const difference_type &d) const
1243  {
1244  // if the index pointed to is
1245  // still within the block we
1246  // currently point into, then we
1247  // can save the computation of
1248  // the block
1249  if ((global_index - d >= next_break_backward) &&
1250  (global_index - d <= next_break_forward))
1251  return Iterator(*parent,
1252  global_index - d,
1253  current_block,
1254  index_within_block - d,
1255  next_break_forward,
1256  next_break_backward);
1257  else
1258  // outside present block, so
1259  // have to seek new block
1260  // anyway
1261  return Iterator(*parent, global_index - d);
1262  }
1263 
1264 
1265 
1266  template <class BlockVectorType, bool Constness>
1267  inline Iterator<BlockVectorType, Constness> &
1268  Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1269  {
1270  // if the index pointed to is
1271  // still within the block we
1272  // currently point into, then we
1273  // can save the computation of
1274  // the block
1275  if ((global_index + d >= next_break_backward) &&
1276  (global_index + d <= next_break_forward))
1277  {
1278  global_index += d;
1279  index_within_block += d;
1280  }
1281  else
1282  // outside present block, so
1283  // have to seek new block
1284  // anyway
1285  *this = Iterator(*parent, global_index + d);
1286 
1287  return *this;
1288  }
1289 
1290 
1291 
1292  template <class BlockVectorType, bool Constness>
1293  inline Iterator<BlockVectorType, Constness> &
1294  Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1295  {
1296  // if the index pointed to is
1297  // still within the block we
1298  // currently point into, then we
1299  // can save the computation of
1300  // the block
1301  if ((global_index - d >= next_break_backward) &&
1302  (global_index - d <= next_break_forward))
1303  {
1304  global_index -= d;
1305  index_within_block -= d;
1306  }
1307  else
1308  // outside present block, so
1309  // have to seek new block
1310  // anyway
1311  *this = Iterator(*parent, global_index - d);
1312 
1313  return *this;
1314  }
1315 
1316 
1317  template <class BlockVectorType, bool Constness>
1318  Iterator<BlockVectorType, Constness>::Iterator(BlockVector & parent,
1319  const size_type global_index)
1320  : parent(&parent)
1321  , global_index(global_index)
1322  {
1323  // find which block we are
1324  // in. for this, take into
1325  // account that it happens at
1326  // times that people want to
1327  // initialize iterators
1328  // past-the-end
1329  if (global_index < parent.size())
1330  {
1331  const std::pair<size_type, size_type> indices =
1332  parent.block_indices.global_to_local(global_index);
1333  current_block = indices.first;
1334  index_within_block = indices.second;
1335 
1336  next_break_backward =
1337  parent.block_indices.local_to_global(current_block, 0);
1338  next_break_forward =
1339  (parent.block_indices.local_to_global(current_block, 0) +
1340  parent.block_indices.block_size(current_block) - 1);
1341  }
1342  else
1343  // past the end. only have one
1344  // value for this
1345  {
1346  this->global_index = parent.size();
1347  current_block = parent.n_blocks();
1348  index_within_block = 0;
1349  next_break_backward = global_index;
1350  next_break_forward = numbers::invalid_size_type;
1351  };
1352  }
1353 
1354 
1355 
1356  template <class BlockVectorType, bool Constness>
1357  void
1358  Iterator<BlockVectorType, Constness>::move_forward()
1359  {
1360  if (global_index != next_break_forward)
1361  ++index_within_block;
1362  else
1363  {
1364  // ok, we traverse a boundary
1365  // between blocks:
1366  index_within_block = 0;
1367  ++current_block;
1368 
1369  // break backwards is now old
1370  // break forward
1371  next_break_backward = next_break_forward + 1;
1372 
1373  // compute new break forward
1374  if (current_block < parent->block_indices.size())
1375  next_break_forward +=
1376  parent->block_indices.block_size(current_block);
1377  else
1378  // if we are beyond the end,
1379  // then move the next
1380  // boundary arbitrarily far
1381  // away
1382  next_break_forward = numbers::invalid_size_type;
1383  };
1384 
1385  ++global_index;
1386  }
1387 
1388 
1389 
1390  template <class BlockVectorType, bool Constness>
1391  void
1392  Iterator<BlockVectorType, Constness>::move_backward()
1393  {
1394  if (global_index != next_break_backward)
1395  --index_within_block;
1396  else if (current_block != 0)
1397  {
1398  // ok, we traverse a boundary
1399  // between blocks:
1400  --current_block;
1401  index_within_block =
1402  parent->block_indices.block_size(current_block) - 1;
1403 
1404  // break forwards is now old
1405  // break backward
1406  next_break_forward = next_break_backward - 1;
1407 
1408  // compute new break forward
1409  next_break_backward -=
1410  parent->block_indices.block_size(current_block);
1411  }
1412  else
1413  // current block was 0, we now
1414  // get into unspecified terrain
1415  {
1416  --current_block;
1417  index_within_block = numbers::invalid_size_type;
1418  next_break_forward = 0;
1419  next_break_backward = 0;
1420  };
1421 
1422  --global_index;
1423  }
1424 
1425 
1426  } // namespace BlockVectorIterators
1427 
1428 } // namespace internal
1429 
1430 
1431 
1432 template <class VectorType>
1433 inline std::size_t
1435 {
1436  return block_indices.total_size();
1437 }
1438 
1439 
1440 
1441 template <class VectorType>
1442 inline std::size_t
1444 {
1445  std::size_t local_size = 0;
1446  for (unsigned int b = 0; b < n_blocks(); ++b)
1447  local_size += block(b).locally_owned_size();
1448  return local_size;
1449 }
1450 
1451 
1452 
1453 template <class VectorType>
1454 inline IndexSet
1456 {
1457  IndexSet is(size());
1458 
1459  // copy index sets from blocks into the global one, shifted
1460  // by the appropriate amount for each block
1461  for (unsigned int b = 0; b < n_blocks(); ++b)
1462  {
1463  IndexSet x = block(b).locally_owned_elements();
1464  is.add_indices(x, block_indices.block_start(b));
1465  }
1466 
1467  is.compress();
1468 
1469  return is;
1470 }
1471 
1472 
1473 
1474 template <class VectorType>
1475 inline unsigned int
1477 {
1478  return block_indices.size();
1479 }
1480 
1481 
1482 template <class VectorType>
1484 BlockVectorBase<VectorType>::block(const unsigned int i)
1485 {
1486  AssertIndexRange(i, n_blocks());
1487 
1488  return components[i];
1489 }
1490 
1491 
1492 
1493 template <class VectorType>
1494 inline const typename BlockVectorBase<VectorType>::BlockType &
1495 BlockVectorBase<VectorType>::block(const unsigned int i) const
1496 {
1497  AssertIndexRange(i, n_blocks());
1498 
1499  return components[i];
1500 }
1501 
1502 
1503 
1504 template <class VectorType>
1505 inline const BlockIndices &
1507 {
1508  return block_indices;
1509 }
1510 
1511 
1512 template <class VectorType>
1513 inline void
1515 {
1516  std::vector<size_type> sizes(n_blocks());
1517 
1518  for (size_type i = 0; i < n_blocks(); ++i)
1519  sizes[i] = block(i).size();
1520 
1521  block_indices.reinit(sizes);
1522 }
1523 
1524 
1525 
1526 template <class VectorType>
1527 inline void
1529  ::VectorOperation::values operation)
1530 {
1531  for (unsigned int i = 0; i < n_blocks(); ++i)
1532  block(i).compress(operation);
1533 }
1534 
1535 
1536 
1537 template <class VectorType>
1540 {
1541  return iterator(*this, 0U);
1542 }
1543 
1544 
1545 
1546 template <class VectorType>
1549 {
1550  return const_iterator(*this, 0U);
1551 }
1552 
1553 
1554 template <class VectorType>
1557 {
1558  return iterator(*this, size());
1559 }
1560 
1561 
1562 
1563 template <class VectorType>
1566 {
1567  return const_iterator(*this, size());
1568 }
1569 
1570 
1571 template <class VectorType>
1572 inline bool
1574 {
1575  const std::pair<size_type, size_type> local_index =
1576  block_indices.global_to_local(global_index);
1577 
1578  return components[local_index.first].in_local_range(global_index);
1579 }
1580 
1581 
1582 template <class VectorType>
1583 bool
1585 {
1586  for (size_type i = 0; i < n_blocks(); ++i)
1587  if (components[i].all_zero() == false)
1588  return false;
1589 
1590  return true;
1591 }
1592 
1593 
1594 
1595 template <class VectorType>
1596 bool
1598 {
1599  for (size_type i = 0; i < n_blocks(); ++i)
1600  if (components[i].is_non_negative() == false)
1601  return false;
1602 
1603  return true;
1604 }
1605 
1606 
1607 
1608 template <class VectorType>
1611 {
1612  Assert(n_blocks() == v.n_blocks(),
1614 
1615  value_type sum = 0.;
1616  for (size_type i = 0; i < n_blocks(); ++i)
1617  sum += components[i] * v.components[i];
1618 
1619  return sum;
1620 }
1621 
1622 
1623 template <class VectorType>
1626 {
1627  real_type sum = 0.;
1628  for (size_type i = 0; i < n_blocks(); ++i)
1629  sum += components[i].norm_sqr();
1630 
1631  return sum;
1632 }
1633 
1634 
1635 
1636 template <class VectorType>
1639 {
1640  value_type sum = 0.;
1641  // need to do static_cast as otherwise it won't work with
1642  // value_type=complex<T>
1643  for (size_type i = 0; i < n_blocks(); ++i)
1644  sum += components[i].mean_value() *
1646  components[i].size()));
1647 
1648  return sum / (typename numbers::NumberTraits<value_type>::real_type(size()));
1649 }
1650 
1651 
1652 
1653 template <class VectorType>
1656 {
1657  real_type sum = 0.;
1658  for (size_type i = 0; i < n_blocks(); ++i)
1659  sum += components[i].l1_norm();
1660 
1661  return sum;
1662 }
1663 
1664 
1665 
1666 template <class VectorType>
1669 {
1670  return std::sqrt(norm_sqr());
1671 }
1672 
1673 
1674 
1675 template <class VectorType>
1678 {
1679  real_type sum = 0.;
1680  for (size_type i = 0; i < n_blocks(); ++i)
1681  {
1682  value_type newval = components[i].linfty_norm();
1683  if (sum < newval)
1684  sum = newval;
1685  }
1686  return sum;
1687 }
1688 
1689 
1690 
1691 template <class VectorType>
1694  const typename BlockVectorBase<VectorType>::value_type a,
1695  const BlockVectorBase<VectorType> & V,
1696  const BlockVectorBase<VectorType> & W)
1697 {
1700 
1701  value_type sum = 0.;
1702  for (size_type i = 0; i < n_blocks(); ++i)
1703  sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1704 
1705  return sum;
1706 }
1707 
1708 
1709 
1710 template <class VectorType>
1713 {
1714  Assert(n_blocks() == v.n_blocks(),
1716 
1717  for (size_type i = 0; i < n_blocks(); ++i)
1718  {
1719  components[i] += v.components[i];
1720  }
1721 
1722  return *this;
1723 }
1724 
1725 
1726 
1727 template <class VectorType>
1730 {
1731  Assert(n_blocks() == v.n_blocks(),
1733 
1734  for (size_type i = 0; i < n_blocks(); ++i)
1735  {
1736  components[i] -= v.components[i];
1737  }
1738  return *this;
1739 }
1740 
1741 
1742 
1743 template <class VectorType>
1744 template <typename Number>
1745 inline void
1746 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1747  const std::vector<Number> & values)
1748 {
1749  Assert(indices.size() == values.size(),
1750  ExcDimensionMismatch(indices.size(), values.size()));
1751  add(indices.size(), indices.data(), values.data());
1752 }
1753 
1754 
1755 
1756 template <class VectorType>
1757 template <typename Number>
1758 inline void
1759 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1760  const Vector<Number> & values)
1761 {
1762  Assert(indices.size() == values.size(),
1763  ExcDimensionMismatch(indices.size(), values.size()));
1764  const size_type n_indices = indices.size();
1765  for (size_type i = 0; i < n_indices; ++i)
1766  (*this)(indices[i]) += values(i);
1767 }
1768 
1769 
1770 
1771 template <class VectorType>
1772 template <typename Number>
1773 inline void
1775  const size_type *indices,
1776  const Number * values)
1777 {
1778  for (size_type i = 0; i < n_indices; ++i)
1779  (*this)(indices[i]) += values[i];
1780 }
1781 
1782 
1783 
1784 template <class VectorType>
1785 void
1787 {
1788  AssertIsFinite(a);
1789 
1790  for (size_type i = 0; i < n_blocks(); ++i)
1791  {
1792  components[i].add(a);
1793  }
1794 }
1795 
1796 
1797 
1798 template <class VectorType>
1799 void
1801  const BlockVectorBase<VectorType> &v)
1802 {
1803  AssertIsFinite(a);
1804 
1805  Assert(n_blocks() == v.n_blocks(),
1807 
1808  for (size_type i = 0; i < n_blocks(); ++i)
1809  {
1810  components[i].add(a, v.components[i]);
1811  }
1812 }
1813 
1814 
1815 
1816 template <class VectorType>
1817 void
1819  const BlockVectorBase<VectorType> &v,
1820  const value_type b,
1822 {
1823  AssertIsFinite(a);
1824  AssertIsFinite(b);
1825 
1826  Assert(n_blocks() == v.n_blocks(),
1828  Assert(n_blocks() == w.n_blocks(),
1830 
1831 
1832  for (size_type i = 0; i < n_blocks(); ++i)
1833  {
1834  components[i].add(a, v.components[i], b, w.components[i]);
1835  }
1836 }
1837 
1838 
1839 
1840 template <class VectorType>
1841 void
1843  const BlockVectorBase<VectorType> &v)
1844 {
1845  AssertIsFinite(x);
1846 
1847  Assert(n_blocks() == v.n_blocks(),
1849 
1850  for (size_type i = 0; i < n_blocks(); ++i)
1851  {
1852  components[i].sadd(x, v.components[i]);
1853  }
1854 }
1855 
1856 
1857 
1858 template <class VectorType>
1859 void
1861  const value_type a,
1862  const BlockVectorBase<VectorType> &v)
1863 {
1864  AssertIsFinite(x);
1865  AssertIsFinite(a);
1866 
1867  Assert(n_blocks() == v.n_blocks(),
1869 
1870  for (size_type i = 0; i < n_blocks(); ++i)
1871  {
1872  components[i].sadd(x, a, v.components[i]);
1873  }
1874 }
1875 
1876 
1877 
1878 template <class VectorType>
1879 void
1881  const value_type a,
1882  const BlockVectorBase<VectorType> &v,
1883  const value_type b,
1884  const BlockVectorBase<VectorType> &w)
1885 {
1886  AssertIsFinite(x);
1887  AssertIsFinite(a);
1888  AssertIsFinite(b);
1889 
1890  Assert(n_blocks() == v.n_blocks(),
1892  Assert(n_blocks() == w.n_blocks(),
1894 
1895  for (size_type i = 0; i < n_blocks(); ++i)
1896  {
1897  components[i].sadd(x, a, v.components[i], b, w.components[i]);
1898  }
1899 }
1900 
1901 
1902 
1903 template <class VectorType>
1904 void
1906  const value_type a,
1907  const BlockVectorBase<VectorType> &v,
1908  const value_type b,
1909  const BlockVectorBase<VectorType> &w,
1910  const value_type c,
1911  const BlockVectorBase<VectorType> &y)
1912 {
1913  AssertIsFinite(x);
1914  AssertIsFinite(a);
1915  AssertIsFinite(b);
1916  AssertIsFinite(c);
1917 
1918  Assert(n_blocks() == v.n_blocks(),
1920  Assert(n_blocks() == w.n_blocks(),
1922  Assert(n_blocks() == y.n_blocks(),
1924 
1925  for (size_type i = 0; i < n_blocks(); ++i)
1926  {
1927  components[i].sadd(
1928  x, a, v.components[i], b, w.components[i], c, y.components[i]);
1929  }
1930 }
1931 
1932 
1933 
1934 template <class VectorType>
1935 template <class BlockVector2>
1936 void
1937 BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
1938 {
1939  Assert(n_blocks() == v.n_blocks(),
1940  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1941  for (size_type i = 0; i < n_blocks(); ++i)
1942  components[i].scale(v.block(i));
1943 }
1944 
1945 
1946 
1947 template <class VectorType>
1948 std::size_t
1950 {
1951  return (MemoryConsumption::memory_consumption(this->block_indices) +
1952  MemoryConsumption::memory_consumption(this->components));
1953 }
1954 
1955 
1956 
1957 template <class VectorType>
1958 template <class BlockVector2>
1959 void
1960 BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
1961 {
1962  AssertIsFinite(a);
1963 
1964  Assert(n_blocks() == v.n_blocks(),
1965  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1966 
1967  for (size_type i = 0; i < n_blocks(); ++i)
1968  components[i].equ(a, v.components[i]);
1969 }
1970 
1971 
1972 
1973 template <class VectorType>
1974 void
1976 {
1977  for (size_type i = 0; i < n_blocks(); ++i)
1978  block(i).update_ghost_values();
1979 }
1980 
1981 
1982 
1983 template <class VectorType>
1986 {
1987  AssertIsFinite(s);
1988 
1989  for (size_type i = 0; i < n_blocks(); ++i)
1990  components[i] = s;
1991 
1992  return *this;
1993 }
1994 
1995 
1996 template <class VectorType>
1999 {
2001 
2002  for (size_type i = 0; i < n_blocks(); ++i)
2003  components[i] = v.components[i];
2004 
2005  return *this;
2006 }
2007 
2008 
2009 template <class VectorType>
2010 template <class VectorType2>
2013 {
2015 
2016  for (size_type i = 0; i < n_blocks(); ++i)
2017  components[i] = v.components[i];
2018 
2019  return *this;
2020 }
2021 
2022 
2023 
2024 template <class VectorType>
2027 {
2028  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2029 
2030  size_type index_v = 0;
2031  for (size_type b = 0; b < n_blocks(); ++b)
2032  for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2033  block(b)(i) = v(index_v);
2034 
2035  return *this;
2036 }
2037 
2038 
2039 
2040 template <class VectorType>
2041 template <class VectorType2>
2042 inline bool
2045 {
2046  Assert(block_indices == v.block_indices, ExcDifferentBlockIndices());
2047 
2048  for (size_type i = 0; i < n_blocks(); ++i)
2049  if (!(components[i] == v.components[i]))
2050  return false;
2051 
2052  return true;
2053 }
2054 
2055 
2056 
2057 template <class VectorType>
2060 {
2061  AssertIsFinite(factor);
2062 
2063  for (size_type i = 0; i < n_blocks(); ++i)
2064  components[i] *= factor;
2065 
2066  return *this;
2067 }
2068 
2069 
2070 
2071 template <class VectorType>
2074 {
2075  AssertIsFinite(factor);
2076  Assert(factor != 0., ExcDivideByZero());
2077 
2078  const value_type inverse_factor = value_type(1.) / factor;
2079 
2080  for (size_type i = 0; i < n_blocks(); ++i)
2081  components[i] *= inverse_factor;
2082 
2083  return *this;
2084 }
2085 
2086 
2087 template <class VectorType>
2090 {
2091  const std::pair<unsigned int, size_type> local_index =
2092  block_indices.global_to_local(i);
2093  return components[local_index.first](local_index.second);
2094 }
2095 
2096 
2097 
2098 template <class VectorType>
2101 {
2102  const std::pair<unsigned int, size_type> local_index =
2103  block_indices.global_to_local(i);
2104  return components[local_index.first](local_index.second);
2105 }
2106 
2107 
2108 
2109 template <class VectorType>
2112 {
2113  return operator()(i);
2114 }
2115 
2116 
2117 
2118 template <class VectorType>
2121 {
2122  return operator()(i);
2123 }
2124 
2125 
2126 
2127 template <typename VectorType>
2128 template <typename OtherNumber>
2129 inline void
2131  const std::vector<size_type> &indices,
2132  std::vector<OtherNumber> & values) const
2133 {
2134  for (size_type i = 0; i < indices.size(); ++i)
2135  values[i] = operator()(indices[i]);
2136 }
2137 
2138 
2139 
2140 template <typename VectorType>
2141 template <typename ForwardIterator, typename OutputIterator>
2142 inline void
2144  ForwardIterator indices_begin,
2145  const ForwardIterator indices_end,
2146  OutputIterator values_begin) const
2147 {
2148  while (indices_begin != indices_end)
2149  {
2150  *values_begin = operator()(*indices_begin);
2151  indices_begin++;
2152  values_begin++;
2153  }
2154 }
2155 
2156 #endif // DOXYGEN
2157 
2159 
2160 #endif
const types::global_dof_index invalid_size_type
Definition: types.h:205
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
value_type operator[](const size_type i) const
void add(const std::vector< size_type > &indices, const std::vector< Number > &values)
real_type l2_norm() const
value_type operator*(const BlockVectorBase &V) const
bool operator==(const BlockVectorBase< VectorType2 > &v) const
void compress(::VectorOperation::values operation)
Expression operator>=(const Expression &lhs, const Expression &rhs)
typename numbers::NumberTraits< number >::real_type real_type
Definition: vector.h:140
value_type & reference
Definition: vector.h:127
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
std::vector< value_type > l2_norm(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2043
bool in_local_range(const size_type global_index) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1703
typename std::conditional< Constness, value_type, typename BlockVectorType::BlockType::reference >::type dereference_type
__global__ void add_and_dot(Number *res, Number *v1, const Number *v2, const Number *v3, const Number a, const size_type N)
static const char V
void collect_sizes()
STL namespace.
static const char U
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
std::enable_if< IsBlockVector< VectorType >::value, void >::type collect_sizes(VectorType &vector)
Definition: operators.h:97
size_type block_size(const unsigned int i) const
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
static ::ExceptionBase & ExcDivideByZero()
typename std::conditional< Constness, const typename BlockVectorType::value_type, typename BlockVectorType::value_type >::type value_type
BlockVectorBase & operator/=(const value_type factor)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
value_type add_and_dot(const value_type a, const BlockVectorBase &V, const BlockVectorBase &W)
size_type total_size() const
bool is_non_negative() const
void scale(const BlockVector2 &v)
static std::true_type check_for_block_vector(const BlockVectorBase< T > *)
std::random_access_iterator_tag iterator_category
BlockIndices block_indices
IndexSet locally_owned_elements() const
real_type norm_sqr() const
std::size_t size() const
std::string compress(const std::string &input)
Definition: utilities.cc:392
number value_type
Definition: vector.h:122
std::size_t locally_owned_size() const
T sum(const T &t, const MPI_Comm &mpi_communicator)
Expression operator<=(const Expression &lhs, const Expression &rhs)
const BlockIndices & get_block_indices() const
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void reinit(const unsigned int n_blocks, const size_type n_elements_per_block)
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2817
real_type l1_norm() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
VectorType::value_type * end(VectorType &V)
void update_ghost_values() const
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
size_type local_to_global(const unsigned int block, const size_type index) const
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
Expression operator>(const Expression &lhs, const Expression &rhs)
static ::ExceptionBase & ExcDifferentBlockIndices()
static const bool value
__global__ void equ(Number *val, const Number a, const Number *V_val, const size_type N)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
BlockVectorBase & operator+=(const BlockVectorBase &V)
std::vector< VectorType > components
__global__ void sadd(const Number s, Number *val, const Number a, const Number *V_val, const size_type N)
BlockVectorBase & operator*=(const value_type factor)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
bool all_zero() const
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int global_dof_index
Definition: types.h:76
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< OtherNumber > &values) const
void sadd(const value_type s, const BlockVectorBase &V)
std::size_t memory_consumption() const
void compress() const
Definition: index_set.h:1642
unsigned int n_blocks() const
BarycentricPolynomial< dim, Number1 > operator+(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
BlockVectorBase & operator-=(const BlockVectorBase &V)
typename BlockType::value_type value_type
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
VectorType::value_type * begin(VectorType &V)
typename BlockType::real_type real_type
iterator end()
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2791
static const char N
typename BlockType::const_reference const_reference
typename BlockVectorType::reference reference
size_type block_start(const unsigned int i) const
value_type operator()(const size_type i) const
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
real_type linfty_norm() const
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)
unsigned int size() const
BlockType & block(const unsigned int i)
iterator begin()
BlockVectorBase & operator=(const value_type s)
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
#define AssertIsFinite(number)
Definition: exceptions.h:1721
void equ(const value_type a, const BlockVector2 &V)
typename BlockType::reference reference
const value_type & const_reference
Definition: vector.h:128
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
value_type mean_value() const