Reference documentation for deal.II version 9.2.0
\(\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:
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 
220  namespace MPI
221  {
222  class BlockVector;
223 
400  class Vector : public Subscriptor
401  {
402  public:
411  using iterator = value_type *;
412  using const_iterator = const value_type *;
413  using reference = internal::VectorReference;
414  using const_reference = const internal::VectorReference;
415 
425  Vector();
426 
430  Vector(const Vector &v);
431 
450  explicit Vector(const IndexSet &parallel_partitioning,
451  const MPI_Comm &communicator = MPI_COMM_WORLD);
452 
464  Vector(const IndexSet &local,
465  const IndexSet &ghost,
466  const MPI_Comm &communicator = MPI_COMM_WORLD);
467 
482  Vector(const IndexSet &parallel_partitioning,
483  const Vector & v,
484  const MPI_Comm &communicator = MPI_COMM_WORLD);
485 
498  template <typename Number>
499  Vector(const IndexSet & parallel_partitioning,
500  const ::Vector<Number> &v,
501  const MPI_Comm & communicator = MPI_COMM_WORLD);
502 
507  Vector(Vector &&v) noexcept;
508 
512  ~Vector() override = default;
513 
518  void
519  clear();
520 
544  void
545  reinit(const Vector &v,
546  const bool omit_zeroing_entries = false,
547  const bool allow_different_maps = false);
548 
571  void
572  reinit(const IndexSet &parallel_partitioning,
573  const MPI_Comm &communicator = MPI_COMM_WORLD,
574  const bool omit_zeroing_entries = false);
575 
601  void
602  reinit(const IndexSet &locally_owned_entries,
603  const IndexSet &ghost_entries,
604  const MPI_Comm &communicator = MPI_COMM_WORLD,
605  const bool vector_writable = false);
606 
610  void
611  reinit(const BlockVector &v, const bool import_data = false);
612 
629  void
630  compress(::VectorOperation::values operation);
631 
644  Vector &
645  operator=(const TrilinosScalar s);
646 
652  Vector &
653  operator=(const Vector &v);
654 
659  Vector &
660  operator=(Vector &&v) noexcept;
661 
669  template <typename Number>
670  Vector &
671  operator=(const ::Vector<Number> &v);
672 
690  void
691  import_nonlocal_data_for_fe(
693  const Vector & vector);
694 
701  void
702  import(const LinearAlgebra::ReadWriteVector<double> &rwv,
703  const VectorOperation::values operation);
704 
705 
711  bool
712  operator==(const Vector &v) const;
713 
719  bool
720  operator!=(const Vector &v) const;
721 
725  size_type
726  size() const;
727 
739  size_type
740  local_size() const;
741 
763  std::pair<size_type, size_type>
764  local_range() const;
765 
773  bool
774  in_local_range(const size_type index) const;
775 
789  IndexSet
790  locally_owned_elements() const;
791 
799  bool
800  has_ghost_elements() const;
801 
808  void
809  update_ghost_values() const;
810 
815  TrilinosScalar operator*(const Vector &vec) const;
816 
820  real_type
821  norm_sqr() const;
822 
827  mean_value() const;
828 
833  min() const;
834 
839  max() const;
840 
844  real_type
845  l1_norm() const;
846 
851  real_type
852  l2_norm() const;
853 
858  real_type
859  lp_norm(const TrilinosScalar p) const;
860 
864  real_type
865  linfty_norm() const;
866 
887  add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
888 
894  bool
895  all_zero() const;
896 
902  bool
903  is_non_negative() const;
905 
906 
911 
919  reference
920  operator()(const size_type index);
921 
930  operator()(const size_type index) const;
931 
937  reference operator[](const size_type index);
938 
944  TrilinosScalar operator[](const size_type index) const;
945 
961  void
962  extract_subvector_to(const std::vector<size_type> &indices,
963  std::vector<TrilinosScalar> & values) const;
964 
992  template <typename ForwardIterator, typename OutputIterator>
993  void
994  extract_subvector_to(ForwardIterator indices_begin,
995  const ForwardIterator indices_end,
996  OutputIterator values_begin) const;
997 
1008  iterator
1009  begin();
1010 
1016  begin() const;
1017 
1022  iterator
1023  end();
1024 
1030  end() const;
1031 
1033 
1034 
1039 
1046  void
1047  set(const std::vector<size_type> & indices,
1048  const std::vector<TrilinosScalar> &values);
1049 
1054  void
1055  set(const std::vector<size_type> & indices,
1056  const ::Vector<TrilinosScalar> &values);
1057 
1063  void
1064  set(const size_type n_elements,
1065  const size_type * indices,
1066  const TrilinosScalar *values);
1067 
1072  void
1073  add(const std::vector<size_type> & indices,
1074  const std::vector<TrilinosScalar> &values);
1075 
1080  void
1081  add(const std::vector<size_type> & indices,
1082  const ::Vector<TrilinosScalar> &values);
1083 
1089  void
1090  add(const size_type n_elements,
1091  const size_type * indices,
1092  const TrilinosScalar *values);
1093 
1097  Vector &
1098  operator*=(const TrilinosScalar factor);
1099 
1103  Vector &
1104  operator/=(const TrilinosScalar factor);
1105 
1109  Vector &
1110  operator+=(const Vector &V);
1111 
1115  Vector &
1116  operator-=(const Vector &V);
1117 
1122  void
1123  add(const TrilinosScalar s);
1124 
1137  void
1138  add(const Vector &V, const bool allow_different_maps = false);
1139 
1143  void
1144  add(const TrilinosScalar a, const Vector &V);
1145 
1149  void
1150  add(const TrilinosScalar a,
1151  const Vector & V,
1152  const TrilinosScalar b,
1153  const Vector & W);
1154 
1159  void
1160  sadd(const TrilinosScalar s, const Vector &V);
1161 
1165  void
1166  sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1167 
1173  void
1174  scale(const Vector &scaling_factors);
1175 
1179  void
1180  equ(const TrilinosScalar a, const Vector &V);
1182 
1187 
1192  const Epetra_MultiVector &
1193  trilinos_vector() const;
1194 
1199  Epetra_FEVector &
1200  trilinos_vector();
1201 
1206  const Epetra_BlockMap &
1207  trilinos_partitioner() const;
1208 
1216  void
1217  print(std::ostream & out,
1218  const unsigned int precision = 3,
1219  const bool scientific = true,
1220  const bool across = true) const;
1221 
1235  void
1236  swap(Vector &v);
1237 
1241  std::size_t
1242  memory_consumption() const;
1243 
1248  const MPI_Comm &
1249  get_mpi_communicator() const;
1251 
1255  DeclException0(ExcDifferentParallelPartitioning);
1256 
1261  int,
1262  << "An error with error number " << arg1
1263  << " occurred while calling a Trilinos function");
1264 
1269  ExcAccessToNonLocalElement,
1270  size_type,
1271  size_type,
1272  size_type,
1273  size_type,
1274  << "You tried to access element " << arg1
1275  << " of a distributed vector, but this element is not stored "
1276  << "on the current processor. Note: There are " << arg2
1277  << " elements stored "
1278  << "on the current processor from within the range " << arg3
1279  << " through " << arg4
1280  << " but Trilinos vectors need not store contiguous "
1281  << "ranges on each processor, and not every element in "
1282  << "this range may in fact be stored locally.");
1283 
1284  private:
1296  Epetra_CombineMode last_action;
1297 
1303 
1309 
1315  std::unique_ptr<Epetra_FEVector> vector;
1316 
1322  std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1323 
1328 
1329  // Make the reference class a friend.
1330  friend class internal::VectorReference;
1331  };
1332 
1333 
1334 
1335  // ------------------- inline and template functions --------------
1336 
1337 
1346  inline void
1348  {
1349  u.swap(v);
1350  }
1351  } // namespace MPI
1352 
1353 # ifndef DOXYGEN
1354 
1355  namespace internal
1356  {
1357  inline VectorReference::VectorReference(MPI::Vector & vector,
1358  const size_type index)
1359  : vector(vector)
1360  , index(index)
1361  {}
1362 
1363 
1364  inline const VectorReference &
1365  VectorReference::operator=(const VectorReference &r) const
1366  {
1367  // as explained in the class
1368  // documentation, this is not the copy
1369  // operator. so simply pass on to the
1370  // "correct" assignment operator
1371  *this = static_cast<TrilinosScalar>(r);
1372 
1373  return *this;
1374  }
1375 
1376 
1377 
1378  inline VectorReference &
1379  VectorReference::operator=(const VectorReference &r)
1380  {
1381  // as above
1382  *this = static_cast<TrilinosScalar>(r);
1383 
1384  return *this;
1385  }
1386 
1387 
1388  inline const VectorReference &
1389  VectorReference::operator=(const TrilinosScalar &value) const
1390  {
1391  vector.set(1, &index, &value);
1392  return *this;
1393  }
1394 
1395 
1396 
1397  inline const VectorReference &
1398  VectorReference::operator+=(const TrilinosScalar &value) const
1399  {
1400  vector.add(1, &index, &value);
1401  return *this;
1402  }
1403 
1404 
1405 
1406  inline const VectorReference &
1407  VectorReference::operator-=(const TrilinosScalar &value) const
1408  {
1409  TrilinosScalar new_value = -value;
1410  vector.add(1, &index, &new_value);
1411  return *this;
1412  }
1413 
1414 
1415 
1416  inline const VectorReference &
1417  VectorReference::operator*=(const TrilinosScalar &value) const
1418  {
1419  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1420  vector.set(1, &index, &new_value);
1421  return *this;
1422  }
1423 
1424 
1425 
1426  inline const VectorReference &
1427  VectorReference::operator/=(const TrilinosScalar &value) const
1428  {
1429  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1430  vector.set(1, &index, &new_value);
1431  return *this;
1432  }
1433  } // namespace internal
1434 
1435  namespace MPI
1436  {
1437  inline bool
1438  Vector::in_local_range(const size_type index) const
1439  {
1440  std::pair<size_type, size_type> range = local_range();
1441 
1442  return ((index >= range.first) && (index < range.second));
1443  }
1444 
1445 
1446 
1447  inline IndexSet
1449  {
1450  Assert(owned_elements.size() == size(),
1451  ExcMessage(
1452  "The locally owned elements have not been properly initialized!"
1453  " This happens for example if this object has been initialized"
1454  " with exactly one overlapping IndexSet."));
1455  return owned_elements;
1456  }
1457 
1458 
1459 
1460  inline bool
1462  {
1463  return has_ghosts;
1464  }
1465 
1466 
1467 
1468  inline void
1470  {}
1471 
1472 
1473 
1474  inline internal::VectorReference
1475  Vector::operator()(const size_type index)
1476  {
1477  return internal::VectorReference(*this, index);
1478  }
1479 
1480 
1481 
1482  inline internal::VectorReference Vector::operator[](const size_type index)
1483  {
1484  return operator()(index);
1485  }
1486 
1487 
1488 
1489  inline TrilinosScalar Vector::operator[](const size_type index) const
1490  {
1491  return operator()(index);
1492  }
1493 
1494 
1495 
1496  inline void
1497  Vector::extract_subvector_to(const std::vector<size_type> &indices,
1498  std::vector<TrilinosScalar> & values) const
1499  {
1500  for (size_type i = 0; i < indices.size(); ++i)
1501  values[i] = operator()(indices[i]);
1502  }
1503 
1504 
1505 
1506  template <typename ForwardIterator, typename OutputIterator>
1507  inline void
1508  Vector::extract_subvector_to(ForwardIterator indices_begin,
1509  const ForwardIterator indices_end,
1510  OutputIterator values_begin) const
1511  {
1512  while (indices_begin != indices_end)
1513  {
1514  *values_begin = operator()(*indices_begin);
1515  indices_begin++;
1516  values_begin++;
1517  }
1518  }
1519 
1520 
1521 
1522  inline Vector::iterator
1523  Vector::begin()
1524  {
1525  return (*vector)[0];
1526  }
1527 
1528 
1529 
1530  inline Vector::iterator
1531  Vector::end()
1532  {
1533  return (*vector)[0] + local_size();
1534  }
1535 
1536 
1537 
1538  inline Vector::const_iterator
1539  Vector::begin() const
1540  {
1541  return (*vector)[0];
1542  }
1543 
1544 
1545 
1546  inline Vector::const_iterator
1547  Vector::end() const
1548  {
1549  return (*vector)[0] + local_size();
1550  }
1551 
1552 
1553 
1554  inline void
1555  Vector::set(const std::vector<size_type> & indices,
1556  const std::vector<TrilinosScalar> &values)
1557  {
1558  // if we have ghost values, do not allow
1559  // writing to this vector at all.
1560  Assert(!has_ghost_elements(), ExcGhostsPresent());
1561 
1562  Assert(indices.size() == values.size(),
1563  ExcDimensionMismatch(indices.size(), values.size()));
1564 
1565  set(indices.size(), indices.data(), values.data());
1566  }
1567 
1568 
1569 
1570  inline void
1571  Vector::set(const std::vector<size_type> & indices,
1572  const ::Vector<TrilinosScalar> &values)
1573  {
1574  // if we have ghost values, do not allow
1575  // writing to this vector at all.
1576  Assert(!has_ghost_elements(), ExcGhostsPresent());
1577 
1578  Assert(indices.size() == values.size(),
1579  ExcDimensionMismatch(indices.size(), values.size()));
1580 
1581  set(indices.size(), indices.data(), values.begin());
1582  }
1583 
1584 
1585 
1586  inline void
1587  Vector::set(const size_type n_elements,
1588  const size_type * indices,
1589  const TrilinosScalar *values)
1590  {
1591  // if we have ghost values, do not allow
1592  // writing to this vector at all.
1593  Assert(!has_ghost_elements(), ExcGhostsPresent());
1594 
1595  if (last_action == Add)
1596  {
1597  const int ierr = vector->GlobalAssemble(Add);
1598  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1599  }
1600 
1601  if (last_action != Insert)
1602  last_action = Insert;
1603 
1604  for (size_type i = 0; i < n_elements; ++i)
1605  {
1606  const TrilinosWrappers::types::int_type row = indices[i];
1607  const TrilinosWrappers::types::int_type local_row =
1608  vector->Map().LID(row);
1609  if (local_row != -1)
1610  (*vector)[0][local_row] = values[i];
1611  else
1612  {
1613  const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1614  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1615  compressed = false;
1616  }
1617  // in set operation, do not use the pre-allocated vector for nonlocal
1618  // entries even if it exists. This is to ensure that we really only
1619  // set the elements touched by the set() method and not all contained
1620  // in the nonlocal entries vector (there is no way to distinguish them
1621  // on the receiving processor)
1622  }
1623  }
1624 
1625 
1626 
1627  inline void
1628  Vector::add(const std::vector<size_type> & indices,
1629  const std::vector<TrilinosScalar> &values)
1630  {
1631  // if we have ghost values, do not allow
1632  // writing to this vector at all.
1633  Assert(!has_ghost_elements(), ExcGhostsPresent());
1634  Assert(indices.size() == values.size(),
1635  ExcDimensionMismatch(indices.size(), values.size()));
1636 
1637  add(indices.size(), indices.data(), values.data());
1638  }
1639 
1640 
1641 
1642  inline void
1643  Vector::add(const std::vector<size_type> & indices,
1644  const ::Vector<TrilinosScalar> &values)
1645  {
1646  // if we have ghost values, do not allow
1647  // writing to this vector at all.
1648  Assert(!has_ghost_elements(), ExcGhostsPresent());
1649  Assert(indices.size() == values.size(),
1650  ExcDimensionMismatch(indices.size(), values.size()));
1651 
1652  add(indices.size(), indices.data(), values.begin());
1653  }
1654 
1655 
1656 
1657  inline void
1658  Vector::add(const size_type n_elements,
1659  const size_type * indices,
1660  const TrilinosScalar *values)
1661  {
1662  // if we have ghost values, do not allow
1663  // writing to this vector at all.
1664  Assert(!has_ghost_elements(), ExcGhostsPresent());
1665 
1666  if (last_action != Add)
1667  {
1668  if (last_action == Insert)
1669  {
1670  const int ierr = vector->GlobalAssemble(Insert);
1671  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1672  }
1673  last_action = Add;
1674  }
1675 
1676  for (size_type i = 0; i < n_elements; ++i)
1677  {
1678  const size_type row = indices[i];
1679  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1680  static_cast<TrilinosWrappers::types::int_type>(row));
1681  if (local_row != -1)
1682  (*vector)[0][local_row] += values[i];
1683  else if (nonlocal_vector.get() == nullptr)
1684  {
1685  const int ierr = vector->SumIntoGlobalValues(
1686  1,
1687  reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1688  &row),
1689  &values[i]);
1690  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1691  compressed = false;
1692  }
1693  else
1694  {
1695  // use pre-allocated vector for non-local entries if it exists for
1696  // addition operation
1697  const TrilinosWrappers::types::int_type my_row =
1698  nonlocal_vector->Map().LID(
1699  static_cast<TrilinosWrappers::types::int_type>(row));
1700  Assert(my_row != -1,
1701  ExcMessage(
1702  "Attempted to write into off-processor vector entry "
1703  "that has not be specified as being writable upon "
1704  "initialization"));
1705  (*nonlocal_vector)[0][my_row] += values[i];
1706  compressed = false;
1707  }
1708  }
1709  }
1710 
1711 
1712 
1713  inline Vector::size_type
1714  Vector::size() const
1715  {
1716 # ifndef DEAL_II_WITH_64BIT_INDICES
1717  return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1718 # else
1719  return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1720 # endif
1721  }
1722 
1723 
1724 
1725  inline Vector::size_type
1726  Vector::local_size() const
1727  {
1728  return vector->Map().NumMyElements();
1729  }
1730 
1731 
1732 
1733  inline std::pair<Vector::size_type, Vector::size_type>
1734  Vector::local_range() const
1735  {
1736 # ifndef DEAL_II_WITH_64BIT_INDICES
1737  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1739  vector->Map().MaxMyGID() + 1;
1740 # else
1741  const TrilinosWrappers::types::int_type begin =
1742  vector->Map().MinMyGID64();
1744  vector->Map().MaxMyGID64() + 1;
1745 # endif
1746 
1747  Assert(
1748  end - begin == vector->Map().NumMyElements(),
1749  ExcMessage(
1750  "This function only makes sense if the elements that this "
1751  "vector stores on the current processor form a contiguous range. "
1752  "This does not appear to be the case for the current vector."));
1753 
1754  return std::make_pair(begin, end);
1755  }
1756 
1757 
1758 
1759  inline TrilinosScalar Vector::operator*(const Vector &vec) const
1760  {
1761  Assert(vector->Map().SameAs(vec.vector->Map()),
1762  ExcDifferentParallelPartitioning());
1763  Assert(!has_ghost_elements(), ExcGhostsPresent());
1764 
1765  TrilinosScalar result;
1766 
1767  const int ierr = vector->Dot(*(vec.vector), &result);
1768  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1769 
1770  return result;
1771  }
1772 
1773 
1774 
1775  inline Vector::real_type
1776  Vector::norm_sqr() const
1777  {
1778  const TrilinosScalar d = l2_norm();
1779  return d * d;
1780  }
1781 
1782 
1783 
1784  inline TrilinosScalar
1785  Vector::mean_value() const
1786  {
1787  Assert(!has_ghost_elements(), ExcGhostsPresent());
1788 
1790  const int ierr = vector->MeanValue(&mean);
1791  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1792 
1793  return mean;
1794  }
1795 
1796 
1797 
1798  inline TrilinosScalar
1799  Vector::min() const
1800  {
1801  TrilinosScalar min_value;
1802  const int ierr = vector->MinValue(&min_value);
1803  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1804 
1805  return min_value;
1806  }
1807 
1808 
1809 
1810  inline TrilinosScalar
1811  Vector::max() const
1812  {
1813  TrilinosScalar max_value;
1814  const int ierr = vector->MaxValue(&max_value);
1815  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1816 
1817  return max_value;
1818  }
1819 
1820 
1821 
1822  inline Vector::real_type
1823  Vector::l1_norm() const
1824  {
1825  Assert(!has_ghost_elements(), ExcGhostsPresent());
1826 
1827  TrilinosScalar d;
1828  const int ierr = vector->Norm1(&d);
1829  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1830 
1831  return d;
1832  }
1833 
1834 
1835 
1836  inline Vector::real_type
1837  Vector::l2_norm() const
1838  {
1839  Assert(!has_ghost_elements(), ExcGhostsPresent());
1840 
1841  TrilinosScalar d;
1842  const int ierr = vector->Norm2(&d);
1843  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1844 
1845  return d;
1846  }
1847 
1848 
1849 
1850  inline Vector::real_type
1851  Vector::lp_norm(const TrilinosScalar p) const
1852  {
1853  Assert(!has_ghost_elements(), ExcGhostsPresent());
1854 
1855  TrilinosScalar norm = 0;
1856  TrilinosScalar sum = 0;
1857  const size_type n_local = local_size();
1858 
1859  // loop over all the elements because
1860  // Trilinos does not support lp norms
1861  for (size_type i = 0; i < n_local; ++i)
1862  sum += std::pow(std::fabs((*vector)[0][i]), p);
1863 
1864  norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1865 
1866  return norm;
1867  }
1868 
1869 
1870 
1871  inline Vector::real_type
1872  Vector::linfty_norm() const
1873  {
1874  // while we disallow the other
1875  // norm operations on ghosted
1876  // vectors, this particular norm
1877  // is safe to run even in the
1878  // presence of ghost elements
1879  TrilinosScalar d;
1880  const int ierr = vector->NormInf(&d);
1881  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1882 
1883  return d;
1884  }
1885 
1886 
1887 
1888  inline TrilinosScalar
1890  const Vector & V,
1891  const Vector & W)
1892  {
1893  this->add(a, V);
1894  return *this * W;
1895  }
1896 
1897 
1898 
1899  // inline also scalar products, vector
1900  // additions etc. since they are all
1901  // representable by a single Trilinos
1902  // call. This reduces the overhead of the
1903  // wrapper class.
1904  inline Vector &
1906  {
1907  AssertIsFinite(a);
1908 
1909  const int ierr = vector->Scale(a);
1910  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1911 
1912  return *this;
1913  }
1914 
1915 
1916 
1917  inline Vector &
1919  {
1920  AssertIsFinite(a);
1921 
1922  const TrilinosScalar factor = 1. / a;
1923 
1924  AssertIsFinite(factor);
1925 
1926  const int ierr = vector->Scale(factor);
1927  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1928 
1929  return *this;
1930  }
1931 
1932 
1933 
1934  inline Vector &
1935  Vector::operator+=(const Vector &v)
1936  {
1937  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
1938  Assert(vector->Map().SameAs(v.vector->Map()),
1939  ExcDifferentParallelPartitioning());
1940 
1941  const int ierr = vector->Update(1.0, *(v.vector), 1.0);
1942  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1943 
1944  return *this;
1945  }
1946 
1947 
1948 
1949  inline Vector &
1950  Vector::operator-=(const Vector &v)
1951  {
1952  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
1953  Assert(vector->Map().SameAs(v.vector->Map()),
1954  ExcDifferentParallelPartitioning());
1955 
1956  const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
1957  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1958 
1959  return *this;
1960  }
1961 
1962 
1963 
1964  inline void
1965  Vector::add(const TrilinosScalar s)
1966  {
1967  // if we have ghost values, do not allow
1968  // writing to this vector at all.
1969  Assert(!has_ghost_elements(), ExcGhostsPresent());
1970  AssertIsFinite(s);
1971 
1972  size_type n_local = local_size();
1973  for (size_type i = 0; i < n_local; ++i)
1974  (*vector)[0][i] += s;
1975  }
1976 
1977 
1978 
1979  inline void
1980  Vector::add(const TrilinosScalar a, const Vector &v)
1981  {
1982  // if we have ghost values, do not allow
1983  // writing to this vector at all.
1984  Assert(!has_ghost_elements(), ExcGhostsPresent());
1985  Assert(local_size() == v.local_size(),
1986  ExcDimensionMismatch(local_size(), v.local_size()));
1987 
1988  AssertIsFinite(a);
1989 
1990  const int ierr = vector->Update(a, *(v.vector), 1.);
1991  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1992  }
1993 
1994 
1995 
1996  inline void
1997  Vector::add(const TrilinosScalar a,
1998  const Vector & v,
1999  const TrilinosScalar b,
2000  const Vector & w)
2001  {
2002  // if we have ghost values, do not allow
2003  // writing to this vector at all.
2004  Assert(!has_ghost_elements(), ExcGhostsPresent());
2005  Assert(local_size() == v.local_size(),
2006  ExcDimensionMismatch(local_size(), v.local_size()));
2007  Assert(local_size() == w.local_size(),
2008  ExcDimensionMismatch(local_size(), w.local_size()));
2009 
2010  AssertIsFinite(a);
2011  AssertIsFinite(b);
2012 
2013  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2014 
2015  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2016  }
2017 
2018 
2019 
2020  inline void
2021  Vector::sadd(const TrilinosScalar s, const Vector &v)
2022  {
2023  // if we have ghost values, do not allow
2024  // writing to this vector at all.
2025  Assert(!has_ghost_elements(), ExcGhostsPresent());
2026  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2027 
2028  AssertIsFinite(s);
2029 
2030  // We assume that the vectors have the same Map
2031  // if the local size is the same and if the vectors are not ghosted
2032  if (local_size() == v.local_size() && !v.has_ghost_elements())
2033  {
2034  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2035  ExcDifferentParallelPartitioning());
2036  const int ierr = vector->Update(1., *(v.vector), s);
2037  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2038  }
2039  else
2040  {
2041  (*this) *= s;
2042  this->add(v, true);
2043  }
2044  }
2045 
2046 
2047 
2048  inline void
2049  Vector::sadd(const TrilinosScalar s,
2050  const TrilinosScalar a,
2051  const Vector & v)
2052  {
2053  // if we have ghost values, do not allow
2054  // writing to this vector at all.
2055  Assert(!has_ghost_elements(), ExcGhostsPresent());
2056  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2057  AssertIsFinite(s);
2058  AssertIsFinite(a);
2059 
2060  // We assume that the vectors have the same Map
2061  // if the local size is the same and if the vectors are not ghosted
2062  if (local_size() == v.local_size() && !v.has_ghost_elements())
2063  {
2064  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2065  ExcDifferentParallelPartitioning());
2066  const int ierr = vector->Update(a, *(v.vector), s);
2067  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2068  }
2069  else
2070  {
2071  (*this) *= s;
2072  Vector tmp = v;
2073  tmp *= a;
2074  this->add(tmp, true);
2075  }
2076  }
2077 
2078 
2079 
2080  inline void
2081  Vector::scale(const Vector &factors)
2082  {
2083  // if we have ghost values, do not allow
2084  // writing to this vector at all.
2085  Assert(!has_ghost_elements(), ExcGhostsPresent());
2086  Assert(local_size() == factors.local_size(),
2087  ExcDimensionMismatch(local_size(), factors.local_size()));
2088 
2089  const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2090  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2091  }
2092 
2093 
2094 
2095  inline void
2096  Vector::equ(const TrilinosScalar a, const Vector &v)
2097  {
2098  // if we have ghost values, do not allow
2099  // writing to this vector at all.
2100  Assert(!has_ghost_elements(), ExcGhostsPresent());
2101  AssertIsFinite(a);
2102 
2103  // If we don't have the same map, copy.
2104  if (vector->Map().SameAs(v.vector->Map()) == false)
2105  {
2106  this->sadd(0., a, v);
2107  }
2108  else
2109  {
2110  // Otherwise, just update
2111  int ierr = vector->Update(a, *v.vector, 0.0);
2112  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2113 
2114  last_action = Zero;
2115  }
2116  }
2117 
2118 
2119 
2120  inline const Epetra_MultiVector &
2121  Vector::trilinos_vector() const
2122  {
2123  return static_cast<const Epetra_MultiVector &>(*vector);
2124  }
2125 
2126 
2127 
2128  inline Epetra_FEVector &
2129  Vector::trilinos_vector()
2130  {
2131  return *vector;
2132  }
2133 
2134 
2135 
2136  inline const Epetra_BlockMap &
2137  Vector::trilinos_partitioner() const
2138  {
2139  return vector->Map();
2140  }
2141 
2142 
2143 
2144  inline const MPI_Comm &
2145  Vector::get_mpi_communicator() const
2146  {
2147  static MPI_Comm comm;
2148 
2149 # ifdef DEAL_II_WITH_MPI
2150 
2151  const Epetra_MpiComm *mpi_comm =
2152  dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2153  comm = mpi_comm->Comm();
2154 
2155 # else
2156 
2157  comm = MPI_COMM_SELF;
2158 
2159 # endif
2160 
2161  return comm;
2162  }
2163 
2164  template <typename number>
2165  Vector::Vector(const IndexSet & parallel_partitioner,
2166  const ::Vector<number> &v,
2167  const MPI_Comm & communicator)
2168  {
2169  *this =
2170  Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2171  owned_elements = parallel_partitioner;
2172  }
2173 
2174 
2175 
2176  inline Vector &
2178  {
2179  AssertIsFinite(s);
2180 
2181  int ierr = vector->PutScalar(s);
2182  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2183 
2184  if (nonlocal_vector.get() != nullptr)
2185  {
2186  ierr = nonlocal_vector->PutScalar(0.);
2187  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2188  }
2189 
2190  return *this;
2191  }
2192  } /* end of namespace MPI */
2193 
2194 # endif /* DOXYGEN */
2195 
2196 } /* end of namespace TrilinosWrappers */
2197 
2201 namespace internal
2202 {
2203  namespace LinearOperatorImplementation
2204  {
2205  template <typename>
2206  class ReinitHelper;
2207 
2212  template <>
2214  {
2215  public:
2216  template <typename Matrix>
2217  static void
2218  reinit_range_vector(const Matrix & matrix,
2220  bool omit_zeroing_entries)
2221  {
2222  v.reinit(matrix.locally_owned_range_indices(),
2223  matrix.get_mpi_communicator(),
2224  omit_zeroing_entries);
2225  }
2226 
2227  template <typename Matrix>
2228  static void
2231  bool omit_zeroing_entries)
2232  {
2233  v.reinit(matrix.locally_owned_domain_indices(),
2234  matrix.get_mpi_communicator(),
2235  omit_zeroing_entries);
2236  }
2237  };
2238 
2239  } // namespace LinearOperatorImplementation
2240 } /* namespace internal */
2241 
2242 
2243 
2249 template <>
2251 {};
2252 
2253 
2255 
2256 #endif // DEAL_II_WITH_TRILINOS
2257 
2258 /*---------------------------- trilinos_vector.h ---------------------------*/
2259 
2260 #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:621
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:837
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:548
Vector< Number > & operator=(const Vector< Number > &in_vector)
virtual Vector< Number > & operator*=(const Number factor) override
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
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:467
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:518
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1419
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
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:2732
#define DeclException0(Exception0)
Definition: exceptions.h:473
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
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:158
virtual size_type size() const override
Definition: la_vector.h:458
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
internal::VectorReference reference
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
VectorType::value_type * begin(VectorType &V)
T min(const T &t, const MPI_Comm &mpi_communicator)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:587
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2705
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
static const bool value
bool has_ghost_elements() const
#define AssertIsFinite(number)
Definition: exceptions.h:1681
T max(const T &t, const MPI_Comm &mpi_communicator)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)