Reference documentation for deal.II version Git 2618e0f 2017-11-23 17:25:26 +0100
block_vector_base.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2017 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 
827  value_type add_and_dot (const value_type a,
828  const BlockVectorBase &V,
829  const BlockVectorBase &W);
830 
835  bool in_local_range (const size_type global_index) const;
836 
842  bool all_zero () const;
843 
849  bool is_non_negative () const;
850 
855  operator += (const BlockVectorBase &V);
856 
861  operator -= (const BlockVectorBase &V);
862 
863 
868  template <typename Number>
869  void add (const std::vector<size_type> &indices,
870  const std::vector<Number> &values);
871 
876  template <typename Number>
877  void add (const std::vector<size_type> &indices,
878  const Vector<Number> &values);
879 
885  template <typename Number>
886  void add (const size_type n_elements,
887  const size_type *indices,
888  const Number *values);
889 
894  void add (const value_type s);
895 
899  void add (const value_type a, const BlockVectorBase &V);
900 
904  void add (const value_type a, const BlockVectorBase &V,
905  const value_type b, const BlockVectorBase &W);
906 
910  void sadd (const value_type s, const BlockVectorBase &V);
911 
915  void sadd (const value_type s, const value_type a, const BlockVectorBase &V);
916 
920  void sadd (const value_type s, const value_type a,
921  const BlockVectorBase &V,
922  const value_type b, const BlockVectorBase &W);
923 
927  void sadd (const value_type s, const value_type a,
928  const BlockVectorBase &V,
929  const value_type b, const BlockVectorBase &W,
930  const value_type c, const BlockVectorBase &X);
931 
935  BlockVectorBase &operator *= (const value_type factor);
936 
940  BlockVectorBase &operator /= (const value_type factor);
941 
946  template <class BlockVector2>
947  void scale (const BlockVector2 &v);
948 
952  template <class BlockVector2>
953  void equ (const value_type a, const BlockVector2 &V);
954 
958  void equ (const value_type a, const BlockVectorBase &V,
959  const value_type b, const BlockVectorBase &W);
960 
965  void update_ghost_values () const;
966 
971  std::size_t memory_consumption () const;
972 
973 protected:
977  std::vector<VectorType> components;
978 
984 
988  template <typename N, bool C>
989  friend class ::internal::BlockVectorIterators::Iterator;
990 
991  template <typename> friend class BlockVectorBase;
992 };
993 
994 
997 /*----------------------- Inline functions ----------------------------------*/
998 
999 
1000 #ifndef DOXYGEN
1001 namespace internal
1002 {
1003  namespace BlockVectorIterators
1004  {
1005  template <class BlockVectorType, bool Constness>
1006  inline
1008  Iterator (const Iterator<BlockVectorType,Constness> &c)
1009  :
1010  parent (c.parent),
1011  global_index (c.global_index),
1012  current_block (c.current_block),
1013  index_within_block (c.index_within_block),
1014  next_break_forward (c.next_break_forward),
1015  next_break_backward (c.next_break_backward)
1016  {}
1017 
1018 
1019 
1020  template <class BlockVectorType, bool Constness>
1021  inline
1022  Iterator<BlockVectorType,Constness>::
1023  Iterator (const Iterator<BlockVectorType,!Constness> &c)
1024  :
1025  parent (c.parent),
1027  current_block (c.current_block),
1028  index_within_block (c.index_within_block),
1029  next_break_forward (c.next_break_forward),
1030  next_break_backward (c.next_break_backward)
1031  {
1032  // Only permit copy-constructing const iterators from non-const
1033  // iterators, and not vice versa (i.e., Constness must always be
1034  // true). As mentioned above, try to check this at compile time if C++11
1035  // support is available.
1036  static_assert(Constness == true,
1037  "Constructing a non-const iterator from a const iterator "
1038  "does not make sense.");
1039  }
1040 
1041 
1042 
1043  template <class BlockVectorType, bool Constness>
1044  inline
1045  Iterator<BlockVectorType,Constness>::
1046  Iterator (BlockVector &parent,
1047  const size_type global_index,
1048  const size_type current_block,
1049  const size_type index_within_block,
1050  const size_type next_break_forward,
1051  const size_type next_break_backward)
1052  :
1053  parent (&parent),
1054  global_index (global_index),
1055  current_block (current_block),
1056  index_within_block (index_within_block),
1057  next_break_forward (next_break_forward),
1058  next_break_backward (next_break_backward)
1059  {
1060  }
1061 
1062 
1063 
1064  template <class BlockVectorType, bool Constness>
1065  inline
1066  Iterator<BlockVectorType,Constness> &
1067  Iterator<BlockVectorType,Constness>::
1068  operator = (const Iterator &c)
1069  {
1070  parent = c.parent;
1071  global_index = c.global_index;
1072  index_within_block = c.index_within_block;
1073  current_block = c.current_block;
1074  next_break_forward = c.next_break_forward;
1075  next_break_backward = c.next_break_backward;
1076 
1077  return *this;
1078  }
1079 
1080 
1081 
1082  template <class BlockVectorType, bool Constness>
1083  inline
1084  typename Iterator<BlockVectorType,Constness>::dereference_type
1085  Iterator<BlockVectorType,Constness>::operator * () const
1086  {
1087  return parent->block(current_block)(index_within_block);
1088  }
1089 
1090 
1091 
1092  template <class BlockVectorType, bool Constness>
1093  inline
1094  typename Iterator<BlockVectorType,Constness>::dereference_type
1095  Iterator<BlockVectorType,Constness>::operator [] (const difference_type d) const
1096  {
1097  // if the index pointed to is
1098  // still within the block we
1099  // currently point into, then we
1100  // can save the computation of
1101  // the block
1102  if ((global_index+d >= next_break_backward) &&
1103  (global_index+d <= next_break_forward))
1104  return parent->block(current_block)(index_within_block + d);
1105 
1106  // if the index is not within the
1107  // block of the block vector into
1108  // which we presently point, then
1109  // there is no way: we have to
1110  // search for the block. this can
1111  // be done through the parent
1112  // class as well.
1113  return (*parent)(global_index+d);
1114  }
1115 
1116 
1117 
1118  template <class BlockVectorType, bool Constness>
1119  inline
1120  Iterator<BlockVectorType,Constness> &
1121  Iterator<BlockVectorType,Constness>::operator ++ ()
1122  {
1123  move_forward ();
1124  return *this;
1125  }
1126 
1127 
1128 
1129  template <class BlockVectorType, bool Constness>
1130  inline
1131  Iterator<BlockVectorType,Constness>
1132  Iterator<BlockVectorType,Constness>::operator ++ (int)
1133  {
1134  const Iterator old_value = *this;
1135  move_forward ();
1136  return old_value;
1137  }
1138 
1139 
1140 
1141  template <class BlockVectorType, bool Constness>
1142  inline
1143  Iterator<BlockVectorType,Constness> &
1144  Iterator<BlockVectorType,Constness>::operator -- ()
1145  {
1146  move_backward ();
1147  return *this;
1148  }
1149 
1150 
1151 
1152  template <class BlockVectorType, bool Constness>
1153  inline
1154  Iterator<BlockVectorType,Constness>
1155  Iterator<BlockVectorType,Constness>::operator -- (int)
1156  {
1157  const Iterator old_value = *this;
1158  move_backward ();
1159  return old_value;
1160  }
1161 
1162 
1163 
1164  template <class BlockVectorType, bool Constness>
1165  template <bool OtherConstness>
1166  inline
1167  bool
1168  Iterator<BlockVectorType,Constness>::
1169  operator == (const Iterator<BlockVectorType, OtherConstness> &i) const
1170  {
1171  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1172 
1173  return (global_index == i.global_index);
1174  }
1175 
1176 
1177 
1178  template <class BlockVectorType, bool Constness>
1179  template <bool OtherConstness>
1180  inline
1181  bool
1182  Iterator<BlockVectorType,Constness>::
1183  operator != (const Iterator<BlockVectorType, OtherConstness> &i) const
1184  {
1185  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1186 
1187  return (global_index != i.global_index);
1188  }
1189 
1190 
1191 
1192  template <class BlockVectorType, bool Constness>
1193  template <bool OtherConstness>
1194  inline
1195  bool
1196  Iterator<BlockVectorType,Constness>::
1197  operator < (const Iterator<BlockVectorType, OtherConstness> &i) const
1198  {
1199  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1200 
1201  return (global_index < i.global_index);
1202  }
1203 
1204 
1205 
1206  template <class BlockVectorType, bool Constness>
1207  template <bool OtherConstness>
1208  inline
1209  bool
1210  Iterator<BlockVectorType,Constness>::
1211  operator <= (const Iterator<BlockVectorType, OtherConstness> &i) const
1212  {
1213  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1214 
1215  return (global_index <= i.global_index);
1216  }
1217 
1218 
1219 
1220  template <class BlockVectorType, bool Constness>
1221  template <bool OtherConstness>
1222  inline
1223  bool
1224  Iterator<BlockVectorType,Constness>::
1225  operator > (const Iterator<BlockVectorType, OtherConstness> &i) const
1226  {
1227  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1228 
1229  return (global_index > i.global_index);
1230  }
1231 
1232 
1233 
1234  template <class BlockVectorType, bool Constness>
1235  template <bool OtherConstness>
1236  inline
1237  bool
1238  Iterator<BlockVectorType,Constness>::
1239  operator >= (const Iterator<BlockVectorType, OtherConstness> &i) const
1240  {
1241  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1242 
1243  return (global_index >= i.global_index);
1244  }
1245 
1246 
1247 
1248  template <class BlockVectorType, bool Constness>
1249  template <bool OtherConstness>
1250  inline
1251  typename Iterator<BlockVectorType,Constness>::difference_type
1252  Iterator<BlockVectorType,Constness>::
1253  operator - (const Iterator<BlockVectorType, OtherConstness> &i) const
1254  {
1255  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1256 
1257  return (static_cast<signed int>(global_index) -
1258  static_cast<signed int>(i.global_index));
1259  }
1260 
1261 
1262 
1263  template <class BlockVectorType, bool Constness>
1264  inline
1265  Iterator<BlockVectorType,Constness>
1266  Iterator<BlockVectorType,Constness>::
1267  operator + (const difference_type &d) const
1268  {
1269  // if the index pointed to is
1270  // still within the block we
1271  // currently point into, then we
1272  // can save the computation of
1273  // the block
1274  if ((global_index+d >= next_break_backward) &&
1275  (global_index+d <= next_break_forward))
1276  return Iterator (*parent, global_index+d, current_block,
1277  index_within_block+d,
1278  next_break_forward, next_break_backward);
1279  else
1280  // outside present block, so
1281  // have to seek new block
1282  // anyway
1283  return Iterator (*parent, global_index+d);
1284  }
1285 
1286 
1287 
1288  template <class BlockVectorType, bool Constness>
1289  inline
1290  Iterator<BlockVectorType,Constness>
1291  Iterator<BlockVectorType,Constness>::
1292  operator - (const difference_type &d) const
1293  {
1294  // if the index pointed to is
1295  // still within the block we
1296  // currently point into, then we
1297  // can save the computation of
1298  // the block
1299  if ((global_index-d >= next_break_backward) &&
1300  (global_index-d <= next_break_forward))
1301  return Iterator (*parent, global_index-d, current_block,
1302  index_within_block-d,
1303  next_break_forward, next_break_backward);
1304  else
1305  // outside present block, so
1306  // have to seek new block
1307  // anyway
1308  return Iterator (*parent, global_index-d);
1309  }
1310 
1311 
1312 
1313  template <class BlockVectorType, bool Constness>
1314  inline
1315  Iterator<BlockVectorType,Constness> &
1316  Iterator<BlockVectorType,Constness>::
1317  operator += (const difference_type &d)
1318  {
1319  // if the index pointed to is
1320  // still within the block we
1321  // currently point into, then we
1322  // can save the computation of
1323  // the block
1324  if ((global_index+d >= next_break_backward) &&
1325  (global_index+d <= next_break_forward))
1326  {
1327  global_index += d;
1328  index_within_block += d;
1329  }
1330  else
1331  // outside present block, so
1332  // have to seek new block
1333  // anyway
1334  *this = Iterator (*parent, global_index+d);
1335 
1336  return *this;
1337  }
1338 
1339 
1340 
1341  template <class BlockVectorType, bool Constness>
1342  inline
1343  Iterator<BlockVectorType,Constness> &
1344  Iterator<BlockVectorType,Constness>::
1345  operator -= (const difference_type &d)
1346  {
1347  // if the index pointed to is
1348  // still within the block we
1349  // currently point into, then we
1350  // can save the computation of
1351  // the block
1352  if ((global_index-d >= next_break_backward) &&
1353  (global_index-d <= next_break_forward))
1354  {
1355  global_index -= d;
1356  index_within_block -= d;
1357  }
1358  else
1359  // outside present block, so
1360  // have to seek new block
1361  // anyway
1362  *this = Iterator (*parent, global_index-d);
1363 
1364  return *this;
1365  }
1366 
1367 
1368  template <class BlockVectorType, bool Constness>
1369  Iterator<BlockVectorType,Constness>::
1370  Iterator (BlockVector &parent,
1371  const size_type global_index)
1372  :
1373  parent (&parent),
1374  global_index (global_index)
1375  {
1376  // find which block we are
1377  // in. for this, take into
1378  // account that it happens at
1379  // times that people want to
1380  // initialize iterators
1381  // past-the-end
1382  if (global_index < parent.size())
1383  {
1384  const std::pair<size_type, size_type>
1385  indices = parent.block_indices.global_to_local(global_index);
1386  current_block = indices.first;
1387  index_within_block = indices.second;
1388 
1389  next_break_backward
1390  = parent.block_indices.local_to_global (current_block, 0);
1391  next_break_forward
1392  = (parent.block_indices.local_to_global (current_block, 0)
1393  +parent.block_indices.block_size(current_block)-1);
1394  }
1395  else
1396  // past the end. only have one
1397  // value for this
1398  {
1399  this->global_index = parent.size ();
1400  current_block = parent.n_blocks();
1401  index_within_block = 0;
1402  next_break_backward = global_index;
1403  next_break_forward = numbers::invalid_size_type;
1404  };
1405  }
1406 
1407 
1408 
1409  template <class BlockVectorType, bool Constness>
1410  void
1411  Iterator<BlockVectorType,Constness>::move_forward ()
1412  {
1413  if (global_index != next_break_forward)
1414  ++index_within_block;
1415  else
1416  {
1417  // ok, we traverse a boundary
1418  // between blocks:
1419  index_within_block = 0;
1420  ++current_block;
1421 
1422  // break backwards is now old
1423  // break forward
1424  next_break_backward = next_break_forward+1;
1425 
1426  // compute new break forward
1427  if (current_block < parent->block_indices.size())
1428  next_break_forward
1429  += parent->block_indices.block_size(current_block);
1430  else
1431  // if we are beyond the end,
1432  // then move the next
1433  // boundary arbitrarily far
1434  // away
1435  next_break_forward = numbers::invalid_size_type;
1436  };
1437 
1438  ++global_index;
1439  }
1440 
1441 
1442 
1443  template <class BlockVectorType, bool Constness>
1444  void
1445  Iterator<BlockVectorType,Constness>::move_backward ()
1446  {
1447  if (global_index != next_break_backward)
1448  --index_within_block;
1449  else if (current_block != 0)
1450  {
1451  // ok, we traverse a boundary
1452  // between blocks:
1453  --current_block;
1454  index_within_block = parent->block_indices.block_size(current_block)-1;
1455 
1456  // break forwards is now old
1457  // break backward
1458  next_break_forward = next_break_backward-1;
1459 
1460  // compute new break forward
1461  next_break_backward
1462  -= parent->block_indices.block_size (current_block);
1463  }
1464  else
1465  // current block was 0, we now
1466  // get into unspecified terrain
1467  {
1468  --current_block;
1469  index_within_block = numbers::invalid_size_type;
1470  next_break_forward = 0;
1471  next_break_backward = 0;
1472  };
1473 
1474  --global_index;
1475  }
1476 
1477 
1478  } // namespace BlockVectorIterators
1479 
1480 } //namespace internal
1481 
1482 
1483 
1484 template <class VectorType>
1485 inline
1486 std::size_t
1488 {
1489  return block_indices.total_size();
1490 }
1491 
1492 
1493 
1494 template <class VectorType>
1495 inline
1496 IndexSet
1498 {
1499  IndexSet is (size());
1500 
1501  // copy index sets from blocks into the global one, shifted
1502  // by the appropriate amount for each block
1503  for (unsigned int b=0; b<n_blocks(); ++b)
1504  {
1505  IndexSet x = block(b).locally_owned_elements();
1506  is.add_indices(x, block_indices.block_start(b));
1507  }
1508 
1509  is.compress();
1510 
1511  return is;
1512 }
1513 
1514 
1515 
1516 template <class VectorType>
1517 inline
1518 unsigned int
1520 {
1521  return block_indices.size();
1522 }
1523 
1524 
1525 template <class VectorType>
1526 inline
1528 BlockVectorBase<VectorType>::block (const unsigned int i)
1529 {
1530  Assert(i<n_blocks(), ExcIndexRange(i,0,n_blocks()));
1531 
1532  return components[i];
1533 }
1534 
1535 
1536 
1537 template <class VectorType>
1538 inline
1540 BlockVectorBase<VectorType>::block (const unsigned int i) const
1541 {
1542  Assert(i<n_blocks(), ExcIndexRange(i,0,n_blocks()));
1543 
1544  return components[i];
1545 }
1546 
1547 
1548 
1549 template <class VectorType>
1550 inline
1551 const BlockIndices &
1553 {
1554  return block_indices;
1555 }
1556 
1557 
1558 template <class VectorType>
1559 inline
1560 void
1562 {
1563  std::vector<size_type> sizes (n_blocks());
1564 
1565  for (size_type i=0; i<n_blocks(); ++i)
1566  sizes[i] = block(i).size();
1567 
1568  block_indices.reinit(sizes);
1569 }
1570 
1571 
1572 
1573 template <class VectorType>
1574 inline
1575 void
1577 {
1578  for (unsigned int i=0; i<n_blocks(); ++i)
1579  block(i).compress (operation);
1580 }
1581 
1582 
1583 
1584 template <class VectorType>
1585 inline
1588 {
1589  return iterator(*this, 0U);
1590 }
1591 
1592 
1593 
1594 template <class VectorType>
1595 inline
1598 {
1599  return const_iterator(*this, 0U);
1600 }
1601 
1602 
1603 template <class VectorType>
1604 inline
1607 {
1608  return iterator(*this, size());
1609 }
1610 
1611 
1612 
1613 template <class VectorType>
1614 inline
1617 {
1618  return const_iterator(*this, size());
1619 }
1620 
1621 
1622 template <class VectorType>
1623 inline
1624 bool
1626 (const size_type global_index) const
1627 {
1628  const std::pair<size_type,size_type> local_index
1629  = block_indices.global_to_local (global_index);
1630 
1631  return components[local_index.first].in_local_range (global_index);
1632 }
1633 
1634 
1635 template <class VectorType>
1636 bool
1638 {
1639  for (size_type i=0; i<n_blocks(); ++i)
1640  if (components[i].all_zero() == false)
1641  return false;
1642 
1643  return true;
1644 }
1645 
1646 
1647 
1648 template <class VectorType>
1649 bool
1651 {
1652  for (size_type i=0; i<n_blocks(); ++i)
1653  if (components[i].is_non_negative() == false)
1654  return false;
1655 
1656  return true;
1657 }
1658 
1659 
1660 
1661 template <class VectorType>
1662 typename BlockVectorBase<VectorType>::value_type
1665 {
1666  Assert (n_blocks() == v.n_blocks(),
1667  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1668 
1669  value_type sum = 0.;
1670  for (size_type i=0; i<n_blocks(); ++i)
1671  sum += components[i]*v.components[i];
1672 
1673  return sum;
1674 }
1675 
1676 
1677 template <class VectorType>
1680 {
1681  real_type sum = 0.;
1682  for (size_type i=0; i<n_blocks(); ++i)
1683  sum += components[i].norm_sqr();
1684 
1685  return sum;
1686 }
1687 
1688 
1689 
1690 template <class VectorType>
1691 typename BlockVectorBase<VectorType>::value_type
1693 {
1694  value_type sum = 0.;
1695  // need to do static_cast as otherwise it won't work with value_type=complex<T>
1696  for (size_type i=0; i<n_blocks(); ++i)
1697  sum += components[i].mean_value() * (typename numbers::NumberTraits<value_type>::real_type(components[i].size()));
1698 
1699  return sum/(typename numbers::NumberTraits<value_type>::real_type(size()));
1700 }
1701 
1702 
1703 
1704 template <class VectorType>
1707 {
1708  real_type sum = 0.;
1709  for (size_type i=0; i<n_blocks(); ++i)
1710  sum += components[i].l1_norm();
1711 
1712  return sum;
1713 }
1714 
1715 
1716 
1717 template <class VectorType>
1720 {
1721  return std::sqrt(norm_sqr());
1722 }
1723 
1724 
1725 
1726 template <class VectorType>
1729 {
1730  real_type sum = 0.;
1731  for (size_type i=0; i<n_blocks(); ++i)
1732  {
1733  value_type newval = components[i].linfty_norm();
1734  if (sum<newval)
1735  sum = newval;
1736  }
1737  return sum;
1738 }
1739 
1740 
1741 
1742 template <class VectorType>
1743 typename BlockVectorBase<VectorType>::value_type
1745 add_and_dot (const typename BlockVectorBase<VectorType>::value_type a,
1746  const BlockVectorBase<VectorType> &V,
1747  const BlockVectorBase<VectorType> &W)
1748 {
1749  AssertDimension (n_blocks(), V.n_blocks());
1750  AssertDimension (n_blocks(), W.n_blocks());
1751 
1752  value_type sum = 0.;
1753  for (size_type i=0; i<n_blocks(); ++i)
1754  sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1755 
1756  return sum;
1757 }
1758 
1759 
1760 
1761 template <class VectorType>
1764 {
1765  Assert (n_blocks() == v.n_blocks(),
1766  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1767 
1768  for (size_type i=0; i<n_blocks(); ++i)
1769  {
1770  components[i] += v.components[i];
1771  }
1772 
1773  return *this;
1774 }
1775 
1776 
1777 
1778 template <class VectorType>
1781 {
1782  Assert (n_blocks() == v.n_blocks(),
1783  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1784 
1785  for (size_type i=0; i<n_blocks(); ++i)
1786  {
1787  components[i] -= v.components[i];
1788  }
1789  return *this;
1790 }
1791 
1792 
1793 
1794 template <class VectorType>
1795 template <typename Number>
1796 inline
1797 void
1798 BlockVectorBase<VectorType>::add (const std::vector<size_type> &indices,
1799  const std::vector<Number> &values)
1800 {
1801  Assert (indices.size() == values.size(),
1802  ExcDimensionMismatch(indices.size(), values.size()));
1803  add (indices.size(), indices.data(), values.data());
1804 }
1805 
1806 
1807 
1808 template <class VectorType>
1809 template <typename Number>
1810 inline
1811 void
1812 BlockVectorBase<VectorType>::add (const std::vector<size_type> &indices,
1813  const Vector<Number> &values)
1814 {
1815  Assert (indices.size() == values.size(),
1816  ExcDimensionMismatch(indices.size(), values.size()));
1817  const size_type n_indices = indices.size();
1818  for (size_type i=0; i<n_indices; ++i)
1819  (*this)(indices[i]) += values(i);
1820 }
1821 
1822 
1823 
1824 template <class VectorType>
1825 template <typename Number>
1826 inline
1827 void
1828 BlockVectorBase<VectorType>::add (const size_type n_indices,
1829  const size_type *indices,
1830  const Number *values)
1831 {
1832  for (size_type i=0; i<n_indices; ++i)
1833  (*this)(indices[i]) += values[i];
1834 }
1835 
1836 
1837 
1838 template <class VectorType>
1839 void BlockVectorBase<VectorType>::add (const value_type a)
1840 {
1841  AssertIsFinite(a);
1842 
1843  for (size_type i=0; i<n_blocks(); ++i)
1844  {
1845  components[i].add(a);
1846  }
1847 }
1848 
1849 
1850 
1851 template <class VectorType>
1852 void BlockVectorBase<VectorType>::add (const value_type a,
1853  const BlockVectorBase<VectorType> &v)
1854 {
1855 
1856  AssertIsFinite(a);
1857 
1858  Assert (n_blocks() == v.n_blocks(),
1859  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1860 
1861  for (size_type i=0; i<n_blocks(); ++i)
1862  {
1863  components[i].add(a, v.components[i]);
1864  }
1865 }
1866 
1867 
1868 
1869 template <class VectorType>
1870 void BlockVectorBase<VectorType>::add (const value_type a,
1871  const BlockVectorBase<VectorType> &v,
1872  const value_type b,
1873  const BlockVectorBase<VectorType> &w)
1874 {
1875 
1876  AssertIsFinite(a);
1877  AssertIsFinite(b);
1878 
1879  Assert (n_blocks() == v.n_blocks(),
1880  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1881  Assert (n_blocks() == w.n_blocks(),
1882  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1883 
1884 
1885  for (size_type i=0; i<n_blocks(); ++i)
1886  {
1887  components[i].add(a, v.components[i], b, w.components[i]);
1888  }
1889 }
1890 
1891 
1892 
1893 template <class VectorType>
1894 void BlockVectorBase<VectorType>::sadd (const value_type x,
1895  const BlockVectorBase<VectorType> &v)
1896 {
1897 
1898  AssertIsFinite(x);
1899 
1900  Assert (n_blocks() == v.n_blocks(),
1901  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1902 
1903  for (size_type i=0; i<n_blocks(); ++i)
1904  {
1905  components[i].sadd(x, v.components[i]);
1906  }
1907 }
1908 
1909 
1910 
1911 template <class VectorType>
1912 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
1913  const BlockVectorBase<VectorType> &v)
1914 {
1915 
1916  AssertIsFinite(x);
1917  AssertIsFinite(a);
1918 
1919  Assert (n_blocks() == v.n_blocks(),
1920  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1921 
1922  for (size_type i=0; i<n_blocks(); ++i)
1923  {
1924  components[i].sadd(x, a, v.components[i]);
1925  }
1926 }
1927 
1928 
1929 
1930 template <class VectorType>
1931 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
1932  const BlockVectorBase<VectorType> &v,
1933  const value_type b,
1934  const BlockVectorBase<VectorType> &w)
1935 {
1936 
1937  AssertIsFinite(x);
1938  AssertIsFinite(a);
1939  AssertIsFinite(b);
1940 
1941  Assert (n_blocks() == v.n_blocks(),
1942  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1943  Assert (n_blocks() == w.n_blocks(),
1944  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1945 
1946  for (size_type i=0; i<n_blocks(); ++i)
1947  {
1948  components[i].sadd(x, a, v.components[i], b, w.components[i]);
1949  }
1950 }
1951 
1952 
1953 
1954 template <class VectorType>
1955 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
1956  const BlockVectorBase<VectorType> &v,
1957  const value_type b,
1958  const BlockVectorBase<VectorType> &w,
1959  const value_type c,
1960  const BlockVectorBase<VectorType> &y)
1961 {
1962 
1963  AssertIsFinite(x);
1964  AssertIsFinite(a);
1965  AssertIsFinite(b);
1966  AssertIsFinite(c);
1967 
1968  Assert (n_blocks() == v.n_blocks(),
1969  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1970  Assert (n_blocks() == w.n_blocks(),
1971  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1972  Assert (n_blocks() == y.n_blocks(),
1973  ExcDimensionMismatch(n_blocks(), y.n_blocks()));
1974 
1975  for (size_type i=0; i<n_blocks(); ++i)
1976  {
1977  components[i].sadd(x, a, v.components[i],
1978  b, w.components[i], c, y.components[i]);
1979  }
1980 }
1981 
1982 
1983 
1984 template <class VectorType>
1985 template <class BlockVector2>
1986 void BlockVectorBase<VectorType>::scale (const BlockVector2 &v)
1987 {
1988  Assert (n_blocks() == v.n_blocks(),
1989  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1990  for (size_type i=0; i<n_blocks(); ++i)
1991  components[i].scale(v.block(i));
1992 }
1993 
1994 
1995 
1996 template <class VectorType>
1997 void BlockVectorBase<VectorType>::equ (const value_type a,
1998  const BlockVectorBase<VectorType> &v,
1999  const value_type b,
2000  const BlockVectorBase<VectorType> &w)
2001 {
2002 
2003  AssertIsFinite(a);
2004  AssertIsFinite(b);
2005 
2006  Assert (n_blocks() == v.n_blocks(),
2007  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
2008  Assert (n_blocks() == w.n_blocks(),
2009  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
2010 
2011  for (size_type i=0; i<n_blocks(); ++i)
2012  {
2013  components[i].equ( a, v.components[i], b, w.components[i]);
2014  }
2015 }
2016 
2017 
2018 
2019 template <class VectorType>
2020 std::size_t
2022 {
2023  return (MemoryConsumption::memory_consumption (this->block_indices)
2024  +
2025  MemoryConsumption::memory_consumption (this->components));
2026 }
2027 
2028 
2029 
2030 template <class VectorType>
2031 template <class BlockVector2>
2032 void BlockVectorBase<VectorType>::equ (const value_type a,
2033  const BlockVector2 &v)
2034 {
2035 
2036  AssertIsFinite(a);
2037 
2038  Assert (n_blocks() == v.n_blocks(),
2039  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
2040 
2041  for (size_type i=0; i<n_blocks(); ++i)
2042  components[i].equ( a, v.components[i]);
2043 }
2044 
2045 
2046 
2047 template <class VectorType>
2049 {
2050  for (size_type i=0; i<n_blocks(); ++i)
2051  block(i).update_ghost_values ();
2052 }
2053 
2054 
2055 
2056 template <class VectorType>
2058 BlockVectorBase<VectorType>::operator = (const value_type s)
2059 {
2060 
2061  AssertIsFinite(s);
2062 
2063  for (size_type i=0; i<n_blocks(); ++i)
2064  components[i] = s;
2065 
2066  return *this;
2067 }
2068 
2069 
2070 template <class VectorType>
2073 {
2074  AssertDimension(n_blocks(), v.n_blocks());
2075 
2076  for (size_type i=0; i<n_blocks(); ++i)
2077  components[i] = v.components[i];
2078 
2079  return *this;
2080 }
2081 
2082 
2083 template <class VectorType>
2084 template <class VectorType2>
2087 {
2088  AssertDimension(n_blocks(), v.n_blocks());
2089 
2090  for (size_type i=0; i<n_blocks(); ++i)
2091  components[i] = v.components[i];
2092 
2093  return *this;
2094 }
2095 
2096 
2097 
2098 template <class VectorType>
2100 BlockVectorBase<VectorType>::operator = (const VectorType &v)
2101 {
2102  Assert (size() == v.size(),
2103  ExcDimensionMismatch(size(), v.size()));
2104 
2105  size_type index_v = 0;
2106  for (size_type b=0; b<n_blocks(); ++b)
2107  for (size_type i=0; i<block(b).size(); ++i, ++index_v)
2108  block(b)(i) = v(index_v);
2109 
2110  return *this;
2111 }
2112 
2113 
2114 
2115 template <class VectorType>
2116 template <class VectorType2>
2117 inline
2118 bool
2121 {
2122  Assert (block_indices == v.block_indices, ExcDifferentBlockIndices());
2123 
2124  for (size_type i=0; i<n_blocks(); ++i)
2125  if ( ! (components[i] == v.components[i]))
2126  return false;
2127 
2128  return true;
2129 }
2130 
2131 
2132 
2133 template <class VectorType>
2134 inline
2136 BlockVectorBase<VectorType>::operator *= (const value_type factor)
2137 {
2138 
2139  AssertIsFinite(factor);
2140 
2141  for (size_type i=0; i<n_blocks(); ++i)
2142  components[i] *= factor;
2143 
2144  return *this;
2145 }
2146 
2147 
2148 
2149 template <class VectorType>
2150 inline
2152 BlockVectorBase<VectorType>::operator /= (const value_type factor)
2153 {
2154 
2155  AssertIsFinite(factor);
2156  Assert (factor != 0., ExcDivideByZero() );
2157 
2158  for (size_type i=0; i<n_blocks(); ++i)
2159  components[i] /= factor;
2160 
2161  return *this;
2162 }
2163 
2164 
2165 template <class VectorType>
2166 inline
2167 typename BlockVectorBase<VectorType>::value_type
2168 BlockVectorBase<VectorType>::operator() (const size_type i) const
2169 {
2170  const std::pair<unsigned int,size_type> local_index
2171  = block_indices.global_to_local (i);
2172  return components[local_index.first](local_index.second);
2173 }
2174 
2175 
2176 
2177 template <class VectorType>
2178 inline
2179 typename BlockVectorBase<VectorType>::reference
2181 {
2182  const std::pair<unsigned int,size_type> local_index
2183  = block_indices.global_to_local (i);
2184  return components[local_index.first](local_index.second);
2185 }
2186 
2187 
2188 
2189 template <class VectorType>
2190 inline
2191 typename BlockVectorBase<VectorType>::value_type
2192 BlockVectorBase<VectorType>::operator[] (const size_type i) const
2193 {
2194  return operator()(i);
2195 }
2196 
2197 
2198 
2199 template <class VectorType>
2200 inline
2201 typename BlockVectorBase<VectorType>::reference
2203 {
2204  return operator()(i);
2205 }
2206 
2207 
2208 
2209 template <typename VectorType>
2210 template <typename OtherNumber>
2211 inline
2212 void BlockVectorBase<VectorType>::extract_subvector_to (const std::vector<size_type> &indices,
2213  std::vector<OtherNumber> &values) const
2214 {
2215  for (size_type i = 0; i < indices.size(); ++i)
2216  values[i] = operator()(indices[i]);
2217 }
2218 
2219 
2220 
2221 template <typename VectorType>
2222 template <typename ForwardIterator, typename OutputIterator>
2223 inline
2224 void BlockVectorBase<VectorType>::extract_subvector_to (ForwardIterator indices_begin,
2225  const ForwardIterator indices_end,
2226  OutputIterator values_begin) const
2227 {
2228  while (indices_begin != indices_end)
2229  {
2230  *values_begin = operator()(*indices_begin);
2231  indices_begin++;
2232  values_begin++;
2233  }
2234 }
2235 
2236 #endif // DOXYGEN
2237 
2238 DEAL_II_NAMESPACE_CLOSE
2239 
2240 #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:1183
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:672
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2028
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:337
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:582
bool operator>=(const Iterator< BlockVectorType, OtherConstness > &i) const
bool operator!=(const Iterator< BlockVectorType, OtherConstness > &i) const
std::size_t size() 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:1219
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