Reference documentation for deal.II version Git e7bb9ce7b3 2020-09-18 12:07:32 -0400
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
block_vector_base.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_block_vector_base_h
17 #define dealii_block_vector_base_h
18 
19 
20 #include <deal.II/base/config.h>
21 
24 #include <deal.II/base/numbers.h>
26 
28 #include <deal.II/lac/exceptions.h>
29 #include <deal.II/lac/vector.h>
31 
32 #include <cmath>
33 #include <cstddef>
34 #include <iterator>
35 #include <vector>
36 
38 
39 
44 // Forward declaration
45 #ifndef DOXYGEN
46 template <typename>
47 class BlockVectorBase;
48 #endif
49 
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 {
115  namespace BlockVectorIterators
116  {
132  template <class BlockVectorType, bool Constness>
133  class Iterator
134  {
135  public:
140 
146  using value_type =
147  typename std::conditional<Constness,
148  const typename BlockVectorType::value_type,
149  typename BlockVectorType::value_type>::type;
150 
156  using iterator_category = std::random_access_iterator_tag;
157  using difference_type = std::ptrdiff_t;
158  using reference = typename BlockVectorType::reference;
159  using pointer = value_type *;
160 
161  using dereference_type = typename std::conditional<
162  Constness,
163  value_type,
164  typename BlockVectorType::BlockType::reference>::type;
165 
170  using BlockVector = typename std::
171  conditional<Constness, const BlockVectorType, BlockVectorType>::type;
172 
181  Iterator(BlockVector &parent, const size_type global_index);
182 
191 
192 
196  Iterator(const Iterator &c);
197 
198  private:
203  Iterator(BlockVector & parent,
204  const size_type global_index,
205  const size_type current_block,
206  const size_type index_within_block,
207  const size_type next_break_forward,
208  const size_type next_break_backward);
209 
210  public:
214  Iterator &
215  operator=(const Iterator &c);
216 
222  dereference_type operator*() const;
223 
228  dereference_type operator[](const difference_type d) const;
229 
234  Iterator &
235  operator++();
236 
242  Iterator
243  operator++(int);
244 
249  Iterator &
250  operator--();
251 
257  Iterator
258  operator--(int);
259 
264  template <bool OtherConstness>
265  bool
267 
272  template <bool OtherConstness>
273  bool
275 
281  template <bool OtherConstness>
282  bool
283  operator<(const Iterator<BlockVectorType, OtherConstness> &i) const;
284 
288  template <bool OtherConstness>
289  bool
290  operator<=(const Iterator<BlockVectorType, OtherConstness> &i) const;
291 
295  template <bool OtherConstness>
296  bool
298 
302  template <bool OtherConstness>
303  bool
305 
309  template <bool OtherConstness>
312 
317  Iterator
318  operator+(const difference_type &d) const;
319 
324  Iterator
325  operator-(const difference_type &d) const;
326 
331  Iterator &
332  operator+=(const difference_type &d);
333 
338  Iterator &
339  operator-=(const difference_type &d);
340 
350  DeclExceptionMsg(ExcPointerToDifferentVectors,
351  "Your program tried to compare iterators pointing to "
352  "different block vectors. There is no reasonable way "
353  "to do this.");
354 
356  private:
363 
368 
373  unsigned int current_block;
375 
384 
388  void
389  move_forward();
390 
394  void
395  move_backward();
396 
397  // Mark all other instances of this template as friends.
398  template <typename, bool>
399  friend class Iterator;
400  };
401  } // namespace BlockVectorIterators
402 } // namespace internal
403 
404 
442 template <class VectorType>
444 {
445 public:
450 
451  /*
452  * Declare standard types used in
453  * all containers. These types
454  * parallel those in the
455  * <tt>C++</tt> standard
456  * libraries
457  * <tt>std::vector<...></tt>
458  * class. This includes iterator
459  * types.
460  */
462  using pointer = value_type *;
463  using const_pointer = const value_type *;
464  using iterator =
466  using const_iterator =
468  using reference = typename BlockType::reference;
471 
481  using real_type = typename BlockType::real_type;
482 
486  BlockVectorBase() = default;
487 
491  BlockVectorBase(const BlockVectorBase & /*V*/) = default;
492 
498  BlockVectorBase(BlockVectorBase && /*V*/) noexcept = default;
499 
506  void
507  collect_sizes();
508 
519  void
520  compress(::VectorOperation::values operation);
521 
525  BlockType &
526  block(const unsigned int i);
527 
531  const BlockType &
532  block(const unsigned int i) const;
533 
539  const BlockIndices &
540  get_block_indices() const;
541 
545  unsigned int
546  n_blocks() const;
547 
552  std::size_t
553  size() const;
554 
571  IndexSet
572  locally_owned_elements() const;
573 
577  iterator
578  begin();
579 
585  begin() const;
586 
590  iterator
591  end();
592 
598  end() const;
599 
603  value_type
604  operator()(const size_type i) const;
605 
609  reference
610  operator()(const size_type i);
611 
617  value_type operator[](const size_type i) const;
618 
624  reference operator[](const size_type i);
625 
641  template <typename OtherNumber>
642  void
643  extract_subvector_to(const std::vector<size_type> &indices,
644  std::vector<OtherNumber> & values) const;
645 
673  template <typename ForwardIterator, typename OutputIterator>
674  void
675  extract_subvector_to(ForwardIterator indices_begin,
676  const ForwardIterator indices_end,
677  OutputIterator values_begin) const;
678 
684  operator=(const value_type s);
685 
690  operator=(const BlockVectorBase &V);
691 
698  operator=(BlockVectorBase && /*V*/) = default; // NOLINT
699 
703  template <class VectorType2>
705  operator=(const BlockVectorBase<VectorType2> &V);
706 
711  operator=(const VectorType &v);
712 
717  template <class VectorType2>
718  bool
719  operator==(const BlockVectorBase<VectorType2> &v) const;
720 
724  value_type operator*(const BlockVectorBase &V) const;
725 
729  real_type
730  norm_sqr() const;
731 
735  value_type
736  mean_value() const;
737 
741  real_type
742  l1_norm() const;
743 
748  real_type
749  l2_norm() const;
750 
755  real_type
756  linfty_norm() const;
757 
780  value_type
781  add_and_dot(const value_type a,
782  const BlockVectorBase &V,
783  const BlockVectorBase &W);
784 
789  bool
790  in_local_range(const size_type global_index) const;
791 
797  bool
798  all_zero() const;
799 
805  bool
806  is_non_negative() const;
807 
812  operator+=(const BlockVectorBase &V);
813 
818  operator-=(const BlockVectorBase &V);
819 
820 
825  template <typename Number>
826  void
827  add(const std::vector<size_type> &indices, const std::vector<Number> &values);
828 
833  template <typename Number>
834  void
835  add(const std::vector<size_type> &indices, const Vector<Number> &values);
836 
842  template <typename Number>
843  void
844  add(const size_type n_elements,
845  const size_type *indices,
846  const Number * values);
847 
852  void
853  add(const value_type s);
854 
858  void
859  add(const value_type a, const BlockVectorBase &V);
860 
864  void
865  add(const value_type a,
866  const BlockVectorBase &V,
867  const value_type b,
868  const BlockVectorBase &W);
869 
873  void
874  sadd(const value_type s, const BlockVectorBase &V);
875 
879  void
880  sadd(const value_type s, const value_type a, const BlockVectorBase &V);
881 
885  void
886  sadd(const value_type s,
887  const value_type a,
888  const BlockVectorBase &V,
889  const value_type b,
890  const BlockVectorBase &W);
891 
895  void
896  sadd(const value_type s,
897  const value_type a,
898  const BlockVectorBase &V,
899  const value_type b,
900  const BlockVectorBase &W,
901  const value_type c,
902  const BlockVectorBase &X);
903 
908  operator*=(const value_type factor);
909 
914  operator/=(const value_type factor);
915 
920  template <class BlockVector2>
921  void
922  scale(const BlockVector2 &v);
923 
927  template <class BlockVector2>
928  void
929  equ(const value_type a, const BlockVector2 &V);
930 
935  void
936  update_ghost_values() const;
937 
942  std::size_t
943  memory_consumption() const;
944 
945 protected:
949  std::vector<VectorType> components;
950 
955  BlockIndices block_indices;
956 
957  // Make the iterator class a friend.
958  template <typename N, bool C>
959  friend class ::internal::BlockVectorIterators::Iterator;
960 
961  template <typename>
962  friend class BlockVectorBase;
963 };
964 
965 
968 /*----------------------- Inline functions ----------------------------------*/
969 
970 
971 #ifndef DOXYGEN
972 namespace internal
973 {
974  namespace BlockVectorIterators
975  {
976  template <class BlockVectorType, bool Constness>
977  inline Iterator<BlockVectorType, Constness>::Iterator(
978  const Iterator<BlockVectorType, Constness> &c)
979  : parent(c.parent)
980  , global_index(c.global_index)
981  , current_block(c.current_block)
982  , index_within_block(c.index_within_block)
983  , next_break_forward(c.next_break_forward)
984  , next_break_backward(c.next_break_backward)
985  {}
986 
987 
988 
989  template <class BlockVectorType, bool Constness>
990  inline Iterator<BlockVectorType, Constness>::Iterator(
991  const Iterator<BlockVectorType, !Constness> &c)
992  : parent(c.parent)
993  , global_index(c.global_index)
994  , current_block(c.current_block)
995  , index_within_block(c.index_within_block)
996  , next_break_forward(c.next_break_forward)
997  , next_break_backward(c.next_break_backward)
998  {
999  // Only permit copy-constructing const iterators from non-const
1000  // iterators, and not vice versa (i.e., Constness must always be
1001  // true).
1002  static_assert(Constness == true,
1003  "Constructing a non-const iterator from a const iterator "
1004  "does not make sense.");
1005  }
1006 
1007 
1008 
1009  template <class BlockVectorType, bool Constness>
1010  inline Iterator<BlockVectorType, Constness>::Iterator(
1011  BlockVector & parent,
1012  const size_type global_index,
1013  const size_type current_block,
1014  const size_type index_within_block,
1015  const size_type next_break_forward,
1016  const size_type next_break_backward)
1017  : parent(&parent)
1018  , global_index(global_index)
1019  , current_block(current_block)
1020  , index_within_block(index_within_block)
1021  , next_break_forward(next_break_forward)
1022  , next_break_backward(next_break_backward)
1023  {}
1024 
1025 
1026 
1027  template <class BlockVectorType, bool Constness>
1028  inline Iterator<BlockVectorType, Constness> &
1029  Iterator<BlockVectorType, Constness>::operator=(const Iterator &c)
1030  {
1031  parent = c.parent;
1032  global_index = c.global_index;
1033  index_within_block = c.index_within_block;
1034  current_block = c.current_block;
1035  next_break_forward = c.next_break_forward;
1036  next_break_backward = c.next_break_backward;
1037 
1038  return *this;
1039  }
1040 
1041 
1042 
1043  template <class BlockVectorType, bool Constness>
1044  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1046  {
1047  return parent->block(current_block)(index_within_block);
1048  }
1049 
1050 
1051 
1052  template <class BlockVectorType, bool Constness>
1053  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1054  Iterator<BlockVectorType, Constness>::
1055  operator[](const difference_type d) const
1056  {
1057  // if the index pointed to is
1058  // still within the block we
1059  // currently point into, then we
1060  // can save the computation of
1061  // the block
1062  if ((global_index + d >= next_break_backward) &&
1063  (global_index + d <= next_break_forward))
1064  return parent->block(current_block)(index_within_block + d);
1065 
1066  // if the index is not within the
1067  // block of the block vector into
1068  // which we presently point, then
1069  // there is no way: we have to
1070  // search for the block. this can
1071  // be done through the parent
1072  // class as well.
1073  return (*parent)(global_index + d);
1074  }
1075 
1076 
1077 
1078  template <class BlockVectorType, bool Constness>
1079  inline Iterator<BlockVectorType, Constness> &
1081  {
1082  move_forward();
1083  return *this;
1084  }
1085 
1086 
1087 
1088  template <class BlockVectorType, bool Constness>
1089  inline Iterator<BlockVectorType, Constness>
1091  {
1092  const Iterator old_value = *this;
1093  move_forward();
1094  return old_value;
1095  }
1096 
1097 
1098 
1099  template <class BlockVectorType, bool Constness>
1100  inline Iterator<BlockVectorType, Constness> &
1101  Iterator<BlockVectorType, Constness>::operator--()
1102  {
1103  move_backward();
1104  return *this;
1105  }
1106 
1107 
1108 
1109  template <class BlockVectorType, bool Constness>
1110  inline Iterator<BlockVectorType, Constness>
1111  Iterator<BlockVectorType, Constness>::operator--(int)
1112  {
1113  const Iterator old_value = *this;
1114  move_backward();
1115  return old_value;
1116  }
1117 
1118 
1119 
1120  template <class BlockVectorType, bool Constness>
1121  template <bool OtherConstness>
1122  inline bool
1124  operator==(const Iterator<BlockVectorType, OtherConstness> &i) const
1125  {
1126  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1127 
1128  return (global_index == i.global_index);
1129  }
1130 
1131 
1132 
1133  template <class BlockVectorType, bool Constness>
1134  template <bool OtherConstness>
1135  inline bool
1137  operator!=(const Iterator<BlockVectorType, OtherConstness> &i) const
1138  {
1139  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1140 
1141  return (global_index != i.global_index);
1142  }
1143 
1144 
1145 
1146  template <class BlockVectorType, bool Constness>
1147  template <bool OtherConstness>
1148  inline bool
1150  operator<(const Iterator<BlockVectorType, OtherConstness> &i) const
1151  {
1152  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1153 
1154  return (global_index < i.global_index);
1155  }
1156 
1157 
1158 
1159  template <class BlockVectorType, bool Constness>
1160  template <bool OtherConstness>
1161  inline bool
1163  operator<=(const Iterator<BlockVectorType, OtherConstness> &i) const
1164  {
1165  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1166 
1167  return (global_index <= i.global_index);
1168  }
1169 
1170 
1171 
1172  template <class BlockVectorType, bool Constness>
1173  template <bool OtherConstness>
1174  inline bool
1176  operator>(const Iterator<BlockVectorType, OtherConstness> &i) const
1177  {
1178  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1179 
1180  return (global_index > i.global_index);
1181  }
1182 
1183 
1184 
1185  template <class BlockVectorType, bool Constness>
1186  template <bool OtherConstness>
1187  inline bool
1189  operator>=(const Iterator<BlockVectorType, OtherConstness> &i) const
1190  {
1191  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1192 
1193  return (global_index >= i.global_index);
1194  }
1195 
1196 
1197 
1198  template <class BlockVectorType, bool Constness>
1199  template <bool OtherConstness>
1200  inline typename Iterator<BlockVectorType, Constness>::difference_type
1202  operator-(const Iterator<BlockVectorType, OtherConstness> &i) const
1203  {
1204  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1205 
1206  return (static_cast<signed int>(global_index) -
1207  static_cast<signed int>(i.global_index));
1208  }
1209 
1210 
1211 
1212  template <class BlockVectorType, bool Constness>
1213  inline Iterator<BlockVectorType, Constness>
1215  operator+(const difference_type &d) const
1216  {
1217  // if the index pointed to is
1218  // still within the block we
1219  // currently point into, then we
1220  // can save the computation of
1221  // the block
1222  if ((global_index + d >= next_break_backward) &&
1223  (global_index + d <= next_break_forward))
1224  return Iterator(*parent,
1225  global_index + d,
1226  current_block,
1227  index_within_block + d,
1228  next_break_forward,
1229  next_break_backward);
1230  else
1231  // outside present block, so
1232  // have to seek new block
1233  // anyway
1234  return Iterator(*parent, global_index + d);
1235  }
1236 
1237 
1238 
1239  template <class BlockVectorType, bool Constness>
1240  inline Iterator<BlockVectorType, Constness>
1242  operator-(const difference_type &d) const
1243  {
1244  // if the index pointed to is
1245  // still within the block we
1246  // currently point into, then we
1247  // can save the computation of
1248  // the block
1249  if ((global_index - d >= next_break_backward) &&
1250  (global_index - d <= next_break_forward))
1251  return Iterator(*parent,
1252  global_index - d,
1253  current_block,
1254  index_within_block - d,
1255  next_break_forward,
1256  next_break_backward);
1257  else
1258  // outside present block, so
1259  // have to seek new block
1260  // anyway
1261  return Iterator(*parent, global_index - d);
1262  }
1263 
1264 
1265 
1266  template <class BlockVectorType, bool Constness>
1267  inline Iterator<BlockVectorType, Constness> &
1268  Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1269  {
1270  // if the index pointed to is
1271  // still within the block we
1272  // currently point into, then we
1273  // can save the computation of
1274  // the block
1275  if ((global_index + d >= next_break_backward) &&
1276  (global_index + d <= next_break_forward))
1277  {
1278  global_index += d;
1279  index_within_block += d;
1280  }
1281  else
1282  // outside present block, so
1283  // have to seek new block
1284  // anyway
1285  *this = Iterator(*parent, global_index + d);
1286 
1287  return *this;
1288  }
1289 
1290 
1291 
1292  template <class BlockVectorType, bool Constness>
1293  inline Iterator<BlockVectorType, Constness> &
1294  Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1295  {
1296  // if the index pointed to is
1297  // still within the block we
1298  // currently point into, then we
1299  // can save the computation of
1300  // the block
1301  if ((global_index - d >= next_break_backward) &&
1302  (global_index - d <= next_break_forward))
1303  {
1304  global_index -= d;
1305  index_within_block -= d;
1306  }
1307  else
1308  // outside present block, so
1309  // have to seek new block
1310  // anyway
1311  *this = Iterator(*parent, global_index - d);
1312 
1313  return *this;
1314  }
1315 
1316 
1317  template <class BlockVectorType, bool Constness>
1318  Iterator<BlockVectorType, Constness>::Iterator(BlockVector & parent,
1319  const size_type global_index)
1320  : parent(&parent)
1321  , global_index(global_index)
1322  {
1323  // find which block we are
1324  // in. for this, take into
1325  // account that it happens at
1326  // times that people want to
1327  // initialize iterators
1328  // past-the-end
1329  if (global_index < parent.size())
1330  {
1331  const std::pair<size_type, size_type> indices =
1332  parent.block_indices.global_to_local(global_index);
1333  current_block = indices.first;
1334  index_within_block = indices.second;
1335 
1336  next_break_backward =
1337  parent.block_indices.local_to_global(current_block, 0);
1338  next_break_forward =
1339  (parent.block_indices.local_to_global(current_block, 0) +
1340  parent.block_indices.block_size(current_block) - 1);
1341  }
1342  else
1343  // past the end. only have one
1344  // value for this
1345  {
1346  this->global_index = parent.size();
1347  current_block = parent.n_blocks();
1348  index_within_block = 0;
1349  next_break_backward = global_index;
1350  next_break_forward = numbers::invalid_size_type;
1351  };
1352  }
1353 
1354 
1355 
1356  template <class BlockVectorType, bool Constness>
1357  void
1358  Iterator<BlockVectorType, Constness>::move_forward()
1359  {
1360  if (global_index != next_break_forward)
1361  ++index_within_block;
1362  else
1363  {
1364  // ok, we traverse a boundary
1365  // between blocks:
1366  index_within_block = 0;
1367  ++current_block;
1368 
1369  // break backwards is now old
1370  // break forward
1371  next_break_backward = next_break_forward + 1;
1372 
1373  // compute new break forward
1374  if (current_block < parent->block_indices.size())
1375  next_break_forward +=
1376  parent->block_indices.block_size(current_block);
1377  else
1378  // if we are beyond the end,
1379  // then move the next
1380  // boundary arbitrarily far
1381  // away
1382  next_break_forward = numbers::invalid_size_type;
1383  };
1384 
1385  ++global_index;
1386  }
1387 
1388 
1389 
1390  template <class BlockVectorType, bool Constness>
1391  void
1392  Iterator<BlockVectorType, Constness>::move_backward()
1393  {
1394  if (global_index != next_break_backward)
1395  --index_within_block;
1396  else if (current_block != 0)
1397  {
1398  // ok, we traverse a boundary
1399  // between blocks:
1400  --current_block;
1401  index_within_block =
1402  parent->block_indices.block_size(current_block) - 1;
1403 
1404  // break forwards is now old
1405  // break backward
1406  next_break_forward = next_break_backward - 1;
1407 
1408  // compute new break forward
1409  next_break_backward -=
1410  parent->block_indices.block_size(current_block);
1411  }
1412  else
1413  // current block was 0, we now
1414  // get into unspecified terrain
1415  {
1416  --current_block;
1417  index_within_block = numbers::invalid_size_type;
1418  next_break_forward = 0;
1419  next_break_backward = 0;
1420  };
1421 
1422  --global_index;
1423  }
1424 
1425 
1426  } // namespace BlockVectorIterators
1427 
1428 } // namespace internal
1429 
1430 
1431 
1432 template <class VectorType>
1433 inline std::size_t
1435 {
1436  return block_indices.total_size();
1437 }
1438 
1439 
1440 
1441 template <class VectorType>
1442 inline IndexSet
1444 {
1445  IndexSet is(size());
1446 
1447  // copy index sets from blocks into the global one, shifted
1448  // by the appropriate amount for each block
1449  for (unsigned int b = 0; b < n_blocks(); ++b)
1450  {
1451  IndexSet x = block(b).locally_owned_elements();
1452  is.add_indices(x, block_indices.block_start(b));
1453  }
1454 
1455  is.compress();
1456 
1457  return is;
1458 }
1459 
1460 
1461 
1462 template <class VectorType>
1463 inline unsigned int
1465 {
1466  return block_indices.size();
1467 }
1468 
1469 
1470 template <class VectorType>
1472 BlockVectorBase<VectorType>::block(const unsigned int i)
1473 {
1474  AssertIndexRange(i, n_blocks());
1475 
1476  return components[i];
1477 }
1478 
1479 
1480 
1481 template <class VectorType>
1482 inline const typename BlockVectorBase<VectorType>::BlockType &
1483 BlockVectorBase<VectorType>::block(const unsigned int i) const
1484 {
1485  AssertIndexRange(i, n_blocks());
1486 
1487  return components[i];
1488 }
1489 
1490 
1491 
1492 template <class VectorType>
1493 inline const BlockIndices &
1495 {
1496  return block_indices;
1497 }
1498 
1499 
1500 template <class VectorType>
1501 inline void
1503 {
1504  std::vector<size_type> sizes(n_blocks());
1505 
1506  for (size_type i = 0; i < n_blocks(); ++i)
1507  sizes[i] = block(i).size();
1508 
1509  block_indices.reinit(sizes);
1510 }
1511 
1512 
1513 
1514 template <class VectorType>
1515 inline void
1517  ::VectorOperation::values operation)
1518 {
1519  for (unsigned int i = 0; i < n_blocks(); ++i)
1520  block(i).compress(operation);
1521 }
1522 
1523 
1524 
1525 template <class VectorType>
1528 {
1529  return iterator(*this, 0U);
1530 }
1531 
1532 
1533 
1534 template <class VectorType>
1537 {
1538  return const_iterator(*this, 0U);
1539 }
1540 
1541 
1542 template <class VectorType>
1545 {
1546  return iterator(*this, size());
1547 }
1548 
1549 
1550 
1551 template <class VectorType>
1554 {
1555  return const_iterator(*this, size());
1556 }
1557 
1558 
1559 template <class VectorType>
1560 inline bool
1562 {
1563  const std::pair<size_type, size_type> local_index =
1564  block_indices.global_to_local(global_index);
1565 
1566  return components[local_index.first].in_local_range(global_index);
1567 }
1568 
1569 
1570 template <class VectorType>
1571 bool
1573 {
1574  for (size_type i = 0; i < n_blocks(); ++i)
1575  if (components[i].all_zero() == false)
1576  return false;
1577 
1578  return true;
1579 }
1580 
1581 
1582 
1583 template <class VectorType>
1584 bool
1586 {
1587  for (size_type i = 0; i < n_blocks(); ++i)
1588  if (components[i].is_non_negative() == false)
1589  return false;
1590 
1591  return true;
1592 }
1593 
1594 
1595 
1596 template <class VectorType>
1599 {
1600  Assert(n_blocks() == v.n_blocks(),
1602 
1603  value_type sum = 0.;
1604  for (size_type i = 0; i < n_blocks(); ++i)
1605  sum += components[i] * v.components[i];
1606 
1607  return sum;
1608 }
1609 
1610 
1611 template <class VectorType>
1614 {
1615  real_type sum = 0.;
1616  for (size_type i = 0; i < n_blocks(); ++i)
1617  sum += components[i].norm_sqr();
1618 
1619  return sum;
1620 }
1621 
1622 
1623 
1624 template <class VectorType>
1627 {
1628  value_type sum = 0.;
1629  // need to do static_cast as otherwise it won't work with
1630  // value_type=complex<T>
1631  for (size_type i = 0; i < n_blocks(); ++i)
1632  sum += components[i].mean_value() *
1634  components[i].size()));
1635 
1636  return sum / (typename numbers::NumberTraits<value_type>::real_type(size()));
1637 }
1638 
1639 
1640 
1641 template <class VectorType>
1644 {
1645  real_type sum = 0.;
1646  for (size_type i = 0; i < n_blocks(); ++i)
1647  sum += components[i].l1_norm();
1648 
1649  return sum;
1650 }
1651 
1652 
1653 
1654 template <class VectorType>
1657 {
1658  return std::sqrt(norm_sqr());
1659 }
1660 
1661 
1662 
1663 template <class VectorType>
1666 {
1667  real_type sum = 0.;
1668  for (size_type i = 0; i < n_blocks(); ++i)
1669  {
1670  value_type newval = components[i].linfty_norm();
1671  if (sum < newval)
1672  sum = newval;
1673  }
1674  return sum;
1675 }
1676 
1677 
1678 
1679 template <class VectorType>
1682  const typename BlockVectorBase<VectorType>::value_type a,
1683  const BlockVectorBase<VectorType> & V,
1684  const BlockVectorBase<VectorType> & W)
1685 {
1688 
1689  value_type sum = 0.;
1690  for (size_type i = 0; i < n_blocks(); ++i)
1691  sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1692 
1693  return sum;
1694 }
1695 
1696 
1697 
1698 template <class VectorType>
1701 {
1702  Assert(n_blocks() == v.n_blocks(),
1704 
1705  for (size_type i = 0; i < n_blocks(); ++i)
1706  {
1707  components[i] += v.components[i];
1708  }
1709 
1710  return *this;
1711 }
1712 
1713 
1714 
1715 template <class VectorType>
1718 {
1719  Assert(n_blocks() == v.n_blocks(),
1721 
1722  for (size_type i = 0; i < n_blocks(); ++i)
1723  {
1724  components[i] -= v.components[i];
1725  }
1726  return *this;
1727 }
1728 
1729 
1730 
1731 template <class VectorType>
1732 template <typename Number>
1733 inline void
1734 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1735  const std::vector<Number> & values)
1736 {
1737  Assert(indices.size() == values.size(),
1738  ExcDimensionMismatch(indices.size(), values.size()));
1739  add(indices.size(), indices.data(), values.data());
1740 }
1741 
1742 
1743 
1744 template <class VectorType>
1745 template <typename Number>
1746 inline void
1747 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1748  const Vector<Number> & values)
1749 {
1750  Assert(indices.size() == values.size(),
1751  ExcDimensionMismatch(indices.size(), values.size()));
1752  const size_type n_indices = indices.size();
1753  for (size_type i = 0; i < n_indices; ++i)
1754  (*this)(indices[i]) += values(i);
1755 }
1756 
1757 
1758 
1759 template <class VectorType>
1760 template <typename Number>
1761 inline void
1763  const size_type *indices,
1764  const Number * values)
1765 {
1766  for (size_type i = 0; i < n_indices; ++i)
1767  (*this)(indices[i]) += values[i];
1768 }
1769 
1770 
1771 
1772 template <class VectorType>
1773 void
1775 {
1776  AssertIsFinite(a);
1777 
1778  for (size_type i = 0; i < n_blocks(); ++i)
1779  {
1780  components[i].add(a);
1781  }
1782 }
1783 
1784 
1785 
1786 template <class VectorType>
1787 void
1789  const BlockVectorBase<VectorType> &v)
1790 {
1791  AssertIsFinite(a);
1792 
1793  Assert(n_blocks() == v.n_blocks(),
1795 
1796  for (size_type i = 0; i < n_blocks(); ++i)
1797  {
1798  components[i].add(a, v.components[i]);
1799  }
1800 }
1801 
1802 
1803 
1804 template <class VectorType>
1805 void
1807  const BlockVectorBase<VectorType> &v,
1808  const value_type b,
1810 {
1811  AssertIsFinite(a);
1812  AssertIsFinite(b);
1813 
1814  Assert(n_blocks() == v.n_blocks(),
1816  Assert(n_blocks() == w.n_blocks(),
1818 
1819 
1820  for (size_type i = 0; i < n_blocks(); ++i)
1821  {
1822  components[i].add(a, v.components[i], b, w.components[i]);
1823  }
1824 }
1825 
1826 
1827 
1828 template <class VectorType>
1829 void
1831  const BlockVectorBase<VectorType> &v)
1832 {
1833  AssertIsFinite(x);
1834 
1835  Assert(n_blocks() == v.n_blocks(),
1837 
1838  for (size_type i = 0; i < n_blocks(); ++i)
1839  {
1840  components[i].sadd(x, v.components[i]);
1841  }
1842 }
1843 
1844 
1845 
1846 template <class VectorType>
1847 void
1849  const value_type a,
1850  const BlockVectorBase<VectorType> &v)
1851 {
1852  AssertIsFinite(x);
1853  AssertIsFinite(a);
1854 
1855  Assert(n_blocks() == v.n_blocks(),
1857 
1858  for (size_type i = 0; i < n_blocks(); ++i)
1859  {
1860  components[i].sadd(x, a, v.components[i]);
1861  }
1862 }
1863 
1864 
1865 
1866 template <class VectorType>
1867 void
1869  const value_type a,
1870  const BlockVectorBase<VectorType> &v,
1871  const value_type b,
1872  const BlockVectorBase<VectorType> &w)
1873 {
1874  AssertIsFinite(x);
1875  AssertIsFinite(a);
1876  AssertIsFinite(b);
1877 
1878  Assert(n_blocks() == v.n_blocks(),
1880  Assert(n_blocks() == w.n_blocks(),
1882 
1883  for (size_type i = 0; i < n_blocks(); ++i)
1884  {
1885  components[i].sadd(x, a, v.components[i], b, w.components[i]);
1886  }
1887 }
1888 
1889 
1890 
1891 template <class VectorType>
1892 void
1894  const value_type a,
1895  const BlockVectorBase<VectorType> &v,
1896  const value_type b,
1897  const BlockVectorBase<VectorType> &w,
1898  const value_type c,
1899  const BlockVectorBase<VectorType> &y)
1900 {
1901  AssertIsFinite(x);
1902  AssertIsFinite(a);
1903  AssertIsFinite(b);
1904  AssertIsFinite(c);
1905 
1906  Assert(n_blocks() == v.n_blocks(),
1908  Assert(n_blocks() == w.n_blocks(),
1910  Assert(n_blocks() == y.n_blocks(),
1912 
1913  for (size_type i = 0; i < n_blocks(); ++i)
1914  {
1915  components[i].sadd(
1916  x, a, v.components[i], b, w.components[i], c, y.components[i]);
1917  }
1918 }
1919 
1920 
1921 
1922 template <class VectorType>
1923 template <class BlockVector2>
1924 void
1925 BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
1926 {
1927  Assert(n_blocks() == v.n_blocks(),
1928  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1929  for (size_type i = 0; i < n_blocks(); ++i)
1930  components[i].scale(v.block(i));
1931 }
1932 
1933 
1934 
1935 template <class VectorType>
1936 std::size_t
1938 {
1939  return (MemoryConsumption::memory_consumption(this->block_indices) +
1940  MemoryConsumption::memory_consumption(this->components));
1941 }
1942 
1943 
1944 
1945 template <class VectorType>
1946 template <class BlockVector2>
1947 void
1948 BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
1949 {
1950  AssertIsFinite(a);
1951 
1952  Assert(n_blocks() == v.n_blocks(),
1953  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1954 
1955  for (size_type i = 0; i < n_blocks(); ++i)
1956  components[i].equ(a, v.components[i]);
1957 }
1958 
1959 
1960 
1961 template <class VectorType>
1962 void
1964 {
1965  for (size_type i = 0; i < n_blocks(); ++i)
1966  block(i).update_ghost_values();
1967 }
1968 
1969 
1970 
1971 template <class VectorType>
1974 {
1975  AssertIsFinite(s);
1976 
1977  for (size_type i = 0; i < n_blocks(); ++i)
1978  components[i] = s;
1979 
1980  return *this;
1981 }
1982 
1983 
1984 template <class VectorType>
1987 {
1989 
1990  for (size_type i = 0; i < n_blocks(); ++i)
1991  components[i] = v.components[i];
1992 
1993  return *this;
1994 }
1995 
1996 
1997 template <class VectorType>
1998 template <class VectorType2>
2001 {
2003 
2004  for (size_type i = 0; i < n_blocks(); ++i)
2005  components[i] = v.components[i];
2006 
2007  return *this;
2008 }
2009 
2010 
2011 
2012 template <class VectorType>
2015 {
2016  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2017 
2018  size_type index_v = 0;
2019  for (size_type b = 0; b < n_blocks(); ++b)
2020  for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2021  block(b)(i) = v(index_v);
2022 
2023  return *this;
2024 }
2025 
2026 
2027 
2028 template <class VectorType>
2029 template <class VectorType2>
2030 inline bool
2033 {
2034  Assert(block_indices == v.block_indices, ExcDifferentBlockIndices());
2035 
2036  for (size_type i = 0; i < n_blocks(); ++i)
2037  if (!(components[i] == v.components[i]))
2038  return false;
2039 
2040  return true;
2041 }
2042 
2043 
2044 
2045 template <class VectorType>
2048 {
2049  AssertIsFinite(factor);
2050 
2051  for (size_type i = 0; i < n_blocks(); ++i)
2052  components[i] *= factor;
2053 
2054  return *this;
2055 }
2056 
2057 
2058 
2059 template <class VectorType>
2062 {
2063  AssertIsFinite(factor);
2064  Assert(factor != 0., ExcDivideByZero());
2065 
2066  const value_type inverse_factor = value_type(1.) / factor;
2067 
2068  for (size_type i = 0; i < n_blocks(); ++i)
2069  components[i] *= inverse_factor;
2070 
2071  return *this;
2072 }
2073 
2074 
2075 template <class VectorType>
2078 {
2079  const std::pair<unsigned int, size_type> local_index =
2080  block_indices.global_to_local(i);
2081  return components[local_index.first](local_index.second);
2082 }
2083 
2084 
2085 
2086 template <class VectorType>
2089 {
2090  const std::pair<unsigned int, size_type> local_index =
2091  block_indices.global_to_local(i);
2092  return components[local_index.first](local_index.second);
2093 }
2094 
2095 
2096 
2097 template <class VectorType>
2100 {
2101  return operator()(i);
2102 }
2103 
2104 
2105 
2106 template <class VectorType>
2109 {
2110  return operator()(i);
2111 }
2112 
2113 
2114 
2115 template <typename VectorType>
2116 template <typename OtherNumber>
2117 inline void
2119  const std::vector<size_type> &indices,
2120  std::vector<OtherNumber> & values) const
2121 {
2122  for (size_type i = 0; i < indices.size(); ++i)
2123  values[i] = operator()(indices[i]);
2124 }
2125 
2126 
2127 
2128 template <typename VectorType>
2129 template <typename ForwardIterator, typename OutputIterator>
2130 inline void
2132  ForwardIterator indices_begin,
2133  const ForwardIterator indices_end,
2134  OutputIterator values_begin) const
2135 {
2136  while (indices_begin != indices_end)
2137  {
2138  *values_begin = operator()(*indices_begin);
2139  indices_begin++;
2140  values_begin++;
2141  }
2142 }
2143 
2144 #endif // DOXYGEN
2145 
2147 
2148 #endif
const types::global_dof_index invalid_size_type
Definition: types.h:205
static yes_type check_for_block_vector(const BlockVectorBase< T > *)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1568
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:835
bool in_local_range(const size_type global_index) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1636
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1701
typename std::conditional< Constness, value_type, typename BlockVectorType::BlockType::reference >::type dereference_type
__global__ void add_and_dot(Number *res, Number *v1, const Number *v2, const Number *v3, const Number a, const size_type N)
static const char V
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
void collect_sizes()
STL namespace.
static const char U
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
std::enable_if< IsBlockVector< VectorType >::value, void >::type collect_sizes(VectorType &vector)
Definition: operators.h:97
size_type block_size(const unsigned int i) const
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
static ::ExceptionBase & ExcDivideByZero()
typename std::conditional< Constness, const typename BlockVectorType::value_type, typename BlockVectorType::value_type >::type value_type
BlockVectorBase & operator/=(const value_type factor)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
value_type add_and_dot(const value_type a, const BlockVectorBase &V, const BlockVectorBase &W)
size_type total_size() const
bool is_non_negative() const
void scale(const BlockVector2 &v)
std::random_access_iterator_tag iterator_category
BlockIndices block_indices
IndexSet locally_owned_elements() const
real_type norm_sqr() const
std::size_t size() const
std::string compress(const std::string &input)
Definition: utilities.cc:392
number value_type
Definition: vector.h:122
T sum(const T &t, const MPI_Comm &mpi_communicator)
Expression operator<=(const Expression &lhs, const Expression &rhs)
const BlockIndices & get_block_indices() const
#define Assert(cond, exc)
Definition: exceptions.h:1411
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
void reinit(const unsigned int n_blocks, const size_type n_elements_per_block)
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2811
real_type l1_norm() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
VectorType::value_type * end(VectorType &V)
void update_ghost_values() const
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
size_type local_to_global(const unsigned int block, const size_type index) const
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
Expression operator>(const Expression &lhs, const Expression &rhs)
static ::ExceptionBase & ExcDifferentBlockIndices()
static const bool value
__global__ void equ(Number *val, const Number a, const Number *V_val, const size_type N)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
BlockVectorBase & operator+=(const BlockVectorBase &V)
std::vector< VectorType > components
__global__ void sadd(const Number s, Number *val, const Number a, const Number *V_val, const size_type N)
BlockVectorBase & operator*=(const value_type factor)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
bool all_zero() const
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int global_dof_index
Definition: types.h:76
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< OtherNumber > &values) const
void sadd(const value_type s, const BlockVectorBase &V)
std::size_t memory_consumption() const
void compress() const
Definition: index_set.h:1640
unsigned int n_blocks() const
BlockVectorBase & operator-=(const BlockVectorBase &V)
typename BlockType::value_type value_type
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
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:2785
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:1667
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