Reference documentation for deal.II version Git 497f915867 2021-09-17 22:46:48 +0200
\(\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 
149  using iterator_category = std::random_access_iterator_tag;
150  using difference_type = std::ptrdiff_t;
151  using reference = typename BlockVectorType::reference;
152  using pointer = value_type *;
153 
154  using dereference_type = typename std::conditional<
155  Constness,
156  value_type,
157  typename BlockVectorType::BlockType::reference>::type;
158 
163  using BlockVector = typename std::
164  conditional<Constness, const BlockVectorType, BlockVectorType>::type;
165 
174  Iterator(BlockVector &parent, const size_type global_index);
175 
184 
185 
189  Iterator(const Iterator &c);
190 
191  private:
196  Iterator(BlockVector & parent,
197  const size_type global_index,
198  const size_type current_block,
199  const size_type index_within_block,
200  const size_type next_break_forward,
201  const size_type next_break_backward);
202 
203  public:
207  Iterator &
208  operator=(const Iterator &c);
209 
216  operator*() const;
217 
223  operator[](const difference_type d) const;
224 
229  Iterator &
230  operator++();
231 
237  Iterator
238  operator++(int);
239 
244  Iterator &
245  operator--();
246 
252  Iterator
253  operator--(int);
254 
259  template <bool OtherConstness>
260  bool
262 
267  template <bool OtherConstness>
268  bool
270 
276  template <bool OtherConstness>
277  bool
278  operator<(const Iterator<BlockVectorType, OtherConstness> &i) const;
279 
283  template <bool OtherConstness>
284  bool
285  operator<=(const Iterator<BlockVectorType, OtherConstness> &i) const;
286 
290  template <bool OtherConstness>
291  bool
293 
297  template <bool OtherConstness>
298  bool
300 
304  template <bool OtherConstness>
307 
312  Iterator
313  operator+(const difference_type &d) const;
314 
319  Iterator
320  operator-(const difference_type &d) const;
321 
326  Iterator &
327  operator+=(const difference_type &d);
328 
333  Iterator &
334  operator-=(const difference_type &d);
335 
345  DeclExceptionMsg(ExcPointerToDifferentVectors,
346  "Your program tried to compare iterators pointing to "
347  "different block vectors. There is no reasonable way "
348  "to do this.");
349 
351  private:
358 
363 
368  unsigned int current_block;
370 
379 
383  void
384  move_forward();
385 
389  void
390  move_backward();
391 
392  // Mark all other instances of this template as friends.
393  template <typename, bool>
394  friend class Iterator;
395  };
396  } // namespace BlockVectorIterators
397 } // namespace internal
398 
399 
437 template <class VectorType>
439 {
440 public:
445 
446  /*
447  * Declare standard types used in
448  * all containers. These types
449  * parallel those in the
450  * <tt>C++</tt> standard
451  * libraries
452  * <tt>std::vector<...></tt>
453  * class. This includes iterator
454  * types.
455  */
457  using pointer = value_type *;
458  using const_pointer = const value_type *;
459  using iterator =
461  using const_iterator =
463  using reference = typename BlockType::reference;
466 
476  using real_type = typename BlockType::real_type;
477 
481  BlockVectorBase() = default;
482 
486  BlockVectorBase(const BlockVectorBase & /*V*/) = default;
487 
493  BlockVectorBase(BlockVectorBase && /*V*/) noexcept = default;
494 
501  void
502  collect_sizes();
503 
514  void
515  compress(::VectorOperation::values operation);
516 
520  BlockType &
521  block(const unsigned int i);
522 
526  const BlockType &
527  block(const unsigned int i) const;
528 
534  const BlockIndices &
535  get_block_indices() const;
536 
540  unsigned int
541  n_blocks() const;
542 
547  std::size_t
548  size() const;
549 
555  std::size_t
556  locally_owned_size() const;
557 
574  IndexSet
575  locally_owned_elements() const;
576 
580  iterator
581  begin();
582 
588  begin() const;
589 
593  iterator
594  end();
595 
601  end() const;
602 
606  value_type
607  operator()(const size_type i) const;
608 
612  reference
613  operator()(const size_type i);
614 
620  value_type
621  operator[](const size_type i) const;
622 
628  reference
629  operator[](const size_type i);
630 
646  template <typename OtherNumber>
647  void
648  extract_subvector_to(const std::vector<size_type> &indices,
649  std::vector<OtherNumber> & values) const;
650 
678  template <typename ForwardIterator, typename OutputIterator>
679  void
680  extract_subvector_to(ForwardIterator indices_begin,
681  const ForwardIterator indices_end,
682  OutputIterator values_begin) const;
683 
689  operator=(const value_type s);
690 
695  operator=(const BlockVectorBase &V);
696 
703  operator=(BlockVectorBase && /*V*/) = default; // NOLINT
704 
708  template <class VectorType2>
710  operator=(const BlockVectorBase<VectorType2> &V);
711 
716  operator=(const VectorType &v);
717 
722  template <class VectorType2>
723  bool
724  operator==(const BlockVectorBase<VectorType2> &v) const;
725 
729  value_type
730  operator*(const BlockVectorBase &V) const;
731 
735  real_type
736  norm_sqr() const;
737 
741  value_type
742  mean_value() const;
743 
747  real_type
748  l1_norm() const;
749 
754  real_type
755  l2_norm() const;
756 
761  real_type
762  linfty_norm() const;
763 
786  value_type
787  add_and_dot(const value_type a,
788  const BlockVectorBase &V,
789  const BlockVectorBase &W);
790 
795  bool
796  in_local_range(const size_type global_index) const;
797 
803  bool
804  all_zero() const;
805 
811  bool
812  is_non_negative() const;
813 
818  operator+=(const BlockVectorBase &V);
819 
824  operator-=(const BlockVectorBase &V);
825 
826 
831  template <typename Number>
832  void
833  add(const std::vector<size_type> &indices, const std::vector<Number> &values);
834 
839  template <typename Number>
840  void
841  add(const std::vector<size_type> &indices, const Vector<Number> &values);
842 
848  template <typename Number>
849  void
850  add(const size_type n_elements,
851  const size_type *indices,
852  const Number * values);
853 
858  void
859  add(const value_type s);
860 
864  void
865  add(const value_type a, const BlockVectorBase &V);
866 
870  void
871  add(const value_type a,
872  const BlockVectorBase &V,
873  const value_type b,
874  const BlockVectorBase &W);
875 
879  void
880  sadd(const value_type s, const BlockVectorBase &V);
881 
885  void
886  sadd(const value_type s, const value_type a, const BlockVectorBase &V);
887 
891  void
892  sadd(const value_type s,
893  const value_type a,
894  const BlockVectorBase &V,
895  const value_type b,
896  const BlockVectorBase &W);
897 
901  void
902  sadd(const value_type s,
903  const value_type a,
904  const BlockVectorBase &V,
905  const value_type b,
906  const BlockVectorBase &W,
907  const value_type c,
908  const BlockVectorBase &X);
909 
914  operator*=(const value_type factor);
915 
920  operator/=(const value_type factor);
921 
926  template <class BlockVector2>
927  void
928  scale(const BlockVector2 &v);
929 
933  template <class BlockVector2>
934  void
935  equ(const value_type a, const BlockVector2 &V);
936 
941  void
942  update_ghost_values() const;
943 
948  std::size_t
949  memory_consumption() const;
950 
951 protected:
955  std::vector<VectorType> components;
956 
961  BlockIndices block_indices;
962 
963  // Make the iterator class a friend.
964  template <typename N, bool C>
965  friend class ::internal::BlockVectorIterators::Iterator;
966 
967  template <typename>
968  friend class BlockVectorBase;
969 };
970 
971 
974 /*----------------------- Inline functions ----------------------------------*/
975 
976 
977 #ifndef DOXYGEN
978 namespace internal
979 {
980  namespace BlockVectorIterators
981  {
982  template <class BlockVectorType, bool Constness>
983  inline Iterator<BlockVectorType, Constness>::Iterator(
984  const Iterator<BlockVectorType, Constness> &c)
985  : parent(c.parent)
986  , global_index(c.global_index)
987  , current_block(c.current_block)
988  , index_within_block(c.index_within_block)
989  , next_break_forward(c.next_break_forward)
990  , next_break_backward(c.next_break_backward)
991  {}
992 
993 
994 
995  template <class BlockVectorType, bool Constness>
996  inline Iterator<BlockVectorType, Constness>::Iterator(
997  const Iterator<BlockVectorType, !Constness> &c)
998  : parent(c.parent)
999  , global_index(c.global_index)
1000  , current_block(c.current_block)
1001  , index_within_block(c.index_within_block)
1002  , next_break_forward(c.next_break_forward)
1003  , next_break_backward(c.next_break_backward)
1004  {
1005  // Only permit copy-constructing const iterators from non-const
1006  // iterators, and not vice versa (i.e., Constness must always be
1007  // true).
1008  static_assert(Constness == true,
1009  "Constructing a non-const iterator from a const iterator "
1010  "does not make sense.");
1011  }
1012 
1013 
1014 
1015  template <class BlockVectorType, bool Constness>
1016  inline Iterator<BlockVectorType, Constness>::Iterator(
1017  BlockVector & parent,
1018  const size_type global_index,
1019  const size_type current_block,
1020  const size_type index_within_block,
1021  const size_type next_break_forward,
1022  const size_type next_break_backward)
1023  : parent(&parent)
1024  , global_index(global_index)
1025  , current_block(current_block)
1026  , index_within_block(index_within_block)
1027  , next_break_forward(next_break_forward)
1028  , next_break_backward(next_break_backward)
1029  {}
1030 
1031 
1032 
1033  template <class BlockVectorType, bool Constness>
1034  inline Iterator<BlockVectorType, Constness> &
1035  Iterator<BlockVectorType, Constness>::operator=(const Iterator &c)
1036  {
1037  parent = c.parent;
1038  global_index = c.global_index;
1039  index_within_block = c.index_within_block;
1040  current_block = c.current_block;
1041  next_break_forward = c.next_break_forward;
1042  next_break_backward = c.next_break_backward;
1043 
1044  return *this;
1045  }
1046 
1047 
1048 
1049  template <class BlockVectorType, bool Constness>
1050  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1052  {
1053  return parent->block(current_block)(index_within_block);
1054  }
1055 
1056 
1057 
1058  template <class BlockVectorType, bool Constness>
1059  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1060  Iterator<BlockVectorType, Constness>::operator[](
1061  const difference_type d) const
1062  {
1063  // if the index pointed to is
1064  // still within the block we
1065  // currently point into, then we
1066  // can save the computation of
1067  // the block
1068  if ((global_index + d >= next_break_backward) &&
1069  (global_index + d <= next_break_forward))
1070  return parent->block(current_block)(index_within_block + d);
1071 
1072  // if the index is not within the
1073  // block of the block vector into
1074  // which we presently point, then
1075  // there is no way: we have to
1076  // search for the block. this can
1077  // be done through the parent
1078  // class as well.
1079  return (*parent)(global_index + d);
1080  }
1081 
1082 
1083 
1084  template <class BlockVectorType, bool Constness>
1085  inline Iterator<BlockVectorType, Constness> &
1087  {
1088  move_forward();
1089  return *this;
1090  }
1091 
1092 
1093 
1094  template <class BlockVectorType, bool Constness>
1095  inline Iterator<BlockVectorType, Constness>
1097  {
1098  const Iterator old_value = *this;
1099  move_forward();
1100  return old_value;
1101  }
1102 
1103 
1104 
1105  template <class BlockVectorType, bool Constness>
1106  inline Iterator<BlockVectorType, Constness> &
1107  Iterator<BlockVectorType, Constness>::operator--()
1108  {
1109  move_backward();
1110  return *this;
1111  }
1112 
1113 
1114 
1115  template <class BlockVectorType, bool Constness>
1116  inline Iterator<BlockVectorType, Constness>
1117  Iterator<BlockVectorType, Constness>::operator--(int)
1118  {
1119  const Iterator old_value = *this;
1120  move_backward();
1121  return old_value;
1122  }
1123 
1124 
1125 
1126  template <class BlockVectorType, bool Constness>
1127  template <bool OtherConstness>
1128  inline bool
1130  const Iterator<BlockVectorType, OtherConstness> &i) const
1131  {
1132  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1133 
1134  return (global_index == i.global_index);
1135  }
1136 
1137 
1138 
1139  template <class BlockVectorType, bool Constness>
1140  template <bool OtherConstness>
1141  inline bool
1143  const Iterator<BlockVectorType, OtherConstness> &i) const
1144  {
1145  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1146 
1147  return (global_index != i.global_index);
1148  }
1149 
1150 
1151 
1152  template <class BlockVectorType, bool Constness>
1153  template <bool OtherConstness>
1154  inline bool
1156  const Iterator<BlockVectorType, OtherConstness> &i) const
1157  {
1158  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1159 
1160  return (global_index < i.global_index);
1161  }
1162 
1163 
1164 
1165  template <class BlockVectorType, bool Constness>
1166  template <bool OtherConstness>
1167  inline bool
1169  const Iterator<BlockVectorType, OtherConstness> &i) const
1170  {
1171  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1172 
1173  return (global_index <= i.global_index);
1174  }
1175 
1176 
1177 
1178  template <class BlockVectorType, bool Constness>
1179  template <bool OtherConstness>
1180  inline bool
1182  const Iterator<BlockVectorType, OtherConstness> &i) const
1183  {
1184  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1185 
1186  return (global_index > i.global_index);
1187  }
1188 
1189 
1190 
1191  template <class BlockVectorType, bool Constness>
1192  template <bool OtherConstness>
1193  inline bool
1195  const Iterator<BlockVectorType, OtherConstness> &i) const
1196  {
1197  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1198 
1199  return (global_index >= i.global_index);
1200  }
1201 
1202 
1203 
1204  template <class BlockVectorType, bool Constness>
1205  template <bool OtherConstness>
1206  inline typename Iterator<BlockVectorType, Constness>::difference_type
1208  const Iterator<BlockVectorType, OtherConstness> &i) const
1209  {
1210  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1211 
1212  return (static_cast<signed int>(global_index) -
1213  static_cast<signed int>(i.global_index));
1214  }
1215 
1216 
1217 
1218  template <class BlockVectorType, bool Constness>
1219  inline Iterator<BlockVectorType, Constness>
1221  const difference_type &d) const
1222  {
1223  // if the index pointed to is
1224  // still within the block we
1225  // currently point into, then we
1226  // can save the computation of
1227  // the block
1228  if ((global_index + d >= next_break_backward) &&
1229  (global_index + d <= next_break_forward))
1230  return Iterator(*parent,
1231  global_index + d,
1232  current_block,
1233  index_within_block + d,
1234  next_break_forward,
1235  next_break_backward);
1236  else
1237  // outside present block, so
1238  // have to seek new block
1239  // anyway
1240  return Iterator(*parent, global_index + d);
1241  }
1242 
1243 
1244 
1245  template <class BlockVectorType, bool Constness>
1246  inline Iterator<BlockVectorType, Constness>
1248  const difference_type &d) const
1249  {
1250  // if the index pointed to is
1251  // still within the block we
1252  // currently point into, then we
1253  // can save the computation of
1254  // the block
1255  if ((global_index - d >= next_break_backward) &&
1256  (global_index - d <= next_break_forward))
1257  return Iterator(*parent,
1258  global_index - d,
1259  current_block,
1260  index_within_block - d,
1261  next_break_forward,
1262  next_break_backward);
1263  else
1264  // outside present block, so
1265  // have to seek new block
1266  // anyway
1267  return Iterator(*parent, global_index - d);
1268  }
1269 
1270 
1271 
1272  template <class BlockVectorType, bool Constness>
1273  inline Iterator<BlockVectorType, Constness> &
1274  Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1275  {
1276  // if the index pointed to is
1277  // still within the block we
1278  // currently point into, then we
1279  // can save the computation of
1280  // the block
1281  if ((global_index + d >= next_break_backward) &&
1282  (global_index + d <= next_break_forward))
1283  {
1284  global_index += d;
1285  index_within_block += d;
1286  }
1287  else
1288  // outside present block, so
1289  // have to seek new block
1290  // anyway
1291  *this = Iterator(*parent, global_index + d);
1292 
1293  return *this;
1294  }
1295 
1296 
1297 
1298  template <class BlockVectorType, bool Constness>
1299  inline Iterator<BlockVectorType, Constness> &
1300  Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1301  {
1302  // if the index pointed to is
1303  // still within the block we
1304  // currently point into, then we
1305  // can save the computation of
1306  // the block
1307  if ((global_index - d >= next_break_backward) &&
1308  (global_index - d <= next_break_forward))
1309  {
1310  global_index -= d;
1311  index_within_block -= d;
1312  }
1313  else
1314  // outside present block, so
1315  // have to seek new block
1316  // anyway
1317  *this = Iterator(*parent, global_index - d);
1318 
1319  return *this;
1320  }
1321 
1322 
1323  template <class BlockVectorType, bool Constness>
1324  Iterator<BlockVectorType, Constness>::Iterator(BlockVector & parent,
1325  const size_type global_index)
1326  : parent(&parent)
1327  , global_index(global_index)
1328  {
1329  // find which block we are
1330  // in. for this, take into
1331  // account that it happens at
1332  // times that people want to
1333  // initialize iterators
1334  // past-the-end
1335  if (global_index < parent.size())
1336  {
1337  const std::pair<size_type, size_type> indices =
1338  parent.block_indices.global_to_local(global_index);
1339  current_block = indices.first;
1340  index_within_block = indices.second;
1341 
1342  next_break_backward =
1343  parent.block_indices.local_to_global(current_block, 0);
1344  next_break_forward =
1345  (parent.block_indices.local_to_global(current_block, 0) +
1346  parent.block_indices.block_size(current_block) - 1);
1347  }
1348  else
1349  // past the end. only have one
1350  // value for this
1351  {
1352  this->global_index = parent.size();
1353  current_block = parent.n_blocks();
1354  index_within_block = 0;
1355  next_break_backward = global_index;
1356  next_break_forward = numbers::invalid_size_type;
1357  };
1358  }
1359 
1360 
1361 
1362  template <class BlockVectorType, bool Constness>
1363  void
1364  Iterator<BlockVectorType, Constness>::move_forward()
1365  {
1366  if (global_index != next_break_forward)
1367  ++index_within_block;
1368  else
1369  {
1370  // ok, we traverse a boundary
1371  // between blocks:
1372  index_within_block = 0;
1373  ++current_block;
1374 
1375  // break backwards is now old
1376  // break forward
1377  next_break_backward = next_break_forward + 1;
1378 
1379  // compute new break forward
1380  if (current_block < parent->block_indices.size())
1381  next_break_forward +=
1382  parent->block_indices.block_size(current_block);
1383  else
1384  // if we are beyond the end,
1385  // then move the next
1386  // boundary arbitrarily far
1387  // away
1388  next_break_forward = numbers::invalid_size_type;
1389  };
1390 
1391  ++global_index;
1392  }
1393 
1394 
1395 
1396  template <class BlockVectorType, bool Constness>
1397  void
1398  Iterator<BlockVectorType, Constness>::move_backward()
1399  {
1400  if (global_index != next_break_backward)
1401  --index_within_block;
1402  else if (current_block != 0)
1403  {
1404  // ok, we traverse a boundary
1405  // between blocks:
1406  --current_block;
1407  index_within_block =
1408  parent->block_indices.block_size(current_block) - 1;
1409 
1410  // break forwards is now old
1411  // break backward
1412  next_break_forward = next_break_backward - 1;
1413 
1414  // compute new break forward
1415  next_break_backward -=
1416  parent->block_indices.block_size(current_block);
1417  }
1418  else
1419  // current block was 0, we now
1420  // get into unspecified terrain
1421  {
1422  --current_block;
1423  index_within_block = numbers::invalid_size_type;
1424  next_break_forward = 0;
1425  next_break_backward = 0;
1426  };
1427 
1428  --global_index;
1429  }
1430 
1431 
1432  } // namespace BlockVectorIterators
1433 
1434 } // namespace internal
1435 
1436 
1437 
1438 template <class VectorType>
1439 inline std::size_t
1441 {
1442  return block_indices.total_size();
1443 }
1444 
1445 
1446 
1447 template <class VectorType>
1448 inline std::size_t
1450 {
1451  std::size_t local_size = 0;
1452  for (unsigned int b = 0; b < n_blocks(); ++b)
1453  local_size += block(b).locally_owned_size();
1454  return local_size;
1455 }
1456 
1457 
1458 
1459 template <class VectorType>
1460 inline IndexSet
1462 {
1463  IndexSet is(size());
1464 
1465  // copy index sets from blocks into the global one, shifted
1466  // by the appropriate amount for each block
1467  for (unsigned int b = 0; b < n_blocks(); ++b)
1468  {
1469  IndexSet x = block(b).locally_owned_elements();
1470  is.add_indices(x, block_indices.block_start(b));
1471  }
1472 
1473  is.compress();
1474 
1475  return is;
1476 }
1477 
1478 
1479 
1480 template <class VectorType>
1481 inline unsigned int
1483 {
1484  return block_indices.size();
1485 }
1486 
1487 
1488 template <class VectorType>
1490 BlockVectorBase<VectorType>::block(const unsigned int i)
1491 {
1492  AssertIndexRange(i, n_blocks());
1493 
1494  return components[i];
1495 }
1496 
1497 
1498 
1499 template <class VectorType>
1500 inline const typename BlockVectorBase<VectorType>::BlockType &
1501 BlockVectorBase<VectorType>::block(const unsigned int i) const
1502 {
1503  AssertIndexRange(i, n_blocks());
1504 
1505  return components[i];
1506 }
1507 
1508 
1509 
1510 template <class VectorType>
1511 inline const BlockIndices &
1513 {
1514  return block_indices;
1515 }
1516 
1517 
1518 template <class VectorType>
1519 inline void
1521 {
1522  std::vector<size_type> sizes(n_blocks());
1523 
1524  for (size_type i = 0; i < n_blocks(); ++i)
1525  sizes[i] = block(i).size();
1526 
1527  block_indices.reinit(sizes);
1528 }
1529 
1530 
1531 
1532 template <class VectorType>
1533 inline void
1535  ::VectorOperation::values operation)
1536 {
1537  for (unsigned int i = 0; i < n_blocks(); ++i)
1538  block(i).compress(operation);
1539 }
1540 
1541 
1542 
1543 template <class VectorType>
1546 {
1547  return iterator(*this, 0U);
1548 }
1549 
1550 
1551 
1552 template <class VectorType>
1555 {
1556  return const_iterator(*this, 0U);
1557 }
1558 
1559 
1560 template <class VectorType>
1563 {
1564  return iterator(*this, size());
1565 }
1566 
1567 
1568 
1569 template <class VectorType>
1572 {
1573  return const_iterator(*this, size());
1574 }
1575 
1576 
1577 template <class VectorType>
1578 inline bool
1580 {
1581  const std::pair<size_type, size_type> local_index =
1582  block_indices.global_to_local(global_index);
1583 
1584  return components[local_index.first].in_local_range(global_index);
1585 }
1586 
1587 
1588 template <class VectorType>
1589 bool
1591 {
1592  for (size_type i = 0; i < n_blocks(); ++i)
1593  if (components[i].all_zero() == false)
1594  return false;
1595 
1596  return true;
1597 }
1598 
1599 
1600 
1601 template <class VectorType>
1602 bool
1604 {
1605  for (size_type i = 0; i < n_blocks(); ++i)
1606  if (components[i].is_non_negative() == false)
1607  return false;
1608 
1609  return true;
1610 }
1611 
1612 
1613 
1614 template <class VectorType>
1617  const BlockVectorBase<VectorType> &v) const
1618 {
1619  Assert(n_blocks() == v.n_blocks(),
1621 
1622  value_type sum = 0.;
1623  for (size_type i = 0; i < n_blocks(); ++i)
1624  sum += components[i] * v.components[i];
1625 
1626  return sum;
1627 }
1628 
1629 
1630 template <class VectorType>
1633 {
1634  real_type sum = 0.;
1635  for (size_type i = 0; i < n_blocks(); ++i)
1636  sum += components[i].norm_sqr();
1637 
1638  return sum;
1639 }
1640 
1641 
1642 
1643 template <class VectorType>
1646 {
1647  value_type sum = 0.;
1648  // need to do static_cast as otherwise it won't work with
1649  // value_type=complex<T>
1650  for (size_type i = 0; i < n_blocks(); ++i)
1651  sum += components[i].mean_value() *
1653  components[i].size()));
1654 
1655  return sum / (typename numbers::NumberTraits<value_type>::real_type(size()));
1656 }
1657 
1658 
1659 
1660 template <class VectorType>
1663 {
1664  real_type sum = 0.;
1665  for (size_type i = 0; i < n_blocks(); ++i)
1666  sum += components[i].l1_norm();
1667 
1668  return sum;
1669 }
1670 
1671 
1672 
1673 template <class VectorType>
1676 {
1677  return std::sqrt(norm_sqr());
1678 }
1679 
1680 
1681 
1682 template <class VectorType>
1685 {
1686  real_type sum = 0.;
1687  for (size_type i = 0; i < n_blocks(); ++i)
1688  {
1689  value_type newval = components[i].linfty_norm();
1690  if (sum < newval)
1691  sum = newval;
1692  }
1693  return sum;
1694 }
1695 
1696 
1697 
1698 template <class VectorType>
1701  const typename BlockVectorBase<VectorType>::value_type a,
1702  const BlockVectorBase<VectorType> & V,
1703  const BlockVectorBase<VectorType> & W)
1704 {
1707 
1708  value_type sum = 0.;
1709  for (size_type i = 0; i < n_blocks(); ++i)
1710  sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1711 
1712  return sum;
1713 }
1714 
1715 
1716 
1717 template <class VectorType>
1720 {
1721  Assert(n_blocks() == v.n_blocks(),
1723 
1724  for (size_type i = 0; i < n_blocks(); ++i)
1725  {
1726  components[i] += v.components[i];
1727  }
1728 
1729  return *this;
1730 }
1731 
1732 
1733 
1734 template <class VectorType>
1737 {
1738  Assert(n_blocks() == v.n_blocks(),
1740 
1741  for (size_type i = 0; i < n_blocks(); ++i)
1742  {
1743  components[i] -= v.components[i];
1744  }
1745  return *this;
1746 }
1747 
1748 
1749 
1750 template <class VectorType>
1751 template <typename Number>
1752 inline void
1753 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1754  const std::vector<Number> & values)
1755 {
1756  Assert(indices.size() == values.size(),
1757  ExcDimensionMismatch(indices.size(), values.size()));
1758  add(indices.size(), indices.data(), values.data());
1759 }
1760 
1761 
1762 
1763 template <class VectorType>
1764 template <typename Number>
1765 inline void
1766 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1767  const Vector<Number> & values)
1768 {
1769  Assert(indices.size() == values.size(),
1770  ExcDimensionMismatch(indices.size(), values.size()));
1771  const size_type n_indices = indices.size();
1772  for (size_type i = 0; i < n_indices; ++i)
1773  (*this)(indices[i]) += values(i);
1774 }
1775 
1776 
1777 
1778 template <class VectorType>
1779 template <typename Number>
1780 inline void
1782  const size_type *indices,
1783  const Number * values)
1784 {
1785  for (size_type i = 0; i < n_indices; ++i)
1786  (*this)(indices[i]) += values[i];
1787 }
1788 
1789 
1790 
1791 template <class VectorType>
1792 void
1794 {
1795  AssertIsFinite(a);
1796 
1797  for (size_type i = 0; i < n_blocks(); ++i)
1798  {
1799  components[i].add(a);
1800  }
1801 }
1802 
1803 
1804 
1805 template <class VectorType>
1806 void
1808  const BlockVectorBase<VectorType> &v)
1809 {
1810  AssertIsFinite(a);
1811 
1812  Assert(n_blocks() == v.n_blocks(),
1814 
1815  for (size_type i = 0; i < n_blocks(); ++i)
1816  {
1817  components[i].add(a, v.components[i]);
1818  }
1819 }
1820 
1821 
1822 
1823 template <class VectorType>
1824 void
1826  const BlockVectorBase<VectorType> &v,
1827  const value_type b,
1829 {
1830  AssertIsFinite(a);
1831  AssertIsFinite(b);
1832 
1833  Assert(n_blocks() == v.n_blocks(),
1835  Assert(n_blocks() == w.n_blocks(),
1837 
1838 
1839  for (size_type i = 0; i < n_blocks(); ++i)
1840  {
1841  components[i].add(a, v.components[i], b, w.components[i]);
1842  }
1843 }
1844 
1845 
1846 
1847 template <class VectorType>
1848 void
1850  const BlockVectorBase<VectorType> &v)
1851 {
1852  AssertIsFinite(x);
1853 
1854  Assert(n_blocks() == v.n_blocks(),
1856 
1857  for (size_type i = 0; i < n_blocks(); ++i)
1858  {
1859  components[i].sadd(x, v.components[i]);
1860  }
1861 }
1862 
1863 
1864 
1865 template <class VectorType>
1866 void
1868  const value_type a,
1869  const BlockVectorBase<VectorType> &v)
1870 {
1871  AssertIsFinite(x);
1872  AssertIsFinite(a);
1873 
1874  Assert(n_blocks() == v.n_blocks(),
1876 
1877  for (size_type i = 0; i < n_blocks(); ++i)
1878  {
1879  components[i].sadd(x, a, v.components[i]);
1880  }
1881 }
1882 
1883 
1884 
1885 template <class VectorType>
1886 void
1888  const value_type a,
1889  const BlockVectorBase<VectorType> &v,
1890  const value_type b,
1891  const BlockVectorBase<VectorType> &w)
1892 {
1893  AssertIsFinite(x);
1894  AssertIsFinite(a);
1895  AssertIsFinite(b);
1896 
1897  Assert(n_blocks() == v.n_blocks(),
1899  Assert(n_blocks() == w.n_blocks(),
1901 
1902  for (size_type i = 0; i < n_blocks(); ++i)
1903  {
1904  components[i].sadd(x, a, v.components[i], b, w.components[i]);
1905  }
1906 }
1907 
1908 
1909 
1910 template <class VectorType>
1911 void
1913  const value_type a,
1914  const BlockVectorBase<VectorType> &v,
1915  const value_type b,
1916  const BlockVectorBase<VectorType> &w,
1917  const value_type c,
1918  const BlockVectorBase<VectorType> &y)
1919 {
1920  AssertIsFinite(x);
1921  AssertIsFinite(a);
1922  AssertIsFinite(b);
1923  AssertIsFinite(c);
1924 
1925  Assert(n_blocks() == v.n_blocks(),
1927  Assert(n_blocks() == w.n_blocks(),
1929  Assert(n_blocks() == y.n_blocks(),
1931 
1932  for (size_type i = 0; i < n_blocks(); ++i)
1933  {
1934  components[i].sadd(
1935  x, a, v.components[i], b, w.components[i], c, y.components[i]);
1936  }
1937 }
1938 
1939 
1940 
1941 template <class VectorType>
1942 template <class BlockVector2>
1943 void
1944 BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
1945 {
1946  Assert(n_blocks() == v.n_blocks(),
1947  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1948  for (size_type i = 0; i < n_blocks(); ++i)
1949  components[i].scale(v.block(i));
1950 }
1951 
1952 
1953 
1954 template <class VectorType>
1955 std::size_t
1957 {
1958  return (MemoryConsumption::memory_consumption(this->block_indices) +
1959  MemoryConsumption::memory_consumption(this->components));
1960 }
1961 
1962 
1963 
1964 template <class VectorType>
1965 template <class BlockVector2>
1966 void
1967 BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
1968 {
1969  AssertIsFinite(a);
1970 
1971  Assert(n_blocks() == v.n_blocks(),
1972  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1973 
1974  for (size_type i = 0; i < n_blocks(); ++i)
1975  components[i].equ(a, v.components[i]);
1976 }
1977 
1978 
1979 
1980 template <class VectorType>
1981 void
1983 {
1984  for (size_type i = 0; i < n_blocks(); ++i)
1985  block(i).update_ghost_values();
1986 }
1987 
1988 
1989 
1990 template <class VectorType>
1993 {
1994  AssertIsFinite(s);
1995 
1996  for (size_type i = 0; i < n_blocks(); ++i)
1997  components[i] = s;
1998 
1999  return *this;
2000 }
2001 
2002 
2003 template <class VectorType>
2006 {
2008 
2009  for (size_type i = 0; i < n_blocks(); ++i)
2010  components[i] = v.components[i];
2011 
2012  return *this;
2013 }
2014 
2015 
2016 template <class VectorType>
2017 template <class VectorType2>
2020 {
2022 
2023  for (size_type i = 0; i < n_blocks(); ++i)
2024  components[i] = v.components[i];
2025 
2026  return *this;
2027 }
2028 
2029 
2030 
2031 template <class VectorType>
2034 {
2035  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2036 
2037  size_type index_v = 0;
2038  for (size_type b = 0; b < n_blocks(); ++b)
2039  for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2040  block(b)(i) = v(index_v);
2041 
2042  return *this;
2043 }
2044 
2045 
2046 
2047 template <class VectorType>
2048 template <class VectorType2>
2049 inline bool
2051  const BlockVectorBase<VectorType2> &v) const
2052 {
2053  Assert(block_indices == v.block_indices, ExcDifferentBlockIndices());
2054 
2055  for (size_type i = 0; i < n_blocks(); ++i)
2056  if (!(components[i] == v.components[i]))
2057  return false;
2058 
2059  return true;
2060 }
2061 
2062 
2063 
2064 template <class VectorType>
2067 {
2068  AssertIsFinite(factor);
2069 
2070  for (size_type i = 0; i < n_blocks(); ++i)
2071  components[i] *= factor;
2072 
2073  return *this;
2074 }
2075 
2076 
2077 
2078 template <class VectorType>
2081 {
2082  AssertIsFinite(factor);
2083  Assert(factor != 0., ExcDivideByZero());
2084 
2085  const value_type inverse_factor = value_type(1.) / factor;
2086 
2087  for (size_type i = 0; i < n_blocks(); ++i)
2088  components[i] *= inverse_factor;
2089 
2090  return *this;
2091 }
2092 
2093 
2094 template <class VectorType>
2097 {
2098  const std::pair<unsigned int, size_type> local_index =
2099  block_indices.global_to_local(i);
2100  return components[local_index.first](local_index.second);
2101 }
2102 
2103 
2104 
2105 template <class VectorType>
2108 {
2109  const std::pair<unsigned int, size_type> local_index =
2110  block_indices.global_to_local(i);
2111  return components[local_index.first](local_index.second);
2112 }
2113 
2114 
2115 
2116 template <class VectorType>
2119 {
2120  return operator()(i);
2121 }
2122 
2123 
2124 
2125 template <class VectorType>
2128 {
2129  return operator()(i);
2130 }
2131 
2132 
2133 
2134 template <typename VectorType>
2135 template <typename OtherNumber>
2136 inline void
2138  const std::vector<size_type> &indices,
2139  std::vector<OtherNumber> & values) const
2140 {
2141  for (size_type i = 0; i < indices.size(); ++i)
2142  values[i] = operator()(indices[i]);
2143 }
2144 
2145 
2146 
2147 template <typename VectorType>
2148 template <typename ForwardIterator, typename OutputIterator>
2149 inline void
2151  ForwardIterator indices_begin,
2152  const ForwardIterator indices_end,
2153  OutputIterator values_begin) const
2154 {
2155  while (indices_begin != indices_end)
2156  {
2157  *values_begin = operator()(*indices_begin);
2158  indices_begin++;
2159  values_begin++;
2160  }
2161 }
2162 
2163 #endif // DOXYGEN
2164 
2166 
2167 #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:1653
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:2041
bool in_local_range(const size_type global_index) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1718
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1708
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:98
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:1461
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:487
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3025
real_type l1_norm() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
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:50
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:1647
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:400
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:2999
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:1744
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