Reference documentation for deal.II version Git 1ad260a 2018-02-22 16:12:58 +0100
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 at
12 // the top level of the deal.II distribution.
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 #include <deal.II/base/numbers.h>
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/subscriptor.h>
24 #include <deal.II/base/memory_consumption.h>
25 #include <deal.II/lac/exceptions.h>
26 #include <deal.II/lac/block_indices.h>
27 #include <deal.II/lac/vector.h>
28 #include <deal.II/lac/vector_operation.h>
29 
30 #include <vector>
31 #include <iterator>
32 #include <cmath>
33 #include <cstddef>
34 
35 DEAL_II_NAMESPACE_OPEN
36 
37 
42 template <typename> class BlockVectorBase;
43 
44 
61 template <typename VectorType>
63 {
64 private:
65  struct yes_type
66  {
67  char c[1];
68  };
69  struct no_type
70  {
71  char c[2];
72  };
73 
78  template <typename T>
79  static yes_type check_for_block_vector (const BlockVectorBase<T> *);
80 
85  static no_type check_for_block_vector (...);
86 
87 public:
93  static const bool value = (sizeof(check_for_block_vector
94  ((VectorType *)nullptr))
95  ==
96  sizeof(yes_type));
97 };
98 
99 
100 // instantiation of the static member
101 template <typename VectorType>
103 
104 
105 
106 
107 namespace internal
108 {
109 
115  namespace BlockVectorIterators
116  {
121  template <class BlockVectorType, bool Constness>
122  struct Types
123  {
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  public std::iterator<std::random_access_iterator_tag,
219  typename Types<BlockVectorType,Constness>::value_type>
220  {
221  public:
226 
232  typedef
235 
241  typedef std::random_access_iterator_tag iterator_type;
242  typedef std::ptrdiff_t difference_type;
243  typedef typename BlockVectorType::reference reference;
244  typedef value_type *pointer;
245 
246  typedef
248  dereference_type;
249 
254  typedef
257 
267  const size_type global_index);
268 
277 
278 
282  Iterator (const Iterator &c);
283 
284  private:
290  const size_type global_index,
291  const size_type current_block,
292  const size_type index_within_block,
293  const size_type next_break_forward,
294  const size_type next_break_backward);
295 
296  public:
297 
301  Iterator &operator = (const Iterator &c);
302 
308  dereference_type operator * () const;
309 
314  dereference_type operator [] (const difference_type d) const;
315 
321 
327  Iterator operator ++ (int);
328 
334 
340  Iterator operator -- (int);
341 
346  template <bool OtherConstness>
348 
353  template <bool OtherConstness>
355 
361  template <bool OtherConstness>
362  bool operator < (const Iterator<BlockVectorType, OtherConstness> &i) const;
363 
367  template <bool OtherConstness>
368  bool operator <= (const Iterator<BlockVectorType, OtherConstness> &i) const;
369 
373  template <bool OtherConstness>
375 
379  template <bool OtherConstness>
381 
385  template <bool OtherConstness>
386  difference_type operator - (const Iterator<BlockVectorType, OtherConstness> &i) const;
387 
392  Iterator operator + (const difference_type &d) const;
393 
398  Iterator operator - (const difference_type &d) const;
399 
404  Iterator &operator += (const difference_type &d);
405 
410  Iterator &operator -= (const difference_type &d);
411 
422  "Your program tried to compare iterators pointing to "
423  "different block vectors. There is no reasonable way "
424  "to do this.");
425 
427  private:
434 
438  size_type global_index;
439 
444  unsigned int current_block;
445  size_type index_within_block;
446 
454  size_type next_break_backward;
455 
459  void move_forward ();
460 
464  void move_backward ();
465 
466 
470  template <typename, bool>
471  friend class Iterator;
472  };
473  } // namespace BlockVectorIterators
474 } // namespace internal
475 
476 
515 template <class VectorType>
516 class BlockVectorBase : public Subscriptor
517 {
518 public:
522  typedef VectorType BlockType;
523 
524  /*
525  * Declare standard types used in
526  * all containers. These types
527  * parallel those in the
528  * <tt>C++</tt> standard
529  * libraries
530  * <tt>std::vector<...></tt>
531  * class. This includes iterator
532  * types.
533  */
534  typedef typename BlockType::value_type value_type;
535  typedef value_type *pointer;
536  typedef const value_type *const_pointer;
537  typedef ::internal::BlockVectorIterators::Iterator<BlockVectorBase,false> iterator;
538  typedef ::internal::BlockVectorIterators::Iterator<BlockVectorBase,true> const_iterator;
539  typedef typename BlockType::reference reference;
540  typedef typename BlockType::const_reference const_reference;
541  typedef types::global_dof_index size_type;
542 
552  typedef typename BlockType::real_type real_type;
553 
557  BlockVectorBase () = default;
558 
562  BlockVectorBase (const BlockVectorBase &/*V*/) = default;
563 
569  BlockVectorBase (BlockVectorBase &&/*V*/) = default;
570 
577  void collect_sizes ();
578 
589  void compress (::VectorOperation::values operation);
590 
594  BlockType &
595  block (const unsigned int i);
596 
600  const BlockType &
601  block (const unsigned int i) const;
602 
608  const BlockIndices &
609  get_block_indices () const;
610 
614  unsigned int n_blocks () const;
615 
620  std::size_t size () const;
621 
639 
643  iterator begin ();
644 
649  const_iterator begin () const;
650 
654  iterator end ();
655 
660  const_iterator end () const;
661 
665  value_type operator() (const size_type i) const;
666 
670  reference operator() (const size_type i);
671 
677  value_type operator[] (const size_type i) const;
678 
684  reference operator[] (const size_type i);
685 
701  template <typename OtherNumber>
702  void extract_subvector_to (const std::vector<size_type> &indices,
703  std::vector<OtherNumber> &values) const;
704 
732  template <typename ForwardIterator, typename OutputIterator>
733  void extract_subvector_to (ForwardIterator indices_begin,
734  const ForwardIterator indices_end,
735  OutputIterator values_begin) const;
736 
741  BlockVectorBase &operator = (const value_type s);
742 
747  operator= (const BlockVectorBase &V);
748 
754  BlockVectorBase &operator= (BlockVectorBase &&/*V*/) = default;
755 
759  template <class VectorType2>
762 
767  operator = (const VectorType &v);
768 
773  template <class VectorType2>
774  bool
776 
780  value_type operator* (const BlockVectorBase &V) const;
781 
785  real_type norm_sqr () const;
786 
790  value_type mean_value () const;
791 
795  real_type l1_norm () const;
796 
801  real_type l2_norm () const;
802 
807  real_type linfty_norm () const;
808 
830  value_type add_and_dot (const value_type a,
831  const BlockVectorBase &V,
832  const BlockVectorBase &W);
833 
838  bool in_local_range (const size_type global_index) const;
839 
845  bool all_zero () const;
846 
852  bool is_non_negative () const;
853 
858  operator += (const BlockVectorBase &V);
859 
864  operator -= (const BlockVectorBase &V);
865 
866 
871  template <typename Number>
872  void add (const std::vector<size_type> &indices,
873  const std::vector<Number> &values);
874 
879  template <typename Number>
880  void add (const std::vector<size_type> &indices,
881  const Vector<Number> &values);
882 
888  template <typename Number>
889  void add (const size_type n_elements,
890  const size_type *indices,
891  const Number *values);
892 
897  void add (const value_type s);
898 
902  void add (const value_type a, const BlockVectorBase &V);
903 
907  void add (const value_type a, const BlockVectorBase &V,
908  const value_type b, const BlockVectorBase &W);
909 
913  void sadd (const value_type s, const BlockVectorBase &V);
914 
918  void sadd (const value_type s, const value_type a, const BlockVectorBase &V);
919 
923  void sadd (const value_type s, const value_type a,
924  const BlockVectorBase &V,
925  const value_type b, const BlockVectorBase &W);
926 
930  void sadd (const value_type s, const value_type a,
931  const BlockVectorBase &V,
932  const value_type b, const BlockVectorBase &W,
933  const value_type c, const BlockVectorBase &X);
934 
938  BlockVectorBase &operator *= (const value_type factor);
939 
943  BlockVectorBase &operator /= (const value_type factor);
944 
949  template <class BlockVector2>
950  void scale (const BlockVector2 &v);
951 
955  template <class BlockVector2>
956  void equ (const value_type a, const BlockVector2 &V);
957 
961  void equ (const value_type a, const BlockVectorBase &V,
962  const value_type b, const BlockVectorBase &W);
963 
968  void update_ghost_values () const;
969 
974  std::size_t memory_consumption () const;
975 
976 protected:
980  std::vector<VectorType> components;
981 
987 
991  template <typename N, bool C>
992  friend class ::internal::BlockVectorIterators::Iterator;
993 
994  template <typename> friend class BlockVectorBase;
995 };
996 
997 
1000 /*----------------------- Inline functions ----------------------------------*/
1001 
1002 
1003 #ifndef DOXYGEN
1004 namespace internal
1005 {
1006  namespace BlockVectorIterators
1007  {
1008  template <class BlockVectorType, bool Constness>
1009  inline
1011  Iterator (const Iterator<BlockVectorType,Constness> &c)
1012  :
1013  parent (c.parent),
1014  global_index (c.global_index),
1015  current_block (c.current_block),
1016  index_within_block (c.index_within_block),
1017  next_break_forward (c.next_break_forward),
1018  next_break_backward (c.next_break_backward)
1019  {}
1020 
1021 
1022 
1023  template <class BlockVectorType, bool Constness>
1024  inline
1025  Iterator<BlockVectorType,Constness>::
1026  Iterator (const Iterator<BlockVectorType,!Constness> &c)
1027  :
1028  parent (c.parent),
1030  current_block (c.current_block),
1031  index_within_block (c.index_within_block),
1032  next_break_forward (c.next_break_forward),
1033  next_break_backward (c.next_break_backward)
1034  {
1035  // Only permit copy-constructing const iterators from non-const
1036  // iterators, and not vice versa (i.e., Constness must always be
1037  // true). As mentioned above, try to check this at compile time if C++11
1038  // support is available.
1039  static_assert(Constness == true,
1040  "Constructing a non-const iterator from a const iterator "
1041  "does not make sense.");
1042  }
1043 
1044 
1045 
1046  template <class BlockVectorType, bool Constness>
1047  inline
1048  Iterator<BlockVectorType,Constness>::
1049  Iterator (BlockVector &parent,
1050  const size_type global_index,
1051  const size_type current_block,
1052  const size_type index_within_block,
1053  const size_type next_break_forward,
1054  const size_type next_break_backward)
1055  :
1056  parent (&parent),
1057  global_index (global_index),
1058  current_block (current_block),
1059  index_within_block (index_within_block),
1060  next_break_forward (next_break_forward),
1061  next_break_backward (next_break_backward)
1062  {
1063  }
1064 
1065 
1066 
1067  template <class BlockVectorType, bool Constness>
1068  inline
1069  Iterator<BlockVectorType,Constness> &
1070  Iterator<BlockVectorType,Constness>::
1071  operator = (const Iterator &c)
1072  {
1073  parent = c.parent;
1074  global_index = c.global_index;
1075  index_within_block = c.index_within_block;
1076  current_block = c.current_block;
1077  next_break_forward = c.next_break_forward;
1078  next_break_backward = c.next_break_backward;
1079 
1080  return *this;
1081  }
1082 
1083 
1084 
1085  template <class BlockVectorType, bool Constness>
1086  inline
1087  typename Iterator<BlockVectorType,Constness>::dereference_type
1088  Iterator<BlockVectorType,Constness>::operator * () const
1089  {
1090  return parent->block(current_block)(index_within_block);
1091  }
1092 
1093 
1094 
1095  template <class BlockVectorType, bool Constness>
1096  inline
1097  typename Iterator<BlockVectorType,Constness>::dereference_type
1098  Iterator<BlockVectorType,Constness>::operator [] (const difference_type d) const
1099  {
1100  // if the index pointed to is
1101  // still within the block we
1102  // currently point into, then we
1103  // can save the computation of
1104  // the block
1105  if ((global_index+d >= next_break_backward) &&
1106  (global_index+d <= next_break_forward))
1107  return parent->block(current_block)(index_within_block + d);
1108 
1109  // if the index is not within the
1110  // block of the block vector into
1111  // which we presently point, then
1112  // there is no way: we have to
1113  // search for the block. this can
1114  // be done through the parent
1115  // class as well.
1116  return (*parent)(global_index+d);
1117  }
1118 
1119 
1120 
1121  template <class BlockVectorType, bool Constness>
1122  inline
1123  Iterator<BlockVectorType,Constness> &
1124  Iterator<BlockVectorType,Constness>::operator ++ ()
1125  {
1126  move_forward ();
1127  return *this;
1128  }
1129 
1130 
1131 
1132  template <class BlockVectorType, bool Constness>
1133  inline
1134  Iterator<BlockVectorType,Constness>
1135  Iterator<BlockVectorType,Constness>::operator ++ (int)
1136  {
1137  const Iterator old_value = *this;
1138  move_forward ();
1139  return old_value;
1140  }
1141 
1142 
1143 
1144  template <class BlockVectorType, bool Constness>
1145  inline
1146  Iterator<BlockVectorType,Constness> &
1147  Iterator<BlockVectorType,Constness>::operator -- ()
1148  {
1149  move_backward ();
1150  return *this;
1151  }
1152 
1153 
1154 
1155  template <class BlockVectorType, bool Constness>
1156  inline
1157  Iterator<BlockVectorType,Constness>
1158  Iterator<BlockVectorType,Constness>::operator -- (int)
1159  {
1160  const Iterator old_value = *this;
1161  move_backward ();
1162  return old_value;
1163  }
1164 
1165 
1166 
1167  template <class BlockVectorType, bool Constness>
1168  template <bool OtherConstness>
1169  inline
1170  bool
1171  Iterator<BlockVectorType,Constness>::
1172  operator == (const Iterator<BlockVectorType, OtherConstness> &i) const
1173  {
1174  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1175 
1176  return (global_index == i.global_index);
1177  }
1178 
1179 
1180 
1181  template <class BlockVectorType, bool Constness>
1182  template <bool OtherConstness>
1183  inline
1184  bool
1185  Iterator<BlockVectorType,Constness>::
1186  operator != (const Iterator<BlockVectorType, OtherConstness> &i) const
1187  {
1188  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1189 
1190  return (global_index != i.global_index);
1191  }
1192 
1193 
1194 
1195  template <class BlockVectorType, bool Constness>
1196  template <bool OtherConstness>
1197  inline
1198  bool
1199  Iterator<BlockVectorType,Constness>::
1200  operator < (const Iterator<BlockVectorType, OtherConstness> &i) const
1201  {
1202  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1203 
1204  return (global_index < i.global_index);
1205  }
1206 
1207 
1208 
1209  template <class BlockVectorType, bool Constness>
1210  template <bool OtherConstness>
1211  inline
1212  bool
1213  Iterator<BlockVectorType,Constness>::
1214  operator <= (const Iterator<BlockVectorType, OtherConstness> &i) const
1215  {
1216  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1217 
1218  return (global_index <= i.global_index);
1219  }
1220 
1221 
1222 
1223  template <class BlockVectorType, bool Constness>
1224  template <bool OtherConstness>
1225  inline
1226  bool
1227  Iterator<BlockVectorType,Constness>::
1228  operator > (const Iterator<BlockVectorType, OtherConstness> &i) const
1229  {
1230  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1231 
1232  return (global_index > i.global_index);
1233  }
1234 
1235 
1236 
1237  template <class BlockVectorType, bool Constness>
1238  template <bool OtherConstness>
1239  inline
1240  bool
1241  Iterator<BlockVectorType,Constness>::
1242  operator >= (const Iterator<BlockVectorType, OtherConstness> &i) const
1243  {
1244  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1245 
1246  return (global_index >= i.global_index);
1247  }
1248 
1249 
1250 
1251  template <class BlockVectorType, bool Constness>
1252  template <bool OtherConstness>
1253  inline
1254  typename Iterator<BlockVectorType,Constness>::difference_type
1255  Iterator<BlockVectorType,Constness>::
1256  operator - (const Iterator<BlockVectorType, OtherConstness> &i) const
1257  {
1258  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1259 
1260  return (static_cast<signed int>(global_index) -
1261  static_cast<signed int>(i.global_index));
1262  }
1263 
1264 
1265 
1266  template <class BlockVectorType, bool Constness>
1267  inline
1268  Iterator<BlockVectorType,Constness>
1269  Iterator<BlockVectorType,Constness>::
1270  operator + (const difference_type &d) const
1271  {
1272  // if the index pointed to is
1273  // still within the block we
1274  // currently point into, then we
1275  // can save the computation of
1276  // the block
1277  if ((global_index+d >= next_break_backward) &&
1278  (global_index+d <= next_break_forward))
1279  return Iterator (*parent, global_index+d, current_block,
1280  index_within_block+d,
1281  next_break_forward, next_break_backward);
1282  else
1283  // outside present block, so
1284  // have to seek new block
1285  // anyway
1286  return Iterator (*parent, global_index+d);
1287  }
1288 
1289 
1290 
1291  template <class BlockVectorType, bool Constness>
1292  inline
1293  Iterator<BlockVectorType,Constness>
1294  Iterator<BlockVectorType,Constness>::
1295  operator - (const difference_type &d) const
1296  {
1297  // if the index pointed to is
1298  // still within the block we
1299  // currently point into, then we
1300  // can save the computation of
1301  // the block
1302  if ((global_index-d >= next_break_backward) &&
1303  (global_index-d <= next_break_forward))
1304  return Iterator (*parent, global_index-d, current_block,
1305  index_within_block-d,
1306  next_break_forward, next_break_backward);
1307  else
1308  // outside present block, so
1309  // have to seek new block
1310  // anyway
1311  return Iterator (*parent, global_index-d);
1312  }
1313 
1314 
1315 
1316  template <class BlockVectorType, bool Constness>
1317  inline
1318  Iterator<BlockVectorType,Constness> &
1319  Iterator<BlockVectorType,Constness>::
1320  operator += (const difference_type &d)
1321  {
1322  // if the index pointed to is
1323  // still within the block we
1324  // currently point into, then we
1325  // can save the computation of
1326  // the block
1327  if ((global_index+d >= next_break_backward) &&
1328  (global_index+d <= next_break_forward))
1329  {
1330  global_index += d;
1331  index_within_block += d;
1332  }
1333  else
1334  // outside present block, so
1335  // have to seek new block
1336  // anyway
1337  *this = Iterator (*parent, global_index+d);
1338 
1339  return *this;
1340  }
1341 
1342 
1343 
1344  template <class BlockVectorType, bool Constness>
1345  inline
1346  Iterator<BlockVectorType,Constness> &
1347  Iterator<BlockVectorType,Constness>::
1348  operator -= (const difference_type &d)
1349  {
1350  // if the index pointed to is
1351  // still within the block we
1352  // currently point into, then we
1353  // can save the computation of
1354  // the block
1355  if ((global_index-d >= next_break_backward) &&
1356  (global_index-d <= next_break_forward))
1357  {
1358  global_index -= d;
1359  index_within_block -= d;
1360  }
1361  else
1362  // outside present block, so
1363  // have to seek new block
1364  // anyway
1365  *this = Iterator (*parent, global_index-d);
1366 
1367  return *this;
1368  }
1369 
1370 
1371  template <class BlockVectorType, bool Constness>
1372  Iterator<BlockVectorType,Constness>::
1373  Iterator (BlockVector &parent,
1374  const size_type global_index)
1375  :
1376  parent (&parent),
1377  global_index (global_index)
1378  {
1379  // find which block we are
1380  // in. for this, take into
1381  // account that it happens at
1382  // times that people want to
1383  // initialize iterators
1384  // past-the-end
1385  if (global_index < parent.size())
1386  {
1387  const std::pair<size_type, size_type>
1388  indices = parent.block_indices.global_to_local(global_index);
1389  current_block = indices.first;
1390  index_within_block = indices.second;
1391 
1392  next_break_backward
1393  = parent.block_indices.local_to_global (current_block, 0);
1394  next_break_forward
1395  = (parent.block_indices.local_to_global (current_block, 0)
1396  +parent.block_indices.block_size(current_block)-1);
1397  }
1398  else
1399  // past the end. only have one
1400  // value for this
1401  {
1402  this->global_index = parent.size ();
1403  current_block = parent.n_blocks();
1404  index_within_block = 0;
1405  next_break_backward = global_index;
1406  next_break_forward = numbers::invalid_size_type;
1407  };
1408  }
1409 
1410 
1411 
1412  template <class BlockVectorType, bool Constness>
1413  void
1414  Iterator<BlockVectorType,Constness>::move_forward ()
1415  {
1416  if (global_index != next_break_forward)
1417  ++index_within_block;
1418  else
1419  {
1420  // ok, we traverse a boundary
1421  // between blocks:
1422  index_within_block = 0;
1423  ++current_block;
1424 
1425  // break backwards is now old
1426  // break forward
1427  next_break_backward = next_break_forward+1;
1428 
1429  // compute new break forward
1430  if (current_block < parent->block_indices.size())
1431  next_break_forward
1432  += parent->block_indices.block_size(current_block);
1433  else
1434  // if we are beyond the end,
1435  // then move the next
1436  // boundary arbitrarily far
1437  // away
1438  next_break_forward = numbers::invalid_size_type;
1439  };
1440 
1441  ++global_index;
1442  }
1443 
1444 
1445 
1446  template <class BlockVectorType, bool Constness>
1447  void
1448  Iterator<BlockVectorType,Constness>::move_backward ()
1449  {
1450  if (global_index != next_break_backward)
1451  --index_within_block;
1452  else if (current_block != 0)
1453  {
1454  // ok, we traverse a boundary
1455  // between blocks:
1456  --current_block;
1457  index_within_block = parent->block_indices.block_size(current_block)-1;
1458 
1459  // break forwards is now old
1460  // break backward
1461  next_break_forward = next_break_backward-1;
1462 
1463  // compute new break forward
1464  next_break_backward
1465  -= parent->block_indices.block_size (current_block);
1466  }
1467  else
1468  // current block was 0, we now
1469  // get into unspecified terrain
1470  {
1471  --current_block;
1472  index_within_block = numbers::invalid_size_type;
1473  next_break_forward = 0;
1474  next_break_backward = 0;
1475  };
1476 
1477  --global_index;
1478  }
1479 
1480 
1481  } // namespace BlockVectorIterators
1482 
1483 } //namespace internal
1484 
1485 
1486 
1487 template <class VectorType>
1488 inline
1489 std::size_t
1491 {
1492  return block_indices.total_size();
1493 }
1494 
1495 
1496 
1497 template <class VectorType>
1498 inline
1499 IndexSet
1501 {
1502  IndexSet is (size());
1503 
1504  // copy index sets from blocks into the global one, shifted
1505  // by the appropriate amount for each block
1506  for (unsigned int b=0; b<n_blocks(); ++b)
1507  {
1508  IndexSet x = block(b).locally_owned_elements();
1509  is.add_indices(x, block_indices.block_start(b));
1510  }
1511 
1512  is.compress();
1513 
1514  return is;
1515 }
1516 
1517 
1518 
1519 template <class VectorType>
1520 inline
1521 unsigned int
1523 {
1524  return block_indices.size();
1525 }
1526 
1527 
1528 template <class VectorType>
1529 inline
1531 BlockVectorBase<VectorType>::block (const unsigned int i)
1532 {
1533  Assert(i<n_blocks(), ExcIndexRange(i,0,n_blocks()));
1534 
1535  return components[i];
1536 }
1537 
1538 
1539 
1540 template <class VectorType>
1541 inline
1543 BlockVectorBase<VectorType>::block (const unsigned int i) const
1544 {
1545  Assert(i<n_blocks(), ExcIndexRange(i,0,n_blocks()));
1546 
1547  return components[i];
1548 }
1549 
1550 
1551 
1552 template <class VectorType>
1553 inline
1554 const BlockIndices &
1556 {
1557  return block_indices;
1558 }
1559 
1560 
1561 template <class VectorType>
1562 inline
1563 void
1565 {
1566  std::vector<size_type> sizes (n_blocks());
1567 
1568  for (size_type i=0; i<n_blocks(); ++i)
1569  sizes[i] = block(i).size();
1570 
1571  block_indices.reinit(sizes);
1572 }
1573 
1574 
1575 
1576 template <class VectorType>
1577 inline
1578 void
1580 {
1581  for (unsigned int i=0; i<n_blocks(); ++i)
1582  block(i).compress (operation);
1583 }
1584 
1585 
1586 
1587 template <class VectorType>
1588 inline
1591 {
1592  return iterator(*this, 0U);
1593 }
1594 
1595 
1596 
1597 template <class VectorType>
1598 inline
1601 {
1602  return const_iterator(*this, 0U);
1603 }
1604 
1605 
1606 template <class VectorType>
1607 inline
1610 {
1611  return iterator(*this, size());
1612 }
1613 
1614 
1615 
1616 template <class VectorType>
1617 inline
1620 {
1621  return const_iterator(*this, size());
1622 }
1623 
1624 
1625 template <class VectorType>
1626 inline
1627 bool
1629 (const size_type global_index) const
1630 {
1631  const std::pair<size_type,size_type> local_index
1632  = block_indices.global_to_local (global_index);
1633 
1634  return components[local_index.first].in_local_range (global_index);
1635 }
1636 
1637 
1638 template <class VectorType>
1639 bool
1641 {
1642  for (size_type i=0; i<n_blocks(); ++i)
1643  if (components[i].all_zero() == false)
1644  return false;
1645 
1646  return true;
1647 }
1648 
1649 
1650 
1651 template <class VectorType>
1652 bool
1654 {
1655  for (size_type i=0; i<n_blocks(); ++i)
1656  if (components[i].is_non_negative() == false)
1657  return false;
1658 
1659  return true;
1660 }
1661 
1662 
1663 
1664 template <class VectorType>
1665 typename BlockVectorBase<VectorType>::value_type
1668 {
1669  Assert (n_blocks() == v.n_blocks(),
1670  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1671 
1672  value_type sum = 0.;
1673  for (size_type i=0; i<n_blocks(); ++i)
1674  sum += components[i]*v.components[i];
1675 
1676  return sum;
1677 }
1678 
1679 
1680 template <class VectorType>
1683 {
1684  real_type sum = 0.;
1685  for (size_type i=0; i<n_blocks(); ++i)
1686  sum += components[i].norm_sqr();
1687 
1688  return sum;
1689 }
1690 
1691 
1692 
1693 template <class VectorType>
1694 typename BlockVectorBase<VectorType>::value_type
1696 {
1697  value_type sum = 0.;
1698  // need to do static_cast as otherwise it won't work with value_type=complex<T>
1699  for (size_type i=0; i<n_blocks(); ++i)
1700  sum += components[i].mean_value() * (typename numbers::NumberTraits<value_type>::real_type(components[i].size()));
1701 
1702  return sum/(typename numbers::NumberTraits<value_type>::real_type(size()));
1703 }
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].l1_norm();
1714 
1715  return sum;
1716 }
1717 
1718 
1719 
1720 template <class VectorType>
1723 {
1724  return std::sqrt(norm_sqr());
1725 }
1726 
1727 
1728 
1729 template <class VectorType>
1732 {
1733  real_type sum = 0.;
1734  for (size_type i=0; i<n_blocks(); ++i)
1735  {
1736  value_type newval = components[i].linfty_norm();
1737  if (sum<newval)
1738  sum = newval;
1739  }
1740  return sum;
1741 }
1742 
1743 
1744 
1745 template <class VectorType>
1746 typename BlockVectorBase<VectorType>::value_type
1748 add_and_dot (const typename BlockVectorBase<VectorType>::value_type a,
1749  const BlockVectorBase<VectorType> &V,
1750  const BlockVectorBase<VectorType> &W)
1751 {
1752  AssertDimension (n_blocks(), V.n_blocks());
1753  AssertDimension (n_blocks(), W.n_blocks());
1754 
1755  value_type sum = 0.;
1756  for (size_type i=0; i<n_blocks(); ++i)
1757  sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1758 
1759  return sum;
1760 }
1761 
1762 
1763 
1764 template <class VectorType>
1767 {
1768  Assert (n_blocks() == v.n_blocks(),
1769  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1770 
1771  for (size_type i=0; i<n_blocks(); ++i)
1772  {
1773  components[i] += v.components[i];
1774  }
1775 
1776  return *this;
1777 }
1778 
1779 
1780 
1781 template <class VectorType>
1784 {
1785  Assert (n_blocks() == v.n_blocks(),
1786  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1787 
1788  for (size_type i=0; i<n_blocks(); ++i)
1789  {
1790  components[i] -= v.components[i];
1791  }
1792  return *this;
1793 }
1794 
1795 
1796 
1797 template <class VectorType>
1798 template <typename Number>
1799 inline
1800 void
1801 BlockVectorBase<VectorType>::add (const std::vector<size_type> &indices,
1802  const std::vector<Number> &values)
1803 {
1804  Assert (indices.size() == values.size(),
1805  ExcDimensionMismatch(indices.size(), values.size()));
1806  add (indices.size(), indices.data(), values.data());
1807 }
1808 
1809 
1810 
1811 template <class VectorType>
1812 template <typename Number>
1813 inline
1814 void
1815 BlockVectorBase<VectorType>::add (const std::vector<size_type> &indices,
1816  const Vector<Number> &values)
1817 {
1818  Assert (indices.size() == values.size(),
1819  ExcDimensionMismatch(indices.size(), values.size()));
1820  const size_type n_indices = indices.size();
1821  for (size_type i=0; i<n_indices; ++i)
1822  (*this)(indices[i]) += values(i);
1823 }
1824 
1825 
1826 
1827 template <class VectorType>
1828 template <typename Number>
1829 inline
1830 void
1831 BlockVectorBase<VectorType>::add (const size_type n_indices,
1832  const size_type *indices,
1833  const Number *values)
1834 {
1835  for (size_type i=0; i<n_indices; ++i)
1836  (*this)(indices[i]) += values[i];
1837 }
1838 
1839 
1840 
1841 template <class VectorType>
1842 void BlockVectorBase<VectorType>::add (const value_type a)
1843 {
1844  AssertIsFinite(a);
1845 
1846  for (size_type i=0; i<n_blocks(); ++i)
1847  {
1848  components[i].add(a);
1849  }
1850 }
1851 
1852 
1853 
1854 template <class VectorType>
1855 void BlockVectorBase<VectorType>::add (const value_type a,
1856  const BlockVectorBase<VectorType> &v)
1857 {
1858 
1859  AssertIsFinite(a);
1860 
1861  Assert (n_blocks() == v.n_blocks(),
1862  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1863 
1864  for (size_type i=0; i<n_blocks(); ++i)
1865  {
1866  components[i].add(a, v.components[i]);
1867  }
1868 }
1869 
1870 
1871 
1872 template <class VectorType>
1873 void BlockVectorBase<VectorType>::add (const value_type a,
1874  const BlockVectorBase<VectorType> &v,
1875  const value_type b,
1876  const BlockVectorBase<VectorType> &w)
1877 {
1878 
1879  AssertIsFinite(a);
1880  AssertIsFinite(b);
1881 
1882  Assert (n_blocks() == v.n_blocks(),
1883  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1884  Assert (n_blocks() == w.n_blocks(),
1885  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1886 
1887 
1888  for (size_type i=0; i<n_blocks(); ++i)
1889  {
1890  components[i].add(a, v.components[i], b, w.components[i]);
1891  }
1892 }
1893 
1894 
1895 
1896 template <class VectorType>
1897 void BlockVectorBase<VectorType>::sadd (const value_type x,
1898  const BlockVectorBase<VectorType> &v)
1899 {
1900 
1901  AssertIsFinite(x);
1902 
1903  Assert (n_blocks() == v.n_blocks(),
1904  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1905 
1906  for (size_type i=0; i<n_blocks(); ++i)
1907  {
1908  components[i].sadd(x, v.components[i]);
1909  }
1910 }
1911 
1912 
1913 
1914 template <class VectorType>
1915 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
1916  const BlockVectorBase<VectorType> &v)
1917 {
1918 
1919  AssertIsFinite(x);
1920  AssertIsFinite(a);
1921 
1922  Assert (n_blocks() == v.n_blocks(),
1923  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1924 
1925  for (size_type i=0; i<n_blocks(); ++i)
1926  {
1927  components[i].sadd(x, a, v.components[i]);
1928  }
1929 }
1930 
1931 
1932 
1933 template <class VectorType>
1934 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
1935  const BlockVectorBase<VectorType> &v,
1936  const value_type b,
1937  const BlockVectorBase<VectorType> &w)
1938 {
1939 
1940  AssertIsFinite(x);
1941  AssertIsFinite(a);
1942  AssertIsFinite(b);
1943 
1944  Assert (n_blocks() == v.n_blocks(),
1945  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1946  Assert (n_blocks() == w.n_blocks(),
1947  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1948 
1949  for (size_type i=0; i<n_blocks(); ++i)
1950  {
1951  components[i].sadd(x, a, v.components[i], b, w.components[i]);
1952  }
1953 }
1954 
1955 
1956 
1957 template <class VectorType>
1958 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
1959  const BlockVectorBase<VectorType> &v,
1960  const value_type b,
1961  const BlockVectorBase<VectorType> &w,
1962  const value_type c,
1963  const BlockVectorBase<VectorType> &y)
1964 {
1965 
1966  AssertIsFinite(x);
1967  AssertIsFinite(a);
1968  AssertIsFinite(b);
1969  AssertIsFinite(c);
1970 
1971  Assert (n_blocks() == v.n_blocks(),
1972  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1973  Assert (n_blocks() == w.n_blocks(),
1974  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1975  Assert (n_blocks() == y.n_blocks(),
1976  ExcDimensionMismatch(n_blocks(), y.n_blocks()));
1977 
1978  for (size_type i=0; i<n_blocks(); ++i)
1979  {
1980  components[i].sadd(x, a, v.components[i],
1981  b, w.components[i], c, y.components[i]);
1982  }
1983 }
1984 
1985 
1986 
1987 template <class VectorType>
1988 template <class BlockVector2>
1989 void BlockVectorBase<VectorType>::scale (const BlockVector2 &v)
1990 {
1991  Assert (n_blocks() == v.n_blocks(),
1992  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1993  for (size_type i=0; i<n_blocks(); ++i)
1994  components[i].scale(v.block(i));
1995 }
1996 
1997 
1998 
1999 template <class VectorType>
2000 void BlockVectorBase<VectorType>::equ (const value_type a,
2001  const BlockVectorBase<VectorType> &v,
2002  const value_type b,
2003  const BlockVectorBase<VectorType> &w)
2004 {
2005 
2006  AssertIsFinite(a);
2007  AssertIsFinite(b);
2008 
2009  Assert (n_blocks() == v.n_blocks(),
2010  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
2011  Assert (n_blocks() == w.n_blocks(),
2012  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
2013 
2014  for (size_type i=0; i<n_blocks(); ++i)
2015  {
2016  components[i].equ( a, v.components[i], b, w.components[i]);
2017  }
2018 }
2019 
2020 
2021 
2022 template <class VectorType>
2023 std::size_t
2025 {
2026  return (MemoryConsumption::memory_consumption (this->block_indices)
2027  +
2028  MemoryConsumption::memory_consumption (this->components));
2029 }
2030 
2031 
2032 
2033 template <class VectorType>
2034 template <class BlockVector2>
2035 void BlockVectorBase<VectorType>::equ (const value_type a,
2036  const BlockVector2 &v)
2037 {
2038 
2039  AssertIsFinite(a);
2040 
2041  Assert (n_blocks() == v.n_blocks(),
2042  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
2043 
2044  for (size_type i=0; i<n_blocks(); ++i)
2045  components[i].equ( a, v.components[i]);
2046 }
2047 
2048 
2049 
2050 template <class VectorType>
2052 {
2053  for (size_type i=0; i<n_blocks(); ++i)
2054  block(i).update_ghost_values ();
2055 }
2056 
2057 
2058 
2059 template <class VectorType>
2061 BlockVectorBase<VectorType>::operator = (const value_type s)
2062 {
2063 
2064  AssertIsFinite(s);
2065 
2066  for (size_type i=0; i<n_blocks(); ++i)
2067  components[i] = s;
2068 
2069  return *this;
2070 }
2071 
2072 
2073 template <class VectorType>
2076 {
2077  AssertDimension(n_blocks(), v.n_blocks());
2078 
2079  for (size_type i=0; i<n_blocks(); ++i)
2080  components[i] = v.components[i];
2081 
2082  return *this;
2083 }
2084 
2085 
2086 template <class VectorType>
2087 template <class VectorType2>
2090 {
2091  AssertDimension(n_blocks(), v.n_blocks());
2092 
2093  for (size_type i=0; i<n_blocks(); ++i)
2094  components[i] = v.components[i];
2095 
2096  return *this;
2097 }
2098 
2099 
2100 
2101 template <class VectorType>
2103 BlockVectorBase<VectorType>::operator = (const VectorType &v)
2104 {
2105  Assert (size() == v.size(),
2106  ExcDimensionMismatch(size(), v.size()));
2107 
2108  size_type index_v = 0;
2109  for (size_type b=0; b<n_blocks(); ++b)
2110  for (size_type i=0; i<block(b).size(); ++i, ++index_v)
2111  block(b)(i) = v(index_v);
2112 
2113  return *this;
2114 }
2115 
2116 
2117 
2118 template <class VectorType>
2119 template <class VectorType2>
2120 inline
2121 bool
2124 {
2125  Assert (block_indices == v.block_indices, ExcDifferentBlockIndices());
2126 
2127  for (size_type i=0; i<n_blocks(); ++i)
2128  if ( ! (components[i] == v.components[i]))
2129  return false;
2130 
2131  return true;
2132 }
2133 
2134 
2135 
2136 template <class VectorType>
2137 inline
2139 BlockVectorBase<VectorType>::operator *= (const value_type factor)
2140 {
2141 
2142  AssertIsFinite(factor);
2143 
2144  for (size_type i=0; i<n_blocks(); ++i)
2145  components[i] *= factor;
2146 
2147  return *this;
2148 }
2149 
2150 
2151 
2152 template <class VectorType>
2153 inline
2155 BlockVectorBase<VectorType>::operator /= (const value_type factor)
2156 {
2157 
2158  AssertIsFinite(factor);
2159  Assert (factor != 0., ExcDivideByZero() );
2160 
2161  for (size_type i=0; i<n_blocks(); ++i)
2162  components[i] /= factor;
2163 
2164  return *this;
2165 }
2166 
2167 
2168 template <class VectorType>
2169 inline
2170 typename BlockVectorBase<VectorType>::value_type
2171 BlockVectorBase<VectorType>::operator() (const size_type i) const
2172 {
2173  const std::pair<unsigned int,size_type> local_index
2174  = block_indices.global_to_local (i);
2175  return components[local_index.first](local_index.second);
2176 }
2177 
2178 
2179 
2180 template <class VectorType>
2181 inline
2182 typename BlockVectorBase<VectorType>::reference
2184 {
2185  const std::pair<unsigned int,size_type> local_index
2186  = block_indices.global_to_local (i);
2187  return components[local_index.first](local_index.second);
2188 }
2189 
2190 
2191 
2192 template <class VectorType>
2193 inline
2194 typename BlockVectorBase<VectorType>::value_type
2195 BlockVectorBase<VectorType>::operator[] (const size_type i) const
2196 {
2197  return operator()(i);
2198 }
2199 
2200 
2201 
2202 template <class VectorType>
2203 inline
2204 typename BlockVectorBase<VectorType>::reference
2206 {
2207  return operator()(i);
2208 }
2209 
2210 
2211 
2212 template <typename VectorType>
2213 template <typename OtherNumber>
2214 inline
2215 void BlockVectorBase<VectorType>::extract_subvector_to (const std::vector<size_type> &indices,
2216  std::vector<OtherNumber> &values) const
2217 {
2218  for (size_type i = 0; i < indices.size(); ++i)
2219  values[i] = operator()(indices[i]);
2220 }
2221 
2222 
2223 
2224 template <typename VectorType>
2225 template <typename ForwardIterator, typename OutputIterator>
2226 inline
2227 void BlockVectorBase<VectorType>::extract_subvector_to (ForwardIterator indices_begin,
2228  const ForwardIterator indices_end,
2229  OutputIterator values_begin) const
2230 {
2231  while (indices_begin != indices_end)
2232  {
2233  *values_begin = operator()(*indices_begin);
2234  indices_begin++;
2235  values_begin++;
2236  }
2237 }
2238 
2239 #endif // DOXYGEN
2240 
2241 DEAL_II_NAMESPACE_CLOSE
2242 
2243 #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:82
std::random_access_iterator_tag iterator_type
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1202
bool operator==(const BlockVectorBase< VectorType2 > &v) const
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)
bool operator==(const Iterator< BlockVectorType, OtherConstness > &i) const
real_type linfty_norm() const
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:679
size_type size() const
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2205
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1597
real_type norm_sqr() const
void collect_sizes()
unsigned int n_blocks() const
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
static::ExceptionBase & ExcDivideByZero()
value_type operator[](const size_type i) const
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
bool operator>(const Iterator< BlockVectorType, OtherConstness > &i) const
bool in_local_range(const size_type global_index) const
unsigned int global_dof_index
Definition: types.h:88
#define Assert(cond, exc)
Definition: exceptions.h:349
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
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:594
bool operator>=(const Iterator< BlockVectorType, OtherConstness > &i) 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
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
BlockVectorBase & operator*=(const value_type factor)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Types< BlockVectorType, Constness >::value_type value_type
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
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
BlockType & block(const unsigned int i)
size_type local_to_global(const unsigned int block, const size_type index) const
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:1253
real_type l1_norm() const
void equ(const value_type a, const BlockVector2 &V)
Iterator & operator=(const Iterator &c)
Iterator & operator-=(const difference_type &d)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
std::size_t memory_consumption() const
dereference_type operator[](const difference_type d) const
Types< BlockVectorType, Constness >::BlockVector BlockVector