Reference documentation for deal.II version Git b1a5775265 2021-07-23 17:28:58 -0400
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_vector.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_trilinos_vector_h
17 #define dealii_trilinos_vector_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #ifdef DEAL_II_WITH_TRILINOS
23 # include <deal.II/base/index_set.h>
24 # include <deal.II/base/mpi.h>
26 # include <deal.II/base/utilities.h>
27 
28 # include <deal.II/lac/exceptions.h>
29 # include <deal.II/lac/vector.h>
32 
33 # include <Epetra_ConfigDefs.h>
34 # include <Epetra_FEVector.h>
35 # include <Epetra_LocalMap.h>
36 # include <Epetra_Map.h>
37 # include <Epetra_MpiComm.h>
38 # include <mpi.h>
39 
40 # include <memory>
41 # include <utility>
42 # include <vector>
43 
45 
46 // Forward declarations
47 # ifndef DOXYGEN
48 namespace LinearAlgebra
49 {
50  // Forward declaration
51  template <typename Number>
52  class ReadWriteVector;
53 } // namespace LinearAlgebra
54 # endif
55 
66 namespace TrilinosWrappers
67 {
68  class SparseMatrix;
69 
80  namespace internal
81  {
86 
96  class VectorReference
97  {
98  private:
103  VectorReference(MPI::Vector &vector, const size_type index);
104 
105  public:
109  VectorReference(const VectorReference &) = default;
110 
122  const VectorReference &
123  operator=(const VectorReference &r) const;
124 
128  VectorReference &
129  operator=(const VectorReference &r);
130 
134  const VectorReference &
135  operator=(const TrilinosScalar &s) const;
136 
140  const VectorReference &
141  operator+=(const TrilinosScalar &s) const;
142 
146  const VectorReference &
147  operator-=(const TrilinosScalar &s) const;
148 
152  const VectorReference &
153  operator*=(const TrilinosScalar &s) const;
154 
158  const VectorReference &
159  operator/=(const TrilinosScalar &s) const;
160 
165  operator TrilinosScalar() const;
166 
171  int,
172  << "An error with error number " << arg1
173  << " occurred while calling a Trilinos function");
174 
175  private:
179  MPI::Vector &vector;
180 
184  const size_type index;
185 
186  // Make the vector class a friend, so that it can create objects of the
187  // present type.
189  };
190  } // namespace internal
195 # ifndef DEAL_II_WITH_64BIT_INDICES
196  // define a helper function that queries the global ID of local ID of
197  // an Epetra_BlockMap object by calling either the 32- or 64-bit
198  // function necessary.
199  inline int
200  gid(const Epetra_BlockMap &map, int i)
201  {
202  return map.GID(i);
203  }
204 # else
205  // define a helper function that queries the global ID of local ID of
206  // an Epetra_BlockMap object by calling either the 32- or 64-bit
207  // function necessary.
208  inline long long int
209  gid(const Epetra_BlockMap &map, int i)
210  {
211  return map.GID64(i);
212  }
213 # endif
214 
220  namespace MPI
221  {
222  class BlockVector;
223 
399  class Vector : public Subscriptor
400  {
401  public:
410  using iterator = value_type *;
411  using const_iterator = const value_type *;
412  using reference = internal::VectorReference;
413  using const_reference = const internal::VectorReference;
414 
424  Vector();
425 
429  Vector(const Vector &v);
430 
449  explicit Vector(const IndexSet &parallel_partitioning,
450  const MPI_Comm &communicator = MPI_COMM_WORLD);
451 
463  Vector(const IndexSet &local,
464  const IndexSet &ghost,
465  const MPI_Comm &communicator = MPI_COMM_WORLD);
466 
481  Vector(const IndexSet &parallel_partitioning,
482  const Vector & v,
483  const MPI_Comm &communicator = MPI_COMM_WORLD);
484 
497  template <typename Number>
498  Vector(const IndexSet & parallel_partitioning,
499  const ::Vector<Number> &v,
500  const MPI_Comm & communicator = MPI_COMM_WORLD);
501 
506  Vector(Vector &&v) noexcept;
507 
511  ~Vector() override = default;
512 
517  void
518  clear();
519 
543  void
544  reinit(const Vector &v,
545  const bool omit_zeroing_entries = false,
546  const bool allow_different_maps = false);
547 
570  void
571  reinit(const IndexSet &parallel_partitioning,
572  const MPI_Comm &communicator = MPI_COMM_WORLD,
573  const bool omit_zeroing_entries = false);
574 
600  void
601  reinit(const IndexSet &locally_owned_entries,
602  const IndexSet &ghost_entries,
603  const MPI_Comm &communicator = MPI_COMM_WORLD,
604  const bool vector_writable = false);
605 
609  void
610  reinit(const BlockVector &v, const bool import_data = false);
611 
628  void
629  compress(::VectorOperation::values operation);
630 
643  Vector &
644  operator=(const TrilinosScalar s);
645 
651  Vector &
652  operator=(const Vector &v);
653 
658  Vector &
659  operator=(Vector &&v) noexcept;
660 
668  template <typename Number>
669  Vector &
670  operator=(const ::Vector<Number> &v);
671 
689  void
690  import_nonlocal_data_for_fe(
692  const Vector & vector);
693 
700  void
701  import(const LinearAlgebra::ReadWriteVector<double> &rwv,
702  const VectorOperation::values operation);
703 
704 
710  bool
711  operator==(const Vector &v) const;
712 
718  bool
719  operator!=(const Vector &v) const;
720 
724  size_type
725  size() const;
726 
741  size_type
742  local_size() const;
743 
748  size_type
749  locally_owned_size() const;
750 
772  std::pair<size_type, size_type>
773  local_range() const;
774 
782  bool
783  in_local_range(const size_type index) const;
784 
798  IndexSet
799  locally_owned_elements() const;
800 
808  bool
809  has_ghost_elements() const;
810 
817  void
818  update_ghost_values() const;
819 
824  TrilinosScalar operator*(const Vector &vec) const;
825 
829  real_type
830  norm_sqr() const;
831 
836  mean_value() const;
837 
842  min() const;
843 
848  max() const;
849 
853  real_type
854  l1_norm() const;
855 
860  real_type
861  l2_norm() const;
862 
867  real_type
868  lp_norm(const TrilinosScalar p) const;
869 
873  real_type
874  linfty_norm() const;
875 
896  add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
897 
903  bool
904  all_zero() const;
905 
911  bool
912  is_non_negative() const;
914 
915 
920 
928  reference
929  operator()(const size_type index);
930 
939  operator()(const size_type index) const;
940 
946  reference operator[](const size_type index);
947 
953  TrilinosScalar operator[](const size_type index) const;
954 
970  void
971  extract_subvector_to(const std::vector<size_type> &indices,
972  std::vector<TrilinosScalar> & values) const;
973 
1001  template <typename ForwardIterator, typename OutputIterator>
1002  void
1003  extract_subvector_to(ForwardIterator indices_begin,
1004  const ForwardIterator indices_end,
1005  OutputIterator values_begin) const;
1006 
1017  iterator
1018  begin();
1019 
1025  begin() const;
1026 
1031  iterator
1032  end();
1033 
1039  end() const;
1040 
1042 
1043 
1048 
1055  void
1056  set(const std::vector<size_type> & indices,
1057  const std::vector<TrilinosScalar> &values);
1058 
1063  void
1064  set(const std::vector<size_type> & indices,
1065  const ::Vector<TrilinosScalar> &values);
1066 
1072  void
1073  set(const size_type n_elements,
1074  const size_type * indices,
1075  const TrilinosScalar *values);
1076 
1081  void
1082  add(const std::vector<size_type> & indices,
1083  const std::vector<TrilinosScalar> &values);
1084 
1089  void
1090  add(const std::vector<size_type> & indices,
1091  const ::Vector<TrilinosScalar> &values);
1092 
1098  void
1099  add(const size_type n_elements,
1100  const size_type * indices,
1101  const TrilinosScalar *values);
1102 
1106  Vector &
1107  operator*=(const TrilinosScalar factor);
1108 
1112  Vector &
1113  operator/=(const TrilinosScalar factor);
1114 
1118  Vector &
1119  operator+=(const Vector &V);
1120 
1124  Vector &
1125  operator-=(const Vector &V);
1126 
1131  void
1132  add(const TrilinosScalar s);
1133 
1146  void
1147  add(const Vector &V, const bool allow_different_maps = false);
1148 
1152  void
1153  add(const TrilinosScalar a, const Vector &V);
1154 
1158  void
1159  add(const TrilinosScalar a,
1160  const Vector & V,
1161  const TrilinosScalar b,
1162  const Vector & W);
1163 
1168  void
1169  sadd(const TrilinosScalar s, const Vector &V);
1170 
1174  void
1175  sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1176 
1182  void
1183  scale(const Vector &scaling_factors);
1184 
1188  void
1189  equ(const TrilinosScalar a, const Vector &V);
1191 
1196 
1201  const Epetra_MultiVector &
1202  trilinos_vector() const;
1203 
1208  Epetra_FEVector &
1209  trilinos_vector();
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 are trying 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  << arg4 << "] but Trilinos vectors need not store contiguous "
1289  << "ranges on each processor, and not every element in "
1290  << "this range may in fact be stored locally."
1291  << "\n\n"
1292  << "A common source for this kind of problem is that you "
1293  << "are passing a 'fully distributed' vector into a function "
1294  << "that needs read access to vector elements that correspond "
1295  << "to degrees of freedom on ghost cells (or at least to "
1296  << "'locally active' degrees of freedom that are not also "
1297  << "'locally owned'). You need to pass a vector that has these "
1298  << "elements as ghost entries.");
1299 
1300  private:
1312  Epetra_CombineMode last_action;
1313 
1319 
1325 
1331  std::unique_ptr<Epetra_FEVector> vector;
1332 
1338  std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1339 
1344 
1345  // Make the reference class a friend.
1346  friend class internal::VectorReference;
1347  };
1348 
1349 
1350 
1351  // ------------------- inline and template functions --------------
1352 
1353 
1361  inline void
1363  {
1364  u.swap(v);
1365  }
1366  } // namespace MPI
1367 
1368 # ifndef DOXYGEN
1369 
1370  namespace internal
1371  {
1372  inline VectorReference::VectorReference(MPI::Vector & vector,
1373  const size_type index)
1374  : vector(vector)
1375  , index(index)
1376  {}
1377 
1378 
1379  inline const VectorReference &
1380  VectorReference::operator=(const VectorReference &r) const
1381  {
1382  // as explained in the class
1383  // documentation, this is not the copy
1384  // operator. so simply pass on to the
1385  // "correct" assignment operator
1386  *this = static_cast<TrilinosScalar>(r);
1387 
1388  return *this;
1389  }
1390 
1391 
1392 
1393  inline VectorReference &
1394  VectorReference::operator=(const VectorReference &r)
1395  {
1396  // as above
1397  *this = static_cast<TrilinosScalar>(r);
1398 
1399  return *this;
1400  }
1401 
1402 
1403  inline const VectorReference &
1404  VectorReference::operator=(const TrilinosScalar &value) const
1405  {
1406  vector.set(1, &index, &value);
1407  return *this;
1408  }
1409 
1410 
1411 
1412  inline const VectorReference &
1413  VectorReference::operator+=(const TrilinosScalar &value) const
1414  {
1415  vector.add(1, &index, &value);
1416  return *this;
1417  }
1418 
1419 
1420 
1421  inline const VectorReference &
1422  VectorReference::operator-=(const TrilinosScalar &value) const
1423  {
1424  TrilinosScalar new_value = -value;
1425  vector.add(1, &index, &new_value);
1426  return *this;
1427  }
1428 
1429 
1430 
1431  inline const VectorReference &
1432  VectorReference::operator*=(const TrilinosScalar &value) const
1433  {
1434  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1435  vector.set(1, &index, &new_value);
1436  return *this;
1437  }
1438 
1439 
1440 
1441  inline const VectorReference &
1442  VectorReference::operator/=(const TrilinosScalar &value) const
1443  {
1444  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1445  vector.set(1, &index, &new_value);
1446  return *this;
1447  }
1448  } // namespace internal
1449 
1450  namespace MPI
1451  {
1452  inline bool
1453  Vector::in_local_range(const size_type index) const
1454  {
1455  std::pair<size_type, size_type> range = local_range();
1456 
1457  return ((index >= range.first) && (index < range.second));
1458  }
1459 
1460 
1461 
1462  inline IndexSet
1464  {
1465  Assert(owned_elements.size() == size(),
1466  ExcMessage(
1467  "The locally owned elements have not been properly initialized!"
1468  " This happens for example if this object has been initialized"
1469  " with exactly one overlapping IndexSet."));
1470  return owned_elements;
1471  }
1472 
1473 
1474 
1475  inline bool
1477  {
1478  return has_ghosts;
1479  }
1480 
1481 
1482 
1483  inline void
1485  {}
1486 
1487 
1488 
1489  inline internal::VectorReference
1490  Vector::operator()(const size_type index)
1491  {
1492  return internal::VectorReference(*this, index);
1493  }
1494 
1495 
1496 
1497  inline internal::VectorReference Vector::operator[](const size_type index)
1498  {
1499  return operator()(index);
1500  }
1501 
1502 
1503 
1504  inline TrilinosScalar Vector::operator[](const size_type index) const
1505  {
1506  return operator()(index);
1507  }
1508 
1509 
1510 
1511  inline void
1512  Vector::extract_subvector_to(const std::vector<size_type> &indices,
1513  std::vector<TrilinosScalar> & values) const
1514  {
1515  for (size_type i = 0; i < indices.size(); ++i)
1516  values[i] = operator()(indices[i]);
1517  }
1518 
1519 
1520 
1521  template <typename ForwardIterator, typename OutputIterator>
1522  inline void
1523  Vector::extract_subvector_to(ForwardIterator indices_begin,
1524  const ForwardIterator indices_end,
1525  OutputIterator values_begin) const
1526  {
1527  while (indices_begin != indices_end)
1528  {
1529  *values_begin = operator()(*indices_begin);
1530  indices_begin++;
1531  values_begin++;
1532  }
1533  }
1534 
1535 
1536 
1537  inline Vector::iterator
1538  Vector::begin()
1539  {
1540  return (*vector)[0];
1541  }
1542 
1543 
1544 
1545  inline Vector::iterator
1546  Vector::end()
1547  {
1548  return (*vector)[0] + locally_owned_size();
1549  }
1550 
1551 
1552 
1553  inline Vector::const_iterator
1554  Vector::begin() const
1555  {
1556  return (*vector)[0];
1557  }
1558 
1559 
1560 
1561  inline Vector::const_iterator
1562  Vector::end() const
1563  {
1564  return (*vector)[0] + locally_owned_size();
1565  }
1566 
1567 
1568 
1569  inline void
1570  Vector::set(const std::vector<size_type> & indices,
1571  const std::vector<TrilinosScalar> &values)
1572  {
1573  // if we have ghost values, do not allow
1574  // writing to this vector at all.
1575  Assert(!has_ghost_elements(), ExcGhostsPresent());
1576 
1577  AssertDimension(indices.size(), values.size());
1578 
1579  set(indices.size(), indices.data(), values.data());
1580  }
1581 
1582 
1583 
1584  inline void
1585  Vector::set(const std::vector<size_type> & indices,
1586  const ::Vector<TrilinosScalar> &values)
1587  {
1588  // if we have ghost values, do not allow
1589  // writing to this vector at all.
1590  Assert(!has_ghost_elements(), ExcGhostsPresent());
1591 
1592  AssertDimension(indices.size(), values.size());
1593 
1594  set(indices.size(), indices.data(), values.begin());
1595  }
1596 
1597 
1598 
1599  inline void
1600  Vector::set(const size_type n_elements,
1601  const size_type * indices,
1602  const TrilinosScalar *values)
1603  {
1604  // if we have ghost values, do not allow
1605  // writing to this vector at all.
1606  Assert(!has_ghost_elements(), ExcGhostsPresent());
1607 
1608  if (last_action == Add)
1609  {
1610  const int ierr = vector->GlobalAssemble(Add);
1611  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1612  }
1613 
1614  if (last_action != Insert)
1615  last_action = Insert;
1616 
1617  for (size_type i = 0; i < n_elements; ++i)
1618  {
1619  const TrilinosWrappers::types::int_type row = indices[i];
1620  const TrilinosWrappers::types::int_type local_row =
1621  vector->Map().LID(row);
1622  if (local_row != -1)
1623  (*vector)[0][local_row] = values[i];
1624  else
1625  {
1626  const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1627  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1628  compressed = false;
1629  }
1630  // in set operation, do not use the pre-allocated vector for nonlocal
1631  // entries even if it exists. This is to ensure that we really only
1632  // set the elements touched by the set() method and not all contained
1633  // in the nonlocal entries vector (there is no way to distinguish them
1634  // on the receiving processor)
1635  }
1636  }
1637 
1638 
1639 
1640  inline void
1641  Vector::add(const std::vector<size_type> & indices,
1642  const std::vector<TrilinosScalar> &values)
1643  {
1644  // if we have ghost values, do not allow
1645  // writing to this vector at all.
1646  Assert(!has_ghost_elements(), ExcGhostsPresent());
1647  AssertDimension(indices.size(), values.size());
1648 
1649  add(indices.size(), indices.data(), values.data());
1650  }
1651 
1652 
1653 
1654  inline void
1655  Vector::add(const std::vector<size_type> & indices,
1656  const ::Vector<TrilinosScalar> &values)
1657  {
1658  // if we have ghost values, do not allow
1659  // writing to this vector at all.
1660  Assert(!has_ghost_elements(), ExcGhostsPresent());
1661  AssertDimension(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 Vector::size_type
1746  {
1747  return owned_elements.n_elements();
1748  }
1749 
1750 
1751 
1752  inline std::pair<Vector::size_type, Vector::size_type>
1753  Vector::local_range() const
1754  {
1755 # ifndef DEAL_II_WITH_64BIT_INDICES
1756  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1758  vector->Map().MaxMyGID() + 1;
1759 # else
1760  const TrilinosWrappers::types::int_type begin =
1761  vector->Map().MinMyGID64();
1763  vector->Map().MaxMyGID64() + 1;
1764 # endif
1765 
1766  Assert(
1767  end - begin == vector->Map().NumMyElements(),
1768  ExcMessage(
1769  "This function only makes sense if the elements that this "
1770  "vector stores on the current processor form a contiguous range. "
1771  "This does not appear to be the case for the current vector."));
1772 
1773  return std::make_pair(begin, end);
1774  }
1775 
1776 
1777 
1778  inline TrilinosScalar Vector::operator*(const Vector &vec) const
1779  {
1780  Assert(vector->Map().SameAs(vec.vector->Map()),
1781  ExcDifferentParallelPartitioning());
1782  Assert(!has_ghost_elements(), ExcGhostsPresent());
1783 
1784  TrilinosScalar result;
1785 
1786  const int ierr = vector->Dot(*(vec.vector), &result);
1787  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1788 
1789  return result;
1790  }
1791 
1792 
1793 
1794  inline Vector::real_type
1795  Vector::norm_sqr() const
1796  {
1797  const TrilinosScalar d = l2_norm();
1798  return d * d;
1799  }
1800 
1801 
1802 
1803  inline TrilinosScalar
1804  Vector::mean_value() const
1805  {
1806  Assert(!has_ghost_elements(), ExcGhostsPresent());
1807 
1809  const int ierr = vector->MeanValue(&mean);
1810  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1811 
1812  return mean;
1813  }
1814 
1815 
1816 
1817  inline TrilinosScalar
1818  Vector::min() const
1819  {
1820  TrilinosScalar min_value;
1821  const int ierr = vector->MinValue(&min_value);
1822  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1823 
1824  return min_value;
1825  }
1826 
1827 
1828 
1829  inline TrilinosScalar
1830  Vector::max() const
1831  {
1832  TrilinosScalar max_value;
1833  const int ierr = vector->MaxValue(&max_value);
1834  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1835 
1836  return max_value;
1837  }
1838 
1839 
1840 
1841  inline Vector::real_type
1842  Vector::l1_norm() const
1843  {
1844  Assert(!has_ghost_elements(), ExcGhostsPresent());
1845 
1846  TrilinosScalar d;
1847  const int ierr = vector->Norm1(&d);
1848  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1849 
1850  return d;
1851  }
1852 
1853 
1854 
1855  inline Vector::real_type
1856  Vector::l2_norm() const
1857  {
1858  Assert(!has_ghost_elements(), ExcGhostsPresent());
1859 
1860  TrilinosScalar d;
1861  const int ierr = vector->Norm2(&d);
1862  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1863 
1864  return d;
1865  }
1866 
1867 
1868 
1869  inline Vector::real_type
1870  Vector::lp_norm(const TrilinosScalar p) const
1871  {
1872  Assert(!has_ghost_elements(), ExcGhostsPresent());
1873 
1874  TrilinosScalar norm = 0;
1875  TrilinosScalar sum = 0;
1876  const size_type n_local = locally_owned_size();
1877 
1878  // loop over all the elements because
1879  // Trilinos does not support lp norms
1880  for (size_type i = 0; i < n_local; ++i)
1881  sum += std::pow(std::fabs((*vector)[0][i]), p);
1882 
1883  norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1884 
1885  return norm;
1886  }
1887 
1888 
1889 
1890  inline Vector::real_type
1891  Vector::linfty_norm() const
1892  {
1893  // while we disallow the other
1894  // norm operations on ghosted
1895  // vectors, this particular norm
1896  // is safe to run even in the
1897  // presence of ghost elements
1898  TrilinosScalar d;
1899  const int ierr = vector->NormInf(&d);
1900  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1901 
1902  return d;
1903  }
1904 
1905 
1906 
1907  inline TrilinosScalar
1909  const Vector & V,
1910  const Vector & W)
1911  {
1912  this->add(a, V);
1913  return *this * W;
1914  }
1915 
1916 
1917 
1918  // inline also scalar products, vector
1919  // additions etc. since they are all
1920  // representable by a single Trilinos
1921  // call. This reduces the overhead of the
1922  // wrapper class.
1923  inline Vector &
1925  {
1926  AssertIsFinite(a);
1927 
1928  const int ierr = vector->Scale(a);
1929  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1930 
1931  return *this;
1932  }
1933 
1934 
1935 
1936  inline Vector &
1938  {
1939  AssertIsFinite(a);
1940 
1941  const TrilinosScalar factor = 1. / a;
1942 
1943  AssertIsFinite(factor);
1944 
1945  const int ierr = vector->Scale(factor);
1946  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1947 
1948  return *this;
1949  }
1950 
1951 
1952 
1953  inline Vector &
1954  Vector::operator+=(const Vector &v)
1955  {
1956  AssertDimension(size(), v.size());
1957  Assert(vector->Map().SameAs(v.vector->Map()),
1958  ExcDifferentParallelPartitioning());
1959 
1960  const int ierr = vector->Update(1.0, *(v.vector), 1.0);
1961  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1962 
1963  return *this;
1964  }
1965 
1966 
1967 
1968  inline Vector &
1969  Vector::operator-=(const Vector &v)
1970  {
1971  AssertDimension(size(), v.size());
1972  Assert(vector->Map().SameAs(v.vector->Map()),
1973  ExcDifferentParallelPartitioning());
1974 
1975  const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
1976  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1977 
1978  return *this;
1979  }
1980 
1981 
1982 
1983  inline void
1984  Vector::add(const TrilinosScalar s)
1985  {
1986  // if we have ghost values, do not allow
1987  // writing to this vector at all.
1988  Assert(!has_ghost_elements(), ExcGhostsPresent());
1989  AssertIsFinite(s);
1990 
1991  size_type n_local = locally_owned_size();
1992  for (size_type i = 0; i < n_local; ++i)
1993  (*vector)[0][i] += s;
1994  }
1995 
1996 
1997 
1998  inline void
1999  Vector::add(const TrilinosScalar a, const Vector &v)
2000  {
2001  // if we have ghost values, do not allow
2002  // writing to this vector at all.
2003  Assert(!has_ghost_elements(), ExcGhostsPresent());
2004  AssertDimension(locally_owned_size(), v.locally_owned_size());
2005 
2006  AssertIsFinite(a);
2007 
2008  const int ierr = vector->Update(a, *(v.vector), 1.);
2009  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2010  }
2011 
2012 
2013 
2014  inline void
2015  Vector::add(const TrilinosScalar a,
2016  const Vector & v,
2017  const TrilinosScalar b,
2018  const Vector & w)
2019  {
2020  // if we have ghost values, do not allow
2021  // writing to this vector at all.
2022  Assert(!has_ghost_elements(), ExcGhostsPresent());
2023  AssertDimension(locally_owned_size(), v.locally_owned_size());
2024  AssertDimension(locally_owned_size(), w.locally_owned_size());
2025 
2026  AssertIsFinite(a);
2027  AssertIsFinite(b);
2028 
2029  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2030 
2031  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2032  }
2033 
2034 
2035 
2036  inline void
2037  Vector::sadd(const TrilinosScalar s, const Vector &v)
2038  {
2039  // if we have ghost values, do not allow
2040  // writing to this vector at all.
2041  Assert(!has_ghost_elements(), ExcGhostsPresent());
2042  AssertDimension(size(), v.size());
2043 
2044  AssertIsFinite(s);
2045 
2046  // We assume that the vectors have the same Map
2047  // if the local size is the same and if the vectors are not ghosted
2048  if (locally_owned_size() == v.locally_owned_size() &&
2049  !v.has_ghost_elements())
2050  {
2051  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2052  ExcDifferentParallelPartitioning());
2053  const int ierr = vector->Update(1., *(v.vector), s);
2054  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2055  }
2056  else
2057  {
2058  (*this) *= s;
2059  this->add(v, true);
2060  }
2061  }
2062 
2063 
2064 
2065  inline void
2066  Vector::sadd(const TrilinosScalar s,
2067  const TrilinosScalar a,
2068  const Vector & v)
2069  {
2070  // if we have ghost values, do not allow
2071  // writing to this vector at all.
2072  Assert(!has_ghost_elements(), ExcGhostsPresent());
2073  AssertDimension(size(), v.size());
2074  AssertIsFinite(s);
2075  AssertIsFinite(a);
2076 
2077  // We assume that the vectors have the same Map
2078  // if the local size is the same and if the vectors are not ghosted
2079  if (locally_owned_size() == v.locally_owned_size() &&
2080  !v.has_ghost_elements())
2081  {
2082  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2083  ExcDifferentParallelPartitioning());
2084  const int ierr = vector->Update(a, *(v.vector), s);
2085  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2086  }
2087  else
2088  {
2089  (*this) *= s;
2090  Vector tmp = v;
2091  tmp *= a;
2092  this->add(tmp, true);
2093  }
2094  }
2095 
2096 
2097 
2098  inline void
2099  Vector::scale(const Vector &factors)
2100  {
2101  // if we have ghost values, do not allow
2102  // writing to this vector at all.
2103  Assert(!has_ghost_elements(), ExcGhostsPresent());
2104  AssertDimension(locally_owned_size(), factors.locally_owned_size());
2105 
2106  const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2107  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2108  }
2109 
2110 
2111 
2112  inline void
2113  Vector::equ(const TrilinosScalar a, const Vector &v)
2114  {
2115  // if we have ghost values, do not allow
2116  // writing to this vector at all.
2117  Assert(!has_ghost_elements(), ExcGhostsPresent());
2118  AssertIsFinite(a);
2119 
2120  // If we don't have the same map, copy.
2121  if (vector->Map().SameAs(v.vector->Map()) == false)
2122  {
2123  this->sadd(0., a, v);
2124  }
2125  else
2126  {
2127  // Otherwise, just update
2128  int ierr = vector->Update(a, *v.vector, 0.0);
2129  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2130 
2131  last_action = Zero;
2132  }
2133  }
2134 
2135 
2136 
2137  inline const Epetra_MultiVector &
2138  Vector::trilinos_vector() const
2139  {
2140  return static_cast<const Epetra_MultiVector &>(*vector);
2141  }
2142 
2143 
2144 
2145  inline Epetra_FEVector &
2146  Vector::trilinos_vector()
2147  {
2148  return *vector;
2149  }
2150 
2151 
2152 
2153  inline const Epetra_BlockMap &
2154  Vector::trilinos_partitioner() const
2155  {
2156  return vector->Map();
2157  }
2158 
2159 
2160 
2161  inline const MPI_Comm &
2162  Vector::get_mpi_communicator() const
2163  {
2164  static MPI_Comm comm;
2165 
2166  const Epetra_MpiComm *mpi_comm =
2167  dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2168  comm = mpi_comm->Comm();
2169 
2170  return comm;
2171  }
2172 
2173  template <typename number>
2174  Vector::Vector(const IndexSet & parallel_partitioner,
2175  const ::Vector<number> &v,
2176  const MPI_Comm & communicator)
2177  {
2178  *this =
2179  Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2180  owned_elements = parallel_partitioner;
2181  }
2182 
2183 
2184 
2185  inline Vector &
2187  {
2188  AssertIsFinite(s);
2189 
2190  int ierr = vector->PutScalar(s);
2191  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2192 
2193  if (nonlocal_vector.get() != nullptr)
2194  {
2195  ierr = nonlocal_vector->PutScalar(0.);
2196  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2197  }
2198 
2199  return *this;
2200  }
2201  } /* end of namespace MPI */
2202 
2203 # endif /* DOXYGEN */
2204 
2205 } /* end of namespace TrilinosWrappers */
2206 
2210 namespace internal
2211 {
2212  namespace LinearOperatorImplementation
2213  {
2214  template <typename>
2215  class ReinitHelper;
2216 
2221  template <>
2223  {
2224  public:
2225  template <typename Matrix>
2226  static void
2227  reinit_range_vector(const Matrix & matrix,
2229  bool omit_zeroing_entries)
2230  {
2231  v.reinit(matrix.locally_owned_range_indices(),
2232  matrix.get_mpi_communicator(),
2233  omit_zeroing_entries);
2234  }
2235 
2236  template <typename Matrix>
2237  static void
2240  bool omit_zeroing_entries)
2241  {
2242  v.reinit(matrix.locally_owned_domain_indices(),
2243  matrix.get_mpi_communicator(),
2244  omit_zeroing_entries);
2245  }
2246  };
2247 
2248  } // namespace LinearOperatorImplementation
2249 } /* namespace internal */
2250 
2251 
2252 
2256 template <>
2258 {};
2259 
2260 
2262 
2263 #endif // DEAL_II_WITH_TRILINOS
2264 
2265 /*---------------------------- trilinos_vector.h ---------------------------*/
2266 
2267 #endif // dealii_trilinos_vector_h
std::unique_ptr< Epetra_MultiVector > nonlocal_vector
static ::ExceptionBase & ExcTrilinosError(int arg1)
void update_ghost_values() const
int gid(const Epetra_BlockMap &map, int i)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1657
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:2040
void reinit(const Vector &v, const bool omit_zeroing_entries=false, const bool allow_different_maps=false)
virtual void add(const Number a) override
void add(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
__global__ void add_and_dot(Number *res, Number *v1, const Number *v2, const Number *v3, const Number a, const size_type N)
const internal::VectorReference const_reference
static const char V
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
Vector< Number > & operator=(const Vector< Number > &in_vector)
virtual Vector< Number > & operator*=(const Number factor) override
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
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:603
virtual ::IndexSet locally_owned_elements() const override
Definition: la_vector.h:462
static void reinit_domain_vector(const Matrix &matrix, TrilinosWrappers::MPI::Vector &v, bool omit_zeroing_entries)
virtual Vector< Number > & operator-=(const VectorSpaceVector< Number > &V) override
virtual void sadd(const Number s, const Number a, const VectorSpaceVector< Number > &V) override
static ::ExceptionBase & ExcMessage(std::string arg1)
std::string compress(const std::string &input)
Definition: utilities.cc:392
bool in_local_range(const size_type global_index) const
void set(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1465
std::unique_ptr< Epetra_FEVector > vector
real_type lp_norm(const real_type p) const
size_type locally_owned_size() const
Number operator()(const size_type global_index) const
void swap(Vector &u, Vector &v)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2939
#define DeclException0(Exception0)
Definition: exceptions.h:470
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
virtual Vector< Number > & operator/=(const Number factor) override
virtual void scale(const VectorSpaceVector< Number > &scaling_factors) override
VectorType::value_type * end(VectorType &V)
double TrilinosScalar
Definition: types.h:163
virtual size_type size() const override
Definition: la_vector.h:453
Expression fabs(const Expression &x)
__global__ void equ(Number *val, const Number a, const Number *V_val, const size_type N)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SparseMatrix< double > SparseMatrix
__global__ void sadd(const Number s, Number *val, const Number a, const Number *V_val, const size_type N)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
size_type locally_owned_size() const
static ::ExceptionBase & ExcGhostsPresent()
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int global_dof_index
Definition: types.h:76
__global__ void set(Number *val, const Number s, const size_type N)
virtual VectorSpaceVector< Number >::real_type l1_norm() const override
virtual VectorSpaceVector< Number >::real_type l2_norm() const override
*braid_SplitCommworld & comm
internal::VectorReference reference
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
VectorType::value_type * begin(VectorType &V)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:584
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2913
virtual value_type mean_value() const override
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< Number2 > &values) const
real_type norm_sqr() const
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
virtual Number add_and_dot(const Number a, const VectorSpaceVector< Number > &V, const VectorSpaceVector< Number > &W) override
#define DEAL_II_DEPRECATED
Definition: config.h:160
bool has_ghost_elements() const
#define AssertIsFinite(number)
Definition: exceptions.h:1748
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)