Reference documentation for deal.II version GIT 6dd81a2340 2022-11-29 05: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 
504  Vector(Vector &&v) noexcept;
505 
509  ~Vector() override = default;
510 
515  void
516  clear();
517 
541  void
542  reinit(const Vector &v,
543  const bool omit_zeroing_entries = false,
544  const bool allow_different_maps = false);
545 
568  void
569  reinit(const IndexSet &parallel_partitioning,
570  const MPI_Comm &communicator = MPI_COMM_WORLD,
571  const bool omit_zeroing_entries = false);
572 
598  void
599  reinit(const IndexSet &locally_owned_entries,
600  const IndexSet &ghost_entries,
601  const MPI_Comm &communicator = MPI_COMM_WORLD,
602  const bool vector_writable = false);
603 
608  void
609  reinit(
610  const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
611  const bool vector_writable = false);
612 
616  void
617  reinit(const BlockVector &v, const bool import_data = false);
618 
635  void
636  compress(::VectorOperation::values operation);
637 
650  Vector &
652 
658  Vector &
659  operator=(const Vector &v);
660 
665  Vector &
666  operator=(Vector &&v) noexcept;
667 
675  template <typename Number>
676  Vector &
677  operator=(const ::Vector<Number> &v);
678 
696  void
699  const Vector & vector);
700 
707  void
708  import(const LinearAlgebra::ReadWriteVector<double> &rwv,
709  const VectorOperation::values operation);
710 
711 
717  bool
718  operator==(const Vector &v) const;
719 
725  bool
726  operator!=(const Vector &v) const;
727 
731  size_type
732  size() const;
733 
748  size_type
749  local_size() const;
750 
755  size_type
757 
779  std::pair<size_type, size_type>
780  local_range() const;
781 
789  bool
790  in_local_range(const size_type index) const;
791 
805  IndexSet
807 
815  bool
817 
824  void
826 
832  operator*(const Vector &vec) const;
833 
837  real_type
838  norm_sqr() const;
839 
844  mean_value() const;
845 
850  min() const;
851 
856  max() const;
857 
861  real_type
862  l1_norm() const;
863 
868  real_type
869  l2_norm() const;
870 
875  real_type
876  lp_norm(const TrilinosScalar p) const;
877 
881  real_type
882  linfty_norm() const;
883 
904  add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
905 
911  bool
912  all_zero() const;
913 
919  bool
920  is_non_negative() const;
936  reference
937  operator()(const size_type index);
938 
947  operator()(const size_type index) const;
948 
954  reference
955  operator[](const size_type index);
956 
963  operator[](const size_type index) const;
964 
980  void
981  extract_subvector_to(const std::vector<size_type> &indices,
982  std::vector<TrilinosScalar> & values) const;
983 
1011  template <typename ForwardIterator, typename OutputIterator>
1012  void
1013  extract_subvector_to(ForwardIterator indices_begin,
1014  const ForwardIterator indices_end,
1015  OutputIterator values_begin) const;
1016 
1027  iterator
1029 
1035  begin() const;
1036 
1041  iterator
1042  end();
1043 
1049  end() const;
1050 
1065  void
1066  set(const std::vector<size_type> & indices,
1067  const std::vector<TrilinosScalar> &values);
1068 
1073  void
1074  set(const std::vector<size_type> & indices,
1075  const ::Vector<TrilinosScalar> &values);
1076 
1082  void
1083  set(const size_type n_elements,
1084  const size_type * indices,
1085  const TrilinosScalar *values);
1086 
1091  void
1092  add(const std::vector<size_type> & indices,
1093  const std::vector<TrilinosScalar> &values);
1094 
1099  void
1100  add(const std::vector<size_type> & indices,
1101  const ::Vector<TrilinosScalar> &values);
1102 
1108  void
1109  add(const size_type n_elements,
1110  const size_type * indices,
1111  const TrilinosScalar *values);
1112 
1116  Vector &
1118 
1122  Vector &
1124 
1128  Vector &
1130 
1134  Vector &
1136 
1141  void
1143 
1156  void
1157  add(const Vector &V, const bool allow_different_maps = false);
1158 
1162  void
1163  add(const TrilinosScalar a, const Vector &V);
1164 
1168  void
1170  const Vector & V,
1171  const TrilinosScalar b,
1172  const Vector & W);
1173 
1178  void
1179  sadd(const TrilinosScalar s, const Vector &V);
1180 
1184  void
1185  sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1186 
1192  void
1193  scale(const Vector &scaling_factors);
1194 
1198  void
1199  equ(const TrilinosScalar a, const Vector &V);
1211  const Epetra_MultiVector &
1213 
1218  Epetra_FEVector &
1220 
1225  const Epetra_BlockMap &
1227 
1235  void
1236  print(std::ostream & out,
1237  const unsigned int precision = 3,
1238  const bool scientific = true,
1239  const bool across = true) const;
1240 
1254  void
1255  swap(Vector &v);
1256 
1260  std::size_t
1261  memory_consumption() const;
1262 
1267  const MPI_Comm &
1275 
1280  int,
1281  << "An error with error number " << arg1
1282  << " occurred while calling a Trilinos function");
1283 
1289  size_type,
1290  size_type,
1291  size_type,
1292  size_type,
1293  << "You are trying to access element " << arg1
1294  << " of a distributed vector, but this element is not stored "
1295  << "on the current processor. Note: There are " << arg2
1296  << " elements stored "
1297  << "on the current processor from within the range [" << arg3 << ','
1298  << arg4 << "] but Trilinos vectors need not store contiguous "
1299  << "ranges on each processor, and not every element in "
1300  << "this range may in fact be stored locally."
1301  << "\n\n"
1302  << "A common source for this kind of problem is that you "
1303  << "are passing a 'fully distributed' vector into a function "
1304  << "that needs read access to vector elements that correspond "
1305  << "to degrees of freedom on ghost cells (or at least to "
1306  << "'locally active' degrees of freedom that are not also "
1307  << "'locally owned'). You need to pass a vector that has these "
1308  << "elements as ghost entries.");
1309 
1310  private:
1322  Epetra_CombineMode last_action;
1323 
1329 
1335 
1341  std::unique_ptr<Epetra_FEVector> vector;
1342 
1348  std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1349 
1354 
1355  // Make the reference class a friend.
1357  };
1358 
1359 
1360 
1361  // ------------------- inline and template functions --------------
1362 
1363 
1371  inline void
1373  {
1374  u.swap(v);
1375  }
1376  } // namespace MPI
1377 
1378 # ifndef DOXYGEN
1379 
1380  namespace internal
1381  {
1382  inline VectorReference::VectorReference(MPI::Vector & vector,
1383  const size_type index)
1384  : vector(vector)
1385  , index(index)
1386  {}
1387 
1388 
1389  inline const VectorReference &
1390  VectorReference::operator=(const VectorReference &r) const
1391  {
1392  // as explained in the class
1393  // documentation, this is not the copy
1394  // operator. so simply pass on to the
1395  // "correct" assignment operator
1396  *this = static_cast<TrilinosScalar>(r);
1397 
1398  return *this;
1399  }
1400 
1401 
1402 
1403  inline VectorReference &
1404  VectorReference::operator=(const VectorReference &r)
1405  {
1406  // as above
1407  *this = static_cast<TrilinosScalar>(r);
1408 
1409  return *this;
1410  }
1411 
1412 
1413  inline const VectorReference &
1414  VectorReference::operator=(const TrilinosScalar &value) const
1415  {
1416  vector.set(1, &index, &value);
1417  return *this;
1418  }
1419 
1420 
1421 
1422  inline const VectorReference &
1423  VectorReference::operator+=(const TrilinosScalar &value) const
1424  {
1425  vector.add(1, &index, &value);
1426  return *this;
1427  }
1428 
1429 
1430 
1431  inline const VectorReference &
1432  VectorReference::operator-=(const TrilinosScalar &value) const
1433  {
1434  TrilinosScalar new_value = -value;
1435  vector.add(1, &index, &new_value);
1436  return *this;
1437  }
1438 
1439 
1440 
1441  inline const VectorReference &
1442  VectorReference::operator*=(const TrilinosScalar &value) const
1443  {
1444  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1445  vector.set(1, &index, &new_value);
1446  return *this;
1447  }
1448 
1449 
1450 
1451  inline const VectorReference &
1452  VectorReference::operator/=(const TrilinosScalar &value) const
1453  {
1454  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1455  vector.set(1, &index, &new_value);
1456  return *this;
1457  }
1458  } // namespace internal
1459 
1460  namespace MPI
1461  {
1462  inline bool
1463  Vector::in_local_range(const size_type index) const
1464  {
1465  std::pair<size_type, size_type> range = local_range();
1466 
1467  return ((index >= range.first) && (index < range.second));
1468  }
1469 
1470 
1471 
1472  inline IndexSet
1474  {
1475  Assert(owned_elements.size() == size(),
1476  ExcMessage(
1477  "The locally owned elements have not been properly initialized!"
1478  " This happens for example if this object has been initialized"
1479  " with exactly one overlapping IndexSet."));
1480  return owned_elements;
1481  }
1482 
1483 
1484 
1485  inline bool
1487  {
1488  return has_ghosts;
1489  }
1490 
1491 
1492 
1493  inline void
1495  {}
1496 
1497 
1498 
1499  inline internal::VectorReference
1500  Vector::operator()(const size_type index)
1501  {
1502  return internal::VectorReference(*this, index);
1503  }
1504 
1505 
1506 
1507  inline internal::VectorReference
1508  Vector::operator[](const size_type index)
1509  {
1510  return operator()(index);
1511  }
1512 
1513 
1514 
1515  inline TrilinosScalar
1516  Vector::operator[](const size_type index) const
1517  {
1518  return operator()(index);
1519  }
1520 
1521 
1522 
1523  inline void
1524  Vector::extract_subvector_to(const std::vector<size_type> &indices,
1525  std::vector<TrilinosScalar> & values) const
1526  {
1527  for (size_type i = 0; i < indices.size(); ++i)
1528  values[i] = operator()(indices[i]);
1529  }
1530 
1531 
1532 
1533  template <typename ForwardIterator, typename OutputIterator>
1534  inline void
1535  Vector::extract_subvector_to(ForwardIterator indices_begin,
1536  const ForwardIterator indices_end,
1537  OutputIterator values_begin) const
1538  {
1539  while (indices_begin != indices_end)
1540  {
1541  *values_begin = operator()(*indices_begin);
1542  indices_begin++;
1543  values_begin++;
1544  }
1545  }
1546 
1547 
1548 
1549  inline Vector::iterator
1550  Vector::begin()
1551  {
1552  return (*vector)[0];
1553  }
1554 
1555 
1556 
1557  inline Vector::iterator
1558  Vector::end()
1559  {
1560  return (*vector)[0] + locally_owned_size();
1561  }
1562 
1563 
1564 
1565  inline Vector::const_iterator
1566  Vector::begin() const
1567  {
1568  return (*vector)[0];
1569  }
1570 
1571 
1572 
1573  inline Vector::const_iterator
1574  Vector::end() const
1575  {
1576  return (*vector)[0] + locally_owned_size();
1577  }
1578 
1579 
1580 
1581  inline void
1582  Vector::set(const std::vector<size_type> & indices,
1583  const std::vector<TrilinosScalar> &values)
1584  {
1585  // if we have ghost values, do not allow
1586  // writing to this vector at all.
1588 
1589  AssertDimension(indices.size(), values.size());
1590 
1591  set(indices.size(), indices.data(), values.data());
1592  }
1593 
1594 
1595 
1596  inline void
1597  Vector::set(const std::vector<size_type> & indices,
1598  const ::Vector<TrilinosScalar> &values)
1599  {
1600  // if we have ghost values, do not allow
1601  // writing to this vector at all.
1603 
1604  AssertDimension(indices.size(), values.size());
1605 
1606  set(indices.size(), indices.data(), values.begin());
1607  }
1608 
1609 
1610 
1611  inline void
1612  Vector::set(const size_type n_elements,
1613  const size_type * indices,
1614  const TrilinosScalar *values)
1615  {
1616  // if we have ghost values, do not allow
1617  // writing to this vector at all.
1619 
1620  if (last_action == Add)
1621  {
1622  const int ierr = vector->GlobalAssemble(Add);
1623  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1624  }
1625 
1626  if (last_action != Insert)
1627  last_action = Insert;
1628 
1629  for (size_type i = 0; i < n_elements; ++i)
1630  {
1631  const TrilinosWrappers::types::int_type row = indices[i];
1632  const TrilinosWrappers::types::int_type local_row =
1633  vector->Map().LID(row);
1634  if (local_row != -1)
1635  (*vector)[0][local_row] = values[i];
1636  else
1637  {
1638  const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1639  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1640  compressed = false;
1641  }
1642  // in set operation, do not use the pre-allocated vector for nonlocal
1643  // entries even if it exists. This is to ensure that we really only
1644  // set the elements touched by the set() method and not all contained
1645  // in the nonlocal entries vector (there is no way to distinguish them
1646  // on the receiving processor)
1647  }
1648  }
1649 
1650 
1651 
1652  inline void
1653  Vector::add(const std::vector<size_type> & indices,
1654  const std::vector<TrilinosScalar> &values)
1655  {
1656  // if we have ghost values, do not allow
1657  // writing to this vector at all.
1659  AssertDimension(indices.size(), values.size());
1660 
1661  add(indices.size(), indices.data(), values.data());
1662  }
1663 
1664 
1665 
1666  inline void
1667  Vector::add(const std::vector<size_type> & indices,
1668  const ::Vector<TrilinosScalar> &values)
1669  {
1670  // if we have ghost values, do not allow
1671  // writing to this vector at all.
1673  AssertDimension(indices.size(), values.size());
1674 
1675  add(indices.size(), indices.data(), values.begin());
1676  }
1677 
1678 
1679 
1680  inline void
1681  Vector::add(const size_type n_elements,
1682  const size_type * indices,
1683  const TrilinosScalar *values)
1684  {
1685  // if we have ghost values, do not allow
1686  // writing to this vector at all.
1688 
1689  if (last_action != Add)
1690  {
1691  if (last_action == Insert)
1692  {
1693  const int ierr = vector->GlobalAssemble(Insert);
1694  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1695  }
1696  last_action = Add;
1697  }
1698 
1699  for (size_type i = 0; i < n_elements; ++i)
1700  {
1701  const size_type row = indices[i];
1702  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1703  static_cast<TrilinosWrappers::types::int_type>(row));
1704  if (local_row != -1)
1705  (*vector)[0][local_row] += values[i];
1706  else if (nonlocal_vector.get() == nullptr)
1707  {
1708  const int ierr = vector->SumIntoGlobalValues(
1709  1,
1710  reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1711  &row),
1712  &values[i]);
1713  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1714  compressed = false;
1715  }
1716  else
1717  {
1718  // use pre-allocated vector for non-local entries if it exists for
1719  // addition operation
1720  const TrilinosWrappers::types::int_type my_row =
1721  nonlocal_vector->Map().LID(
1722  static_cast<TrilinosWrappers::types::int_type>(row));
1723  Assert(my_row != -1,
1724  ExcMessage(
1725  "Attempted to write into off-processor vector entry "
1726  "that has not be specified as being writable upon "
1727  "initialization"));
1728  (*nonlocal_vector)[0][my_row] += values[i];
1729  compressed = false;
1730  }
1731  }
1732  }
1733 
1734 
1735 
1736  inline Vector::size_type
1737  Vector::size() const
1738  {
1739 # ifndef DEAL_II_WITH_64BIT_INDICES
1740  return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1741 # else
1742  return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1743 # endif
1744  }
1745 
1746 
1747 
1748  inline Vector::size_type
1749  Vector::local_size() const
1750  {
1751  return vector->Map().NumMyElements();
1752  }
1753 
1754 
1755 
1756  inline Vector::size_type
1758  {
1759  return owned_elements.n_elements();
1760  }
1761 
1762 
1763 
1764  inline std::pair<Vector::size_type, Vector::size_type>
1765  Vector::local_range() const
1766  {
1767 # ifndef DEAL_II_WITH_64BIT_INDICES
1768  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1770  vector->Map().MaxMyGID() + 1;
1771 # else
1773  vector->Map().MinMyGID64();
1775  vector->Map().MaxMyGID64() + 1;
1776 # endif
1777 
1778  Assert(
1779  end - begin == vector->Map().NumMyElements(),
1780  ExcMessage(
1781  "This function only makes sense if the elements that this "
1782  "vector stores on the current processor form a contiguous range. "
1783  "This does not appear to be the case for the current vector."));
1784 
1785  return std::make_pair(begin, end);
1786  }
1787 
1788 
1789 
1790  inline TrilinosScalar
1791  Vector::operator*(const Vector &vec) const
1792  {
1793  Assert(vector->Map().SameAs(vec.vector->Map()),
1796 
1797  TrilinosScalar result;
1798 
1799  const int ierr = vector->Dot(*(vec.vector), &result);
1800  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1801 
1802  return result;
1803  }
1804 
1805 
1806 
1807  inline Vector::real_type
1808  Vector::norm_sqr() const
1809  {
1810  const TrilinosScalar d = l2_norm();
1811  return d * d;
1812  }
1813 
1814 
1815 
1816  inline TrilinosScalar
1817  Vector::mean_value() const
1818  {
1820 
1822  const int ierr = vector->MeanValue(&mean);
1823  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1824 
1825  return mean;
1826  }
1827 
1828 
1829 
1830  inline TrilinosScalar
1831  Vector::min() const
1832  {
1833  TrilinosScalar min_value;
1834  const int ierr = vector->MinValue(&min_value);
1835  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1836 
1837  return min_value;
1838  }
1839 
1840 
1841 
1842  inline TrilinosScalar
1843  Vector::max() const
1844  {
1845  TrilinosScalar max_value;
1846  const int ierr = vector->MaxValue(&max_value);
1847  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1848 
1849  return max_value;
1850  }
1851 
1852 
1853 
1854  inline Vector::real_type
1855  Vector::l1_norm() const
1856  {
1858 
1859  TrilinosScalar d;
1860  const int ierr = vector->Norm1(&d);
1861  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1862 
1863  return d;
1864  }
1865 
1866 
1867 
1868  inline Vector::real_type
1869  Vector::l2_norm() const
1870  {
1872 
1873  TrilinosScalar d;
1874  const int ierr = vector->Norm2(&d);
1875  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1876 
1877  return d;
1878  }
1879 
1880 
1881 
1882  inline Vector::real_type
1883  Vector::lp_norm(const TrilinosScalar p) const
1884  {
1886 
1887  TrilinosScalar norm = 0;
1888  TrilinosScalar sum = 0;
1889  const size_type n_local = locally_owned_size();
1890 
1891  // loop over all the elements because
1892  // Trilinos does not support lp norms
1893  for (size_type i = 0; i < n_local; ++i)
1894  sum += std::pow(std::fabs((*vector)[0][i]), p);
1895 
1896  norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1897 
1898  return norm;
1899  }
1900 
1901 
1902 
1903  inline Vector::real_type
1904  Vector::linfty_norm() const
1905  {
1906  // while we disallow the other
1907  // norm operations on ghosted
1908  // vectors, this particular norm
1909  // is safe to run even in the
1910  // presence of ghost elements
1911  TrilinosScalar d;
1912  const int ierr = vector->NormInf(&d);
1913  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1914 
1915  return d;
1916  }
1917 
1918 
1919 
1920  inline TrilinosScalar
1922  const Vector & V,
1923  const Vector & W)
1924  {
1925  this->add(a, V);
1926  return *this * W;
1927  }
1928 
1929 
1930 
1931  // inline also scalar products, vector
1932  // additions etc. since they are all
1933  // representable by a single Trilinos
1934  // call. This reduces the overhead of the
1935  // wrapper class.
1936  inline Vector &
1938  {
1939  AssertIsFinite(a);
1940 
1941  const int ierr = vector->Scale(a);
1942  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1943 
1944  return *this;
1945  }
1946 
1947 
1948 
1949  inline Vector &
1951  {
1952  AssertIsFinite(a);
1953 
1954  const TrilinosScalar factor = 1. / a;
1955 
1956  AssertIsFinite(factor);
1957 
1958  const int ierr = vector->Scale(factor);
1959  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1960 
1961  return *this;
1962  }
1963 
1964 
1965 
1966  inline Vector &
1967  Vector::operator+=(const Vector &v)
1968  {
1969  AssertDimension(size(), v.size());
1970  Assert(vector->Map().SameAs(v.vector->Map()),
1972 
1973  const int ierr = vector->Update(1.0, *(v.vector), 1.0);
1974  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1975 
1976  return *this;
1977  }
1978 
1979 
1980 
1981  inline Vector &
1982  Vector::operator-=(const Vector &v)
1983  {
1984  AssertDimension(size(), v.size());
1985  Assert(vector->Map().SameAs(v.vector->Map()),
1987 
1988  const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
1989  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1990 
1991  return *this;
1992  }
1993 
1994 
1995 
1996  inline void
1997  Vector::add(const TrilinosScalar s)
1998  {
1999  // if we have ghost values, do not allow
2000  // writing to this vector at all.
2002  AssertIsFinite(s);
2003 
2004  size_type n_local = locally_owned_size();
2005  for (size_type i = 0; i < n_local; ++i)
2006  (*vector)[0][i] += s;
2007  }
2008 
2009 
2010 
2011  inline void
2012  Vector::add(const TrilinosScalar a, const Vector &v)
2013  {
2014  // if we have ghost values, do not allow
2015  // writing to this vector at all.
2018 
2019  AssertIsFinite(a);
2020 
2021  const int ierr = vector->Update(a, *(v.vector), 1.);
2022  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2023  }
2024 
2025 
2026 
2027  inline void
2028  Vector::add(const TrilinosScalar a,
2029  const Vector & v,
2030  const TrilinosScalar b,
2031  const Vector & w)
2032  {
2033  // if we have ghost values, do not allow
2034  // writing to this vector at all.
2037  AssertDimension(locally_owned_size(), w.locally_owned_size());
2038 
2039  AssertIsFinite(a);
2040  AssertIsFinite(b);
2041 
2042  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2043 
2044  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2045  }
2046 
2047 
2048 
2049  inline void
2050  Vector::sadd(const TrilinosScalar s, const Vector &v)
2051  {
2052  // if we have ghost values, do not allow
2053  // writing to this vector at all.
2055  AssertDimension(size(), v.size());
2056 
2057  AssertIsFinite(s);
2058 
2059  // We assume that the vectors have the same Map
2060  // if the local size is the same and if the vectors are not ghosted
2061  if (locally_owned_size() == v.locally_owned_size() &&
2062  !v.has_ghost_elements())
2063  {
2064  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2066  const int ierr = vector->Update(1., *(v.vector), s);
2067  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2068  }
2069  else
2070  {
2071  (*this) *= s;
2072  this->add(v, true);
2073  }
2074  }
2075 
2076 
2077 
2078  inline void
2079  Vector::sadd(const TrilinosScalar s,
2080  const TrilinosScalar a,
2081  const Vector & v)
2082  {
2083  // if we have ghost values, do not allow
2084  // writing to this vector at all.
2086  AssertDimension(size(), v.size());
2087  AssertIsFinite(s);
2088  AssertIsFinite(a);
2089 
2090  // We assume that the vectors have the same Map
2091  // if the local size is the same and if the vectors are not ghosted
2092  if (locally_owned_size() == v.locally_owned_size() &&
2093  !v.has_ghost_elements())
2094  {
2095  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2097  const int ierr = vector->Update(a, *(v.vector), s);
2098  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2099  }
2100  else
2101  {
2102  (*this) *= s;
2103  Vector tmp = v;
2104  tmp *= a;
2105  this->add(tmp, true);
2106  }
2107  }
2108 
2109 
2110 
2111  inline void
2112  Vector::scale(const Vector &factors)
2113  {
2114  // if we have ghost values, do not allow
2115  // writing to this vector at all.
2118 
2119  const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2120  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2121  }
2122 
2123 
2124 
2125  inline void
2126  Vector::equ(const TrilinosScalar a, const Vector &v)
2127  {
2128  // if we have ghost values, do not allow
2129  // writing to this vector at all.
2131  AssertIsFinite(a);
2132 
2133  // If we don't have the same map, copy.
2134  if (vector->Map().SameAs(v.vector->Map()) == false)
2135  {
2136  this->sadd(0., a, v);
2137  }
2138  else
2139  {
2140  // Otherwise, just update
2141  int ierr = vector->Update(a, *v.vector, 0.0);
2142  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2143 
2144  last_action = Zero;
2145  }
2146  }
2147 
2148 
2149 
2150  inline const Epetra_MultiVector &
2151  Vector::trilinos_vector() const
2152  {
2153  return static_cast<const Epetra_MultiVector &>(*vector);
2154  }
2155 
2156 
2157 
2158  inline Epetra_FEVector &
2160  {
2161  return *vector;
2162  }
2163 
2164 
2165 
2166  inline const Epetra_BlockMap &
2168  {
2169  return vector->Map();
2170  }
2171 
2172 
2173 
2174  inline const MPI_Comm &
2176  {
2177  static MPI_Comm comm;
2178 
2179  const Epetra_MpiComm *mpi_comm =
2180  dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2181  comm = mpi_comm->Comm();
2182 
2183  return comm;
2184  }
2185 
2186  template <typename number>
2187  Vector::Vector(const IndexSet & parallel_partitioner,
2188  const ::Vector<number> &v,
2189  const MPI_Comm & communicator)
2190  {
2191  *this =
2192  Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2193  owned_elements = parallel_partitioner;
2194  }
2195 
2196 
2197 
2198  inline Vector &
2200  {
2201  AssertIsFinite(s);
2202 
2203  int ierr = vector->PutScalar(s);
2204  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2205 
2206  if (nonlocal_vector.get() != nullptr)
2207  {
2208  ierr = nonlocal_vector->PutScalar(0.);
2209  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2210  }
2211 
2212  return *this;
2213  }
2214  } /* end of namespace MPI */
2215 
2216 # endif /* DOXYGEN */
2217 
2218 } /* end of namespace TrilinosWrappers */
2219 
2223 namespace internal
2224 {
2225  namespace LinearOperatorImplementation
2226  {
2227  template <typename>
2228  class ReinitHelper;
2229 
2234  template <>
2236  {
2237  public:
2238  template <typename Matrix>
2239  static void
2240  reinit_range_vector(const Matrix & matrix,
2242  bool omit_zeroing_entries)
2243  {
2244  v.reinit(matrix.locally_owned_range_indices(),
2245  matrix.get_mpi_communicator(),
2246  omit_zeroing_entries);
2247  }
2248 
2249  template <typename Matrix>
2250  static void
2253  bool omit_zeroing_entries)
2254  {
2255  v.reinit(matrix.locally_owned_domain_indices(),
2256  matrix.get_mpi_communicator(),
2257  omit_zeroing_entries);
2258  }
2259  };
2260 
2261  } // namespace LinearOperatorImplementation
2262 } /* namespace internal */
2263 
2264 
2265 
2269 template <>
2271 {};
2272 
2273 
2275 
2276 #endif
2277 
2278 #endif
size_type size() const
Definition: index_set.h:1626
size_type n_elements() const
Definition: index_set.h:1774
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:777
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)
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)
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:164
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
#define DeclException0(Exception0)
Definition: exceptions.h:464
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:578
#define Assert(cond, exc)
Definition: exceptions.h:1501
#define AssertIsFinite(number)
Definition: exceptions.h:1786
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
static ::ExceptionBase & ExcDifferentParallelPartitioning()
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
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:1611
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:81
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const MPI_Comm & comm
double TrilinosScalar
Definition: types.h:168