Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\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 // Forward declaration
45 #ifndef DOXYGEN
46 template <typename>
47 class BlockVectorBase;
48 #endif
49 
66 template <typename VectorType>
68 {
69 private:
70  struct yes_type
71  {
72  char c[1];
73  };
74  struct no_type
75  {
76  char c[2];
77  };
78 
83  template <typename T>
84  static yes_type
86 
91  static no_type
93 
94 public:
100  static const bool value =
101  (sizeof(check_for_block_vector(static_cast<VectorType *>(nullptr))) ==
102  sizeof(yes_type));
103 };
104 
105 
106 // instantiation of the static member
107 template <typename VectorType>
109 
110 
111 
112 namespace internal
113 {
119  namespace BlockVectorIterators
120  {
138  template <class BlockVectorType, bool Constness>
139  class Iterator
140  {
141  public:
146 
152  using value_type =
153  typename std::conditional<Constness,
154  const typename BlockVectorType::value_type,
155  typename BlockVectorType::value_type>::type;
156 
162  using iterator_category = std::random_access_iterator_tag;
163  using difference_type = std::ptrdiff_t;
164  using reference = typename BlockVectorType::reference;
165  using pointer = value_type *;
166 
167  using dereference_type = typename std::conditional<
168  Constness,
169  value_type,
170  typename BlockVectorType::BlockType::reference>::type;
171 
176  using BlockVector = typename std::
177  conditional<Constness, const BlockVectorType, BlockVectorType>::type;
178 
187  Iterator(BlockVector &parent, const size_type global_index);
188 
197 
198 
202  Iterator(const Iterator &c);
203 
204  private:
209  Iterator(BlockVector & parent,
210  const size_type global_index,
211  const size_type current_block,
212  const size_type index_within_block,
213  const size_type next_break_forward,
214  const size_type next_break_backward);
215 
216  public:
220  Iterator &
221  operator=(const Iterator &c);
222 
228  dereference_type operator*() const;
229 
234  dereference_type operator[](const difference_type d) const;
235 
240  Iterator &
241  operator++();
242 
248  Iterator
249  operator++(int);
250 
255  Iterator &
256  operator--();
257 
263  Iterator
264  operator--(int);
265 
270  template <bool OtherConstness>
271  bool
272  operator==(const Iterator<BlockVectorType, OtherConstness> &i) const;
273 
278  template <bool OtherConstness>
279  bool
280  operator!=(const Iterator<BlockVectorType, OtherConstness> &i) const;
281 
287  template <bool OtherConstness>
288  bool
289  operator<(const Iterator<BlockVectorType, OtherConstness> &i) const;
290 
294  template <bool OtherConstness>
295  bool
296  operator<=(const Iterator<BlockVectorType, OtherConstness> &i) const;
297 
301  template <bool OtherConstness>
302  bool
303  operator>(const Iterator<BlockVectorType, OtherConstness> &i) const;
304 
308  template <bool OtherConstness>
309  bool
310  operator>=(const Iterator<BlockVectorType, OtherConstness> &i) const;
311 
315  template <bool OtherConstness>
316  difference_type
318 
323  Iterator
324  operator+(const difference_type &d) const;
325 
330  Iterator
331  operator-(const difference_type &d) const;
332 
337  Iterator &
338  operator+=(const difference_type &d);
339 
344  Iterator &
345  operator-=(const difference_type &d);
346 
356  DeclExceptionMsg(ExcPointerToDifferentVectors,
357  "Your program tried to compare iterators pointing to "
358  "different block vectors. There is no reasonable way "
359  "to do this.");
360 
362  private:
369 
374 
379  unsigned int current_block;
380  size_type index_within_block;
381 
389  size_type next_break_backward;
390 
394  void
395  move_forward();
396 
400  void
401  move_backward();
402 
403  // Mark all other instances of this template as friends.
404  template <typename, bool>
405  friend class Iterator;
406  };
407  } // namespace BlockVectorIterators
408 } // namespace internal
409 
410 
449 template <class VectorType>
451 {
452 public:
456  using BlockType = VectorType;
457 
458  /*
459  * Declare standard types used in
460  * all containers. These types
461  * parallel those in the
462  * <tt>C++</tt> standard
463  * libraries
464  * <tt>std::vector<...></tt>
465  * class. This includes iterator
466  * types.
467  */
468  using value_type = typename BlockType::value_type;
469  using pointer = value_type *;
470  using const_pointer = const value_type *;
471  using iterator =
473  using const_iterator =
475  using reference = typename BlockType::reference;
476  using const_reference = typename BlockType::const_reference;
477  using size_type = types::global_dof_index;
478 
488  using real_type = typename BlockType::real_type;
489 
493  BlockVectorBase() = default;
494 
498  BlockVectorBase(const BlockVectorBase & /*V*/) = default;
499 
505  BlockVectorBase(BlockVectorBase && /*V*/) noexcept = default;
506 
513  void
514  collect_sizes();
515 
526  void
527  compress(::VectorOperation::values operation);
528 
532  BlockType &
533  block(const unsigned int i);
534 
538  const BlockType &
539  block(const unsigned int i) const;
540 
546  const BlockIndices &
547  get_block_indices() const;
548 
552  unsigned int
553  n_blocks() const;
554 
559  std::size_t
560  size() const;
561 
578  IndexSet
579  locally_owned_elements() const;
580 
584  iterator
585  begin();
586 
592  begin() const;
593 
597  iterator
598  end();
599 
605  end() const;
606 
610  value_type
611  operator()(const size_type i) const;
612 
616  reference
617  operator()(const size_type i);
618 
624  value_type operator[](const size_type i) const;
625 
631  reference operator[](const size_type i);
632 
648  template <typename OtherNumber>
649  void
650  extract_subvector_to(const std::vector<size_type> &indices,
651  std::vector<OtherNumber> & values) const;
652 
680  template <typename ForwardIterator, typename OutputIterator>
681  void
682  extract_subvector_to(ForwardIterator indices_begin,
683  const ForwardIterator indices_end,
684  OutputIterator values_begin) const;
685 
691  operator=(const value_type s);
692 
697  operator=(const BlockVectorBase &V);
698 
705  operator=(BlockVectorBase && /*V*/) = default; // NOLINT
706 
710  template <class VectorType2>
712  operator=(const BlockVectorBase<VectorType2> &V);
713 
718  operator=(const VectorType &v);
719 
724  template <class VectorType2>
725  bool
726  operator==(const BlockVectorBase<VectorType2> &v) const;
727 
731  value_type operator*(const BlockVectorBase &V) const;
732 
736  real_type
737  norm_sqr() const;
738 
742  value_type
743  mean_value() const;
744 
748  real_type
749  l1_norm() const;
750 
755  real_type
756  l2_norm() const;
757 
762  real_type
763  linfty_norm() const;
764 
787  value_type
788  add_and_dot(const value_type a,
789  const BlockVectorBase &V,
790  const BlockVectorBase &W);
791 
796  bool
797  in_local_range(const size_type global_index) const;
798 
804  bool
805  all_zero() const;
806 
812  bool
813  is_non_negative() const;
814 
819  operator+=(const BlockVectorBase &V);
820 
825  operator-=(const BlockVectorBase &V);
826 
827 
832  template <typename Number>
833  void
834  add(const std::vector<size_type> &indices, const std::vector<Number> &values);
835 
840  template <typename Number>
841  void
842  add(const std::vector<size_type> &indices, const Vector<Number> &values);
843 
849  template <typename Number>
850  void
851  add(const size_type n_elements,
852  const size_type *indices,
853  const Number * values);
854 
859  void
860  add(const value_type s);
861 
865  void
866  add(const value_type a, const BlockVectorBase &V);
867 
871  void
872  add(const value_type a,
873  const BlockVectorBase &V,
874  const value_type b,
875  const BlockVectorBase &W);
876 
880  void
881  sadd(const value_type s, const BlockVectorBase &V);
882 
886  void
887  sadd(const value_type s, const value_type a, const BlockVectorBase &V);
888 
892  void
893  sadd(const value_type s,
894  const value_type a,
895  const BlockVectorBase &V,
896  const value_type b,
897  const BlockVectorBase &W);
898 
902  void
903  sadd(const value_type s,
904  const value_type a,
905  const BlockVectorBase &V,
906  const value_type b,
907  const BlockVectorBase &W,
908  const value_type c,
909  const BlockVectorBase &X);
910 
915  operator*=(const value_type factor);
916 
921  operator/=(const value_type factor);
922 
927  template <class BlockVector2>
928  void
929  scale(const BlockVector2 &v);
930 
934  template <class BlockVector2>
935  void
936  equ(const value_type a, const BlockVector2 &V);
937 
942  void
943  update_ghost_values() const;
944 
949  std::size_t
950  memory_consumption() const;
951 
952 protected:
956  std::vector<VectorType> components;
957 
962  BlockIndices block_indices;
963 
964  // Make the iterator class a friend.
965  template <typename N, bool C>
966  friend class ::internal::BlockVectorIterators::Iterator;
967 
968  template <typename>
969  friend class BlockVectorBase;
970 };
971 
972 
975 /*----------------------- Inline functions ----------------------------------*/
976 
977 
978 #ifndef DOXYGEN
979 namespace internal
980 {
981  namespace BlockVectorIterators
982  {
983  template <class BlockVectorType, bool Constness>
984  inline Iterator<BlockVectorType, Constness>::Iterator(
985  const Iterator<BlockVectorType, Constness> &c)
986  : parent(c.parent)
987  , global_index(c.global_index)
988  , current_block(c.current_block)
989  , index_within_block(c.index_within_block)
990  , next_break_forward(c.next_break_forward)
991  , next_break_backward(c.next_break_backward)
992  {}
993 
994 
995 
996  template <class BlockVectorType, bool Constness>
997  inline Iterator<BlockVectorType, Constness>::Iterator(
998  const Iterator<BlockVectorType, !Constness> &c)
999  : parent(c.parent)
1000  , global_index(c.global_index)
1001  , current_block(c.current_block)
1002  , index_within_block(c.index_within_block)
1003  , next_break_forward(c.next_break_forward)
1004  , next_break_backward(c.next_break_backward)
1005  {
1006  // Only permit copy-constructing const iterators from non-const
1007  // iterators, and not vice versa (i.e., Constness must always be
1008  // true).
1009  static_assert(Constness == true,
1010  "Constructing a non-const iterator from a const iterator "
1011  "does not make sense.");
1012  }
1013 
1014 
1015 
1016  template <class BlockVectorType, bool Constness>
1017  inline Iterator<BlockVectorType, Constness>::Iterator(
1018  BlockVector & parent,
1019  const size_type global_index,
1020  const size_type current_block,
1021  const size_type index_within_block,
1022  const size_type next_break_forward,
1023  const size_type next_break_backward)
1024  : parent(&parent)
1025  , global_index(global_index)
1026  , current_block(current_block)
1027  , index_within_block(index_within_block)
1028  , next_break_forward(next_break_forward)
1029  , next_break_backward(next_break_backward)
1030  {}
1031 
1032 
1033 
1034  template <class BlockVectorType, bool Constness>
1035  inline Iterator<BlockVectorType, Constness> &
1036  Iterator<BlockVectorType, Constness>::operator=(const Iterator &c)
1037  {
1038  parent = c.parent;
1039  global_index = c.global_index;
1040  index_within_block = c.index_within_block;
1041  current_block = c.current_block;
1042  next_break_forward = c.next_break_forward;
1043  next_break_backward = c.next_break_backward;
1044 
1045  return *this;
1046  }
1047 
1048 
1049 
1050  template <class BlockVectorType, bool Constness>
1051  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1053  {
1054  return parent->block(current_block)(index_within_block);
1055  }
1056 
1057 
1058 
1059  template <class BlockVectorType, bool Constness>
1060  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1061  Iterator<BlockVectorType, Constness>::
1062  operator[](const difference_type d) const
1063  {
1064  // if the index pointed to is
1065  // still within the block we
1066  // currently point into, then we
1067  // can save the computation of
1068  // the block
1069  if ((global_index + d >= next_break_backward) &&
1070  (global_index + d <= next_break_forward))
1071  return parent->block(current_block)(index_within_block + d);
1072 
1073  // if the index is not within the
1074  // block of the block vector into
1075  // which we presently point, then
1076  // there is no way: we have to
1077  // search for the block. this can
1078  // be done through the parent
1079  // class as well.
1080  return (*parent)(global_index + d);
1081  }
1082 
1083 
1084 
1085  template <class BlockVectorType, bool Constness>
1086  inline Iterator<BlockVectorType, Constness> &
1087  Iterator<BlockVectorType, Constness>::operator++()
1088  {
1089  move_forward();
1090  return *this;
1091  }
1092 
1093 
1094 
1095  template <class BlockVectorType, bool Constness>
1096  inline Iterator<BlockVectorType, Constness>
1097  Iterator<BlockVectorType, Constness>::operator++(int)
1098  {
1099  const Iterator old_value = *this;
1100  move_forward();
1101  return old_value;
1102  }
1103 
1104 
1105 
1106  template <class BlockVectorType, bool Constness>
1107  inline Iterator<BlockVectorType, Constness> &
1108  Iterator<BlockVectorType, Constness>::operator--()
1109  {
1110  move_backward();
1111  return *this;
1112  }
1113 
1114 
1115 
1116  template <class BlockVectorType, bool Constness>
1117  inline Iterator<BlockVectorType, Constness>
1118  Iterator<BlockVectorType, Constness>::operator--(int)
1119  {
1120  const Iterator old_value = *this;
1121  move_backward();
1122  return old_value;
1123  }
1124 
1125 
1126 
1127  template <class BlockVectorType, bool Constness>
1128  template <bool OtherConstness>
1129  inline bool
1130  Iterator<BlockVectorType, Constness>::
1131  operator==(const Iterator<BlockVectorType, OtherConstness> &i) const
1132  {
1133  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1134 
1135  return (global_index == i.global_index);
1136  }
1137 
1138 
1139 
1140  template <class BlockVectorType, bool Constness>
1141  template <bool OtherConstness>
1142  inline bool
1143  Iterator<BlockVectorType, Constness>::
1144  operator!=(const Iterator<BlockVectorType, OtherConstness> &i) const
1145  {
1146  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1147 
1148  return (global_index != i.global_index);
1149  }
1150 
1151 
1152 
1153  template <class BlockVectorType, bool Constness>
1154  template <bool OtherConstness>
1155  inline bool
1156  Iterator<BlockVectorType, Constness>::
1157  operator<(const Iterator<BlockVectorType, OtherConstness> &i) const
1158  {
1159  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1160 
1161  return (global_index < i.global_index);
1162  }
1163 
1164 
1165 
1166  template <class BlockVectorType, bool Constness>
1167  template <bool OtherConstness>
1168  inline bool
1169  Iterator<BlockVectorType, Constness>::
1170  operator<=(const Iterator<BlockVectorType, OtherConstness> &i) const
1171  {
1172  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1173 
1174  return (global_index <= i.global_index);
1175  }
1176 
1177 
1178 
1179  template <class BlockVectorType, bool Constness>
1180  template <bool OtherConstness>
1181  inline bool
1182  Iterator<BlockVectorType, Constness>::
1183  operator>(const Iterator<BlockVectorType, OtherConstness> &i) const
1184  {
1185  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1186 
1187  return (global_index > i.global_index);
1188  }
1189 
1190 
1191 
1192  template <class BlockVectorType, bool Constness>
1193  template <bool OtherConstness>
1194  inline bool
1195  Iterator<BlockVectorType, Constness>::
1196  operator>=(const Iterator<BlockVectorType, OtherConstness> &i) const
1197  {
1198  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1199 
1200  return (global_index >= i.global_index);
1201  }
1202 
1203 
1204 
1205  template <class BlockVectorType, bool Constness>
1206  template <bool OtherConstness>
1207  inline typename Iterator<BlockVectorType, Constness>::difference_type
1209  operator-(const Iterator<BlockVectorType, OtherConstness> &i) const
1210  {
1211  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1212 
1213  return (static_cast<signed int>(global_index) -
1214  static_cast<signed int>(i.global_index));
1215  }
1216 
1217 
1218 
1219  template <class BlockVectorType, bool Constness>
1220  inline Iterator<BlockVectorType, Constness>
1222  operator+(const difference_type &d) const
1223  {
1224  // if the index pointed to is
1225  // still within the block we
1226  // currently point into, then we
1227  // can save the computation of
1228  // the block
1229  if ((global_index + d >= next_break_backward) &&
1230  (global_index + d <= next_break_forward))
1231  return Iterator(*parent,
1232  global_index + d,
1233  current_block,
1234  index_within_block + d,
1235  next_break_forward,
1236  next_break_backward);
1237  else
1238  // outside present block, so
1239  // have to seek new block
1240  // anyway
1241  return Iterator(*parent, global_index + d);
1242  }
1243 
1244 
1245 
1246  template <class BlockVectorType, bool Constness>
1247  inline Iterator<BlockVectorType, Constness>
1249  operator-(const difference_type &d) const
1250  {
1251  // if the index pointed to is
1252  // still within the block we
1253  // currently point into, then we
1254  // can save the computation of
1255  // the block
1256  if ((global_index - d >= next_break_backward) &&
1257  (global_index - d <= next_break_forward))
1258  return Iterator(*parent,
1259  global_index - d,
1260  current_block,
1261  index_within_block - d,
1262  next_break_forward,
1263  next_break_backward);
1264  else
1265  // outside present block, so
1266  // have to seek new block
1267  // anyway
1268  return Iterator(*parent, global_index - d);
1269  }
1270 
1271 
1272 
1273  template <class BlockVectorType, bool Constness>
1274  inline Iterator<BlockVectorType, Constness> &
1275  Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1276  {
1277  // if the index pointed to is
1278  // still within the block we
1279  // currently point into, then we
1280  // can save the computation of
1281  // the block
1282  if ((global_index + d >= next_break_backward) &&
1283  (global_index + d <= next_break_forward))
1284  {
1285  global_index += d;
1286  index_within_block += d;
1287  }
1288  else
1289  // outside present block, so
1290  // have to seek new block
1291  // anyway
1292  *this = Iterator(*parent, global_index + d);
1293 
1294  return *this;
1295  }
1296 
1297 
1298 
1299  template <class BlockVectorType, bool Constness>
1300  inline Iterator<BlockVectorType, Constness> &
1301  Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1302  {
1303  // if the index pointed to is
1304  // still within the block we
1305  // currently point into, then we
1306  // can save the computation of
1307  // the block
1308  if ((global_index - d >= next_break_backward) &&
1309  (global_index - d <= next_break_forward))
1310  {
1311  global_index -= d;
1312  index_within_block -= d;
1313  }
1314  else
1315  // outside present block, so
1316  // have to seek new block
1317  // anyway
1318  *this = Iterator(*parent, global_index - d);
1319 
1320  return *this;
1321  }
1322 
1323 
1324  template <class BlockVectorType, bool Constness>
1325  Iterator<BlockVectorType, Constness>::Iterator(BlockVector & parent,
1326  const size_type global_index)
1327  : parent(&parent)
1328  , global_index(global_index)
1329  {
1330  // find which block we are
1331  // in. for this, take into
1332  // account that it happens at
1333  // times that people want to
1334  // initialize iterators
1335  // past-the-end
1336  if (global_index < parent.size())
1337  {
1338  const std::pair<size_type, size_type> indices =
1339  parent.block_indices.global_to_local(global_index);
1340  current_block = indices.first;
1341  index_within_block = indices.second;
1342 
1343  next_break_backward =
1344  parent.block_indices.local_to_global(current_block, 0);
1345  next_break_forward =
1346  (parent.block_indices.local_to_global(current_block, 0) +
1347  parent.block_indices.block_size(current_block) - 1);
1348  }
1349  else
1350  // past the end. only have one
1351  // value for this
1352  {
1353  this->global_index = parent.size();
1354  current_block = parent.n_blocks();
1355  index_within_block = 0;
1356  next_break_backward = global_index;
1357  next_break_forward = numbers::invalid_size_type;
1358  };
1359  }
1360 
1361 
1362 
1363  template <class BlockVectorType, bool Constness>
1364  void
1365  Iterator<BlockVectorType, Constness>::move_forward()
1366  {
1367  if (global_index != next_break_forward)
1368  ++index_within_block;
1369  else
1370  {
1371  // ok, we traverse a boundary
1372  // between blocks:
1373  index_within_block = 0;
1374  ++current_block;
1375 
1376  // break backwards is now old
1377  // break forward
1378  next_break_backward = next_break_forward + 1;
1379 
1380  // compute new break forward
1381  if (current_block < parent->block_indices.size())
1382  next_break_forward +=
1383  parent->block_indices.block_size(current_block);
1384  else
1385  // if we are beyond the end,
1386  // then move the next
1387  // boundary arbitrarily far
1388  // away
1389  next_break_forward = numbers::invalid_size_type;
1390  };
1391 
1392  ++global_index;
1393  }
1394 
1395 
1396 
1397  template <class BlockVectorType, bool Constness>
1398  void
1399  Iterator<BlockVectorType, Constness>::move_backward()
1400  {
1401  if (global_index != next_break_backward)
1402  --index_within_block;
1403  else if (current_block != 0)
1404  {
1405  // ok, we traverse a boundary
1406  // between blocks:
1407  --current_block;
1408  index_within_block =
1409  parent->block_indices.block_size(current_block) - 1;
1410 
1411  // break forwards is now old
1412  // break backward
1413  next_break_forward = next_break_backward - 1;
1414 
1415  // compute new break forward
1416  next_break_backward -=
1417  parent->block_indices.block_size(current_block);
1418  }
1419  else
1420  // current block was 0, we now
1421  // get into unspecified terrain
1422  {
1423  --current_block;
1424  index_within_block = numbers::invalid_size_type;
1425  next_break_forward = 0;
1426  next_break_backward = 0;
1427  };
1428 
1429  --global_index;
1430  }
1431 
1432 
1433  } // namespace BlockVectorIterators
1434 
1435 } // namespace internal
1436 
1437 
1438 
1439 template <class VectorType>
1440 inline std::size_t
1442 {
1443  return block_indices.total_size();
1444 }
1445 
1446 
1447 
1448 template <class VectorType>
1449 inline IndexSet
1451 {
1452  IndexSet is(size());
1453 
1454  // copy index sets from blocks into the global one, shifted
1455  // by the appropriate amount for each block
1456  for (unsigned int b = 0; b < n_blocks(); ++b)
1457  {
1458  IndexSet x = block(b).locally_owned_elements();
1459  is.add_indices(x, block_indices.block_start(b));
1460  }
1461 
1462  is.compress();
1463 
1464  return is;
1465 }
1466 
1467 
1468 
1469 template <class VectorType>
1470 inline unsigned int
1472 {
1473  return block_indices.size();
1474 }
1475 
1476 
1477 template <class VectorType>
1479 BlockVectorBase<VectorType>::block(const unsigned int i)
1480 {
1481  Assert(i < n_blocks(), ExcIndexRange(i, 0, n_blocks()));
1482 
1483  return components[i];
1484 }
1485 
1486 
1487 
1488 template <class VectorType>
1489 inline const typename BlockVectorBase<VectorType>::BlockType &
1490 BlockVectorBase<VectorType>::block(const unsigned int i) const
1491 {
1492  Assert(i < n_blocks(), ExcIndexRange(i, 0, n_blocks()));
1493 
1494  return components[i];
1495 }
1496 
1497 
1498 
1499 template <class VectorType>
1500 inline const BlockIndices &
1502 {
1503  return block_indices;
1504 }
1505 
1506 
1507 template <class VectorType>
1508 inline void
1510 {
1511  std::vector<size_type> sizes(n_blocks());
1512 
1513  for (size_type i = 0; i < n_blocks(); ++i)
1514  sizes[i] = block(i).size();
1515 
1516  block_indices.reinit(sizes);
1517 }
1518 
1519 
1520 
1521 template <class VectorType>
1522 inline void
1524  ::VectorOperation::values operation)
1525 {
1526  for (unsigned int i = 0; i < n_blocks(); ++i)
1527  block(i).compress(operation);
1528 }
1529 
1530 
1531 
1532 template <class VectorType>
1535 {
1536  return iterator(*this, 0U);
1537 }
1538 
1539 
1540 
1541 template <class VectorType>
1544 {
1545  return const_iterator(*this, 0U);
1546 }
1547 
1548 
1549 template <class VectorType>
1552 {
1553  return iterator(*this, size());
1554 }
1555 
1556 
1557 
1558 template <class VectorType>
1561 {
1562  return const_iterator(*this, size());
1563 }
1564 
1565 
1566 template <class VectorType>
1567 inline bool
1568 BlockVectorBase<VectorType>::in_local_range(const size_type global_index) const
1569 {
1570  const std::pair<size_type, size_type> local_index =
1571  block_indices.global_to_local(global_index);
1572 
1573  return components[local_index.first].in_local_range(global_index);
1574 }
1575 
1576 
1577 template <class VectorType>
1578 bool
1580 {
1581  for (size_type i = 0; i < n_blocks(); ++i)
1582  if (components[i].all_zero() == false)
1583  return false;
1584 
1585  return true;
1586 }
1587 
1588 
1589 
1590 template <class VectorType>
1591 bool
1593 {
1594  for (size_type i = 0; i < n_blocks(); ++i)
1595  if (components[i].is_non_negative() == false)
1596  return false;
1597 
1598  return true;
1599 }
1600 
1601 
1602 
1603 template <class VectorType>
1604 typename BlockVectorBase<VectorType>::value_type BlockVectorBase<VectorType>::
1606 {
1607  Assert(n_blocks() == v.n_blocks(),
1608  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1609 
1610  value_type sum = 0.;
1611  for (size_type i = 0; i < n_blocks(); ++i)
1612  sum += components[i] * v.components[i];
1613 
1614  return sum;
1615 }
1616 
1617 
1618 template <class VectorType>
1621 {
1622  real_type sum = 0.;
1623  for (size_type i = 0; i < n_blocks(); ++i)
1624  sum += components[i].norm_sqr();
1625 
1626  return sum;
1627 }
1628 
1629 
1630 
1631 template <class VectorType>
1632 typename BlockVectorBase<VectorType>::value_type
1634 {
1635  value_type sum = 0.;
1636  // need to do static_cast as otherwise it won't work with
1637  // value_type=complex<T>
1638  for (size_type i = 0; i < n_blocks(); ++i)
1639  sum += components[i].mean_value() *
1641  components[i].size()));
1642 
1643  return sum / (typename numbers::NumberTraits<value_type>::real_type(size()));
1644 }
1645 
1646 
1647 
1648 template <class VectorType>
1651 {
1652  real_type sum = 0.;
1653  for (size_type i = 0; i < n_blocks(); ++i)
1654  sum += components[i].l1_norm();
1655 
1656  return sum;
1657 }
1658 
1659 
1660 
1661 template <class VectorType>
1664 {
1665  return std::sqrt(norm_sqr());
1666 }
1667 
1668 
1669 
1670 template <class VectorType>
1673 {
1674  real_type sum = 0.;
1675  for (size_type i = 0; i < n_blocks(); ++i)
1676  {
1677  value_type newval = components[i].linfty_norm();
1678  if (sum < newval)
1679  sum = newval;
1680  }
1681  return sum;
1682 }
1683 
1684 
1685 
1686 template <class VectorType>
1687 typename BlockVectorBase<VectorType>::value_type
1689  const typename BlockVectorBase<VectorType>::value_type a,
1690  const BlockVectorBase<VectorType> & V,
1691  const BlockVectorBase<VectorType> & W)
1692 {
1693  AssertDimension(n_blocks(), V.n_blocks());
1694  AssertDimension(n_blocks(), W.n_blocks());
1695 
1696  value_type sum = 0.;
1697  for (size_type i = 0; i < n_blocks(); ++i)
1698  sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1699 
1700  return sum;
1701 }
1702 
1703 
1704 
1705 template <class VectorType>
1708 {
1709  Assert(n_blocks() == v.n_blocks(),
1710  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1711 
1712  for (size_type i = 0; i < n_blocks(); ++i)
1713  {
1714  components[i] += v.components[i];
1715  }
1716 
1717  return *this;
1718 }
1719 
1720 
1721 
1722 template <class VectorType>
1725 {
1726  Assert(n_blocks() == v.n_blocks(),
1727  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1728 
1729  for (size_type i = 0; i < n_blocks(); ++i)
1730  {
1731  components[i] -= v.components[i];
1732  }
1733  return *this;
1734 }
1735 
1736 
1737 
1738 template <class VectorType>
1739 template <typename Number>
1740 inline void
1741 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1742  const std::vector<Number> & values)
1743 {
1744  Assert(indices.size() == values.size(),
1745  ExcDimensionMismatch(indices.size(), values.size()));
1746  add(indices.size(), indices.data(), values.data());
1747 }
1748 
1749 
1750 
1751 template <class VectorType>
1752 template <typename Number>
1753 inline void
1754 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1755  const Vector<Number> & values)
1756 {
1757  Assert(indices.size() == values.size(),
1758  ExcDimensionMismatch(indices.size(), values.size()));
1759  const size_type n_indices = indices.size();
1760  for (size_type i = 0; i < n_indices; ++i)
1761  (*this)(indices[i]) += values(i);
1762 }
1763 
1764 
1765 
1766 template <class VectorType>
1767 template <typename Number>
1768 inline void
1769 BlockVectorBase<VectorType>::add(const size_type n_indices,
1770  const size_type *indices,
1771  const Number * values)
1772 {
1773  for (size_type i = 0; i < n_indices; ++i)
1774  (*this)(indices[i]) += values[i];
1775 }
1776 
1777 
1778 
1779 template <class VectorType>
1780 void
1781 BlockVectorBase<VectorType>::add(const value_type a)
1782 {
1783  AssertIsFinite(a);
1784 
1785  for (size_type i = 0; i < n_blocks(); ++i)
1786  {
1787  components[i].add(a);
1788  }
1789 }
1790 
1791 
1792 
1793 template <class VectorType>
1794 void
1795 BlockVectorBase<VectorType>::add(const value_type a,
1796  const BlockVectorBase<VectorType> &v)
1797 {
1798  AssertIsFinite(a);
1799 
1800  Assert(n_blocks() == v.n_blocks(),
1801  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1802 
1803  for (size_type i = 0; i < n_blocks(); ++i)
1804  {
1805  components[i].add(a, v.components[i]);
1806  }
1807 }
1808 
1809 
1810 
1811 template <class VectorType>
1812 void
1813 BlockVectorBase<VectorType>::add(const value_type a,
1814  const BlockVectorBase<VectorType> &v,
1815  const value_type b,
1816  const BlockVectorBase<VectorType> &w)
1817 {
1818  AssertIsFinite(a);
1819  AssertIsFinite(b);
1820 
1821  Assert(n_blocks() == v.n_blocks(),
1822  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1823  Assert(n_blocks() == w.n_blocks(),
1824  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1825 
1826 
1827  for (size_type i = 0; i < n_blocks(); ++i)
1828  {
1829  components[i].add(a, v.components[i], b, w.components[i]);
1830  }
1831 }
1832 
1833 
1834 
1835 template <class VectorType>
1836 void
1837 BlockVectorBase<VectorType>::sadd(const value_type x,
1838  const BlockVectorBase<VectorType> &v)
1839 {
1840  AssertIsFinite(x);
1841 
1842  Assert(n_blocks() == v.n_blocks(),
1843  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1844 
1845  for (size_type i = 0; i < n_blocks(); ++i)
1846  {
1847  components[i].sadd(x, v.components[i]);
1848  }
1849 }
1850 
1851 
1852 
1853 template <class VectorType>
1854 void
1855 BlockVectorBase<VectorType>::sadd(const value_type x,
1856  const value_type a,
1857  const BlockVectorBase<VectorType> &v)
1858 {
1859  AssertIsFinite(x);
1860  AssertIsFinite(a);
1861 
1862  Assert(n_blocks() == v.n_blocks(),
1863  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1864 
1865  for (size_type i = 0; i < n_blocks(); ++i)
1866  {
1867  components[i].sadd(x, a, v.components[i]);
1868  }
1869 }
1870 
1871 
1872 
1873 template <class VectorType>
1874 void
1875 BlockVectorBase<VectorType>::sadd(const value_type x,
1876  const value_type a,
1877  const BlockVectorBase<VectorType> &v,
1878  const value_type b,
1879  const BlockVectorBase<VectorType> &w)
1880 {
1881  AssertIsFinite(x);
1882  AssertIsFinite(a);
1883  AssertIsFinite(b);
1884 
1885  Assert(n_blocks() == v.n_blocks(),
1886  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1887  Assert(n_blocks() == w.n_blocks(),
1888  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1889 
1890  for (size_type i = 0; i < n_blocks(); ++i)
1891  {
1892  components[i].sadd(x, a, v.components[i], b, w.components[i]);
1893  }
1894 }
1895 
1896 
1897 
1898 template <class VectorType>
1899 void
1900 BlockVectorBase<VectorType>::sadd(const value_type x,
1901  const value_type a,
1902  const BlockVectorBase<VectorType> &v,
1903  const value_type b,
1904  const BlockVectorBase<VectorType> &w,
1905  const value_type c,
1906  const BlockVectorBase<VectorType> &y)
1907 {
1908  AssertIsFinite(x);
1909  AssertIsFinite(a);
1910  AssertIsFinite(b);
1911  AssertIsFinite(c);
1912 
1913  Assert(n_blocks() == v.n_blocks(),
1914  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1915  Assert(n_blocks() == w.n_blocks(),
1916  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1917  Assert(n_blocks() == y.n_blocks(),
1918  ExcDimensionMismatch(n_blocks(), y.n_blocks()));
1919 
1920  for (size_type i = 0; i < n_blocks(); ++i)
1921  {
1922  components[i].sadd(
1923  x, a, v.components[i], b, w.components[i], c, y.components[i]);
1924  }
1925 }
1926 
1927 
1928 
1929 template <class VectorType>
1930 template <class BlockVector2>
1931 void
1932 BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
1933 {
1934  Assert(n_blocks() == v.n_blocks(),
1935  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1936  for (size_type i = 0; i < n_blocks(); ++i)
1937  components[i].scale(v.block(i));
1938 }
1939 
1940 
1941 
1942 template <class VectorType>
1943 std::size_t
1945 {
1946  return (MemoryConsumption::memory_consumption(this->block_indices) +
1947  MemoryConsumption::memory_consumption(this->components));
1948 }
1949 
1950 
1951 
1952 template <class VectorType>
1953 template <class BlockVector2>
1954 void
1955 BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
1956 {
1957  AssertIsFinite(a);
1958 
1959  Assert(n_blocks() == v.n_blocks(),
1960  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1961 
1962  for (size_type i = 0; i < n_blocks(); ++i)
1963  components[i].equ(a, v.components[i]);
1964 }
1965 
1966 
1967 
1968 template <class VectorType>
1969 void
1971 {
1972  for (size_type i = 0; i < n_blocks(); ++i)
1973  block(i).update_ghost_values();
1974 }
1975 
1976 
1977 
1978 template <class VectorType>
1980 BlockVectorBase<VectorType>::operator=(const value_type s)
1981 {
1982  AssertIsFinite(s);
1983 
1984  for (size_type i = 0; i < n_blocks(); ++i)
1985  components[i] = s;
1986 
1987  return *this;
1988 }
1989 
1990 
1991 template <class VectorType>
1994 {
1995  AssertDimension(n_blocks(), v.n_blocks());
1996 
1997  for (size_type i = 0; i < n_blocks(); ++i)
1998  components[i] = v.components[i];
1999 
2000  return *this;
2001 }
2002 
2003 
2004 template <class VectorType>
2005 template <class VectorType2>
2008 {
2009  AssertDimension(n_blocks(), v.n_blocks());
2010 
2011  for (size_type i = 0; i < n_blocks(); ++i)
2012  components[i] = v.components[i];
2013 
2014  return *this;
2015 }
2016 
2017 
2018 
2019 template <class VectorType>
2021 BlockVectorBase<VectorType>::operator=(const VectorType &v)
2022 {
2023  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2024 
2025  size_type index_v = 0;
2026  for (size_type b = 0; b < n_blocks(); ++b)
2027  for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2028  block(b)(i) = v(index_v);
2029 
2030  return *this;
2031 }
2032 
2033 
2034 
2035 template <class VectorType>
2036 template <class VectorType2>
2037 inline bool
2040 {
2041  Assert(block_indices == v.block_indices, ExcDifferentBlockIndices());
2042 
2043  for (size_type i = 0; i < n_blocks(); ++i)
2044  if (!(components[i] == v.components[i]))
2045  return false;
2046 
2047  return true;
2048 }
2049 
2050 
2051 
2052 template <class VectorType>
2054 BlockVectorBase<VectorType>::operator*=(const value_type factor)
2055 {
2056  AssertIsFinite(factor);
2057 
2058  for (size_type i = 0; i < n_blocks(); ++i)
2059  components[i] *= factor;
2060 
2061  return *this;
2062 }
2063 
2064 
2065 
2066 template <class VectorType>
2068 BlockVectorBase<VectorType>::operator/=(const value_type factor)
2069 {
2070  AssertIsFinite(factor);
2071  Assert(factor != 0., ExcDivideByZero());
2072 
2073  const value_type inverse_factor = value_type(1.) / factor;
2074 
2075  for (size_type i = 0; i < n_blocks(); ++i)
2076  components[i] *= inverse_factor;
2077 
2078  return *this;
2079 }
2080 
2081 
2082 template <class VectorType>
2083 inline typename BlockVectorBase<VectorType>::value_type
2084 BlockVectorBase<VectorType>::operator()(const size_type i) const
2085 {
2086  const std::pair<unsigned int, size_type> local_index =
2087  block_indices.global_to_local(i);
2088  return components[local_index.first](local_index.second);
2089 }
2090 
2091 
2092 
2093 template <class VectorType>
2094 inline typename BlockVectorBase<VectorType>::reference
2095 BlockVectorBase<VectorType>::operator()(const size_type i)
2096 {
2097  const std::pair<unsigned int, size_type> local_index =
2098  block_indices.global_to_local(i);
2099  return components[local_index.first](local_index.second);
2100 }
2101 
2102 
2103 
2104 template <class VectorType>
2105 inline typename BlockVectorBase<VectorType>::value_type
2106  BlockVectorBase<VectorType>::operator[](const size_type i) const
2107 {
2108  return operator()(i);
2109 }
2110 
2111 
2112 
2113 template <class VectorType>
2114 inline typename BlockVectorBase<VectorType>::reference
2115  BlockVectorBase<VectorType>::operator[](const size_type i)
2116 {
2117  return operator()(i);
2118 }
2119 
2120 
2121 
2122 template <typename VectorType>
2123 template <typename OtherNumber>
2124 inline void
2126  const std::vector<size_type> &indices,
2127  std::vector<OtherNumber> & values) const
2128 {
2129  for (size_type i = 0; i < indices.size(); ++i)
2130  values[i] = operator()(indices[i]);
2131 }
2132 
2133 
2134 
2135 template <typename VectorType>
2136 template <typename ForwardIterator, typename OutputIterator>
2137 inline void
2139  ForwardIterator indices_begin,
2140  const ForwardIterator indices_end,
2141  OutputIterator values_begin) const
2142 {
2143  while (indices_begin != indices_end)
2144  {
2145  *values_begin = operator()(*indices_begin);
2146  indices_begin++;
2147  values_begin++;
2148  }
2149 }
2150 
2151 #endif // DOXYGEN
2152 
2153 DEAL_II_NAMESPACE_CLOSE
2154 
2155 #endif
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
const types::global_dof_index invalid_size_type
Definition: types.h:196
static yes_type check_for_block_vector(const BlockVectorBase< T > *)
#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
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
bool operator==(const BlockVectorBase< VectorType2 > &v) const
void compress(::VectorOperation::values operation)
typename numbers::NumberTraits< number >::real_type real_type
Definition: vector.h:143
bool in_local_range(const size_type global_index) const
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1670
__global__ void add_and_dot(Number *res, Number *v1, const Number *v2, const Number *v3, const Number a, const size_type N)
void collect_sizes()
STL namespace.
size_type block_size(const unsigned int i) const
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
static ::ExceptionBase & ExcDivideByZero()
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
__global__ void scale(Number *val, const Number *V_val, const size_type N)
void scale(const BlockVector2 &v)
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
number value_type
Definition: vector.h:125
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
static ::ExceptionBase & ExcDifferentBlockIndices()
static const bool value
__global__ void equ(Number *val, const Number a, const Number *V_val, const size_type N)
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)
bool all_zero() const
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
void compress() const
Definition: index_set.h:1608
unsigned int n_blocks() const
BlockVectorBase & operator-=(const BlockVectorBase &V)
typename BlockType::real_type real_type
iterator end()
size_type block_start(const unsigned int i) const
value_type operator()(const size_type i) const
real_type linfty_norm() const
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:1669
void equ(const value_type a, const BlockVector2 &V)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
value_type mean_value() const