Reference documentation for deal.II version 9.1.1
\(\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 namespace LinearAlgebra
51 {
52  // Forward declaration
53  template <typename Number>
54  class ReadWriteVector;
55 } // namespace LinearAlgebra
56 
67 namespace TrilinosWrappers
68 {
69  class SparseMatrix;
70 
81  namespace internal
82  {
86  using size_type = ::types::global_dof_index;
87 
97  class VectorReference
98  {
99  private:
104  VectorReference(MPI::Vector &vector, const size_type index);
105 
106  public:
118  const VectorReference &
119  operator=(const VectorReference &r) const;
120 
124  VectorReference &
125  operator=(const VectorReference &r);
126 
130  const VectorReference &
131  operator=(const TrilinosScalar &s) const;
132 
136  const VectorReference &
137  operator+=(const TrilinosScalar &s) const;
138 
142  const VectorReference &
143  operator-=(const TrilinosScalar &s) const;
144 
148  const VectorReference &
149  operator*=(const TrilinosScalar &s) const;
150 
154  const VectorReference &
155  operator/=(const TrilinosScalar &s) const;
156 
161  operator TrilinosScalar() const;
162 
167  int,
168  << "An error with error number " << arg1
169  << " occurred while calling a Trilinos function");
170 
171  private:
175  MPI::Vector &vector;
176 
180  const size_type index;
181 
187  };
188  } // namespace internal
193 # ifndef DEAL_II_WITH_64BIT_INDICES
194  // define a helper function that queries the global ID of local ID of
195  // an Epetra_BlockMap object by calling either the 32- or 64-bit
196  // function necessary.
197  inline int
198  gid(const Epetra_BlockMap &map, int i)
199  {
200  return map.GID(i);
201  }
202 # else
203  // define a helper function that queries the global ID of local ID of
204  // an Epetra_BlockMap object by calling either the 32- or 64-bit
205  // function necessary.
206  inline long long int
207  gid(const Epetra_BlockMap &map, int i)
208  {
209  return map.GID64(i);
210  }
211 # endif
212 
219  namespace MPI
220  {
221  class BlockVector;
222 
399  class Vector : public Subscriptor
400  {
401  public:
407  using value_type = TrilinosScalar;
408  using real_type = TrilinosScalar;
409  using size_type = ::types::global_dof_index;
410  using iterator = value_type *;
411  using const_iterator = const value_type *;
412  using reference = internal::VectorReference;
413  using const_reference = const internal::VectorReference;
414 
424  Vector();
425 
429  Vector(const Vector &v);
430 
449  explicit Vector(const IndexSet &parallel_partitioning,
450  const MPI_Comm &communicator = MPI_COMM_WORLD);
451 
463  Vector(const IndexSet &local,
464  const IndexSet &ghost,
465  const MPI_Comm &communicator = MPI_COMM_WORLD);
466 
481  Vector(const IndexSet &parallel_partitioning,
482  const Vector & v,
483  const MPI_Comm &communicator = MPI_COMM_WORLD);
484 
497  template <typename Number>
498  Vector(const IndexSet & parallel_partitioning,
499  const ::Vector<Number> &v,
500  const MPI_Comm & communicator = MPI_COMM_WORLD);
501 
506  Vector(Vector &&v) noexcept;
507 
511  ~Vector() override = default;
512 
517  void
518  clear();
519 
543  void
544  reinit(const Vector &v,
545  const bool omit_zeroing_entries = false,
546  const bool allow_different_maps = false);
547 
570  void
571  reinit(const IndexSet &parallel_partitioning,
572  const MPI_Comm &communicator = MPI_COMM_WORLD,
573  const bool omit_zeroing_entries = false);
574 
600  void
601  reinit(const IndexSet &locally_owned_entries,
602  const IndexSet &ghost_entries,
603  const MPI_Comm &communicator = MPI_COMM_WORLD,
604  const bool vector_writable = false);
605 
609  void
610  reinit(const BlockVector &v, const bool import_data = false);
611 
628  void
629  compress(::VectorOperation::values operation);
630 
643  Vector &
644  operator=(const TrilinosScalar s);
645 
651  Vector &
652  operator=(const Vector &v);
653 
658  Vector &
659  operator=(Vector &&v) noexcept;
660 
668  template <typename Number>
669  Vector &
670  operator=(const ::Vector<Number> &v);
671 
689  void
690  import_nonlocal_data_for_fe(
691  const ::TrilinosWrappers::SparseMatrix &matrix,
692  const Vector & vector);
693 
700  void
701  import(const LinearAlgebra::ReadWriteVector<double> &rwv,
702  const VectorOperation::values operation);
703 
704 
710  bool
711  operator==(const Vector &v) const;
712 
718  bool
719  operator!=(const Vector &v) const;
720 
724  size_type
725  size() const;
726 
738  size_type
739  local_size() const;
740 
762  std::pair<size_type, size_type>
763  local_range() const;
764 
772  bool
773  in_local_range(const size_type index) const;
774 
788  IndexSet
789  locally_owned_elements() const;
790 
798  bool
799  has_ghost_elements() const;
800 
807  void
808  update_ghost_values() const;
809 
814  TrilinosScalar operator*(const Vector &vec) const;
815 
819  real_type
820  norm_sqr() const;
821 
825  TrilinosScalar
826  mean_value() const;
827 
831  TrilinosScalar
832  min() const;
833 
837  TrilinosScalar
838  max() const;
839 
843  real_type
844  l1_norm() const;
845 
850  real_type
851  l2_norm() const;
852 
857  real_type
858  lp_norm(const TrilinosScalar p) const;
859 
863  real_type
864  linfty_norm() const;
865 
885  TrilinosScalar
886  add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
887 
893  bool
894  all_zero() const;
895 
901  bool
902  is_non_negative() const;
904 
905 
910 
918  reference
919  operator()(const size_type index);
920 
928  TrilinosScalar
929  operator()(const size_type index) const;
930 
936  reference operator[](const size_type index);
937 
943  TrilinosScalar operator[](const size_type index) const;
944 
960  void
961  extract_subvector_to(const std::vector<size_type> &indices,
962  std::vector<TrilinosScalar> & values) const;
963 
991  template <typename ForwardIterator, typename OutputIterator>
992  void
993  extract_subvector_to(ForwardIterator indices_begin,
994  const ForwardIterator indices_end,
995  OutputIterator values_begin) const;
996 
1007  iterator
1008  begin();
1009 
1014  const_iterator
1015  begin() const;
1016 
1021  iterator
1022  end();
1023 
1028  const_iterator
1029  end() const;
1030 
1032 
1033 
1038 
1045  void
1046  set(const std::vector<size_type> & indices,
1047  const std::vector<TrilinosScalar> &values);
1048 
1053  void
1054  set(const std::vector<size_type> & indices,
1055  const ::Vector<TrilinosScalar> &values);
1056 
1062  void
1063  set(const size_type n_elements,
1064  const size_type * indices,
1065  const TrilinosScalar *values);
1066 
1071  void
1072  add(const std::vector<size_type> & indices,
1073  const std::vector<TrilinosScalar> &values);
1074 
1079  void
1080  add(const std::vector<size_type> & indices,
1081  const ::Vector<TrilinosScalar> &values);
1082 
1088  void
1089  add(const size_type n_elements,
1090  const size_type * indices,
1091  const TrilinosScalar *values);
1092 
1096  Vector &
1097  operator*=(const TrilinosScalar factor);
1098 
1102  Vector &
1103  operator/=(const TrilinosScalar factor);
1104 
1108  Vector &
1109  operator+=(const Vector &V);
1110 
1114  Vector &
1115  operator-=(const Vector &V);
1116 
1121  void
1122  add(const TrilinosScalar s);
1123 
1136  void
1137  add(const Vector &V, const bool allow_different_maps = false);
1138 
1142  void
1143  add(const TrilinosScalar a, const Vector &V);
1144 
1148  void
1149  add(const TrilinosScalar a,
1150  const Vector & V,
1151  const TrilinosScalar b,
1152  const Vector & W);
1153 
1158  void
1159  sadd(const TrilinosScalar s, const Vector &V);
1160 
1164  void
1165  sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1166 
1172  void
1173  scale(const Vector &scaling_factors);
1174 
1178  void
1179  equ(const TrilinosScalar a, const Vector &V);
1181 
1186 
1191  const Epetra_MultiVector &
1192  trilinos_vector() const;
1193 
1198  Epetra_FEVector &
1199  trilinos_vector();
1200 
1207  DEAL_II_DEPRECATED
1208  const Epetra_Map &
1209  vector_partitioner() const;
1210 
1215  const Epetra_BlockMap &
1216  trilinos_partitioner() const;
1217 
1225  void
1226  print(std::ostream & out,
1227  const unsigned int precision = 3,
1228  const bool scientific = true,
1229  const bool across = true) const;
1230 
1244  void
1245  swap(Vector &v);
1246 
1250  std::size_t
1251  memory_consumption() const;
1252 
1257  const MPI_Comm &
1258  get_mpi_communicator() const;
1260 
1264  DeclException0(ExcDifferentParallelPartitioning);
1265 
1270  int,
1271  << "An error with error number " << arg1
1272  << " occurred while calling a Trilinos function");
1273 
1278  ExcAccessToNonLocalElement,
1279  size_type,
1280  size_type,
1281  size_type,
1282  size_type,
1283  << "You tried to access element " << arg1
1284  << " of a distributed vector, but this element is not stored "
1285  << "on the current processor. Note: There are " << arg2
1286  << " elements stored "
1287  << "on the current processor from within the range " << arg3
1288  << " through " << arg4
1289  << " but Trilinos vectors need not store contiguous "
1290  << "ranges on each processor, and not every element in "
1291  << "this range may in fact be stored locally.");
1292 
1293  private:
1305  Epetra_CombineMode last_action;
1306 
1312 
1318 
1324  std::unique_ptr<Epetra_FEVector> vector;
1325 
1331  std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1332 
1337 
1341  friend class internal::VectorReference;
1342  };
1343 
1344 
1345 
1346  // ------------------- inline and template functions --------------
1347 
1348 
1357  inline void
1359  {
1360  u.swap(v);
1361  }
1362  } // namespace MPI
1363 
1364 # ifndef DOXYGEN
1365 
1366  namespace internal
1367  {
1368  inline VectorReference::VectorReference(MPI::Vector & vector,
1369  const size_type index)
1370  : vector(vector)
1371  , index(index)
1372  {}
1373 
1374 
1375  inline const VectorReference &
1376  VectorReference::operator=(const VectorReference &r) const
1377  {
1378  // as explained in the class
1379  // documentation, this is not the copy
1380  // operator. so simply pass on to the
1381  // "correct" assignment operator
1382  *this = static_cast<TrilinosScalar>(r);
1383 
1384  return *this;
1385  }
1386 
1387 
1388 
1389  inline VectorReference &
1390  VectorReference::operator=(const VectorReference &r)
1391  {
1392  // as above
1393  *this = static_cast<TrilinosScalar>(r);
1394 
1395  return *this;
1396  }
1397 
1398 
1399  inline const VectorReference &
1400  VectorReference::operator=(const TrilinosScalar &value) const
1401  {
1402  vector.set(1, &index, &value);
1403  return *this;
1404  }
1405 
1406 
1407 
1408  inline const VectorReference &
1409  VectorReference::operator+=(const TrilinosScalar &value) const
1410  {
1411  vector.add(1, &index, &value);
1412  return *this;
1413  }
1414 
1415 
1416 
1417  inline const VectorReference &
1418  VectorReference::operator-=(const TrilinosScalar &value) const
1419  {
1420  TrilinosScalar new_value = -value;
1421  vector.add(1, &index, &new_value);
1422  return *this;
1423  }
1424 
1425 
1426 
1427  inline const VectorReference &
1428  VectorReference::operator*=(const TrilinosScalar &value) const
1429  {
1430  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1431  vector.set(1, &index, &new_value);
1432  return *this;
1433  }
1434 
1435 
1436 
1437  inline const VectorReference &
1438  VectorReference::operator/=(const TrilinosScalar &value) const
1439  {
1440  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1441  vector.set(1, &index, &new_value);
1442  return *this;
1443  }
1444  } // namespace internal
1445 
1446  namespace MPI
1447  {
1448  inline bool
1449  Vector::in_local_range(const size_type index) const
1450  {
1451  std::pair<size_type, size_type> range = local_range();
1452 
1453  return ((index >= range.first) && (index < range.second));
1454  }
1455 
1456 
1457 
1458  inline IndexSet
1460  {
1461  Assert(owned_elements.size() == size(),
1462  ExcMessage(
1463  "The locally owned elements have not been properly initialized!"
1464  " This happens for example if this object has been initialized"
1465  " with exactly one overlapping IndexSet."));
1466  return owned_elements;
1467  }
1468 
1469 
1470 
1471  inline bool
1472  Vector::has_ghost_elements() const
1473  {
1474  return has_ghosts;
1475  }
1476 
1477 
1478 
1479  inline void
1481  {}
1482 
1483 
1484 
1485  inline internal::VectorReference
1486  Vector::operator()(const size_type index)
1487  {
1488  return internal::VectorReference(*this, index);
1489  }
1490 
1491 
1492 
1493  inline internal::VectorReference Vector::operator[](const size_type index)
1494  {
1495  return operator()(index);
1496  }
1497 
1498 
1499 
1500  inline TrilinosScalar Vector::operator[](const size_type index) const
1501  {
1502  return operator()(index);
1503  }
1504 
1505 
1506 
1507  inline void
1508  Vector::extract_subvector_to(const std::vector<size_type> &indices,
1509  std::vector<TrilinosScalar> & values) const
1510  {
1511  for (size_type i = 0; i < indices.size(); ++i)
1512  values[i] = operator()(indices[i]);
1513  }
1514 
1515 
1516 
1517  template <typename ForwardIterator, typename OutputIterator>
1518  inline void
1519  Vector::extract_subvector_to(ForwardIterator indices_begin,
1520  const ForwardIterator indices_end,
1521  OutputIterator values_begin) const
1522  {
1523  while (indices_begin != indices_end)
1524  {
1525  *values_begin = operator()(*indices_begin);
1526  indices_begin++;
1527  values_begin++;
1528  }
1529  }
1530 
1531 
1532 
1533  inline Vector::iterator
1534  Vector::begin()
1535  {
1536  return (*vector)[0];
1537  }
1538 
1539 
1540 
1541  inline Vector::iterator
1542  Vector::end()
1543  {
1544  return (*vector)[0] + local_size();
1545  }
1546 
1547 
1548 
1549  inline Vector::const_iterator
1550  Vector::begin() const
1551  {
1552  return (*vector)[0];
1553  }
1554 
1555 
1556 
1557  inline Vector::const_iterator
1558  Vector::end() const
1559  {
1560  return (*vector)[0] + local_size();
1561  }
1562 
1563 
1564 
1565  inline void
1566  Vector::set(const std::vector<size_type> & indices,
1567  const std::vector<TrilinosScalar> &values)
1568  {
1569  // if we have ghost values, do not allow
1570  // writing to this vector at all.
1571  Assert(!has_ghost_elements(), ExcGhostsPresent());
1572 
1573  Assert(indices.size() == values.size(),
1574  ExcDimensionMismatch(indices.size(), values.size()));
1575 
1576  set(indices.size(), indices.data(), values.data());
1577  }
1578 
1579 
1580 
1581  inline void
1582  Vector::set(const std::vector<size_type> & indices,
1583  const ::Vector<TrilinosScalar> &values)
1584  {
1585  // if we have ghost values, do not allow
1586  // writing to this vector at all.
1587  Assert(!has_ghost_elements(), ExcGhostsPresent());
1588 
1589  Assert(indices.size() == values.size(),
1590  ExcDimensionMismatch(indices.size(), values.size()));
1591 
1592  set(indices.size(), indices.data(), values.begin());
1593  }
1594 
1595 
1596 
1597  inline void
1598  Vector::set(const size_type n_elements,
1599  const size_type * indices,
1600  const TrilinosScalar *values)
1601  {
1602  // if we have ghost values, do not allow
1603  // writing to this vector at all.
1604  Assert(!has_ghost_elements(), ExcGhostsPresent());
1605 
1606  if (last_action == Add)
1607  {
1608  const int ierr = vector->GlobalAssemble(Add);
1609  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1610  }
1611 
1612  if (last_action != Insert)
1613  last_action = Insert;
1614 
1615  for (size_type i = 0; i < n_elements; ++i)
1616  {
1617  const TrilinosWrappers::types::int_type row = indices[i];
1618  const TrilinosWrappers::types::int_type local_row =
1619  vector->Map().LID(row);
1620  if (local_row != -1)
1621  (*vector)[0][local_row] = values[i];
1622  else
1623  {
1624  const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1625  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1626  compressed = false;
1627  }
1628  // in set operation, do not use the pre-allocated vector for nonlocal
1629  // entries even if it exists. This is to ensure that we really only
1630  // set the elements touched by the set() method and not all contained
1631  // in the nonlocal entries vector (there is no way to distinguish them
1632  // on the receiving processor)
1633  }
1634  }
1635 
1636 
1637 
1638  inline void
1639  Vector::add(const std::vector<size_type> & indices,
1640  const std::vector<TrilinosScalar> &values)
1641  {
1642  // if we have ghost values, do not allow
1643  // writing to this vector at all.
1644  Assert(!has_ghost_elements(), ExcGhostsPresent());
1645  Assert(indices.size() == values.size(),
1646  ExcDimensionMismatch(indices.size(), values.size()));
1647 
1648  add(indices.size(), indices.data(), values.data());
1649  }
1650 
1651 
1652 
1653  inline void
1654  Vector::add(const std::vector<size_type> & indices,
1655  const ::Vector<TrilinosScalar> &values)
1656  {
1657  // if we have ghost values, do not allow
1658  // writing to this vector at all.
1659  Assert(!has_ghost_elements(), ExcGhostsPresent());
1660  Assert(indices.size() == values.size(),
1661  ExcDimensionMismatch(indices.size(), values.size()));
1662 
1663  add(indices.size(), indices.data(), values.begin());
1664  }
1665 
1666 
1667 
1668  inline void
1669  Vector::add(const size_type n_elements,
1670  const size_type * indices,
1671  const TrilinosScalar *values)
1672  {
1673  // if we have ghost values, do not allow
1674  // writing to this vector at all.
1675  Assert(!has_ghost_elements(), ExcGhostsPresent());
1676 
1677  if (last_action != Add)
1678  {
1679  if (last_action == Insert)
1680  {
1681  const int ierr = vector->GlobalAssemble(Insert);
1682  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1683  }
1684  last_action = Add;
1685  }
1686 
1687  for (size_type i = 0; i < n_elements; ++i)
1688  {
1689  const size_type row = indices[i];
1690  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1691  static_cast<TrilinosWrappers::types::int_type>(row));
1692  if (local_row != -1)
1693  (*vector)[0][local_row] += values[i];
1694  else if (nonlocal_vector.get() == nullptr)
1695  {
1696  const int ierr = vector->SumIntoGlobalValues(
1697  1,
1698  reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1699  &row),
1700  &values[i]);
1701  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1702  compressed = false;
1703  }
1704  else
1705  {
1706  // use pre-allocated vector for non-local entries if it exists for
1707  // addition operation
1708  const TrilinosWrappers::types::int_type my_row =
1709  nonlocal_vector->Map().LID(
1710  static_cast<TrilinosWrappers::types::int_type>(row));
1711  Assert(my_row != -1,
1712  ExcMessage(
1713  "Attempted to write into off-processor vector entry "
1714  "that has not be specified as being writable upon "
1715  "initialization"));
1716  (*nonlocal_vector)[0][my_row] += values[i];
1717  compressed = false;
1718  }
1719  }
1720  }
1721 
1722 
1723 
1724  inline Vector::size_type
1725  Vector::size() const
1726  {
1727 # ifndef DEAL_II_WITH_64BIT_INDICES
1728  return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1729 # else
1730  return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1731 # endif
1732  }
1733 
1734 
1735 
1736  inline Vector::size_type
1737  Vector::local_size() const
1738  {
1739  return vector->Map().NumMyElements();
1740  }
1741 
1742 
1743 
1744  inline std::pair<Vector::size_type, Vector::size_type>
1745  Vector::local_range() const
1746  {
1747 # ifndef DEAL_II_WITH_64BIT_INDICES
1748  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1749  const TrilinosWrappers::types::int_type end =
1750  vector->Map().MaxMyGID() + 1;
1751 # else
1752  const TrilinosWrappers::types::int_type begin =
1753  vector->Map().MinMyGID64();
1754  const TrilinosWrappers::types::int_type end =
1755  vector->Map().MaxMyGID64() + 1;
1756 # endif
1757 
1758  Assert(
1759  end - begin == vector->Map().NumMyElements(),
1760  ExcMessage(
1761  "This function only makes sense if the elements that this "
1762  "vector stores on the current processor form a contiguous range. "
1763  "This does not appear to be the case for the current vector."));
1764 
1765  return std::make_pair(begin, end);
1766  }
1767 
1768 
1769 
1770  inline TrilinosScalar Vector::operator*(const Vector &vec) const
1771  {
1772  Assert(vector->Map().SameAs(vec.vector->Map()),
1773  ExcDifferentParallelPartitioning());
1774  Assert(!has_ghost_elements(), ExcGhostsPresent());
1775 
1776  TrilinosScalar result;
1777 
1778  const int ierr = vector->Dot(*(vec.vector), &result);
1779  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1780 
1781  return result;
1782  }
1783 
1784 
1785 
1786  inline Vector::real_type
1787  Vector::norm_sqr() const
1788  {
1789  const TrilinosScalar d = l2_norm();
1790  return d * d;
1791  }
1792 
1793 
1794 
1795  inline TrilinosScalar
1796  Vector::mean_value() const
1797  {
1798  Assert(!has_ghost_elements(), ExcGhostsPresent());
1799 
1800  TrilinosScalar mean;
1801  const int ierr = vector->MeanValue(&mean);
1802  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1803 
1804  return mean;
1805  }
1806 
1807 
1808 
1809  inline TrilinosScalar
1810  Vector::min() const
1811  {
1812  TrilinosScalar min_value;
1813  const int ierr = vector->MinValue(&min_value);
1814  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1815 
1816  return min_value;
1817  }
1818 
1819 
1820 
1821  inline TrilinosScalar
1822  Vector::max() const
1823  {
1824  TrilinosScalar max_value;
1825  const int ierr = vector->MaxValue(&max_value);
1826  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1827 
1828  return max_value;
1829  }
1830 
1831 
1832 
1833  inline Vector::real_type
1834  Vector::l1_norm() const
1835  {
1836  Assert(!has_ghost_elements(), ExcGhostsPresent());
1837 
1838  TrilinosScalar d;
1839  const int ierr = vector->Norm1(&d);
1840  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1841 
1842  return d;
1843  }
1844 
1845 
1846 
1847  inline Vector::real_type
1848  Vector::l2_norm() const
1849  {
1850  Assert(!has_ghost_elements(), ExcGhostsPresent());
1851 
1852  TrilinosScalar d;
1853  const int ierr = vector->Norm2(&d);
1854  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1855 
1856  return d;
1857  }
1858 
1859 
1860 
1861  inline Vector::real_type
1862  Vector::lp_norm(const TrilinosScalar p) const
1863  {
1864  Assert(!has_ghost_elements(), ExcGhostsPresent());
1865 
1866  TrilinosScalar norm = 0;
1867  TrilinosScalar sum = 0;
1868  const size_type n_local = local_size();
1869 
1870  // loop over all the elements because
1871  // Trilinos does not support lp norms
1872  for (size_type i = 0; i < n_local; ++i)
1873  sum += std::pow(std::fabs((*vector)[0][i]), p);
1874 
1875  norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1876 
1877  return norm;
1878  }
1879 
1880 
1881 
1882  inline Vector::real_type
1883  Vector::linfty_norm() const
1884  {
1885  // while we disallow the other
1886  // norm operations on ghosted
1887  // vectors, this particular norm
1888  // is safe to run even in the
1889  // presence of ghost elements
1890  TrilinosScalar d;
1891  const int ierr = vector->NormInf(&d);
1892  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1893 
1894  return d;
1895  }
1896 
1897 
1898 
1899  inline TrilinosScalar
1900  Vector::add_and_dot(const TrilinosScalar a,
1901  const Vector & V,
1902  const Vector & W)
1903  {
1904  this->add(a, V);
1905  return *this * W;
1906  }
1907 
1908 
1909 
1910  // inline also scalar products, vector
1911  // additions etc. since they are all
1912  // representable by a single Trilinos
1913  // call. This reduces the overhead of the
1914  // wrapper class.
1915  inline Vector &
1916  Vector::operator*=(const TrilinosScalar a)
1917  {
1918  AssertIsFinite(a);
1919 
1920  const int ierr = vector->Scale(a);
1921  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1922 
1923  return *this;
1924  }
1925 
1926 
1927 
1928  inline Vector &
1929  Vector::operator/=(const TrilinosScalar a)
1930  {
1931  AssertIsFinite(a);
1932 
1933  const TrilinosScalar factor = 1. / a;
1934 
1935  AssertIsFinite(factor);
1936 
1937  const int ierr = vector->Scale(factor);
1938  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1939 
1940  return *this;
1941  }
1942 
1943 
1944 
1945  inline Vector &
1946  Vector::operator+=(const Vector &v)
1947  {
1948  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
1949  Assert(vector->Map().SameAs(v.vector->Map()),
1950  ExcDifferentParallelPartitioning());
1951 
1952  const int ierr = vector->Update(1.0, *(v.vector), 1.0);
1953  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1954 
1955  return *this;
1956  }
1957 
1958 
1959 
1960  inline Vector &
1961  Vector::operator-=(const Vector &v)
1962  {
1963  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
1964  Assert(vector->Map().SameAs(v.vector->Map()),
1965  ExcDifferentParallelPartitioning());
1966 
1967  const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
1968  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1969 
1970  return *this;
1971  }
1972 
1973 
1974 
1975  inline void
1976  Vector::add(const TrilinosScalar s)
1977  {
1978  // if we have ghost values, do not allow
1979  // writing to this vector at all.
1980  Assert(!has_ghost_elements(), ExcGhostsPresent());
1981  AssertIsFinite(s);
1982 
1983  size_type n_local = local_size();
1984  for (size_type i = 0; i < n_local; ++i)
1985  (*vector)[0][i] += s;
1986  }
1987 
1988 
1989 
1990  inline void
1991  Vector::add(const TrilinosScalar a, const Vector &v)
1992  {
1993  // if we have ghost values, do not allow
1994  // writing to this vector at all.
1995  Assert(!has_ghost_elements(), ExcGhostsPresent());
1996  Assert(local_size() == v.local_size(),
1997  ExcDimensionMismatch(local_size(), v.local_size()));
1998 
1999  AssertIsFinite(a);
2000 
2001  const int ierr = vector->Update(a, *(v.vector), 1.);
2002  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2003  }
2004 
2005 
2006 
2007  inline void
2008  Vector::add(const TrilinosScalar a,
2009  const Vector & v,
2010  const TrilinosScalar b,
2011  const Vector & w)
2012  {
2013  // if we have ghost values, do not allow
2014  // writing to this vector at all.
2015  Assert(!has_ghost_elements(), ExcGhostsPresent());
2016  Assert(local_size() == v.local_size(),
2017  ExcDimensionMismatch(local_size(), v.local_size()));
2018  Assert(local_size() == w.local_size(),
2019  ExcDimensionMismatch(local_size(), w.local_size()));
2020 
2021  AssertIsFinite(a);
2022  AssertIsFinite(b);
2023 
2024  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2025 
2026  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2027  }
2028 
2029 
2030 
2031  inline void
2032  Vector::sadd(const TrilinosScalar s, const Vector &v)
2033  {
2034  // if we have ghost values, do not allow
2035  // writing to this vector at all.
2036  Assert(!has_ghost_elements(), ExcGhostsPresent());
2037  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2038 
2039  AssertIsFinite(s);
2040 
2041  // We assume that the vectors have the same Map
2042  // if the local size is the same and if the vectors are not ghosted
2043  if (local_size() == v.local_size() && !v.has_ghost_elements())
2044  {
2045  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2046  ExcDifferentParallelPartitioning());
2047  const int ierr = vector->Update(1., *(v.vector), s);
2048  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2049  }
2050  else
2051  {
2052  (*this) *= s;
2053  this->add(v, true);
2054  }
2055  }
2056 
2057 
2058 
2059  inline void
2060  Vector::sadd(const TrilinosScalar s,
2061  const TrilinosScalar a,
2062  const Vector & v)
2063  {
2064  // if we have ghost values, do not allow
2065  // writing to this vector at all.
2066  Assert(!has_ghost_elements(), ExcGhostsPresent());
2067  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2068  AssertIsFinite(s);
2069  AssertIsFinite(a);
2070 
2071  // We assume that the vectors have the same Map
2072  // if the local size is the same and if the vectors are not ghosted
2073  if (local_size() == v.local_size() && !v.has_ghost_elements())
2074  {
2075  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2076  ExcDifferentParallelPartitioning());
2077  const int ierr = vector->Update(a, *(v.vector), s);
2078  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2079  }
2080  else
2081  {
2082  (*this) *= s;
2083  Vector tmp = v;
2084  tmp *= a;
2085  this->add(tmp, true);
2086  }
2087  }
2088 
2089 
2090 
2091  inline void
2092  Vector::scale(const Vector &factors)
2093  {
2094  // if we have ghost values, do not allow
2095  // writing to this vector at all.
2096  Assert(!has_ghost_elements(), ExcGhostsPresent());
2097  Assert(local_size() == factors.local_size(),
2098  ExcDimensionMismatch(local_size(), factors.local_size()));
2099 
2100  const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2101  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2102  }
2103 
2104 
2105 
2106  inline void
2107  Vector::equ(const TrilinosScalar a, const Vector &v)
2108  {
2109  // if we have ghost values, do not allow
2110  // writing to this vector at all.
2111  Assert(!has_ghost_elements(), ExcGhostsPresent());
2112  AssertIsFinite(a);
2113 
2114  // If we don't have the same map, copy.
2115  if (vector->Map().SameAs(v.vector->Map()) == false)
2116  {
2117  this->sadd(0., a, v);
2118  }
2119  else
2120  {
2121  // Otherwise, just update
2122  int ierr = vector->Update(a, *v.vector, 0.0);
2123  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2124 
2125  last_action = Zero;
2126  }
2127  }
2128 
2129 
2130 
2131  inline const Epetra_MultiVector &
2132  Vector::trilinos_vector() const
2133  {
2134  return static_cast<const Epetra_MultiVector &>(*vector);
2135  }
2136 
2137 
2138 
2139  inline Epetra_FEVector &
2140  Vector::trilinos_vector()
2141  {
2142  return *vector;
2143  }
2144 
2145 
2146 
2147  inline const Epetra_Map &
2148  Vector::vector_partitioner() const
2149  {
2150  // TODO A dynamic_cast fails here. This is suspicious.
2151  return static_cast<const Epetra_Map &>(vector->Map()); // NOLINT
2152  }
2153 
2154 
2155 
2156  inline const Epetra_BlockMap &
2157  Vector::trilinos_partitioner() const
2158  {
2159  return vector->Map();
2160  }
2161 
2162 
2163 
2164  inline const MPI_Comm &
2165  Vector::get_mpi_communicator() const
2166  {
2167  static MPI_Comm comm;
2168 
2169 # ifdef DEAL_II_WITH_MPI
2170 
2171  const Epetra_MpiComm *mpi_comm =
2172  dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2173  comm = mpi_comm->Comm();
2174 
2175 # else
2176 
2177  comm = MPI_COMM_SELF;
2178 
2179 # endif
2180 
2181  return comm;
2182  }
2183 
2184  template <typename number>
2185  Vector::Vector(const IndexSet & parallel_partitioner,
2186  const ::Vector<number> &v,
2187  const MPI_Comm & communicator)
2188  {
2189  *this =
2190  Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2191  owned_elements = parallel_partitioner;
2192  }
2193 
2194 
2195 
2196  inline Vector &
2197  Vector::operator=(const TrilinosScalar s)
2198  {
2199  AssertIsFinite(s);
2200 
2201  int ierr = vector->PutScalar(s);
2202  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2203 
2204  if (nonlocal_vector.get() != nullptr)
2205  {
2206  ierr = nonlocal_vector->PutScalar(0.);
2207  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2208  }
2209 
2210  return *this;
2211  }
2212  } /* end of namespace MPI */
2213 
2214 # endif /* DOXYGEN */
2215 
2216 } /* end of namespace TrilinosWrappers */
2217 
2221 namespace internal
2222 {
2223  namespace LinearOperatorImplementation
2224  {
2225  template <typename>
2226  class ReinitHelper;
2227 
2232  template <>
2234  {
2235  public:
2236  template <typename Matrix>
2237  static void
2238  reinit_range_vector(const Matrix & matrix,
2240  bool omit_zeroing_entries)
2241  {
2242  v.reinit(matrix.locally_owned_range_indices(),
2243  matrix.get_mpi_communicator(),
2244  omit_zeroing_entries);
2245  }
2246 
2247  template <typename Matrix>
2248  static void
2249  reinit_domain_vector(const Matrix & matrix,
2251  bool omit_zeroing_entries)
2252  {
2253  v.reinit(matrix.locally_owned_domain_indices(),
2254  matrix.get_mpi_communicator(),
2255  omit_zeroing_entries);
2256  }
2257  };
2258 
2259  } // namespace LinearOperatorImplementation
2260 } /* namespace internal */
2261 
2262 
2263 
2269 template <>
2271 {};
2272 
2273 
2274 DEAL_II_NAMESPACE_CLOSE
2275 
2276 #endif // DEAL_II_WITH_TRILINOS
2277 
2278 /*---------------------------- trilinos_vector.h ---------------------------*/
2279 
2280 #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
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:1519
virtual Number operator*(const VectorSpaceVector< Number > &V) const override
virtual Vector< Number > & operator+=(const VectorSpaceVector< Number > &V) override
SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
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:594
__global__ void scale(Number *val, const Number *V_val, const size_type N)
virtual ::IndexSet locally_owned_elements() const override
Definition: la_vector.h:448
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:1407
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:439
__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
#define AssertIsFinite(number)
Definition: exceptions.h:1669