Reference documentation for deal.II version Git 95a79d5 2018-06-20 15:34:22 +0200
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 
61 namespace TrilinosWrappers
62 {
63  class SparseMatrix;
64 
75  namespace internal
76  {
80  typedef ::types::global_dof_index size_type;
81 
91  class VectorReference
92  {
93  private:
98  VectorReference(MPI::Vector &vector, const size_type index);
99 
100  public:
112  const VectorReference &
113  operator=(const VectorReference &r) const;
114 
118  VectorReference &
119  operator=(const VectorReference &r);
120 
124  const VectorReference &
125  operator=(const TrilinosScalar &s) const;
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 
155  operator TrilinosScalar() const;
156 
161  int,
162  << "An error with error number " << arg1
163  << " occurred while calling a Trilinos function");
164 
165  private:
169  MPI::Vector &vector;
170 
174  const size_type index;
175 
181  };
182  } // namespace internal
187 # ifndef DEAL_II_WITH_64BIT_INDICES
188  // define a helper function that queries the global ID of local ID of
189  // an Epetra_BlockMap object by calling either the 32- or 64-bit
190  // function necessary.
191  inline int
192  gid(const Epetra_BlockMap &map, int i)
193  {
194  return map.GID(i);
195  }
196 # else
197  // define a helper function that queries the global ID of local ID of
198  // an Epetra_BlockMap object by calling either the 32- or 64-bit
199  // function necessary.
200  inline long long int
201  gid(const Epetra_BlockMap &map, int i)
202  {
203  return map.GID64(i);
204  }
205 # endif
206 
213  namespace MPI
214  {
215  class BlockVector;
216 
393  class Vector : public Subscriptor
394  {
395  public:
401  typedef TrilinosScalar value_type;
402  typedef TrilinosScalar real_type;
403  typedef ::types::global_dof_index size_type;
404  typedef value_type * iterator;
405  typedef const value_type * const_iterator;
406  typedef internal::VectorReference reference;
407  typedef const internal::VectorReference const_reference;
408 
418  Vector();
419 
423  Vector(const Vector &v);
424 
443  explicit Vector(const IndexSet &parallel_partitioning,
444  const MPI_Comm &communicator = MPI_COMM_WORLD);
445 
457  Vector(const IndexSet &local,
458  const IndexSet &ghost,
459  const MPI_Comm &communicator = MPI_COMM_WORLD);
460 
475  Vector(const IndexSet &parallel_partitioning,
476  const Vector & v,
477  const MPI_Comm &communicator = MPI_COMM_WORLD);
478 
491  template <typename Number>
492  Vector(const IndexSet & parallel_partitioning,
493  const ::Vector<Number> &v,
494  const MPI_Comm & communicator = MPI_COMM_WORLD);
495 
500  Vector(Vector &&v) noexcept;
501 
505  ~Vector() override = default;
506 
511  void
512  clear();
513 
537  void
538  reinit(const Vector &v,
539  const bool omit_zeroing_entries = false,
540  const bool allow_different_maps = false);
541 
564  void
565  reinit(const IndexSet &parallel_partitioning,
566  const MPI_Comm &communicator = MPI_COMM_WORLD,
567  const bool omit_zeroing_entries = false);
568 
594  void
595  reinit(const IndexSet &locally_owned_entries,
596  const IndexSet &ghost_entries,
597  const MPI_Comm &communicator = MPI_COMM_WORLD,
598  const bool vector_writable = false);
599 
603  void
604  reinit(const BlockVector &v, const bool import_data = false);
605 
622  void
623  compress(::VectorOperation::values operation);
624 
637  Vector &
638  operator=(const TrilinosScalar s);
639 
645  Vector &
646  operator=(const Vector &v);
647 
652  Vector &
653  operator=(Vector &&v) noexcept;
654 
662  template <typename Number>
663  Vector &
664  operator=(const ::Vector<Number> &v);
665 
683  void
685  const ::TrilinosWrappers::SparseMatrix &matrix,
686  const Vector & vector);
687 
693  bool
694  operator==(const Vector &v) const;
695 
701  bool
702  operator!=(const Vector &v) const;
703 
707  size_type
708  size() const;
709 
721  size_type
722  local_size() const;
723 
745  std::pair<size_type, size_type>
746  local_range() const;
747 
755  bool
756  in_local_range(const size_type index) const;
757 
771  IndexSet
772  locally_owned_elements() const;
773 
781  bool
782  has_ghost_elements() const;
783 
790  void
791  update_ghost_values() const;
792 
797  TrilinosScalar operator*(const Vector &vec) const;
798 
802  real_type
803  norm_sqr() const;
804 
808  TrilinosScalar
809  mean_value() const;
810 
814  TrilinosScalar
815  min() const;
816 
820  TrilinosScalar
821  max() const;
822 
826  real_type
827  l1_norm() const;
828 
833  real_type
834  l2_norm() const;
835 
840  real_type
841  lp_norm(const TrilinosScalar p) const;
842 
846  real_type
847  linfty_norm() const;
848 
868  TrilinosScalar
869  add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
870 
876  bool
877  all_zero() const;
878 
884  bool
885  is_non_negative() const;
887 
888 
893 
901  reference
902  operator()(const size_type index);
903 
911  TrilinosScalar
912  operator()(const size_type index) const;
913 
919  reference operator[](const size_type index);
920 
926  TrilinosScalar operator[](const size_type index) const;
927 
943  void
944  extract_subvector_to(const std::vector<size_type> &indices,
945  std::vector<TrilinosScalar> & values) const;
946 
974  template <typename ForwardIterator, typename OutputIterator>
975  void
976  extract_subvector_to(ForwardIterator indices_begin,
977  const ForwardIterator indices_end,
978  OutputIterator values_begin) const;
979 
990  iterator
991  begin();
992 
997  const_iterator
998  begin() const;
999 
1004  iterator
1005  end();
1006 
1011  const_iterator
1012  end() const;
1013 
1015 
1016 
1021 
1028  void
1029  set(const std::vector<size_type> & indices,
1030  const std::vector<TrilinosScalar> &values);
1031 
1036  void
1037  set(const std::vector<size_type> & indices,
1038  const ::Vector<TrilinosScalar> &values);
1039 
1045  void
1046  set(const size_type n_elements,
1047  const size_type * indices,
1048  const TrilinosScalar *values);
1049 
1054  void
1055  add(const std::vector<size_type> & indices,
1056  const std::vector<TrilinosScalar> &values);
1057 
1062  void
1063  add(const std::vector<size_type> & indices,
1064  const ::Vector<TrilinosScalar> &values);
1065 
1071  void
1072  add(const size_type n_elements,
1073  const size_type * indices,
1074  const TrilinosScalar *values);
1075 
1079  Vector &
1080  operator*=(const TrilinosScalar factor);
1081 
1085  Vector &
1086  operator/=(const TrilinosScalar factor);
1087 
1091  Vector &
1092  operator+=(const Vector &V);
1093 
1097  Vector &
1098  operator-=(const Vector &V);
1099 
1104  void
1105  add(const TrilinosScalar s);
1106 
1119  void
1120  add(const Vector &V, const bool allow_different_maps = false);
1121 
1125  void
1126  add(const TrilinosScalar a, const Vector &V);
1127 
1131  void
1132  add(const TrilinosScalar a,
1133  const Vector & V,
1134  const TrilinosScalar b,
1135  const Vector & W);
1136 
1141  void
1142  sadd(const TrilinosScalar s, const Vector &V);
1143 
1147  void
1148  sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1149 
1155  void
1156  scale(const Vector &scaling_factors);
1157 
1161  void
1162  equ(const TrilinosScalar a, const Vector &V);
1164 
1169 
1174  const Epetra_MultiVector &
1175  trilinos_vector() const;
1176 
1181  Epetra_FEVector &
1182  trilinos_vector();
1183 
1188  const Epetra_Map &
1189  vector_partitioner() const;
1190 
1198  void
1199  print(std::ostream & out,
1200  const unsigned int precision = 3,
1201  const bool scientific = true,
1202  const bool across = true) const;
1203 
1217  void
1218  swap(Vector &v);
1219 
1223  std::size_t
1224  memory_consumption() const;
1225 
1230  const MPI_Comm &
1231  get_mpi_communicator() const;
1233 
1238 
1243  int,
1244  << "An error with error number " << arg1
1245  << " occurred while calling a Trilinos function");
1246 
1252  size_type,
1253  size_type,
1254  size_type,
1255  size_type,
1256  << "You tried to access element " << arg1
1257  << " of a distributed vector, but this element is not stored "
1258  << "on the current processor. Note: There are " << arg2
1259  << " elements stored "
1260  << "on the current processor from within the range " << arg3
1261  << " through " << arg4
1262  << " but Trilinos vectors need not store contiguous "
1263  << "ranges on each processor, and not every element in "
1264  << "this range may in fact be stored locally.");
1265 
1266  private:
1278  Epetra_CombineMode last_action;
1279 
1285 
1291 
1297  std::unique_ptr<Epetra_FEVector> vector;
1298 
1304  std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1305 
1310 
1314  friend class internal::VectorReference;
1315  };
1316 
1317 
1318 
1319  // ------------------- inline and template functions --------------
1320 
1321 
1330  inline void
1332  {
1333  u.swap(v);
1334  }
1335  } // namespace MPI
1336 
1337 # ifndef DOXYGEN
1338 
1339  namespace internal
1340  {
1341  inline VectorReference::VectorReference(MPI::Vector & vector,
1342  const size_type index)
1343  : vector(vector)
1344  , index(index)
1345  {}
1346 
1347 
1348  inline const VectorReference &
1349  VectorReference::operator=(const VectorReference &r) const
1350  {
1351  // as explained in the class
1352  // documentation, this is not the copy
1353  // operator. so simply pass on to the
1354  // "correct" assignment operator
1355  *this = static_cast<TrilinosScalar>(r);
1356 
1357  return *this;
1358  }
1359 
1360 
1361 
1362  inline VectorReference &
1363  VectorReference::operator=(const VectorReference &r)
1364  {
1365  // as above
1366  *this = static_cast<TrilinosScalar>(r);
1367 
1368  return *this;
1369  }
1370 
1371 
1372  inline const VectorReference &
1373  VectorReference::operator=(const TrilinosScalar &value) const
1374  {
1375  vector.set(1, &index, &value);
1376  return *this;
1377  }
1378 
1379 
1380 
1381  inline const VectorReference &
1382  VectorReference::operator+=(const TrilinosScalar &value) const
1383  {
1384  vector.add(1, &index, &value);
1385  return *this;
1386  }
1387 
1388 
1389 
1390  inline const VectorReference &
1391  VectorReference::operator-=(const TrilinosScalar &value) const
1392  {
1393  TrilinosScalar new_value = -value;
1394  vector.add(1, &index, &new_value);
1395  return *this;
1396  }
1397 
1398 
1399 
1400  inline const VectorReference &
1401  VectorReference::operator*=(const TrilinosScalar &value) const
1402  {
1403  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1404  vector.set(1, &index, &new_value);
1405  return *this;
1406  }
1407 
1408 
1409 
1410  inline const VectorReference &
1411  VectorReference::operator/=(const TrilinosScalar &value) const
1412  {
1413  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1414  vector.set(1, &index, &new_value);
1415  return *this;
1416  }
1417  } // namespace internal
1418 
1419  namespace MPI
1420  {
1421  inline bool
1422  Vector::in_local_range(const size_type index) const
1423  {
1424  std::pair<size_type, size_type> range = local_range();
1425 
1426  return ((index >= range.first) && (index < range.second));
1427  }
1428 
1429 
1430 
1431  inline IndexSet
1433  {
1434  Assert(owned_elements.size() == size(),
1435  ExcMessage(
1436  "The locally owned elements have not been properly initialized!"
1437  " This happens for example if this object has been initialized"
1438  " with exactly one overlapping IndexSet."));
1439  return owned_elements;
1440  }
1441 
1442 
1443 
1444  inline bool
1445  Vector::has_ghost_elements() const
1446  {
1447  return has_ghosts;
1448  }
1449 
1450 
1451 
1452  inline void
1454  {}
1455 
1456 
1457 
1458  inline internal::VectorReference
1459  Vector::operator()(const size_type index)
1460  {
1461  return internal::VectorReference(*this, index);
1462  }
1463 
1464 
1465 
1466  inline internal::VectorReference Vector::operator[](const size_type index)
1467  {
1468  return operator()(index);
1469  }
1470 
1471 
1472 
1473  inline TrilinosScalar Vector::operator[](const size_type index) const
1474  {
1475  return operator()(index);
1476  }
1477 
1478 
1479 
1480  inline void
1481  Vector::extract_subvector_to(const std::vector<size_type> &indices,
1482  std::vector<TrilinosScalar> & values) const
1483  {
1484  for (size_type i = 0; i < indices.size(); ++i)
1485  values[i] = operator()(indices[i]);
1486  }
1487 
1488 
1489 
1490  template <typename ForwardIterator, typename OutputIterator>
1491  inline void
1492  Vector::extract_subvector_to(ForwardIterator indices_begin,
1493  const ForwardIterator indices_end,
1494  OutputIterator values_begin) const
1495  {
1496  while (indices_begin != indices_end)
1497  {
1498  *values_begin = operator()(*indices_begin);
1499  indices_begin++;
1500  values_begin++;
1501  }
1502  }
1503 
1504 
1505 
1506  inline Vector::iterator
1507  Vector::begin()
1508  {
1509  return (*vector)[0];
1510  }
1511 
1512 
1513 
1514  inline Vector::iterator
1515  Vector::end()
1516  {
1517  return (*vector)[0] + local_size();
1518  }
1519 
1520 
1521 
1522  inline Vector::const_iterator
1523  Vector::begin() const
1524  {
1525  return (*vector)[0];
1526  }
1527 
1528 
1529 
1530  inline Vector::const_iterator
1531  Vector::end() const
1532  {
1533  return (*vector)[0] + local_size();
1534  }
1535 
1536 
1537 
1538  inline void
1539  Vector::set(const std::vector<size_type> & indices,
1540  const std::vector<TrilinosScalar> &values)
1541  {
1542  // if we have ghost values, do not allow
1543  // writing to this vector at all.
1544  Assert(!has_ghost_elements(), ExcGhostsPresent());
1545 
1546  Assert(indices.size() == values.size(),
1547  ExcDimensionMismatch(indices.size(), values.size()));
1548 
1549  set(indices.size(), indices.data(), values.data());
1550  }
1551 
1552 
1553 
1554  inline void
1555  Vector::set(const std::vector<size_type> & indices,
1556  const ::Vector<TrilinosScalar> &values)
1557  {
1558  // if we have ghost values, do not allow
1559  // writing to this vector at all.
1560  Assert(!has_ghost_elements(), ExcGhostsPresent());
1561 
1562  Assert(indices.size() == values.size(),
1563  ExcDimensionMismatch(indices.size(), values.size()));
1564 
1565  set(indices.size(), indices.data(), values.begin());
1566  }
1567 
1568 
1569 
1570  inline void
1571  Vector::set(const size_type n_elements,
1572  const size_type * indices,
1573  const TrilinosScalar *values)
1574  {
1575  // if we have ghost values, do not allow
1576  // writing to this vector at all.
1577  Assert(!has_ghost_elements(), ExcGhostsPresent());
1578 
1579  if (last_action == Add)
1580  {
1581  const int ierr = vector->GlobalAssemble(Add);
1582  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1583  }
1584 
1585  if (last_action != Insert)
1586  last_action = Insert;
1587 
1588  for (size_type i = 0; i < n_elements; ++i)
1589  {
1590  const size_type row = indices[i];
1591  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1592  static_cast<TrilinosWrappers::types::int_type>(row));
1593  if (local_row != -1)
1594  (*vector)[0][local_row] = values[i];
1595  else
1596  {
1597  const int ierr = vector->ReplaceGlobalValues(
1598  1,
1599  (const TrilinosWrappers::types::int_type *)(&row),
1600  &values[i]);
1601  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1602  compressed = false;
1603  }
1604  // in set operation, do not use the pre-allocated vector for nonlocal
1605  // entries even if it exists. This is to ensure that we really only
1606  // set the elements touched by the set() method and not all contained
1607  // in the nonlocal entries vector (there is no way to distinguish them
1608  // on the receiving processor)
1609  }
1610  }
1611 
1612 
1613 
1614  inline void
1615  Vector::add(const std::vector<size_type> & indices,
1616  const std::vector<TrilinosScalar> &values)
1617  {
1618  // if we have ghost values, do not allow
1619  // writing to this vector at all.
1620  Assert(!has_ghost_elements(), ExcGhostsPresent());
1621  Assert(indices.size() == values.size(),
1622  ExcDimensionMismatch(indices.size(), values.size()));
1623 
1624  add(indices.size(), indices.data(), values.data());
1625  }
1626 
1627 
1628 
1629  inline void
1630  Vector::add(const std::vector<size_type> & indices,
1631  const ::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.begin());
1640  }
1641 
1642 
1643 
1644  inline void
1645  Vector::add(const size_type n_elements,
1646  const size_type * indices,
1647  const TrilinosScalar *values)
1648  {
1649  // if we have ghost values, do not allow
1650  // writing to this vector at all.
1651  Assert(!has_ghost_elements(), ExcGhostsPresent());
1652 
1653  if (last_action != Add)
1654  {
1655  if (last_action == Insert)
1656  {
1657  const int ierr = vector->GlobalAssemble(Insert);
1658  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1659  }
1660  last_action = Add;
1661  }
1662 
1663  for (size_type i = 0; i < n_elements; ++i)
1664  {
1665  const size_type row = indices[i];
1666  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1667  static_cast<TrilinosWrappers::types::int_type>(row));
1668  if (local_row != -1)
1669  (*vector)[0][local_row] += values[i];
1670  else if (nonlocal_vector.get() == nullptr)
1671  {
1672  const int ierr = vector->SumIntoGlobalValues(
1673  1,
1674  (const TrilinosWrappers::types::int_type *)(&row),
1675  &values[i]);
1676  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1677  compressed = false;
1678  }
1679  else
1680  {
1681  // use pre-allocated vector for non-local entries if it exists for
1682  // addition operation
1683  const TrilinosWrappers::types::int_type my_row =
1684  nonlocal_vector->Map().LID(
1685  static_cast<TrilinosWrappers::types::int_type>(row));
1686  Assert(my_row != -1,
1687  ExcMessage(
1688  "Attempted to write into off-processor vector entry "
1689  "that has not be specified as being writable upon "
1690  "initialization"));
1691  (*nonlocal_vector)[0][my_row] += values[i];
1692  compressed = false;
1693  }
1694  }
1695  }
1696 
1697 
1698 
1699  inline Vector::size_type
1700  Vector::size() const
1701  {
1702 # ifndef DEAL_II_WITH_64BIT_INDICES
1703  return (size_type)(vector->Map().MaxAllGID() + 1 -
1704  vector->Map().MinAllGID());
1705 # else
1706  return (size_type)(vector->Map().MaxAllGID64() + 1 -
1707  vector->Map().MinAllGID64());
1708 # endif
1709  }
1710 
1711 
1712 
1713  inline Vector::size_type
1714  Vector::local_size() const
1715  {
1716  return (size_type)vector->Map().NumMyElements();
1717  }
1718 
1719 
1720 
1721  inline std::pair<Vector::size_type, Vector::size_type>
1722  Vector::local_range() const
1723  {
1724 # ifndef DEAL_II_WITH_64BIT_INDICES
1725  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1726  const TrilinosWrappers::types::int_type end =
1727  vector->Map().MaxMyGID() + 1;
1728 # else
1729  const TrilinosWrappers::types::int_type begin =
1730  vector->Map().MinMyGID64();
1731  const TrilinosWrappers::types::int_type end =
1732  vector->Map().MaxMyGID64() + 1;
1733 # endif
1734 
1735  Assert(
1736  end - begin == vector->Map().NumMyElements(),
1737  ExcMessage(
1738  "This function only makes sense if the elements that this "
1739  "vector stores on the current processor form a contiguous range. "
1740  "This does not appear to be the case for the current vector."));
1741 
1742  return std::make_pair(begin, end);
1743  }
1744 
1745 
1746 
1747  inline TrilinosScalar Vector::operator*(const Vector &vec) const
1748  {
1749  Assert(vector->Map().SameAs(vec.vector->Map()),
1750  ExcDifferentParallelPartitioning());
1751  Assert(!has_ghost_elements(), ExcGhostsPresent());
1752 
1753  TrilinosScalar result;
1754 
1755  const int ierr = vector->Dot(*(vec.vector), &result);
1756  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1757 
1758  return result;
1759  }
1760 
1761 
1762 
1763  inline Vector::real_type
1764  Vector::norm_sqr() const
1765  {
1766  const TrilinosScalar d = l2_norm();
1767  return d * d;
1768  }
1769 
1770 
1771 
1772  inline TrilinosScalar
1773  Vector::mean_value() const
1774  {
1775  Assert(!has_ghost_elements(), ExcGhostsPresent());
1776 
1777  TrilinosScalar mean;
1778  const int ierr = vector->MeanValue(&mean);
1779  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1780 
1781  return mean;
1782  }
1783 
1784 
1785 
1786  inline TrilinosScalar
1787  Vector::min() const
1788  {
1789  TrilinosScalar min_value;
1790  const int ierr = vector->MinValue(&min_value);
1791  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1792 
1793  return min_value;
1794  }
1795 
1796 
1797 
1798  inline TrilinosScalar
1799  Vector::max() const
1800  {
1801  TrilinosScalar max_value;
1802  const int ierr = vector->MaxValue(&max_value);
1803  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1804 
1805  return max_value;
1806  }
1807 
1808 
1809 
1810  inline Vector::real_type
1811  Vector::l1_norm() const
1812  {
1813  Assert(!has_ghost_elements(), ExcGhostsPresent());
1814 
1815  TrilinosScalar d;
1816  const int ierr = vector->Norm1(&d);
1817  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1818 
1819  return d;
1820  }
1821 
1822 
1823 
1824  inline Vector::real_type
1825  Vector::l2_norm() const
1826  {
1827  Assert(!has_ghost_elements(), ExcGhostsPresent());
1828 
1829  TrilinosScalar d;
1830  const int ierr = vector->Norm2(&d);
1831  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1832 
1833  return d;
1834  }
1835 
1836 
1837 
1838  inline Vector::real_type
1839  Vector::lp_norm(const TrilinosScalar p) const
1840  {
1841  Assert(!has_ghost_elements(), ExcGhostsPresent());
1842 
1843  TrilinosScalar norm = 0;
1844  TrilinosScalar sum = 0;
1845  const size_type n_local = local_size();
1846 
1847  // loop over all the elements because
1848  // Trilinos does not support lp norms
1849  for (size_type i = 0; i < n_local; i++)
1850  sum += std::pow(std::fabs((*vector)[0][i]), p);
1851 
1852  norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1853 
1854  return norm;
1855  }
1856 
1857 
1858 
1859  inline Vector::real_type
1860  Vector::linfty_norm() const
1861  {
1862  // while we disallow the other
1863  // norm operations on ghosted
1864  // vectors, this particular norm
1865  // is safe to run even in the
1866  // presence of ghost elements
1867  TrilinosScalar d;
1868  const int ierr = vector->NormInf(&d);
1869  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1870 
1871  return d;
1872  }
1873 
1874 
1875 
1876  inline TrilinosScalar
1877  Vector::add_and_dot(const TrilinosScalar a,
1878  const Vector & V,
1879  const Vector & W)
1880  {
1881  this->add(a, V);
1882  return *this * W;
1883  }
1884 
1885 
1886 
1887  // inline also scalar products, vector
1888  // additions etc. since they are all
1889  // representable by a single Trilinos
1890  // call. This reduces the overhead of the
1891  // wrapper class.
1892  inline Vector &
1893  Vector::operator*=(const TrilinosScalar a)
1894  {
1895  AssertIsFinite(a);
1896 
1897  const int ierr = vector->Scale(a);
1898  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1899 
1900  return *this;
1901  }
1902 
1903 
1904 
1905  inline Vector &
1906  Vector::operator/=(const TrilinosScalar a)
1907  {
1908  AssertIsFinite(a);
1909 
1910  const TrilinosScalar factor = 1. / a;
1911 
1912  AssertIsFinite(factor);
1913 
1914  const int ierr = vector->Scale(factor);
1915  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1916 
1917  return *this;
1918  }
1919 
1920 
1921 
1922  inline Vector &
1923  Vector::operator+=(const Vector &v)
1924  {
1925  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
1926  Assert(vector->Map().SameAs(v.vector->Map()),
1927  ExcDifferentParallelPartitioning());
1928 
1929  const int ierr = vector->Update(1.0, *(v.vector), 1.0);
1930  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1931 
1932  return *this;
1933  }
1934 
1935 
1936 
1937  inline Vector &
1938  Vector::operator-=(const Vector &v)
1939  {
1940  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
1941  Assert(vector->Map().SameAs(v.vector->Map()),
1942  ExcDifferentParallelPartitioning());
1943 
1944  const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
1945  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1946 
1947  return *this;
1948  }
1949 
1950 
1951 
1952  inline void
1953  Vector::add(const TrilinosScalar s)
1954  {
1955  // if we have ghost values, do not allow
1956  // writing to this vector at all.
1957  Assert(!has_ghost_elements(), ExcGhostsPresent());
1958  AssertIsFinite(s);
1959 
1960  size_type n_local = local_size();
1961  for (size_type i = 0; i < n_local; i++)
1962  (*vector)[0][i] += s;
1963  }
1964 
1965 
1966 
1967  inline void
1968  Vector::add(const TrilinosScalar a, const Vector &v)
1969  {
1970  // if we have ghost values, do not allow
1971  // writing to this vector at all.
1972  Assert(!has_ghost_elements(), ExcGhostsPresent());
1973  Assert(local_size() == v.local_size(),
1974  ExcDimensionMismatch(local_size(), v.local_size()));
1975 
1976  AssertIsFinite(a);
1977 
1978  const int ierr = vector->Update(a, *(v.vector), 1.);
1979  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1980  }
1981 
1982 
1983 
1984  inline void
1985  Vector::add(const TrilinosScalar a,
1986  const Vector & v,
1987  const TrilinosScalar b,
1988  const Vector & w)
1989  {
1990  // if we have ghost values, do not allow
1991  // writing to this vector at all.
1992  Assert(!has_ghost_elements(), ExcGhostsPresent());
1993  Assert(local_size() == v.local_size(),
1994  ExcDimensionMismatch(local_size(), v.local_size()));
1995  Assert(local_size() == w.local_size(),
1996  ExcDimensionMismatch(local_size(), w.local_size()));
1997 
1998  AssertIsFinite(a);
1999  AssertIsFinite(b);
2000 
2001  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2002 
2003  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2004  }
2005 
2006 
2007 
2008  inline void
2009  Vector::sadd(const TrilinosScalar s, const Vector &v)
2010  {
2011  // if we have ghost values, do not allow
2012  // writing to this vector at all.
2013  Assert(!has_ghost_elements(), ExcGhostsPresent());
2014  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2015 
2016  AssertIsFinite(s);
2017 
2018  // We assume that the vectors have the same Map
2019  // if the local size is the same and if the vectors are not ghosted
2020  if (local_size() == v.local_size() && !v.has_ghost_elements())
2021  {
2022  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2023  ExcDifferentParallelPartitioning());
2024  const int ierr = vector->Update(1., *(v.vector), s);
2025  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2026  }
2027  else
2028  {
2029  (*this) *= s;
2030  this->add(v, true);
2031  }
2032  }
2033 
2034 
2035 
2036  inline void
2037  Vector::sadd(const TrilinosScalar s,
2038  const TrilinosScalar a,
2039  const Vector & v)
2040  {
2041  // if we have ghost values, do not allow
2042  // writing to this vector at all.
2043  Assert(!has_ghost_elements(), ExcGhostsPresent());
2044  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
2045  AssertIsFinite(s);
2046  AssertIsFinite(a);
2047 
2048  // We assume that the vectors have the same Map
2049  // if the local size is the same and if the vectors are not ghosted
2050  if (local_size() == v.local_size() && !v.has_ghost_elements())
2051  {
2052  Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2053  ExcDifferentParallelPartitioning());
2054  const int ierr = vector->Update(a, *(v.vector), s);
2055  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2056  }
2057  else
2058  {
2059  (*this) *= s;
2060  Vector tmp = v;
2061  tmp *= a;
2062  this->add(tmp, true);
2063  }
2064  }
2065 
2066 
2067 
2068  inline void
2069  Vector::scale(const Vector &factors)
2070  {
2071  // if we have ghost values, do not allow
2072  // writing to this vector at all.
2073  Assert(!has_ghost_elements(), ExcGhostsPresent());
2074  Assert(local_size() == factors.local_size(),
2075  ExcDimensionMismatch(local_size(), factors.local_size()));
2076 
2077  const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2078  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2079  }
2080 
2081 
2082 
2083  inline void
2084  Vector::equ(const TrilinosScalar a, const Vector &v)
2085  {
2086  // if we have ghost values, do not allow
2087  // writing to this vector at all.
2088  Assert(!has_ghost_elements(), ExcGhostsPresent());
2089  AssertIsFinite(a);
2090 
2091  // If we don't have the same map, copy.
2092  if (vector->Map().SameAs(v.vector->Map()) == false)
2093  {
2094  this->sadd(0., a, v);
2095  }
2096  else
2097  {
2098  // Otherwise, just update
2099  int ierr = vector->Update(a, *v.vector, 0.0);
2100  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2101 
2102  last_action = Zero;
2103  }
2104  }
2105 
2106 
2107 
2108  inline const Epetra_MultiVector &
2109  Vector::trilinos_vector() const
2110  {
2111  return static_cast<const Epetra_MultiVector &>(*vector);
2112  }
2113 
2114 
2115 
2116  inline Epetra_FEVector &
2117  Vector::trilinos_vector()
2118  {
2119  return *vector;
2120  }
2121 
2122 
2123 
2124  inline const Epetra_Map &
2125  Vector::vector_partitioner() const
2126  {
2127  return static_cast<const Epetra_Map &>(vector->Map());
2128  }
2129 
2130 
2131 
2132  inline const MPI_Comm &
2133  Vector::get_mpi_communicator() const
2134  {
2135  static MPI_Comm comm;
2136 
2137 # ifdef DEAL_II_WITH_MPI
2138 
2139  const Epetra_MpiComm *mpi_comm =
2140  dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2141  comm = mpi_comm->Comm();
2142 
2143 # else
2144 
2145  comm = MPI_COMM_SELF;
2146 
2147 # endif
2148 
2149  return comm;
2150  }
2151 
2152  template <typename number>
2153  Vector::Vector(const IndexSet & parallel_partitioner,
2154  const ::Vector<number> &v,
2155  const MPI_Comm & communicator)
2156  {
2157  *this =
2158  Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2159  owned_elements = parallel_partitioner;
2160  }
2161 
2162 
2163 
2164  inline Vector &
2165  Vector::operator=(const TrilinosScalar s)
2166  {
2167  AssertIsFinite(s);
2168 
2169  int ierr = vector->PutScalar(s);
2170  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2171 
2172  if (nonlocal_vector.get() != nullptr)
2173  {
2174  ierr = nonlocal_vector->PutScalar(0.);
2175  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2176  }
2177 
2178  return *this;
2179  }
2180  } /* end of namespace MPI */
2181 
2182 # endif /* DOXYGEN */
2183 
2184 } /* end of namespace TrilinosWrappers */
2185 
2189 namespace internal
2190 {
2191  namespace LinearOperatorImplementation
2192  {
2193  template <typename>
2194  class ReinitHelper;
2195 
2200  template <>
2202  {
2203  public:
2204  template <typename Matrix>
2205  static void
2206  reinit_range_vector(const Matrix & matrix,
2208  bool omit_zeroing_entries)
2209  {
2210  v.reinit(matrix.locally_owned_range_indices(),
2211  matrix.get_mpi_communicator(),
2212  omit_zeroing_entries);
2213  }
2214 
2215  template <typename Matrix>
2216  static void
2217  reinit_domain_vector(const Matrix & matrix,
2219  bool omit_zeroing_entries)
2220  {
2221  v.reinit(matrix.locally_owned_domain_indices(),
2222  matrix.get_mpi_communicator(),
2223  omit_zeroing_entries);
2224  }
2225  };
2226 
2227  } // namespace LinearOperatorImplementation
2228 } /* namespace internal */
2229 
2230 
2231 
2237 template <>
2239 {};
2240 
2241 
2242 DEAL_II_NAMESPACE_CLOSE
2243 
2244 #endif // DEAL_II_WITH_TRILINOS
2245 
2246 /*---------------------------- trilinos_vector.h ---------------------------*/
2247 
2248 #endif // dealii_trilinos_vector_h
real_type lp_norm(const real_type p) const
std::unique_ptr< Epetra_MultiVector > nonlocal_vector
static::ExceptionBase & ExcTrilinosError(int arg1)
real_type norm_sqr() const
real_type linfty_norm() const
void sadd(const TrilinosScalar s, const Vector &V)
IndexSet locally_owned_elements() const
real_type l2_norm() const
Number add_and_dot(const Number a, const Vector< Number > &V, const Vector< Number > &W)
TrilinosScalar max() const
Number mean_value() const
void sadd(const Number s, const Vector< Number > &V)
void reinit(const Vector &v, const bool omit_zeroing_entries=false, const bool allow_different_maps=false)
void add(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
size_type size() const
void equ(const Number a, const Vector< Number > &u)
real_type l1_norm() const
void equ(const TrilinosScalar a, const Vector &V)
bool in_local_range(const size_type global_index) const
STL namespace.
std::pair< size_type, size_type > local_range() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:1329
TrilinosScalar min() const
void update_ghost_values() const
iterator end()
Vector< Number > & operator=(const Number s)
LinearAlgebra::distributed::Vector< Number > Vector
static void reinit_range_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
Number operator*(const Vector< Number2 > &V) const
static void reinit_domain_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
real_type linfty_norm() const
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
static::ExceptionBase & ExcMessage(std::string arg1)
numbers::NumberTraits< Number >::real_type real_type
Definition: vector.h:133
void scale(const Vector< Number > &scaling_factors)
void set(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:408
Vector< Number > & operator+=(const Vector< Number > &V)
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1227
real_type l2_norm() const
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::unique_ptr< Epetra_FEVector > vector
void compress(::VectorOperation::values operation)
#define DeclException0(Exception0)
Definition: exceptions.h:385
Number operator[](const size_type i) const
IndexSet locally_owned_elements() const
real_type norm_sqr() const
Vector< Number > & operator-=(const Vector< Number > &V)
const Epetra_Map & vector_partitioner() const
Vector< Number > & operator*=(const Number factor)
const MPI_Comm & get_mpi_communicator() const
friend class internal::VectorReference
void import_nonlocal_data_for_fe(const ::TrilinosWrappers::SparseMatrix &matrix, const Vector &vector)
iterator begin()
real_type l1_norm() const
TrilinosScalar mean_value() const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:523
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
static::ExceptionBase & ExcGhostsPresent()
std::size_t memory_consumption() const
Vector< Number > & operator/=(const Number factor)
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< TrilinosScalar > &values) const
real_type lp_norm(const TrilinosScalar p) const
bool in_local_range(const size_type index) const
void scale(const Vector &scaling_factors)
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< OtherNumber > &values) const
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:444
size_type local_size() const
double norm(const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim >>>> &Du)
Definition: divergence.h:553
static::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
TrilinosScalar add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W)
static::ExceptionBase & ExcDifferentParallelPartitioning()
Number operator()(const size_type i) const
static::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertIsFinite(number)
Definition: exceptions.h:1423
const Epetra_MultiVector & trilinos_vector() const