Reference documentation for deal.II version GIT 5dcc62ab57 2022-07-04 21:05:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_vector.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2022 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.md at
12 // the top level directory of deal.II.
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/index_set.h>
24 # include <deal.II/base/mpi.h>
26 # include <deal.II/base/utilities.h>
27 
28 # include <deal.II/lac/exceptions.h>
29 # include <deal.II/lac/vector.h>
32 
33 # include <Epetra_ConfigDefs.h>
34 # include <Epetra_FEVector.h>
35 # include <Epetra_LocalMap.h>
36 # include <Epetra_Map.h>
37 # include <Epetra_MpiComm.h>
38 # include <mpi.h>
39 
40 # include <memory>
41 # include <utility>
42 # include <vector>
43 
45 
46 // Forward declarations
47 # ifndef DOXYGEN
48 namespace LinearAlgebra
49 {
50  // Forward declaration
51  template <typename Number>
52  class ReadWriteVector;
53 } // namespace LinearAlgebra
54 # endif
55 
66 namespace TrilinosWrappers
67 {
68  class SparseMatrix;
69 
80  namespace internal
81  {
86 
96  class VectorReference
97  {
98  private:
103  VectorReference(MPI::Vector &vector, const size_type index);
104 
105  public:
109  VectorReference(const VectorReference &) = default;
110 
122  const VectorReference &
123  operator=(const VectorReference &r) const;
124 
128  VectorReference &
129  operator=(const VectorReference &r);
130 
134  const VectorReference &
135  operator=(const TrilinosScalar &s) const;
136 
140  const VectorReference &
141  operator+=(const TrilinosScalar &s) const;
142 
146  const VectorReference &
147  operator-=(const TrilinosScalar &s) const;
148 
152  const VectorReference &
153  operator*=(const TrilinosScalar &s) const;
154 
158  const VectorReference &
159  operator/=(const TrilinosScalar &s) const;
160 
165  operator TrilinosScalar() const;
166 
171  int,
172  << "An error with error number " << arg1
173  << " occurred while calling a Trilinos function");
174 
175  private:
179  MPI::Vector &vector;
180 
184  const size_type index;
185 
186  // Make the vector class a friend, so that it can create objects of the
187  // present type.
189  };
190  } // namespace internal
195 # ifndef DEAL_II_WITH_64BIT_INDICES
196  // define a helper function that queries the global ID of local ID of
197  // an Epetra_BlockMap object by calling either the 32- or 64-bit
198  // function necessary.
199  inline int
200  gid(const Epetra_BlockMap &map, int i)
201  {
202  return map.GID(i);
203  }
204 # else
205  // define a helper function that queries the global ID of local ID of
206  // an Epetra_BlockMap object by calling either the 32- or 64-bit
207  // function necessary.
208  inline long long int
209  gid(const Epetra_BlockMap &map, int i)
210  {
211  return map.GID64(i);
212  }
213 # endif
214 
220  namespace MPI
221  {
222  class BlockVector;
223 
398  class Vector : public Subscriptor
399  {
400  public:
409  using iterator = value_type *;
410  using const_iterator = const value_type *;
413 
423  Vector();
424 
428  Vector(const Vector &v);
429 
448  explicit Vector(const IndexSet &parallel_partitioning,
449  const MPI_Comm &communicator = MPI_COMM_WORLD);
450 
462  Vector(const IndexSet &local,
463  const IndexSet &ghost,
464  const MPI_Comm &communicator = MPI_COMM_WORLD);
465 
480  Vector(const IndexSet &parallel_partitioning,
481  const Vector & v,
482  const MPI_Comm &communicator = MPI_COMM_WORLD);
483 
496  template <typename Number>
497  Vector(const IndexSet & parallel_partitioning,
498  const ::Vector<Number> &v,
499  const MPI_Comm & communicator = MPI_COMM_WORLD);
500 
505  Vector(Vector &&v) noexcept;
506 
510  ~Vector() override = default;
511 
516  void
517  clear();
518 
542  void
543  reinit(const Vector &v,
544  const bool omit_zeroing_entries = false,
545  const bool allow_different_maps = false);
546 
569  void
570  reinit(const IndexSet &parallel_partitioning,
571  const MPI_Comm &communicator = MPI_COMM_WORLD,
572  const bool omit_zeroing_entries = false);
573 
599  void
600  reinit(const IndexSet &locally_owned_entries,
601  const IndexSet &ghost_entries,
602  const MPI_Comm &communicator = MPI_COMM_WORLD,
603  const bool vector_writable = false);
604 
608  void
609  reinit(const BlockVector &v, const bool import_data = false);
610 
627  void
628  compress(::VectorOperation::values operation);
629 
642  Vector &
644 
650  Vector &
651  operator=(const Vector &v);
652 
657  Vector &
658  operator=(Vector &&v) noexcept;
659 
667  template <typename Number>
668  Vector &
669  operator=(const ::Vector<Number> &v);
670 
688  void
691  const Vector & vector);
692 
699  void
700  import(const LinearAlgebra::ReadWriteVector<double> &rwv,
701  const VectorOperation::values operation);
702 
703 
709  bool
710  operator==(const Vector &v) const;
711 
717  bool
718  operator!=(const Vector &v) const;
719 
723  size_type
724  size() const;
725 
740  size_type
741  local_size() const;
742 
747  size_type
749 
771  std::pair<size_type, size_type>
772  local_range() const;
773 
781  bool
782  in_local_range(const size_type index) const;
783 
797  IndexSet
799 
807  bool
809 
816  void
818 
824  operator*(const Vector &vec) const;
825 
829  real_type
830  norm_sqr() const;
831 
836  mean_value() const;
837 
842  min() const;
843 
848  max() const;
849 
853  real_type
854  l1_norm() const;
855 
860  real_type
861  l2_norm() const;
862 
867  real_type
868  lp_norm(const TrilinosScalar p) const;
869 
873  real_type
874  linfty_norm() const;
875 
896  add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
897 
903  bool
904  all_zero() const;
905 
911  bool
912  is_non_negative() const;
914 
915 
920 
928  reference
929  operator()(const size_type index);
930 
939  operator()(const size_type index) const;
940 
946  reference
947  operator[](const size_type index);
948 
955  operator[](const size_type index) const;
956 
972  void
973  extract_subvector_to(const std::vector<size_type> &indices,
974  std::vector<TrilinosScalar> & values) const;
975 
1003  template <typename ForwardIterator, typename OutputIterator>
1004  void
1005  extract_subvector_to(ForwardIterator indices_begin,
1006  const ForwardIterator indices_end,
1007  OutputIterator values_begin) const;
1008 
1019  iterator
1021 
1027  begin() const;
1028 
1033  iterator
1034  end();
1035 
1041  end() const;
1042 
1044 
1045 
1050 
1057  void
1058  set(const std::vector<size_type> & indices,
1059  const std::vector<TrilinosScalar> &values);
1060 
1065  void
1066  set(const std::vector<size_type> & indices,
1067  const ::Vector<TrilinosScalar> &values);
1068 
1074  void
1075  set(const size_type n_elements,
1076  const size_type * indices,
1077  const TrilinosScalar *values);
1078 
1083  void
1084  add(const std::vector<size_type> & indices,
1085  const std::vector<TrilinosScalar> &values);
1086 
1091  void
1092  add(const std::vector<size_type> & indices,
1093  const ::Vector<TrilinosScalar> &values);
1094 
1100  void
1101  add(const size_type n_elements,
1102  const size_type * indices,
1103  const TrilinosScalar *values);
1104 
1108  Vector &
1110 
1114  Vector &
1116 
1120  Vector &
1122 
1126  Vector &
1128 
1133  void
1135 
1148  void
1149  add(const Vector &V, const bool allow_different_maps = false);
1150 
1154  void
1155  add(const TrilinosScalar a, const Vector &V);
1156 
1160  void
1162  const Vector & V,
1163  const TrilinosScalar b,
1164  const Vector & W);
1165 
1170  void
1171  sadd(const TrilinosScalar s, const Vector &V);
1172 
1176  void
1177  sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1178 
1184  void
1185  scale(const Vector &scaling_factors);
1186 
1190  void
1191  equ(const TrilinosScalar a, const Vector &V);
1193 
1198 
1203  const Epetra_MultiVector &
1205 
1210  Epetra_FEVector &
1212 
1217  const Epetra_BlockMap &
1219 
1227  void
1228  print(std::ostream & out,
1229  const unsigned int precision = 3,
1230  const bool scientific = true,
1231  const bool across = true) const;
1232 
1246  void
1247  swap(Vector &v);
1248 
1252  std::size_t
1253  memory_consumption() const;
1254 
1259  const MPI_Comm &
1262 
1267 
1272  int,
1273  << "An error with error number " << arg1
1274  << " occurred while calling a Trilinos function");
1275 
1281  size_type,
1282  size_type,
1283  size_type,
1284  size_type,
1285  << "You are trying to access element " << arg1
1286  << " of a distributed vector, but this element is not stored "
1287  << "on the current processor. Note: There are " << arg2
1288  << " elements stored "
1289  << "on the current processor from within the range [" << arg3 << ','
1290  << arg4 << "] but Trilinos vectors need not store contiguous "
1291  << "ranges on each processor, and not every element in "
1292  << "this range may in fact be stored locally."
1293  << "\n\n"
1294  << "A common source for this kind of problem is that you "
1295  << "are passing a 'fully distributed' vector into a function "
1296  << "that needs read access to vector elements that correspond "
1297  << "to degrees of freedom on ghost cells (or at least to "
1298  << "'locally active' degrees of freedom that are not also "
1299  << "'locally owned'). You need to pass a vector that has these "
1300  << "elements as ghost entries.");
1301 
1302  private:
1314  Epetra_CombineMode last_action;
1315 
1321 
1327 
1333  std::unique_ptr<Epetra_FEVector> vector;
1334 
1340  std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1341 
1346 
1347  // Make the reference class a friend.
1349  };
1350 
1351 
1352 
1353  // ------------------- inline and template functions --------------
1354 
1355 
1363  inline void
1365  {
1366  u.swap(v);
1367  }
1368  } // namespace MPI
1369 
1370 # ifndef DOXYGEN
1371 
1372  namespace internal
1373  {
1374  inline VectorReference::VectorReference(MPI::Vector & vector,
1375  const size_type index)
1376  : vector(vector)
1377  , index(index)
1378  {}
1379 
1380 
1381  inline const VectorReference &
1382  VectorReference::operator=(const VectorReference &r) const
1383  {
1384  // as explained in the class
1385  // documentation, this is not the copy
1386  // operator. so simply pass on to the
1387  // "correct" assignment operator
1388  *this = static_cast<TrilinosScalar>(r);
1389 
1390  return *this;
1391  }
1392 
1393 
1394 
1395  inline VectorReference &
1396  VectorReference::operator=(const VectorReference &r)
1397  {
1398  // as above
1399  *this = static_cast<TrilinosScalar>(r);
1400 
1401  return *this;
1402  }
1403 
1404 
1405  inline const VectorReference &
1406  VectorReference::operator=(const TrilinosScalar &value) const
1407  {
1408  vector.set(1, &index, &value);
1409  return *this;
1410  }
1411 
1412 
1413 
1414  inline const VectorReference &
1415  VectorReference::operator+=(const TrilinosScalar &value) const
1416  {
1417  vector.add(1, &index, &value);
1418  return *this;
1419  }
1420 
1421 
1422 
1423  inline const VectorReference &
1424  VectorReference::operator-=(const TrilinosScalar &value) const
1425  {
1426  TrilinosScalar new_value = -value;
1427  vector.add(1, &index, &new_value);
1428  return *this;
1429  }
1430 
1431 
1432 
1433  inline const VectorReference &
1434  VectorReference::operator*=(const TrilinosScalar &value) const
1435  {
1436  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1437  vector.set(1, &index, &new_value);
1438  return *this;
1439  }
1440 
1441 
1442 
1443  inline const VectorReference &
1444  VectorReference::operator/=(const TrilinosScalar &value) const
1445  {
1446  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1447  vector.set(1, &index, &new_value);
1448  return *this;
1449  }
1450  } // namespace internal
1451 
1452  namespace MPI
1453  {
1454  inline bool
1455  Vector::in_local_range(const size_type index) const
1456  {
1457  std::pair<size_type, size_type> range = local_range();
1458 
1459  return ((index >= range.first) && (index < range.second));
1460  }
1461 
1462 
1463 
1464  inline IndexSet
1466  {
1467  Assert(owned_elements.size() == size(),
1468  ExcMessage(
1469  "The locally owned elements have not been properly initialized!"
1470  " This happens for example if this object has been initialized"
1471  " with exactly one overlapping IndexSet."));
1472  return owned_elements;
1473  }
1474 
1475 
1476 
1477  inline bool
1479  {
1480  return has_ghosts;
1481  }
1482 
1483 
1484 
1485  inline void
1487  {}
1488 
1489 
1490 
1491  inline internal::VectorReference
1492  Vector::operator()(const size_type index)
1493  {
1494  return internal::VectorReference(*this, index);
1495  }
1496 
1497 
1498 
1499  inline internal::VectorReference
1500  Vector::operator[](const size_type index)
1501  {
1502  return operator()(index);
1503  }
1504 
1505 
1506 
1507  inline TrilinosScalar
1508  Vector::operator[](const size_type index) const
1509  {
1510  return operator()(index);
1511  }
1512 
1513 
1514 
1515  inline void
1516  Vector::extract_subvector_to(const std::vector<size_type> &indices,
1517  std::vector<TrilinosScalar> & values) const
1518  {
1519  for (size_type i = 0; i < indices.size(); ++i)
1520  values[i] = operator()(indices[i]);
1521  }
1522 
1523 
1524 
1525  template <typename ForwardIterator, typename OutputIterator>
1526  inline void
1527  Vector::extract_subvector_to(ForwardIterator indices_begin,
1528  const ForwardIterator indices_end,
1529  OutputIterator values_begin) const
1530  {
1531  while (indices_begin != indices_end)
1532  {
1533  *values_begin = operator()(*indices_begin);
1534  indices_begin++;
1535  values_begin++;
1536  }
1537  }
1538 
1539 
1540 
1541  inline Vector::iterator
1542  Vector::begin()
1543  {
1544  return (*vector)[0];
1545  }
1546 
1547 
1548 
1549  inline Vector::iterator
1550  Vector::end()
1551  {
1552  return (*vector)[0] + locally_owned_size();
1553  }
1554 
1555 
1556 
1557  inline Vector::const_iterator
1558  Vector::begin() const
1559  {
1560  return (*vector)[0];
1561  }
1562 
1563 
1564 
1565  inline Vector::const_iterator
1566  Vector::end() const
1567  {
1568  return (*vector)[0] + locally_owned_size();
1569  }
1570 
1571 
1572 
1573  inline void
1574  Vector::set(const std::vector<size_type> & indices,
1575  const std::vector<TrilinosScalar> &values)
1576  {
1577  // if we have ghost values, do not allow
1578  // writing to this vector at all.
1580 
1581  AssertDimension(indices.size(), values.size());
1582 
1583  set(indices.size(), indices.data(), values.data());
1584  }
1585 
1586 
1587 
1588  inline void
1589  Vector::set(const std::vector<size_type> & indices,
1590  const ::Vector<TrilinosScalar> &values)
1591  {
1592  // if we have ghost values, do not allow
1593  // writing to this vector at all.
1595 
1596  AssertDimension(indices.size(), values.size());
1597 
1598  set(indices.size(), indices.data(), values.begin());
1599  }
1600 
1601 
1602 
1603  inline void
1604  Vector::set(const size_type n_elements,
1605  const size_type * indices,
1606  const TrilinosScalar *values)
1607  {
1608  // if we have ghost values, do not allow
1609  // writing to this vector at all.
1611 
1612  if (last_action == Add)
1613  {
1614  const int ierr = vector->GlobalAssemble(Add);
1615  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1616  }
1617 
1618  if (last_action != Insert)
1619  last_action = Insert;
1620 
1621  for (size_type i = 0; i < n_elements; ++i)
1622  {
1623  const TrilinosWrappers::types::int_type row = indices[i];
1624  const TrilinosWrappers::types::int_type local_row =
1625  vector->Map().LID(row);
1626  if (local_row != -1)
1627  (*vector)[0][local_row] = values[i];
1628  else
1629  {
1630  const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1631  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1632  compressed = false;
1633  }
1634  // in set operation, do not use the pre-allocated vector for nonlocal
1635  // entries even if it exists. This is to ensure that we really only
1636  // set the elements touched by the set() method and not all contained
1637  // in the nonlocal entries vector (there is no way to distinguish them
1638  // on the receiving processor)
1639  }
1640  }
1641 
1642 
1643 
1644  inline void
1645  Vector::add(const std::vector<size_type> & indices,
1646  const std::vector<TrilinosScalar> &values)
1647  {
1648  // if we have ghost values, do not allow
1649  // writing to this vector at all.
1651  AssertDimension(indices.size(), values.size());
1652 
1653  add(indices.size(), indices.data(), values.data());
1654  }
1655 
1656 
1657 
1658  inline void
1659  Vector::add(const std::vector<size_type> & indices,
1660  const ::Vector<TrilinosScalar> &values)
1661  {
1662  // if we have ghost values, do not allow
1663  // writing to this vector at all.
1665  AssertDimension(indices.size(), values.size());
1666 
1667  add(indices.size(), indices.data(), values.begin());
1668  }
1669 
1670 
1671 
1672  inline void
1673  Vector::add(const size_type n_elements,
1674  const size_type * indices,
1675  const TrilinosScalar *values)
1676  {
1677  // if we have ghost values, do not allow
1678  // writing to this vector at all.
1680 
1681  if (last_action != Add)
1682  {
1683  if (last_action == Insert)
1684  {
1685  const int ierr = vector->GlobalAssemble(Insert);
1686  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1687  }
1688  last_action = Add;
1689  }
1690 
1691  for (size_type i = 0; i < n_elements; ++i)
1692  {
1693  const size_type row = indices[i];
1694  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1695  static_cast<TrilinosWrappers::types::int_type>(row));
1696  if (local_row != -1)
1697  (*vector)[0][local_row] += values[i];
1698  else if (nonlocal_vector.get() == nullptr)
1699  {
1700  const int ierr = vector->SumIntoGlobalValues(
1701  1,
1702  reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1703  &row),
1704  &values[i]);
1705  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1706  compressed = false;
1707  }
1708  else
1709  {
1710  // use pre-allocated vector for non-local entries if it exists for
1711  // addition operation
1712  const TrilinosWrappers::types::int_type my_row =
1713  nonlocal_vector->Map().LID(
1714  static_cast<TrilinosWrappers::types::int_type>(row));
1715  Assert(my_row != -1,
1716  ExcMessage(
1717  "Attempted to write into off-processor vector entry "
1718  "that has not be specified as being writable upon "
1719  "initialization"));
1720  (*nonlocal_vector)[0][my_row] += values[i];
1721  compressed = false;
1722  }
1723  }
1724  }
1725 
1726 
1727 
1728  inline Vector::size_type
1729  Vector::size() const
1730  {
1731 # ifndef DEAL_II_WITH_64BIT_INDICES
1732  return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1733 # else
1734  return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1735 # endif
1736  }
1737 
1738 
1739 
1740  inline Vector::size_type
1741  Vector::local_size() const
1742  {
1743  return vector->Map().NumMyElements();
1744  }
1745 
1746 
1747 
1748  inline Vector::size_type
1750  {
1751  return owned_elements.n_elements();
1752  }
1753 
1754 
1755 
1756  inline std::pair<Vector::size_type, Vector::size_type>
1757  Vector::local_range() const
1758  {
1759 # ifndef DEAL_II_WITH_64BIT_INDICES
1760  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1762  vector->Map().MaxMyGID() + 1;
1763 # else
1765  vector->Map().MinMyGID64();
1767  vector->Map().MaxMyGID64() + 1;
1768 # endif
1769 
1770  Assert(
1771  end - begin == vector->Map().NumMyElements(),
1772  ExcMessage(
1773  "This function only makes sense if the elements that this "
1774  "vector stores on the current processor form a contiguous range. "
1775  "This does not appear to be the case for the current vector."));
1776 
1777  return std::make_pair(begin, end);
1778  }
1779 
1780 
1781 
1782  inline TrilinosScalar
1783  Vector::operator*(const Vector &vec) const
1784  {
1785  Assert(vector->Map().SameAs(vec.vector->Map()),
1788 
1789  TrilinosScalar result;
1790 
1791  const int ierr = vector->Dot(*(vec.vector), &result);
1792  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1793 
1794  return result;
1795  }
1796 
1797 
1798 
1799  inline Vector::real_type
1800  Vector::norm_sqr() const
1801  {
1802  const TrilinosScalar d = l2_norm();
1803  return d * d;
1804  }
1805 
1806 
1807 
1808  inline TrilinosScalar
1809  Vector::mean_value() const
1810  {
1812 
1814  const int ierr = vector->MeanValue(&mean);
1815  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1816 
1817  return mean;
1818  }
1819 
1820 
1821 
1822  inline TrilinosScalar
1823  Vector::min() const
1824  {
1825  TrilinosScalar min_value;
1826  const int ierr = vector->MinValue(&min_value);
1827  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1828 
1829  return min_value;
1830  }
1831 
1832 
1833 
1834  inline TrilinosScalar
1835  Vector::max() const
1836  {
1837  TrilinosScalar max_value;
1838  const int ierr = vector->MaxValue(&max_value);
1839  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1840 
1841  return max_value;
1842  }
1843 
1844 
1845 
1846  inline Vector::real_type
1847  Vector::l1_norm() const
1848  {
1850 
1851  TrilinosScalar d;
1852  const int ierr = vector->Norm1(&d);
1853  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1854 
1855  return d;
1856  }
1857 
1858 
1859 
1860  inline Vector::real_type
1861  Vector::l2_norm() const
1862  {
1864 
1865  TrilinosScalar d;
1866  const int ierr = vector->Norm2(&d);
1867  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1868 
1869  return d;
1870  }
1871 
1872 
1873 
1874  inline Vector::real_type
1875  Vector::lp_norm(const TrilinosScalar p) const
1876  {
1878 
1879  TrilinosScalar norm = 0;
1880  TrilinosScalar sum = 0;
1881  const size_type n_local = locally_owned_size();
1882 
1883  // loop over all the elements because
1884  // Trilinos does not support lp norms
1885  for (size_type i = 0; i < n_local; ++i)
1886  sum += std::pow(std::fabs((*vector)[0][i]), p);
1887 
1888  norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1889 
1890  return norm;
1891  }
1892 
1893 
1894 
1895  inline Vector::real_type
1896  Vector::linfty_norm() const
1897  {
1898  // while we disallow the other
1899  // norm operations on ghosted
1900  // vectors, this particular norm
1901  // is safe to run even in the
1902  // presence of ghost elements
1903  TrilinosScalar d;
1904  const int ierr = vector->NormInf(&d);
1905  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1906 
1907  return d;
1908  }
1909 
1910 
1911 
1912  inline TrilinosScalar
1914  const Vector & V,
1915  const Vector & W)
1916  {
1917  this->add(a, V);
1918  return *this * W;
1919  }
1920 
1921 
1922 
1923  // inline also scalar products, vector
1924  // additions etc. since they are all
1925  // representable by a single Trilinos
1926  // call. This reduces the overhead of the
1927  // wrapper class.
1928  inline Vector &
1930  {
1931  AssertIsFinite(a);
1932 
1933  const int ierr = vector->Scale(a);
1934  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1935 
1936  return *this;
1937  }
1938 
1939 
1940 
1941  inline Vector &
1943  {
1944  AssertIsFinite(a);
1945 
1946  const TrilinosScalar factor = 1. / a;
1947 
1948  AssertIsFinite(factor);
1949 
1950  const int ierr = vector->Scale(factor);
1951  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1952 
1953  return *this;
1954  }
1955 
1956 
1957 
1958  inline Vector &
1959  Vector::operator+=(const Vector &v)
1960  {
1961  AssertDimension(size(), v.size());
1962  Assert(vector->Map().SameAs(v.vector->Map()),
1964 
1965  const int ierr = vector->Update(1.0, *(v.vector), 1.0);
1966  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1967 
1968  return *this;
1969  }
1970 
1971 
1972 
1973  inline Vector &
1974  Vector::operator-=(const Vector &v)
1975  {
1976  AssertDimension(size(), v.size());
1977  Assert(vector->Map().SameAs(v.vector->Map()),
1979 
1980  const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
1981  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1982 
1983  return *this;
1984  }
1985 
1986 
1987 
1988  inline void
1989  Vector::add(const TrilinosScalar s)
1990  {
1991  // if we have ghost values, do not allow
1992  // writing to this vector at all.
1994  AssertIsFinite(s);
1995 
1996  size_type n_local = locally_owned_size();
1997  for (size_type i = 0; i < n_local; ++i)
1998  (*vector)[0][i] += s;
1999  }
2000 
2001 
2002 
2003  inline void
2004  Vector::add(const TrilinosScalar a, const Vector &v)
2005  {
2006  // if we have ghost values, do not allow
2007  // writing to this vector at all.
2010 
2011  AssertIsFinite(a);
2012 
2013  const int ierr = vector->Update(a, *(v.vector), 1.);
2014  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2015  }
2016 
2017 
2018 
2019  inline void
2020  Vector::add(const TrilinosScalar a,
2021  const Vector & v,
2022  const TrilinosScalar b,
2023  const Vector & w)
2024  {
2025  // if we have ghost values, do not allow
2026  // writing to this vector at all.
2029  AssertDimension(locally_owned_size(), w.locally_owned_size());
2030 
2031  AssertIsFinite(a);
2032  AssertIsFinite(b);
2033 
2034  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2035 
2036  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2037  }
2038 
2039 
2040 
2041  inline void
2042  Vector::sadd(const TrilinosScalar s, const Vector &v)
2043  {
2044  // if we have ghost values, do not allow
2045  // writing to this vector at all.
2047  AssertDimension(size(), v.size());
2048 
2049  AssertIsFinite(s);
2050 
2051  // We assume that the vectors have the same Map
2052  // if the local size is the same and if the vectors are not ghosted
2053  if (locally_owned_size() == v.locally_owned_size() &&
2054  !v.has_ghost_elements())
2055  {
2056  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2058  const int ierr = vector->Update(1., *(v.vector), s);
2059  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2060  }
2061  else
2062  {
2063  (*this) *= s;
2064  this->add(v, true);
2065  }
2066  }
2067 
2068 
2069 
2070  inline void
2071  Vector::sadd(const TrilinosScalar s,
2072  const TrilinosScalar a,
2073  const Vector & v)
2074  {
2075  // if we have ghost values, do not allow
2076  // writing to this vector at all.
2078  AssertDimension(size(), v.size());
2079  AssertIsFinite(s);
2080  AssertIsFinite(a);
2081 
2082  // We assume that the vectors have the same Map
2083  // if the local size is the same and if the vectors are not ghosted
2084  if (locally_owned_size() == v.locally_owned_size() &&
2085  !v.has_ghost_elements())
2086  {
2087  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2089  const int ierr = vector->Update(a, *(v.vector), s);
2090  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2091  }
2092  else
2093  {
2094  (*this) *= s;
2095  Vector tmp = v;
2096  tmp *= a;
2097  this->add(tmp, true);
2098  }
2099  }
2100 
2101 
2102 
2103  inline void
2104  Vector::scale(const Vector &factors)
2105  {
2106  // if we have ghost values, do not allow
2107  // writing to this vector at all.
2110 
2111  const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2112  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2113  }
2114 
2115 
2116 
2117  inline void
2118  Vector::equ(const TrilinosScalar a, const Vector &v)
2119  {
2120  // if we have ghost values, do not allow
2121  // writing to this vector at all.
2123  AssertIsFinite(a);
2124 
2125  // If we don't have the same map, copy.
2126  if (vector->Map().SameAs(v.vector->Map()) == false)
2127  {
2128  this->sadd(0., a, v);
2129  }
2130  else
2131  {
2132  // Otherwise, just update
2133  int ierr = vector->Update(a, *v.vector, 0.0);
2134  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2135 
2136  last_action = Zero;
2137  }
2138  }
2139 
2140 
2141 
2142  inline const Epetra_MultiVector &
2143  Vector::trilinos_vector() const
2144  {
2145  return static_cast<const Epetra_MultiVector &>(*vector);
2146  }
2147 
2148 
2149 
2150  inline Epetra_FEVector &
2152  {
2153  return *vector;
2154  }
2155 
2156 
2157 
2158  inline const Epetra_BlockMap &
2160  {
2161  return vector->Map();
2162  }
2163 
2164 
2165 
2166  inline const MPI_Comm &
2168  {
2169  static MPI_Comm comm;
2170 
2171  const Epetra_MpiComm *mpi_comm =
2172  dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2173  comm = mpi_comm->Comm();
2174 
2175  return comm;
2176  }
2177 
2178  template <typename number>
2179  Vector::Vector(const IndexSet & parallel_partitioner,
2180  const ::Vector<number> &v,
2181  const MPI_Comm & communicator)
2182  {
2183  *this =
2184  Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2185  owned_elements = parallel_partitioner;
2186  }
2187 
2188 
2189 
2190  inline Vector &
2192  {
2193  AssertIsFinite(s);
2194 
2195  int ierr = vector->PutScalar(s);
2196  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2197 
2198  if (nonlocal_vector.get() != nullptr)
2199  {
2200  ierr = nonlocal_vector->PutScalar(0.);
2201  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2202  }
2203 
2204  return *this;
2205  }
2206  } /* end of namespace MPI */
2207 
2208 # endif /* DOXYGEN */
2209 
2210 } /* end of namespace TrilinosWrappers */
2211 
2215 namespace internal
2216 {
2217  namespace LinearOperatorImplementation
2218  {
2219  template <typename>
2220  class ReinitHelper;
2221 
2226  template <>
2228  {
2229  public:
2230  template <typename Matrix>
2231  static void
2232  reinit_range_vector(const Matrix & matrix,
2234  bool omit_zeroing_entries)
2235  {
2236  v.reinit(matrix.locally_owned_range_indices(),
2237  matrix.get_mpi_communicator(),
2238  omit_zeroing_entries);
2239  }
2240 
2241  template <typename Matrix>
2242  static void
2245  bool omit_zeroing_entries)
2246  {
2247  v.reinit(matrix.locally_owned_domain_indices(),
2248  matrix.get_mpi_communicator(),
2249  omit_zeroing_entries);
2250  }
2251  };
2252 
2253  } // namespace LinearOperatorImplementation
2254 } /* namespace internal */
2255 
2256 
2257 
2261 template <>
2263 {};
2264 
2265 
2267 
2268 #endif // DEAL_II_WITH_TRILINOS
2269 
2270 /*---------------------------- trilinos_vector.h ---------------------------*/
2271 
2272 #endif // dealii_trilinos_vector_h
size_type size() const
Definition: index_set.h:1635
size_type n_elements() const
Definition: index_set.h:1783
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:777
Definition: vector.h:109
#define DEAL_II_DEPRECATED
Definition: config.h:164
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DeclException0(Exception0)
Definition: exceptions.h:464
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:578
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertIsFinite(number)
Definition: exceptions.h:1758
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcGhostsPresent()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
void sadd(const TrilinosScalar s, const Vector &V)
TrilinosScalar mean_value() const
void add(const size_type n_elements, const size_type *indices, const TrilinosScalar *values)
static ::ExceptionBase & ExcTrilinosError(int arg1)
const Epetra_BlockMap & trilinos_partitioner() const
void add(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
const MPI_Comm & get_mpi_communicator() const
static void reinit_domain_vector(const Matrix &matrix, TrilinosWrappers::MPI::Vector &v, bool omit_zeroing_entries)
std::unique_ptr< Epetra_MultiVector > nonlocal_vector
Epetra_FEVector & trilinos_vector()
void add(const TrilinosScalar s)
int gid(const Epetra_BlockMap &map, int i)
reference operator()(const size_type index)
void import_nonlocal_data_for_fe(const ::TrilinosWrappers::SparseMatrix &matrix, const Vector &vector)
std::unique_ptr< Epetra_FEVector > vector
bool in_local_range(const size_type index) const
size_type local_size() const
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
friend class internal::VectorReference
real_type lp_norm(const TrilinosScalar p) const
IndexSet locally_owned_elements() const
const internal::VectorReference const_reference
void reinit(const Vector &v, const bool omit_zeroing_entries=false, const bool allow_different_maps=false)
Vector & operator+=(const Vector &V)
void compress(::VectorOperation::values operation)
real_type norm_sqr() const
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< TrilinosScalar > &values) const
bool operator!=(const Vector &v) const
~Vector() override=default
const_iterator begin() const
real_type linfty_norm() const
internal::VectorReference reference
TrilinosScalar min() const
void set(const std::vector< size_type > &indices, const ::Vector< TrilinosScalar > &values)
void set(const size_type n_elements, const size_type *indices, const TrilinosScalar *values)
static void reinit_range_vector(const Matrix &matrix, TrilinosWrappers::MPI::Vector &v, bool omit_zeroing_entries)
void scale(const Vector &scaling_factors)
void equ(const TrilinosScalar a, const Vector &V)
const Epetra_MultiVector & trilinos_vector() const
bool operator==(const Vector &v) const
Vector & operator/=(const TrilinosScalar factor)
void sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V)
reference operator[](const size_type index)
static ::ExceptionBase & ExcDifferentParallelPartitioning()
size_type locally_owned_size() const
std::pair< size_type, size_type > local_range() const
TrilinosScalar operator[](const size_type index) const
Vector & operator*=(const TrilinosScalar factor)
Vector & operator=(const TrilinosScalar s)
Vector & operator=(const ::Vector< Number > &v)
TrilinosScalar add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W)
::types::global_dof_index size_type
TrilinosScalar operator*(const Vector &vec) const
void add(const std::vector< size_type > &indices, const ::Vector< TrilinosScalar > &values)
Vector & operator-=(const Vector &V)
void extract_subvector_to(ForwardIterator indices_begin, const ForwardIterator indices_end, OutputIterator values_begin) const
std::size_t memory_consumption() const
Vector(const IndexSet &parallel_partitioning, const ::Vector< Number > &v, const MPI_Comm &communicator=MPI_COMM_WORLD)
void add(const TrilinosScalar a, const Vector &V, const TrilinosScalar b, const Vector &W)
void add(const TrilinosScalar a, const Vector &V)
void set(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
void swap(Vector &u, Vector &v)
TrilinosScalar max() const
const_iterator end() const
typename numbers::NumberTraits< Number >::real_type real_type
Definition: vector.h:139
bool has_ghost_elements() const
const value_type * const_iterator
Definition: vector.h:125
types::global_dof_index size_type
Definition: vector.h:128
size_type size() const
value_type * iterator
Definition: vector.h:124
size_type locally_owned_size() const
Expression fabs(const Expression &x)
@ matrix
Contents is actually a matrix.
static const char V
SparseMatrix< double > SparseMatrix
types::global_dof_index size_type
Definition: cuda_kernels.h:45
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int global_dof_index
Definition: types.h:76
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const MPI_Comm & comm
double TrilinosScalar
Definition: types.h:163