Reference documentation for deal.II version Git 4abc4a1666 2020-07-04 19:58:34 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_vector.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_trilinos_vector_h
17 #define dealii_trilinos_vector_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #ifdef DEAL_II_WITH_TRILINOS
23 # include <deal.II/base/index_set.h>
24 # include <deal.II/base/mpi.h>
26 # include <deal.II/base/utilities.h>
27 
28 # include <deal.II/lac/exceptions.h>
29 # include <deal.II/lac/vector.h>
32 
33 # include <Epetra_ConfigDefs.h>
34 
35 # include <memory>
36 # include <utility>
37 # include <vector>
38 # ifdef DEAL_II_WITH_MPI // only if MPI is installed
39 # include <Epetra_MpiComm.h>
40 # include <mpi.h>
41 # else
42 # include <Epetra_SerialComm.h>
43 # endif
44 # include <Epetra_FEVector.h>
45 # include <Epetra_LocalMap.h>
46 # include <Epetra_Map.h>
47 
49 
50 // Forward declarations
51 # ifndef DOXYGEN
52 namespace LinearAlgebra
53 {
54  // Forward declaration
55  template <typename Number>
56  class ReadWriteVector;
57 } // namespace LinearAlgebra
58 # endif
59 
70 namespace TrilinosWrappers
71 {
72  class SparseMatrix;
73 
84  namespace internal
85  {
90 
100  class VectorReference
101  {
102  private:
107  VectorReference(MPI::Vector &vector, const size_type index);
108 
109  public:
113  VectorReference(const VectorReference &) = default;
114 
126  const VectorReference &
127  operator=(const VectorReference &r) const;
128 
132  VectorReference &
133  operator=(const VectorReference &r);
134 
138  const VectorReference &
139  operator=(const TrilinosScalar &s) const;
140 
144  const VectorReference &
145  operator+=(const TrilinosScalar &s) const;
146 
150  const VectorReference &
151  operator-=(const TrilinosScalar &s) const;
152 
156  const VectorReference &
157  operator*=(const TrilinosScalar &s) const;
158 
162  const VectorReference &
163  operator/=(const TrilinosScalar &s) const;
164 
169  operator TrilinosScalar() const;
170 
175  int,
176  << "An error with error number " << arg1
177  << " occurred while calling a Trilinos function");
178 
179  private:
183  MPI::Vector &vector;
184 
188  const size_type index;
189 
190  // Make the vector class a friend, so that it can create objects of the
191  // present type.
193  };
194  } // namespace internal
199 # ifndef DEAL_II_WITH_64BIT_INDICES
200  // define a helper function that queries the global ID of local ID of
201  // an Epetra_BlockMap object by calling either the 32- or 64-bit
202  // function necessary.
203  inline int
204  gid(const Epetra_BlockMap &map, int i)
205  {
206  return map.GID(i);
207  }
208 # else
209  // define a helper function that queries the global ID of local ID of
210  // an Epetra_BlockMap object by calling either the 32- or 64-bit
211  // function necessary.
212  inline long long int
213  gid(const Epetra_BlockMap &map, int i)
214  {
215  return map.GID64(i);
216  }
217 # endif
218 
224  namespace MPI
225  {
226  class BlockVector;
227 
403  class Vector : public Subscriptor
404  {
405  public:
414  using iterator = value_type *;
415  using const_iterator = const value_type *;
416  using reference = internal::VectorReference;
417  using const_reference = const internal::VectorReference;
418 
428  Vector();
429 
433  Vector(const Vector &v);
434 
453  explicit Vector(const IndexSet &parallel_partitioning,
454  const MPI_Comm &communicator = MPI_COMM_WORLD);
455 
467  Vector(const IndexSet &local,
468  const IndexSet &ghost,
469  const MPI_Comm &communicator = MPI_COMM_WORLD);
470 
485  Vector(const IndexSet &parallel_partitioning,
486  const Vector & v,
487  const MPI_Comm &communicator = MPI_COMM_WORLD);
488 
501  template <typename Number>
502  Vector(const IndexSet & parallel_partitioning,
503  const ::Vector<Number> &v,
504  const MPI_Comm & communicator = MPI_COMM_WORLD);
505 
510  Vector(Vector &&v) noexcept;
511 
515  ~Vector() override = default;
516 
521  void
522  clear();
523 
547  void
548  reinit(const Vector &v,
549  const bool omit_zeroing_entries = false,
550  const bool allow_different_maps = false);
551 
574  void
575  reinit(const IndexSet &parallel_partitioning,
576  const MPI_Comm &communicator = MPI_COMM_WORLD,
577  const bool omit_zeroing_entries = false);
578 
604  void
605  reinit(const IndexSet &locally_owned_entries,
606  const IndexSet &ghost_entries,
607  const MPI_Comm &communicator = MPI_COMM_WORLD,
608  const bool vector_writable = false);
609 
613  void
614  reinit(const BlockVector &v, const bool import_data = false);
615 
632  void
633  compress(::VectorOperation::values operation);
634 
647  Vector &
648  operator=(const TrilinosScalar s);
649 
655  Vector &
656  operator=(const Vector &v);
657 
662  Vector &
663  operator=(Vector &&v) noexcept;
664 
672  template <typename Number>
673  Vector &
674  operator=(const ::Vector<Number> &v);
675 
693  void
694  import_nonlocal_data_for_fe(
696  const Vector & vector);
697 
704  void
705  import(const LinearAlgebra::ReadWriteVector<double> &rwv,
706  const VectorOperation::values operation);
707 
708 
714  bool
715  operator==(const Vector &v) const;
716 
722  bool
723  operator!=(const Vector &v) const;
724 
728  size_type
729  size() const;
730 
742  size_type
743  local_size() const;
744 
766  std::pair<size_type, size_type>
767  local_range() const;
768 
776  bool
777  in_local_range(const size_type index) const;
778 
792  IndexSet
793  locally_owned_elements() const;
794 
802  bool
803  has_ghost_elements() const;
804 
811  void
812  update_ghost_values() const;
813 
818  TrilinosScalar operator*(const Vector &vec) const;
819 
823  real_type
824  norm_sqr() const;
825 
830  mean_value() const;
831 
836  min() const;
837 
842  max() const;
843 
847  real_type
848  l1_norm() const;
849 
854  real_type
855  l2_norm() const;
856 
861  real_type
862  lp_norm(const TrilinosScalar p) const;
863 
867  real_type
868  linfty_norm() const;
869 
890  add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
891 
897  bool
898  all_zero() const;
899 
905  bool
906  is_non_negative() const;
908 
909 
914 
922  reference
923  operator()(const size_type index);
924 
933  operator()(const size_type index) const;
934 
940  reference operator[](const size_type index);
941 
947  TrilinosScalar operator[](const size_type index) const;
948 
964  void
965  extract_subvector_to(const std::vector<size_type> &indices,
966  std::vector<TrilinosScalar> & values) const;
967 
995  template <typename ForwardIterator, typename OutputIterator>
996  void
997  extract_subvector_to(ForwardIterator indices_begin,
998  const ForwardIterator indices_end,
999  OutputIterator values_begin) const;
1000 
1011  iterator
1012  begin();
1013 
1019  begin() const;
1020 
1025  iterator
1026  end();
1027 
1033  end() const;
1034 
1036 
1037 
1042 
1049  void
1050  set(const std::vector<size_type> & indices,
1051  const std::vector<TrilinosScalar> &values);
1052 
1057  void
1058  set(const std::vector<size_type> & indices,
1059  const ::Vector<TrilinosScalar> &values);
1060 
1066  void
1067  set(const size_type n_elements,
1068  const size_type * indices,
1069  const TrilinosScalar *values);
1070 
1075  void
1076  add(const std::vector<size_type> & indices,
1077  const std::vector<TrilinosScalar> &values);
1078 
1083  void
1084  add(const std::vector<size_type> & indices,
1085  const ::Vector<TrilinosScalar> &values);
1086 
1092  void
1093  add(const size_type n_elements,
1094  const size_type * indices,
1095  const TrilinosScalar *values);
1096 
1100  Vector &
1101  operator*=(const TrilinosScalar factor);
1102 
1106  Vector &
1107  operator/=(const TrilinosScalar factor);
1108 
1112  Vector &
1113  operator+=(const Vector &V);
1114 
1118  Vector &
1119  operator-=(const Vector &V);
1120 
1125  void
1126  add(const TrilinosScalar s);
1127 
1140  void
1141  add(const Vector &V, const bool allow_different_maps = false);
1142 
1146  void
1147  add(const TrilinosScalar a, const Vector &V);
1148 
1152  void
1153  add(const TrilinosScalar a,
1154  const Vector & V,
1155  const TrilinosScalar b,
1156  const Vector & W);
1157 
1162  void
1163  sadd(const TrilinosScalar s, const Vector &V);
1164 
1168  void
1169  sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1170 
1176  void
1177  scale(const Vector &scaling_factors);
1178 
1182  void
1183  equ(const TrilinosScalar a, const Vector &V);
1185 
1190 
1195  const Epetra_MultiVector &
1196  trilinos_vector() const;
1197 
1202  Epetra_FEVector &
1203  trilinos_vector();
1204 
1209  const Epetra_BlockMap &
1210  trilinos_partitioner() const;
1211 
1219  void
1220  print(std::ostream & out,
1221  const unsigned int precision = 3,
1222  const bool scientific = true,
1223  const bool across = true) const;
1224 
1238  void
1239  swap(Vector &v);
1240 
1244  std::size_t
1245  memory_consumption() const;
1246 
1251  const MPI_Comm &
1252  get_mpi_communicator() const;
1254 
1258  DeclException0(ExcDifferentParallelPartitioning);
1259 
1264  int,
1265  << "An error with error number " << arg1
1266  << " occurred while calling a Trilinos function");
1267 
1272  ExcAccessToNonLocalElement,
1273  size_type,
1274  size_type,
1275  size_type,
1276  size_type,
1277  << "You tried to access element " << arg1
1278  << " of a distributed vector, but this element is not stored "
1279  << "on the current processor. Note: There are " << arg2
1280  << " elements stored "
1281  << "on the current processor from within the range " << arg3
1282  << " through " << arg4
1283  << " but Trilinos vectors need not store contiguous "
1284  << "ranges on each processor, and not every element in "
1285  << "this range may in fact be stored locally.");
1286 
1287  private:
1299  Epetra_CombineMode last_action;
1300 
1306 
1312 
1318  std::unique_ptr<Epetra_FEVector> vector;
1319 
1325  std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1326 
1331 
1332  // Make the reference class a friend.
1333  friend class internal::VectorReference;
1334  };
1335 
1336 
1337 
1338  // ------------------- inline and template functions --------------
1339 
1340 
1348  inline void
1350  {
1351  u.swap(v);
1352  }
1353  } // namespace MPI
1354 
1355 # ifndef DOXYGEN
1356 
1357  namespace internal
1358  {
1359  inline VectorReference::VectorReference(MPI::Vector & vector,
1360  const size_type index)
1361  : vector(vector)
1362  , index(index)
1363  {}
1364 
1365 
1366  inline const VectorReference &
1367  VectorReference::operator=(const VectorReference &r) const
1368  {
1369  // as explained in the class
1370  // documentation, this is not the copy
1371  // operator. so simply pass on to the
1372  // "correct" assignment operator
1373  *this = static_cast<TrilinosScalar>(r);
1374 
1375  return *this;
1376  }
1377 
1378 
1379 
1380  inline VectorReference &
1381  VectorReference::operator=(const VectorReference &r)
1382  {
1383  // as above
1384  *this = static_cast<TrilinosScalar>(r);
1385 
1386  return *this;
1387  }
1388 
1389 
1390  inline const VectorReference &
1391  VectorReference::operator=(const TrilinosScalar &value) const
1392  {
1393  vector.set(1, &index, &value);
1394  return *this;
1395  }
1396 
1397 
1398 
1399  inline const VectorReference &
1400  VectorReference::operator+=(const TrilinosScalar &value) const
1401  {
1402  vector.add(1, &index, &value);
1403  return *this;
1404  }
1405 
1406 
1407 
1408  inline const VectorReference &
1409  VectorReference::operator-=(const TrilinosScalar &value) const
1410  {
1411  TrilinosScalar new_value = -value;
1412  vector.add(1, &index, &new_value);
1413  return *this;
1414  }
1415 
1416 
1417 
1418  inline const VectorReference &
1419  VectorReference::operator*=(const TrilinosScalar &value) const
1420  {
1421  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1422  vector.set(1, &index, &new_value);
1423  return *this;
1424  }
1425 
1426 
1427 
1428  inline const VectorReference &
1429  VectorReference::operator/=(const TrilinosScalar &value) const
1430  {
1431  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1432  vector.set(1, &index, &new_value);
1433  return *this;
1434  }
1435  } // namespace internal
1436 
1437  namespace MPI
1438  {
1439  inline bool
1440  Vector::in_local_range(const size_type index) const
1441  {
1442  std::pair<size_type, size_type> range = local_range();
1443 
1444  return ((index >= range.first) && (index < range.second));
1445  }
1446 
1447 
1448 
1449  inline IndexSet
1451  {
1452  Assert(owned_elements.size() == size(),
1453  ExcMessage(
1454  "The locally owned elements have not been properly initialized!"
1455  " This happens for example if this object has been initialized"
1456  " with exactly one overlapping IndexSet."));
1457  return owned_elements;
1458  }
1459 
1460 
1461 
1462  inline bool
1464  {
1465  return has_ghosts;
1466  }
1467 
1468 
1469 
1470  inline void
1472  {}
1473 
1474 
1475 
1476  inline internal::VectorReference
1477  Vector::operator()(const size_type index)
1478  {
1479  return internal::VectorReference(*this, index);
1480  }
1481 
1482 
1483 
1484  inline internal::VectorReference Vector::operator[](const size_type index)
1485  {
1486  return operator()(index);
1487  }
1488 
1489 
1490 
1491  inline TrilinosScalar Vector::operator[](const size_type index) const
1492  {
1493  return operator()(index);
1494  }
1495 
1496 
1497 
1498  inline void
1499  Vector::extract_subvector_to(const std::vector<size_type> &indices,
1500  std::vector<TrilinosScalar> & values) const
1501  {
1502  for (size_type i = 0; i < indices.size(); ++i)
1503  values[i] = operator()(indices[i]);
1504  }
1505 
1506 
1507 
1508  template <typename ForwardIterator, typename OutputIterator>
1509  inline void
1510  Vector::extract_subvector_to(ForwardIterator indices_begin,
1511  const ForwardIterator indices_end,
1512  OutputIterator values_begin) const
1513  {
1514  while (indices_begin != indices_end)
1515  {
1516  *values_begin = operator()(*indices_begin);
1517  indices_begin++;
1518  values_begin++;
1519  }
1520  }
1521 
1522 
1523 
1524  inline Vector::iterator
1525  Vector::begin()
1526  {
1527  return (*vector)[0];
1528  }
1529 
1530 
1531 
1532  inline Vector::iterator
1533  Vector::end()
1534  {
1535  return (*vector)[0] + local_size();
1536  }
1537 
1538 
1539 
1540  inline Vector::const_iterator
1541  Vector::begin() const
1542  {
1543  return (*vector)[0];
1544  }
1545 
1546 
1547 
1548  inline Vector::const_iterator
1549  Vector::end() const
1550  {
1551  return (*vector)[0] + local_size();
1552  }
1553 
1554 
1555 
1556  inline void
1557  Vector::set(const std::vector<size_type> & indices,
1558  const std::vector<TrilinosScalar> &values)
1559  {
1560  // if we have ghost values, do not allow
1561  // writing to this vector at all.
1562  Assert(!has_ghost_elements(), ExcGhostsPresent());
1563 
1564  Assert(indices.size() == values.size(),
1565  ExcDimensionMismatch(indices.size(), values.size()));
1566 
1567  set(indices.size(), indices.data(), values.data());
1568  }
1569 
1570 
1571 
1572  inline void
1573  Vector::set(const std::vector<size_type> & indices,
1574  const ::Vector<TrilinosScalar> &values)
1575  {
1576  // if we have ghost values, do not allow
1577  // writing to this vector at all.
1578  Assert(!has_ghost_elements(), ExcGhostsPresent());
1579 
1580  Assert(indices.size() == values.size(),
1581  ExcDimensionMismatch(indices.size(), values.size()));
1582 
1583  set(indices.size(), indices.data(), values.begin());
1584  }
1585 
1586 
1587 
1588  inline void
1589  Vector::set(const size_type n_elements,
1590  const size_type * indices,
1591  const TrilinosScalar *values)
1592  {
1593  // if we have ghost values, do not allow
1594  // writing to this vector at all.
1595  Assert(!has_ghost_elements(), ExcGhostsPresent());
1596 
1597  if (last_action == Add)
1598  {
1599  const int ierr = vector->GlobalAssemble(Add);
1600  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1601  }
1602 
1603  if (last_action != Insert)
1604  last_action = Insert;
1605 
1606  for (size_type i = 0; i < n_elements; ++i)
1607  {
1608  const TrilinosWrappers::types::int_type row = indices[i];
1609  const TrilinosWrappers::types::int_type local_row =
1610  vector->Map().LID(row);
1611  if (local_row != -1)
1612  (*vector)[0][local_row] = values[i];
1613  else
1614  {
1615  const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1616  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1617  compressed = false;
1618  }
1619  // in set operation, do not use the pre-allocated vector for nonlocal
1620  // entries even if it exists. This is to ensure that we really only
1621  // set the elements touched by the set() method and not all contained
1622  // in the nonlocal entries vector (there is no way to distinguish them
1623  // on the receiving processor)
1624  }
1625  }
1626 
1627 
1628 
1629  inline void
1630  Vector::add(const std::vector<size_type> & indices,
1631  const std::vector<TrilinosScalar> &values)
1632  {
1633  // if we have ghost values, do not allow
1634  // writing to this vector at all.
1635  Assert(!has_ghost_elements(), ExcGhostsPresent());
1636  Assert(indices.size() == values.size(),
1637  ExcDimensionMismatch(indices.size(), values.size()));
1638 
1639  add(indices.size(), indices.data(), values.data());
1640  }
1641 
1642 
1643 
1644  inline void
1645  Vector::add(const std::vector<size_type> & indices,
1646  const ::Vector<TrilinosScalar> &values)
1647  {
1648  // if we have ghost values, do not allow
1649  // writing to this vector at all.
1650  Assert(!has_ghost_elements(), ExcGhostsPresent());
1651  Assert(indices.size() == values.size(),
1652  ExcDimensionMismatch(indices.size(), values.size()));
1653 
1654  add(indices.size(), indices.data(), values.begin());
1655  }
1656 
1657 
1658 
1659  inline void
1660  Vector::add(const size_type n_elements,
1661  const size_type * indices,
1662  const TrilinosScalar *values)
1663  {
1664  // if we have ghost values, do not allow
1665  // writing to this vector at all.
1666  Assert(!has_ghost_elements(), ExcGhostsPresent());
1667 
1668  if (last_action != Add)
1669  {
1670  if (last_action == Insert)
1671  {
1672  const int ierr = vector->GlobalAssemble(Insert);
1673  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1674  }
1675  last_action = Add;
1676  }
1677 
1678  for (size_type i = 0; i < n_elements; ++i)
1679  {
1680  const size_type row = indices[i];
1681  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1682  static_cast<TrilinosWrappers::types::int_type>(row));
1683  if (local_row != -1)
1684  (*vector)[0][local_row] += values[i];
1685  else if (nonlocal_vector.get() == nullptr)
1686  {
1687  const int ierr = vector->SumIntoGlobalValues(
1688  1,
1689  reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1690  &row),
1691  &values[i]);
1692  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1693  compressed = false;
1694  }
1695  else
1696  {
1697  // use pre-allocated vector for non-local entries if it exists for
1698  // addition operation
1699  const TrilinosWrappers::types::int_type my_row =
1700  nonlocal_vector->Map().LID(
1701  static_cast<TrilinosWrappers::types::int_type>(row));
1702  Assert(my_row != -1,
1703  ExcMessage(
1704  "Attempted to write into off-processor vector entry "
1705  "that has not be specified as being writable upon "
1706  "initialization"));
1707  (*nonlocal_vector)[0][my_row] += values[i];
1708  compressed = false;
1709  }
1710  }
1711  }
1712 
1713 
1714 
1715  inline Vector::size_type
1716  Vector::size() const
1717  {
1718 # ifndef DEAL_II_WITH_64BIT_INDICES
1719  return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1720 # else
1721  return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1722 # endif
1723  }
1724 
1725 
1726 
1727  inline Vector::size_type
1728  Vector::local_size() const
1729  {
1730  return vector->Map().NumMyElements();
1731  }
1732 
1733 
1734 
1735  inline std::pair<Vector::size_type, Vector::size_type>
1736  Vector::local_range() const
1737  {
1738 # ifndef DEAL_II_WITH_64BIT_INDICES
1739  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1741  vector->Map().MaxMyGID() + 1;
1742 # else
1743  const TrilinosWrappers::types::int_type begin =
1744  vector->Map().MinMyGID64();
1746  vector->Map().MaxMyGID64() + 1;
1747 # endif
1748 
1749  Assert(
1750  end - begin == vector->Map().NumMyElements(),
1751  ExcMessage(
1752  "This function only makes sense if the elements that this "
1753  "vector stores on the current processor form a contiguous range. "
1754  "This does not appear to be the case for the current vector."));
1755 
1756  return std::make_pair(begin, end);
1757  }
1758 
1759 
1760 
1761  inline TrilinosScalar Vector::operator*(const Vector &vec) const
1762  {
1763  Assert(vector->Map().SameAs(vec.vector->Map()),
1764  ExcDifferentParallelPartitioning());
1765  Assert(!has_ghost_elements(), ExcGhostsPresent());
1766 
1767  TrilinosScalar result;
1768 
1769  const int ierr = vector->Dot(*(vec.vector), &result);
1770  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1771 
1772  return result;
1773  }
1774 
1775 
1776 
1777  inline Vector::real_type
1778  Vector::norm_sqr() const
1779  {
1780  const TrilinosScalar d = l2_norm();
1781  return d * d;
1782  }
1783 
1784 
1785 
1786  inline TrilinosScalar
1787  Vector::mean_value() const
1788  {
1789  Assert(!has_ghost_elements(), ExcGhostsPresent());
1790 
1792  const int ierr = vector->MeanValue(&mean);
1793  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1794 
1795  return mean;
1796  }
1797 
1798 
1799 
1800  inline TrilinosScalar
1801  Vector::min() const
1802  {
1803  TrilinosScalar min_value;
1804  const int ierr = vector->MinValue(&min_value);
1805  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1806 
1807  return min_value;
1808  }
1809 
1810 
1811 
1812  inline TrilinosScalar
1813  Vector::max() const
1814  {
1815  TrilinosScalar max_value;
1816  const int ierr = vector->MaxValue(&max_value);
1817  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1818 
1819  return max_value;
1820  }
1821 
1822 
1823 
1824  inline Vector::real_type
1825  Vector::l1_norm() const
1826  {
1827  Assert(!has_ghost_elements(), ExcGhostsPresent());
1828 
1829  TrilinosScalar d;
1830  const int ierr = vector->Norm1(&d);
1831  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1832 
1833  return d;
1834  }
1835 
1836 
1837 
1838  inline Vector::real_type
1839  Vector::l2_norm() const
1840  {
1841  Assert(!has_ghost_elements(), ExcGhostsPresent());
1842 
1843  TrilinosScalar d;
1844  const int ierr = vector->Norm2(&d);
1845  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1846 
1847  return d;
1848  }
1849 
1850 
1851 
1852  inline Vector::real_type
1853  Vector::lp_norm(const TrilinosScalar p) const
1854  {
1855  Assert(!has_ghost_elements(), ExcGhostsPresent());
1856 
1857  TrilinosScalar norm = 0;
1858  TrilinosScalar sum = 0;
1859  const size_type n_local = local_size();
1860 
1861  // loop over all the elements because
1862  // Trilinos does not support lp norms
1863  for (size_type i = 0; i < n_local; ++i)
1864  sum += std::pow(std::fabs((*vector)[0][i]), p);
1865 
1866  norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1867 
1868  return norm;
1869  }
1870 
1871 
1872 
1873  inline Vector::real_type
1874  Vector::linfty_norm() const
1875  {
1876  // while we disallow the other
1877  // norm operations on ghosted
1878  // vectors, this particular norm
1879  // is safe to run even in the
1880  // presence of ghost elements
1881  TrilinosScalar d;
1882  const int ierr = vector->NormInf(&d);
1883  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1884 
1885  return d;
1886  }
1887 
1888 
1889 
1890  inline TrilinosScalar
1892  const Vector & V,
1893  const Vector & W)
1894  {
1895  this->add(a, V);
1896  return *this * W;
1897  }
1898 
1899 
1900 
1901  // inline also scalar products, vector
1902  // additions etc. since they are all
1903  // representable by a single Trilinos
1904  // call. This reduces the overhead of the
1905  // wrapper class.
1906  inline Vector &
1908  {
1909  AssertIsFinite(a);
1910 
1911  const int ierr = vector->Scale(a);
1912  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1913 
1914  return *this;
1915  }
1916 
1917 
1918 
1919  inline Vector &
1921  {
1922  AssertIsFinite(a);
1923 
1924  const TrilinosScalar factor = 1. / a;
1925 
1926  AssertIsFinite(factor);
1927 
1928  const int ierr = vector->Scale(factor);
1929  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1930 
1931  return *this;
1932  }
1933 
1934 
1935 
1936  inline Vector &
1937  Vector::operator+=(const Vector &v)
1938  {
1939  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
1940  Assert(vector->Map().SameAs(v.vector->Map()),
1941  ExcDifferentParallelPartitioning());
1942 
1943  const int ierr = vector->Update(1.0, *(v.vector), 1.0);
1944  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1945 
1946  return *this;
1947  }
1948 
1949 
1950 
1951  inline Vector &
1952  Vector::operator-=(const Vector &v)
1953  {
1954  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
1955  Assert(vector->Map().SameAs(v.vector->Map()),
1956  ExcDifferentParallelPartitioning());
1957 
1958  const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
1959  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1960 
1961  return *this;
1962  }
1963 
1964 
1965 
1966  inline void
1967  Vector::add(const TrilinosScalar s)
1968  {
1969  // if we have ghost values, do not allow
1970  // writing to this vector at all.
1971  Assert(!has_ghost_elements(), ExcGhostsPresent());
1972  AssertIsFinite(s);
1973 
1974  size_type n_local = local_size();
1975  for (size_type i = 0; i < n_local; ++i)
1976  (*vector)[0][i] += s;
1977  }
1978 
1979 
1980 
1981  inline void
1982  Vector::add(const TrilinosScalar a, const Vector &v)
1983  {
1984  // if we have ghost values, do not allow
1985  // writing to this vector at all.
1986  Assert(!has_ghost_elements(), ExcGhostsPresent());
1987  Assert(local_size() == v.local_size(),
1988  ExcDimensionMismatch(local_size(), v.local_size()));
1989 
1990  AssertIsFinite(a);
1991 
1992  const int ierr = vector->Update(a, *(v.vector), 1.);
1993  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1994  }
1995 
1996 
1997 
1998  inline void
1999  Vector::add(const TrilinosScalar a,
2000  const Vector & v,
2001  const TrilinosScalar b,
2002  const Vector & w)
2003  {
2004  // if we have ghost values, do not allow
2005  // writing to this vector at all.
2006  Assert(!has_ghost_elements(), ExcGhostsPresent());
2007  Assert(local_size() == v.local_size(),
2008  ExcDimensionMismatch(local_size(), v.local_size()));
2009  Assert(local_size() == w.local_size(),
2010  ExcDimensionMismatch(local_size(), w.local_size()));
2011 
2012  AssertIsFinite(a);
2013  AssertIsFinite(b);
2014 
2015  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2016 
2017  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2018  }
2019 
2020 
2021 
2022  inline void
2023  Vector::sadd(const TrilinosScalar s, const Vector &v)
2024  {
2025  // if we have ghost values, do not allow
2026  // writing to this vector at all.
2027  Assert(!has_ghost_elements(), ExcGhostsPresent());
2028  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2029 
2030  AssertIsFinite(s);
2031 
2032  // We assume that the vectors have the same Map
2033  // if the local size is the same and if the vectors are not ghosted
2034  if (local_size() == v.local_size() && !v.has_ghost_elements())
2035  {
2036  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2037  ExcDifferentParallelPartitioning());
2038  const int ierr = vector->Update(1., *(v.vector), s);
2039  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2040  }
2041  else
2042  {
2043  (*this) *= s;
2044  this->add(v, true);
2045  }
2046  }
2047 
2048 
2049 
2050  inline void
2051  Vector::sadd(const TrilinosScalar s,
2052  const TrilinosScalar a,
2053  const Vector & v)
2054  {
2055  // if we have ghost values, do not allow
2056  // writing to this vector at all.
2057  Assert(!has_ghost_elements(), ExcGhostsPresent());
2058  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2059  AssertIsFinite(s);
2060  AssertIsFinite(a);
2061 
2062  // We assume that the vectors have the same Map
2063  // if the local size is the same and if the vectors are not ghosted
2064  if (local_size() == v.local_size() && !v.has_ghost_elements())
2065  {
2066  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2067  ExcDifferentParallelPartitioning());
2068  const int ierr = vector->Update(a, *(v.vector), s);
2069  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2070  }
2071  else
2072  {
2073  (*this) *= s;
2074  Vector tmp = v;
2075  tmp *= a;
2076  this->add(tmp, true);
2077  }
2078  }
2079 
2080 
2081 
2082  inline void
2083  Vector::scale(const Vector &factors)
2084  {
2085  // if we have ghost values, do not allow
2086  // writing to this vector at all.
2087  Assert(!has_ghost_elements(), ExcGhostsPresent());
2088  Assert(local_size() == factors.local_size(),
2089  ExcDimensionMismatch(local_size(), factors.local_size()));
2090 
2091  const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2092  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2093  }
2094 
2095 
2096 
2097  inline void
2098  Vector::equ(const TrilinosScalar a, const Vector &v)
2099  {
2100  // if we have ghost values, do not allow
2101  // writing to this vector at all.
2102  Assert(!has_ghost_elements(), ExcGhostsPresent());
2103  AssertIsFinite(a);
2104 
2105  // If we don't have the same map, copy.
2106  if (vector->Map().SameAs(v.vector->Map()) == false)
2107  {
2108  this->sadd(0., a, v);
2109  }
2110  else
2111  {
2112  // Otherwise, just update
2113  int ierr = vector->Update(a, *v.vector, 0.0);
2114  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2115 
2116  last_action = Zero;
2117  }
2118  }
2119 
2120 
2121 
2122  inline const Epetra_MultiVector &
2123  Vector::trilinos_vector() const
2124  {
2125  return static_cast<const Epetra_MultiVector &>(*vector);
2126  }
2127 
2128 
2129 
2130  inline Epetra_FEVector &
2131  Vector::trilinos_vector()
2132  {
2133  return *vector;
2134  }
2135 
2136 
2137 
2138  inline const Epetra_BlockMap &
2139  Vector::trilinos_partitioner() const
2140  {
2141  return vector->Map();
2142  }
2143 
2144 
2145 
2146  inline const MPI_Comm &
2147  Vector::get_mpi_communicator() const
2148  {
2149  static MPI_Comm comm;
2150 
2151 # ifdef DEAL_II_WITH_MPI
2152 
2153  const Epetra_MpiComm *mpi_comm =
2154  dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2155  comm = mpi_comm->Comm();
2156 
2157 # else
2158 
2159  comm = MPI_COMM_SELF;
2160 
2161 # endif
2162 
2163  return comm;
2164  }
2165 
2166  template <typename number>
2167  Vector::Vector(const IndexSet & parallel_partitioner,
2168  const ::Vector<number> &v,
2169  const MPI_Comm & communicator)
2170  {
2171  *this =
2172  Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2173  owned_elements = parallel_partitioner;
2174  }
2175 
2176 
2177 
2178  inline Vector &
2180  {
2181  AssertIsFinite(s);
2182 
2183  int ierr = vector->PutScalar(s);
2184  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2185 
2186  if (nonlocal_vector.get() != nullptr)
2187  {
2188  ierr = nonlocal_vector->PutScalar(0.);
2189  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2190  }
2191 
2192  return *this;
2193  }
2194  } /* end of namespace MPI */
2195 
2196 # endif /* DOXYGEN */
2197 
2198 } /* end of namespace TrilinosWrappers */
2199 
2203 namespace internal
2204 {
2205  namespace LinearOperatorImplementation
2206  {
2207  template <typename>
2208  class ReinitHelper;
2209 
2214  template <>
2216  {
2217  public:
2218  template <typename Matrix>
2219  static void
2220  reinit_range_vector(const Matrix & matrix,
2222  bool omit_zeroing_entries)
2223  {
2224  v.reinit(matrix.locally_owned_range_indices(),
2225  matrix.get_mpi_communicator(),
2226  omit_zeroing_entries);
2227  }
2228 
2229  template <typename Matrix>
2230  static void
2233  bool omit_zeroing_entries)
2234  {
2235  v.reinit(matrix.locally_owned_domain_indices(),
2236  matrix.get_mpi_communicator(),
2237  omit_zeroing_entries);
2238  }
2239  };
2240 
2241  } // namespace LinearOperatorImplementation
2242 } /* namespace internal */
2243 
2244 
2245 
2249 template <>
2251 {};
2252 
2253 
2255 
2256 #endif // DEAL_II_WITH_TRILINOS
2257 
2258 /*---------------------------- trilinos_vector.h ---------------------------*/
2259 
2260 #endif // dealii_trilinos_vector_h
std::unique_ptr< Epetra_MultiVector > nonlocal_vector
static ::ExceptionBase & ExcTrilinosError(int arg1)
void update_ghost_values() const
int gid(const Epetra_BlockMap &map, int i)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
Number operator[](const size_type global_index) const
virtual VectorSpaceVector< Number >::real_type linfty_norm() const override
::types::global_dof_index size_type
Contents is actually a matrix.
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:839
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:1513
virtual Number operator*(const VectorSpaceVector< Number > &V) const override
virtual Vector< Number > & operator+=(const VectorSpaceVector< Number > &V) override
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
virtual void equ(const Number a, const VectorSpaceVector< Number > &V) override
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:596
virtual ::IndexSet locally_owned_elements() const override
Definition: la_vector.h:465
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:1403
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::unique_ptr< Epetra_FEVector > vector
real_type lp_norm(const real_type p) const
Number operator()(const size_type global_index) const
void swap(Vector &u, Vector &v)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2779
#define DeclException0(Exception0)
Definition: exceptions.h:470
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
virtual Vector< Number > & operator/=(const Number factor) override
virtual void scale(const VectorSpaceVector< Number > &scaling_factors) override
VectorType::value_type * end(VectorType &V)
double TrilinosScalar
Definition: types.h:158
virtual size_type size() const override
Definition: la_vector.h:456
Expression fabs(const Expression &x)
__global__ void equ(Number *val, const Number a, const Number *V_val, const size_type N)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SparseMatrix< double > SparseMatrix
__global__ void sadd(const Number s, Number *val, const Number a, const Number *V_val, const size_type N)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
static ::ExceptionBase & ExcGhostsPresent()
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int global_dof_index
Definition: types.h:76
__global__ void set(Number *val, const Number s, const size_type N)
virtual VectorSpaceVector< Number >::real_type l1_norm() const override
size_type local_size() const
virtual VectorSpaceVector< Number >::real_type l2_norm() const override
*braid_SplitCommworld & comm
internal::VectorReference reference
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
VectorType::value_type * begin(VectorType &V)
T min(const T &t, const MPI_Comm &mpi_communicator)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:584
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2753
virtual value_type mean_value() const override
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< Number2 > &values) const
real_type norm_sqr() const
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
virtual Number add_and_dot(const Number a, const VectorSpaceVector< Number > &V, const VectorSpaceVector< Number > &W) override
bool has_ghost_elements() const
#define AssertIsFinite(number)
Definition: exceptions.h:1659
T max(const T &t, const MPI_Comm &mpi_communicator)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)