Reference documentation for deal.II version Git 1dc1051882 2021-04-22 23:57:03 +0200
\(\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 - 2020 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 
35 # include <memory>
36 # include <utility>
37 # include <vector>
38 # ifdef DEAL_II_WITH_MPI // only if MPI is installed
39 # include <Epetra_MpiComm.h>
40 # include <mpi.h>
41 # else
42 # include <Epetra_SerialComm.h>
43 # endif
44 # include <Epetra_FEVector.h>
45 # include <Epetra_LocalMap.h>
46 # include <Epetra_Map.h>
47 
49 
50 // Forward declarations
51 # ifndef DOXYGEN
52 namespace LinearAlgebra
53 {
54  // Forward declaration
55  template <typename Number>
56  class ReadWriteVector;
57 } // namespace LinearAlgebra
58 # endif
59 
70 namespace TrilinosWrappers
71 {
72  class SparseMatrix;
73 
84  namespace internal
85  {
90 
100  class VectorReference
101  {
102  private:
107  VectorReference(MPI::Vector &vector, const size_type index);
108 
109  public:
113  VectorReference(const VectorReference &) = default;
114 
126  const VectorReference &
127  operator=(const VectorReference &r) const;
128 
132  VectorReference &
133  operator=(const VectorReference &r);
134 
138  const VectorReference &
139  operator=(const TrilinosScalar &s) const;
140 
144  const VectorReference &
145  operator+=(const TrilinosScalar &s) const;
146 
150  const VectorReference &
151  operator-=(const TrilinosScalar &s) const;
152 
156  const VectorReference &
157  operator*=(const TrilinosScalar &s) const;
158 
162  const VectorReference &
163  operator/=(const TrilinosScalar &s) const;
164 
169  operator TrilinosScalar() const;
170 
175  int,
176  << "An error with error number " << arg1
177  << " occurred while calling a Trilinos function");
178 
179  private:
183  MPI::Vector &vector;
184 
188  const size_type index;
189 
190  // Make the vector class a friend, so that it can create objects of the
191  // present type.
193  };
194  } // namespace internal
199 # ifndef DEAL_II_WITH_64BIT_INDICES
200  // define a helper function that queries the global ID of local ID of
201  // an Epetra_BlockMap object by calling either the 32- or 64-bit
202  // function necessary.
203  inline int
204  gid(const Epetra_BlockMap &map, int i)
205  {
206  return map.GID(i);
207  }
208 # else
209  // define a helper function that queries the global ID of local ID of
210  // an Epetra_BlockMap object by calling either the 32- or 64-bit
211  // function necessary.
212  inline long long int
213  gid(const Epetra_BlockMap &map, int i)
214  {
215  return map.GID64(i);
216  }
217 # endif
218 
224  namespace MPI
225  {
226  class BlockVector;
227 
403  class Vector : public Subscriptor
404  {
405  public:
414  using iterator = value_type *;
415  using const_iterator = const value_type *;
416  using reference = internal::VectorReference;
417  using const_reference = const internal::VectorReference;
418 
428  Vector();
429 
433  Vector(const Vector &v);
434 
453  explicit Vector(const IndexSet &parallel_partitioning,
454  const MPI_Comm &communicator = MPI_COMM_WORLD);
455 
467  Vector(const IndexSet &local,
468  const IndexSet &ghost,
469  const MPI_Comm &communicator = MPI_COMM_WORLD);
470 
485  Vector(const IndexSet &parallel_partitioning,
486  const Vector & v,
487  const MPI_Comm &communicator = MPI_COMM_WORLD);
488 
501  template <typename Number>
502  Vector(const IndexSet & parallel_partitioning,
503  const ::Vector<Number> &v,
504  const MPI_Comm & communicator = MPI_COMM_WORLD);
505 
510  Vector(Vector &&v) noexcept;
511 
515  ~Vector() override = default;
516 
521  void
522  clear();
523 
547  void
548  reinit(const Vector &v,
549  const bool omit_zeroing_entries = false,
550  const bool allow_different_maps = false);
551 
574  void
575  reinit(const IndexSet &parallel_partitioning,
576  const MPI_Comm &communicator = MPI_COMM_WORLD,
577  const bool omit_zeroing_entries = false);
578 
604  void
605  reinit(const IndexSet &locally_owned_entries,
606  const IndexSet &ghost_entries,
607  const MPI_Comm &communicator = MPI_COMM_WORLD,
608  const bool vector_writable = false);
609 
613  void
614  reinit(const BlockVector &v, const bool import_data = false);
615 
632  void
633  compress(::VectorOperation::values operation);
634 
647  Vector &
648  operator=(const TrilinosScalar s);
649 
655  Vector &
656  operator=(const Vector &v);
657 
662  Vector &
663  operator=(Vector &&v) noexcept;
664 
672  template <typename Number>
673  Vector &
674  operator=(const ::Vector<Number> &v);
675 
693  void
694  import_nonlocal_data_for_fe(
696  const Vector & vector);
697 
704  void
705  import(const LinearAlgebra::ReadWriteVector<double> &rwv,
706  const VectorOperation::values operation);
707 
708 
714  bool
715  operator==(const Vector &v) const;
716 
722  bool
723  operator!=(const Vector &v) const;
724 
728  size_type
729  size() const;
730 
745  size_type
746  local_size() const;
747 
752  size_type
753  locally_owned_size() const;
754 
776  std::pair<size_type, size_type>
777  local_range() const;
778 
786  bool
787  in_local_range(const size_type index) const;
788 
802  IndexSet
803  locally_owned_elements() const;
804 
812  bool
813  has_ghost_elements() const;
814 
821  void
822  update_ghost_values() const;
823 
828  TrilinosScalar operator*(const Vector &vec) const;
829 
833  real_type
834  norm_sqr() const;
835 
840  mean_value() const;
841 
846  min() const;
847 
852  max() const;
853 
857  real_type
858  l1_norm() const;
859 
864  real_type
865  l2_norm() const;
866 
871  real_type
872  lp_norm(const TrilinosScalar p) const;
873 
877  real_type
878  linfty_norm() const;
879 
900  add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
901 
907  bool
908  all_zero() const;
909 
915  bool
916  is_non_negative() const;
918 
919 
924 
932  reference
933  operator()(const size_type index);
934 
943  operator()(const size_type index) const;
944 
950  reference operator[](const size_type index);
951 
957  TrilinosScalar operator[](const size_type index) const;
958 
974  void
975  extract_subvector_to(const std::vector<size_type> &indices,
976  std::vector<TrilinosScalar> & values) const;
977 
1005  template <typename ForwardIterator, typename OutputIterator>
1006  void
1007  extract_subvector_to(ForwardIterator indices_begin,
1008  const ForwardIterator indices_end,
1009  OutputIterator values_begin) const;
1010 
1021  iterator
1022  begin();
1023 
1029  begin() const;
1030 
1035  iterator
1036  end();
1037 
1043  end() const;
1044 
1046 
1047 
1052 
1059  void
1060  set(const std::vector<size_type> & indices,
1061  const std::vector<TrilinosScalar> &values);
1062 
1067  void
1068  set(const std::vector<size_type> & indices,
1069  const ::Vector<TrilinosScalar> &values);
1070 
1076  void
1077  set(const size_type n_elements,
1078  const size_type * indices,
1079  const TrilinosScalar *values);
1080 
1085  void
1086  add(const std::vector<size_type> & indices,
1087  const std::vector<TrilinosScalar> &values);
1088 
1093  void
1094  add(const std::vector<size_type> & indices,
1095  const ::Vector<TrilinosScalar> &values);
1096 
1102  void
1103  add(const size_type n_elements,
1104  const size_type * indices,
1105  const TrilinosScalar *values);
1106 
1110  Vector &
1111  operator*=(const TrilinosScalar factor);
1112 
1116  Vector &
1117  operator/=(const TrilinosScalar factor);
1118 
1122  Vector &
1123  operator+=(const Vector &V);
1124 
1128  Vector &
1129  operator-=(const Vector &V);
1130 
1135  void
1136  add(const TrilinosScalar s);
1137 
1150  void
1151  add(const Vector &V, const bool allow_different_maps = false);
1152 
1156  void
1157  add(const TrilinosScalar a, const Vector &V);
1158 
1162  void
1163  add(const TrilinosScalar a,
1164  const Vector & V,
1165  const TrilinosScalar b,
1166  const Vector & W);
1167 
1172  void
1173  sadd(const TrilinosScalar s, const Vector &V);
1174 
1178  void
1179  sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1180 
1186  void
1187  scale(const Vector &scaling_factors);
1188 
1192  void
1193  equ(const TrilinosScalar a, const Vector &V);
1195 
1200 
1205  const Epetra_MultiVector &
1206  trilinos_vector() const;
1207 
1212  Epetra_FEVector &
1213  trilinos_vector();
1214 
1219  const Epetra_BlockMap &
1220  trilinos_partitioner() const;
1221 
1229  void
1230  print(std::ostream & out,
1231  const unsigned int precision = 3,
1232  const bool scientific = true,
1233  const bool across = true) const;
1234 
1248  void
1249  swap(Vector &v);
1250 
1254  std::size_t
1255  memory_consumption() const;
1256 
1261  const MPI_Comm &
1262  get_mpi_communicator() const;
1264 
1268  DeclException0(ExcDifferentParallelPartitioning);
1269 
1274  int,
1275  << "An error with error number " << arg1
1276  << " occurred while calling a Trilinos function");
1277 
1282  ExcAccessToNonLocalElement,
1283  size_type,
1284  size_type,
1285  size_type,
1286  size_type,
1287  << "You are trying to access element " << arg1
1288  << " of a distributed vector, but this element is not stored "
1289  << "on the current processor. Note: There are " << arg2
1290  << " elements stored "
1291  << "on the current processor from within the range [" << arg3 << ","
1292  << arg4 << "] but Trilinos vectors need not store contiguous "
1293  << "ranges on each processor, and not every element in "
1294  << "this range may in fact be stored locally."
1295  << "\n\n"
1296  << "A common source for this kind of problem is that you "
1297  << "are passing a 'fully distributed' vector into a function "
1298  << "that needs read access to vector elements that correspond "
1299  << "to degrees of freedom on ghost cells (or at least to "
1300  << "'locally active' degrees of freedom that are not also "
1301  << "'locally owned'). You need to pass a vector that has these "
1302  << "elements as ghost entries.");
1303 
1304  private:
1316  Epetra_CombineMode last_action;
1317 
1323 
1329 
1335  std::unique_ptr<Epetra_FEVector> vector;
1336 
1342  std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1343 
1348 
1349  // Make the reference class a friend.
1350  friend class internal::VectorReference;
1351  };
1352 
1353 
1354 
1355  // ------------------- inline and template functions --------------
1356 
1357 
1365  inline void
1367  {
1368  u.swap(v);
1369  }
1370  } // namespace MPI
1371 
1372 # ifndef DOXYGEN
1373 
1374  namespace internal
1375  {
1376  inline VectorReference::VectorReference(MPI::Vector & vector,
1377  const size_type index)
1378  : vector(vector)
1379  , index(index)
1380  {}
1381 
1382 
1383  inline const VectorReference &
1384  VectorReference::operator=(const VectorReference &r) const
1385  {
1386  // as explained in the class
1387  // documentation, this is not the copy
1388  // operator. so simply pass on to the
1389  // "correct" assignment operator
1390  *this = static_cast<TrilinosScalar>(r);
1391 
1392  return *this;
1393  }
1394 
1395 
1396 
1397  inline VectorReference &
1398  VectorReference::operator=(const VectorReference &r)
1399  {
1400  // as above
1401  *this = static_cast<TrilinosScalar>(r);
1402 
1403  return *this;
1404  }
1405 
1406 
1407  inline const VectorReference &
1408  VectorReference::operator=(const TrilinosScalar &value) const
1409  {
1410  vector.set(1, &index, &value);
1411  return *this;
1412  }
1413 
1414 
1415 
1416  inline const VectorReference &
1417  VectorReference::operator+=(const TrilinosScalar &value) const
1418  {
1419  vector.add(1, &index, &value);
1420  return *this;
1421  }
1422 
1423 
1424 
1425  inline const VectorReference &
1426  VectorReference::operator-=(const TrilinosScalar &value) const
1427  {
1428  TrilinosScalar new_value = -value;
1429  vector.add(1, &index, &new_value);
1430  return *this;
1431  }
1432 
1433 
1434 
1435  inline const VectorReference &
1436  VectorReference::operator*=(const TrilinosScalar &value) const
1437  {
1438  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1439  vector.set(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  } // namespace internal
1453 
1454  namespace MPI
1455  {
1456  inline bool
1457  Vector::in_local_range(const size_type index) const
1458  {
1459  std::pair<size_type, size_type> range = local_range();
1460 
1461  return ((index >= range.first) && (index < range.second));
1462  }
1463 
1464 
1465 
1466  inline IndexSet
1468  {
1469  Assert(owned_elements.size() == size(),
1470  ExcMessage(
1471  "The locally owned elements have not been properly initialized!"
1472  " This happens for example if this object has been initialized"
1473  " with exactly one overlapping IndexSet."));
1474  return owned_elements;
1475  }
1476 
1477 
1478 
1479  inline bool
1481  {
1482  return has_ghosts;
1483  }
1484 
1485 
1486 
1487  inline void
1489  {}
1490 
1491 
1492 
1493  inline internal::VectorReference
1494  Vector::operator()(const size_type index)
1495  {
1496  return internal::VectorReference(*this, index);
1497  }
1498 
1499 
1500 
1501  inline internal::VectorReference Vector::operator[](const size_type index)
1502  {
1503  return operator()(index);
1504  }
1505 
1506 
1507 
1508  inline TrilinosScalar Vector::operator[](const size_type index) const
1509  {
1510  return operator()(index);
1511  }
1512 
1513 
1514 
1515  inline void
1516  Vector::extract_subvector_to(const std::vector<size_type> &indices,
1517  std::vector<TrilinosScalar> & values) const
1518  {
1519  for (size_type i = 0; i < indices.size(); ++i)
1520  values[i] = operator()(indices[i]);
1521  }
1522 
1523 
1524 
1525  template <typename ForwardIterator, typename OutputIterator>
1526  inline void
1527  Vector::extract_subvector_to(ForwardIterator indices_begin,
1528  const ForwardIterator indices_end,
1529  OutputIterator values_begin) const
1530  {
1531  while (indices_begin != indices_end)
1532  {
1533  *values_begin = operator()(*indices_begin);
1534  indices_begin++;
1535  values_begin++;
1536  }
1537  }
1538 
1539 
1540 
1541  inline Vector::iterator
1542  Vector::begin()
1543  {
1544  return (*vector)[0];
1545  }
1546 
1547 
1548 
1549  inline Vector::iterator
1550  Vector::end()
1551  {
1552  return (*vector)[0] + local_size();
1553  }
1554 
1555 
1556 
1557  inline Vector::const_iterator
1558  Vector::begin() const
1559  {
1560  return (*vector)[0];
1561  }
1562 
1563 
1564 
1565  inline Vector::const_iterator
1566  Vector::end() const
1567  {
1568  return (*vector)[0] + local_size();
1569  }
1570 
1571 
1572 
1573  inline void
1574  Vector::set(const std::vector<size_type> & indices,
1575  const std::vector<TrilinosScalar> &values)
1576  {
1577  // if we have ghost values, do not allow
1578  // writing to this vector at all.
1579  Assert(!has_ghost_elements(), ExcGhostsPresent());
1580 
1581  Assert(indices.size() == values.size(),
1582  ExcDimensionMismatch(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  Assert(indices.size() == values.size(),
1598  ExcDimensionMismatch(indices.size(), values.size()));
1599 
1600  set(indices.size(), indices.data(), values.begin());
1601  }
1602 
1603 
1604 
1605  inline void
1606  Vector::set(const size_type n_elements,
1607  const size_type * indices,
1608  const TrilinosScalar *values)
1609  {
1610  // if we have ghost values, do not allow
1611  // writing to this vector at all.
1612  Assert(!has_ghost_elements(), ExcGhostsPresent());
1613 
1614  if (last_action == Add)
1615  {
1616  const int ierr = vector->GlobalAssemble(Add);
1617  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1618  }
1619 
1620  if (last_action != Insert)
1621  last_action = Insert;
1622 
1623  for (size_type i = 0; i < n_elements; ++i)
1624  {
1625  const TrilinosWrappers::types::int_type row = indices[i];
1626  const TrilinosWrappers::types::int_type local_row =
1627  vector->Map().LID(row);
1628  if (local_row != -1)
1629  (*vector)[0][local_row] = values[i];
1630  else
1631  {
1632  const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1633  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1634  compressed = false;
1635  }
1636  // in set operation, do not use the pre-allocated vector for nonlocal
1637  // entries even if it exists. This is to ensure that we really only
1638  // set the elements touched by the set() method and not all contained
1639  // in the nonlocal entries vector (there is no way to distinguish them
1640  // on the receiving processor)
1641  }
1642  }
1643 
1644 
1645 
1646  inline void
1647  Vector::add(const std::vector<size_type> & indices,
1648  const std::vector<TrilinosScalar> &values)
1649  {
1650  // if we have ghost values, do not allow
1651  // writing to this vector at all.
1652  Assert(!has_ghost_elements(), ExcGhostsPresent());
1653  Assert(indices.size() == values.size(),
1654  ExcDimensionMismatch(indices.size(), values.size()));
1655 
1656  add(indices.size(), indices.data(), values.data());
1657  }
1658 
1659 
1660 
1661  inline void
1662  Vector::add(const std::vector<size_type> & indices,
1663  const ::Vector<TrilinosScalar> &values)
1664  {
1665  // if we have ghost values, do not allow
1666  // writing to this vector at all.
1667  Assert(!has_ghost_elements(), ExcGhostsPresent());
1668  Assert(indices.size() == values.size(),
1669  ExcDimensionMismatch(indices.size(), values.size()));
1670 
1671  add(indices.size(), indices.data(), values.begin());
1672  }
1673 
1674 
1675 
1676  inline void
1677  Vector::add(const size_type n_elements,
1678  const size_type * indices,
1679  const TrilinosScalar *values)
1680  {
1681  // if we have ghost values, do not allow
1682  // writing to this vector at all.
1683  Assert(!has_ghost_elements(), ExcGhostsPresent());
1684 
1685  if (last_action != Add)
1686  {
1687  if (last_action == Insert)
1688  {
1689  const int ierr = vector->GlobalAssemble(Insert);
1690  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1691  }
1692  last_action = Add;
1693  }
1694 
1695  for (size_type i = 0; i < n_elements; ++i)
1696  {
1697  const size_type row = indices[i];
1698  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1699  static_cast<TrilinosWrappers::types::int_type>(row));
1700  if (local_row != -1)
1701  (*vector)[0][local_row] += values[i];
1702  else if (nonlocal_vector.get() == nullptr)
1703  {
1704  const int ierr = vector->SumIntoGlobalValues(
1705  1,
1706  reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1707  &row),
1708  &values[i]);
1709  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1710  compressed = false;
1711  }
1712  else
1713  {
1714  // use pre-allocated vector for non-local entries if it exists for
1715  // addition operation
1716  const TrilinosWrappers::types::int_type my_row =
1717  nonlocal_vector->Map().LID(
1718  static_cast<TrilinosWrappers::types::int_type>(row));
1719  Assert(my_row != -1,
1720  ExcMessage(
1721  "Attempted to write into off-processor vector entry "
1722  "that has not be specified as being writable upon "
1723  "initialization"));
1724  (*nonlocal_vector)[0][my_row] += values[i];
1725  compressed = false;
1726  }
1727  }
1728  }
1729 
1730 
1731 
1732  inline Vector::size_type
1733  Vector::size() const
1734  {
1735 # ifndef DEAL_II_WITH_64BIT_INDICES
1736  return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1737 # else
1738  return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1739 # endif
1740  }
1741 
1742 
1743 
1744  inline Vector::size_type
1745  Vector::local_size() const
1746  {
1747  return vector->Map().NumMyElements();
1748  }
1749 
1750 
1751 
1752  inline Vector::size_type
1754  {
1755  return owned_elements.n_elements();
1756  }
1757 
1758 
1759 
1760  inline std::pair<Vector::size_type, Vector::size_type>
1761  Vector::local_range() const
1762  {
1763 # ifndef DEAL_II_WITH_64BIT_INDICES
1764  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1766  vector->Map().MaxMyGID() + 1;
1767 # else
1768  const TrilinosWrappers::types::int_type begin =
1769  vector->Map().MinMyGID64();
1771  vector->Map().MaxMyGID64() + 1;
1772 # endif
1773 
1774  Assert(
1775  end - begin == vector->Map().NumMyElements(),
1776  ExcMessage(
1777  "This function only makes sense if the elements that this "
1778  "vector stores on the current processor form a contiguous range. "
1779  "This does not appear to be the case for the current vector."));
1780 
1781  return std::make_pair(begin, end);
1782  }
1783 
1784 
1785 
1786  inline TrilinosScalar Vector::operator*(const Vector &vec) const
1787  {
1788  Assert(vector->Map().SameAs(vec.vector->Map()),
1789  ExcDifferentParallelPartitioning());
1790  Assert(!has_ghost_elements(), ExcGhostsPresent());
1791 
1792  TrilinosScalar result;
1793 
1794  const int ierr = vector->Dot(*(vec.vector), &result);
1795  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1796 
1797  return result;
1798  }
1799 
1800 
1801 
1802  inline Vector::real_type
1803  Vector::norm_sqr() const
1804  {
1805  const TrilinosScalar d = l2_norm();
1806  return d * d;
1807  }
1808 
1809 
1810 
1811  inline TrilinosScalar
1812  Vector::mean_value() const
1813  {
1814  Assert(!has_ghost_elements(), ExcGhostsPresent());
1815 
1817  const int ierr = vector->MeanValue(&mean);
1818  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1819 
1820  return mean;
1821  }
1822 
1823 
1824 
1825  inline TrilinosScalar
1826  Vector::min() const
1827  {
1828  TrilinosScalar min_value;
1829  const int ierr = vector->MinValue(&min_value);
1830  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1831 
1832  return min_value;
1833  }
1834 
1835 
1836 
1837  inline TrilinosScalar
1838  Vector::max() const
1839  {
1840  TrilinosScalar max_value;
1841  const int ierr = vector->MaxValue(&max_value);
1842  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1843 
1844  return max_value;
1845  }
1846 
1847 
1848 
1849  inline Vector::real_type
1850  Vector::l1_norm() const
1851  {
1852  Assert(!has_ghost_elements(), ExcGhostsPresent());
1853 
1854  TrilinosScalar d;
1855  const int ierr = vector->Norm1(&d);
1856  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1857 
1858  return d;
1859  }
1860 
1861 
1862 
1863  inline Vector::real_type
1864  Vector::l2_norm() const
1865  {
1866  Assert(!has_ghost_elements(), ExcGhostsPresent());
1867 
1868  TrilinosScalar d;
1869  const int ierr = vector->Norm2(&d);
1870  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1871 
1872  return d;
1873  }
1874 
1875 
1876 
1877  inline Vector::real_type
1878  Vector::lp_norm(const TrilinosScalar p) const
1879  {
1880  Assert(!has_ghost_elements(), ExcGhostsPresent());
1881 
1882  TrilinosScalar norm = 0;
1883  TrilinosScalar sum = 0;
1884  const size_type n_local = local_size();
1885 
1886  // loop over all the elements because
1887  // Trilinos does not support lp norms
1888  for (size_type i = 0; i < n_local; ++i)
1889  sum += std::pow(std::fabs((*vector)[0][i]), p);
1890 
1891  norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1892 
1893  return norm;
1894  }
1895 
1896 
1897 
1898  inline Vector::real_type
1899  Vector::linfty_norm() const
1900  {
1901  // while we disallow the other
1902  // norm operations on ghosted
1903  // vectors, this particular norm
1904  // is safe to run even in the
1905  // presence of ghost elements
1906  TrilinosScalar d;
1907  const int ierr = vector->NormInf(&d);
1908  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1909 
1910  return d;
1911  }
1912 
1913 
1914 
1915  inline TrilinosScalar
1917  const Vector & V,
1918  const Vector & W)
1919  {
1920  this->add(a, V);
1921  return *this * W;
1922  }
1923 
1924 
1925 
1926  // inline also scalar products, vector
1927  // additions etc. since they are all
1928  // representable by a single Trilinos
1929  // call. This reduces the overhead of the
1930  // wrapper class.
1931  inline Vector &
1933  {
1934  AssertIsFinite(a);
1935 
1936  const int ierr = vector->Scale(a);
1937  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1938 
1939  return *this;
1940  }
1941 
1942 
1943 
1944  inline Vector &
1946  {
1947  AssertIsFinite(a);
1948 
1949  const TrilinosScalar factor = 1. / a;
1950 
1951  AssertIsFinite(factor);
1952 
1953  const int ierr = vector->Scale(factor);
1954  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1955 
1956  return *this;
1957  }
1958 
1959 
1960 
1961  inline Vector &
1962  Vector::operator+=(const Vector &v)
1963  {
1964  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
1965  Assert(vector->Map().SameAs(v.vector->Map()),
1966  ExcDifferentParallelPartitioning());
1967 
1968  const int ierr = vector->Update(1.0, *(v.vector), 1.0);
1969  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1970 
1971  return *this;
1972  }
1973 
1974 
1975 
1976  inline Vector &
1977  Vector::operator-=(const Vector &v)
1978  {
1979  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
1980  Assert(vector->Map().SameAs(v.vector->Map()),
1981  ExcDifferentParallelPartitioning());
1982 
1983  const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
1984  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1985 
1986  return *this;
1987  }
1988 
1989 
1990 
1991  inline void
1992  Vector::add(const TrilinosScalar s)
1993  {
1994  // if we have ghost values, do not allow
1995  // writing to this vector at all.
1996  Assert(!has_ghost_elements(), ExcGhostsPresent());
1997  AssertIsFinite(s);
1998 
1999  size_type n_local = local_size();
2000  for (size_type i = 0; i < n_local; ++i)
2001  (*vector)[0][i] += s;
2002  }
2003 
2004 
2005 
2006  inline void
2007  Vector::add(const TrilinosScalar a, const Vector &v)
2008  {
2009  // if we have ghost values, do not allow
2010  // writing to this vector at all.
2011  Assert(!has_ghost_elements(), ExcGhostsPresent());
2012  Assert(local_size() == v.local_size(),
2013  ExcDimensionMismatch(local_size(), v.local_size()));
2014 
2015  AssertIsFinite(a);
2016 
2017  const int ierr = vector->Update(a, *(v.vector), 1.);
2018  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2019  }
2020 
2021 
2022 
2023  inline void
2024  Vector::add(const TrilinosScalar a,
2025  const Vector & v,
2026  const TrilinosScalar b,
2027  const Vector & w)
2028  {
2029  // if we have ghost values, do not allow
2030  // writing to this vector at all.
2031  Assert(!has_ghost_elements(), ExcGhostsPresent());
2032  Assert(local_size() == v.local_size(),
2033  ExcDimensionMismatch(local_size(), v.local_size()));
2034  Assert(local_size() == w.local_size(),
2035  ExcDimensionMismatch(local_size(), w.local_size()));
2036 
2037  AssertIsFinite(a);
2038  AssertIsFinite(b);
2039 
2040  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2041 
2042  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2043  }
2044 
2045 
2046 
2047  inline void
2048  Vector::sadd(const TrilinosScalar s, const Vector &v)
2049  {
2050  // if we have ghost values, do not allow
2051  // writing to this vector at all.
2052  Assert(!has_ghost_elements(), ExcGhostsPresent());
2053  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2054 
2055  AssertIsFinite(s);
2056 
2057  // We assume that the vectors have the same Map
2058  // if the local size is the same and if the vectors are not ghosted
2059  if (local_size() == v.local_size() && !v.has_ghost_elements())
2060  {
2061  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2062  ExcDifferentParallelPartitioning());
2063  const int ierr = vector->Update(1., *(v.vector), s);
2064  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2065  }
2066  else
2067  {
2068  (*this) *= s;
2069  this->add(v, true);
2070  }
2071  }
2072 
2073 
2074 
2075  inline void
2076  Vector::sadd(const TrilinosScalar s,
2077  const TrilinosScalar a,
2078  const Vector & v)
2079  {
2080  // if we have ghost values, do not allow
2081  // writing to this vector at all.
2082  Assert(!has_ghost_elements(), ExcGhostsPresent());
2083  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2084  AssertIsFinite(s);
2085  AssertIsFinite(a);
2086 
2087  // We assume that the vectors have the same Map
2088  // if the local size is the same and if the vectors are not ghosted
2089  if (local_size() == v.local_size() && !v.has_ghost_elements())
2090  {
2091  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2092  ExcDifferentParallelPartitioning());
2093  const int ierr = vector->Update(a, *(v.vector), s);
2094  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2095  }
2096  else
2097  {
2098  (*this) *= s;
2099  Vector tmp = v;
2100  tmp *= a;
2101  this->add(tmp, true);
2102  }
2103  }
2104 
2105 
2106 
2107  inline void
2108  Vector::scale(const Vector &factors)
2109  {
2110  // if we have ghost values, do not allow
2111  // writing to this vector at all.
2112  Assert(!has_ghost_elements(), ExcGhostsPresent());
2113  Assert(local_size() == factors.local_size(),
2114  ExcDimensionMismatch(local_size(), factors.local_size()));
2115 
2116  const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2117  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2118  }
2119 
2120 
2121 
2122  inline void
2123  Vector::equ(const TrilinosScalar a, const Vector &v)
2124  {
2125  // if we have ghost values, do not allow
2126  // writing to this vector at all.
2127  Assert(!has_ghost_elements(), ExcGhostsPresent());
2128  AssertIsFinite(a);
2129 
2130  // If we don't have the same map, copy.
2131  if (vector->Map().SameAs(v.vector->Map()) == false)
2132  {
2133  this->sadd(0., a, v);
2134  }
2135  else
2136  {
2137  // Otherwise, just update
2138  int ierr = vector->Update(a, *v.vector, 0.0);
2139  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2140 
2141  last_action = Zero;
2142  }
2143  }
2144 
2145 
2146 
2147  inline const Epetra_MultiVector &
2148  Vector::trilinos_vector() const
2149  {
2150  return static_cast<const Epetra_MultiVector &>(*vector);
2151  }
2152 
2153 
2154 
2155  inline Epetra_FEVector &
2156  Vector::trilinos_vector()
2157  {
2158  return *vector;
2159  }
2160 
2161 
2162 
2163  inline const Epetra_BlockMap &
2164  Vector::trilinos_partitioner() const
2165  {
2166  return vector->Map();
2167  }
2168 
2169 
2170 
2171  inline const MPI_Comm &
2172  Vector::get_mpi_communicator() const
2173  {
2174  static MPI_Comm comm;
2175 
2176 # ifdef DEAL_II_WITH_MPI
2177 
2178  const Epetra_MpiComm *mpi_comm =
2179  dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2180  comm = mpi_comm->Comm();
2181 
2182 # else
2183 
2184  comm = MPI_COMM_SELF;
2185 
2186 # endif
2187 
2188  return comm;
2189  }
2190 
2191  template <typename number>
2192  Vector::Vector(const IndexSet & parallel_partitioner,
2193  const ::Vector<number> &v,
2194  const MPI_Comm & communicator)
2195  {
2196  *this =
2197  Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2198  owned_elements = parallel_partitioner;
2199  }
2200 
2201 
2202 
2203  inline Vector &
2205  {
2206  AssertIsFinite(s);
2207 
2208  int ierr = vector->PutScalar(s);
2209  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2210 
2211  if (nonlocal_vector.get() != nullptr)
2212  {
2213  ierr = nonlocal_vector->PutScalar(0.);
2214  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2215  }
2216 
2217  return *this;
2218  }
2219  } /* end of namespace MPI */
2220 
2221 # endif /* DOXYGEN */
2222 
2223 } /* end of namespace TrilinosWrappers */
2224 
2228 namespace internal
2229 {
2230  namespace LinearOperatorImplementation
2231  {
2232  template <typename>
2233  class ReinitHelper;
2234 
2239  template <>
2241  {
2242  public:
2243  template <typename Matrix>
2244  static void
2245  reinit_range_vector(const Matrix & matrix,
2247  bool omit_zeroing_entries)
2248  {
2249  v.reinit(matrix.locally_owned_range_indices(),
2250  matrix.get_mpi_communicator(),
2251  omit_zeroing_entries);
2252  }
2253 
2254  template <typename Matrix>
2255  static void
2258  bool omit_zeroing_entries)
2259  {
2260  v.reinit(matrix.locally_owned_domain_indices(),
2261  matrix.get_mpi_communicator(),
2262  omit_zeroing_entries);
2263  }
2264  };
2265 
2266  } // namespace LinearOperatorImplementation
2267 } /* namespace internal */
2268 
2269 
2270 
2274 template <>
2276 {};
2277 
2278 
2280 
2281 #endif // DEAL_II_WITH_TRILINOS
2282 
2283 /*---------------------------- trilinos_vector.h ---------------------------*/
2284 
2285 #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
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.
#define DEAL_II_DEPRECATED_EARLY
Definition: config.h:166
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:1025
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:1583
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:596
virtual ::IndexSet locally_owned_elements() const override
Definition: la_vector.h:462
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:515
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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:2817
#define DeclException0(Exception0)
Definition: exceptions.h:470
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
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:453
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)
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
size_type local_size() const
virtual VectorSpaceVector< Number >::real_type l2_norm() const override
*braid_SplitCommworld & comm
internal::VectorReference reference
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
VectorType::value_type * begin(VectorType &V)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:584
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2791
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
bool has_ghost_elements() const
#define AssertIsFinite(number)
Definition: exceptions.h:1729
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)