Reference documentation for deal.II version Git 05e4468 2017-09-21 10:18:23 +0200
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/exceptions.h>
22 #include <deal.II/base/subscriptor.h>
23 #include <deal.II/base/memory_consumption.h>
24 #include <deal.II/lac/exceptions.h>
25 #include <deal.II/lac/block_indices.h>
26 #include <deal.II/lac/vector.h>
27 #include <vector>
28 #include <iterator>
29 #include <cmath>
30 #include <cstddef>
31 
32 DEAL_II_NAMESPACE_OPEN
33 
34 
39 template <typename> class BlockVectorBase;
40 
41 
58 template <typename VectorType>
60 {
61 private:
62  struct yes_type
63  {
64  char c[1];
65  };
66  struct no_type
67  {
68  char c[2];
69  };
70 
75  template <typename T>
76  static yes_type check_for_block_vector (const BlockVectorBase<T> *);
77 
82  static no_type check_for_block_vector (...);
83 
84 public:
90  static const bool value = (sizeof(check_for_block_vector
91  ((VectorType *)nullptr))
92  ==
93  sizeof(yes_type));
94 };
95 
96 
97 // instantiation of the static member
98 template <typename VectorType>
100 
101 
102 
103 
104 namespace internal
105 {
106 
112  namespace BlockVectorIterators
113  {
118  template <class BlockVectorType, bool Constness>
119  struct Types
120  {
121  };
122 
123 
124 
131  template <class BlockVectorType>
132  struct Types<BlockVectorType,false>
133  {
138  typedef typename BlockVectorType::BlockType Vector;
139 
144  typedef BlockVectorType BlockVector;
145 
150 
155  typedef typename Vector::reference dereference_type;
156  };
157 
158 
159 
166  template <class BlockVectorType>
167  struct Types<BlockVectorType,true>
168  {
173  typedef const typename BlockVectorType::BlockType Vector;
174 
179  typedef const BlockVectorType BlockVector;
180 
185  typedef const typename BlockVector::value_type value_type;
186 
192  typedef value_type dereference_type;
193  };
194 
195 
213  template <class BlockVectorType, bool Constness>
214  class Iterator :
215  public std::iterator<std::random_access_iterator_tag,
216  typename Types<BlockVectorType,Constness>::value_type>
217  {
218  public:
223 
229  typedef
232 
238  typedef std::random_access_iterator_tag iterator_type;
239  typedef std::ptrdiff_t difference_type;
240  typedef typename BlockVectorType::reference reference;
241  typedef value_type *pointer;
242 
243  typedef
245  dereference_type;
246 
251  typedef
254 
264  const size_type global_index);
265 
274 
275 
279  Iterator (const Iterator &c);
280 
281  private:
287  const size_type global_index,
288  const size_type current_block,
289  const size_type index_within_block,
290  const size_type next_break_forward,
291  const size_type next_break_backward);
292 
293  public:
294 
298  Iterator &operator = (const Iterator &c);
299 
305  dereference_type operator * () const;
306 
311  dereference_type operator [] (const difference_type d) const;
312 
318 
324  Iterator operator ++ (int);
325 
331 
337  Iterator operator -- (int);
338 
343  template <bool OtherConstness>
345 
350  template <bool OtherConstness>
352 
358  template <bool OtherConstness>
359  bool operator < (const Iterator<BlockVectorType, OtherConstness> &i) const;
360 
364  template <bool OtherConstness>
365  bool operator <= (const Iterator<BlockVectorType, OtherConstness> &i) const;
366 
370  template <bool OtherConstness>
372 
376  template <bool OtherConstness>
378 
382  template <bool OtherConstness>
383  difference_type operator - (const Iterator<BlockVectorType, OtherConstness> &i) const;
384 
389  Iterator operator + (const difference_type &d) const;
390 
395  Iterator operator - (const difference_type &d) const;
396 
401  Iterator &operator += (const difference_type &d);
402 
407  Iterator &operator -= (const difference_type &d);
408 
419  "Your program tried to compare iterators pointing to "
420  "different block vectors. There is no reasonable way "
421  "to do this.");
422 
424  private:
431 
435  size_type global_index;
436 
441  unsigned int current_block;
442  size_type index_within_block;
443 
451  size_type next_break_backward;
452 
456  void move_forward ();
457 
461  void move_backward ();
462 
463 
467  template <typename, bool>
468  friend class Iterator;
469  };
470  } // namespace BlockVectorIterators
471 } // namespace internal
472 
473 
512 template <class VectorType>
513 class BlockVectorBase : public Subscriptor
514 {
515 public:
519  typedef VectorType BlockType;
520 
521  /*
522  * Declare standard types used in
523  * all containers. These types
524  * parallel those in the
525  * <tt>C++</tt> standard
526  * libraries
527  * <tt>std::vector<...></tt>
528  * class. This includes iterator
529  * types.
530  */
531  typedef typename BlockType::value_type value_type;
532  typedef value_type *pointer;
533  typedef const value_type *const_pointer;
534  typedef ::internal::BlockVectorIterators::Iterator<BlockVectorBase,false> iterator;
535  typedef ::internal::BlockVectorIterators::Iterator<BlockVectorBase,true> const_iterator;
536  typedef typename BlockType::reference reference;
537  typedef typename BlockType::const_reference const_reference;
538  typedef types::global_dof_index size_type;
539 
549  typedef typename BlockType::real_type real_type;
550 
554  BlockVectorBase () = default;
555 
559  BlockVectorBase (const BlockVectorBase &/*V*/) = default;
560 
566  BlockVectorBase (BlockVectorBase &&/*V*/) = default;
567 
574  void collect_sizes ();
575 
586  void compress (::VectorOperation::values operation);
587 
591  BlockType &
592  block (const unsigned int i);
593 
597  const BlockType &
598  block (const unsigned int i) const;
599 
605  const BlockIndices &
606  get_block_indices () const;
607 
611  unsigned int n_blocks () const;
612 
617  std::size_t size () const;
618 
636 
640  iterator begin ();
641 
646  const_iterator begin () const;
647 
651  iterator end ();
652 
657  const_iterator end () const;
658 
662  value_type operator() (const size_type i) const;
663 
667  reference operator() (const size_type i);
668 
674  value_type operator[] (const size_type i) const;
675 
681  reference operator[] (const size_type i);
682 
698  template <typename OtherNumber>
699  void extract_subvector_to (const std::vector<size_type> &indices,
700  std::vector<OtherNumber> &values) const;
701 
729  template <typename ForwardIterator, typename OutputIterator>
730  void extract_subvector_to (ForwardIterator indices_begin,
731  const ForwardIterator indices_end,
732  OutputIterator values_begin) const;
733 
738  BlockVectorBase &operator = (const value_type s);
739 
744  operator= (const BlockVectorBase &V);
745 
751  BlockVectorBase &operator= (BlockVectorBase &&/*V*/) = default;
752 
756  template <class VectorType2>
759 
764  operator = (const VectorType &v);
765 
770  template <class VectorType2>
771  bool
773 
777  value_type operator* (const BlockVectorBase &V) const;
778 
782  real_type norm_sqr () const;
783 
787  value_type mean_value () const;
788 
792  real_type l1_norm () const;
793 
798  real_type l2_norm () const;
799 
804  real_type linfty_norm () const;
805 
824  value_type add_and_dot (const value_type a,
825  const BlockVectorBase &V,
826  const BlockVectorBase &W);
827 
832  bool in_local_range (const size_type global_index) const;
833 
839  bool all_zero () const;
840 
846  bool is_non_negative () const;
847 
852  operator += (const BlockVectorBase &V);
853 
858  operator -= (const BlockVectorBase &V);
859 
860 
865  template <typename Number>
866  void add (const std::vector<size_type> &indices,
867  const std::vector<Number> &values);
868 
873  template <typename Number>
874  void add (const std::vector<size_type> &indices,
875  const Vector<Number> &values);
876 
882  template <typename Number>
883  void add (const size_type n_elements,
884  const size_type *indices,
885  const Number *values);
886 
891  void add (const value_type s);
892 
896  void add (const value_type a, const BlockVectorBase &V);
897 
901  void add (const value_type a, const BlockVectorBase &V,
902  const value_type b, const BlockVectorBase &W);
903 
907  void sadd (const value_type s, const BlockVectorBase &V);
908 
912  void sadd (const value_type s, const value_type a, const BlockVectorBase &V);
913 
917  void sadd (const value_type s, const value_type a,
918  const BlockVectorBase &V,
919  const value_type b, const BlockVectorBase &W);
920 
924  void sadd (const value_type s, const value_type a,
925  const BlockVectorBase &V,
926  const value_type b, const BlockVectorBase &W,
927  const value_type c, const BlockVectorBase &X);
928 
932  BlockVectorBase &operator *= (const value_type factor);
933 
937  BlockVectorBase &operator /= (const value_type factor);
938 
943  template <class BlockVector2>
944  void scale (const BlockVector2 &v);
945 
949  template <class BlockVector2>
950  void equ (const value_type a, const BlockVector2 &V);
951 
955  void equ (const value_type a, const BlockVectorBase &V,
956  const value_type b, const BlockVectorBase &W);
957 
962  void update_ghost_values () const;
963 
968  std::size_t memory_consumption () const;
969 
970 protected:
974  std::vector<VectorType> components;
975 
981 
985  template <typename N, bool C>
986  friend class ::internal::BlockVectorIterators::Iterator;
987 
988  template <typename> friend class BlockVectorBase;
989 };
990 
991 
994 /*----------------------- Inline functions ----------------------------------*/
995 
996 
997 #ifndef DOXYGEN
998 namespace internal
999 {
1000  namespace BlockVectorIterators
1001  {
1002  template <class BlockVectorType, bool Constness>
1003  inline
1005  Iterator (const Iterator<BlockVectorType,Constness> &c)
1006  :
1007  parent (c.parent),
1008  global_index (c.global_index),
1009  current_block (c.current_block),
1010  index_within_block (c.index_within_block),
1011  next_break_forward (c.next_break_forward),
1012  next_break_backward (c.next_break_backward)
1013  {}
1014 
1015 
1016 
1017  template <class BlockVectorType, bool Constness>
1018  inline
1019  Iterator<BlockVectorType,Constness>::
1020  Iterator (const Iterator<BlockVectorType,!Constness> &c)
1021  :
1022  parent (c.parent),
1024  current_block (c.current_block),
1025  index_within_block (c.index_within_block),
1026  next_break_forward (c.next_break_forward),
1027  next_break_backward (c.next_break_backward)
1028  {
1029  // Only permit copy-constructing const iterators from non-const
1030  // iterators, and not vice versa (i.e., Constness must always be
1031  // true). As mentioned above, try to check this at compile time if C++11
1032  // support is available.
1033  static_assert(Constness == true,
1034  "Constructing a non-const iterator from a const iterator "
1035  "does not make sense.");
1036  }
1037 
1038 
1039 
1040  template <class BlockVectorType, bool Constness>
1041  inline
1042  Iterator<BlockVectorType,Constness>::
1043  Iterator (BlockVector &parent,
1044  const size_type global_index,
1045  const size_type current_block,
1046  const size_type index_within_block,
1047  const size_type next_break_forward,
1048  const size_type next_break_backward)
1049  :
1050  parent (&parent),
1051  global_index (global_index),
1052  current_block (current_block),
1053  index_within_block (index_within_block),
1054  next_break_forward (next_break_forward),
1055  next_break_backward (next_break_backward)
1056  {
1057  }
1058 
1059 
1060 
1061  template <class BlockVectorType, bool Constness>
1062  inline
1063  Iterator<BlockVectorType,Constness> &
1064  Iterator<BlockVectorType,Constness>::
1065  operator = (const Iterator &c)
1066  {
1067  parent = c.parent;
1068  global_index = c.global_index;
1069  index_within_block = c.index_within_block;
1070  current_block = c.current_block;
1071  next_break_forward = c.next_break_forward;
1072  next_break_backward = c.next_break_backward;
1073 
1074  return *this;
1075  }
1076 
1077 
1078 
1079  template <class BlockVectorType, bool Constness>
1080  inline
1081  typename Iterator<BlockVectorType,Constness>::dereference_type
1082  Iterator<BlockVectorType,Constness>::operator * () const
1083  {
1084  return parent->block(current_block)(index_within_block);
1085  }
1086 
1087 
1088 
1089  template <class BlockVectorType, bool Constness>
1090  inline
1091  typename Iterator<BlockVectorType,Constness>::dereference_type
1092  Iterator<BlockVectorType,Constness>::operator [] (const difference_type d) const
1093  {
1094  // if the index pointed to is
1095  // still within the block we
1096  // currently point into, then we
1097  // can save the computation of
1098  // the block
1099  if ((global_index+d >= next_break_backward) &&
1100  (global_index+d <= next_break_forward))
1101  return parent->block(current_block)(index_within_block + d);
1102 
1103  // if the index is not within the
1104  // block of the block vector into
1105  // which we presently point, then
1106  // there is no way: we have to
1107  // search for the block. this can
1108  // be done through the parent
1109  // class as well.
1110  return (*parent)(global_index+d);
1111  }
1112 
1113 
1114 
1115  template <class BlockVectorType, bool Constness>
1116  inline
1117  Iterator<BlockVectorType,Constness> &
1118  Iterator<BlockVectorType,Constness>::operator ++ ()
1119  {
1120  move_forward ();
1121  return *this;
1122  }
1123 
1124 
1125 
1126  template <class BlockVectorType, bool Constness>
1127  inline
1128  Iterator<BlockVectorType,Constness>
1129  Iterator<BlockVectorType,Constness>::operator ++ (int)
1130  {
1131  const Iterator old_value = *this;
1132  move_forward ();
1133  return old_value;
1134  }
1135 
1136 
1137 
1138  template <class BlockVectorType, bool Constness>
1139  inline
1140  Iterator<BlockVectorType,Constness> &
1141  Iterator<BlockVectorType,Constness>::operator -- ()
1142  {
1143  move_backward ();
1144  return *this;
1145  }
1146 
1147 
1148 
1149  template <class BlockVectorType, bool Constness>
1150  inline
1151  Iterator<BlockVectorType,Constness>
1152  Iterator<BlockVectorType,Constness>::operator -- (int)
1153  {
1154  const Iterator old_value = *this;
1155  move_backward ();
1156  return old_value;
1157  }
1158 
1159 
1160 
1161  template <class BlockVectorType, bool Constness>
1162  template <bool OtherConstness>
1163  inline
1164  bool
1165  Iterator<BlockVectorType,Constness>::
1166  operator == (const Iterator<BlockVectorType, OtherConstness> &i) const
1167  {
1168  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1169 
1170  return (global_index == i.global_index);
1171  }
1172 
1173 
1174 
1175  template <class BlockVectorType, bool Constness>
1176  template <bool OtherConstness>
1177  inline
1178  bool
1179  Iterator<BlockVectorType,Constness>::
1180  operator != (const Iterator<BlockVectorType, OtherConstness> &i) const
1181  {
1182  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1183 
1184  return (global_index != i.global_index);
1185  }
1186 
1187 
1188 
1189  template <class BlockVectorType, bool Constness>
1190  template <bool OtherConstness>
1191  inline
1192  bool
1193  Iterator<BlockVectorType,Constness>::
1194  operator < (const Iterator<BlockVectorType, OtherConstness> &i) const
1195  {
1196  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1197 
1198  return (global_index < i.global_index);
1199  }
1200 
1201 
1202 
1203  template <class BlockVectorType, bool Constness>
1204  template <bool OtherConstness>
1205  inline
1206  bool
1207  Iterator<BlockVectorType,Constness>::
1208  operator <= (const Iterator<BlockVectorType, OtherConstness> &i) const
1209  {
1210  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1211 
1212  return (global_index <= i.global_index);
1213  }
1214 
1215 
1216 
1217  template <class BlockVectorType, bool Constness>
1218  template <bool OtherConstness>
1219  inline
1220  bool
1221  Iterator<BlockVectorType,Constness>::
1222  operator > (const Iterator<BlockVectorType, OtherConstness> &i) const
1223  {
1224  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1225 
1226  return (global_index > i.global_index);
1227  }
1228 
1229 
1230 
1231  template <class BlockVectorType, bool Constness>
1232  template <bool OtherConstness>
1233  inline
1234  bool
1235  Iterator<BlockVectorType,Constness>::
1236  operator >= (const Iterator<BlockVectorType, OtherConstness> &i) const
1237  {
1238  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1239 
1240  return (global_index >= i.global_index);
1241  }
1242 
1243 
1244 
1245  template <class BlockVectorType, bool Constness>
1246  template <bool OtherConstness>
1247  inline
1248  typename Iterator<BlockVectorType,Constness>::difference_type
1249  Iterator<BlockVectorType,Constness>::
1250  operator - (const Iterator<BlockVectorType, OtherConstness> &i) const
1251  {
1252  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1253 
1254  return (static_cast<signed int>(global_index) -
1255  static_cast<signed int>(i.global_index));
1256  }
1257 
1258 
1259 
1260  template <class BlockVectorType, bool Constness>
1261  inline
1262  Iterator<BlockVectorType,Constness>
1263  Iterator<BlockVectorType,Constness>::
1264  operator + (const difference_type &d) const
1265  {
1266  // if the index pointed to is
1267  // still within the block we
1268  // currently point into, then we
1269  // can save the computation of
1270  // the block
1271  if ((global_index+d >= next_break_backward) &&
1272  (global_index+d <= next_break_forward))
1273  return Iterator (*parent, global_index+d, current_block,
1274  index_within_block+d,
1275  next_break_forward, next_break_backward);
1276  else
1277  // outside present block, so
1278  // have to seek new block
1279  // anyway
1280  return Iterator (*parent, global_index+d);
1281  }
1282 
1283 
1284 
1285  template <class BlockVectorType, bool Constness>
1286  inline
1287  Iterator<BlockVectorType,Constness>
1288  Iterator<BlockVectorType,Constness>::
1289  operator - (const difference_type &d) const
1290  {
1291  // if the index pointed to is
1292  // still within the block we
1293  // currently point into, then we
1294  // can save the computation of
1295  // the block
1296  if ((global_index-d >= next_break_backward) &&
1297  (global_index-d <= next_break_forward))
1298  return Iterator (*parent, global_index-d, current_block,
1299  index_within_block-d,
1300  next_break_forward, next_break_backward);
1301  else
1302  // outside present block, so
1303  // have to seek new block
1304  // anyway
1305  return Iterator (*parent, global_index-d);
1306  }
1307 
1308 
1309 
1310  template <class BlockVectorType, bool Constness>
1311  inline
1312  Iterator<BlockVectorType,Constness> &
1313  Iterator<BlockVectorType,Constness>::
1314  operator += (const difference_type &d)
1315  {
1316  // if the index pointed to is
1317  // still within the block we
1318  // currently point into, then we
1319  // can save the computation of
1320  // the block
1321  if ((global_index+d >= next_break_backward) &&
1322  (global_index+d <= next_break_forward))
1323  {
1324  global_index += d;
1325  index_within_block += d;
1326  }
1327  else
1328  // outside present block, so
1329  // have to seek new block
1330  // anyway
1331  *this = Iterator (*parent, global_index+d);
1332 
1333  return *this;
1334  }
1335 
1336 
1337 
1338  template <class BlockVectorType, bool Constness>
1339  inline
1340  Iterator<BlockVectorType,Constness> &
1341  Iterator<BlockVectorType,Constness>::
1342  operator -= (const difference_type &d)
1343  {
1344  // if the index pointed to is
1345  // still within the block we
1346  // currently point into, then we
1347  // can save the computation of
1348  // the block
1349  if ((global_index-d >= next_break_backward) &&
1350  (global_index-d <= next_break_forward))
1351  {
1352  global_index -= d;
1353  index_within_block -= d;
1354  }
1355  else
1356  // outside present block, so
1357  // have to seek new block
1358  // anyway
1359  *this = Iterator (*parent, global_index-d);
1360 
1361  return *this;
1362  }
1363 
1364 
1365  template <class BlockVectorType, bool Constness>
1366  Iterator<BlockVectorType,Constness>::
1367  Iterator (BlockVector &parent,
1368  const size_type global_index)
1369  :
1370  parent (&parent),
1371  global_index (global_index)
1372  {
1373  // find which block we are
1374  // in. for this, take into
1375  // account that it happens at
1376  // times that people want to
1377  // initialize iterators
1378  // past-the-end
1379  if (global_index < parent.size())
1380  {
1381  const std::pair<size_type, size_type>
1382  indices = parent.block_indices.global_to_local(global_index);
1383  current_block = indices.first;
1384  index_within_block = indices.second;
1385 
1386  next_break_backward
1387  = parent.block_indices.local_to_global (current_block, 0);
1388  next_break_forward
1389  = (parent.block_indices.local_to_global (current_block, 0)
1390  +parent.block_indices.block_size(current_block)-1);
1391  }
1392  else
1393  // past the end. only have one
1394  // value for this
1395  {
1396  this->global_index = parent.size ();
1397  current_block = parent.n_blocks();
1398  index_within_block = 0;
1399  next_break_backward = global_index;
1400  next_break_forward = numbers::invalid_size_type;
1401  };
1402  }
1403 
1404 
1405 
1406  template <class BlockVectorType, bool Constness>
1407  void
1408  Iterator<BlockVectorType,Constness>::move_forward ()
1409  {
1410  if (global_index != next_break_forward)
1411  ++index_within_block;
1412  else
1413  {
1414  // ok, we traverse a boundary
1415  // between blocks:
1416  index_within_block = 0;
1417  ++current_block;
1418 
1419  // break backwards is now old
1420  // break forward
1421  next_break_backward = next_break_forward+1;
1422 
1423  // compute new break forward
1424  if (current_block < parent->block_indices.size())
1425  next_break_forward
1426  += parent->block_indices.block_size(current_block);
1427  else
1428  // if we are beyond the end,
1429  // then move the next
1430  // boundary arbitrarily far
1431  // away
1432  next_break_forward = numbers::invalid_size_type;
1433  };
1434 
1435  ++global_index;
1436  }
1437 
1438 
1439 
1440  template <class BlockVectorType, bool Constness>
1441  void
1442  Iterator<BlockVectorType,Constness>::move_backward ()
1443  {
1444  if (global_index != next_break_backward)
1445  --index_within_block;
1446  else if (current_block != 0)
1447  {
1448  // ok, we traverse a boundary
1449  // between blocks:
1450  --current_block;
1451  index_within_block = parent->block_indices.block_size(current_block)-1;
1452 
1453  // break forwards is now old
1454  // break backward
1455  next_break_forward = next_break_backward-1;
1456 
1457  // compute new break forward
1458  next_break_backward
1459  -= parent->block_indices.block_size (current_block);
1460  }
1461  else
1462  // current block was 0, we now
1463  // get into unspecified terrain
1464  {
1465  --current_block;
1466  index_within_block = numbers::invalid_size_type;
1467  next_break_forward = 0;
1468  next_break_backward = 0;
1469  };
1470 
1471  --global_index;
1472  }
1473 
1474 
1475  } // namespace BlockVectorIterators
1476 
1477 } //namespace internal
1478 
1479 
1480 
1481 template <class VectorType>
1482 inline
1483 std::size_t
1485 {
1486  return block_indices.total_size();
1487 }
1488 
1489 
1490 
1491 template <class VectorType>
1492 inline
1493 IndexSet
1495 {
1496  IndexSet is (size());
1497 
1498  // copy index sets from blocks into the global one, shifted
1499  // by the appropriate amount for each block
1500  for (unsigned int b=0; b<n_blocks(); ++b)
1501  {
1502  IndexSet x = block(b).locally_owned_elements();
1503  is.add_indices(x, block_indices.block_start(b));
1504  }
1505 
1506  is.compress();
1507 
1508  return is;
1509 }
1510 
1511 
1512 
1513 template <class VectorType>
1514 inline
1515 unsigned int
1517 {
1518  return block_indices.size();
1519 }
1520 
1521 
1522 template <class VectorType>
1523 inline
1525 BlockVectorBase<VectorType>::block (const unsigned int i)
1526 {
1527  Assert(i<n_blocks(), ExcIndexRange(i,0,n_blocks()));
1528 
1529  return components[i];
1530 }
1531 
1532 
1533 
1534 template <class VectorType>
1535 inline
1537 BlockVectorBase<VectorType>::block (const unsigned int i) const
1538 {
1539  Assert(i<n_blocks(), ExcIndexRange(i,0,n_blocks()));
1540 
1541  return components[i];
1542 }
1543 
1544 
1545 
1546 template <class VectorType>
1547 inline
1548 const BlockIndices &
1550 {
1551  return block_indices;
1552 }
1553 
1554 
1555 template <class VectorType>
1556 inline
1557 void
1559 {
1560  std::vector<size_type> sizes (n_blocks());
1561 
1562  for (size_type i=0; i<n_blocks(); ++i)
1563  sizes[i] = block(i).size();
1564 
1565  block_indices.reinit(sizes);
1566 }
1567 
1568 
1569 
1570 template <class VectorType>
1571 inline
1572 void
1574 {
1575  for (unsigned int i=0; i<n_blocks(); ++i)
1576  block(i).compress (operation);
1577 }
1578 
1579 
1580 
1581 template <class VectorType>
1582 inline
1585 {
1586  return iterator(*this, 0U);
1587 }
1588 
1589 
1590 
1591 template <class VectorType>
1592 inline
1595 {
1596  return const_iterator(*this, 0U);
1597 }
1598 
1599 
1600 template <class VectorType>
1601 inline
1604 {
1605  return iterator(*this, size());
1606 }
1607 
1608 
1609 
1610 template <class VectorType>
1611 inline
1614 {
1615  return const_iterator(*this, size());
1616 }
1617 
1618 
1619 template <class VectorType>
1620 inline
1621 bool
1623 (const size_type global_index) const
1624 {
1625  const std::pair<size_type,size_type> local_index
1626  = block_indices.global_to_local (global_index);
1627 
1628  return components[local_index.first].in_local_range (global_index);
1629 }
1630 
1631 
1632 template <class VectorType>
1633 bool
1635 {
1636  for (size_type i=0; i<n_blocks(); ++i)
1637  if (components[i].all_zero() == false)
1638  return false;
1639 
1640  return true;
1641 }
1642 
1643 
1644 
1645 template <class VectorType>
1646 bool
1648 {
1649  for (size_type i=0; i<n_blocks(); ++i)
1650  if (components[i].is_non_negative() == false)
1651  return false;
1652 
1653  return true;
1654 }
1655 
1656 
1657 
1658 template <class VectorType>
1659 typename BlockVectorBase<VectorType>::value_type
1662 {
1663  Assert (n_blocks() == v.n_blocks(),
1664  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1665 
1666  value_type sum = 0.;
1667  for (size_type i=0; i<n_blocks(); ++i)
1668  sum += components[i]*v.components[i];
1669 
1670  return sum;
1671 }
1672 
1673 
1674 template <class VectorType>
1677 {
1678  real_type sum = 0.;
1679  for (size_type i=0; i<n_blocks(); ++i)
1680  sum += components[i].norm_sqr();
1681 
1682  return sum;
1683 }
1684 
1685 
1686 
1687 template <class VectorType>
1688 typename BlockVectorBase<VectorType>::value_type
1690 {
1691  value_type sum = 0.;
1692  // need to do static_cast as otherwise it won't work with value_type=complex<T>
1693  for (size_type i=0; i<n_blocks(); ++i)
1694  sum += components[i].mean_value() * static_cast<double>(components[i].size());
1695 
1696  return sum/static_cast<double>(size());
1697 }
1698 
1699 
1700 
1701 template <class VectorType>
1704 {
1705  real_type sum = 0.;
1706  for (size_type i=0; i<n_blocks(); ++i)
1707  sum += components[i].l1_norm();
1708 
1709  return sum;
1710 }
1711 
1712 
1713 
1714 template <class VectorType>
1717 {
1718  return std::sqrt(norm_sqr());
1719 }
1720 
1721 
1722 
1723 template <class VectorType>
1726 {
1727  real_type sum = 0.;
1728  for (size_type i=0; i<n_blocks(); ++i)
1729  {
1730  value_type newval = components[i].linfty_norm();
1731  if (sum<newval)
1732  sum = newval;
1733  }
1734  return sum;
1735 }
1736 
1737 
1738 
1739 template <class VectorType>
1740 typename BlockVectorBase<VectorType>::value_type
1742 add_and_dot (const typename BlockVectorBase<VectorType>::value_type a,
1743  const BlockVectorBase<VectorType> &V,
1744  const BlockVectorBase<VectorType> &W)
1745 {
1746  AssertDimension (n_blocks(), V.n_blocks());
1747  AssertDimension (n_blocks(), W.n_blocks());
1748 
1749  value_type sum = 0.;
1750  for (size_type i=0; i<n_blocks(); ++i)
1751  sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
1752 
1753  return sum;
1754 }
1755 
1756 
1757 
1758 template <class VectorType>
1761 {
1762  Assert (n_blocks() == v.n_blocks(),
1763  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1764 
1765  for (size_type i=0; i<n_blocks(); ++i)
1766  {
1767  components[i] += v.components[i];
1768  }
1769 
1770  return *this;
1771 }
1772 
1773 
1774 
1775 template <class VectorType>
1778 {
1779  Assert (n_blocks() == v.n_blocks(),
1780  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1781 
1782  for (size_type i=0; i<n_blocks(); ++i)
1783  {
1784  components[i] -= v.components[i];
1785  }
1786  return *this;
1787 }
1788 
1789 
1790 
1791 template <class VectorType>
1792 template <typename Number>
1793 inline
1794 void
1795 BlockVectorBase<VectorType>::add (const std::vector<size_type> &indices,
1796  const std::vector<Number> &values)
1797 {
1798  Assert (indices.size() == values.size(),
1799  ExcDimensionMismatch(indices.size(), values.size()));
1800  add (indices.size(), &indices[0], &values[0]);
1801 }
1802 
1803 
1804 
1805 template <class VectorType>
1806 template <typename Number>
1807 inline
1808 void
1809 BlockVectorBase<VectorType>::add (const std::vector<size_type> &indices,
1810  const Vector<Number> &values)
1811 {
1812  Assert (indices.size() == values.size(),
1813  ExcDimensionMismatch(indices.size(), values.size()));
1814  const size_type n_indices = indices.size();
1815  for (size_type i=0; i<n_indices; ++i)
1816  (*this)(indices[i]) += values(i);
1817 }
1818 
1819 
1820 
1821 template <class VectorType>
1822 template <typename Number>
1823 inline
1824 void
1825 BlockVectorBase<VectorType>::add (const size_type n_indices,
1826  const size_type *indices,
1827  const Number *values)
1828 {
1829  for (size_type i=0; i<n_indices; ++i)
1830  (*this)(indices[i]) += values[i];
1831 }
1832 
1833 
1834 
1835 template <class VectorType>
1836 void BlockVectorBase<VectorType>::add (const value_type a)
1837 {
1838  AssertIsFinite(a);
1839 
1840  for (size_type i=0; i<n_blocks(); ++i)
1841  {
1842  components[i].add(a);
1843  }
1844 }
1845 
1846 
1847 
1848 template <class VectorType>
1849 void BlockVectorBase<VectorType>::add (const value_type a,
1850  const BlockVectorBase<VectorType> &v)
1851 {
1852 
1853  AssertIsFinite(a);
1854 
1855  Assert (n_blocks() == v.n_blocks(),
1856  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1857 
1858  for (size_type i=0; i<n_blocks(); ++i)
1859  {
1860  components[i].add(a, v.components[i]);
1861  }
1862 }
1863 
1864 
1865 
1866 template <class VectorType>
1867 void BlockVectorBase<VectorType>::add (const value_type a,
1868  const BlockVectorBase<VectorType> &v,
1869  const value_type b,
1870  const BlockVectorBase<VectorType> &w)
1871 {
1872 
1873  AssertIsFinite(a);
1874  AssertIsFinite(b);
1875 
1876  Assert (n_blocks() == v.n_blocks(),
1877  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1878  Assert (n_blocks() == w.n_blocks(),
1879  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1880 
1881 
1882  for (size_type i=0; i<n_blocks(); ++i)
1883  {
1884  components[i].add(a, v.components[i], b, w.components[i]);
1885  }
1886 }
1887 
1888 
1889 
1890 template <class VectorType>
1891 void BlockVectorBase<VectorType>::sadd (const value_type x,
1892  const BlockVectorBase<VectorType> &v)
1893 {
1894 
1895  AssertIsFinite(x);
1896 
1897  Assert (n_blocks() == v.n_blocks(),
1898  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1899 
1900  for (size_type i=0; i<n_blocks(); ++i)
1901  {
1902  components[i].sadd(x, v.components[i]);
1903  }
1904 }
1905 
1906 
1907 
1908 template <class VectorType>
1909 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
1910  const BlockVectorBase<VectorType> &v)
1911 {
1912 
1913  AssertIsFinite(x);
1914  AssertIsFinite(a);
1915 
1916  Assert (n_blocks() == v.n_blocks(),
1917  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1918 
1919  for (size_type i=0; i<n_blocks(); ++i)
1920  {
1921  components[i].sadd(x, a, v.components[i]);
1922  }
1923 }
1924 
1925 
1926 
1927 template <class VectorType>
1928 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
1929  const BlockVectorBase<VectorType> &v,
1930  const value_type b,
1931  const BlockVectorBase<VectorType> &w)
1932 {
1933 
1934  AssertIsFinite(x);
1935  AssertIsFinite(a);
1936  AssertIsFinite(b);
1937 
1938  Assert (n_blocks() == v.n_blocks(),
1939  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1940  Assert (n_blocks() == w.n_blocks(),
1941  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1942 
1943  for (size_type i=0; i<n_blocks(); ++i)
1944  {
1945  components[i].sadd(x, a, v.components[i], b, w.components[i]);
1946  }
1947 }
1948 
1949 
1950 
1951 template <class VectorType>
1952 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
1953  const BlockVectorBase<VectorType> &v,
1954  const value_type b,
1955  const BlockVectorBase<VectorType> &w,
1956  const value_type c,
1957  const BlockVectorBase<VectorType> &y)
1958 {
1959 
1960  AssertIsFinite(x);
1961  AssertIsFinite(a);
1962  AssertIsFinite(b);
1963  AssertIsFinite(c);
1964 
1965  Assert (n_blocks() == v.n_blocks(),
1966  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1967  Assert (n_blocks() == w.n_blocks(),
1968  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
1969  Assert (n_blocks() == y.n_blocks(),
1970  ExcDimensionMismatch(n_blocks(), y.n_blocks()));
1971 
1972  for (size_type i=0; i<n_blocks(); ++i)
1973  {
1974  components[i].sadd(x, a, v.components[i],
1975  b, w.components[i], c, y.components[i]);
1976  }
1977 }
1978 
1979 
1980 
1981 template <class VectorType>
1982 template <class BlockVector2>
1983 void BlockVectorBase<VectorType>::scale (const BlockVector2 &v)
1984 {
1985  Assert (n_blocks() == v.n_blocks(),
1986  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
1987  for (size_type i=0; i<n_blocks(); ++i)
1988  components[i].scale(v.block(i));
1989 }
1990 
1991 
1992 
1993 template <class VectorType>
1994 void BlockVectorBase<VectorType>::equ (const value_type a,
1995  const BlockVectorBase<VectorType> &v,
1996  const value_type b,
1997  const BlockVectorBase<VectorType> &w)
1998 {
1999 
2000  AssertIsFinite(a);
2001  AssertIsFinite(b);
2002 
2003  Assert (n_blocks() == v.n_blocks(),
2004  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
2005  Assert (n_blocks() == w.n_blocks(),
2006  ExcDimensionMismatch(n_blocks(), w.n_blocks()));
2007 
2008  for (size_type i=0; i<n_blocks(); ++i)
2009  {
2010  components[i].equ( a, v.components[i], b, w.components[i]);
2011  }
2012 }
2013 
2014 
2015 
2016 template <class VectorType>
2017 std::size_t
2019 {
2020  return (MemoryConsumption::memory_consumption (this->block_indices)
2021  +
2022  MemoryConsumption::memory_consumption (this->components));
2023 }
2024 
2025 
2026 
2027 template <class VectorType>
2028 template <class BlockVector2>
2029 void BlockVectorBase<VectorType>::equ (const value_type a,
2030  const BlockVector2 &v)
2031 {
2032 
2033  AssertIsFinite(a);
2034 
2035  Assert (n_blocks() == v.n_blocks(),
2036  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
2037 
2038  for (size_type i=0; i<n_blocks(); ++i)
2039  components[i].equ( a, v.components[i]);
2040 }
2041 
2042 
2043 
2044 template <class VectorType>
2046 {
2047  for (size_type i=0; i<n_blocks(); ++i)
2048  block(i).update_ghost_values ();
2049 }
2050 
2051 
2052 
2053 template <class VectorType>
2055 BlockVectorBase<VectorType>::operator = (const value_type s)
2056 {
2057 
2058  AssertIsFinite(s);
2059 
2060  for (size_type i=0; i<n_blocks(); ++i)
2061  components[i] = s;
2062 
2063  return *this;
2064 }
2065 
2066 
2067 template <class VectorType>
2070 {
2071  AssertDimension(n_blocks(), v.n_blocks());
2072 
2073  for (size_type i=0; i<n_blocks(); ++i)
2074  components[i] = v.components[i];
2075 
2076  return *this;
2077 }
2078 
2079 
2080 template <class VectorType>
2081 template <class VectorType2>
2084 {
2085  AssertDimension(n_blocks(), v.n_blocks());
2086 
2087  for (size_type i=0; i<n_blocks(); ++i)
2088  components[i] = v.components[i];
2089 
2090  return *this;
2091 }
2092 
2093 
2094 
2095 template <class VectorType>
2097 BlockVectorBase<VectorType>::operator = (const VectorType &v)
2098 {
2099  Assert (size() == v.size(),
2100  ExcDimensionMismatch(size(), v.size()));
2101 
2102  size_type index_v = 0;
2103  for (size_type b=0; b<n_blocks(); ++b)
2104  for (size_type i=0; i<block(b).size(); ++i, ++index_v)
2105  block(b)(i) = v(index_v);
2106 
2107  return *this;
2108 }
2109 
2110 
2111 
2112 template <class VectorType>
2113 template <class VectorType2>
2114 inline
2115 bool
2118 {
2119  Assert (block_indices == v.block_indices, ExcDifferentBlockIndices());
2120 
2121  for (size_type i=0; i<n_blocks(); ++i)
2122  if ( ! (components[i] == v.components[i]))
2123  return false;
2124 
2125  return true;
2126 }
2127 
2128 
2129 
2130 template <class VectorType>
2131 inline
2133 BlockVectorBase<VectorType>::operator *= (const value_type factor)
2134 {
2135 
2136  AssertIsFinite(factor);
2137 
2138  for (size_type i=0; i<n_blocks(); ++i)
2139  components[i] *= factor;
2140 
2141  return *this;
2142 }
2143 
2144 
2145 
2146 template <class VectorType>
2147 inline
2149 BlockVectorBase<VectorType>::operator /= (const value_type factor)
2150 {
2151 
2152  AssertIsFinite(factor);
2153  Assert (factor != 0., ExcDivideByZero() );
2154 
2155  for (size_type i=0; i<n_blocks(); ++i)
2156  components[i] /= factor;
2157 
2158  return *this;
2159 }
2160 
2161 
2162 template <class VectorType>
2163 inline
2164 typename BlockVectorBase<VectorType>::value_type
2165 BlockVectorBase<VectorType>::operator() (const size_type i) const
2166 {
2167  const std::pair<unsigned int,size_type> local_index
2168  = block_indices.global_to_local (i);
2169  return components[local_index.first](local_index.second);
2170 }
2171 
2172 
2173 
2174 template <class VectorType>
2175 inline
2176 typename BlockVectorBase<VectorType>::reference
2178 {
2179  const std::pair<unsigned int,size_type> local_index
2180  = block_indices.global_to_local (i);
2181  return components[local_index.first](local_index.second);
2182 }
2183 
2184 
2185 
2186 template <class VectorType>
2187 inline
2188 typename BlockVectorBase<VectorType>::value_type
2189 BlockVectorBase<VectorType>::operator[] (const size_type i) const
2190 {
2191  return operator()(i);
2192 }
2193 
2194 
2195 
2196 template <class VectorType>
2197 inline
2198 typename BlockVectorBase<VectorType>::reference
2200 {
2201  return operator()(i);
2202 }
2203 
2204 
2205 
2206 template <typename VectorType>
2207 template <typename OtherNumber>
2208 inline
2209 void BlockVectorBase<VectorType>::extract_subvector_to (const std::vector<size_type> &indices,
2210  std::vector<OtherNumber> &values) const
2211 {
2212  for (size_type i = 0; i < indices.size(); ++i)
2213  values[i] = operator()(indices[i]);
2214 }
2215 
2216 
2217 
2218 template <typename VectorType>
2219 template <typename ForwardIterator, typename OutputIterator>
2220 inline
2221 void BlockVectorBase<VectorType>::extract_subvector_to (ForwardIterator indices_begin,
2222  const ForwardIterator indices_end,
2223  OutputIterator values_begin) const
2224 {
2225  while (indices_begin != indices_end)
2226  {
2227  *values_begin = operator()(*indices_begin);
2228  indices_begin++;
2229  values_begin++;
2230  }
2231 }
2232 
2233 #endif // DOXYGEN
2234 
2235 DEAL_II_NAMESPACE_CLOSE
2236 
2237 #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:81
std::random_access_iterator_tag iterator_type
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1175
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:654
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:1211
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