Reference documentation for deal.II version Git 2375602 2018-06-23 12:24:17 -0500
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((VectorType *)nullptr)) == sizeof(yes_type));
100 };
101 
102 
103 // instantiation of the static member
104 template <typename VectorType>
106 
107 
108 
109 namespace internal
110 {
116  namespace BlockVectorIterators
117  {
122  template <class BlockVectorType, bool Constness>
123  struct Types
124  {};
125 
126 
127 
134  template <class BlockVectorType>
135  struct Types<BlockVectorType, false>
136  {
141  typedef typename BlockVectorType::BlockType Vector;
142 
147  typedef BlockVectorType BlockVector;
148 
153 
158  typedef typename Vector::reference dereference_type;
159  };
160 
161 
162 
169  template <class BlockVectorType>
170  struct Types<BlockVectorType, true>
171  {
176  typedef const typename BlockVectorType::BlockType Vector;
177 
182  typedef const BlockVectorType BlockVector;
183 
188  typedef const typename BlockVector::value_type value_type;
189 
195  typedef value_type dereference_type;
196  };
197 
198 
216  template <class BlockVectorType, bool Constness>
217  class Iterator
218  {
219  public:
224 
231 
237  typedef std::random_access_iterator_tag iterator_category;
238  typedef std::ptrdiff_t difference_type;
239  typedef typename BlockVectorType::reference reference;
240  typedef value_type * pointer;
241 
243  dereference_type;
244 
249  typedef
251 
260  Iterator(BlockVector &parent, const size_type global_index);
261 
270 
271 
275  Iterator(const Iterator &c);
276 
277  private:
282  Iterator(BlockVector & parent,
283  const size_type global_index,
284  const size_type current_block,
285  const size_type index_within_block,
286  const size_type next_break_forward,
287  const size_type next_break_backward);
288 
289  public:
293  Iterator &
294  operator=(const Iterator &c);
295 
301  dereference_type operator*() const;
302 
307  dereference_type operator[](const difference_type d) const;
308 
313  Iterator &
314  operator++();
315 
321  Iterator
322  operator++(int);
323 
328  Iterator &
329  operator--();
330 
336  Iterator
337  operator--(int);
338 
343  template <bool OtherConstness>
344  bool
346 
351  template <bool OtherConstness>
352  bool
354 
360  template <bool OtherConstness>
361  bool
362  operator<(const Iterator<BlockVectorType, OtherConstness> &i) const;
363 
367  template <bool OtherConstness>
368  bool
369  operator<=(const Iterator<BlockVectorType, OtherConstness> &i) const;
370 
374  template <bool OtherConstness>
375  bool
377 
381  template <bool OtherConstness>
382  bool
384 
388  template <bool OtherConstness>
389  difference_type
391 
396  Iterator
397  operator+(const difference_type &d) const;
398 
403  Iterator
404  operator-(const difference_type &d) const;
405 
410  Iterator &
411  operator+=(const difference_type &d);
412 
417  Iterator &
418  operator-=(const difference_type &d);
419 
430  "Your program tried to compare iterators pointing to "
431  "different block vectors. There is no reasonable way "
432  "to do this.");
433 
435  private:
441  BlockVector *parent;
442 
446  size_type global_index;
447 
452  unsigned int current_block;
453  size_type index_within_block;
454 
462  size_type next_break_backward;
463 
467  void
468  move_forward();
469 
473  void
474  move_backward();
475 
476 
480  template <typename, bool>
481  friend class Iterator;
482  };
483  } // namespace BlockVectorIterators
484 } // namespace internal
485 
486 
525 template <class VectorType>
526 class BlockVectorBase : public Subscriptor
527 {
528 public:
532  typedef VectorType BlockType;
533 
534  /*
535  * Declare standard types used in
536  * all containers. These types
537  * parallel those in the
538  * <tt>C++</tt> standard
539  * libraries
540  * <tt>std::vector<...></tt>
541  * class. This includes iterator
542  * types.
543  */
544  typedef typename BlockType::value_type value_type;
545  typedef value_type * pointer;
546  typedef const value_type * const_pointer;
547  typedef ::internal::BlockVectorIterators::Iterator<BlockVectorBase,
548  false>
549  iterator;
550  typedef ::internal::BlockVectorIterators::Iterator<BlockVectorBase,
551  true>
553  typedef typename BlockType::reference reference;
554  typedef typename BlockType::const_reference const_reference;
555  typedef types::global_dof_index size_type;
556 
566  typedef typename BlockType::real_type real_type;
567 
571  BlockVectorBase() = default;
572 
576  BlockVectorBase(const BlockVectorBase & /*V*/) = default;
577 
583  BlockVectorBase(BlockVectorBase && /*V*/) noexcept = default;
584 
591  void
592  collect_sizes();
593 
604  void
605  compress(::VectorOperation::values operation);
606 
610  BlockType &
611  block(const unsigned int i);
612 
616  const BlockType &
617  block(const unsigned int i) const;
618 
624  const BlockIndices &
625  get_block_indices() const;
626 
630  unsigned int
631  n_blocks() const;
632 
637  std::size_t
638  size() const;
639 
656  IndexSet
657  locally_owned_elements() const;
658 
662  iterator
663  begin();
664 
669  const_iterator
670  begin() const;
671 
675  iterator
676  end();
677 
682  const_iterator
683  end() const;
684 
688  value_type
689  operator()(const size_type i) const;
690 
694  reference
695  operator()(const size_type i);
696 
702  value_type operator[](const size_type i) const;
703 
709  reference operator[](const size_type i);
710 
726  template <typename OtherNumber>
727  void
728  extract_subvector_to(const std::vector<size_type> &indices,
729  std::vector<OtherNumber> & values) const;
730 
758  template <typename ForwardIterator, typename OutputIterator>
759  void
760  extract_subvector_to(ForwardIterator indices_begin,
761  const ForwardIterator indices_end,
762  OutputIterator values_begin) const;
763 
768  BlockVectorBase &
769  operator=(const value_type s);
770 
774  BlockVectorBase &
775  operator=(const BlockVectorBase &V);
776 
782  BlockVectorBase &
783  operator=(BlockVectorBase && /*V*/) = default; // NOLINT
784 
788  template <class VectorType2>
789  BlockVectorBase &
790  operator=(const BlockVectorBase<VectorType2> &V);
791 
795  BlockVectorBase &
796  operator=(const VectorType &v);
797 
802  template <class VectorType2>
803  bool
804  operator==(const BlockVectorBase<VectorType2> &v) const;
805 
809  value_type operator*(const BlockVectorBase &V) const;
810 
814  real_type
815  norm_sqr() const;
816 
820  value_type
821  mean_value() const;
822 
826  real_type
827  l1_norm() const;
828 
833  real_type
834  l2_norm() const;
835 
840  real_type
841  linfty_norm() const;
842 
865  value_type
866  add_and_dot(const value_type a,
867  const BlockVectorBase &V,
868  const BlockVectorBase &W);
869 
874  bool
875  in_local_range(const size_type global_index) const;
876 
882  bool
883  all_zero() const;
884 
890  bool
891  is_non_negative() const;
892 
896  BlockVectorBase &
897  operator+=(const BlockVectorBase &V);
898 
902  BlockVectorBase &
903  operator-=(const BlockVectorBase &V);
904 
905 
910  template <typename Number>
911  void
912  add(const std::vector<size_type> &indices, const std::vector<Number> &values);
913 
918  template <typename Number>
919  void
920  add(const std::vector<size_type> &indices, const Vector<Number> &values);
921 
927  template <typename Number>
928  void
929  add(const size_type n_elements,
930  const size_type *indices,
931  const Number * values);
932 
937  void
938  add(const value_type s);
939 
943  void
944  add(const value_type a, const BlockVectorBase &V);
945 
949  void
950  add(const value_type a,
951  const BlockVectorBase &V,
952  const value_type b,
953  const BlockVectorBase &W);
954 
958  void
959  sadd(const value_type s, const BlockVectorBase &V);
960 
964  void
965  sadd(const value_type s, const value_type a, const BlockVectorBase &V);
966 
970  void
971  sadd(const value_type s,
972  const value_type a,
973  const BlockVectorBase &V,
974  const value_type b,
975  const BlockVectorBase &W);
976 
980  void
981  sadd(const value_type s,
982  const value_type a,
983  const BlockVectorBase &V,
984  const value_type b,
985  const BlockVectorBase &W,
986  const value_type c,
987  const BlockVectorBase &X);
988 
992  BlockVectorBase &
993  operator*=(const value_type factor);
994 
998  BlockVectorBase &
999  operator/=(const value_type factor);
1000 
1005  template <class BlockVector2>
1006  void
1007  scale(const BlockVector2 &v);
1008 
1012  template <class BlockVector2>
1013  void
1014  equ(const value_type a, const BlockVector2 &V);
1015 
1019  void
1020  equ(const value_type a,
1021  const BlockVectorBase &V,
1022  const value_type b,
1023  const BlockVectorBase &W);
1024 
1029  void
1030  update_ghost_values() const;
1031 
1036  std::size_t
1037  memory_consumption() const;
1038 
1039 protected:
1043  std::vector<VectorType> components;
1044 
1050 
1054  template <typename N, bool C>
1055  friend class ::internal::BlockVectorIterators::Iterator;
1056 
1057  template <typename>
1058  friend class BlockVectorBase;
1059 };
1060 
1061 
1064 /*----------------------- Inline functions ----------------------------------*/
1065 
1066 
1067 #ifndef DOXYGEN
1068 namespace internal
1069 {
1070  namespace BlockVectorIterators
1071  {
1072  template <class BlockVectorType, bool Constness>
1073  inline Iterator<BlockVectorType, Constness>::Iterator(
1074  const Iterator<BlockVectorType, Constness> &c)
1075  : parent(c.parent)
1076  , global_index(c.global_index)
1077  , current_block(c.current_block)
1078  , index_within_block(c.index_within_block)
1079  , next_break_forward(c.next_break_forward)
1080  , next_break_backward(c.next_break_backward)
1081  {}
1082 
1083 
1084 
1085  template <class BlockVectorType, bool Constness>
1086  inline Iterator<BlockVectorType, Constness>::Iterator(
1087  const Iterator<BlockVectorType, !Constness> &c)
1088  : parent(c.parent)
1089  , global_index(c.global_index)
1090  , current_block(c.current_block)
1091  , index_within_block(c.index_within_block)
1092  , next_break_forward(c.next_break_forward)
1093  , next_break_backward(c.next_break_backward)
1094  {
1095  // Only permit copy-constructing const iterators from non-const
1096  // iterators, and not vice versa (i.e., Constness must always be
1097  // true).
1098  static_assert(Constness == true,
1099  "Constructing a non-const iterator from a const iterator "
1100  "does not make sense.");
1101  }
1102 
1103 
1104 
1105  template <class BlockVectorType, bool Constness>
1106  inline Iterator<BlockVectorType, Constness>::Iterator(
1107  BlockVector & parent,
1108  const size_type global_index,
1109  const size_type current_block,
1110  const size_type index_within_block,
1111  const size_type next_break_forward,
1112  const size_type next_break_backward)
1113  : parent(&parent)
1114  , global_index(global_index)
1115  , current_block(current_block)
1116  , index_within_block(index_within_block)
1117  , next_break_forward(next_break_forward)
1118  , next_break_backward(next_break_backward)
1119  {}
1120 
1121 
1122 
1123  template <class BlockVectorType, bool Constness>
1124  inline Iterator<BlockVectorType, Constness> &
1125  Iterator<BlockVectorType, Constness>::operator=(const Iterator &c)
1126  {
1127  parent = c.parent;
1128  global_index = c.global_index;
1129  index_within_block = c.index_within_block;
1130  current_block = c.current_block;
1131  next_break_forward = c.next_break_forward;
1132  next_break_backward = c.next_break_backward;
1133 
1134  return *this;
1135  }
1136 
1137 
1138 
1139  template <class BlockVectorType, bool Constness>
1140  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1142  {
1143  return parent->block(current_block)(index_within_block);
1144  }
1145 
1146 
1147 
1148  template <class BlockVectorType, bool Constness>
1149  inline typename Iterator<BlockVectorType, Constness>::dereference_type
1150  Iterator<BlockVectorType, Constness>::
1151  operator[](const difference_type d) const
1152  {
1153  // if the index pointed to is
1154  // still within the block we
1155  // currently point into, then we
1156  // can save the computation of
1157  // the block
1158  if ((global_index + d >= next_break_backward) &&
1159  (global_index + d <= next_break_forward))
1160  return parent->block(current_block)(index_within_block + d);
1161 
1162  // if the index is not within the
1163  // block of the block vector into
1164  // which we presently point, then
1165  // there is no way: we have to
1166  // search for the block. this can
1167  // be done through the parent
1168  // class as well.
1169  return (*parent)(global_index + d);
1170  }
1171 
1172 
1173 
1174  template <class BlockVectorType, bool Constness>
1175  inline Iterator<BlockVectorType, Constness> &
1176  Iterator<BlockVectorType, Constness>::operator++()
1177  {
1178  move_forward();
1179  return *this;
1180  }
1181 
1182 
1183 
1184  template <class BlockVectorType, bool Constness>
1185  inline Iterator<BlockVectorType, Constness>
1186  Iterator<BlockVectorType, Constness>::operator++(int)
1187  {
1188  const Iterator old_value = *this;
1189  move_forward();
1190  return old_value;
1191  }
1192 
1193 
1194 
1195  template <class BlockVectorType, bool Constness>
1196  inline Iterator<BlockVectorType, Constness> &
1197  Iterator<BlockVectorType, Constness>::operator--()
1198  {
1199  move_backward();
1200  return *this;
1201  }
1202 
1203 
1204 
1205  template <class BlockVectorType, bool Constness>
1206  inline Iterator<BlockVectorType, Constness>
1207  Iterator<BlockVectorType, Constness>::operator--(int)
1208  {
1209  const Iterator old_value = *this;
1210  move_backward();
1211  return old_value;
1212  }
1213 
1214 
1215 
1216  template <class BlockVectorType, bool Constness>
1217  template <bool OtherConstness>
1218  inline bool
1219  Iterator<BlockVectorType, Constness>::
1220  operator==(const Iterator<BlockVectorType, OtherConstness> &i) const
1221  {
1222  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1223 
1224  return (global_index == i.global_index);
1225  }
1226 
1227 
1228 
1229  template <class BlockVectorType, bool Constness>
1230  template <bool OtherConstness>
1231  inline bool
1232  Iterator<BlockVectorType, Constness>::
1233  operator!=(const Iterator<BlockVectorType, OtherConstness> &i) const
1234  {
1235  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1236 
1237  return (global_index != i.global_index);
1238  }
1239 
1240 
1241 
1242  template <class BlockVectorType, bool Constness>
1243  template <bool OtherConstness>
1244  inline bool
1245  Iterator<BlockVectorType, Constness>::
1246  operator<(const Iterator<BlockVectorType, OtherConstness> &i) const
1247  {
1248  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1249 
1250  return (global_index < i.global_index);
1251  }
1252 
1253 
1254 
1255  template <class BlockVectorType, bool Constness>
1256  template <bool OtherConstness>
1257  inline bool
1258  Iterator<BlockVectorType, Constness>::
1259  operator<=(const Iterator<BlockVectorType, OtherConstness> &i) const
1260  {
1261  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1262 
1263  return (global_index <= i.global_index);
1264  }
1265 
1266 
1267 
1268  template <class BlockVectorType, bool Constness>
1269  template <bool OtherConstness>
1270  inline bool
1271  Iterator<BlockVectorType, Constness>::
1272  operator>(const Iterator<BlockVectorType, OtherConstness> &i) const
1273  {
1274  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1275 
1276  return (global_index > i.global_index);
1277  }
1278 
1279 
1280 
1281  template <class BlockVectorType, bool Constness>
1282  template <bool OtherConstness>
1283  inline bool
1284  Iterator<BlockVectorType, Constness>::
1285  operator>=(const Iterator<BlockVectorType, OtherConstness> &i) const
1286  {
1287  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1288 
1289  return (global_index >= i.global_index);
1290  }
1291 
1292 
1293 
1294  template <class BlockVectorType, bool Constness>
1295  template <bool OtherConstness>
1296  inline typename Iterator<BlockVectorType, Constness>::difference_type
1298  operator-(const Iterator<BlockVectorType, OtherConstness> &i) const
1299  {
1300  Assert(parent == i.parent, ExcPointerToDifferentVectors());
1301 
1302  return (static_cast<signed int>(global_index) -
1303  static_cast<signed int>(i.global_index));
1304  }
1305 
1306 
1307 
1308  template <class BlockVectorType, bool Constness>
1309  inline Iterator<BlockVectorType, Constness>
1311  operator+(const difference_type &d) const
1312  {
1313  // if the index pointed to is
1314  // still within the block we
1315  // currently point into, then we
1316  // can save the computation of
1317  // the block
1318  if ((global_index + d >= next_break_backward) &&
1319  (global_index + d <= next_break_forward))
1320  return Iterator(*parent,
1321  global_index + d,
1322  current_block,
1323  index_within_block + d,
1324  next_break_forward,
1325  next_break_backward);
1326  else
1327  // outside present block, so
1328  // have to seek new block
1329  // anyway
1330  return Iterator(*parent, global_index + d);
1331  }
1332 
1333 
1334 
1335  template <class BlockVectorType, bool Constness>
1336  inline Iterator<BlockVectorType, Constness>
1338  operator-(const difference_type &d) const
1339  {
1340  // if the index pointed to is
1341  // still within the block we
1342  // currently point into, then we
1343  // can save the computation of
1344  // the block
1345  if ((global_index - d >= next_break_backward) &&
1346  (global_index - d <= next_break_forward))
1347  return Iterator(*parent,
1348  global_index - d,
1349  current_block,
1350  index_within_block - d,
1351  next_break_forward,
1352  next_break_backward);
1353  else
1354  // outside present block, so
1355  // have to seek new block
1356  // anyway
1357  return Iterator(*parent, global_index - d);
1358  }
1359 
1360 
1361 
1362  template <class BlockVectorType, bool Constness>
1363  inline Iterator<BlockVectorType, Constness> &
1364  Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
1365  {
1366  // if the index pointed to is
1367  // still within the block we
1368  // currently point into, then we
1369  // can save the computation of
1370  // the block
1371  if ((global_index + d >= next_break_backward) &&
1372  (global_index + d <= next_break_forward))
1373  {
1374  global_index += d;
1375  index_within_block += d;
1376  }
1377  else
1378  // outside present block, so
1379  // have to seek new block
1380  // anyway
1381  *this = Iterator(*parent, global_index + d);
1382 
1383  return *this;
1384  }
1385 
1386 
1387 
1388  template <class BlockVectorType, bool Constness>
1389  inline Iterator<BlockVectorType, Constness> &
1390  Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
1391  {
1392  // if the index pointed to is
1393  // still within the block we
1394  // currently point into, then we
1395  // can save the computation of
1396  // the block
1397  if ((global_index - d >= next_break_backward) &&
1398  (global_index - d <= next_break_forward))
1399  {
1400  global_index -= d;
1401  index_within_block -= d;
1402  }
1403  else
1404  // outside present block, so
1405  // have to seek new block
1406  // anyway
1407  *this = Iterator(*parent, global_index - d);
1408 
1409  return *this;
1410  }
1411 
1412 
1413  template <class BlockVectorType, bool Constness>
1414  Iterator<BlockVectorType, Constness>::Iterator(BlockVector & parent,
1415  const size_type global_index)
1416  : parent(&parent)
1417  , global_index(global_index)
1418  {
1419  // find which block we are
1420  // in. for this, take into
1421  // account that it happens at
1422  // times that people want to
1423  // initialize iterators
1424  // past-the-end
1425  if (global_index < parent.size())
1426  {
1427  const std::pair<size_type, size_type> indices =
1428  parent.block_indices.global_to_local(global_index);
1429  current_block = indices.first;
1430  index_within_block = indices.second;
1431 
1432  next_break_backward =
1433  parent.block_indices.local_to_global(current_block, 0);
1434  next_break_forward =
1435  (parent.block_indices.local_to_global(current_block, 0) +
1436  parent.block_indices.block_size(current_block) - 1);
1437  }
1438  else
1439  // past the end. only have one
1440  // value for this
1441  {
1442  this->global_index = parent.size();
1443  current_block = parent.n_blocks();
1444  index_within_block = 0;
1445  next_break_backward = global_index;
1446  next_break_forward = numbers::invalid_size_type;
1447  };
1448  }
1449 
1450 
1451 
1452  template <class BlockVectorType, bool Constness>
1453  void
1454  Iterator<BlockVectorType, Constness>::move_forward()
1455  {
1456  if (global_index != next_break_forward)
1457  ++index_within_block;
1458  else
1459  {
1460  // ok, we traverse a boundary
1461  // between blocks:
1462  index_within_block = 0;
1463  ++current_block;
1464 
1465  // break backwards is now old
1466  // break forward
1467  next_break_backward = next_break_forward + 1;
1468 
1469  // compute new break forward
1470  if (current_block < parent->block_indices.size())
1471  next_break_forward +=
1472  parent->block_indices.block_size(current_block);
1473  else
1474  // if we are beyond the end,
1475  // then move the next
1476  // boundary arbitrarily far
1477  // away
1478  next_break_forward = numbers::invalid_size_type;
1479  };
1480 
1481  ++global_index;
1482  }
1483 
1484 
1485 
1486  template <class BlockVectorType, bool Constness>
1487  void
1488  Iterator<BlockVectorType, Constness>::move_backward()
1489  {
1490  if (global_index != next_break_backward)
1491  --index_within_block;
1492  else if (current_block != 0)
1493  {
1494  // ok, we traverse a boundary
1495  // between blocks:
1496  --current_block;
1497  index_within_block =
1498  parent->block_indices.block_size(current_block) - 1;
1499 
1500  // break forwards is now old
1501  // break backward
1502  next_break_forward = next_break_backward - 1;
1503 
1504  // compute new break forward
1505  next_break_backward -=
1506  parent->block_indices.block_size(current_block);
1507  }
1508  else
1509  // current block was 0, we now
1510  // get into unspecified terrain
1511  {
1512  --current_block;
1513  index_within_block = numbers::invalid_size_type;
1514  next_break_forward = 0;
1515  next_break_backward = 0;
1516  };
1517 
1518  --global_index;
1519  }
1520 
1521 
1522  } // namespace BlockVectorIterators
1523 
1524 } // namespace internal
1525 
1526 
1527 
1528 template <class VectorType>
1529 inline std::size_t
1531 {
1532  return block_indices.total_size();
1533 }
1534 
1535 
1536 
1537 template <class VectorType>
1538 inline IndexSet
1540 {
1541  IndexSet is(size());
1542 
1543  // copy index sets from blocks into the global one, shifted
1544  // by the appropriate amount for each block
1545  for (unsigned int b = 0; b < n_blocks(); ++b)
1546  {
1547  IndexSet x = block(b).locally_owned_elements();
1548  is.add_indices(x, block_indices.block_start(b));
1549  }
1550 
1551  is.compress();
1552 
1553  return is;
1554 }
1555 
1556 
1557 
1558 template <class VectorType>
1559 inline unsigned int
1561 {
1562  return block_indices.size();
1563 }
1564 
1565 
1566 template <class VectorType>
1568 BlockVectorBase<VectorType>::block(const unsigned int i)
1569 {
1570  Assert(i < n_blocks(), ExcIndexRange(i, 0, n_blocks()));
1571 
1572  return components[i];
1573 }
1574 
1575 
1576 
1577 template <class VectorType>
1578 inline const typename BlockVectorBase<VectorType>::BlockType &
1579 BlockVectorBase<VectorType>::block(const unsigned int i) const
1580 {
1581  Assert(i < n_blocks(), ExcIndexRange(i, 0, n_blocks()));
1582 
1583  return components[i];
1584 }
1585 
1586 
1587 
1588 template <class VectorType>
1589 inline const BlockIndices &
1591 {
1592  return block_indices;
1593 }
1594 
1595 
1596 template <class VectorType>
1597 inline void
1599 {
1600  std::vector<size_type> sizes(n_blocks());
1601 
1602  for (size_type i = 0; i < n_blocks(); ++i)
1603  sizes[i] = block(i).size();
1604 
1605  block_indices.reinit(sizes);
1606 }
1607 
1608 
1609 
1610 template <class VectorType>
1611 inline void
1613  ::VectorOperation::values operation)
1614 {
1615  for (unsigned int i = 0; i < n_blocks(); ++i)
1616  block(i).compress(operation);
1617 }
1618 
1619 
1620 
1621 template <class VectorType>
1624 {
1625  return iterator(*this, 0U);
1626 }
1627 
1628 
1629 
1630 template <class VectorType>
1633 {
1634  return const_iterator(*this, 0U);
1635 }
1636 
1637 
1638 template <class VectorType>
1641 {
1642  return iterator(*this, size());
1643 }
1644 
1645 
1646 
1647 template <class VectorType>
1650 {
1651  return const_iterator(*this, size());
1652 }
1653 
1654 
1655 template <class VectorType>
1656 inline bool
1657 BlockVectorBase<VectorType>::in_local_range(const size_type global_index) const
1658 {
1659  const std::pair<size_type, size_type> local_index =
1660  block_indices.global_to_local(global_index);
1661 
1662  return components[local_index.first].in_local_range(global_index);
1663 }
1664 
1665 
1666 template <class VectorType>
1667 bool
1669 {
1670  for (size_type i = 0; i < n_blocks(); ++i)
1671  if (components[i].all_zero() == false)
1672  return false;
1673 
1674  return true;
1675 }
1676 
1677 
1678 
1679 template <class VectorType>
1680 bool
1682 {
1683  for (size_type i = 0; i < n_blocks(); ++i)
1684  if (components[i].is_non_negative() == false)
1685  return false;
1686 
1687  return true;
1688 }
1689 
1690 
1691 
1692 template <class VectorType>
1693 typename BlockVectorBase<VectorType>::value_type BlockVectorBase<VectorType>::
1695 {
1696  Assert(n_blocks() == v.n_blocks(),
1698 
1699  value_type sum = 0.;
1700  for (size_type i = 0; i < n_blocks(); ++i)
1701  sum += components[i] * v.components[i];
1702 
1703  return sum;
1704 }
1705 
1706 
1707 template <class VectorType>
1710 {
1711  real_type sum = 0.;
1712  for (size_type i = 0; i < n_blocks(); ++i)
1713  sum += components[i].norm_sqr();
1714 
1715  return sum;
1716 }
1717 
1718 
1719 
1720 template <class VectorType>
1721 typename BlockVectorBase<VectorType>::value_type
1723 {
1724  value_type sum = 0.;
1725  // need to do static_cast as otherwise it won't work with
1726  // value_type=complex<T>
1727  for (size_type i = 0; i < n_blocks(); ++i)
1728  sum += components[i].mean_value() *
1730  components[i].size()));
1731 
1732  return sum / (typename numbers::NumberTraits<value_type>::real_type(size()));
1733 }
1734 
1735 
1736 
1737 template <class VectorType>
1740 {
1741  real_type sum = 0.;
1742  for (size_type i = 0; i < n_blocks(); ++i)
1743  sum += components[i].l1_norm();
1744 
1745  return sum;
1746 }
1747 
1748 
1749 
1750 template <class VectorType>
1753 {
1754  return std::sqrt(norm_sqr());
1755 }
1756 
1757 
1758 
1759 template <class VectorType>
1762 {
1763  real_type sum = 0.;
1764  for (size_type i = 0; i < n_blocks(); ++i)
1765  {
1766  value_type newval = components[i].linfty_norm();
1767  if (sum < newval)
1768  sum = newval;
1769  }
1770  return sum;
1771 }
1772 
1773 
1774 
1775 template <class VectorType>
1776 typename BlockVectorBase<VectorType>::value_type
1778  const typename BlockVectorBase<VectorType>::value_type a,
1779  const BlockVectorBase<VectorType> & V,
1780  const BlockVectorBase<VectorType> & W)
1781 {
1784 
1785  value_type sum = 0.;
1786  for (size_type i = 0; i < n_blocks(); ++i)
1787  sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1788 
1789  return sum;
1790 }
1791 
1792 
1793 
1794 template <class VectorType>
1797 {
1798  Assert(n_blocks() == v.n_blocks(),
1800 
1801  for (size_type i = 0; i < n_blocks(); ++i)
1802  {
1803  components[i] += v.components[i];
1804  }
1805 
1806  return *this;
1807 }
1808 
1809 
1810 
1811 template <class VectorType>
1814 {
1815  Assert(n_blocks() == v.n_blocks(),
1817 
1818  for (size_type i = 0; i < n_blocks(); ++i)
1819  {
1820  components[i] -= v.components[i];
1821  }
1822  return *this;
1823 }
1824 
1825 
1826 
1827 template <class VectorType>
1828 template <typename Number>
1829 inline void
1830 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1831  const std::vector<Number> & values)
1832 {
1833  Assert(indices.size() == values.size(),
1834  ExcDimensionMismatch(indices.size(), values.size()));
1835  add(indices.size(), indices.data(), values.data());
1836 }
1837 
1838 
1839 
1840 template <class VectorType>
1841 template <typename Number>
1842 inline void
1843 BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
1844  const Vector<Number> & values)
1845 {
1846  Assert(indices.size() == values.size(),
1847  ExcDimensionMismatch(indices.size(), values.size()));
1848  const size_type n_indices = indices.size();
1849  for (size_type i = 0; i < n_indices; ++i)
1850  (*this)(indices[i]) += values(i);
1851 }
1852 
1853 
1854 
1855 template <class VectorType>
1856 template <typename Number>
1857 inline void
1858 BlockVectorBase<VectorType>::add(const size_type n_indices,
1859  const size_type *indices,
1860  const Number * values)
1861 {
1862  for (size_type i = 0; i < n_indices; ++i)
1863  (*this)(indices[i]) += values[i];
1864 }
1865 
1866 
1867 
1868 template <class VectorType>
1869 void
1870 BlockVectorBase<VectorType>::add(const value_type a)
1871 {
1872  AssertIsFinite(a);
1873 
1874  for (size_type i = 0; i < n_blocks(); ++i)
1875  {
1876  components[i].add(a);
1877  }
1878 }
1879 
1880 
1881 
1882 template <class VectorType>
1883 void
1884 BlockVectorBase<VectorType>::add(const value_type a,
1885  const BlockVectorBase<VectorType> &v)
1886 {
1887  AssertIsFinite(a);
1888 
1889  Assert(n_blocks() == v.n_blocks(),
1891 
1892  for (size_type i = 0; i < n_blocks(); ++i)
1893  {
1894  components[i].add(a, v.components[i]);
1895  }
1896 }
1897 
1898 
1899 
1900 template <class VectorType>
1901 void
1902 BlockVectorBase<VectorType>::add(const value_type a,
1903  const BlockVectorBase<VectorType> &v,
1904  const value_type b,
1905  const BlockVectorBase<VectorType> &w)
1906 {
1907  AssertIsFinite(a);
1908  AssertIsFinite(b);
1909 
1910  Assert(n_blocks() == v.n_blocks(),
1912  Assert(n_blocks() == w.n_blocks(),
1914 
1915 
1916  for (size_type i = 0; i < n_blocks(); ++i)
1917  {
1918  components[i].add(a, v.components[i], b, w.components[i]);
1919  }
1920 }
1921 
1922 
1923 
1924 template <class VectorType>
1925 void
1926 BlockVectorBase<VectorType>::sadd(const value_type x,
1927  const BlockVectorBase<VectorType> &v)
1928 {
1929  AssertIsFinite(x);
1930 
1931  Assert(n_blocks() == v.n_blocks(),
1933 
1934  for (size_type i = 0; i < n_blocks(); ++i)
1935  {
1936  components[i].sadd(x, v.components[i]);
1937  }
1938 }
1939 
1940 
1941 
1942 template <class VectorType>
1943 void
1944 BlockVectorBase<VectorType>::sadd(const value_type x,
1945  const value_type a,
1946  const BlockVectorBase<VectorType> &v)
1947 {
1948  AssertIsFinite(x);
1949  AssertIsFinite(a);
1950 
1951  Assert(n_blocks() == v.n_blocks(),
1953 
1954  for (size_type i = 0; i < n_blocks(); ++i)
1955  {
1956  components[i].sadd(x, a, v.components[i]);
1957  }
1958 }
1959 
1960 
1961 
1962 template <class VectorType>
1963 void
1964 BlockVectorBase<VectorType>::sadd(const value_type x,
1965  const value_type a,
1966  const BlockVectorBase<VectorType> &v,
1967  const value_type b,
1968  const BlockVectorBase<VectorType> &w)
1969 {
1970  AssertIsFinite(x);
1971  AssertIsFinite(a);
1972  AssertIsFinite(b);
1973 
1974  Assert(n_blocks() == v.n_blocks(),
1976  Assert(n_blocks() == w.n_blocks(),
1978 
1979  for (size_type i = 0; i < n_blocks(); ++i)
1980  {
1981  components[i].sadd(x, a, v.components[i], b, w.components[i]);
1982  }
1983 }
1984 
1985 
1986 
1987 template <class VectorType>
1988 void
1989 BlockVectorBase<VectorType>::sadd(const value_type x,
1990  const value_type a,
1991  const BlockVectorBase<VectorType> &v,
1992  const value_type b,
1993  const BlockVectorBase<VectorType> &w,
1994  const value_type c,
1995  const BlockVectorBase<VectorType> &y)
1996 {
1997  AssertIsFinite(x);
1998  AssertIsFinite(a);
1999  AssertIsFinite(b);
2000  AssertIsFinite(c);
2001 
2002  Assert(n_blocks() == v.n_blocks(),
2004  Assert(n_blocks() == w.n_blocks(),
2006  Assert(n_blocks() == y.n_blocks(),
2008 
2009  for (size_type i = 0; i < n_blocks(); ++i)
2010  {
2011  components[i].sadd(
2012  x, a, v.components[i], b, w.components[i], c, y.components[i]);
2013  }
2014 }
2015 
2016 
2017 
2018 template <class VectorType>
2019 template <class BlockVector2>
2020 void
2021 BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
2022 {
2023  Assert(n_blocks() == v.n_blocks(),
2024  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
2025  for (size_type i = 0; i < n_blocks(); ++i)
2026  components[i].scale(v.block(i));
2027 }
2028 
2029 
2030 
2031 template <class VectorType>
2032 void
2033 BlockVectorBase<VectorType>::equ(const value_type a,
2034  const BlockVectorBase<VectorType> &v,
2035  const value_type b,
2036  const BlockVectorBase<VectorType> &w)
2037 {
2038  AssertIsFinite(a);
2039  AssertIsFinite(b);
2040 
2041  Assert(n_blocks() == v.n_blocks(),
2043  Assert(n_blocks() == w.n_blocks(),
2045 
2046  for (size_type i = 0; i < n_blocks(); ++i)
2047  {
2048  components[i].equ(a, v.components[i], b, w.components[i]);
2049  }
2050 }
2051 
2052 
2053 
2054 template <class VectorType>
2055 std::size_t
2057 {
2058  return (MemoryConsumption::memory_consumption(this->block_indices) +
2059  MemoryConsumption::memory_consumption(this->components));
2060 }
2061 
2062 
2063 
2064 template <class VectorType>
2065 template <class BlockVector2>
2066 void
2067 BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
2068 {
2069  AssertIsFinite(a);
2070 
2071  Assert(n_blocks() == v.n_blocks(),
2072  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
2073 
2074  for (size_type i = 0; i < n_blocks(); ++i)
2075  components[i].equ(a, v.components[i]);
2076 }
2077 
2078 
2079 
2080 template <class VectorType>
2081 void
2083 {
2084  for (size_type i = 0; i < n_blocks(); ++i)
2085  block(i).update_ghost_values();
2086 }
2087 
2088 
2089 
2090 template <class VectorType>
2092 BlockVectorBase<VectorType>::operator=(const value_type s)
2093 {
2094  AssertIsFinite(s);
2095 
2096  for (size_type i = 0; i < n_blocks(); ++i)
2097  components[i] = s;
2098 
2099  return *this;
2100 }
2101 
2102 
2103 template <class VectorType>
2106 {
2108 
2109  for (size_type i = 0; i < n_blocks(); ++i)
2110  components[i] = v.components[i];
2111 
2112  return *this;
2113 }
2114 
2115 
2116 template <class VectorType>
2117 template <class VectorType2>
2120 {
2122 
2123  for (size_type i = 0; i < n_blocks(); ++i)
2124  components[i] = v.components[i];
2125 
2126  return *this;
2127 }
2128 
2129 
2130 
2131 template <class VectorType>
2133 BlockVectorBase<VectorType>::operator=(const VectorType &v)
2134 {
2135  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2136 
2137  size_type index_v = 0;
2138  for (size_type b = 0; b < n_blocks(); ++b)
2139  for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
2140  block(b)(i) = v(index_v);
2141 
2142  return *this;
2143 }
2144 
2145 
2146 
2147 template <class VectorType>
2148 template <class VectorType2>
2149 inline bool
2152 {
2153  Assert(block_indices == v.block_indices, ExcDifferentBlockIndices());
2154 
2155  for (size_type i = 0; i < n_blocks(); ++i)
2156  if (!(components[i] == v.components[i]))
2157  return false;
2158 
2159  return true;
2160 }
2161 
2162 
2163 
2164 template <class VectorType>
2166 BlockVectorBase<VectorType>::operator*=(const value_type factor)
2167 {
2168  AssertIsFinite(factor);
2169 
2170  for (size_type i = 0; i < n_blocks(); ++i)
2171  components[i] *= factor;
2172 
2173  return *this;
2174 }
2175 
2176 
2177 
2178 template <class VectorType>
2180 BlockVectorBase<VectorType>::operator/=(const value_type factor)
2181 {
2182  AssertIsFinite(factor);
2183  Assert(factor != 0., ExcDivideByZero());
2184 
2185  for (size_type i = 0; i < n_blocks(); ++i)
2186  components[i] /= factor;
2187 
2188  return *this;
2189 }
2190 
2191 
2192 template <class VectorType>
2193 inline typename BlockVectorBase<VectorType>::value_type
2194 BlockVectorBase<VectorType>::operator()(const size_type i) const
2195 {
2196  const std::pair<unsigned int, size_type> local_index =
2197  block_indices.global_to_local(i);
2198  return components[local_index.first](local_index.second);
2199 }
2200 
2201 
2202 
2203 template <class VectorType>
2204 inline typename BlockVectorBase<VectorType>::reference
2205 BlockVectorBase<VectorType>::operator()(const size_type i)
2206 {
2207  const std::pair<unsigned int, size_type> local_index =
2208  block_indices.global_to_local(i);
2209  return components[local_index.first](local_index.second);
2210 }
2211 
2212 
2213 
2214 template <class VectorType>
2215 inline typename BlockVectorBase<VectorType>::value_type
2216  BlockVectorBase<VectorType>::operator[](const size_type i) const
2217 {
2218  return operator()(i);
2219 }
2220 
2221 
2222 
2223 template <class VectorType>
2224 inline typename BlockVectorBase<VectorType>::reference
2225  BlockVectorBase<VectorType>::operator[](const size_type i)
2226 {
2227  return operator()(i);
2228 }
2229 
2230 
2231 
2232 template <typename VectorType>
2233 template <typename OtherNumber>
2234 inline void
2236  const std::vector<size_type> &indices,
2237  std::vector<OtherNumber> & values) const
2238 {
2239  for (size_type i = 0; i < indices.size(); ++i)
2240  values[i] = operator()(indices[i]);
2241 }
2242 
2243 
2244 
2245 template <typename VectorType>
2246 template <typename ForwardIterator, typename OutputIterator>
2247 inline void
2249  ForwardIterator indices_begin,
2250  const ForwardIterator indices_end,
2251  OutputIterator values_begin) const
2252 {
2253  while (indices_begin != indices_end)
2254  {
2255  *values_begin = operator()(*indices_begin);
2256  indices_begin++;
2257  values_begin++;
2258  }
2259 }
2260 
2261 #endif // DOXYGEN
2262 
2263 DEAL_II_NAMESPACE_CLOSE
2264 
2265 #endif
const types::global_dof_index invalid_size_type
Definition: types.h:182
static yes_type check_for_block_vector(const BlockVectorBase< T > *)
std::size_t size() const
BaseClass::value_type value_type
Definition: block_vector.h:84
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1366
bool operator==(const BlockVectorBase< VectorType2 > &v) const
Types< BlockVectorType, Constness >::value_type value_type
void add(const std::vector< size_type > &indices, const std::vector< Number > &values)
difference_type operator-(const Iterator< BlockVectorType, OtherConstness > &i) const
void compress(::VectorOperation::values operation)
Types< BlockVectorType, Constness >::BlockVector BlockVector
bool operator==(const Iterator< BlockVectorType, OtherConstness > &i) const
real_type linfty_norm() const
real_type norm_sqr() const
void collect_sizes()
STL namespace.
unsigned int n_blocks() 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()
value_type operator[](const size_type i) const
LinearAlgebra::distributed::Vector< Number > Vector
BlockVectorBase & operator/=(const value_type factor)
value_type add_and_dot(const value_type a, const BlockVectorBase &V, const BlockVectorBase &W)
Iterator operator+(const difference_type &d) const
void scale(const BlockVector2 &v)
BlockVectorBase()=default
bool is_non_negative() const
dereference_type operator*() const
value_type operator()(const size_type i) const
BlockIndices block_indices
size_type block_start(const unsigned int i) const
bool operator>(const Iterator< BlockVectorType, OtherConstness > &i) const
bool in_local_range(const size_type global_index) const
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int global_dof_index
Definition: types.h:89
#define Assert(cond, exc)
Definition: exceptions.h:1227
void update_ghost_values() const
BlockType::real_type real_type
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const BlockIndices & get_block_indices() const
void reinit(const unsigned int n_blocks, const size_type n_elements_per_block)
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:397
bool operator>=(const Iterator< BlockVectorType, OtherConstness > &i) const
size_type total_size() const
bool operator!=(const Iterator< BlockVectorType, OtherConstness > &i) const
value_type mean_value() const
bool all_zero() const
static::ExceptionBase & ExcDifferentBlockIndices()
static const bool value
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
BlockVectorBase & operator+=(const BlockVectorBase &V)
std::vector< VectorType > components
size_type block_size(const unsigned int i) const
BlockVectorBase & operator*=(const value_type factor)
Iterator & operator+=(const difference_type &d)
void sadd(const value_type s, const BlockVectorBase &V)
BlockVectorBase & operator-=(const BlockVectorBase &V)
iterator end()
value_type operator*(const BlockVectorBase &V) const
real_type l2_norm() const
IndexSet locally_owned_elements() const
unsigned int size() const
BlockType & block(const unsigned int i)
size_type local_to_global(const unsigned int block, const size_type index) const
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)
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< OtherNumber > &values) const
static::ExceptionBase & ExcPointerToDifferentVectors()
#define AssertIsFinite(number)
Definition: exceptions.h:1423
real_type l1_norm() const
void equ(const value_type a, const BlockVector2 &V)
Iterator & operator=(const Iterator &c)
std::random_access_iterator_tag iterator_category
Iterator & operator-=(const difference_type &d)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
std::size_t memory_consumption() const
dereference_type operator[](const difference_type d) const