Reference documentation for deal.II version GIT 58febcd5cf 2023-09-30 20:00:01+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 - 2023 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_stub.h>
27 
28 # include <deal.II/lac/exceptions.h>
29 # include <deal.II/lac/read_vector.h>
30 # include <deal.II/lac/vector.h>
33 
34 # include <Epetra_ConfigDefs.h>
35 # include <Epetra_FEVector.h>
36 # include <Epetra_LocalMap.h>
37 # include <Epetra_Map.h>
38 # include <Epetra_MpiComm.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, public ReadVector<TrilinosScalar>
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 
509  Vector(Vector &&v); // NOLINT
510 
514  ~Vector() override = default;
515 
520  void
521  clear();
522 
546  void
547  reinit(const Vector &v,
548  const bool omit_zeroing_entries = false,
549  const bool allow_different_maps = false);
550 
573  void
574  reinit(const IndexSet &parallel_partitioning,
575  const MPI_Comm communicator = MPI_COMM_WORLD,
576  const bool omit_zeroing_entries = false);
577 
612  void
613  reinit(const IndexSet &locally_owned_entries,
614  const IndexSet &locally_relevant_or_ghost_entries,
615  const MPI_Comm communicator = MPI_COMM_WORLD,
616  const bool vector_writable = false);
617 
628  void
629  reinit(
630  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
631  const bool make_ghosted = true,
632  const bool vector_writable = false);
633 
637  void
638  reinit(const BlockVector &v, const bool import_data = false);
639 
656  void
658 
671  Vector &
673 
679  Vector &
680  operator=(const Vector &v);
681 
686  Vector &
687  operator=(Vector &&v) noexcept;
688 
696  template <typename Number>
697  Vector &
698  operator=(const ::Vector<Number> &v);
699 
717  void
720  const Vector &vector);
721 
728  void
730  const VectorOperation::values operation);
731 
736  void
738  const VectorOperation::values operation)
739  {
740  import_elements(rwv, operation);
741  }
742 
743 
749  bool
750  operator==(const Vector &v) const;
751 
757  bool
758  operator!=(const Vector &v) const;
759 
763  size_type
764  size() const override;
765 
770  size_type
772 
795  std::pair<size_type, size_type>
796  local_range() const;
797 
805  bool
806  in_local_range(const size_type index) const;
807 
821  IndexSet
823 
831  bool
833 
840  void
842 
848  operator*(const Vector &vec) const;
849 
853  real_type
854  norm_sqr() const;
855 
860  mean_value() const;
861 
866  min() const;
867 
872  max() const;
873 
877  real_type
878  l1_norm() const;
879 
884  real_type
885  l2_norm() const;
886 
891  real_type
892  lp_norm(const TrilinosScalar p) const;
893 
897  real_type
898  linfty_norm() const;
899 
920  add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
921 
927  bool
928  all_zero() const;
929 
935  bool
936  is_non_negative() const;
952  reference
953  operator()(const size_type index);
954 
963  operator()(const size_type index) const;
964 
970  reference
971  operator[](const size_type index);
972 
979  operator[](const size_type index) const;
980 
996  void
997  extract_subvector_to(const std::vector<size_type> &indices,
998  std::vector<TrilinosScalar> &values) const;
999 
1003  virtual void
1005  ArrayView<TrilinosScalar> &elements) const override;
1006 
1034  template <typename ForwardIterator, typename OutputIterator>
1035  void
1036  extract_subvector_to(ForwardIterator indices_begin,
1037  const ForwardIterator indices_end,
1038  OutputIterator values_begin) const;
1039 
1048  iterator
1050 
1056  begin() const;
1057 
1062  iterator
1063  end();
1064 
1070  end() const;
1071 
1086  void
1087  set(const std::vector<size_type> &indices,
1088  const std::vector<TrilinosScalar> &values);
1089 
1094  void
1095  set(const std::vector<size_type> &indices,
1096  const ::Vector<TrilinosScalar> &values);
1097 
1103  void
1104  set(const size_type n_elements,
1105  const size_type *indices,
1106  const TrilinosScalar *values);
1107 
1112  void
1113  add(const std::vector<size_type> &indices,
1114  const std::vector<TrilinosScalar> &values);
1115 
1120  void
1121  add(const std::vector<size_type> &indices,
1122  const ::Vector<TrilinosScalar> &values);
1123 
1129  void
1130  add(const size_type n_elements,
1131  const size_type *indices,
1132  const TrilinosScalar *values);
1133 
1137  Vector &
1139 
1143  Vector &
1145 
1149  Vector &
1151 
1155  Vector &
1157 
1162  void
1164 
1177  void
1178  add(const Vector &V, const bool allow_different_maps = false);
1179 
1183  void
1184  add(const TrilinosScalar a, const Vector &V);
1185 
1189  void
1191  const Vector &V,
1192  const TrilinosScalar b,
1193  const Vector &W);
1194 
1199  void
1200  sadd(const TrilinosScalar s, const Vector &V);
1201 
1205  void
1206  sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1207 
1213  void
1214  scale(const Vector &scaling_factors);
1215 
1219  void
1220  equ(const TrilinosScalar a, const Vector &V);
1232  const Epetra_MultiVector &
1234 
1239  Epetra_FEVector &
1241 
1246  const Epetra_BlockMap &
1248 
1256  void
1257  print(std::ostream &out,
1258  const unsigned int precision = 3,
1259  const bool scientific = true,
1260  const bool across = true) const;
1261 
1275  void
1276  swap(Vector &v);
1277 
1281  std::size_t
1282  memory_consumption() const;
1283 
1287  MPI_Comm
1295 
1300  int,
1301  << "An error with error number " << arg1
1302  << " occurred while calling a Trilinos function");
1303 
1309  size_type,
1310  size_type,
1311  size_type,
1312  size_type,
1313  << "You are trying to access element " << arg1
1314  << " of a distributed vector, but this element is not stored "
1315  << "on the current processor. Note: There are " << arg2
1316  << " elements stored "
1317  << "on the current processor from within the range [" << arg3 << ','
1318  << arg4 << "] but Trilinos vectors need not store contiguous "
1319  << "ranges on each processor, and not every element in "
1320  << "this range may in fact be stored locally."
1321  << "\n\n"
1322  << "A common source for this kind of problem is that you "
1323  << "are passing a 'fully distributed' vector into a function "
1324  << "that needs read access to vector elements that correspond "
1325  << "to degrees of freedom on ghost cells (or at least to "
1326  << "'locally active' degrees of freedom that are not also "
1327  << "'locally owned'). You need to pass a vector that has these "
1328  << "elements as ghost entries.");
1329 
1330  private:
1342  Epetra_CombineMode last_action;
1343 
1349 
1355 
1361  std::unique_ptr<Epetra_FEVector> vector;
1362 
1368  std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1369 
1374 
1375  // Make the reference class a friend.
1377  };
1378 
1379 
1380 
1381  // ------------------- inline and template functions --------------
1382 
1383 
1391  inline void
1393  {
1394  u.swap(v);
1395  }
1396  } // namespace MPI
1397 
1398 # ifndef DOXYGEN
1399 
1400  namespace internal
1401  {
1402  inline VectorReference::VectorReference(MPI::Vector &vector,
1403  const size_type index)
1404  : vector(vector)
1405  , index(index)
1406  {}
1407 
1408 
1409  inline const VectorReference &
1410  VectorReference::operator=(const VectorReference &r) const
1411  {
1412  // as explained in the class
1413  // documentation, this is not the copy
1414  // operator. so simply pass on to the
1415  // "correct" assignment operator
1416  *this = static_cast<TrilinosScalar>(r);
1417 
1418  return *this;
1419  }
1420 
1421 
1422 
1423  inline VectorReference &
1424  VectorReference::operator=(const VectorReference &r)
1425  {
1426  // as above
1427  *this = static_cast<TrilinosScalar>(r);
1428 
1429  return *this;
1430  }
1431 
1432 
1433  inline const VectorReference &
1434  VectorReference::operator=(const TrilinosScalar &value) const
1435  {
1436  vector.set(1, &index, &value);
1437  return *this;
1438  }
1439 
1440 
1441 
1442  inline const VectorReference &
1443  VectorReference::operator+=(const TrilinosScalar &value) const
1444  {
1445  vector.add(1, &index, &value);
1446  return *this;
1447  }
1448 
1449 
1450 
1451  inline const VectorReference &
1452  VectorReference::operator-=(const TrilinosScalar &value) const
1453  {
1454  TrilinosScalar new_value = -value;
1455  vector.add(1, &index, &new_value);
1456  return *this;
1457  }
1458 
1459 
1460 
1461  inline const VectorReference &
1462  VectorReference::operator*=(const TrilinosScalar &value) const
1463  {
1464  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1465  vector.set(1, &index, &new_value);
1466  return *this;
1467  }
1468 
1469 
1470 
1471  inline const VectorReference &
1472  VectorReference::operator/=(const TrilinosScalar &value) const
1473  {
1474  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1475  vector.set(1, &index, &new_value);
1476  return *this;
1477  }
1478  } // namespace internal
1479 
1480  namespace MPI
1481  {
1482  inline bool
1483  Vector::in_local_range(const size_type index) const
1484  {
1485  std::pair<size_type, size_type> range = local_range();
1486 
1487  return ((index >= range.first) && (index < range.second));
1488  }
1489 
1490 
1491 
1492  inline IndexSet
1494  {
1495  Assert(owned_elements.size() == size(),
1496  ExcMessage(
1497  "The locally owned elements have not been properly initialized!"
1498  " This happens for example if this object has been initialized"
1499  " with exactly one overlapping IndexSet."));
1500  return owned_elements;
1501  }
1502 
1503 
1504 
1505  inline bool
1507  {
1508  return has_ghosts;
1509  }
1510 
1511 
1512 
1513  inline void
1515  {}
1516 
1517 
1518 
1519  inline internal::VectorReference
1520  Vector::operator()(const size_type index)
1521  {
1522  return internal::VectorReference(*this, index);
1523  }
1524 
1525 
1526 
1527  inline internal::VectorReference
1528  Vector::operator[](const size_type index)
1529  {
1530  return operator()(index);
1531  }
1532 
1533 
1534 
1535  inline TrilinosScalar
1536  Vector::operator[](const size_type index) const
1537  {
1538  return operator()(index);
1539  }
1540 
1541 
1542 
1543  inline void
1544  Vector::extract_subvector_to(const std::vector<size_type> &indices,
1545  std::vector<TrilinosScalar> &values) const
1546  {
1547  for (size_type i = 0; i < indices.size(); ++i)
1548  values[i] = operator()(indices[i]);
1549  }
1550 
1551 
1552 
1553  inline void
1555  ArrayView<TrilinosScalar> &elements) const
1556  {
1557  AssertDimension(indices.size(), elements.size());
1558  for (unsigned int i = 0; i < indices.size(); ++i)
1559  {
1560  AssertIndexRange(indices[i], size());
1561  elements[i] = (*this)[indices[i]];
1562  }
1563  }
1564 
1565 
1566 
1567  template <typename ForwardIterator, typename OutputIterator>
1568  inline void
1569  Vector::extract_subvector_to(ForwardIterator indices_begin,
1570  const ForwardIterator indices_end,
1571  OutputIterator values_begin) const
1572  {
1573  while (indices_begin != indices_end)
1574  {
1575  *values_begin = operator()(*indices_begin);
1576  indices_begin++;
1577  values_begin++;
1578  }
1579  }
1580 
1581 
1582 
1583  inline Vector::iterator
1584  Vector::begin()
1585  {
1586  return (*vector)[0];
1587  }
1588 
1589 
1590 
1591  inline Vector::iterator
1592  Vector::end()
1593  {
1594  return (*vector)[0] + locally_owned_size();
1595  }
1596 
1597 
1598 
1599  inline Vector::const_iterator
1600  Vector::begin() const
1601  {
1602  return (*vector)[0];
1603  }
1604 
1605 
1606 
1607  inline Vector::const_iterator
1608  Vector::end() const
1609  {
1610  return (*vector)[0] + locally_owned_size();
1611  }
1612 
1613 
1614 
1615  inline void
1616  Vector::set(const std::vector<size_type> &indices,
1617  const std::vector<TrilinosScalar> &values)
1618  {
1619  // if we have ghost values, do not allow
1620  // writing to this vector at all.
1622 
1623  AssertDimension(indices.size(), values.size());
1624 
1625  set(indices.size(), indices.data(), values.data());
1626  }
1627 
1628 
1629 
1630  inline void
1631  Vector::set(const std::vector<size_type> &indices,
1632  const ::Vector<TrilinosScalar> &values)
1633  {
1634  // if we have ghost values, do not allow
1635  // writing to this vector at all.
1637 
1638  AssertDimension(indices.size(), values.size());
1639 
1640  set(indices.size(), indices.data(), values.begin());
1641  }
1642 
1643 
1644 
1645  inline void
1646  Vector::set(const size_type n_elements,
1647  const size_type *indices,
1648  const TrilinosScalar *values)
1649  {
1650  // if we have ghost values, do not allow
1651  // writing to this vector at all.
1653 
1654  if (last_action == Add)
1655  {
1656  const int ierr = vector->GlobalAssemble(Add);
1657  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1658  }
1659 
1660  if (last_action != Insert)
1661  last_action = Insert;
1662 
1663  for (size_type i = 0; i < n_elements; ++i)
1664  {
1665  const TrilinosWrappers::types::int_type row = indices[i];
1666  const TrilinosWrappers::types::int_type local_row =
1667  vector->Map().LID(row);
1668  if (local_row != -1)
1669  (*vector)[0][local_row] = values[i];
1670  else
1671  {
1672  const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1673  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1674  compressed = false;
1675  }
1676  // in set operation, do not use the pre-allocated vector for nonlocal
1677  // entries even if it exists. This is to ensure that we really only
1678  // set the elements touched by the set() method and not all contained
1679  // in the nonlocal entries vector (there is no way to distinguish them
1680  // on the receiving processor)
1681  }
1682  }
1683 
1684 
1685 
1686  inline void
1687  Vector::add(const std::vector<size_type> &indices,
1688  const std::vector<TrilinosScalar> &values)
1689  {
1690  // if we have ghost values, do not allow
1691  // writing to this vector at all.
1693  AssertDimension(indices.size(), values.size());
1694 
1695  add(indices.size(), indices.data(), values.data());
1696  }
1697 
1698 
1699 
1700  inline void
1701  Vector::add(const std::vector<size_type> &indices,
1702  const ::Vector<TrilinosScalar> &values)
1703  {
1704  // if we have ghost values, do not allow
1705  // writing to this vector at all.
1707  AssertDimension(indices.size(), values.size());
1708 
1709  add(indices.size(), indices.data(), values.begin());
1710  }
1711 
1712 
1713 
1714  inline void
1715  Vector::add(const size_type n_elements,
1716  const size_type *indices,
1717  const TrilinosScalar *values)
1718  {
1719  // if we have ghost values, do not allow
1720  // writing to this vector at all.
1722 
1723  if (last_action != Add)
1724  {
1725  if (last_action == Insert)
1726  {
1727  const int ierr = vector->GlobalAssemble(Insert);
1728  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1729  }
1730  last_action = Add;
1731  }
1732 
1733  for (size_type i = 0; i < n_elements; ++i)
1734  {
1735  const size_type row = indices[i];
1736  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1737  static_cast<TrilinosWrappers::types::int_type>(row));
1738  if (local_row != -1)
1739  (*vector)[0][local_row] += values[i];
1740  else if (nonlocal_vector.get() == nullptr)
1741  {
1742  const int ierr = vector->SumIntoGlobalValues(
1743  1,
1744  reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1745  &row),
1746  &values[i]);
1747  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1748  compressed = false;
1749  }
1750  else
1751  {
1752  // use pre-allocated vector for non-local entries if it exists for
1753  // addition operation
1754  const TrilinosWrappers::types::int_type my_row =
1755  nonlocal_vector->Map().LID(
1756  static_cast<TrilinosWrappers::types::int_type>(row));
1757  Assert(my_row != -1,
1758  ExcMessage(
1759  "Attempted to write into off-processor vector entry "
1760  "that has not be specified as being writable upon "
1761  "initialization"));
1762  (*nonlocal_vector)[0][my_row] += values[i];
1763  compressed = false;
1764  }
1765  }
1766  }
1767 
1768 
1769 
1770  inline Vector::size_type
1771  Vector::size() const
1772  {
1773 # ifndef DEAL_II_WITH_64BIT_INDICES
1774  return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1775 # else
1776  return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1777 # endif
1778  }
1779 
1780 
1781 
1782  inline Vector::size_type
1784  {
1785  return owned_elements.n_elements();
1786  }
1787 
1788 
1789 
1790  inline std::pair<Vector::size_type, Vector::size_type>
1791  Vector::local_range() const
1792  {
1793 # ifndef DEAL_II_WITH_64BIT_INDICES
1794  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1796  vector->Map().MaxMyGID() + 1;
1797 # else
1799  vector->Map().MinMyGID64();
1801  vector->Map().MaxMyGID64() + 1;
1802 # endif
1803 
1804  Assert(
1805  end - begin == vector->Map().NumMyElements(),
1806  ExcMessage(
1807  "This function only makes sense if the elements that this "
1808  "vector stores on the current processor form a contiguous range. "
1809  "This does not appear to be the case for the current vector."));
1810 
1811  return std::make_pair(begin, end);
1812  }
1813 
1814 
1815 
1816  inline TrilinosScalar
1817  Vector::operator*(const Vector &vec) const
1818  {
1819  Assert(vector->Map().SameAs(vec.vector->Map()),
1822 
1823  TrilinosScalar result;
1824 
1825  const int ierr = vector->Dot(*(vec.vector), &result);
1826  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1827 
1828  return result;
1829  }
1830 
1831 
1832 
1833  inline Vector::real_type
1834  Vector::norm_sqr() const
1835  {
1836  const TrilinosScalar d = l2_norm();
1837  return d * d;
1838  }
1839 
1840 
1841 
1842  inline TrilinosScalar
1843  Vector::mean_value() const
1844  {
1846 
1848  const int ierr = vector->MeanValue(&mean);
1849  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1850 
1851  return mean;
1852  }
1853 
1854 
1855 
1856  inline TrilinosScalar
1857  Vector::min() const
1858  {
1859  TrilinosScalar min_value;
1860  const int ierr = vector->MinValue(&min_value);
1861  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1862 
1863  return min_value;
1864  }
1865 
1866 
1867 
1868  inline TrilinosScalar
1869  Vector::max() const
1870  {
1871  TrilinosScalar max_value;
1872  const int ierr = vector->MaxValue(&max_value);
1873  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1874 
1875  return max_value;
1876  }
1877 
1878 
1879 
1880  inline Vector::real_type
1881  Vector::l1_norm() const
1882  {
1884 
1885  TrilinosScalar d;
1886  const int ierr = vector->Norm1(&d);
1887  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1888 
1889  return d;
1890  }
1891 
1892 
1893 
1894  inline Vector::real_type
1895  Vector::l2_norm() const
1896  {
1898 
1899  TrilinosScalar d;
1900  const int ierr = vector->Norm2(&d);
1901  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1902 
1903  return d;
1904  }
1905 
1906 
1907 
1908  inline Vector::real_type
1909  Vector::lp_norm(const TrilinosScalar p) const
1910  {
1912 
1913  TrilinosScalar norm = 0;
1914  TrilinosScalar sum = 0;
1915  const size_type n_local = locally_owned_size();
1916 
1917  // loop over all the elements because
1918  // Trilinos does not support lp norms
1919  for (size_type i = 0; i < n_local; ++i)
1920  sum += std::pow(std::fabs((*vector)[0][i]), p);
1921 
1922  norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1923 
1924  return norm;
1925  }
1926 
1927 
1928 
1929  inline Vector::real_type
1930  Vector::linfty_norm() const
1931  {
1932  // while we disallow the other
1933  // norm operations on ghosted
1934  // vectors, this particular norm
1935  // is safe to run even in the
1936  // presence of ghost elements
1937  TrilinosScalar d;
1938  const int ierr = vector->NormInf(&d);
1939  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1940 
1941  return d;
1942  }
1943 
1944 
1945 
1946  inline TrilinosScalar
1948  const Vector &V,
1949  const Vector &W)
1950  {
1951  this->add(a, V);
1952  return *this * W;
1953  }
1954 
1955 
1956 
1957  // inline also scalar products, vector
1958  // additions etc. since they are all
1959  // representable by a single Trilinos
1960  // call. This reduces the overhead of the
1961  // wrapper class.
1962  inline Vector &
1964  {
1965  AssertIsFinite(a);
1966 
1967  const int ierr = vector->Scale(a);
1968  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1969 
1970  return *this;
1971  }
1972 
1973 
1974 
1975  inline Vector &
1977  {
1978  AssertIsFinite(a);
1979 
1980  const TrilinosScalar factor = 1. / a;
1981 
1982  AssertIsFinite(factor);
1983 
1984  const int ierr = vector->Scale(factor);
1985  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1986 
1987  return *this;
1988  }
1989 
1990 
1991 
1992  inline Vector &
1993  Vector::operator+=(const Vector &v)
1994  {
1995  AssertDimension(size(), v.size());
1996  Assert(vector->Map().SameAs(v.vector->Map()),
1998 
1999  const int ierr = vector->Update(1.0, *(v.vector), 1.0);
2000  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2001 
2002  return *this;
2003  }
2004 
2005 
2006 
2007  inline Vector &
2008  Vector::operator-=(const Vector &v)
2009  {
2010  AssertDimension(size(), v.size());
2011  Assert(vector->Map().SameAs(v.vector->Map()),
2013 
2014  const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
2015  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2016 
2017  return *this;
2018  }
2019 
2020 
2021 
2022  inline void
2023  Vector::add(const TrilinosScalar s)
2024  {
2025  // if we have ghost values, do not allow
2026  // writing to this vector at all.
2028  AssertIsFinite(s);
2029 
2030  size_type n_local = locally_owned_size();
2031  for (size_type i = 0; i < n_local; ++i)
2032  (*vector)[0][i] += s;
2033  }
2034 
2035 
2036 
2037  inline void
2038  Vector::add(const TrilinosScalar a, const Vector &v)
2039  {
2040  // if we have ghost values, do not allow
2041  // writing to this vector at all.
2044 
2045  AssertIsFinite(a);
2046 
2047  const int ierr = vector->Update(a, *(v.vector), 1.);
2048  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2049  }
2050 
2051 
2052 
2053  inline void
2054  Vector::add(const TrilinosScalar a,
2055  const Vector &v,
2056  const TrilinosScalar b,
2057  const Vector &w)
2058  {
2059  // if we have ghost values, do not allow
2060  // writing to this vector at all.
2063  AssertDimension(locally_owned_size(), w.locally_owned_size());
2064 
2065  AssertIsFinite(a);
2066  AssertIsFinite(b);
2067 
2068  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2069 
2070  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2071  }
2072 
2073 
2074 
2075  inline void
2076  Vector::sadd(const TrilinosScalar s, const Vector &v)
2077  {
2078  // if we have ghost values, do not allow
2079  // writing to this vector at all.
2081  AssertDimension(size(), v.size());
2082 
2083  AssertIsFinite(s);
2084 
2085  // We assume that the vectors have the same Map
2086  // if the local size is the same and if the vectors are not ghosted
2087  if (locally_owned_size() == v.locally_owned_size() &&
2088  !v.has_ghost_elements())
2089  {
2090  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2092  const int ierr = vector->Update(1., *(v.vector), s);
2093  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2094  }
2095  else
2096  {
2097  (*this) *= s;
2098  this->add(v, true);
2099  }
2100  }
2101 
2102 
2103 
2104  inline void
2105  Vector::sadd(const TrilinosScalar s,
2106  const TrilinosScalar a,
2107  const Vector &v)
2108  {
2109  // if we have ghost values, do not allow
2110  // writing to this vector at all.
2112  AssertDimension(size(), v.size());
2113  AssertIsFinite(s);
2114  AssertIsFinite(a);
2115 
2116  // We assume that the vectors have the same Map
2117  // if the local size is the same and if the vectors are not ghosted
2118  if (locally_owned_size() == v.locally_owned_size() &&
2119  !v.has_ghost_elements())
2120  {
2121  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2123  const int ierr = vector->Update(a, *(v.vector), s);
2124  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2125  }
2126  else
2127  {
2128  (*this) *= s;
2129  Vector tmp = v;
2130  tmp *= a;
2131  this->add(tmp, true);
2132  }
2133  }
2134 
2135 
2136 
2137  inline void
2138  Vector::scale(const Vector &factors)
2139  {
2140  // if we have ghost values, do not allow
2141  // writing to this vector at all.
2144 
2145  const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2146  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2147  }
2148 
2149 
2150 
2151  inline void
2152  Vector::equ(const TrilinosScalar a, const Vector &v)
2153  {
2154  // if we have ghost values, do not allow
2155  // writing to this vector at all.
2157  AssertIsFinite(a);
2158 
2159  // If we don't have the same map, copy.
2160  if (vector->Map().SameAs(v.vector->Map()) == false)
2161  {
2162  this->sadd(0., a, v);
2163  }
2164  else
2165  {
2166  // Otherwise, just update
2167  int ierr = vector->Update(a, *v.vector, 0.0);
2168  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2169 
2170  last_action = Zero;
2171  }
2172  }
2173 
2174 
2175 
2176  inline const Epetra_MultiVector &
2177  Vector::trilinos_vector() const
2178  {
2179  return static_cast<const Epetra_MultiVector &>(*vector);
2180  }
2181 
2182 
2183 
2184  inline Epetra_FEVector &
2186  {
2187  return *vector;
2188  }
2189 
2190 
2191 
2192  inline const Epetra_BlockMap &
2194  {
2195  return vector->Map();
2196  }
2197 
2198 
2199 
2200  inline MPI_Comm
2202  {
2203  const Epetra_MpiComm *mpi_comm =
2204  dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2205  return mpi_comm->Comm();
2206  }
2207 
2208 
2209 
2210  template <typename number>
2211  Vector::Vector(const IndexSet &parallel_partitioner,
2212  const ::Vector<number> &v,
2213  const MPI_Comm communicator)
2214  {
2215  *this =
2216  Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2217  owned_elements = parallel_partitioner;
2218  }
2219 
2220 
2221 
2222  inline Vector &
2224  {
2225  AssertIsFinite(s);
2226 
2227  int ierr = vector->PutScalar(s);
2228  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2229 
2230  if (nonlocal_vector.get() != nullptr)
2231  {
2232  ierr = nonlocal_vector->PutScalar(0.);
2233  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2234  }
2235 
2236  return *this;
2237  }
2238  } /* end of namespace MPI */
2239 
2240 # endif /* DOXYGEN */
2241 
2242 } /* end of namespace TrilinosWrappers */
2243 
2247 namespace internal
2248 {
2249  namespace LinearOperatorImplementation
2250  {
2251  template <typename>
2252  class ReinitHelper;
2253 
2258  template <>
2260  {
2261  public:
2262  template <typename Matrix>
2263  static void
2266  bool omit_zeroing_entries)
2267  {
2268  v.reinit(matrix.locally_owned_range_indices(),
2269  matrix.get_mpi_communicator(),
2270  omit_zeroing_entries);
2271  }
2272 
2273  template <typename Matrix>
2274  static void
2277  bool omit_zeroing_entries)
2278  {
2279  v.reinit(matrix.locally_owned_domain_indices(),
2280  matrix.get_mpi_communicator(),
2281  omit_zeroing_entries);
2282  }
2283  };
2284 
2285  } // namespace LinearOperatorImplementation
2286 } /* namespace internal */
2287 
2288 
2289 
2293 template <>
2295 {};
2296 
2297 
2299 
2300 #endif
2301 
2302 #endif
std::size_t size() const
Definition: array_view.h:573
size_type size() const
Definition: index_set.h:1696
size_type n_elements() const
Definition: index_set.h:1854
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:840
void compress(VectorOperation::values operation)
Vector(const IndexSet &parallel_partitioning, const ::Vector< Number > &v, const MPI_Comm communicator=MPI_COMM_WORLD)
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)
const Epetra_BlockMap & trilinos_partitioner() const
void add(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
std::unique_ptr< Epetra_MultiVector > nonlocal_vector
Epetra_FEVector & trilinos_vector()
void add(const TrilinosScalar s)
void import_elements(const LinearAlgebra::ReadWriteVector< double > &rwv, const VectorOperation::values operation)
MPI_Comm get_mpi_communicator() const
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
size_type size() const override
bool in_local_range(const size_type index) 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
virtual void extract_subvector_to(const ArrayView< const size_type > &indices, ArrayView< TrilinosScalar > &elements) const override
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)
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)
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)
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
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)
void swap(Vector &u, Vector &v)
TrilinosScalar max() const
const_iterator end() const
Definition: vector.h:110
typename numbers::NumberTraits< Number >::real_type real_type
Definition: vector.h:147
bool has_ghost_elements() const
const value_type * const_iterator
Definition: vector.h:133
types::global_dof_index size_type
Definition: vector.h:136
virtual size_type size() const override
value_type * iterator
Definition: vector.h:132
size_type locally_owned_size() const
static void reinit_domain_vector(const Matrix &matrix, TrilinosWrappers::MPI::Vector &v, bool omit_zeroing_entries)
static void reinit_range_vector(const Matrix &matrix, TrilinosWrappers::MPI::Vector &v, bool omit_zeroing_entries)
#define DEAL_II_DEPRECATED
Definition: config.h:178
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define DeclException0(Exception0)
Definition: exceptions.h:467
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:581
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define AssertIsFinite(number)
Definition: exceptions.h:1884
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
static ::ExceptionBase & ExcDifferentParallelPartitioning()
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:512
static ::ExceptionBase & ExcGhostsPresent()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
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)
int gid(const Epetra_BlockMap &map, int i)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int global_dof_index
Definition: types.h:82
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
double TrilinosScalar
Definition: types.h:175