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