Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
block_vector_base.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2018 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 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/memory_consumption.h>
24 #include <deal.II/base/numbers.h>
25 #include <deal.II/base/subscriptor.h>
26 
27 #include <deal.II/lac/block_indices.h>
28 #include <deal.II/lac/exceptions.h>
29 #include <deal.II/lac/vector.h>
30 #include <deal.II/lac/vector_operation.h>
31 
32 #include <cmath>
33 #include <cstddef>
34 #include <iterator>
35 #include <vector>
36 
37 DEAL_II_NAMESPACE_OPEN
38 
39 
44 template <typename>
46 
47 
64 template <typename VectorType>
66 {
67 private:
68  struct yes_type
69  {
70  char c[1];
71  };
72  struct no_type
73  {
74  char c[2];
75  };
76 
81  template <typename T>
82  static yes_type
84 
89  static no_type
91 
92 public:
98  static const bool value =
99  (sizeof(check_for_block_vector(static_cast<VectorType *>(nullptr))) ==
100  sizeof(yes_type));
101 };
102 
103 
104 // instantiation of the static member
105 template <typename VectorType>
107 
108 
109 
110 namespace internal
111 {
117  namespace BlockVectorIterators
118  {
136  template <class BlockVectorType, bool Constness>
137  class Iterator
138  {
139  public:
144 
150  using value_type =
151  typename std::conditional<Constness,
152  const typename BlockVectorType::value_type,
153  typename BlockVectorType::value_type>::type;
154 
160  using iterator_category = std::random_access_iterator_tag;
161  using difference_type = std::ptrdiff_t;
162  using reference = typename BlockVectorType::reference;
163  using pointer = value_type *;
164 
165  using dereference_type = typename std::conditional<
166  Constness,
167  value_type,
168  typename BlockVectorType::BlockType::reference>::type;
169 
174  using BlockVector = typename std::
175  conditional<Constness, const BlockVectorType, BlockVectorType>::type;
176 
186 
195 
196 
200  Iterator(const Iterator &c);
201 
202  private:
208  const size_type global_index,
209  const size_type current_block,
210  const size_type index_within_block,
212  const size_type next_break_backward);
213 
214  public:
218  Iterator &
219  operator=(const Iterator &c);
220 
226  dereference_type operator*() const;
227 
232  dereference_type operator[](const difference_type d) const;
233 
238  Iterator &
239  operator++();
240 
246  Iterator
247  operator++(int);
248 
253  Iterator &
254  operator--();
255 
261  Iterator
262  operator--(int);
263 
268  template <bool OtherConstness>
269  bool
271 
276  template <bool OtherConstness>
277  bool
279 
285  template <bool OtherConstness>
286  bool
287  operator<(const Iterator<BlockVectorType, OtherConstness> &i) const;
288 
292  template <bool OtherConstness>
293  bool
294  operator<=(const Iterator<BlockVectorType, OtherConstness> &i) const;
295 
299  template <bool OtherConstness>
300  bool
302 
306  template <bool OtherConstness>
307  bool
309 
313  template <bool OtherConstness>
314  difference_type
316 
321  Iterator
322  operator+(const difference_type &d) const;
323 
328  Iterator
329  operator-(const difference_type &d) const;
330 
335  Iterator &
336  operator+=(const difference_type &d);
337 
342  Iterator &
343  operator-=(const difference_type &d);
344 
355  "Your program tried to compare iterators pointing to "
356  "different block vectors. There is no reasonable way "
357  "to do this.");
358 
360  private:
367 
372 
377  unsigned int current_block;
378  size_type index_within_block;
379 
387  size_type next_break_backward;
388 
392  void
393  move_forward();
394 
398  void
399  move_backward();
400 
401 
405  template <typename, bool>
406  friend class Iterator;
407  };
408  } // namespace BlockVectorIterators
409 } // namespace internal
410 
411 
450 template <class VectorType>
451 class BlockVectorBase : public Subscriptor
452 {
453 public:
457  using BlockType = VectorType;
458 
459  /*
460  * Declare standard types used in
461  * all containers. These types
462  * parallel those in the
463  * <tt>C++</tt> standard
464  * libraries
465  * <tt>std::vector<...></tt>
466  * class. This includes iterator
467  * types.
468  */
469  using value_type = typename BlockType::value_type;
470  using pointer = value_type *;
471  using const_pointer = const value_type *;
472  using iterator =
474  using const_iterator =
476  using reference = typename BlockType::reference;
477  using const_reference = typename BlockType::const_reference;
478  using size_type = types::global_dof_index;
479 
489  using real_type = typename BlockType::real_type;
490 
494  BlockVectorBase() = default;
495 
499  BlockVectorBase(const BlockVectorBase & /*V*/) = default;
500 
506  BlockVectorBase(BlockVectorBase && /*V*/) noexcept = default;
507 
514  void
515  collect_sizes();
516 
527  void
528  compress(::VectorOperation::values operation);
529 
533  BlockType &
534  block(const unsigned int i);
535 
539  const BlockType &
540  block(const unsigned int i) const;
541 
547  const BlockIndices &
548  get_block_indices() const;
549 
553  unsigned int
554  n_blocks() const;
555 
560  std::size_t
561  size() const;
562 
579  IndexSet
580  locally_owned_elements() const;
581 
585  iterator
586  begin();
587 
593  begin() const;
594 
598  iterator
599  end();
600 
606  end() const;
607 
611  value_type
612  operator()(const size_type i) const;
613 
617  reference
618  operator()(const size_type i);
619 
625  value_type operator[](const size_type i) const;
626 
632  reference operator[](const size_type i);
633 
649  template <typename OtherNumber>
650  void
651  extract_subvector_to(const std::vector<size_type> &indices,
652  std::vector<OtherNumber> & values) const;
653 
681  template <typename ForwardIterator, typename OutputIterator>
682  void
683  extract_subvector_to(ForwardIterator indices_begin,
684  const ForwardIterator indices_end,
685  OutputIterator values_begin) const;
686 
692  operator=(const value_type s);
693 
698  operator=(const BlockVectorBase &V);
699 
706  operator=(BlockVectorBase && /*V*/) = default; // NOLINT
707 
711  template <class VectorType2>
713  operator=(const BlockVectorBase<VectorType2> &V);
714 
719  operator=(const VectorType &v);
720 
725  template <class VectorType2>
726  bool
727  operator==(const BlockVectorBase<VectorType2> &v) const;
728 
732  value_type operator*(const BlockVectorBase &V) const;
733 
737  real_type
738  norm_sqr() const;
739 
743  value_type
744  mean_value() const;
745 
749  real_type
750  l1_norm() const;
751 
756  real_type
757  l2_norm() const;
758 
763  real_type
764  linfty_norm() const;
765 
788  value_type
789  add_and_dot(const value_type a,
790  const BlockVectorBase &V,
791  const BlockVectorBase &W);
792 
797  bool
798  in_local_range(const size_type global_index) const;
799 
805  bool
806  all_zero() const;
807 
813  bool
814  is_non_negative() const;
815 
820  operator+=(const BlockVectorBase &V);
821 
826  operator-=(const BlockVectorBase &V);
827 
828 
833  template <typename Number>
834  void
835  add(const std::vector<size_type> &indices, const std::vector<Number> &values);
836 
841  template <typename Number>
842  void
843  add(const std::vector<size_type> &indices, const Vector<Number> &values);
844 
850  template <typename Number>
851  void
852  add(const size_type n_elements,
853  const size_type *indices,
854  const Number * values);
855 
860  void
861  add(const value_type s);
862 
866  void
867  add(const value_type a, const BlockVectorBase &V);
868 
872  void
873  add(const value_type a,
874  const BlockVectorBase &V,
875  const value_type b,
876  const BlockVectorBase &W);
877 
881  void
882  sadd(const value_type s, const BlockVectorBase &V);
883 
887  void
888  sadd(const value_type s, const value_type a, const BlockVectorBase &V);
889 
893  void
894  sadd(const value_type s,
895  const value_type a,
896  const BlockVectorBase &V,
897  const value_type b,
898  const BlockVectorBase &W);
899 
903  void
904  sadd(const value_type s,
905  const value_type a,
906  const BlockVectorBase &V,
907  const value_type b,
908  const BlockVectorBase &W,
909  const value_type c,
910  const BlockVectorBase &X);
911 
916  operator*=(const value_type factor);
917 
922  operator/=(const value_type factor);
923 
928  template <class BlockVector2>
929  void
930  scale(const BlockVector2 &v);
931 
935  template <class BlockVector2>
936  void
937  equ(const value_type a, const BlockVector2 &V);
938 
942  void
943  equ(const value_type a,
944  const BlockVectorBase &V,
945  const value_type b,
946  const BlockVectorBase &W);
947 
952  void
953  update_ghost_values() const;
954 
959  std::size_t
960  memory_consumption() const;
961 
962 protected:
966  std::vector<VectorType> components;
967 
973 
977  template <typename N, bool C>
978  friend class ::internal::BlockVectorIterators::Iterator;
979 
980  template <typename>
981  friend class BlockVectorBase;
982 };
983 
984 
987 /*----------------------- Inline functions ----------------------------------*/
988 
989 
990 #ifndef DOXYGEN
991 namespace internal
992 {
993  namespace BlockVectorIterators
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 
1006 
1007 
1008  template <class BlockVectorType, bool Constness>
1009  inline Iterator<BlockVectorType, Constness>::Iterator(
1010  const Iterator<BlockVectorType, !Constness> &c)
1011  : parent(c.parent)
1012  , global_index(c.global_index)
1013  , current_block(c.current_block)
1014  , index_within_block(c.index_within_block)
1015  , next_break_forward(c.next_break_forward)
1016  , next_break_backward(c.next_break_backward)
1017  {
1018  // Only permit copy-constructing const iterators from non-const
1019  // iterators, and not vice versa (i.e., Constness must always be
1020  // true).
1021  static_assert(Constness == true,
1022  "Constructing a non-const iterator from a const iterator "
1023  "does not make sense.");
1024  }
1025 
1026 
1027 
1028  template <class BlockVectorType, bool Constness>
1029  inline Iterator<BlockVectorType, Constness>::Iterator(
1030  BlockVector & parent,
1031  const size_type global_index,
1032  const size_type current_block,
1033  const size_type index_within_block,
1034  const size_type next_break_forward,
1035  const size_type next_break_backward)
1036  : parent(&parent)
1037  , global_index(global_index)
1038  , current_block(current_block)
1039  , index_within_block(index_within_block)
1040  , next_break_forward(next_break_forward)
1041  , next_break_backward(next_break_backward)
1042  {}
1043 
1044 
1045 
1046  template <class BlockVectorType, bool Constness>
1047  inline Iterator<BlockVectorType, Constness> &
1048  Iterator<BlockVectorType, Constness>::operator=(const Iterator &c)
1049  {
1050  parent = c.parent;
1051  global_index = c.global_index;
1052  index_within_block = c.index_within_block;
1053  current_block = c.current_block;
1054  next_break_forward = c.next_break_forward;
1055  next_break_backward = c.next_break_backward;
1056 
1057  return *this;
1058  }
1059 
1060 
1061 
1062  template <class BlockVectorType, bool Constness>
1063  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1065  {
1066  return parent->block(current_block)(index_within_block);
1067  }
1068 
1069 
1070 
1071  template <class BlockVectorType, bool Constness>
1072  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1073  Iterator<BlockVectorType, Constness>::
1074  operator[](const difference_type d) const
1075  {
1076  // if the index pointed to is
1077  // still within the block we
1078  // currently point into, then we
1079  // can save the computation of
1080  // the block
1081  if ((global_index + d >= next_break_backward) &&
1082  (global_index + d <= next_break_forward))
1083  return parent->block(current_block)(index_within_block + d);
1084 
1085  // if the index is not within the
1086  // block of the block vector into
1087  // which we presently point, then
1088  // there is no way: we have to
1089  // search for the block. this can
1090  // be done through the parent
1091  // class as well.
1092  return (*parent)(global_index + d);
1093  }
1094 
1095 
1096 
1097  template <class BlockVectorType, bool Constness>
1098  inline Iterator<BlockVectorType, Constness> &
1099  Iterator<BlockVectorType, Constness>::operator++()
1100  {
1101  move_forward();
1102  return *this;
1103  }
1104 
1105 
1106 
1107  template <class BlockVectorType, bool Constness>
1108  inline Iterator<BlockVectorType, Constness>
1109  Iterator<BlockVectorType, Constness>::operator++(int)
1110  {
1111  const Iterator old_value = *this;
1112  move_forward();
1113  return old_value;
1114  }
1115 
1116 
1117 
1118  template <class BlockVectorType, bool Constness>
1119  inline Iterator<BlockVectorType, Constness> &
1120  Iterator<BlockVectorType, Constness>::operator--()
1121  {
1122  move_backward();
1123  return *this;
1124  }
1125 
1126 
1127 
1128  template <class BlockVectorType, bool Constness>
1129  inline Iterator<BlockVectorType, Constness>
1130  Iterator<BlockVectorType, Constness>::operator--(int)
1131  {
1132  const Iterator old_value = *this;
1133  move_backward();
1134  return old_value;
1135  }
1136 
1137 
1138 
1139  template <class BlockVectorType, bool Constness>
1140  template <bool OtherConstness>
1141  inline bool
1142  Iterator<BlockVectorType, Constness>::
1143  operator==(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
1155  Iterator<BlockVectorType, Constness>::
1156  operator!=(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
1168  Iterator<BlockVectorType, Constness>::
1169  operator<(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
1181  Iterator<BlockVectorType, Constness>::
1182  operator<=(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
1194  Iterator<BlockVectorType, Constness>::
1195  operator>(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 bool
1207  Iterator<BlockVectorType, Constness>::
1208  operator>=(const Iterator<BlockVectorType, OtherConstness> &i) const
1209  {
1210  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1211 
1212  return (global_index >= i.global_index);
1213  }
1214 
1215 
1216 
1217  template <class BlockVectorType, bool Constness>
1218  template <bool OtherConstness>
1219  inline typename Iterator<BlockVectorType, Constness>::difference_type
1221  operator-(const Iterator<BlockVectorType, OtherConstness> &i) const
1222  {
1223  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1224 
1225  return (static_cast<signed int>(global_index) -
1226  static_cast<signed int>(i.global_index));
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>
1261  operator-(const difference_type &d) const
1262  {
1263  // if the index pointed to is
1264  // still within the block we
1265  // currently point into, then we
1266  // can save the computation of
1267  // the block
1268  if ((global_index - d >= next_break_backward) &&
1269  (global_index - d <= next_break_forward))
1270  return Iterator(*parent,
1271  global_index - d,
1272  current_block,
1273  index_within_block - d,
1274  next_break_forward,
1275  next_break_backward);
1276  else
1277  // outside present block, so
1278  // have to seek new block
1279  // anyway
1280  return Iterator(*parent, global_index - d);
1281  }
1282 
1283 
1284 
1285  template <class BlockVectorType, bool Constness>
1286  inline Iterator<BlockVectorType, Constness> &
1287  Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1288  {
1289  // if the index pointed to is
1290  // still within the block we
1291  // currently point into, then we
1292  // can save the computation of
1293  // the block
1294  if ((global_index + d >= next_break_backward) &&
1295  (global_index + d <= next_break_forward))
1296  {
1297  global_index += d;
1298  index_within_block += d;
1299  }
1300  else
1301  // outside present block, so
1302  // have to seek new block
1303  // anyway
1304  *this = Iterator(*parent, global_index + d);
1305 
1306  return *this;
1307  }
1308 
1309 
1310 
1311  template <class BlockVectorType, bool Constness>
1312  inline Iterator<BlockVectorType, Constness> &
1313  Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1314  {
1315  // if the index pointed to is
1316  // still within the block we
1317  // currently point into, then we
1318  // can save the computation of
1319  // the block
1320  if ((global_index - d >= next_break_backward) &&
1321  (global_index - d <= next_break_forward))
1322  {
1323  global_index -= d;
1324  index_within_block -= d;
1325  }
1326  else
1327  // outside present block, so
1328  // have to seek new block
1329  // anyway
1330  *this = Iterator(*parent, global_index - d);
1331 
1332  return *this;
1333  }
1334 
1335 
1336  template <class BlockVectorType, bool Constness>
1337  Iterator<BlockVectorType, Constness>::Iterator(BlockVector & parent,
1338  const size_type global_index)
1339  : parent(&parent)
1340  , global_index(global_index)
1341  {
1342  // find which block we are
1343  // in. for this, take into
1344  // account that it happens at
1345  // times that people want to
1346  // initialize iterators
1347  // past-the-end
1348  if (global_index < parent.size())
1349  {
1350  const std::pair<size_type, size_type> indices =
1351  parent.block_indices.global_to_local(global_index);
1352  current_block = indices.first;
1353  index_within_block = indices.second;
1354 
1355  next_break_backward =
1356  parent.block_indices.local_to_global(current_block, 0);
1357  next_break_forward =
1358  (parent.block_indices.local_to_global(current_block, 0) +
1359  parent.block_indices.block_size(current_block) - 1);
1360  }
1361  else
1362  // past the end. only have one
1363  // value for this
1364  {
1365  this->global_index = parent.size();
1366  current_block = parent.n_blocks();
1367  index_within_block = 0;
1368  next_break_backward = global_index;
1369  next_break_forward = numbers::invalid_size_type;
1370  };
1371  }
1372 
1373 
1374 
1375  template <class BlockVectorType, bool Constness>
1376  void
1377  Iterator<BlockVectorType, Constness>::move_forward()
1378  {
1379  if (global_index != next_break_forward)
1380  ++index_within_block;
1381  else
1382  {
1383  // ok, we traverse a boundary
1384  // between blocks:
1385  index_within_block = 0;
1386  ++current_block;
1387 
1388  // break backwards is now old
1389  // break forward
1390  next_break_backward = next_break_forward + 1;
1391 
1392  // compute new break forward
1393  if (current_block < parent->block_indices.size())
1394  next_break_forward +=
1395  parent->block_indices.block_size(current_block);
1396  else
1397  // if we are beyond the end,
1398  // then move the next
1399  // boundary arbitrarily far
1400  // away
1401  next_break_forward = numbers::invalid_size_type;
1402  };
1403 
1404  ++global_index;
1405  }
1406 
1407 
1408 
1409  template <class BlockVectorType, bool Constness>
1410  void
1411  Iterator<BlockVectorType, Constness>::move_backward()
1412  {
1413  if (global_index != next_break_backward)
1414  --index_within_block;
1415  else if (current_block != 0)
1416  {
1417  // ok, we traverse a boundary
1418  // between blocks:
1419  --current_block;
1420  index_within_block =
1421  parent->block_indices.block_size(current_block) - 1;
1422 
1423  // break forwards is now old
1424  // break backward
1425  next_break_forward = next_break_backward - 1;
1426 
1427  // compute new break forward
1428  next_break_backward -=
1429  parent->block_indices.block_size(current_block);
1430  }
1431  else
1432  // current block was 0, we now
1433  // get into unspecified terrain
1434  {
1435  --current_block;
1436  index_within_block = numbers::invalid_size_type;
1437  next_break_forward = 0;
1438  next_break_backward = 0;
1439  };
1440 
1441  --global_index;
1442  }
1443 
1444 
1445  } // namespace BlockVectorIterators
1446 
1447 } // namespace internal
1448 
1449 
1450 
1451 template <class VectorType>
1452 inline std::size_t
1454 {
1455  return block_indices.total_size();
1456 }
1457 
1458 
1459 
1460 template <class VectorType>
1461 inline IndexSet
1463 {
1464  IndexSet is(size());
1465 
1466  // copy index sets from blocks into the global one, shifted
1467  // by the appropriate amount for each block
1468  for (unsigned int b = 0; b < n_blocks(); ++b)
1469  {
1470  IndexSet x = block(b).locally_owned_elements();
1471  is.add_indices(x, block_indices.block_start(b));
1472  }
1473 
1474  is.compress();
1475 
1476  return is;
1477 }
1478 
1479 
1480 
1481 template <class VectorType>
1482 inline unsigned int
1484 {
1485  return block_indices.size();
1486 }
1487 
1488 
1489 template <class VectorType>
1491 BlockVectorBase<VectorType>::block(const unsigned int i)
1492 {
1493  Assert(i < n_blocks(), ExcIndexRange(i, 0, n_blocks()));
1494 
1495  return components[i];
1496 }
1497 
1498 
1499 
1500 template <class VectorType>
1501 inline const typename BlockVectorBase<VectorType>::BlockType &
1502 BlockVectorBase<VectorType>::block(const unsigned int i) const
1503 {
1504  Assert(i < n_blocks(), ExcIndexRange(i, 0, n_blocks()));
1505 
1506  return components[i];
1507 }
1508 
1509 
1510 
1511 template <class VectorType>
1512 inline const BlockIndices &
1514 {
1515  return block_indices;
1516 }
1517 
1518 
1519 template <class VectorType>
1520 inline void
1522 {
1523  std::vector<size_type> sizes(n_blocks());
1524 
1525  for (size_type i = 0; i < n_blocks(); ++i)
1526  sizes[i] = block(i).size();
1527 
1528  block_indices.reinit(sizes);
1529 }
1530 
1531 
1532 
1533 template <class VectorType>
1534 inline void
1536  ::VectorOperation::values operation)
1537 {
1538  for (unsigned int i = 0; i < n_blocks(); ++i)
1539  block(i).compress(operation);
1540 }
1541 
1542 
1543 
1544 template <class VectorType>
1547 {
1548  return iterator(*this, 0U);
1549 }
1550 
1551 
1552 
1553 template <class VectorType>
1556 {
1557  return const_iterator(*this, 0U);
1558 }
1559 
1560 
1561 template <class VectorType>
1564 {
1565  return iterator(*this, size());
1566 }
1567 
1568 
1569 
1570 template <class VectorType>
1573 {
1574  return const_iterator(*this, size());
1575 }
1576 
1577 
1578 template <class VectorType>
1579 inline bool
1580 BlockVectorBase<VectorType>::in_local_range(const size_type global_index) const
1581 {
1582  const std::pair<size_type, size_type> local_index =
1583  block_indices.global_to_local(global_index);
1584 
1585  return components[local_index.first].in_local_range(global_index);
1586 }
1587 
1588 
1589 template <class VectorType>
1590 bool
1592 {
1593  for (size_type i = 0; i < n_blocks(); ++i)
1594  if (components[i].all_zero() == false)
1595  return false;
1596 
1597  return true;
1598 }
1599 
1600 
1601 
1602 template <class VectorType>
1603 bool
1605 {
1606  for (size_type i = 0; i < n_blocks(); ++i)
1607  if (components[i].is_non_negative() == false)
1608  return false;
1609 
1610  return true;
1611 }
1612 
1613 
1614 
1615 template <class VectorType>
1616 typename BlockVectorBase<VectorType>::value_type BlockVectorBase<VectorType>::
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>
1644 typename BlockVectorBase<VectorType>::value_type
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>
1699 typename BlockVectorBase<VectorType>::value_type
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
1781 BlockVectorBase<VectorType>::add(const size_type n_indices,
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
1793 BlockVectorBase<VectorType>::add(const value_type a)
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
1807 BlockVectorBase<VectorType>::add(const value_type a,
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
1825 BlockVectorBase<VectorType>::add(const value_type a,
1826  const BlockVectorBase<VectorType> &v,
1827  const value_type b,
1828  const BlockVectorBase<VectorType> &w)
1829 {
1830  AssertIsFinite(a);
1831  AssertIsFinite(b);
1832 
1833  Assert(n_blocks() == v.n_blocks(),
1835  Assert(n_blocks() == w.n_blocks(),
1836  ExcDimensionMismatch(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
1849 BlockVectorBase<VectorType>::sadd(const value_type x,
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
1867 BlockVectorBase<VectorType>::sadd(const value_type x,
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
1887 BlockVectorBase<VectorType>::sadd(const value_type x,
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(),
1900  ExcDimensionMismatch(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
1912 BlockVectorBase<VectorType>::sadd(const value_type x,
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(),
1928  ExcDimensionMismatch(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 void
1956 BlockVectorBase<VectorType>::equ(const value_type a,
1957  const BlockVectorBase<VectorType> &v,
1958  const value_type b,
1959  const BlockVectorBase<VectorType> &w)
1960 {
1961  AssertIsFinite(a);
1962  AssertIsFinite(b);
1963 
1964  Assert(n_blocks() == v.n_blocks(),
1966  Assert(n_blocks() == w.n_blocks(),
1967  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1968 
1969  for (size_type i = 0; i < n_blocks(); ++i)
1970  {
1971  components[i].equ(a, v.components[i], b, w.components[i]);
1972  }
1973 }
1974 
1975 
1976 
1977 template <class VectorType>
1978 std::size_t
1980 {
1981  return (MemoryConsumption::memory_consumption(this->block_indices) +
1982  MemoryConsumption::memory_consumption(this->components));
1983 }
1984 
1985 
1986 
1987 template <class VectorType>
1988 template <class BlockVector2>
1989 void
1990 BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
1991 {
1992  AssertIsFinite(a);
1993 
1994  Assert(n_blocks() == v.n_blocks(),
1995  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1996 
1997  for (size_type i = 0; i < n_blocks(); ++i)
1998  components[i].equ(a, v.components[i]);
1999 }
2000 
2001 
2002 
2003 template <class VectorType>
2004 void
2006 {
2007  for (size_type i = 0; i < n_blocks(); ++i)
2008  block(i).update_ghost_values();
2009 }
2010 
2011 
2012 
2013 template <class VectorType>
2015 BlockVectorBase<VectorType>::operator=(const value_type s)
2016 {
2017  AssertIsFinite(s);
2018 
2019  for (size_type i = 0; i < n_blocks(); ++i)
2020  components[i] = s;
2021 
2022  return *this;
2023 }
2024 
2025 
2026 template <class VectorType>
2029 {
2031 
2032  for (size_type i = 0; i < n_blocks(); ++i)
2033  components[i] = v.components[i];
2034 
2035  return *this;
2036 }
2037 
2038 
2039 template <class VectorType>
2040 template <class VectorType2>
2043 {
2045 
2046  for (size_type i = 0; i < n_blocks(); ++i)
2047  components[i] = v.components[i];
2048 
2049  return *this;
2050 }
2051 
2052 
2053 
2054 template <class VectorType>
2056 BlockVectorBase<VectorType>::operator=(const VectorType &v)
2057 {
2058  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2059 
2060  size_type index_v = 0;
2061  for (size_type b = 0; b < n_blocks(); ++b)
2062  for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2063  block(b)(i) = v(index_v);
2064 
2065  return *this;
2066 }
2067 
2068 
2069 
2070 template <class VectorType>
2071 template <class VectorType2>
2072 inline bool
2075 {
2077 
2078  for (size_type i = 0; i < n_blocks(); ++i)
2079  if (!(components[i] == v.components[i]))
2080  return false;
2081 
2082  return true;
2083 }
2084 
2085 
2086 
2087 template <class VectorType>
2089 BlockVectorBase<VectorType>::operator*=(const value_type factor)
2090 {
2091  AssertIsFinite(factor);
2092 
2093  for (size_type i = 0; i < n_blocks(); ++i)
2094  components[i] *= factor;
2095 
2096  return *this;
2097 }
2098 
2099 
2100 
2101 template <class VectorType>
2103 BlockVectorBase<VectorType>::operator/=(const value_type factor)
2104 {
2105  AssertIsFinite(factor);
2106  Assert(factor != 0., ExcDivideByZero());
2107 
2108  for (size_type i = 0; i < n_blocks(); ++i)
2109  components[i] /= factor;
2110 
2111  return *this;
2112 }
2113 
2114 
2115 template <class VectorType>
2116 inline typename BlockVectorBase<VectorType>::value_type
2117 BlockVectorBase<VectorType>::operator()(const size_type i) const
2118 {
2119  const std::pair<unsigned int, size_type> local_index =
2121  return components[local_index.first](local_index.second);
2122 }
2123 
2124 
2125 
2126 template <class VectorType>
2127 inline typename BlockVectorBase<VectorType>::reference
2128 BlockVectorBase<VectorType>::operator()(const size_type i)
2129 {
2130  const std::pair<unsigned int, size_type> local_index =
2132  return components[local_index.first](local_index.second);
2133 }
2134 
2135 
2136 
2137 template <class VectorType>
2138 inline typename BlockVectorBase<VectorType>::value_type
2139  BlockVectorBase<VectorType>::operator[](const size_type i) const
2140 {
2141  return operator()(i);
2142 }
2143 
2144 
2145 
2146 template <class VectorType>
2147 inline typename BlockVectorBase<VectorType>::reference
2148  BlockVectorBase<VectorType>::operator[](const size_type i)
2149 {
2150  return operator()(i);
2151 }
2152 
2153 
2154 
2155 template <typename VectorType>
2156 template <typename OtherNumber>
2157 inline void
2159  const std::vector<size_type> &indices,
2160  std::vector<OtherNumber> & values) const
2161 {
2162  for (size_type i = 0; i < indices.size(); ++i)
2163  values[i] = operator()(indices[i]);
2164 }
2165 
2166 
2167 
2168 template <typename VectorType>
2169 template <typename ForwardIterator, typename OutputIterator>
2170 inline void
2172  ForwardIterator indices_begin,
2173  const ForwardIterator indices_end,
2174  OutputIterator values_begin) const
2175 {
2176  while (indices_begin != indices_end)
2177  {
2178  *values_begin = operator()(*indices_begin);
2179  indices_begin++;
2180  values_begin++;
2181  }
2182 }
2183 
2184 #endif // DOXYGEN
2185 
2186 DEAL_II_NAMESPACE_CLOSE
2187 
2188 #endif
const types::global_dof_index invalid_size_type
Definition: types.h:182
static yes_type check_for_block_vector(const BlockVectorBase< T > *)
bool operator>(const Iterator< BlockVectorType, OtherConstness > &i) const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
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
dereference_type operator[](const difference_type d) const
bool operator==(const BlockVectorBase< VectorType2 > &v) const
void compress(::VectorOperation::values operation)
bool in_local_range(const size_type global_index) const
void collect_sizes()
STL namespace.
size_type block_size(const unsigned int i) const
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
static ::ExceptionBase & ExcDivideByZero()
Iterator operator+(const difference_type &d) const
LinearAlgebra::distributed::Vector< Number > Vector
typename std::conditional< Constness, const typename BlockVectorType::value_type, typename BlockVectorType::value_type >::type value_type
BlockVectorBase & operator/=(const value_type factor)
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)
BlockVectorBase()=default
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
T sum(const T &t, const MPI_Comm &mpi_communicator)
const BlockIndices & get_block_indices() const
#define Assert(cond, exc)
Definition: exceptions.h:1407
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:496
real_type l1_norm() const
void update_ghost_values() const
size_type local_to_global(const unsigned int block, const size_type index) const
bool operator==(const Iterator< BlockVectorType, OtherConstness > &i) const
static ::ExceptionBase & ExcDifferentBlockIndices()
static const bool value
BlockVectorBase & operator+=(const BlockVectorBase &V)
std::vector< VectorType > components
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)
Iterator & operator+=(const difference_type &d)
unsigned int global_dof_index
Definition: types.h:89
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
unsigned int n_blocks() const
BlockVectorBase & operator-=(const BlockVectorBase &V)
typename BlockType::real_type real_type
iterator end()
bool operator!=(const Iterator< BlockVectorType, OtherConstness > &i) const
size_type block_start(const unsigned int i) const
value_type operator()(const size_type i) const
difference_type operator-(const Iterator< BlockVectorType, OtherConstness > &i) const
real_type linfty_norm() const
bool operator>=(const Iterator< BlockVectorType, OtherConstness > &i) const
unsigned int size() const
BlockType & block(const unsigned int i)
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
iterator begin()
BlockVectorBase & operator=(const value_type s)
static ::ExceptionBase & ExcPointerToDifferentVectors()
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
#define AssertIsFinite(number)
Definition: exceptions.h:1669
void equ(const value_type a, const BlockVector2 &V)
Iterator & operator=(const Iterator &c)
Iterator & operator-=(const difference_type &d)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
dereference_type operator*() const
value_type mean_value() const