Reference documentation for deal.II version Git 7916fb8 2017-12-10 11:04:15 +0100
trilinos_vector.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 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_trilinos_vector_h
17 #define dealii_trilinos_vector_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #ifdef DEAL_II_WITH_TRILINOS
23 # include <deal.II/base/subscriptor.h>
24 # include <deal.II/base/index_set.h>
25 # include <deal.II/base/utilities.h>
26 # include <deal.II/base/mpi.h>
27 # include <deal.II/lac/exceptions.h>
28 # include <deal.II/lac/vector.h>
29 # include <deal.II/lac/vector_operation.h>
30 # include <deal.II/lac/vector_type_traits.h>
31 
32 # include <vector>
33 # include <utility>
34 # include <memory>
35 
36 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
37 # include <Epetra_ConfigDefs.h>
38 # ifdef DEAL_II_WITH_MPI // only if MPI is installed
39 # include <mpi.h>
40 # include <Epetra_MpiComm.h>
41 # else
42 # include <Epetra_SerialComm.h>
43 # endif
44 # include <Epetra_FEVector.h>
45 # include <Epetra_Map.h>
46 # include <Epetra_LocalMap.h>
47 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
48 
49 DEAL_II_NAMESPACE_OPEN
50 
51 
62 namespace TrilinosWrappers
63 {
64  class SparseMatrix;
65 
76  namespace internal
77  {
81  typedef ::types::global_dof_index size_type;
82 
92  class VectorReference
93  {
94  private:
99  VectorReference (MPI::Vector &vector,
100  const size_type index);
101 
102  public:
103 
115  const VectorReference &
116  operator= (const VectorReference &r) const;
117 
121  VectorReference &
122  operator= (const VectorReference &r);
123 
127  const VectorReference &
128  operator= (const TrilinosScalar &s) const;
129 
133  const VectorReference &
134  operator+= (const TrilinosScalar &s) const;
135 
139  const VectorReference &
140  operator-= (const TrilinosScalar &s) const;
141 
145  const VectorReference &
146  operator*= (const TrilinosScalar &s) const;
147 
151  const VectorReference &
152  operator/= (const TrilinosScalar &s) const;
153 
158  operator TrilinosScalar () const;
159 
164  int,
165  << "An error with error number " << arg1
166  << " occurred while calling a Trilinos function");
167 
168  private:
172  MPI::Vector &vector;
173 
177  const size_type index;
178 
183  friend class ::TrilinosWrappers::MPI::Vector;
184  };
185  }
190  namespace
191  {
192 #ifndef DEAL_II_WITH_64BIT_INDICES
193  // define a helper function that queries the global ID of local ID of
194  // an Epetra_BlockMap object by calling either the 32- or 64-bit
195  // function necessary.
196  inline
197  int gid(const Epetra_BlockMap &map, int i)
198  {
199  return map.GID(i);
200  }
201 #else
202  // define a helper function that queries the global ID of local ID of
203  // an Epetra_BlockMap object by calling either the 32- or 64-bit
204  // function necessary.
205  inline
206  long long int gid(const Epetra_BlockMap &map, int i)
207  {
208  return map.GID64(i);
209  }
210 #endif
211  }
212 
219  namespace MPI
220  {
221  class BlockVector;
222 
399  class Vector : public Subscriptor
400  {
401  public:
407  typedef TrilinosScalar value_type;
408  typedef TrilinosScalar real_type;
409  typedef ::types::global_dof_index size_type;
410  typedef value_type *iterator;
411  typedef const value_type *const_iterator;
412  typedef internal::VectorReference reference;
413  typedef const internal::VectorReference const_reference;
414 
424  Vector ();
425 
429  Vector (const Vector &v);
430 
449  explicit Vector (const IndexSet &parallel_partitioning,
450  const MPI_Comm &communicator = MPI_COMM_WORLD);
451 
463  Vector (const IndexSet &local,
464  const IndexSet &ghost,
465  const MPI_Comm &communicator = MPI_COMM_WORLD);
466 
481  Vector (const IndexSet &parallel_partitioning,
482  const Vector &v,
483  const MPI_Comm &communicator = MPI_COMM_WORLD);
484 
497  template <typename Number>
498  Vector (const IndexSet &parallel_partitioning,
499  const ::Vector<Number> &v,
500  const MPI_Comm &communicator = MPI_COMM_WORLD);
501 
506  Vector (Vector &&v);
507 
511  ~Vector () = default;
512 
517  void clear ();
518 
542  void reinit (const Vector &v,
543  const bool omit_zeroing_entries = false,
544  const bool allow_different_maps = false);
545 
568  void reinit (const IndexSet &parallel_partitioning,
569  const MPI_Comm &communicator = MPI_COMM_WORLD,
570  const bool omit_zeroing_entries = false);
571 
597  void reinit (const IndexSet &locally_owned_entries,
598  const IndexSet &ghost_entries,
599  const MPI_Comm &communicator = MPI_COMM_WORLD,
600  const bool vector_writable = false);
601 
605  void reinit (const BlockVector &v,
606  const bool import_data = false);
607 
624  void compress (::VectorOperation::values operation);
625 
638  Vector &operator= (const TrilinosScalar s);
639 
645  Vector &operator= (const Vector &v);
646 
651  Vector &operator= (Vector &&v);
652 
660  template <typename Number>
661  Vector &operator= (const ::Vector<Number> &v);
662 
681  (const ::TrilinosWrappers::SparseMatrix &matrix,
682  const Vector &vector);
683 
689  bool operator== (const Vector &v) const;
690 
696  bool operator!= (const Vector &v) const;
697 
701  size_type size () const;
702 
714  size_type local_size () const;
715 
737  std::pair<size_type, size_type> local_range () const;
738 
746  bool in_local_range (const size_type index) const;
747 
762 
770  bool has_ghost_elements() const;
771 
778  void update_ghost_values () const;
779 
784  TrilinosScalar operator* (const Vector &vec) const;
785 
789  real_type norm_sqr () const;
790 
794  TrilinosScalar mean_value () const;
795 
799  TrilinosScalar min () const;
800 
804  TrilinosScalar max () const;
805 
809  real_type l1_norm () const;
810 
815  real_type l2_norm () const;
816 
821  real_type lp_norm (const TrilinosScalar p) const;
822 
826  real_type linfty_norm () const;
827 
843  TrilinosScalar add_and_dot (const TrilinosScalar a,
844  const Vector &V,
845  const Vector &W);
846 
852  bool all_zero () const;
853 
859  bool is_non_negative () const;
861 
862 
867 
877  reference
878  operator() (const size_type index);
879 
889  TrilinosScalar
890  operator() (const size_type index) const;
891 
897  reference
898  operator[] (const size_type index);
899 
905  TrilinosScalar
906  operator[] (const size_type index) const;
907 
923  void extract_subvector_to (const std::vector<size_type> &indices,
924  std::vector<TrilinosScalar> &values) const;
925 
953  template <typename ForwardIterator, typename OutputIterator>
954  void extract_subvector_to (ForwardIterator indices_begin,
955  const ForwardIterator indices_end,
956  OutputIterator values_begin) const;
957 
968  iterator begin ();
969 
974  const_iterator begin () const;
975 
980  iterator end ();
981 
986  const_iterator end () const;
987 
989 
990 
995 
1002  void set (const std::vector<size_type> &indices,
1003  const std::vector<TrilinosScalar> &values);
1004 
1009  void set (const std::vector<size_type> &indices,
1010  const ::Vector<TrilinosScalar> &values);
1011 
1017  void set (const size_type n_elements,
1018  const size_type *indices,
1019  const TrilinosScalar *values);
1020 
1025  void add (const std::vector<size_type> &indices,
1026  const std::vector<TrilinosScalar> &values);
1027 
1032  void add (const std::vector<size_type> &indices,
1033  const ::Vector<TrilinosScalar> &values);
1034 
1040  void add (const size_type n_elements,
1041  const size_type *indices,
1042  const TrilinosScalar *values);
1043 
1047  Vector &operator*= (const TrilinosScalar factor);
1048 
1052  Vector &operator/= (const TrilinosScalar factor);
1053 
1057  Vector &operator+= (const Vector &V);
1058 
1062  Vector &operator-= (const Vector &V);
1063 
1068  void add (const TrilinosScalar s);
1069 
1082  void add (const Vector &V,
1083  const bool allow_different_maps = false);
1084 
1088  void add (const TrilinosScalar a,
1089  const Vector &V);
1090 
1094  void add (const TrilinosScalar a,
1095  const Vector &V,
1096  const TrilinosScalar b,
1097  const Vector &W);
1098 
1103  void sadd (const TrilinosScalar s,
1104  const Vector &V);
1105 
1109  void sadd (const TrilinosScalar s,
1110  const TrilinosScalar a,
1111  const Vector &V);
1112 
1118  void scale (const Vector &scaling_factors);
1119 
1123  void equ (const TrilinosScalar a,
1124  const Vector &V);
1126 
1131 
1136  const Epetra_MultiVector &trilinos_vector () const;
1137 
1142  Epetra_FEVector &trilinos_vector ();
1143 
1148  const Epetra_Map &vector_partitioner () const;
1149 
1157  void print (std::ostream &out,
1158  const unsigned int precision = 3,
1159  const bool scientific = true,
1160  const bool across = true) const;
1161 
1175  void swap (Vector &v);
1176 
1180  std::size_t memory_consumption () const;
1181 
1186  const MPI_Comm &get_mpi_communicator () const;
1188 
1193 
1198  int,
1199  << "An error with error number " << arg1
1200  << " occurred while calling a Trilinos function");
1201 
1206  size_type, size_type, size_type, size_type,
1207  << "You tried to access element " << arg1
1208  << " of a distributed vector, but this element is not stored "
1209  << "on the current processor. Note: There are "
1210  << arg2 << " elements stored "
1211  << "on the current processor from within the range "
1212  << arg3 << " through " << arg4
1213  << " but Trilinos vectors need not store contiguous "
1214  << "ranges on each processor, and not every element in "
1215  << "this range may in fact be stored locally.");
1216 
1217  private:
1229  Epetra_CombineMode last_action;
1230 
1236 
1242 
1248  std::shared_ptr<Epetra_FEVector> vector;
1249 
1255  std::shared_ptr<Epetra_MultiVector> nonlocal_vector;
1256 
1261 
1266  };
1267 
1268 
1269 
1270 
1271 // ------------------- inline and template functions --------------
1272 
1273 
1282  inline
1283  void swap (Vector &u, Vector &v)
1284  {
1285  u.swap (v);
1286  }
1287  }
1288 
1289 #ifndef DOXYGEN
1290 
1291  namespace internal
1292  {
1293  inline
1294  VectorReference::VectorReference (MPI::Vector &vector,
1295  const size_type index)
1296  :
1297  vector (vector),
1298  index (index)
1299  {}
1300 
1301 
1302  inline
1303  const VectorReference &
1304  VectorReference::operator= (const VectorReference &r) const
1305  {
1306  // as explained in the class
1307  // documentation, this is not the copy
1308  // operator. so simply pass on to the
1309  // "correct" assignment operator
1310  *this = static_cast<TrilinosScalar> (r);
1311 
1312  return *this;
1313  }
1314 
1315 
1316 
1317  inline
1318  VectorReference &
1319  VectorReference::operator= (const VectorReference &r)
1320  {
1321  // as above
1322  *this = static_cast<TrilinosScalar> (r);
1323 
1324  return *this;
1325  }
1326 
1327 
1328  inline
1329  const VectorReference &
1330  VectorReference::operator= (const TrilinosScalar &value) const
1331  {
1332  vector.set (1, &index, &value);
1333  return *this;
1334  }
1335 
1336 
1337 
1338  inline
1339  const VectorReference &
1340  VectorReference::operator+= (const TrilinosScalar &value) const
1341  {
1342  vector.add (1, &index, &value);
1343  return *this;
1344  }
1345 
1346 
1347 
1348  inline
1349  const VectorReference &
1350  VectorReference::operator-= (const TrilinosScalar &value) const
1351  {
1352  TrilinosScalar new_value = -value;
1353  vector.add (1, &index, &new_value);
1354  return *this;
1355  }
1356 
1357 
1358 
1359  inline
1360  const VectorReference &
1361  VectorReference::operator*= (const TrilinosScalar &value) const
1362  {
1363  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1364  vector.set (1, &index, &new_value);
1365  return *this;
1366  }
1367 
1368 
1369 
1370  inline
1371  const VectorReference &
1372  VectorReference::operator/= (const TrilinosScalar &value) const
1373  {
1374  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1375  vector.set (1, &index, &new_value);
1376  return *this;
1377  }
1378  }
1379 
1380  namespace MPI
1381  {
1382  inline
1383  bool
1384  Vector::in_local_range (const size_type index) const
1385  {
1386  std::pair<size_type, size_type> range = local_range();
1387 
1388  return ((index >= range.first) && (index < range.second));
1389  }
1390 
1391 
1392 
1393  inline
1394  IndexSet
1396  {
1398  ExcMessage("The locally owned elements have not been properly initialized!"
1399  " This happens for example if this object has been initialized"
1400  " with exactly one overlapping IndexSet."));
1401  return owned_elements;
1402  }
1403 
1404 
1405 
1406  inline
1407  bool
1409  {
1410  return has_ghosts;
1411  }
1412 
1413 
1414 
1415  inline
1416  void
1418  {}
1419 
1420 
1421 
1422  inline
1423  internal::VectorReference
1424  Vector::operator() (const size_type index)
1425  {
1426  return internal::VectorReference (*this, index);
1427  }
1428 
1429 
1430 
1431  inline
1432  internal::VectorReference
1433  Vector::operator[] (const size_type index)
1434  {
1435  return operator() (index);
1436  }
1437 
1438 
1439 
1440  inline
1441  TrilinosScalar
1442  Vector::operator[] (const size_type index) const
1443  {
1444  return operator() (index);
1445  }
1446 
1447 
1448 
1449  inline
1450  void Vector::extract_subvector_to (const std::vector<size_type> &indices,
1451  std::vector<TrilinosScalar> &values) const
1452  {
1453  for (size_type i = 0; i < indices.size(); ++i)
1454  values[i] = operator()(indices[i]);
1455  }
1456 
1457 
1458 
1459  template <typename ForwardIterator, typename OutputIterator>
1460  inline
1461  void Vector::extract_subvector_to (ForwardIterator indices_begin,
1462  const ForwardIterator indices_end,
1463  OutputIterator values_begin) const
1464  {
1465  while (indices_begin != indices_end)
1466  {
1467  *values_begin = operator()(*indices_begin);
1468  indices_begin++;
1469  values_begin++;
1470  }
1471  }
1472 
1473 
1474 
1475  inline
1476  Vector::iterator
1477  Vector::begin()
1478  {
1479  return (*vector)[0];
1480  }
1481 
1482 
1483 
1484  inline
1485  Vector::iterator
1486  Vector::end()
1487  {
1488  return (*vector)[0]+local_size();
1489  }
1490 
1491 
1492 
1493  inline
1494  Vector::const_iterator
1495  Vector::begin() const
1496  {
1497  return (*vector)[0];
1498  }
1499 
1500 
1501 
1502  inline
1503  Vector::const_iterator
1504  Vector::end() const
1505  {
1506  return (*vector)[0]+local_size();
1507  }
1508 
1509 
1510 
1511  inline
1512  void
1513  Vector::set (const std::vector<size_type> &indices,
1514  const std::vector<TrilinosScalar> &values)
1515  {
1516  // if we have ghost values, do not allow
1517  // writing to this vector at all.
1519 
1520  Assert (indices.size() == values.size(),
1521  ExcDimensionMismatch(indices.size(),values.size()));
1522 
1523  set (indices.size(), indices.data(), values.data());
1524  }
1525 
1526 
1527 
1528  inline
1529  void
1530  Vector::set (const std::vector<size_type> &indices,
1531  const ::Vector<TrilinosScalar> &values)
1532  {
1533  // if we have ghost values, do not allow
1534  // writing to this vector at all.
1536 
1537  Assert (indices.size() == values.size(),
1538  ExcDimensionMismatch(indices.size(),values.size()));
1539 
1540  set (indices.size(), indices.data(), values.begin());
1541  }
1542 
1543 
1544 
1545  inline
1546  void
1547  Vector::set (const size_type n_elements,
1548  const size_type *indices,
1549  const TrilinosScalar *values)
1550  {
1551  // if we have ghost values, do not allow
1552  // writing to this vector at all.
1554 
1555  if (last_action == Add)
1556  {
1557  const int ierr = vector->GlobalAssemble(Add);
1558  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1559  }
1560 
1561  if (last_action != Insert)
1562  last_action = Insert;
1563 
1564  for (size_type i=0; i<n_elements; ++i)
1565  {
1566  const size_type row = indices[i];
1567  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1568  if (local_row != -1)
1569  (*vector)[0][local_row] = values[i];
1570  else
1571  {
1572  const int ierr = vector->ReplaceGlobalValues (1,
1573  (const TrilinosWrappers::types::int_type *)(&row),
1574  &values[i]);
1575  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1576  compressed = false;
1577  }
1578  // in set operation, do not use the pre-allocated vector for nonlocal
1579  // entries even if it exists. This is to ensure that we really only
1580  // set the elements touched by the set() method and not all contained
1581  // in the nonlocal entries vector (there is no way to distinguish them
1582  // on the receiving processor)
1583  }
1584  }
1585 
1586 
1587 
1588  inline
1589  void
1590  Vector::add (const std::vector<size_type> &indices,
1591  const std::vector<TrilinosScalar> &values)
1592  {
1593  // if we have ghost values, do not allow
1594  // writing to this vector at all.
1596  Assert (indices.size() == values.size(),
1597  ExcDimensionMismatch(indices.size(),values.size()));
1598 
1599  add (indices.size(), indices.data(), values.data());
1600  }
1601 
1602 
1603 
1604  inline
1605  void
1606  Vector::add (const std::vector<size_type> &indices,
1607  const ::Vector<TrilinosScalar> &values)
1608  {
1609  // if we have ghost values, do not allow
1610  // writing to this vector at all.
1612  Assert (indices.size() == values.size(),
1613  ExcDimensionMismatch(indices.size(),values.size()));
1614 
1615  add (indices.size(), indices.data(), values.begin());
1616  }
1617 
1618 
1619 
1620  inline
1621  void
1622  Vector::add (const size_type n_elements,
1623  const size_type *indices,
1624  const TrilinosScalar *values)
1625  {
1626  // if we have ghost values, do not allow
1627  // writing to this vector at all.
1629 
1630  if (last_action != Add)
1631  {
1632  if (last_action == Insert)
1633  {
1634  const int ierr = vector->GlobalAssemble(Insert);
1635  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1636  }
1637  last_action = Add;
1638  }
1639 
1640  for (size_type i=0; i<n_elements; ++i)
1641  {
1642  const size_type row = indices[i];
1643  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1644  if (local_row != -1)
1645  (*vector)[0][local_row] += values[i];
1646  else if (nonlocal_vector.get() == nullptr)
1647  {
1648  const int ierr = vector->SumIntoGlobalValues (1,
1649  (const TrilinosWrappers::types::int_type *)(&row),
1650  &values[i]);
1651  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1652  compressed = false;
1653  }
1654  else
1655  {
1656  // use pre-allocated vector for non-local entries if it exists for
1657  // addition operation
1658  const TrilinosWrappers::types::int_type my_row = nonlocal_vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1659  Assert(my_row != -1,
1660  ExcMessage("Attempted to write into off-processor vector entry "
1661  "that has not be specified as being writable upon "
1662  "initialization"));
1663  (*nonlocal_vector)[0][my_row] += values[i];
1664  compressed = false;
1665  }
1666  }
1667  }
1668 
1669 
1670 
1671  inline
1672  Vector::size_type
1673  Vector::size () const
1674  {
1675 #ifndef DEAL_II_WITH_64BIT_INDICES
1676  return (size_type) (vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID());
1677 #else
1678  return (size_type) (vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64());
1679 #endif
1680  }
1681 
1682 
1683 
1684  inline
1685  Vector::size_type
1686  Vector::local_size () const
1687  {
1688  return (size_type) vector->Map().NumMyElements();
1689  }
1690 
1691 
1692 
1693  inline
1694  std::pair<Vector::size_type, Vector::size_type>
1695  Vector::local_range () const
1696  {
1697 #ifndef DEAL_II_WITH_64BIT_INDICES
1698  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1699  const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID()+1;
1700 #else
1701  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID64();
1702  const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID64()+1;
1703 #endif
1704 
1705  Assert (end-begin == vector->Map().NumMyElements(),
1706  ExcMessage ("This function only makes sense if the elements that this "
1707  "vector stores on the current processor form a contiguous range. "
1708  "This does not appear to be the case for the current vector."));
1709 
1710  return std::make_pair (begin, end);
1711  }
1712 
1713 
1714 
1715  inline
1716  TrilinosScalar
1717  Vector::operator* (const Vector &vec) const
1718  {
1719  Assert (vector->Map().SameAs(vec.vector->Map()),
1722 
1723  TrilinosScalar result;
1724 
1725  const int ierr = vector->Dot(*(vec.vector), &result);
1726  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1727 
1728  return result;
1729  }
1730 
1731 
1732 
1733  inline
1735  Vector::norm_sqr () const
1736  {
1737  const TrilinosScalar d = l2_norm();
1738  return d*d;
1739  }
1740 
1741 
1742 
1743  inline
1744  TrilinosScalar
1745  Vector::mean_value () const
1746  {
1748 
1749  TrilinosScalar mean;
1750  const int ierr = vector->MeanValue (&mean);
1751  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1752 
1753  return mean;
1754  }
1755 
1756 
1757 
1758  inline
1759  TrilinosScalar
1760  Vector::min () const
1761  {
1762  TrilinosScalar min_value;
1763  const int ierr = vector->MinValue (&min_value);
1764  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1765 
1766  return min_value;
1767  }
1768 
1769 
1770 
1771  inline
1772  TrilinosScalar
1773  Vector::max () const
1774  {
1775  TrilinosScalar max_value;
1776  const int ierr = vector->MaxValue (&max_value);
1777  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1778 
1779  return max_value;
1780  }
1781 
1782 
1783 
1784  inline
1786  Vector::l1_norm () const
1787  {
1789 
1790  TrilinosScalar d;
1791  const int ierr = vector->Norm1 (&d);
1792  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1793 
1794  return d;
1795  }
1796 
1797 
1798 
1799  inline
1801  Vector::l2_norm () const
1802  {
1804 
1805  TrilinosScalar d;
1806  const int ierr = vector->Norm2 (&d);
1807  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1808 
1809  return d;
1810  }
1811 
1812 
1813 
1814  inline
1816  Vector::lp_norm (const TrilinosScalar p) const
1817  {
1819 
1820  TrilinosScalar norm = 0;
1821  TrilinosScalar sum=0;
1822  const size_type n_local = local_size();
1823 
1824  // loop over all the elements because
1825  // Trilinos does not support lp norms
1826  for (size_type i=0; i<n_local; i++)
1827  sum += std::pow(std::fabs((*vector)[0][i]), p);
1828 
1829  norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
1830 
1831  return norm;
1832  }
1833 
1834 
1835 
1836  inline
1838  Vector::linfty_norm () const
1839  {
1840  // while we disallow the other
1841  // norm operations on ghosted
1842  // vectors, this particular norm
1843  // is safe to run even in the
1844  // presence of ghost elements
1845  TrilinosScalar d;
1846  const int ierr = vector->NormInf (&d);
1847  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1848 
1849  return d;
1850  }
1851 
1852 
1853 
1854  inline
1855  TrilinosScalar
1856  Vector::add_and_dot (const TrilinosScalar a,
1857  const Vector &V,
1858  const Vector &W)
1859  {
1860  this->add(a, V);
1861  return *this * W;
1862  }
1863 
1864 
1865 
1866  // inline also scalar products, vector
1867  // additions etc. since they are all
1868  // representable by a single Trilinos
1869  // call. This reduces the overhead of the
1870  // wrapper class.
1871  inline
1872  Vector &
1873  Vector::operator*= (const TrilinosScalar a)
1874  {
1875  AssertIsFinite(a);
1876 
1877  const int ierr = vector->Scale(a);
1878  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1879 
1880  return *this;
1881  }
1882 
1883 
1884 
1885  inline
1886  Vector &
1887  Vector::operator/= (const TrilinosScalar a)
1888  {
1889  AssertIsFinite(a);
1890 
1891  const TrilinosScalar factor = 1./a;
1892 
1893  AssertIsFinite(factor);
1894 
1895  const int ierr = vector->Scale(factor);
1896  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1897 
1898  return *this;
1899  }
1900 
1901 
1902 
1903  inline
1904  Vector &
1905  Vector::operator+= (const Vector &v)
1906  {
1907  Assert (size() == v.size(),
1908  ExcDimensionMismatch(size(), v.size()));
1909  Assert (vector->Map().SameAs(v.vector->Map()),
1911 
1912  const int ierr = vector->Update (1.0, *(v.vector), 1.0);
1913  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1914 
1915  return *this;
1916  }
1917 
1918 
1919 
1920  inline
1921  Vector &
1922  Vector::operator-= (const Vector &v)
1923  {
1924  Assert (size() == v.size(),
1925  ExcDimensionMismatch(size(), v.size()));
1926  Assert (vector->Map().SameAs(v.vector->Map()),
1928 
1929  const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
1930  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1931 
1932  return *this;
1933  }
1934 
1935 
1936 
1937  inline
1938  void
1939  Vector::add (const TrilinosScalar s)
1940  {
1941  // if we have ghost values, do not allow
1942  // writing to this vector at all.
1944  AssertIsFinite(s);
1945 
1946  size_type n_local = local_size();
1947  for (size_type i=0; i<n_local; i++)
1948  (*vector)[0][i] += s;
1949  }
1950 
1951 
1952 
1953  inline
1954  void
1955  Vector::add (const TrilinosScalar a,
1956  const Vector &v)
1957  {
1958  // if we have ghost values, do not allow
1959  // writing to this vector at all.
1961  Assert (local_size() == v.local_size(),
1962  ExcDimensionMismatch(local_size(), v.local_size()));
1963 
1964  AssertIsFinite(a);
1965 
1966  const int ierr = vector->Update(a, *(v.vector), 1.);
1967  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1968  }
1969 
1970 
1971 
1972  inline
1973  void
1974  Vector::add (const TrilinosScalar a,
1975  const Vector &v,
1976  const TrilinosScalar b,
1977  const Vector &w)
1978  {
1979  // if we have ghost values, do not allow
1980  // writing to this vector at all.
1982  Assert (local_size() == v.local_size(),
1983  ExcDimensionMismatch(local_size(), v.local_size()));
1984  Assert (local_size() == w.local_size(),
1985  ExcDimensionMismatch(local_size(), w.local_size()));
1986 
1987  AssertIsFinite(a);
1988  AssertIsFinite(b);
1989 
1990  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
1991 
1992  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1993  }
1994 
1995 
1996 
1997  inline
1998  void
1999  Vector::sadd (const TrilinosScalar s,
2000  const Vector &v)
2001  {
2002  // if we have ghost values, do not allow
2003  // writing to this vector at all.
2005  Assert (size() == v.size(),
2006  ExcDimensionMismatch (size(), v.size()));
2007 
2008  AssertIsFinite(s);
2009 
2010  // We assume that the vectors have the same Map
2011  // if the local size is the same and if the vectors are not ghosted
2012  if (local_size() == v.local_size() && !v.has_ghost_elements())
2013  {
2014  Assert (this->vector->Map().SameAs(v.vector->Map())==true,
2016  const int ierr = vector->Update(1., *(v.vector), s);
2017  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2018  }
2019  else
2020  {
2021  (*this) *= s;
2022  this->add(v, true);
2023  }
2024  }
2025 
2026 
2027 
2028  inline
2029  void
2030  Vector::sadd (const TrilinosScalar s,
2031  const TrilinosScalar a,
2032  const Vector &v)
2033  {
2034  // if we have ghost values, do not allow
2035  // writing to this vector at all.
2037  Assert (size() == v.size(),
2038  ExcDimensionMismatch (size(), v.size()));
2039  AssertIsFinite(s);
2040  AssertIsFinite(a);
2041 
2042  // We assume that the vectors have the same Map
2043  // if the local size is the same and if the vectors are not ghosted
2044  if (local_size() == v.local_size() && !v.has_ghost_elements())
2045  {
2046  Assert (this->vector->Map().SameAs(v.vector->Map())==true,
2048  const int ierr = vector->Update(a, *(v.vector), s);
2049  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2050  }
2051  else
2052  {
2053  (*this) *= s;
2054  Vector tmp = v;
2055  tmp *= a;
2056  this->add(tmp, true);
2057  }
2058  }
2059 
2060 
2061 
2062  inline
2063  void
2064  Vector::scale (const Vector &factors)
2065  {
2066  // if we have ghost values, do not allow
2067  // writing to this vector at all.
2069  Assert (local_size() == factors.local_size(),
2070  ExcDimensionMismatch(local_size(), factors.local_size()));
2071 
2072  const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
2073  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2074  }
2075 
2076 
2077 
2078  inline
2079  void
2080  Vector::equ (const TrilinosScalar a,
2081  const Vector &v)
2082  {
2083  // if we have ghost values, do not allow
2084  // writing to this vector at all.
2086  AssertIsFinite(a);
2087 
2088  // If we don't have the same map, copy.
2089  if (vector->Map().SameAs(v.vector->Map())==false)
2090  {
2091  this->sadd(0., a, v);
2092  }
2093  else
2094  {
2095  // Otherwise, just update
2096  int ierr = vector->Update(a, *v.vector, 0.0);
2097  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2098 
2099  last_action = Zero;
2100  }
2101 
2102  }
2103 
2104 
2105 
2106  inline
2107  const Epetra_MultiVector &
2108  Vector::trilinos_vector () const
2109  {
2110  return static_cast<const Epetra_MultiVector &>(*vector);
2111  }
2112 
2113 
2114 
2115  inline
2116  Epetra_FEVector &
2118  {
2119  return *vector;
2120  }
2121 
2122 
2123 
2124  inline
2125  const Epetra_Map &
2127  {
2128  return static_cast<const Epetra_Map &>(vector->Map());
2129  }
2130 
2131 
2132 
2133  inline
2134  const MPI_Comm &
2136  {
2137  static MPI_Comm comm;
2138 
2139 #ifdef DEAL_II_WITH_MPI
2140 
2141  const Epetra_MpiComm *mpi_comm
2142  = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2143  comm = mpi_comm->Comm();
2144 
2145 #else
2146 
2147  comm = MPI_COMM_SELF;
2148 
2149 #endif
2150 
2151  return comm;
2152  }
2153 
2154  template <typename number>
2155  Vector::Vector (const IndexSet &parallel_partitioner,
2156  const ::Vector<number> &v,
2157  const MPI_Comm &communicator)
2158  {
2159  *this = Vector(parallel_partitioner.make_trilinos_map (communicator, true),
2160  v);
2161  owned_elements = parallel_partitioner;
2162  }
2163 
2164 
2165 
2166  inline
2167  Vector &
2168  Vector::operator= (const TrilinosScalar s)
2169  {
2170  AssertIsFinite(s);
2171 
2172  int ierr = vector->PutScalar(s);
2173  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2174 
2175  if (nonlocal_vector.get() != nullptr)
2176  {
2177  ierr = nonlocal_vector->PutScalar(0.);
2178  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2179  }
2180 
2181  return *this;
2182  }
2183  } /* end of namespace MPI */
2184 
2185 #endif /* DOXYGEN */
2186 
2187 } /* end of namespace TrilinosWrappers */
2188 
2192 namespace internal
2193 {
2194  namespace LinearOperator
2195  {
2196  template <typename> class ReinitHelper;
2197 
2202  template <>
2203  class ReinitHelper<TrilinosWrappers::MPI::Vector>
2204  {
2205  public:
2206  template <typename Matrix>
2207  static
2208  void reinit_range_vector (const Matrix &matrix,
2210  bool omit_zeroing_entries)
2211  {
2212  v.reinit(matrix.locally_owned_range_indices(), matrix.get_mpi_communicator(), omit_zeroing_entries);
2213  }
2214 
2215  template <typename Matrix>
2216  static
2217  void reinit_domain_vector(const Matrix &matrix,
2219  bool omit_zeroing_entries)
2220  {
2221  v.reinit(matrix.locally_owned_domain_indices(), matrix.get_mpi_communicator(), omit_zeroing_entries);
2222  }
2223  };
2224 
2225  } /* namespace LinearOperator */
2226 } /* namespace internal */
2227 
2228 
2229 
2235 template <>
2236 struct is_serial_vector< TrilinosWrappers::MPI::Vector > : std::false_type
2237 {
2238 };
2239 
2240 
2241 DEAL_II_NAMESPACE_CLOSE
2242 
2243 #endif // DEAL_II_WITH_TRILINOS
2244 
2245 /*---------------------------- trilinos_vector.h ---------------------------*/
2246 
2247 #endif // dealii_trilinos_vector_h
Vector & operator=(const TrilinosScalar s)
static::ExceptionBase & ExcTrilinosError(int arg1)
real_type linfty_norm() const
void sadd(const TrilinosScalar s, const Vector &V)
TrilinosScalar max() const
reference operator[](const size_type index)
void reinit(const Vector &v, const bool omit_zeroing_entries=false, const bool allow_different_maps=false)
void add(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
real_type l1_norm() const
void equ(const TrilinosScalar a, const Vector &V)
std::pair< size_type, size_type > local_range() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:398
TrilinosScalar min() const
size_type size() const
Definition: index_set.h:1553
static void reinit_domain_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
std::shared_ptr< Epetra_MultiVector > nonlocal_vector
std::shared_ptr< Epetra_FEVector > vector
static::ExceptionBase & ExcMessage(std::string arg1)
numbers::NumberTraits< Number >::real_type real_type
Definition: vector.h:128
bool operator!=(const Vector &v) const
reference operator()(const size_type index)
void set(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:593
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:337
real_type l2_norm() const
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void compress(::VectorOperation::values operation)
#define DeclException0(Exception0)
Definition: exceptions.h:570
IndexSet locally_owned_elements() const
std::size_t size() const
Vector & operator/=(const TrilinosScalar factor)
real_type norm_sqr() const
const Epetra_Map & vector_partitioner() const
void swap(Vector &u, Vector &v)
const MPI_Comm & get_mpi_communicator() const
friend class internal::VectorReference
void import_nonlocal_data_for_fe(const ::TrilinosWrappers::SparseMatrix &matrix, const Vector &vector)
TrilinosScalar mean_value() const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:549
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
static void reinit_range_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
static::ExceptionBase & ExcGhostsPresent()
std::size_t memory_consumption() const
double norm(const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Du)
Definition: divergence.h:534
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< TrilinosScalar > &values) const
real_type lp_norm(const TrilinosScalar p) const
bool operator==(const Vector &v) const
TrilinosScalar operator*(const Vector &vec) const
bool in_local_range(const size_type index) const
void scale(const Vector &scaling_factors)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:629
Vector & operator-=(const Vector &V)
size_type local_size() const
static::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
TrilinosScalar add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W)
Vector & operator+=(const Vector &V)
static::ExceptionBase & ExcDifferentParallelPartitioning()
Vector & operator*=(const TrilinosScalar factor)
static::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertIsFinite(number)
Definition: exceptions.h:1219
const Epetra_MultiVector & trilinos_vector() const