Reference documentation for deal.II version Git ede8f93e86 2020-12-03 14:59:20 -0700
\(\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 - 2020 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 
563  IndexSet
564  locally_owned_elements() const;
565 
569  iterator
570  begin();
571 
577  begin() const;
578 
582  iterator
583  end();
584 
590  end() const;
591 
595  value_type
596  operator()(const size_type i) const;
597 
601  reference
602  operator()(const size_type i);
603 
609  value_type operator[](const size_type i) const;
610 
616  reference operator[](const size_type i);
617 
633  template <typename OtherNumber>
634  void
635  extract_subvector_to(const std::vector<size_type> &indices,
636  std::vector<OtherNumber> & values) const;
637 
665  template <typename ForwardIterator, typename OutputIterator>
666  void
667  extract_subvector_to(ForwardIterator indices_begin,
668  const ForwardIterator indices_end,
669  OutputIterator values_begin) const;
670 
676  operator=(const value_type s);
677 
682  operator=(const BlockVectorBase &V);
683 
690  operator=(BlockVectorBase && /*V*/) = default; // NOLINT
691 
695  template <class VectorType2>
697  operator=(const BlockVectorBase<VectorType2> &V);
698 
703  operator=(const VectorType &v);
704 
709  template <class VectorType2>
710  bool
711  operator==(const BlockVectorBase<VectorType2> &v) const;
712 
716  value_type operator*(const BlockVectorBase &V) const;
717 
721  real_type
722  norm_sqr() const;
723 
727  value_type
728  mean_value() const;
729 
733  real_type
734  l1_norm() const;
735 
740  real_type
741  l2_norm() const;
742 
747  real_type
748  linfty_norm() const;
749 
772  value_type
773  add_and_dot(const value_type a,
774  const BlockVectorBase &V,
775  const BlockVectorBase &W);
776 
781  bool
782  in_local_range(const size_type global_index) const;
783 
789  bool
790  all_zero() const;
791 
797  bool
798  is_non_negative() const;
799 
804  operator+=(const BlockVectorBase &V);
805 
810  operator-=(const BlockVectorBase &V);
811 
812 
817  template <typename Number>
818  void
819  add(const std::vector<size_type> &indices, const std::vector<Number> &values);
820 
825  template <typename Number>
826  void
827  add(const std::vector<size_type> &indices, const Vector<Number> &values);
828 
834  template <typename Number>
835  void
836  add(const size_type n_elements,
837  const size_type *indices,
838  const Number * values);
839 
844  void
845  add(const value_type s);
846 
850  void
851  add(const value_type a, const BlockVectorBase &V);
852 
856  void
857  add(const value_type a,
858  const BlockVectorBase &V,
859  const value_type b,
860  const BlockVectorBase &W);
861 
865  void
866  sadd(const value_type s, const BlockVectorBase &V);
867 
871  void
872  sadd(const value_type s, const value_type a, const BlockVectorBase &V);
873 
877  void
878  sadd(const value_type s,
879  const value_type a,
880  const BlockVectorBase &V,
881  const value_type b,
882  const BlockVectorBase &W);
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  const value_type c,
894  const BlockVectorBase &X);
895 
900  operator*=(const value_type factor);
901 
906  operator/=(const value_type factor);
907 
912  template <class BlockVector2>
913  void
914  scale(const BlockVector2 &v);
915 
919  template <class BlockVector2>
920  void
921  equ(const value_type a, const BlockVector2 &V);
922 
927  void
928  update_ghost_values() const;
929 
934  std::size_t
935  memory_consumption() const;
936 
937 protected:
941  std::vector<VectorType> components;
942 
947  BlockIndices block_indices;
948 
949  // Make the iterator class a friend.
950  template <typename N, bool C>
951  friend class ::internal::BlockVectorIterators::Iterator;
952 
953  template <typename>
954  friend class BlockVectorBase;
955 };
956 
957 
960 /*----------------------- Inline functions ----------------------------------*/
961 
962 
963 #ifndef DOXYGEN
964 namespace internal
965 {
966  namespace BlockVectorIterators
967  {
968  template <class BlockVectorType, bool Constness>
969  inline Iterator<BlockVectorType, Constness>::Iterator(
970  const Iterator<BlockVectorType, Constness> &c)
971  : parent(c.parent)
972  , global_index(c.global_index)
973  , current_block(c.current_block)
974  , index_within_block(c.index_within_block)
975  , next_break_forward(c.next_break_forward)
976  , next_break_backward(c.next_break_backward)
977  {}
978 
979 
980 
981  template <class BlockVectorType, bool Constness>
982  inline Iterator<BlockVectorType, Constness>::Iterator(
983  const Iterator<BlockVectorType, !Constness> &c)
984  : parent(c.parent)
985  , global_index(c.global_index)
986  , current_block(c.current_block)
987  , index_within_block(c.index_within_block)
988  , next_break_forward(c.next_break_forward)
989  , next_break_backward(c.next_break_backward)
990  {
991  // Only permit copy-constructing const iterators from non-const
992  // iterators, and not vice versa (i.e., Constness must always be
993  // true).
994  static_assert(Constness == true,
995  "Constructing a non-const iterator from a const iterator "
996  "does not make sense.");
997  }
998 
999 
1000 
1001  template <class BlockVectorType, bool Constness>
1002  inline Iterator<BlockVectorType, Constness>::Iterator(
1003  BlockVector & parent,
1004  const size_type global_index,
1005  const size_type current_block,
1006  const size_type index_within_block,
1007  const size_type next_break_forward,
1008  const size_type next_break_backward)
1009  : parent(&parent)
1010  , global_index(global_index)
1011  , current_block(current_block)
1012  , index_within_block(index_within_block)
1013  , next_break_forward(next_break_forward)
1014  , next_break_backward(next_break_backward)
1015  {}
1016 
1017 
1018 
1019  template <class BlockVectorType, bool Constness>
1020  inline Iterator<BlockVectorType, Constness> &
1021  Iterator<BlockVectorType, Constness>::operator=(const Iterator &c)
1022  {
1023  parent = c.parent;
1024  global_index = c.global_index;
1025  index_within_block = c.index_within_block;
1026  current_block = c.current_block;
1027  next_break_forward = c.next_break_forward;
1028  next_break_backward = c.next_break_backward;
1029 
1030  return *this;
1031  }
1032 
1033 
1034 
1035  template <class BlockVectorType, bool Constness>
1036  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1038  {
1039  return parent->block(current_block)(index_within_block);
1040  }
1041 
1042 
1043 
1044  template <class BlockVectorType, bool Constness>
1045  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1046  Iterator<BlockVectorType, Constness>::
1047  operator[](const difference_type d) const
1048  {
1049  // if the index pointed to is
1050  // still within the block we
1051  // currently point into, then we
1052  // can save the computation of
1053  // the block
1054  if ((global_index + d >= next_break_backward) &&
1055  (global_index + d <= next_break_forward))
1056  return parent->block(current_block)(index_within_block + d);
1057 
1058  // if the index is not within the
1059  // block of the block vector into
1060  // which we presently point, then
1061  // there is no way: we have to
1062  // search for the block. this can
1063  // be done through the parent
1064  // class as well.
1065  return (*parent)(global_index + d);
1066  }
1067 
1068 
1069 
1070  template <class BlockVectorType, bool Constness>
1071  inline Iterator<BlockVectorType, Constness> &
1073  {
1074  move_forward();
1075  return *this;
1076  }
1077 
1078 
1079 
1080  template <class BlockVectorType, bool Constness>
1081  inline Iterator<BlockVectorType, Constness>
1083  {
1084  const Iterator old_value = *this;
1085  move_forward();
1086  return old_value;
1087  }
1088 
1089 
1090 
1091  template <class BlockVectorType, bool Constness>
1092  inline Iterator<BlockVectorType, Constness> &
1093  Iterator<BlockVectorType, Constness>::operator--()
1094  {
1095  move_backward();
1096  return *this;
1097  }
1098 
1099 
1100 
1101  template <class BlockVectorType, bool Constness>
1102  inline Iterator<BlockVectorType, Constness>
1103  Iterator<BlockVectorType, Constness>::operator--(int)
1104  {
1105  const Iterator old_value = *this;
1106  move_backward();
1107  return old_value;
1108  }
1109 
1110 
1111 
1112  template <class BlockVectorType, bool Constness>
1113  template <bool OtherConstness>
1114  inline bool
1116  operator==(const Iterator<BlockVectorType, OtherConstness> &i) const
1117  {
1118  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1119 
1120  return (global_index == i.global_index);
1121  }
1122 
1123 
1124 
1125  template <class BlockVectorType, bool Constness>
1126  template <bool OtherConstness>
1127  inline bool
1129  operator!=(const Iterator<BlockVectorType, OtherConstness> &i) const
1130  {
1131  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1132 
1133  return (global_index != i.global_index);
1134  }
1135 
1136 
1137 
1138  template <class BlockVectorType, bool Constness>
1139  template <bool OtherConstness>
1140  inline bool
1142  operator<(const Iterator<BlockVectorType, OtherConstness> &i) const
1143  {
1144  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1145 
1146  return (global_index < i.global_index);
1147  }
1148 
1149 
1150 
1151  template <class BlockVectorType, bool Constness>
1152  template <bool OtherConstness>
1153  inline bool
1155  operator<=(const Iterator<BlockVectorType, OtherConstness> &i) const
1156  {
1157  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1158 
1159  return (global_index <= i.global_index);
1160  }
1161 
1162 
1163 
1164  template <class BlockVectorType, bool Constness>
1165  template <bool OtherConstness>
1166  inline bool
1168  operator>(const Iterator<BlockVectorType, OtherConstness> &i) const
1169  {
1170  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1171 
1172  return (global_index > i.global_index);
1173  }
1174 
1175 
1176 
1177  template <class BlockVectorType, bool Constness>
1178  template <bool OtherConstness>
1179  inline bool
1181  operator>=(const Iterator<BlockVectorType, OtherConstness> &i) const
1182  {
1183  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1184 
1185  return (global_index >= i.global_index);
1186  }
1187 
1188 
1189 
1190  template <class BlockVectorType, bool Constness>
1191  template <bool OtherConstness>
1192  inline typename Iterator<BlockVectorType, Constness>::difference_type
1194  operator-(const Iterator<BlockVectorType, OtherConstness> &i) const
1195  {
1196  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1197 
1198  return (static_cast<signed int>(global_index) -
1199  static_cast<signed int>(i.global_index));
1200  }
1201 
1202 
1203 
1204  template <class BlockVectorType, bool Constness>
1205  inline Iterator<BlockVectorType, Constness>
1207  operator+(const difference_type &d) const
1208  {
1209  // if the index pointed to is
1210  // still within the block we
1211  // currently point into, then we
1212  // can save the computation of
1213  // the block
1214  if ((global_index + d >= next_break_backward) &&
1215  (global_index + d <= next_break_forward))
1216  return Iterator(*parent,
1217  global_index + d,
1218  current_block,
1219  index_within_block + d,
1220  next_break_forward,
1221  next_break_backward);
1222  else
1223  // outside present block, so
1224  // have to seek new block
1225  // anyway
1226  return Iterator(*parent, global_index + d);
1227  }
1228 
1229 
1230 
1231  template <class BlockVectorType, bool Constness>
1232  inline Iterator<BlockVectorType, Constness>
1234  operator-(const difference_type &d) const
1235  {
1236  // if the index pointed to is
1237  // still within the block we
1238  // currently point into, then we
1239  // can save the computation of
1240  // the block
1241  if ((global_index - d >= next_break_backward) &&
1242  (global_index - d <= next_break_forward))
1243  return Iterator(*parent,
1244  global_index - d,
1245  current_block,
1246  index_within_block - d,
1247  next_break_forward,
1248  next_break_backward);
1249  else
1250  // outside present block, so
1251  // have to seek new block
1252  // anyway
1253  return Iterator(*parent, global_index - d);
1254  }
1255 
1256 
1257 
1258  template <class BlockVectorType, bool Constness>
1259  inline Iterator<BlockVectorType, Constness> &
1260  Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1261  {
1262  // if the index pointed to is
1263  // still within the block we
1264  // currently point into, then we
1265  // can save the computation of
1266  // the block
1267  if ((global_index + d >= next_break_backward) &&
1268  (global_index + d <= next_break_forward))
1269  {
1270  global_index += d;
1271  index_within_block += d;
1272  }
1273  else
1274  // outside present block, so
1275  // have to seek new block
1276  // anyway
1277  *this = Iterator(*parent, global_index + d);
1278 
1279  return *this;
1280  }
1281 
1282 
1283 
1284  template <class BlockVectorType, bool Constness>
1285  inline Iterator<BlockVectorType, Constness> &
1286  Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1287  {
1288  // if the index pointed to is
1289  // still within the block we
1290  // currently point into, then we
1291  // can save the computation of
1292  // the block
1293  if ((global_index - d >= next_break_backward) &&
1294  (global_index - d <= next_break_forward))
1295  {
1296  global_index -= d;
1297  index_within_block -= d;
1298  }
1299  else
1300  // outside present block, so
1301  // have to seek new block
1302  // anyway
1303  *this = Iterator(*parent, global_index - d);
1304 
1305  return *this;
1306  }
1307 
1308 
1309  template <class BlockVectorType, bool Constness>
1310  Iterator<BlockVectorType, Constness>::Iterator(BlockVector & parent,
1311  const size_type global_index)
1312  : parent(&parent)
1313  , global_index(global_index)
1314  {
1315  // find which block we are
1316  // in. for this, take into
1317  // account that it happens at
1318  // times that people want to
1319  // initialize iterators
1320  // past-the-end
1321  if (global_index < parent.size())
1322  {
1323  const std::pair<size_type, size_type> indices =
1324  parent.block_indices.global_to_local(global_index);
1325  current_block = indices.first;
1326  index_within_block = indices.second;
1327 
1328  next_break_backward =
1329  parent.block_indices.local_to_global(current_block, 0);
1330  next_break_forward =
1331  (parent.block_indices.local_to_global(current_block, 0) +
1332  parent.block_indices.block_size(current_block) - 1);
1333  }
1334  else
1335  // past the end. only have one
1336  // value for this
1337  {
1338  this->global_index = parent.size();
1339  current_block = parent.n_blocks();
1340  index_within_block = 0;
1341  next_break_backward = global_index;
1342  next_break_forward = numbers::invalid_size_type;
1343  };
1344  }
1345 
1346 
1347 
1348  template <class BlockVectorType, bool Constness>
1349  void
1350  Iterator<BlockVectorType, Constness>::move_forward()
1351  {
1352  if (global_index != next_break_forward)
1353  ++index_within_block;
1354  else
1355  {
1356  // ok, we traverse a boundary
1357  // between blocks:
1358  index_within_block = 0;
1359  ++current_block;
1360 
1361  // break backwards is now old
1362  // break forward
1363  next_break_backward = next_break_forward + 1;
1364 
1365  // compute new break forward
1366  if (current_block < parent->block_indices.size())
1367  next_break_forward +=
1368  parent->block_indices.block_size(current_block);
1369  else
1370  // if we are beyond the end,
1371  // then move the next
1372  // boundary arbitrarily far
1373  // away
1374  next_break_forward = numbers::invalid_size_type;
1375  };
1376 
1377  ++global_index;
1378  }
1379 
1380 
1381 
1382  template <class BlockVectorType, bool Constness>
1383  void
1384  Iterator<BlockVectorType, Constness>::move_backward()
1385  {
1386  if (global_index != next_break_backward)
1387  --index_within_block;
1388  else if (current_block != 0)
1389  {
1390  // ok, we traverse a boundary
1391  // between blocks:
1392  --current_block;
1393  index_within_block =
1394  parent->block_indices.block_size(current_block) - 1;
1395 
1396  // break forwards is now old
1397  // break backward
1398  next_break_forward = next_break_backward - 1;
1399 
1400  // compute new break forward
1401  next_break_backward -=
1402  parent->block_indices.block_size(current_block);
1403  }
1404  else
1405  // current block was 0, we now
1406  // get into unspecified terrain
1407  {
1408  --current_block;
1409  index_within_block = numbers::invalid_size_type;
1410  next_break_forward = 0;
1411  next_break_backward = 0;
1412  };
1413 
1414  --global_index;
1415  }
1416 
1417 
1418  } // namespace BlockVectorIterators
1419 
1420 } // namespace internal
1421 
1422 
1423 
1424 template <class VectorType>
1425 inline std::size_t
1427 {
1428  return block_indices.total_size();
1429 }
1430 
1431 
1432 
1433 template <class VectorType>
1434 inline IndexSet
1436 {
1437  IndexSet is(size());
1438 
1439  // copy index sets from blocks into the global one, shifted
1440  // by the appropriate amount for each block
1441  for (unsigned int b = 0; b < n_blocks(); ++b)
1442  {
1443  IndexSet x = block(b).locally_owned_elements();
1444  is.add_indices(x, block_indices.block_start(b));
1445  }
1446 
1447  is.compress();
1448 
1449  return is;
1450 }
1451 
1452 
1453 
1454 template <class VectorType>
1455 inline unsigned int
1457 {
1458  return block_indices.size();
1459 }
1460 
1461 
1462 template <class VectorType>
1464 BlockVectorBase<VectorType>::block(const unsigned int i)
1465 {
1466  AssertIndexRange(i, n_blocks());
1467 
1468  return components[i];
1469 }
1470 
1471 
1472 
1473 template <class VectorType>
1474 inline const typename BlockVectorBase<VectorType>::BlockType &
1475 BlockVectorBase<VectorType>::block(const unsigned int i) const
1476 {
1477  AssertIndexRange(i, n_blocks());
1478 
1479  return components[i];
1480 }
1481 
1482 
1483 
1484 template <class VectorType>
1485 inline const BlockIndices &
1487 {
1488  return block_indices;
1489 }
1490 
1491 
1492 template <class VectorType>
1493 inline void
1495 {
1496  std::vector<size_type> sizes(n_blocks());
1497 
1498  for (size_type i = 0; i < n_blocks(); ++i)
1499  sizes[i] = block(i).size();
1500 
1501  block_indices.reinit(sizes);
1502 }
1503 
1504 
1505 
1506 template <class VectorType>
1507 inline void
1509  ::VectorOperation::values operation)
1510 {
1511  for (unsigned int i = 0; i < n_blocks(); ++i)
1512  block(i).compress(operation);
1513 }
1514 
1515 
1516 
1517 template <class VectorType>
1520 {
1521  return iterator(*this, 0U);
1522 }
1523 
1524 
1525 
1526 template <class VectorType>
1529 {
1530  return const_iterator(*this, 0U);
1531 }
1532 
1533 
1534 template <class VectorType>
1537 {
1538  return iterator(*this, size());
1539 }
1540 
1541 
1542 
1543 template <class VectorType>
1546 {
1547  return const_iterator(*this, size());
1548 }
1549 
1550 
1551 template <class VectorType>
1552 inline bool
1554 {
1555  const std::pair<size_type, size_type> local_index =
1556  block_indices.global_to_local(global_index);
1557 
1558  return components[local_index.first].in_local_range(global_index);
1559 }
1560 
1561 
1562 template <class VectorType>
1563 bool
1565 {
1566  for (size_type i = 0; i < n_blocks(); ++i)
1567  if (components[i].all_zero() == false)
1568  return false;
1569 
1570  return true;
1571 }
1572 
1573 
1574 
1575 template <class VectorType>
1576 bool
1578 {
1579  for (size_type i = 0; i < n_blocks(); ++i)
1580  if (components[i].is_non_negative() == false)
1581  return false;
1582 
1583  return true;
1584 }
1585 
1586 
1587 
1588 template <class VectorType>
1591 {
1592  Assert(n_blocks() == v.n_blocks(),
1594 
1595  value_type sum = 0.;
1596  for (size_type i = 0; i < n_blocks(); ++i)
1597  sum += components[i] * v.components[i];
1598 
1599  return sum;
1600 }
1601 
1602 
1603 template <class VectorType>
1606 {
1607  real_type sum = 0.;
1608  for (size_type i = 0; i < n_blocks(); ++i)
1609  sum += components[i].norm_sqr();
1610 
1611  return sum;
1612 }
1613 
1614 
1615 
1616 template <class VectorType>
1619 {
1620  value_type sum = 0.;
1621  // need to do static_cast as otherwise it won't work with
1622  // value_type=complex<T>
1623  for (size_type i = 0; i < n_blocks(); ++i)
1624  sum += components[i].mean_value() *
1626  components[i].size()));
1627 
1628  return sum / (typename numbers::NumberTraits<value_type>::real_type(size()));
1629 }
1630 
1631 
1632 
1633 template <class VectorType>
1636 {
1637  real_type sum = 0.;
1638  for (size_type i = 0; i < n_blocks(); ++i)
1639  sum += components[i].l1_norm();
1640 
1641  return sum;
1642 }
1643 
1644 
1645 
1646 template <class VectorType>
1649 {
1650  return std::sqrt(norm_sqr());
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  {
1662  value_type newval = components[i].linfty_norm();
1663  if (sum < newval)
1664  sum = newval;
1665  }
1666  return sum;
1667 }
1668 
1669 
1670 
1671 template <class VectorType>
1674  const typename BlockVectorBase<VectorType>::value_type a,
1675  const BlockVectorBase<VectorType> & V,
1676  const BlockVectorBase<VectorType> & W)
1677 {
1680 
1681  value_type sum = 0.;
1682  for (size_type i = 0; i < n_blocks(); ++i)
1683  sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1684 
1685  return sum;
1686 }
1687 
1688 
1689 
1690 template <class VectorType>
1693 {
1694  Assert(n_blocks() == v.n_blocks(),
1696 
1697  for (size_type i = 0; i < n_blocks(); ++i)
1698  {
1699  components[i] += v.components[i];
1700  }
1701 
1702  return *this;
1703 }
1704 
1705 
1706 
1707 template <class VectorType>
1710 {
1711  Assert(n_blocks() == v.n_blocks(),
1713 
1714  for (size_type i = 0; i < n_blocks(); ++i)
1715  {
1716  components[i] -= v.components[i];
1717  }
1718  return *this;
1719 }
1720 
1721 
1722 
1723 template <class VectorType>
1724 template <typename Number>
1725 inline void
1726 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1727  const std::vector<Number> & values)
1728 {
1729  Assert(indices.size() == values.size(),
1730  ExcDimensionMismatch(indices.size(), values.size()));
1731  add(indices.size(), indices.data(), values.data());
1732 }
1733 
1734 
1735 
1736 template <class VectorType>
1737 template <typename Number>
1738 inline void
1739 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1740  const Vector<Number> & values)
1741 {
1742  Assert(indices.size() == values.size(),
1743  ExcDimensionMismatch(indices.size(), values.size()));
1744  const size_type n_indices = indices.size();
1745  for (size_type i = 0; i < n_indices; ++i)
1746  (*this)(indices[i]) += values(i);
1747 }
1748 
1749 
1750 
1751 template <class VectorType>
1752 template <typename Number>
1753 inline void
1755  const size_type *indices,
1756  const Number * values)
1757 {
1758  for (size_type i = 0; i < n_indices; ++i)
1759  (*this)(indices[i]) += values[i];
1760 }
1761 
1762 
1763 
1764 template <class VectorType>
1765 void
1767 {
1768  AssertIsFinite(a);
1769 
1770  for (size_type i = 0; i < n_blocks(); ++i)
1771  {
1772  components[i].add(a);
1773  }
1774 }
1775 
1776 
1777 
1778 template <class VectorType>
1779 void
1781  const BlockVectorBase<VectorType> &v)
1782 {
1783  AssertIsFinite(a);
1784 
1785  Assert(n_blocks() == v.n_blocks(),
1787 
1788  for (size_type i = 0; i < n_blocks(); ++i)
1789  {
1790  components[i].add(a, v.components[i]);
1791  }
1792 }
1793 
1794 
1795 
1796 template <class VectorType>
1797 void
1799  const BlockVectorBase<VectorType> &v,
1800  const value_type b,
1802 {
1803  AssertIsFinite(a);
1804  AssertIsFinite(b);
1805 
1806  Assert(n_blocks() == v.n_blocks(),
1808  Assert(n_blocks() == w.n_blocks(),
1810 
1811 
1812  for (size_type i = 0; i < n_blocks(); ++i)
1813  {
1814  components[i].add(a, v.components[i], b, w.components[i]);
1815  }
1816 }
1817 
1818 
1819 
1820 template <class VectorType>
1821 void
1823  const BlockVectorBase<VectorType> &v)
1824 {
1825  AssertIsFinite(x);
1826 
1827  Assert(n_blocks() == v.n_blocks(),
1829 
1830  for (size_type i = 0; i < n_blocks(); ++i)
1831  {
1832  components[i].sadd(x, v.components[i]);
1833  }
1834 }
1835 
1836 
1837 
1838 template <class VectorType>
1839 void
1841  const value_type a,
1842  const BlockVectorBase<VectorType> &v)
1843 {
1844  AssertIsFinite(x);
1845  AssertIsFinite(a);
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, a, 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  const value_type b,
1864  const BlockVectorBase<VectorType> &w)
1865 {
1866  AssertIsFinite(x);
1867  AssertIsFinite(a);
1868  AssertIsFinite(b);
1869 
1870  Assert(n_blocks() == v.n_blocks(),
1872  Assert(n_blocks() == w.n_blocks(),
1874 
1875  for (size_type i = 0; i < n_blocks(); ++i)
1876  {
1877  components[i].sadd(x, a, v.components[i], b, w.components[i]);
1878  }
1879 }
1880 
1881 
1882 
1883 template <class VectorType>
1884 void
1886  const value_type a,
1887  const BlockVectorBase<VectorType> &v,
1888  const value_type b,
1889  const BlockVectorBase<VectorType> &w,
1890  const value_type c,
1891  const BlockVectorBase<VectorType> &y)
1892 {
1893  AssertIsFinite(x);
1894  AssertIsFinite(a);
1895  AssertIsFinite(b);
1896  AssertIsFinite(c);
1897 
1898  Assert(n_blocks() == v.n_blocks(),
1900  Assert(n_blocks() == w.n_blocks(),
1902  Assert(n_blocks() == y.n_blocks(),
1904 
1905  for (size_type i = 0; i < n_blocks(); ++i)
1906  {
1907  components[i].sadd(
1908  x, a, v.components[i], b, w.components[i], c, y.components[i]);
1909  }
1910 }
1911 
1912 
1913 
1914 template <class VectorType>
1915 template <class BlockVector2>
1916 void
1917 BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
1918 {
1919  Assert(n_blocks() == v.n_blocks(),
1920  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1921  for (size_type i = 0; i < n_blocks(); ++i)
1922  components[i].scale(v.block(i));
1923 }
1924 
1925 
1926 
1927 template <class VectorType>
1928 std::size_t
1930 {
1931  return (MemoryConsumption::memory_consumption(this->block_indices) +
1932  MemoryConsumption::memory_consumption(this->components));
1933 }
1934 
1935 
1936 
1937 template <class VectorType>
1938 template <class BlockVector2>
1939 void
1940 BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
1941 {
1942  AssertIsFinite(a);
1943 
1944  Assert(n_blocks() == v.n_blocks(),
1945  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1946 
1947  for (size_type i = 0; i < n_blocks(); ++i)
1948  components[i].equ(a, v.components[i]);
1949 }
1950 
1951 
1952 
1953 template <class VectorType>
1954 void
1956 {
1957  for (size_type i = 0; i < n_blocks(); ++i)
1958  block(i).update_ghost_values();
1959 }
1960 
1961 
1962 
1963 template <class VectorType>
1966 {
1967  AssertIsFinite(s);
1968 
1969  for (size_type i = 0; i < n_blocks(); ++i)
1970  components[i] = s;
1971 
1972  return *this;
1973 }
1974 
1975 
1976 template <class VectorType>
1979 {
1981 
1982  for (size_type i = 0; i < n_blocks(); ++i)
1983  components[i] = v.components[i];
1984 
1985  return *this;
1986 }
1987 
1988 
1989 template <class VectorType>
1990 template <class VectorType2>
1993 {
1995 
1996  for (size_type i = 0; i < n_blocks(); ++i)
1997  components[i] = v.components[i];
1998 
1999  return *this;
2000 }
2001 
2002 
2003 
2004 template <class VectorType>
2007 {
2008  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2009 
2010  size_type index_v = 0;
2011  for (size_type b = 0; b < n_blocks(); ++b)
2012  for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2013  block(b)(i) = v(index_v);
2014 
2015  return *this;
2016 }
2017 
2018 
2019 
2020 template <class VectorType>
2021 template <class VectorType2>
2022 inline bool
2025 {
2026  Assert(block_indices == v.block_indices, ExcDifferentBlockIndices());
2027 
2028  for (size_type i = 0; i < n_blocks(); ++i)
2029  if (!(components[i] == v.components[i]))
2030  return false;
2031 
2032  return true;
2033 }
2034 
2035 
2036 
2037 template <class VectorType>
2040 {
2041  AssertIsFinite(factor);
2042 
2043  for (size_type i = 0; i < n_blocks(); ++i)
2044  components[i] *= factor;
2045 
2046  return *this;
2047 }
2048 
2049 
2050 
2051 template <class VectorType>
2054 {
2055  AssertIsFinite(factor);
2056  Assert(factor != 0., ExcDivideByZero());
2057 
2058  const value_type inverse_factor = value_type(1.) / factor;
2059 
2060  for (size_type i = 0; i < n_blocks(); ++i)
2061  components[i] *= inverse_factor;
2062 
2063  return *this;
2064 }
2065 
2066 
2067 template <class VectorType>
2070 {
2071  const std::pair<unsigned int, size_type> local_index =
2072  block_indices.global_to_local(i);
2073  return components[local_index.first](local_index.second);
2074 }
2075 
2076 
2077 
2078 template <class VectorType>
2081 {
2082  const std::pair<unsigned int, size_type> local_index =
2083  block_indices.global_to_local(i);
2084  return components[local_index.first](local_index.second);
2085 }
2086 
2087 
2088 
2089 template <class VectorType>
2092 {
2093  return operator()(i);
2094 }
2095 
2096 
2097 
2098 template <class VectorType>
2101 {
2102  return operator()(i);
2103 }
2104 
2105 
2106 
2107 template <typename VectorType>
2108 template <typename OtherNumber>
2109 inline void
2111  const std::vector<size_type> &indices,
2112  std::vector<OtherNumber> & values) const
2113 {
2114  for (size_type i = 0; i < indices.size(); ++i)
2115  values[i] = operator()(indices[i]);
2116 }
2117 
2118 
2119 
2120 template <typename VectorType>
2121 template <typename ForwardIterator, typename OutputIterator>
2122 inline void
2124  ForwardIterator indices_begin,
2125  const ForwardIterator indices_end,
2126  OutputIterator values_begin) const
2127 {
2128  while (indices_begin != indices_end)
2129  {
2130  *values_begin = operator()(*indices_begin);
2131  indices_begin++;
2132  values_begin++;
2133  }
2134 }
2135 
2136 #endif // DOXYGEN
2137 
2139 
2140 #endif
const types::global_dof_index invalid_size_type
Definition: types.h:205
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
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:951
bool in_local_range(const size_type global_index) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1701
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
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
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
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:1466
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
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:2816
real_type l1_norm() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
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:1640
unsigned int n_blocks() const
BlockVectorBase & operator-=(const BlockVectorBase &V)
typename BlockType::value_type value_type
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
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:2790
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:1722
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