Reference documentation for deal.II version Git 41b7f7ec9f 2019-11-14 19:01:26 -0600
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
trilinos_vector.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2018 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>
25 # include <deal.II/base/subscriptor.h>
26 # include <deal.II/base/utilities.h>
27 
28 # include <deal.II/lac/exceptions.h>
29 # include <deal.II/lac/vector.h>
30 # include <deal.II/lac/vector_operation.h>
31 # include <deal.II/lac/vector_type_traits.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 
48 DEAL_II_NAMESPACE_OPEN
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  {
89  using size_type = ::types::global_dof_index;
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:
408  using value_type = TrilinosScalar;
409  using real_type = TrilinosScalar;
410  using size_type = ::types::global_dof_index;
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(
692  const ::TrilinosWrappers::SparseMatrix &matrix,
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 
826  TrilinosScalar
827  mean_value() const;
828 
832  TrilinosScalar
833  min() const;
834 
838  TrilinosScalar
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 
886  TrilinosScalar
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 
929  TrilinosScalar
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 
1015  const_iterator
1016  begin() const;
1017 
1022  iterator
1023  end();
1024 
1029  const_iterator
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 
1208  DEAL_II_DEPRECATED
1209  const Epetra_Map &
1210  vector_partitioner() const;
1211 
1216  const Epetra_BlockMap &
1217  trilinos_partitioner() const;
1218 
1226  void
1227  print(std::ostream & out,
1228  const unsigned int precision = 3,
1229  const bool scientific = true,
1230  const bool across = true) const;
1231 
1245  void
1246  swap(Vector &v);
1247 
1251  std::size_t
1252  memory_consumption() const;
1253 
1258  const MPI_Comm &
1259  get_mpi_communicator() const;
1261 
1265  DeclException0(ExcDifferentParallelPartitioning);
1266 
1271  int,
1272  << "An error with error number " << arg1
1273  << " occurred while calling a Trilinos function");
1274 
1279  ExcAccessToNonLocalElement,
1280  size_type,
1281  size_type,
1282  size_type,
1283  size_type,
1284  << "You tried to access element " << arg1
1285  << " of a distributed vector, but this element is not stored "
1286  << "on the current processor. Note: There are " << arg2
1287  << " elements stored "
1288  << "on the current processor from within the range " << arg3
1289  << " through " << arg4
1290  << " but Trilinos vectors need not store contiguous "
1291  << "ranges on each processor, and not every element in "
1292  << "this range may in fact be stored locally.");
1293 
1294  private:
1306  Epetra_CombineMode last_action;
1307 
1313 
1319 
1325  std::unique_ptr<Epetra_FEVector> vector;
1326 
1332  std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1333 
1338 
1339  // Make the reference class a friend.
1340  friend class internal::VectorReference;
1341  };
1342 
1343 
1344 
1345  // ------------------- inline and template functions --------------
1346 
1347 
1356  inline void
1358  {
1359  u.swap(v);
1360  }
1361  } // namespace MPI
1362 
1363 # ifndef DOXYGEN
1364 
1365  namespace internal
1366  {
1367  inline VectorReference::VectorReference(MPI::Vector & vector,
1368  const size_type index)
1369  : vector(vector)
1370  , index(index)
1371  {}
1372 
1373 
1374  inline const VectorReference &
1375  VectorReference::operator=(const VectorReference &r) const
1376  {
1377  // as explained in the class
1378  // documentation, this is not the copy
1379  // operator. so simply pass on to the
1380  // "correct" assignment operator
1381  *this = static_cast<TrilinosScalar>(r);
1382 
1383  return *this;
1384  }
1385 
1386 
1387 
1388  inline VectorReference &
1389  VectorReference::operator=(const VectorReference &r)
1390  {
1391  // as above
1392  *this = static_cast<TrilinosScalar>(r);
1393 
1394  return *this;
1395  }
1396 
1397 
1398  inline const VectorReference &
1399  VectorReference::operator=(const TrilinosScalar &value) const
1400  {
1401  vector.set(1, &index, &value);
1402  return *this;
1403  }
1404 
1405 
1406 
1407  inline const VectorReference &
1408  VectorReference::operator+=(const TrilinosScalar &value) const
1409  {
1410  vector.add(1, &index, &value);
1411  return *this;
1412  }
1413 
1414 
1415 
1416  inline const VectorReference &
1417  VectorReference::operator-=(const TrilinosScalar &value) const
1418  {
1419  TrilinosScalar new_value = -value;
1420  vector.add(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 
1434 
1435 
1436  inline const VectorReference &
1437  VectorReference::operator/=(const TrilinosScalar &value) const
1438  {
1439  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1440  vector.set(1, &index, &new_value);
1441  return *this;
1442  }
1443  } // namespace internal
1444 
1445  namespace MPI
1446  {
1447  inline bool
1448  Vector::in_local_range(const size_type index) const
1449  {
1450  std::pair<size_type, size_type> range = local_range();
1451 
1452  return ((index >= range.first) && (index < range.second));
1453  }
1454 
1455 
1456 
1457  inline IndexSet
1459  {
1460  Assert(owned_elements.size() == size(),
1461  ExcMessage(
1462  "The locally owned elements have not been properly initialized!"
1463  " This happens for example if this object has been initialized"
1464  " with exactly one overlapping IndexSet."));
1465  return owned_elements;
1466  }
1467 
1468 
1469 
1470  inline bool
1472  {
1473  return has_ghosts;
1474  }
1475 
1476 
1477 
1478  inline void
1480  {}
1481 
1482 
1483 
1484  inline internal::VectorReference
1485  Vector::operator()(const size_type index)
1486  {
1487  return internal::VectorReference(*this, index);
1488  }
1489 
1490 
1491 
1492  inline internal::VectorReference Vector::operator[](const size_type index)
1493  {
1494  return operator()(index);
1495  }
1496 
1497 
1498 
1499  inline TrilinosScalar Vector::operator[](const size_type index) const
1500  {
1501  return operator()(index);
1502  }
1503 
1504 
1505 
1506  inline void
1507  Vector::extract_subvector_to(const std::vector<size_type> &indices,
1508  std::vector<TrilinosScalar> & values) const
1509  {
1510  for (size_type i = 0; i < indices.size(); ++i)
1511  values[i] = operator()(indices[i]);
1512  }
1513 
1514 
1515 
1516  template <typename ForwardIterator, typename OutputIterator>
1517  inline void
1518  Vector::extract_subvector_to(ForwardIterator indices_begin,
1519  const ForwardIterator indices_end,
1520  OutputIterator values_begin) const
1521  {
1522  while (indices_begin != indices_end)
1523  {
1524  *values_begin = operator()(*indices_begin);
1525  indices_begin++;
1526  values_begin++;
1527  }
1528  }
1529 
1530 
1531 
1532  inline Vector::iterator
1533  Vector::begin()
1534  {
1535  return (*vector)[0];
1536  }
1537 
1538 
1539 
1540  inline Vector::iterator
1541  Vector::end()
1542  {
1543  return (*vector)[0] + local_size();
1544  }
1545 
1546 
1547 
1548  inline Vector::const_iterator
1549  Vector::begin() const
1550  {
1551  return (*vector)[0];
1552  }
1553 
1554 
1555 
1556  inline Vector::const_iterator
1557  Vector::end() const
1558  {
1559  return (*vector)[0] + local_size();
1560  }
1561 
1562 
1563 
1564  inline void
1565  Vector::set(const std::vector<size_type> & indices,
1566  const std::vector<TrilinosScalar> &values)
1567  {
1568  // if we have ghost values, do not allow
1569  // writing to this vector at all.
1570  Assert(!has_ghost_elements(), ExcGhostsPresent());
1571 
1572  Assert(indices.size() == values.size(),
1573  ExcDimensionMismatch(indices.size(), values.size()));
1574 
1575  set(indices.size(), indices.data(), values.data());
1576  }
1577 
1578 
1579 
1580  inline void
1581  Vector::set(const std::vector<size_type> & indices,
1582  const ::Vector<TrilinosScalar> &values)
1583  {
1584  // if we have ghost values, do not allow
1585  // writing to this vector at all.
1586  Assert(!has_ghost_elements(), ExcGhostsPresent());
1587 
1588  Assert(indices.size() == values.size(),
1589  ExcDimensionMismatch(indices.size(), values.size()));
1590 
1591  set(indices.size(), indices.data(), values.begin());
1592  }
1593 
1594 
1595 
1596  inline void
1597  Vector::set(const size_type n_elements,
1598  const size_type * indices,
1599  const TrilinosScalar *values)
1600  {
1601  // if we have ghost values, do not allow
1602  // writing to this vector at all.
1603  Assert(!has_ghost_elements(), ExcGhostsPresent());
1604 
1605  if (last_action == Add)
1606  {
1607  const int ierr = vector->GlobalAssemble(Add);
1608  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1609  }
1610 
1611  if (last_action != Insert)
1612  last_action = Insert;
1613 
1614  for (size_type i = 0; i < n_elements; ++i)
1615  {
1616  const TrilinosWrappers::types::int_type row = indices[i];
1617  const TrilinosWrappers::types::int_type local_row =
1618  vector->Map().LID(row);
1619  if (local_row != -1)
1620  (*vector)[0][local_row] = values[i];
1621  else
1622  {
1623  const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1624  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1625  compressed = false;
1626  }
1627  // in set operation, do not use the pre-allocated vector for nonlocal
1628  // entries even if it exists. This is to ensure that we really only
1629  // set the elements touched by the set() method and not all contained
1630  // in the nonlocal entries vector (there is no way to distinguish them
1631  // on the receiving processor)
1632  }
1633  }
1634 
1635 
1636 
1637  inline void
1638  Vector::add(const std::vector<size_type> & indices,
1639  const std::vector<TrilinosScalar> &values)
1640  {
1641  // if we have ghost values, do not allow
1642  // writing to this vector at all.
1643  Assert(!has_ghost_elements(), ExcGhostsPresent());
1644  Assert(indices.size() == values.size(),
1645  ExcDimensionMismatch(indices.size(), values.size()));
1646 
1647  add(indices.size(), indices.data(), values.data());
1648  }
1649 
1650 
1651 
1652  inline void
1653  Vector::add(const std::vector<size_type> & indices,
1654  const ::Vector<TrilinosScalar> &values)
1655  {
1656  // if we have ghost values, do not allow
1657  // writing to this vector at all.
1658  Assert(!has_ghost_elements(), ExcGhostsPresent());
1659  Assert(indices.size() == values.size(),
1660  ExcDimensionMismatch(indices.size(), values.size()));
1661 
1662  add(indices.size(), indices.data(), values.begin());
1663  }
1664 
1665 
1666 
1667  inline void
1668  Vector::add(const size_type n_elements,
1669  const size_type * indices,
1670  const TrilinosScalar *values)
1671  {
1672  // if we have ghost values, do not allow
1673  // writing to this vector at all.
1674  Assert(!has_ghost_elements(), ExcGhostsPresent());
1675 
1676  if (last_action != Add)
1677  {
1678  if (last_action == Insert)
1679  {
1680  const int ierr = vector->GlobalAssemble(Insert);
1681  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1682  }
1683  last_action = Add;
1684  }
1685 
1686  for (size_type i = 0; i < n_elements; ++i)
1687  {
1688  const size_type row = indices[i];
1689  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1690  static_cast<TrilinosWrappers::types::int_type>(row));
1691  if (local_row != -1)
1692  (*vector)[0][local_row] += values[i];
1693  else if (nonlocal_vector.get() == nullptr)
1694  {
1695  const int ierr = vector->SumIntoGlobalValues(
1696  1,
1697  reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1698  &row),
1699  &values[i]);
1700  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1701  compressed = false;
1702  }
1703  else
1704  {
1705  // use pre-allocated vector for non-local entries if it exists for
1706  // addition operation
1707  const TrilinosWrappers::types::int_type my_row =
1708  nonlocal_vector->Map().LID(
1709  static_cast<TrilinosWrappers::types::int_type>(row));
1710  Assert(my_row != -1,
1711  ExcMessage(
1712  "Attempted to write into off-processor vector entry "
1713  "that has not be specified as being writable upon "
1714  "initialization"));
1715  (*nonlocal_vector)[0][my_row] += values[i];
1716  compressed = false;
1717  }
1718  }
1719  }
1720 
1721 
1722 
1723  inline Vector::size_type
1724  Vector::size() const
1725  {
1726 # ifndef DEAL_II_WITH_64BIT_INDICES
1727  return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1728 # else
1729  return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1730 # endif
1731  }
1732 
1733 
1734 
1735  inline Vector::size_type
1736  Vector::local_size() const
1737  {
1738  return vector->Map().NumMyElements();
1739  }
1740 
1741 
1742 
1743  inline std::pair<Vector::size_type, Vector::size_type>
1744  Vector::local_range() const
1745  {
1746 # ifndef DEAL_II_WITH_64BIT_INDICES
1747  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1748  const TrilinosWrappers::types::int_type end =
1749  vector->Map().MaxMyGID() + 1;
1750 # else
1751  const TrilinosWrappers::types::int_type begin =
1752  vector->Map().MinMyGID64();
1753  const TrilinosWrappers::types::int_type end =
1754  vector->Map().MaxMyGID64() + 1;
1755 # endif
1756 
1757  Assert(
1758  end - begin == vector->Map().NumMyElements(),
1759  ExcMessage(
1760  "This function only makes sense if the elements that this "
1761  "vector stores on the current processor form a contiguous range. "
1762  "This does not appear to be the case for the current vector."));
1763 
1764  return std::make_pair(begin, end);
1765  }
1766 
1767 
1768 
1769  inline TrilinosScalar Vector::operator*(const Vector &vec) const
1770  {
1771  Assert(vector->Map().SameAs(vec.vector->Map()),
1772  ExcDifferentParallelPartitioning());
1773  Assert(!has_ghost_elements(), ExcGhostsPresent());
1774 
1775  TrilinosScalar result;
1776 
1777  const int ierr = vector->Dot(*(vec.vector), &result);
1778  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1779 
1780  return result;
1781  }
1782 
1783 
1784 
1785  inline Vector::real_type
1786  Vector::norm_sqr() const
1787  {
1788  const TrilinosScalar d = l2_norm();
1789  return d * d;
1790  }
1791 
1792 
1793 
1794  inline TrilinosScalar
1795  Vector::mean_value() const
1796  {
1797  Assert(!has_ghost_elements(), ExcGhostsPresent());
1798 
1799  TrilinosScalar mean;
1800  const int ierr = vector->MeanValue(&mean);
1801  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1802 
1803  return mean;
1804  }
1805 
1806 
1807 
1808  inline TrilinosScalar
1809  Vector::min() const
1810  {
1811  TrilinosScalar min_value;
1812  const int ierr = vector->MinValue(&min_value);
1813  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1814 
1815  return min_value;
1816  }
1817 
1818 
1819 
1820  inline TrilinosScalar
1821  Vector::max() const
1822  {
1823  TrilinosScalar max_value;
1824  const int ierr = vector->MaxValue(&max_value);
1825  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1826 
1827  return max_value;
1828  }
1829 
1830 
1831 
1832  inline Vector::real_type
1833  Vector::l1_norm() const
1834  {
1835  Assert(!has_ghost_elements(), ExcGhostsPresent());
1836 
1837  TrilinosScalar d;
1838  const int ierr = vector->Norm1(&d);
1839  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1840 
1841  return d;
1842  }
1843 
1844 
1845 
1846  inline Vector::real_type
1847  Vector::l2_norm() const
1848  {
1849  Assert(!has_ghost_elements(), ExcGhostsPresent());
1850 
1851  TrilinosScalar d;
1852  const int ierr = vector->Norm2(&d);
1853  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1854 
1855  return d;
1856  }
1857 
1858 
1859 
1860  inline Vector::real_type
1861  Vector::lp_norm(const TrilinosScalar p) const
1862  {
1863  Assert(!has_ghost_elements(), ExcGhostsPresent());
1864 
1865  TrilinosScalar norm = 0;
1866  TrilinosScalar sum = 0;
1867  const size_type n_local = local_size();
1868 
1869  // loop over all the elements because
1870  // Trilinos does not support lp norms
1871  for (size_type i = 0; i < n_local; ++i)
1872  sum += std::pow(std::fabs((*vector)[0][i]), p);
1873 
1874  norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1875 
1876  return norm;
1877  }
1878 
1879 
1880 
1881  inline Vector::real_type
1882  Vector::linfty_norm() const
1883  {
1884  // while we disallow the other
1885  // norm operations on ghosted
1886  // vectors, this particular norm
1887  // is safe to run even in the
1888  // presence of ghost elements
1889  TrilinosScalar d;
1890  const int ierr = vector->NormInf(&d);
1891  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1892 
1893  return d;
1894  }
1895 
1896 
1897 
1898  inline TrilinosScalar
1899  Vector::add_and_dot(const TrilinosScalar a,
1900  const Vector & V,
1901  const Vector & W)
1902  {
1903  this->add(a, V);
1904  return *this * W;
1905  }
1906 
1907 
1908 
1909  // inline also scalar products, vector
1910  // additions etc. since they are all
1911  // representable by a single Trilinos
1912  // call. This reduces the overhead of the
1913  // wrapper class.
1914  inline Vector &
1915  Vector::operator*=(const TrilinosScalar a)
1916  {
1917  AssertIsFinite(a);
1918 
1919  const int ierr = vector->Scale(a);
1920  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1921 
1922  return *this;
1923  }
1924 
1925 
1926 
1927  inline Vector &
1928  Vector::operator/=(const TrilinosScalar a)
1929  {
1930  AssertIsFinite(a);
1931 
1932  const TrilinosScalar factor = 1. / a;
1933 
1934  AssertIsFinite(factor);
1935 
1936  const int ierr = vector->Scale(factor);
1937  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1938 
1939  return *this;
1940  }
1941 
1942 
1943 
1944  inline Vector &
1945  Vector::operator+=(const Vector &v)
1946  {
1947  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
1948  Assert(vector->Map().SameAs(v.vector->Map()),
1949  ExcDifferentParallelPartitioning());
1950 
1951  const int ierr = vector->Update(1.0, *(v.vector), 1.0);
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  Assert(size() == v.size(), ExcDimensionMismatch(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 void
1975  Vector::add(const TrilinosScalar s)
1976  {
1977  // if we have ghost values, do not allow
1978  // writing to this vector at all.
1979  Assert(!has_ghost_elements(), ExcGhostsPresent());
1980  AssertIsFinite(s);
1981 
1982  size_type n_local = local_size();
1983  for (size_type i = 0; i < n_local; ++i)
1984  (*vector)[0][i] += s;
1985  }
1986 
1987 
1988 
1989  inline void
1990  Vector::add(const TrilinosScalar a, const Vector &v)
1991  {
1992  // if we have ghost values, do not allow
1993  // writing to this vector at all.
1994  Assert(!has_ghost_elements(), ExcGhostsPresent());
1995  Assert(local_size() == v.local_size(),
1996  ExcDimensionMismatch(local_size(), v.local_size()));
1997 
1998  AssertIsFinite(a);
1999 
2000  const int ierr = vector->Update(a, *(v.vector), 1.);
2001  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2002  }
2003 
2004 
2005 
2006  inline void
2007  Vector::add(const TrilinosScalar a,
2008  const Vector & v,
2009  const TrilinosScalar b,
2010  const Vector & w)
2011  {
2012  // if we have ghost values, do not allow
2013  // writing to this vector at all.
2014  Assert(!has_ghost_elements(), ExcGhostsPresent());
2015  Assert(local_size() == v.local_size(),
2016  ExcDimensionMismatch(local_size(), v.local_size()));
2017  Assert(local_size() == w.local_size(),
2018  ExcDimensionMismatch(local_size(), w.local_size()));
2019 
2020  AssertIsFinite(a);
2021  AssertIsFinite(b);
2022 
2023  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2024 
2025  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2026  }
2027 
2028 
2029 
2030  inline void
2031  Vector::sadd(const TrilinosScalar s, const Vector &v)
2032  {
2033  // if we have ghost values, do not allow
2034  // writing to this vector at all.
2035  Assert(!has_ghost_elements(), ExcGhostsPresent());
2036  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2037 
2038  AssertIsFinite(s);
2039 
2040  // We assume that the vectors have the same Map
2041  // if the local size is the same and if the vectors are not ghosted
2042  if (local_size() == v.local_size() && !v.has_ghost_elements())
2043  {
2044  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2045  ExcDifferentParallelPartitioning());
2046  const int ierr = vector->Update(1., *(v.vector), s);
2047  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2048  }
2049  else
2050  {
2051  (*this) *= s;
2052  this->add(v, true);
2053  }
2054  }
2055 
2056 
2057 
2058  inline void
2059  Vector::sadd(const TrilinosScalar s,
2060  const TrilinosScalar a,
2061  const Vector & v)
2062  {
2063  // if we have ghost values, do not allow
2064  // writing to this vector at all.
2065  Assert(!has_ghost_elements(), ExcGhostsPresent());
2066  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2067  AssertIsFinite(s);
2068  AssertIsFinite(a);
2069 
2070  // We assume that the vectors have the same Map
2071  // if the local size is the same and if the vectors are not ghosted
2072  if (local_size() == v.local_size() && !v.has_ghost_elements())
2073  {
2074  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2075  ExcDifferentParallelPartitioning());
2076  const int ierr = vector->Update(a, *(v.vector), s);
2077  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2078  }
2079  else
2080  {
2081  (*this) *= s;
2082  Vector tmp = v;
2083  tmp *= a;
2084  this->add(tmp, true);
2085  }
2086  }
2087 
2088 
2089 
2090  inline void
2091  Vector::scale(const Vector &factors)
2092  {
2093  // if we have ghost values, do not allow
2094  // writing to this vector at all.
2095  Assert(!has_ghost_elements(), ExcGhostsPresent());
2096  Assert(local_size() == factors.local_size(),
2097  ExcDimensionMismatch(local_size(), factors.local_size()));
2098 
2099  const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2100  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2101  }
2102 
2103 
2104 
2105  inline void
2106  Vector::equ(const TrilinosScalar a, const Vector &v)
2107  {
2108  // if we have ghost values, do not allow
2109  // writing to this vector at all.
2110  Assert(!has_ghost_elements(), ExcGhostsPresent());
2111  AssertIsFinite(a);
2112 
2113  // If we don't have the same map, copy.
2114  if (vector->Map().SameAs(v.vector->Map()) == false)
2115  {
2116  this->sadd(0., a, v);
2117  }
2118  else
2119  {
2120  // Otherwise, just update
2121  int ierr = vector->Update(a, *v.vector, 0.0);
2122  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2123 
2124  last_action = Zero;
2125  }
2126  }
2127 
2128 
2129 
2130  inline const Epetra_MultiVector &
2131  Vector::trilinos_vector() const
2132  {
2133  return static_cast<const Epetra_MultiVector &>(*vector);
2134  }
2135 
2136 
2137 
2138  inline Epetra_FEVector &
2139  Vector::trilinos_vector()
2140  {
2141  return *vector;
2142  }
2143 
2144 
2145 
2146  inline const Epetra_Map &
2147  Vector::vector_partitioner() const
2148  {
2149  // TODO A dynamic_cast fails here. This is suspicious.
2150  return static_cast<const Epetra_Map &>(vector->Map()); // NOLINT
2151  }
2152 
2153 
2154 
2155  inline const Epetra_BlockMap &
2156  Vector::trilinos_partitioner() const
2157  {
2158  return vector->Map();
2159  }
2160 
2161 
2162 
2163  inline const MPI_Comm &
2164  Vector::get_mpi_communicator() const
2165  {
2166  static MPI_Comm comm;
2167 
2168 # ifdef DEAL_II_WITH_MPI
2169 
2170  const Epetra_MpiComm *mpi_comm =
2171  dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2172  comm = mpi_comm->Comm();
2173 
2174 # else
2175 
2176  comm = MPI_COMM_SELF;
2177 
2178 # endif
2179 
2180  return comm;
2181  }
2182 
2183  template <typename number>
2184  Vector::Vector(const IndexSet & parallel_partitioner,
2185  const ::Vector<number> &v,
2186  const MPI_Comm & communicator)
2187  {
2188  *this =
2189  Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2190  owned_elements = parallel_partitioner;
2191  }
2192 
2193 
2194 
2195  inline Vector &
2196  Vector::operator=(const TrilinosScalar s)
2197  {
2198  AssertIsFinite(s);
2199 
2200  int ierr = vector->PutScalar(s);
2201  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2202 
2203  if (nonlocal_vector.get() != nullptr)
2204  {
2205  ierr = nonlocal_vector->PutScalar(0.);
2206  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2207  }
2208 
2209  return *this;
2210  }
2211  } /* end of namespace MPI */
2212 
2213 # endif /* DOXYGEN */
2214 
2215 } /* end of namespace TrilinosWrappers */
2216 
2220 namespace internal
2221 {
2222  namespace LinearOperatorImplementation
2223  {
2224  template <typename>
2225  class ReinitHelper;
2226 
2231  template <>
2233  {
2234  public:
2235  template <typename Matrix>
2236  static void
2237  reinit_range_vector(const Matrix & matrix,
2239  bool omit_zeroing_entries)
2240  {
2241  v.reinit(matrix.locally_owned_range_indices(),
2242  matrix.get_mpi_communicator(),
2243  omit_zeroing_entries);
2244  }
2245 
2246  template <typename Matrix>
2247  static void
2248  reinit_domain_vector(const Matrix & matrix,
2250  bool omit_zeroing_entries)
2251  {
2252  v.reinit(matrix.locally_owned_domain_indices(),
2253  matrix.get_mpi_communicator(),
2254  omit_zeroing_entries);
2255  }
2256  };
2257 
2258  } // namespace LinearOperatorImplementation
2259 } /* namespace internal */
2260 
2261 
2262 
2268 template <>
2270 {};
2271 
2272 
2273 DEAL_II_NAMESPACE_CLOSE
2274 
2275 #endif // DEAL_II_WITH_TRILINOS
2276 
2277 /*---------------------------- trilinos_vector.h ---------------------------*/
2278 
2279 #endif // dealii_trilinos_vector_h
std::unique_ptr< Epetra_MultiVector > nonlocal_vector
static ::ExceptionBase & ExcTrilinosError(int arg1)
void update_ghost_values() const
Number operator[](const size_type global_index) const
virtual VectorSpaceVector< Number >::real_type linfty_norm() const override
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
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)
Vector< Number > & operator=(const Vector< Number > &in_vector)
virtual Vector< Number > & operator*=(const Number factor) override
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
virtual Number operator*(const VectorSpaceVector< Number > &V) const override
virtual Vector< Number > & operator+=(const VectorSpaceVector< Number > &V) override
LinearAlgebra::distributed::Vector< Number > Vector
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:565
__global__ void scale(Number *val, const Number *V_val, const size_type N)
virtual ::IndexSet locally_owned_elements() const override
Definition: la_vector.h:467
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)
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
#define Assert(cond, exc)
Definition: exceptions.h:1411
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
#define DeclException0(Exception0)
Definition: exceptions.h:473
virtual Vector< Number > & operator/=(const Number factor) override
virtual void scale(const VectorSpaceVector< Number > &scaling_factors) override
virtual size_type size() const override
Definition: la_vector.h:458
__global__ void equ(Number *val, const Number a, const Number *V_val, const size_type N)
__global__ void sadd(const Number s, Number *val, const Number a, const Number *V_val, const size_type N)
static ::ExceptionBase & ExcGhostsPresent()
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
unsigned int global_dof_index
Definition: types.h:89
__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
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:587
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
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:1673