Reference documentation for deal.II version 8.1.0
block_vector_base.h
1 // ---------------------------------------------------------------------
2 // @f$Id: block_vector_base.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2004 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__block_vector_base_h
18 #define __deal2__block_vector_base_h
19 
20 
21 #include <deal.II/base/config.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 <vector>
29 #include <iterator>
30 #include <cmath>
31 #include <cstddef>
32 
34 
35 
40 template <typename> class BlockVectorBase;
41 
42 
59 template <typename VectorType>
61 {
62 private:
63  struct yes_type
64  {
65  char c[1];
66  };
67  struct no_type
68  {
69  char c[2];
70  };
71 
77  template <typename T>
79 
85  static no_type check_for_block_vector (...);
86 
87 public:
95  static const bool value = (sizeof(check_for_block_vector
96  ((VectorType *)0))
97  ==
98  sizeof(yes_type));
99 };
100 
101 
102 // instantiation of the static member
103 template <typename VectorType>
105 
106 
107 
108 
109 namespace internal
110 {
111 
117  namespace BlockVectorIterators
118  {
126  template <class BlockVectorType, bool constness>
127  struct Types
128  {
129  };
130 
131 
132 
143  template <class BlockVectorType>
144  struct Types<BlockVectorType,false>
145  {
153  typedef typename BlockVectorType::BlockType Vector;
154 
161  typedef BlockVectorType BlockVector;
162 
169 
176  typedef typename Vector::reference dereference_type;
177  };
178 
179 
180 
191  template <class BlockVectorType>
192  struct Types<BlockVectorType,true>
193  {
202  typedef const typename BlockVectorType::BlockType Vector;
203 
211  typedef const BlockVectorType BlockVector;
212 
219  typedef const typename BlockVector::value_type value_type;
220 
230  typedef value_type dereference_type;
231  };
232 
233 
292  template <class BlockVectorType, bool constness>
293  class Iterator :
294  public std::iterator<std::random_access_iterator_tag,
295  typename Types<BlockVectorType,constness>::value_type>
296  {
297  private:
305 
306  public:
311 
321  typedef
324 
333  typedef std::random_access_iterator_tag iterator_type;
334  typedef std::ptrdiff_t difference_type;
335  typedef typename BlockVectorType::reference reference;
336  typedef value_type *pointer;
337 
338  typedef
340  dereference_type;
341 
349  typedef
352 
368  const size_type global_index);
369 
374 
385  Iterator (const InverseConstnessIterator &c);
386 
387  private:
396  const size_type global_index,
397  const size_type current_block,
398  const size_type index_within_block,
399  const size_type next_break_forward,
400  const size_type next_break_backward);
401 
402  public:
403 
407  Iterator &operator = (const Iterator &c);
408 
417  dereference_type operator * () const;
418 
425  dereference_type operator [] (const difference_type d) const;
426 
434 
443  Iterator operator ++ (int);
444 
452 
461  Iterator operator -- (int);
462 
471  bool operator == (const Iterator &i) const;
472 
478  bool operator == (const InverseConstnessIterator &i) const;
479 
488  bool operator != (const Iterator &i) const;
489 
495  bool operator != (const InverseConstnessIterator &i) const;
496 
508  bool operator < (const Iterator &i) const;
509 
515  bool operator < (const InverseConstnessIterator &i) const;
516 
521  bool operator <= (const Iterator &i) const;
522 
528  bool operator <= (const InverseConstnessIterator &i) const;
529 
534  bool operator > (const Iterator &i) const;
535 
541  bool operator > (const InverseConstnessIterator &i) const;
542 
547  bool operator >= (const Iterator &i) const;
548 
554  bool operator >= (const InverseConstnessIterator &i) const;
555 
561  difference_type operator - (const Iterator &i) const;
562 
567  difference_type operator - (const InverseConstnessIterator &i) const;
568 
574  Iterator operator + (const difference_type &d) const;
575 
581  Iterator operator - (const difference_type &d) const;
582 
588  Iterator &operator += (const difference_type &d);
589 
595  Iterator &operator -= (const difference_type &d);
596 
603  DeclException0 (ExcPointerToDifferentVectors);
607  DeclException0 (ExcCastingAwayConstness);
609  private:
621 
627  size_type global_index;
628 
635  unsigned int current_block;
636  size_type index_within_block;
637 
651  size_type next_break_backward;
652 
656  void move_forward ();
657 
661  void move_backward ();
662 
663 
673  template <typename N, bool C>
674  friend class Iterator;
675  };
676  } // namespace BlockVectorIterators
677 } // namespace internal
678 
679 
720 template <class VectorType>
721 class BlockVectorBase : public Subscriptor
722 {
723 public:
728  typedef VectorType BlockType;
729 
730  /*
731  * Declare standard types used in
732  * all containers. These types
733  * parallel those in the
734  * <tt>C++</tt> standard
735  * libraries
736  * <tt>std::vector<...></tt>
737  * class. This includes iterator
738  * types.
739  */
740  typedef typename BlockType::value_type value_type;
741  typedef value_type *pointer;
742  typedef const value_type *const_pointer;
743  typedef ::internal::BlockVectorIterators::Iterator<BlockVectorBase,false> iterator;
744  typedef ::internal::BlockVectorIterators::Iterator<BlockVectorBase,true> const_iterator;
745  typedef typename BlockType::reference reference;
746  typedef typename BlockType::const_reference const_reference;
747  typedef types::global_dof_index size_type;
748 
768  typedef typename BlockType::real_type real_type;
769 
781  static const bool supports_distributed_data = BlockType::supports_distributed_data;
782 
786  BlockVectorBase ();
787 
799  void collect_sizes ();
800 
814  void compress (::VectorOperation::values operation);
815 
820 
824  BlockType &
825  block (const unsigned int i);
826 
830  const BlockType &
831  block (const unsigned int i) const;
832 
843  const BlockIndices &
844  get_block_indices () const;
845 
849  unsigned int n_blocks () const;
850 
856  std::size_t size () const;
857 
877 
882  iterator begin ();
883 
889  const_iterator begin () const;
890 
895  iterator end ();
896 
902  const_iterator end () const;
903 
907  value_type operator() (const size_type i) const;
908 
913  reference operator() (const size_type i);
914 
920  value_type operator[] (const size_type i) const;
921 
928  reference operator[] (const size_type i);
929 
940  template <typename OtherNumber>
941  void extract_subvector_to (const std::vector<size_type> &indices,
942  std::vector<OtherNumber> &values) const;
943 
948  template <typename ForwardIterator, typename OutputIterator>
949  void extract_subvector_to (ForwardIterator indices_begin,
950  const ForwardIterator indices_end,
951  OutputIterator values_begin) const;
952 
958  BlockVectorBase &operator = (const value_type s);
959 
965  operator= (const BlockVectorBase &V);
966 
971  template <class VectorType2>
973  operator= (const BlockVectorBase<VectorType2> &V);
974 
980  operator = (const VectorType &v);
981 
988  template <class VectorType2>
989  bool
990  operator == (const BlockVectorBase<VectorType2> &v) const;
991 
995  value_type operator* (const BlockVectorBase &V) const;
996 
1000  real_type norm_sqr () const;
1001 
1006  value_type mean_value () const;
1007 
1012  real_type l1_norm () const;
1013 
1019  real_type l2_norm () const;
1020 
1026  real_type linfty_norm () const;
1027 
1033  bool in_local_range (const size_type global_index) const;
1034 
1043  bool all_zero () const;
1044 
1053  bool is_non_negative () const;
1054 
1059  BlockVectorBase &
1060  operator += (const BlockVectorBase &V);
1061 
1066  BlockVectorBase &
1067  operator -= (const BlockVectorBase &V);
1068 
1069 
1078  template <typename Number>
1079  void add (const std::vector<size_type> &indices,
1080  const std::vector<Number> &values);
1081 
1089  template <typename Number>
1090  void add (const std::vector<size_type> &indices,
1091  const Vector<Number> &values);
1092 
1102  template <typename Number>
1103  void add (const size_type n_elements,
1104  const size_type *indices,
1105  const Number *values);
1106 
1113  void add (const value_type s);
1114 
1120  void add (const BlockVectorBase &V);
1121 
1126  void add (const value_type a, const BlockVectorBase &V);
1127 
1132  void add (const value_type a, const BlockVectorBase &V,
1133  const value_type b, const BlockVectorBase &W);
1134 
1139  void sadd (const value_type s, const BlockVectorBase &V);
1140 
1145  void sadd (const value_type s, const value_type a, const BlockVectorBase &V);
1146 
1151  void sadd (const value_type s, const value_type a,
1152  const BlockVectorBase &V,
1153  const value_type b, const BlockVectorBase &W);
1154 
1159  void sadd (const value_type s, const value_type a,
1160  const BlockVectorBase &V,
1161  const value_type b, const BlockVectorBase &W,
1162  const value_type c, const BlockVectorBase &X);
1163 
1169  BlockVectorBase &operator *= (const value_type factor);
1170 
1176  BlockVectorBase &operator /= (const value_type factor);
1177 
1183  template <class BlockVector2>
1184  void scale (const BlockVector2 &v);
1185 
1189  template <class BlockVector2>
1190  void equ (const value_type a, const BlockVector2 &V);
1191 
1196  void equ (const value_type a, const BlockVectorBase &V,
1197  const value_type b, const BlockVectorBase &W);
1198 
1214  void update_ghost_values () const;
1215 
1221  std::size_t memory_consumption () const;
1222 
1223 protected:
1227  std::vector<VectorType> components;
1228 
1236 
1241  template <typename N, bool C>
1242  friend class ::internal::BlockVectorIterators::Iterator;
1243 
1244  template <typename> friend class BlockVectorBase;
1245 };
1246 
1247 
1250 /*----------------------- Inline functions ----------------------------------*/
1251 
1252 
1253 #ifndef DOXYGEN
1254 namespace internal
1255 {
1256  namespace BlockVectorIterators
1257  {
1258 
1259  template <class BlockVectorType, bool constness>
1260  inline
1261  Iterator<BlockVectorType,constness>::
1262  Iterator (const Iterator<BlockVectorType,constness> &c)
1263  :
1264  parent (c.parent),
1265  global_index (c.global_index),
1266  current_block (c.current_block),
1267  index_within_block (c.index_within_block),
1268  next_break_forward (c.next_break_forward),
1269  next_break_backward (c.next_break_backward)
1270  {}
1271 
1272 
1273 
1274  template <class BlockVectorType, bool constness>
1275  inline
1276  Iterator<BlockVectorType,constness>::
1277  Iterator (const InverseConstnessIterator &c)
1278  :
1279  parent (const_cast<BlockVectorType *>(c.parent)),
1280  global_index (c.global_index),
1281  current_block (c.current_block),
1282  index_within_block (c.index_within_block),
1283  next_break_forward (c.next_break_forward),
1284  next_break_backward (c.next_break_backward)
1285  {
1286  // if constness==false, then the
1287  // constness of the iterator we
1288  // got is true and we are trying
1289  // to cast away the
1290  // constness. disallow this
1291  Assert (constness==true, ExcCastingAwayConstness());
1292  }
1293 
1294 
1295 
1296  template <class BlockVectorType, bool constness>
1297  inline
1298  Iterator<BlockVectorType,constness>::
1299  Iterator (BlockVector &parent,
1300  const size_type global_index,
1301  const size_type current_block,
1302  const size_type index_within_block,
1303  const size_type next_break_forward,
1304  const size_type next_break_backward)
1305  :
1306  parent (&parent),
1307  global_index (global_index),
1308  current_block (current_block),
1309  index_within_block (index_within_block),
1310  next_break_forward (next_break_forward),
1311  next_break_backward (next_break_backward)
1312  {
1313  }
1314 
1315 
1316 
1317  template <class BlockVectorType, bool constness>
1318  inline
1319  Iterator<BlockVectorType,constness> &
1320  Iterator<BlockVectorType,constness>::
1321  operator = (const Iterator &c)
1322  {
1323  parent = c.parent;
1324  global_index = c.global_index;
1325  index_within_block = c.index_within_block;
1326  current_block = c.current_block;
1327  next_break_forward = c.next_break_forward;
1328  next_break_backward = c.next_break_backward;
1329 
1330  return *this;
1331  }
1332 
1333 
1334 
1335  template <class BlockVectorType, bool constness>
1336  inline
1337  typename Iterator<BlockVectorType,constness>::dereference_type
1338  Iterator<BlockVectorType,constness>::operator * () const
1339  {
1340  return parent->block(current_block)(index_within_block);
1341  }
1342 
1343 
1344 
1345  template <class BlockVectorType, bool constness>
1346  inline
1347  typename Iterator<BlockVectorType,constness>::dereference_type
1348  Iterator<BlockVectorType,constness>::operator [] (const difference_type d) const
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  return parent->block(current_block)(index_within_block + d);
1358 
1359  // if the index is not within the
1360  // block of the block vector into
1361  // which we presently point, then
1362  // there is no way: we have to
1363  // search for the block. this can
1364  // be done through the parent
1365  // class as well.
1366  return (*parent)(global_index+d);
1367  }
1368 
1369 
1370 
1371  template <class BlockVectorType, bool constness>
1372  inline
1373  Iterator<BlockVectorType,constness> &
1374  Iterator<BlockVectorType,constness>::operator ++ ()
1375  {
1376  move_forward ();
1377  return *this;
1378  }
1379 
1380 
1381 
1382  template <class BlockVectorType, bool constness>
1383  inline
1384  Iterator<BlockVectorType,constness>
1385  Iterator<BlockVectorType,constness>::operator ++ (int)
1386  {
1387  const Iterator old_value = *this;
1388  move_forward ();
1389  return old_value;
1390  }
1391 
1392 
1393 
1394  template <class BlockVectorType, bool constness>
1395  inline
1396  Iterator<BlockVectorType,constness> &
1397  Iterator<BlockVectorType,constness>::operator -- ()
1398  {
1399  move_backward ();
1400  return *this;
1401  }
1402 
1403 
1404 
1405  template <class BlockVectorType, bool constness>
1406  inline
1407  Iterator<BlockVectorType,constness>
1408  Iterator<BlockVectorType,constness>::operator -- (int)
1409  {
1410  const Iterator old_value = *this;
1411  move_backward ();
1412  return old_value;
1413  }
1414 
1415 
1416 
1417  template <class BlockVectorType, bool constness>
1418  inline
1419  bool
1420  Iterator<BlockVectorType,constness>::
1421  operator == (const Iterator &i) const
1422  {
1423  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1424 
1425  return (global_index == i.global_index);
1426  }
1427 
1428 
1429 
1430  template <class BlockVectorType, bool constness>
1431  inline
1432  bool
1433  Iterator<BlockVectorType,constness>::
1434  operator == (const InverseConstnessIterator &i) const
1435  {
1436  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1437 
1438  return (global_index == i.global_index);
1439  }
1440 
1441 
1442 
1443  template <class BlockVectorType, bool constness>
1444  inline
1445  bool
1446  Iterator<BlockVectorType,constness>::
1447  operator != (const Iterator &i) const
1448  {
1449  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1450 
1451  return (global_index != i.global_index);
1452  }
1453 
1454 
1455 
1456  template <class BlockVectorType, bool constness>
1457  inline
1458  bool
1459  Iterator<BlockVectorType,constness>::
1460  operator != (const InverseConstnessIterator &i) const
1461  {
1462  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1463 
1464  return (global_index != i.global_index);
1465  }
1466 
1467 
1468 
1469  template <class BlockVectorType, bool constness>
1470  inline
1471  bool
1472  Iterator<BlockVectorType,constness>::
1473  operator < (const Iterator &i) const
1474  {
1475  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1476 
1477  return (global_index < i.global_index);
1478  }
1479 
1480 
1481 
1482  template <class BlockVectorType, bool constness>
1483  inline
1484  bool
1485  Iterator<BlockVectorType,constness>::
1486  operator < (const InverseConstnessIterator &i) const
1487  {
1488  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1489 
1490  return (global_index < i.global_index);
1491  }
1492 
1493 
1494 
1495  template <class BlockVectorType, bool constness>
1496  inline
1497  bool
1498  Iterator<BlockVectorType,constness>::
1499  operator <= (const Iterator &i) const
1500  {
1501  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1502 
1503  return (global_index <= i.global_index);
1504  }
1505 
1506 
1507 
1508  template <class BlockVectorType, bool constness>
1509  inline
1510  bool
1511  Iterator<BlockVectorType,constness>::
1512  operator <= (const InverseConstnessIterator &i) const
1513  {
1514  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1515 
1516  return (global_index <= i.global_index);
1517  }
1518 
1519 
1520 
1521  template <class BlockVectorType, bool constness>
1522  inline
1523  bool
1524  Iterator<BlockVectorType,constness>::
1525  operator > (const Iterator &i) const
1526  {
1527  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1528 
1529  return (global_index > i.global_index);
1530  }
1531 
1532 
1533 
1534  template <class BlockVectorType, bool constness>
1535  inline
1536  bool
1537  Iterator<BlockVectorType,constness>::
1538  operator > (const InverseConstnessIterator &i) const
1539  {
1540  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1541 
1542  return (global_index > i.global_index);
1543  }
1544 
1545 
1546 
1547  template <class BlockVectorType, bool constness>
1548  inline
1549  bool
1550  Iterator<BlockVectorType,constness>::
1551  operator >= (const Iterator &i) const
1552  {
1553  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1554 
1555  return (global_index >= i.global_index);
1556  }
1557 
1558 
1559 
1560  template <class BlockVectorType, bool constness>
1561  inline
1562  bool
1563  Iterator<BlockVectorType,constness>::
1564  operator >= (const InverseConstnessIterator &i) const
1565  {
1566  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1567 
1568  return (global_index >= i.global_index);
1569  }
1570 
1571 
1572 
1573  template <class BlockVectorType, bool constness>
1574  inline
1575  typename Iterator<BlockVectorType,constness>::difference_type
1576  Iterator<BlockVectorType,constness>::
1577  operator - (const Iterator &i) const
1578  {
1579  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1580 
1581  return (static_cast<signed int>(global_index) -
1582  static_cast<signed int>(i.global_index));
1583  }
1584 
1585 
1586 
1587  template <class BlockVectorType, bool constness>
1588  inline
1589  typename Iterator<BlockVectorType,constness>::difference_type
1590  Iterator<BlockVectorType,constness>::
1591  operator - (const InverseConstnessIterator &i) const
1592  {
1593  Assert (parent == i.parent, ExcPointerToDifferentVectors());
1594 
1595  return (static_cast<signed int>(global_index) -
1596  static_cast<signed int>(i.global_index));
1597  }
1598 
1599 
1600 
1601  template <class BlockVectorType, bool constness>
1602  inline
1603  Iterator<BlockVectorType,constness>
1604  Iterator<BlockVectorType,constness>::
1605  operator + (const difference_type &d) const
1606  {
1607  // if the index pointed to is
1608  // still within the block we
1609  // currently point into, then we
1610  // can save the computation of
1611  // the block
1612  if ((global_index+d >= next_break_backward) &&
1613  (global_index+d <= next_break_forward))
1614  return Iterator (*parent, global_index+d, current_block,
1615  index_within_block+d,
1616  next_break_forward, next_break_backward);
1617  else
1618  // outside present block, so
1619  // have to seek new block
1620  // anyway
1621  return Iterator (*parent, global_index+d);
1622  }
1623 
1624 
1625 
1626  template <class BlockVectorType, bool constness>
1627  inline
1628  Iterator<BlockVectorType,constness>
1629  Iterator<BlockVectorType,constness>::
1630  operator - (const difference_type &d) const
1631  {
1632  // if the index pointed to is
1633  // still within the block we
1634  // currently point into, then we
1635  // can save the computation of
1636  // the block
1637  if ((global_index-d >= next_break_backward) &&
1638  (global_index-d <= next_break_forward))
1639  return Iterator (*parent, global_index-d, current_block,
1640  index_within_block-d,
1641  next_break_forward, next_break_backward);
1642  else
1643  // outside present block, so
1644  // have to seek new block
1645  // anyway
1646  return Iterator (*parent, global_index-d);
1647  }
1648 
1649 
1650 
1651  template <class BlockVectorType, bool constness>
1652  inline
1653  Iterator<BlockVectorType,constness> &
1654  Iterator<BlockVectorType,constness>::
1655  operator += (const difference_type &d)
1656  {
1657  // if the index pointed to is
1658  // still within the block we
1659  // currently point into, then we
1660  // can save the computation of
1661  // the block
1662  if ((global_index+d >= next_break_backward) &&
1663  (global_index+d <= next_break_forward))
1664  {
1665  global_index += d;
1666  index_within_block += d;
1667  }
1668  else
1669  // outside present block, so
1670  // have to seek new block
1671  // anyway
1672  *this = Iterator (*parent, global_index+d);
1673 
1674  return *this;
1675  }
1676 
1677 
1678 
1679  template <class BlockVectorType, bool constness>
1680  inline
1681  Iterator<BlockVectorType,constness> &
1682  Iterator<BlockVectorType,constness>::
1683  operator -= (const difference_type &d)
1684  {
1685  // if the index pointed to is
1686  // still within the block we
1687  // currently point into, then we
1688  // can save the computation of
1689  // the block
1690  if ((global_index-d >= next_break_backward) &&
1691  (global_index-d <= next_break_forward))
1692  {
1693  global_index -= d;
1694  index_within_block -= d;
1695  }
1696  else
1697  // outside present block, so
1698  // have to seek new block
1699  // anyway
1700  *this = Iterator (*parent, global_index-d);
1701 
1702  return *this;
1703  }
1704 
1705 
1706  template <class BlockVectorType, bool constness>
1707  Iterator<BlockVectorType,constness>::
1708  Iterator (BlockVector &parent,
1709  const size_type global_index)
1710  :
1711  parent (&parent),
1712  global_index (global_index)
1713  {
1714  // find which block we are
1715  // in. for this, take into
1716  // account that it happens at
1717  // times that people want to
1718  // initialize iterators
1719  // past-the-end
1720  if (global_index < parent.size())
1721  {
1722  const std::pair<size_type, size_type>
1723  indices = parent.block_indices.global_to_local(global_index);
1724  current_block = indices.first;
1725  index_within_block = indices.second;
1726 
1727  next_break_backward
1728  = parent.block_indices.local_to_global (current_block, 0);
1729  next_break_forward
1730  = (parent.block_indices.local_to_global (current_block, 0)
1731  +parent.block_indices.block_size(current_block)-1);
1732  }
1733  else
1734  // past the end. only have one
1735  // value for this
1736  {
1737  this->global_index = parent.size ();
1738  current_block = parent.n_blocks();
1739  index_within_block = 0;
1740  next_break_backward = global_index;
1741  next_break_forward = numbers::invalid_size_type;
1742  };
1743  }
1744 
1745 
1746 
1747  template <class BlockVectorType, bool constness>
1748  void
1749  Iterator<BlockVectorType,constness>::move_forward ()
1750  {
1751  if (global_index != next_break_forward)
1752  ++index_within_block;
1753  else
1754  {
1755  // ok, we traverse a boundary
1756  // between blocks:
1757  index_within_block = 0;
1758  ++current_block;
1759 
1760  // break backwards is now old
1761  // break forward
1762  next_break_backward = next_break_forward+1;
1763 
1764  // compute new break forward
1765  if (current_block < parent->block_indices.size())
1766  next_break_forward
1767  += parent->block_indices.block_size(current_block);
1768  else
1769  // if we are beyond the end,
1770  // then move the next
1771  // boundary arbitrarily far
1772  // away
1773  next_break_forward = numbers::invalid_size_type;
1774  };
1775 
1776  ++global_index;
1777  }
1778 
1779 
1780 
1781  template <class BlockVectorType, bool constness>
1782  void
1783  Iterator<BlockVectorType,constness>::move_backward ()
1784  {
1785  if (global_index != next_break_backward)
1786  --index_within_block;
1787  else if (current_block != 0)
1788  {
1789  // ok, we traverse a boundary
1790  // between blocks:
1791  --current_block;
1792  index_within_block = parent->block_indices.block_size(current_block)-1;
1793 
1794  // break forwards is now old
1795  // break backward
1796  next_break_forward = next_break_backward-1;
1797 
1798  // compute new break forward
1799  next_break_backward
1800  -= parent->block_indices.block_size (current_block);
1801  }
1802  else
1803  // current block was 0, we now
1804  // get into unspecified terrain
1805  {
1806  --current_block;
1807  index_within_block = numbers::invalid_size_type;
1808  next_break_forward = 0;
1809  next_break_backward = 0;
1810  };
1811 
1812  --global_index;
1813  }
1814 
1815 
1816  } // namespace BlockVectorIterators
1817 
1818 } //namespace internal
1819 
1820 
1821 template <class VectorType>
1822 inline
1824 {}
1825 
1826 
1827 
1828 template <class VectorType>
1829 inline
1830 std::size_t
1832 {
1833  return block_indices.total_size();
1834 }
1835 
1836 
1837 
1838 template <class VectorType>
1839 inline
1840 IndexSet
1842 {
1843  IndexSet is (size());
1844 
1845  // copy index sets from blocks into the global one, shifted
1846  // by the appropriate amount for each block
1847  for (unsigned int b=0; b<n_blocks(); ++b)
1848  {
1849  IndexSet x = block(b).locally_owned_elements();
1850  is.add_indices(x, block_indices.block_start(b));
1851  }
1852 
1853  is.compress();
1854 
1855  return is;
1856 }
1857 
1858 
1859 
1860 template <class VectorType>
1861 inline
1862 unsigned int
1864 {
1865  return block_indices.size();
1866 }
1867 
1868 
1869 template <class VectorType>
1870 inline
1872 BlockVectorBase<VectorType>::block (const unsigned int i)
1873 {
1874  Assert(i<n_blocks(), ExcIndexRange(i,0,n_blocks()));
1875 
1876  return components[i];
1877 }
1878 
1879 
1880 
1881 template <class VectorType>
1882 inline
1884 BlockVectorBase<VectorType>::block (const unsigned int i) const
1885 {
1886  Assert(i<n_blocks(), ExcIndexRange(i,0,n_blocks()));
1887 
1888  return components[i];
1889 }
1890 
1891 
1892 
1893 template <class VectorType>
1894 inline
1895 const BlockIndices &
1897 {
1898  return block_indices;
1899 }
1900 
1901 
1902 template <class VectorType>
1903 inline
1904 void
1906 {
1907  std::vector<size_type> sizes (n_blocks());
1908 
1909  for (size_type i=0; i<n_blocks(); ++i)
1910  sizes[i] = block(i).size();
1911 
1912  block_indices.reinit(sizes);
1913 }
1914 
1915 
1916 
1917 template <class VectorType>
1918 inline
1919 void
1920 BlockVectorBase<VectorType>::compress (::VectorOperation::values operation)
1921 {
1922  for (unsigned int i=0; i<n_blocks(); ++i)
1923  block(i).compress (operation);
1924 }
1925 
1926 
1927 
1928 template <class VectorType>
1929 inline
1930 void
1932 {
1933  compress(VectorOperation::unknown);
1934 }
1935 
1936 
1937 
1938 template <class VectorType>
1939 inline
1942 {
1943  return iterator(*this, 0U);
1944 }
1945 
1946 
1947 
1948 template <class VectorType>
1949 inline
1952 {
1953  return const_iterator(*this, 0U);
1954 }
1955 
1956 
1957 template <class VectorType>
1958 inline
1961 {
1962  return iterator(*this, size());
1963 }
1964 
1965 
1966 
1967 template <class VectorType>
1968 inline
1971 {
1972  return const_iterator(*this, size());
1973 }
1974 
1975 
1976 template <class VectorType>
1977 inline
1978 bool
1980 (const size_type global_index) const
1981 {
1982  const std::pair<size_type,size_type> local_index
1983  = block_indices.global_to_local (global_index);
1984 
1985  return components[local_index.first].in_local_range (global_index);
1986 }
1987 
1988 
1989 template <class VectorType>
1990 bool
1992 {
1993  for (size_type i=0; i<n_blocks(); ++i)
1994  if (components[i].all_zero() == false)
1995  return false;
1996 
1997  return true;
1998 }
1999 
2000 
2001 
2002 template <class VectorType>
2003 bool
2005 {
2006  for (size_type i=0; i<n_blocks(); ++i)
2007  if (components[i].is_non_negative() == false)
2008  return false;
2009 
2010  return true;
2011 }
2012 
2013 
2014 
2015 template <class VectorType>
2016 typename BlockVectorBase<VectorType>::value_type
2019 {
2020  Assert (n_blocks() == v.n_blocks(),
2022 
2023  value_type sum = 0.;
2024  for (size_type i=0; i<n_blocks(); ++i)
2025  sum += components[i]*v.components[i];
2026 
2027  return sum;
2028 }
2029 
2030 
2031 template <class VectorType>
2034 {
2035  real_type sum = 0.;
2036  for (size_type i=0; i<n_blocks(); ++i)
2037  sum += components[i].norm_sqr();
2038 
2039  return sum;
2040 }
2041 
2042 
2043 
2044 template <class VectorType>
2045 typename BlockVectorBase<VectorType>::value_type
2047 {
2048  value_type sum = 0.;
2049  for (size_type i=0; i<n_blocks(); ++i)
2050  sum += components[i].mean_value() * components[i].size();
2051 
2052  return sum/size();
2053 }
2054 
2055 
2056 
2057 template <class VectorType>
2060 {
2061  real_type sum = 0.;
2062  for (size_type i=0; i<n_blocks(); ++i)
2063  sum += components[i].l1_norm();
2064 
2065  return sum;
2066 }
2067 
2068 
2069 
2070 template <class VectorType>
2073 {
2074  return std::sqrt(norm_sqr());
2075 }
2076 
2077 
2078 
2079 template <class VectorType>
2082 {
2083  real_type sum = 0.;
2084  for (size_type i=0; i<n_blocks(); ++i)
2085  {
2086  value_type newval = components[i].linfty_norm();
2087  if (sum<newval)
2088  sum = newval;
2089  }
2090  return sum;
2091 }
2092 
2093 
2094 
2095 template <class VectorType>
2098 {
2099  add (v);
2100  return *this;
2101 }
2102 
2103 
2104 
2105 template <class VectorType>
2108 {
2109  Assert (n_blocks() == v.n_blocks(),
2111 
2112  for (size_type i=0; i<n_blocks(); ++i)
2113  {
2114  components[i] -= v.components[i];
2115  }
2116  return *this;
2117 }
2118 
2119 
2120 
2121 template <class VectorType>
2122 template <typename Number>
2123 inline
2124 void
2125 BlockVectorBase<VectorType>::add (const std::vector<size_type> &indices,
2126  const std::vector<Number> &values)
2127 {
2128  Assert (indices.size() == values.size(),
2129  ExcDimensionMismatch(indices.size(), values.size()));
2130  add (indices.size(), &indices[0], &values[0]);
2131 }
2132 
2133 
2134 
2135 template <class VectorType>
2136 template <typename Number>
2137 inline
2138 void
2139 BlockVectorBase<VectorType>::add (const std::vector<size_type> &indices,
2140  const Vector<Number> &values)
2141 {
2142  Assert (indices.size() == values.size(),
2143  ExcDimensionMismatch(indices.size(), values.size()));
2144  const size_type n_indices = indices.size();
2145  for (size_type i=0; i<n_indices; ++i)
2146  (*this)(indices[i]) += values(i);
2147 }
2148 
2149 
2150 
2151 template <class VectorType>
2152 template <typename Number>
2153 inline
2154 void
2155 BlockVectorBase<VectorType>::add (const size_type n_indices,
2156  const size_type *indices,
2157  const Number *values)
2158 {
2159  for (size_type i=0; i<n_indices; ++i)
2160  (*this)(indices[i]) += values[i];
2161 }
2162 
2163 
2164 
2165 template <class VectorType>
2166 void BlockVectorBase<VectorType>::add (const value_type a)
2167 {
2169 
2170  for (size_type i=0; i<n_blocks(); ++i)
2171  {
2172  components[i].add(a);
2173  }
2174 }
2175 
2176 
2177 
2178 template <class VectorType>
2180 {
2181  Assert (n_blocks() == v.n_blocks(),
2183 
2184  for (size_type i=0; i<n_blocks(); ++i)
2185  {
2186  components[i].add(v.components[i]);
2187  }
2188 }
2189 
2190 
2191 
2192 template <class VectorType>
2193 void BlockVectorBase<VectorType>::add (const value_type a,
2194  const BlockVectorBase<VectorType> &v)
2195 {
2196 
2198 
2199  Assert (n_blocks() == v.n_blocks(),
2201 
2202  for (size_type i=0; i<n_blocks(); ++i)
2203  {
2204  components[i].add(a, v.components[i]);
2205  }
2206 }
2207 
2208 
2209 
2210 template <class VectorType>
2211 void BlockVectorBase<VectorType>::add (const value_type a,
2212  const BlockVectorBase<VectorType> &v,
2213  const value_type b,
2214  const BlockVectorBase<VectorType> &w)
2215 {
2216 
2219 
2220  Assert (n_blocks() == v.n_blocks(),
2222  Assert (n_blocks() == w.n_blocks(),
2224 
2225 
2226  for (size_type i=0; i<n_blocks(); ++i)
2227  {
2228  components[i].add(a, v.components[i], b, w.components[i]);
2229  }
2230 }
2231 
2232 
2233 
2234 template <class VectorType>
2235 void BlockVectorBase<VectorType>::sadd (const value_type x,
2236  const BlockVectorBase<VectorType> &v)
2237 {
2238 
2240 
2241  Assert (n_blocks() == v.n_blocks(),
2243 
2244  for (size_type i=0; i<n_blocks(); ++i)
2245  {
2246  components[i].sadd(x, v.components[i]);
2247  }
2248 }
2249 
2250 
2251 
2252 template <class VectorType>
2253 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
2254  const BlockVectorBase<VectorType> &v)
2255 {
2256 
2259 
2260  Assert (n_blocks() == v.n_blocks(),
2262 
2263  for (size_type i=0; i<n_blocks(); ++i)
2264  {
2265  components[i].sadd(x, a, v.components[i]);
2266  }
2267 }
2268 
2269 
2270 
2271 template <class VectorType>
2272 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
2273  const BlockVectorBase<VectorType> &v,
2274  const value_type b,
2275  const BlockVectorBase<VectorType> &w)
2276 {
2277 
2281 
2282  Assert (n_blocks() == v.n_blocks(),
2284  Assert (n_blocks() == w.n_blocks(),
2286 
2287  for (size_type i=0; i<n_blocks(); ++i)
2288  {
2289  components[i].sadd(x, a, v.components[i], b, w.components[i]);
2290  }
2291 }
2292 
2293 
2294 
2295 template <class VectorType>
2296 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
2297  const BlockVectorBase<VectorType> &v,
2298  const value_type b,
2299  const BlockVectorBase<VectorType> &w,
2300  const value_type c,
2301  const BlockVectorBase<VectorType> &y)
2302 {
2303 
2308 
2309  Assert (n_blocks() == v.n_blocks(),
2311  Assert (n_blocks() == w.n_blocks(),
2313  Assert (n_blocks() == y.n_blocks(),
2315 
2316  for (size_type i=0; i<n_blocks(); ++i)
2317  {
2318  components[i].sadd(x, a, v.components[i],
2319  b, w.components[i], c, y.components[i]);
2320  }
2321 }
2322 
2323 
2324 
2325 template <class VectorType>
2326 template <class BlockVector2>
2327 void BlockVectorBase<VectorType>::scale (const BlockVector2 &v)
2328 {
2329  Assert (n_blocks() == v.n_blocks(),
2330  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
2331  for (size_type i=0; i<n_blocks(); ++i)
2332  components[i].scale(v.block(i));
2333 }
2334 
2335 
2336 
2337 template <class VectorType>
2338 void BlockVectorBase<VectorType>::equ (const value_type a,
2339  const BlockVectorBase<VectorType> &v,
2340  const value_type b,
2341  const BlockVectorBase<VectorType> &w)
2342 {
2343 
2346 
2347  Assert (n_blocks() == v.n_blocks(),
2349  Assert (n_blocks() == w.n_blocks(),
2351 
2352  for (size_type i=0; i<n_blocks(); ++i)
2353  {
2354  components[i].equ( a, v.components[i], b, w.components[i]);
2355  }
2356 }
2357 
2358 
2359 
2360 template <class VectorType>
2361 std::size_t
2363 {
2364  std::size_t mem = sizeof(this->n_blocks());
2365  for (size_type i=0; i<this->components.size(); ++i)
2366  mem += MemoryConsumption::memory_consumption (this->components[i]);
2367  mem += MemoryConsumption::memory_consumption (this->block_indices);
2368  return mem;
2369 }
2370 
2371 
2372 
2373 template <class VectorType>
2374 template <class BlockVector2>
2375 void BlockVectorBase<VectorType>::equ (const value_type a,
2376  const BlockVector2 &v)
2377 {
2378 
2380 
2381  Assert (n_blocks() == v.n_blocks(),
2382  ExcDimensionMismatch(n_blocks(), v.n_blocks()));
2383 
2384  for (size_type i=0; i<n_blocks(); ++i)
2385  components[i].equ( a, v.components[i]);
2386 }
2387 
2388 
2389 
2390 template <class VectorType>
2392 {
2393  for (size_type i=0; i<n_blocks(); ++i)
2394  block(i).update_ghost_values ();
2395 }
2396 
2397 
2398 
2399 template <class VectorType>
2401 BlockVectorBase<VectorType>::operator = (const value_type s)
2402 {
2403 
2405 
2406  for (size_type i=0; i<n_blocks(); ++i)
2407  components[i] = s;
2408 
2409  return *this;
2410 }
2411 
2412 
2413 template <class VectorType>
2416 {
2418 
2419  for (size_type i=0; i<n_blocks(); ++i)
2420  components[i] = v.components[i];
2421 
2422  return *this;
2423 }
2424 
2425 
2426 template <class VectorType>
2427 template <class VectorType2>
2430 {
2432 
2433  for (size_type i=0; i<n_blocks(); ++i)
2434  components[i] = v.components[i];
2435 
2436  return *this;
2437 }
2438 
2439 
2440 
2441 template <class VectorType>
2443 BlockVectorBase<VectorType>::operator = (const VectorType &v)
2444 {
2445  Assert (size() == v.size(),
2446  ExcDimensionMismatch(size(), v.size()));
2447 
2448  size_type index_v = 0;
2449  for (size_type b=0; b<n_blocks(); ++b)
2450  for (size_type i=0; i<block(b).size(); ++i, ++index_v)
2451  block(b)(i) = v(index_v);
2452 
2453  return *this;
2454 }
2455 
2456 
2457 
2458 template <class VectorType>
2459 template <class VectorType2>
2460 inline
2461 bool
2464 {
2465  Assert (block_indices == v.block_indices, ExcDifferentBlockIndices());
2466 
2467  for (size_type i=0; i<n_blocks(); ++i)
2468  if ( ! (components[i] == v.components[i]))
2469  return false;
2470 
2471  return true;
2472 }
2473 
2474 
2475 
2476 template <class VectorType>
2477 inline
2479 BlockVectorBase<VectorType>::operator *= (const value_type factor)
2480 {
2481 
2483 
2484  for (size_type i=0; i<n_blocks(); ++i)
2485  components[i] *= factor;
2486 
2487  return *this;
2488 }
2489 
2490 
2491 
2492 template <class VectorType>
2493 inline
2495 BlockVectorBase<VectorType>::operator /= (const value_type factor)
2496 {
2497 
2499  Assert (factor > 0., ExcDivideByZero() );
2500 
2501  for (size_type i=0; i<n_blocks(); ++i)
2502  components[i] /= factor;
2503 
2504  return *this;
2505 }
2506 
2507 
2508 template <class VectorType>
2509 inline
2510 typename BlockVectorBase<VectorType>::value_type
2511 BlockVectorBase<VectorType>::operator() (const size_type i) const
2512 {
2513  const std::pair<size_type,size_type> local_index
2514  = block_indices.global_to_local (i);
2515  return components[local_index.first](local_index.second);
2516 }
2517 
2518 
2519 
2520 template <class VectorType>
2521 inline
2522 typename BlockVectorBase<VectorType>::reference
2524 {
2525  const std::pair<size_type,size_type> local_index
2526  = block_indices.global_to_local (i);
2527  return components[local_index.first](local_index.second);
2528 }
2529 
2530 
2531 
2532 template <class VectorType>
2533 inline
2534 typename BlockVectorBase<VectorType>::value_type
2535 BlockVectorBase<VectorType>::operator[] (const size_type i) const
2536 {
2537  return operator()(i);
2538 }
2539 
2540 
2541 
2542 template <class VectorType>
2543 inline
2544 typename BlockVectorBase<VectorType>::reference
2546 {
2547  return operator()(i);
2548 }
2549 
2550 
2551 
2552 template <typename VectorType>
2553 template <typename OtherNumber>
2554 inline
2555 void BlockVectorBase<VectorType>::extract_subvector_to (const std::vector<size_type> &indices,
2556  std::vector<OtherNumber> &values) const
2557 {
2558  for (size_type i = 0; i < indices.size(); ++i)
2559  values[i] = operator()(indices[i]);
2560 }
2561 
2562 
2563 
2564 template <typename VectorType>
2565 template <typename ForwardIterator, typename OutputIterator>
2566 inline
2567 void BlockVectorBase<VectorType>::extract_subvector_to (ForwardIterator indices_begin,
2568  const ForwardIterator indices_end,
2569  OutputIterator values_begin) const
2570 {
2571  while (indices_begin != indices_end)
2572  {
2573  *values_begin = operator()(*indices_begin);
2574  indices_begin++;
2575  values_begin++;
2576  }
2577 }
2578 
2579 #endif // DOXYGEN
2580 
2581 DEAL_II_NAMESPACE_CLOSE
2582 
2583 #endif
const types::global_dof_index invalid_size_type
Definition: types.h:211
static yes_type check_for_block_vector(const BlockVectorBase< T > *)
std::size_t size() const
BaseClass::value_type value_type
Definition: block_vector.h:80
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
bool operator==(const BlockVectorBase< VectorType2 > &v) const
void add(const std::vector< size_type > &indices, const std::vector< Number > &values)
real_type linfty_norm() const
bool operator>(const Iterator &i) const
Auxiliary class aiding in the handling of block structures like in BlockVector or FESystem...
Definition: block_indices.h:55
real_type norm_sqr() const
void collect_sizes()
STL namespace.
bool operator>=(const Iterator &i) const
unsigned int n_blocks() const
value_type operator[](const size_type i) const
BlockVectorBase & operator/=(const value_type factor)
bool operator==(const Iterator &i) const
bool is_finite(const double x)
void scale(const BlockVector2 &v)
dereference_type operator*() const
Iterator & operator+=(const difference_type &d)
bool is_non_negative() const
value_type operator()(const size_type i) const
BlockIndices block_indices
size_type block_start(const unsigned int i) const
bool in_local_range(const size_type global_index) const
Iterator & operator=(const Iterator &c)
dereference_type operator[](const difference_type d) const
difference_type operator-(const Iterator &i) const
T sum(const T &t, const MPI_Comm &mpi_communicator)
Definition: mpi.h:447
unsigned int global_dof_index
Definition: types.h:100
bool operator<(const Iterator &i) const
#define Assert(cond, exc)
Definition: exceptions.h:299
void update_ghost_values() const
BlockType::real_type real_type
std::size_t memory_consumption(const T &t)
const BlockIndices & get_block_indices() const
void reinit(const unsigned int n_blocks, const size_type n_elements_per_block)
size_type total_size() const
Types< BlockVectorType, constness >::value_type value_type
bool operator!=(const Iterator &i) const
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
value_type mean_value() const
bool all_zero() const
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
DeclException0(ExcPointerToDifferentVectors)
static const bool value
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
BlockVectorBase & operator+=(const BlockVectorBase &V)
std::vector< VectorType > components
size_type block_size(const unsigned int i) const
BlockVectorBase & operator*=(const value_type factor)
void sadd(const value_type s, const BlockVectorBase &V)
::ExceptionBase & ExcNumberNotFinite()
BlockVectorBase & operator-=(const BlockVectorBase &V)
Iterator< BlockVectorType,!constness > InverseConstnessIterator
Iterator operator+(const difference_type &d) const
bool operator<=(const Iterator &i) const
iterator end()
Iterator & operator-=(const difference_type &d)
value_type operator*(const BlockVectorBase &V) const
real_type l2_norm() const
void compress() DEAL_II_DEPRECATED
IndexSet locally_owned_elements() const
unsigned int size() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::random_access_iterator_tag iterator_type
::ExceptionBase & ExcDivideByZero()
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
real_type l1_norm() const
void equ(const value_type a, const BlockVector2 &V)
static const bool supports_distributed_data
Types< BlockVectorType, constness >::BlockVector BlockVector
std::size_t memory_consumption() const