Reference documentation for deal.II version GIT b8135fa6eb 2023-03-25 15:55: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_stub.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 
39 # include <memory>
40 # include <utility>
41 # include <vector>
42 
44 
45 // Forward declarations
46 # ifndef DOXYGEN
47 namespace LinearAlgebra
48 {
49  // Forward declaration
50  template <typename Number>
51  class ReadWriteVector;
52 } // namespace LinearAlgebra
53 # endif
54 
65 namespace TrilinosWrappers
66 {
67  class SparseMatrix;
68 
79  namespace internal
80  {
85 
95  class VectorReference
96  {
97  private:
102  VectorReference(MPI::Vector &vector, const size_type index);
103 
104  public:
108  VectorReference(const VectorReference &) = default;
109 
121  const VectorReference &
122  operator=(const VectorReference &r) const;
123 
127  VectorReference &
128  operator=(const VectorReference &r);
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 
157  const VectorReference &
158  operator/=(const TrilinosScalar &s) const;
159 
164  operator TrilinosScalar() const;
165 
170  int,
171  << "An error with error number " << arg1
172  << " occurred while calling a Trilinos function");
173 
174  private:
178  MPI::Vector &vector;
179 
183  const size_type index;
184 
185  // Make the vector class a friend, so that it can create objects of the
186  // present type.
188  };
189  } // namespace internal
194 # ifndef DEAL_II_WITH_64BIT_INDICES
195  // define a helper function that queries the global ID of local ID of
196  // an Epetra_BlockMap object by calling either the 32- or 64-bit
197  // function necessary.
198  inline int
199  gid(const Epetra_BlockMap &map, int i)
200  {
201  return map.GID(i);
202  }
203 # else
204  // define a helper function that queries the global ID of local ID of
205  // an Epetra_BlockMap object by calling either the 32- or 64-bit
206  // function necessary.
207  inline long long int
208  gid(const Epetra_BlockMap &map, int i)
209  {
210  return map.GID64(i);
211  }
212 # endif
213 
219  namespace MPI
220  {
221  class BlockVector;
222 
397  class Vector : public Subscriptor
398  {
399  public:
408  using iterator = value_type *;
409  using const_iterator = const value_type *;
412 
422  Vector();
423 
427  Vector(const Vector &v);
428 
447  explicit Vector(const IndexSet &parallel_partitioning,
448  const MPI_Comm &communicator = MPI_COMM_WORLD);
449 
461  Vector(const IndexSet &local,
462  const IndexSet &ghost,
463  const MPI_Comm &communicator = MPI_COMM_WORLD);
464 
479  Vector(const IndexSet &parallel_partitioning,
480  const Vector & v,
481  const MPI_Comm &communicator = MPI_COMM_WORLD);
482 
495  template <typename Number>
496  Vector(const IndexSet & parallel_partitioning,
497  const ::Vector<Number> &v,
498  const MPI_Comm & communicator = MPI_COMM_WORLD);
499 
508  Vector(Vector &&v); // NOLINT
509 
513  ~Vector() override = default;
514 
519  void
520  clear();
521 
545  void
546  reinit(const Vector &v,
547  const bool omit_zeroing_entries = false,
548  const bool allow_different_maps = false);
549 
572  void
573  reinit(const IndexSet &parallel_partitioning,
574  const MPI_Comm &communicator = MPI_COMM_WORLD,
575  const bool omit_zeroing_entries = false);
576 
602  void
603  reinit(const IndexSet &locally_owned_entries,
604  const IndexSet &ghost_entries,
605  const MPI_Comm &communicator = MPI_COMM_WORLD,
606  const bool vector_writable = false);
607 
612  void
613  reinit(
614  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
615  const bool vector_writable = false);
616 
620  void
621  reinit(const BlockVector &v, const bool import_data = false);
622 
639  void
641 
654  Vector &
656 
662  Vector &
663  operator=(const Vector &v);
664 
669  Vector &
670  operator=(Vector &&v) noexcept;
671 
679  template <typename Number>
680  Vector &
681  operator=(const ::Vector<Number> &v);
682 
700  void
703  const Vector & vector);
704 
711  void
712  import(const LinearAlgebra::ReadWriteVector<double> &rwv,
713  const VectorOperation::values operation);
714 
715 
721  bool
722  operator==(const Vector &v) const;
723 
729  bool
730  operator!=(const Vector &v) const;
731 
735  size_type
736  size() const;
737 
752  size_type
753  local_size() const;
754 
759  size_type
761 
783  std::pair<size_type, size_type>
784  local_range() const;
785 
793  bool
794  in_local_range(const size_type index) const;
795 
809  IndexSet
811 
819  bool
821 
828  void
830 
836  operator*(const Vector &vec) const;
837 
841  real_type
842  norm_sqr() const;
843 
848  mean_value() const;
849 
854  min() const;
855 
860  max() const;
861 
865  real_type
866  l1_norm() const;
867 
872  real_type
873  l2_norm() const;
874 
879  real_type
880  lp_norm(const TrilinosScalar p) const;
881 
885  real_type
886  linfty_norm() const;
887 
908  add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
909 
915  bool
916  all_zero() const;
917 
923  bool
924  is_non_negative() const;
940  reference
941  operator()(const size_type index);
942 
951  operator()(const size_type index) const;
952 
958  reference
959  operator[](const size_type index);
960 
967  operator[](const size_type index) const;
968 
984  void
985  extract_subvector_to(const std::vector<size_type> &indices,
986  std::vector<TrilinosScalar> & values) const;
987 
1015  template <typename ForwardIterator, typename OutputIterator>
1016  void
1017  extract_subvector_to(ForwardIterator indices_begin,
1018  const ForwardIterator indices_end,
1019  OutputIterator values_begin) const;
1020 
1031  iterator
1033 
1039  begin() const;
1040 
1045  iterator
1046  end();
1047 
1053  end() const;
1054 
1069  void
1070  set(const std::vector<size_type> & indices,
1071  const std::vector<TrilinosScalar> &values);
1072 
1077  void
1078  set(const std::vector<size_type> & indices,
1079  const ::Vector<TrilinosScalar> &values);
1080 
1086  void
1087  set(const size_type n_elements,
1088  const size_type * indices,
1089  const TrilinosScalar *values);
1090 
1095  void
1096  add(const std::vector<size_type> & indices,
1097  const std::vector<TrilinosScalar> &values);
1098 
1103  void
1104  add(const std::vector<size_type> & indices,
1105  const ::Vector<TrilinosScalar> &values);
1106 
1112  void
1113  add(const size_type n_elements,
1114  const size_type * indices,
1115  const TrilinosScalar *values);
1116 
1120  Vector &
1122 
1126  Vector &
1128 
1132  Vector &
1134 
1138  Vector &
1140 
1145  void
1147 
1160  void
1161  add(const Vector &V, const bool allow_different_maps = false);
1162 
1166  void
1167  add(const TrilinosScalar a, const Vector &V);
1168 
1172  void
1174  const Vector & V,
1175  const TrilinosScalar b,
1176  const Vector & W);
1177 
1182  void
1183  sadd(const TrilinosScalar s, const Vector &V);
1184 
1188  void
1189  sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1190 
1196  void
1197  scale(const Vector &scaling_factors);
1198 
1202  void
1203  equ(const TrilinosScalar a, const Vector &V);
1215  const Epetra_MultiVector &
1217 
1222  Epetra_FEVector &
1224 
1229  const Epetra_BlockMap &
1231 
1239  void
1240  print(std::ostream & out,
1241  const unsigned int precision = 3,
1242  const bool scientific = true,
1243  const bool across = true) const;
1244 
1258  void
1259  swap(Vector &v);
1260 
1264  std::size_t
1265  memory_consumption() const;
1266 
1271  const MPI_Comm &
1279 
1284  int,
1285  << "An error with error number " << arg1
1286  << " occurred while calling a Trilinos function");
1287 
1293  size_type,
1294  size_type,
1295  size_type,
1296  size_type,
1297  << "You are trying to access element " << arg1
1298  << " of a distributed vector, but this element is not stored "
1299  << "on the current processor. Note: There are " << arg2
1300  << " elements stored "
1301  << "on the current processor from within the range [" << arg3 << ','
1302  << arg4 << "] but Trilinos vectors need not store contiguous "
1303  << "ranges on each processor, and not every element in "
1304  << "this range may in fact be stored locally."
1305  << "\n\n"
1306  << "A common source for this kind of problem is that you "
1307  << "are passing a 'fully distributed' vector into a function "
1308  << "that needs read access to vector elements that correspond "
1309  << "to degrees of freedom on ghost cells (or at least to "
1310  << "'locally active' degrees of freedom that are not also "
1311  << "'locally owned'). You need to pass a vector that has these "
1312  << "elements as ghost entries.");
1313 
1314  private:
1326  Epetra_CombineMode last_action;
1327 
1333 
1339 
1345  std::unique_ptr<Epetra_FEVector> vector;
1346 
1352  std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1353 
1358 
1359  // Make the reference class a friend.
1361  };
1362 
1363 
1364 
1365  // ------------------- inline and template functions --------------
1366 
1367 
1375  inline void
1377  {
1378  u.swap(v);
1379  }
1380  } // namespace MPI
1381 
1382 # ifndef DOXYGEN
1383 
1384  namespace internal
1385  {
1386  inline VectorReference::VectorReference(MPI::Vector & vector,
1387  const size_type index)
1388  : vector(vector)
1389  , index(index)
1390  {}
1391 
1392 
1393  inline const VectorReference &
1394  VectorReference::operator=(const VectorReference &r) const
1395  {
1396  // as explained in the class
1397  // documentation, this is not the copy
1398  // operator. so simply pass on to the
1399  // "correct" assignment operator
1400  *this = static_cast<TrilinosScalar>(r);
1401 
1402  return *this;
1403  }
1404 
1405 
1406 
1407  inline VectorReference &
1408  VectorReference::operator=(const VectorReference &r)
1409  {
1410  // as above
1411  *this = static_cast<TrilinosScalar>(r);
1412 
1413  return *this;
1414  }
1415 
1416 
1417  inline const VectorReference &
1418  VectorReference::operator=(const TrilinosScalar &value) const
1419  {
1420  vector.set(1, &index, &value);
1421  return *this;
1422  }
1423 
1424 
1425 
1426  inline const VectorReference &
1427  VectorReference::operator+=(const TrilinosScalar &value) const
1428  {
1429  vector.add(1, &index, &value);
1430  return *this;
1431  }
1432 
1433 
1434 
1435  inline const VectorReference &
1436  VectorReference::operator-=(const TrilinosScalar &value) const
1437  {
1438  TrilinosScalar new_value = -value;
1439  vector.add(1, &index, &new_value);
1440  return *this;
1441  }
1442 
1443 
1444 
1445  inline const VectorReference &
1446  VectorReference::operator*=(const TrilinosScalar &value) const
1447  {
1448  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1449  vector.set(1, &index, &new_value);
1450  return *this;
1451  }
1452 
1453 
1454 
1455  inline const VectorReference &
1456  VectorReference::operator/=(const TrilinosScalar &value) const
1457  {
1458  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1459  vector.set(1, &index, &new_value);
1460  return *this;
1461  }
1462  } // namespace internal
1463 
1464  namespace MPI
1465  {
1466  inline bool
1467  Vector::in_local_range(const size_type index) const
1468  {
1469  std::pair<size_type, size_type> range = local_range();
1470 
1471  return ((index >= range.first) && (index < range.second));
1472  }
1473 
1474 
1475 
1476  inline IndexSet
1478  {
1479  Assert(owned_elements.size() == size(),
1480  ExcMessage(
1481  "The locally owned elements have not been properly initialized!"
1482  " This happens for example if this object has been initialized"
1483  " with exactly one overlapping IndexSet."));
1484  return owned_elements;
1485  }
1486 
1487 
1488 
1489  inline bool
1491  {
1492  return has_ghosts;
1493  }
1494 
1495 
1496 
1497  inline void
1499  {}
1500 
1501 
1502 
1503  inline internal::VectorReference
1504  Vector::operator()(const size_type index)
1505  {
1506  return internal::VectorReference(*this, index);
1507  }
1508 
1509 
1510 
1511  inline internal::VectorReference
1512  Vector::operator[](const size_type index)
1513  {
1514  return operator()(index);
1515  }
1516 
1517 
1518 
1519  inline TrilinosScalar
1520  Vector::operator[](const size_type index) const
1521  {
1522  return operator()(index);
1523  }
1524 
1525 
1526 
1527  inline void
1528  Vector::extract_subvector_to(const std::vector<size_type> &indices,
1529  std::vector<TrilinosScalar> & values) const
1530  {
1531  for (size_type i = 0; i < indices.size(); ++i)
1532  values[i] = operator()(indices[i]);
1533  }
1534 
1535 
1536 
1537  template <typename ForwardIterator, typename OutputIterator>
1538  inline void
1539  Vector::extract_subvector_to(ForwardIterator indices_begin,
1540  const ForwardIterator indices_end,
1541  OutputIterator values_begin) const
1542  {
1543  while (indices_begin != indices_end)
1544  {
1545  *values_begin = operator()(*indices_begin);
1546  indices_begin++;
1547  values_begin++;
1548  }
1549  }
1550 
1551 
1552 
1553  inline Vector::iterator
1554  Vector::begin()
1555  {
1556  return (*vector)[0];
1557  }
1558 
1559 
1560 
1561  inline Vector::iterator
1562  Vector::end()
1563  {
1564  return (*vector)[0] + locally_owned_size();
1565  }
1566 
1567 
1568 
1569  inline Vector::const_iterator
1570  Vector::begin() const
1571  {
1572  return (*vector)[0];
1573  }
1574 
1575 
1576 
1577  inline Vector::const_iterator
1578  Vector::end() const
1579  {
1580  return (*vector)[0] + locally_owned_size();
1581  }
1582 
1583 
1584 
1585  inline void
1586  Vector::set(const std::vector<size_type> & indices,
1587  const std::vector<TrilinosScalar> &values)
1588  {
1589  // if we have ghost values, do not allow
1590  // writing to this vector at all.
1592 
1593  AssertDimension(indices.size(), values.size());
1594 
1595  set(indices.size(), indices.data(), values.data());
1596  }
1597 
1598 
1599 
1600  inline void
1601  Vector::set(const std::vector<size_type> & indices,
1602  const ::Vector<TrilinosScalar> &values)
1603  {
1604  // if we have ghost values, do not allow
1605  // writing to this vector at all.
1607 
1608  AssertDimension(indices.size(), values.size());
1609 
1610  set(indices.size(), indices.data(), values.begin());
1611  }
1612 
1613 
1614 
1615  inline void
1616  Vector::set(const size_type n_elements,
1617  const size_type * indices,
1618  const TrilinosScalar *values)
1619  {
1620  // if we have ghost values, do not allow
1621  // writing to this vector at all.
1623 
1624  if (last_action == Add)
1625  {
1626  const int ierr = vector->GlobalAssemble(Add);
1627  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1628  }
1629 
1630  if (last_action != Insert)
1631  last_action = Insert;
1632 
1633  for (size_type i = 0; i < n_elements; ++i)
1634  {
1635  const TrilinosWrappers::types::int_type row = indices[i];
1636  const TrilinosWrappers::types::int_type local_row =
1637  vector->Map().LID(row);
1638  if (local_row != -1)
1639  (*vector)[0][local_row] = values[i];
1640  else
1641  {
1642  const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1643  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1644  compressed = false;
1645  }
1646  // in set operation, do not use the pre-allocated vector for nonlocal
1647  // entries even if it exists. This is to ensure that we really only
1648  // set the elements touched by the set() method and not all contained
1649  // in the nonlocal entries vector (there is no way to distinguish them
1650  // on the receiving processor)
1651  }
1652  }
1653 
1654 
1655 
1656  inline void
1657  Vector::add(const std::vector<size_type> & indices,
1658  const std::vector<TrilinosScalar> &values)
1659  {
1660  // if we have ghost values, do not allow
1661  // writing to this vector at all.
1663  AssertDimension(indices.size(), values.size());
1664 
1665  add(indices.size(), indices.data(), values.data());
1666  }
1667 
1668 
1669 
1670  inline void
1671  Vector::add(const std::vector<size_type> & indices,
1672  const ::Vector<TrilinosScalar> &values)
1673  {
1674  // if we have ghost values, do not allow
1675  // writing to this vector at all.
1677  AssertDimension(indices.size(), values.size());
1678 
1679  add(indices.size(), indices.data(), values.begin());
1680  }
1681 
1682 
1683 
1684  inline void
1685  Vector::add(const size_type n_elements,
1686  const size_type * indices,
1687  const TrilinosScalar *values)
1688  {
1689  // if we have ghost values, do not allow
1690  // writing to this vector at all.
1692 
1693  if (last_action != Add)
1694  {
1695  if (last_action == Insert)
1696  {
1697  const int ierr = vector->GlobalAssemble(Insert);
1698  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1699  }
1700  last_action = Add;
1701  }
1702 
1703  for (size_type i = 0; i < n_elements; ++i)
1704  {
1705  const size_type row = indices[i];
1706  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1707  static_cast<TrilinosWrappers::types::int_type>(row));
1708  if (local_row != -1)
1709  (*vector)[0][local_row] += values[i];
1710  else if (nonlocal_vector.get() == nullptr)
1711  {
1712  const int ierr = vector->SumIntoGlobalValues(
1713  1,
1714  reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1715  &row),
1716  &values[i]);
1717  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1718  compressed = false;
1719  }
1720  else
1721  {
1722  // use pre-allocated vector for non-local entries if it exists for
1723  // addition operation
1724  const TrilinosWrappers::types::int_type my_row =
1725  nonlocal_vector->Map().LID(
1726  static_cast<TrilinosWrappers::types::int_type>(row));
1727  Assert(my_row != -1,
1728  ExcMessage(
1729  "Attempted to write into off-processor vector entry "
1730  "that has not be specified as being writable upon "
1731  "initialization"));
1732  (*nonlocal_vector)[0][my_row] += values[i];
1733  compressed = false;
1734  }
1735  }
1736  }
1737 
1738 
1739 
1740  inline Vector::size_type
1741  Vector::size() const
1742  {
1743 # ifndef DEAL_II_WITH_64BIT_INDICES
1744  return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1745 # else
1746  return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1747 # endif
1748  }
1749 
1750 
1751 
1752  inline Vector::size_type
1753  Vector::local_size() const
1754  {
1755  return vector->Map().NumMyElements();
1756  }
1757 
1758 
1759 
1760  inline Vector::size_type
1762  {
1763  return owned_elements.n_elements();
1764  }
1765 
1766 
1767 
1768  inline std::pair<Vector::size_type, Vector::size_type>
1769  Vector::local_range() const
1770  {
1771 # ifndef DEAL_II_WITH_64BIT_INDICES
1772  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1774  vector->Map().MaxMyGID() + 1;
1775 # else
1777  vector->Map().MinMyGID64();
1779  vector->Map().MaxMyGID64() + 1;
1780 # endif
1781 
1782  Assert(
1783  end - begin == vector->Map().NumMyElements(),
1784  ExcMessage(
1785  "This function only makes sense if the elements that this "
1786  "vector stores on the current processor form a contiguous range. "
1787  "This does not appear to be the case for the current vector."));
1788 
1789  return std::make_pair(begin, end);
1790  }
1791 
1792 
1793 
1794  inline TrilinosScalar
1795  Vector::operator*(const Vector &vec) const
1796  {
1797  Assert(vector->Map().SameAs(vec.vector->Map()),
1800 
1801  TrilinosScalar result;
1802 
1803  const int ierr = vector->Dot(*(vec.vector), &result);
1804  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1805 
1806  return result;
1807  }
1808 
1809 
1810 
1811  inline Vector::real_type
1812  Vector::norm_sqr() const
1813  {
1814  const TrilinosScalar d = l2_norm();
1815  return d * d;
1816  }
1817 
1818 
1819 
1820  inline TrilinosScalar
1821  Vector::mean_value() const
1822  {
1824 
1826  const int ierr = vector->MeanValue(&mean);
1827  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1828 
1829  return mean;
1830  }
1831 
1832 
1833 
1834  inline TrilinosScalar
1835  Vector::min() const
1836  {
1837  TrilinosScalar min_value;
1838  const int ierr = vector->MinValue(&min_value);
1839  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1840 
1841  return min_value;
1842  }
1843 
1844 
1845 
1846  inline TrilinosScalar
1847  Vector::max() const
1848  {
1849  TrilinosScalar max_value;
1850  const int ierr = vector->MaxValue(&max_value);
1851  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1852 
1853  return max_value;
1854  }
1855 
1856 
1857 
1858  inline Vector::real_type
1859  Vector::l1_norm() const
1860  {
1862 
1863  TrilinosScalar d;
1864  const int ierr = vector->Norm1(&d);
1865  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1866 
1867  return d;
1868  }
1869 
1870 
1871 
1872  inline Vector::real_type
1873  Vector::l2_norm() const
1874  {
1876 
1877  TrilinosScalar d;
1878  const int ierr = vector->Norm2(&d);
1879  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1880 
1881  return d;
1882  }
1883 
1884 
1885 
1886  inline Vector::real_type
1887  Vector::lp_norm(const TrilinosScalar p) const
1888  {
1890 
1891  TrilinosScalar norm = 0;
1892  TrilinosScalar sum = 0;
1893  const size_type n_local = locally_owned_size();
1894 
1895  // loop over all the elements because
1896  // Trilinos does not support lp norms
1897  for (size_type i = 0; i < n_local; ++i)
1898  sum += std::pow(std::fabs((*vector)[0][i]), p);
1899 
1900  norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1901 
1902  return norm;
1903  }
1904 
1905 
1906 
1907  inline Vector::real_type
1908  Vector::linfty_norm() const
1909  {
1910  // while we disallow the other
1911  // norm operations on ghosted
1912  // vectors, this particular norm
1913  // is safe to run even in the
1914  // presence of ghost elements
1915  TrilinosScalar d;
1916  const int ierr = vector->NormInf(&d);
1917  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1918 
1919  return d;
1920  }
1921 
1922 
1923 
1924  inline TrilinosScalar
1926  const Vector & V,
1927  const Vector & W)
1928  {
1929  this->add(a, V);
1930  return *this * W;
1931  }
1932 
1933 
1934 
1935  // inline also scalar products, vector
1936  // additions etc. since they are all
1937  // representable by a single Trilinos
1938  // call. This reduces the overhead of the
1939  // wrapper class.
1940  inline Vector &
1942  {
1943  AssertIsFinite(a);
1944 
1945  const int ierr = vector->Scale(a);
1946  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1947 
1948  return *this;
1949  }
1950 
1951 
1952 
1953  inline Vector &
1955  {
1956  AssertIsFinite(a);
1957 
1958  const TrilinosScalar factor = 1. / a;
1959 
1960  AssertIsFinite(factor);
1961 
1962  const int ierr = vector->Scale(factor);
1963  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1964 
1965  return *this;
1966  }
1967 
1968 
1969 
1970  inline Vector &
1971  Vector::operator+=(const Vector &v)
1972  {
1973  AssertDimension(size(), v.size());
1974  Assert(vector->Map().SameAs(v.vector->Map()),
1976 
1977  const int ierr = vector->Update(1.0, *(v.vector), 1.0);
1978  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1979 
1980  return *this;
1981  }
1982 
1983 
1984 
1985  inline Vector &
1986  Vector::operator-=(const Vector &v)
1987  {
1988  AssertDimension(size(), v.size());
1989  Assert(vector->Map().SameAs(v.vector->Map()),
1991 
1992  const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
1993  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1994 
1995  return *this;
1996  }
1997 
1998 
1999 
2000  inline void
2001  Vector::add(const TrilinosScalar s)
2002  {
2003  // if we have ghost values, do not allow
2004  // writing to this vector at all.
2006  AssertIsFinite(s);
2007 
2008  size_type n_local = locally_owned_size();
2009  for (size_type i = 0; i < n_local; ++i)
2010  (*vector)[0][i] += s;
2011  }
2012 
2013 
2014 
2015  inline void
2016  Vector::add(const TrilinosScalar a, const Vector &v)
2017  {
2018  // if we have ghost values, do not allow
2019  // writing to this vector at all.
2022 
2023  AssertIsFinite(a);
2024 
2025  const int ierr = vector->Update(a, *(v.vector), 1.);
2026  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2027  }
2028 
2029 
2030 
2031  inline void
2032  Vector::add(const TrilinosScalar a,
2033  const Vector & v,
2034  const TrilinosScalar b,
2035  const Vector & w)
2036  {
2037  // if we have ghost values, do not allow
2038  // writing to this vector at all.
2041  AssertDimension(locally_owned_size(), w.locally_owned_size());
2042 
2043  AssertIsFinite(a);
2044  AssertIsFinite(b);
2045 
2046  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2047 
2048  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2049  }
2050 
2051 
2052 
2053  inline void
2054  Vector::sadd(const TrilinosScalar s, const Vector &v)
2055  {
2056  // if we have ghost values, do not allow
2057  // writing to this vector at all.
2059  AssertDimension(size(), v.size());
2060 
2061  AssertIsFinite(s);
2062 
2063  // We assume that the vectors have the same Map
2064  // if the local size is the same and if the vectors are not ghosted
2065  if (locally_owned_size() == v.locally_owned_size() &&
2066  !v.has_ghost_elements())
2067  {
2068  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2070  const int ierr = vector->Update(1., *(v.vector), s);
2071  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2072  }
2073  else
2074  {
2075  (*this) *= s;
2076  this->add(v, true);
2077  }
2078  }
2079 
2080 
2081 
2082  inline void
2083  Vector::sadd(const TrilinosScalar s,
2084  const TrilinosScalar a,
2085  const Vector & v)
2086  {
2087  // if we have ghost values, do not allow
2088  // writing to this vector at all.
2090  AssertDimension(size(), v.size());
2091  AssertIsFinite(s);
2092  AssertIsFinite(a);
2093 
2094  // We assume that the vectors have the same Map
2095  // if the local size is the same and if the vectors are not ghosted
2096  if (locally_owned_size() == v.locally_owned_size() &&
2097  !v.has_ghost_elements())
2098  {
2099  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2101  const int ierr = vector->Update(a, *(v.vector), s);
2102  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2103  }
2104  else
2105  {
2106  (*this) *= s;
2107  Vector tmp = v;
2108  tmp *= a;
2109  this->add(tmp, true);
2110  }
2111  }
2112 
2113 
2114 
2115  inline void
2116  Vector::scale(const Vector &factors)
2117  {
2118  // if we have ghost values, do not allow
2119  // writing to this vector at all.
2122 
2123  const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2124  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2125  }
2126 
2127 
2128 
2129  inline void
2130  Vector::equ(const TrilinosScalar a, const Vector &v)
2131  {
2132  // if we have ghost values, do not allow
2133  // writing to this vector at all.
2135  AssertIsFinite(a);
2136 
2137  // If we don't have the same map, copy.
2138  if (vector->Map().SameAs(v.vector->Map()) == false)
2139  {
2140  this->sadd(0., a, v);
2141  }
2142  else
2143  {
2144  // Otherwise, just update
2145  int ierr = vector->Update(a, *v.vector, 0.0);
2146  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2147 
2148  last_action = Zero;
2149  }
2150  }
2151 
2152 
2153 
2154  inline const Epetra_MultiVector &
2155  Vector::trilinos_vector() const
2156  {
2157  return static_cast<const Epetra_MultiVector &>(*vector);
2158  }
2159 
2160 
2161 
2162  inline Epetra_FEVector &
2164  {
2165  return *vector;
2166  }
2167 
2168 
2169 
2170  inline const Epetra_BlockMap &
2172  {
2173  return vector->Map();
2174  }
2175 
2176 
2177 
2178  inline const MPI_Comm &
2180  {
2181  static MPI_Comm comm;
2182 
2183  const Epetra_MpiComm *mpi_comm =
2184  dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2185  comm = mpi_comm->Comm();
2186 
2187  return comm;
2188  }
2189 
2190  template <typename number>
2191  Vector::Vector(const IndexSet & parallel_partitioner,
2192  const ::Vector<number> &v,
2193  const MPI_Comm & communicator)
2194  {
2195  *this =
2196  Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2197  owned_elements = parallel_partitioner;
2198  }
2199 
2200 
2201 
2202  inline Vector &
2204  {
2205  AssertIsFinite(s);
2206 
2207  int ierr = vector->PutScalar(s);
2208  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2209 
2210  if (nonlocal_vector.get() != nullptr)
2211  {
2212  ierr = nonlocal_vector->PutScalar(0.);
2213  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2214  }
2215 
2216  return *this;
2217  }
2218  } /* end of namespace MPI */
2219 
2220 # endif /* DOXYGEN */
2221 
2222 } /* end of namespace TrilinosWrappers */
2223 
2227 namespace internal
2228 {
2229  namespace LinearOperatorImplementation
2230  {
2231  template <typename>
2232  class ReinitHelper;
2233 
2238  template <>
2240  {
2241  public:
2242  template <typename Matrix>
2243  static void
2244  reinit_range_vector(const Matrix & matrix,
2246  bool omit_zeroing_entries)
2247  {
2248  v.reinit(matrix.locally_owned_range_indices(),
2249  matrix.get_mpi_communicator(),
2250  omit_zeroing_entries);
2251  }
2252 
2253  template <typename Matrix>
2254  static void
2257  bool omit_zeroing_entries)
2258  {
2259  v.reinit(matrix.locally_owned_domain_indices(),
2260  matrix.get_mpi_communicator(),
2261  omit_zeroing_entries);
2262  }
2263  };
2264 
2265  } // namespace LinearOperatorImplementation
2266 } /* namespace internal */
2267 
2268 
2269 
2273 template <>
2275 {};
2276 
2277 
2279 
2280 #endif
2281 
2282 #endif
size_type size() const
Definition: index_set.h:1648
size_type n_elements() const
Definition: index_set.h:1796
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:789
void compress(VectorOperation::values operation)
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)
const MPI_Comm & get_mpi_communicator() const
std::unique_ptr< Epetra_MultiVector > nonlocal_vector
Epetra_FEVector & trilinos_vector()
void add(const TrilinosScalar s)
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)
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
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)
void swap(Vector &u, Vector &v)
TrilinosScalar max() const
const_iterator end() const
Definition: vector.h:109
typename numbers::NumberTraits< Number >::real_type real_type
Definition: vector.h:146
bool has_ghost_elements() const
const value_type * const_iterator
Definition: vector.h:132
types::global_dof_index size_type
Definition: vector.h:135
size_type size() const
value_type * iterator
Definition: vector.h:131
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:174
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
#define DeclException0(Exception0)
Definition: exceptions.h:465
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:579
#define Assert(cond, exc)
Definition: exceptions.h:1586
#define AssertIsFinite(number)
Definition: exceptions.h:1854
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
static ::ExceptionBase & ExcDifferentParallelPartitioning()
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:510
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:1675
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)
const MPI_Comm & comm
double TrilinosScalar
Definition: types.h:175