Reference documentation for deal.II version Git 091d4db 2018-04-24 17:37:29 -0600
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 at
12 // the top level of the deal.II distribution.
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/subscriptor.h>
24 # include <deal.II/base/index_set.h>
25 # include <deal.II/base/utilities.h>
26 # include <deal.II/base/mpi.h>
27 # include <deal.II/lac/exceptions.h>
28 # include <deal.II/lac/vector.h>
29 # include <deal.II/lac/vector_operation.h>
30 # include <deal.II/lac/vector_type_traits.h>
31 
32 # include <vector>
33 # include <utility>
34 # include <memory>
35 
36 # include <Epetra_ConfigDefs.h>
37 # ifdef DEAL_II_WITH_MPI // only if MPI is installed
38 # include <mpi.h>
39 # include <Epetra_MpiComm.h>
40 # else
41 # include <Epetra_SerialComm.h>
42 # endif
43 # include <Epetra_FEVector.h>
44 # include <Epetra_Map.h>
45 # include <Epetra_LocalMap.h>
46 
47 DEAL_II_NAMESPACE_OPEN
48 
49 
60 namespace TrilinosWrappers
61 {
62  class SparseMatrix;
63 
74  namespace internal
75  {
79  typedef ::types::global_dof_index size_type;
80 
90  class VectorReference
91  {
92  private:
97  VectorReference (MPI::Vector &vector,
98  const size_type index);
99 
100  public:
101 
113  const VectorReference &
114  operator= (const VectorReference &r) const;
115 
119  VectorReference &
120  operator= (const VectorReference &r);
121 
125  const VectorReference &
126  operator= (const TrilinosScalar &s) const;
127 
131  const VectorReference &
132  operator+= (const TrilinosScalar &s) const;
133 
137  const VectorReference &
138  operator-= (const TrilinosScalar &s) const;
139 
143  const VectorReference &
144  operator*= (const TrilinosScalar &s) const;
145 
149  const VectorReference &
150  operator/= (const TrilinosScalar &s) const;
151 
156  operator TrilinosScalar () const;
157 
162  int,
163  << "An error with error number " << arg1
164  << " occurred while calling a Trilinos function");
165 
166  private:
170  MPI::Vector &vector;
171 
175  const size_type index;
176 
181  friend class ::TrilinosWrappers::MPI::Vector;
182  };
183  }
188  namespace
189  {
190 #ifndef DEAL_II_WITH_64BIT_INDICES
191  // define a helper function that queries the global ID of local ID of
192  // an Epetra_BlockMap object by calling either the 32- or 64-bit
193  // function necessary.
194  inline
195  int gid(const Epetra_BlockMap &map, int i)
196  {
197  return map.GID(i);
198  }
199 #else
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
204  long long int gid(const Epetra_BlockMap &map, int i)
205  {
206  return map.GID64(i);
207  }
208 #endif
209  }
210 
217  namespace MPI
218  {
219  class BlockVector;
220 
397  class Vector : public Subscriptor
398  {
399  public:
405  typedef TrilinosScalar value_type;
406  typedef TrilinosScalar real_type;
407  typedef ::types::global_dof_index size_type;
408  typedef value_type *iterator;
409  typedef const value_type *const_iterator;
410  typedef internal::VectorReference reference;
411  typedef const internal::VectorReference const_reference;
412 
422  Vector ();
423 
427  Vector (const Vector &v);
428 
447  explicit Vector (const IndexSet &parallel_partitioning,
448  const MPI_Comm &communicator = MPI_COMM_WORLD);
449 
461  Vector (const IndexSet &local,
462  const IndexSet &ghost,
463  const MPI_Comm &communicator = MPI_COMM_WORLD);
464 
479  Vector (const IndexSet &parallel_partitioning,
480  const Vector &v,
481  const MPI_Comm &communicator = MPI_COMM_WORLD);
482 
495  template <typename Number>
496  Vector (const IndexSet &parallel_partitioning,
497  const ::Vector<Number> &v,
498  const MPI_Comm &communicator = MPI_COMM_WORLD);
499 
504  Vector (Vector &&v) noexcept;
505 
509  ~Vector () = default;
510 
515  void clear ();
516 
540  void reinit (const Vector &v,
541  const bool omit_zeroing_entries = false,
542  const bool allow_different_maps = false);
543 
566  void reinit (const IndexSet &parallel_partitioning,
567  const MPI_Comm &communicator = MPI_COMM_WORLD,
568  const bool omit_zeroing_entries = false);
569 
595  void 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 reinit (const BlockVector &v,
604  const bool import_data = false);
605 
622  void compress (::VectorOperation::values operation);
623 
636  Vector &operator= (const TrilinosScalar s);
637 
643  Vector &operator= (const Vector &v);
644 
649  Vector &operator= (Vector &&v) noexcept;
650 
658  template <typename Number>
659  Vector &operator= (const ::Vector<Number> &v);
660 
679  (const ::TrilinosWrappers::SparseMatrix &matrix,
680  const Vector &vector);
681 
687  bool operator== (const Vector &v) const;
688 
694  bool operator!= (const Vector &v) const;
695 
699  size_type size () const;
700 
712  size_type local_size () const;
713 
735  std::pair<size_type, size_type> local_range () const;
736 
744  bool in_local_range (const size_type index) const;
745 
760 
768  bool has_ghost_elements() const;
769 
776  void update_ghost_values () const;
777 
782  TrilinosScalar operator* (const Vector &vec) const;
783 
787  real_type norm_sqr () const;
788 
792  TrilinosScalar mean_value () const;
793 
797  TrilinosScalar min () const;
798 
802  TrilinosScalar max () const;
803 
807  real_type l1_norm () const;
808 
813  real_type l2_norm () const;
814 
819  real_type lp_norm (const TrilinosScalar p) const;
820 
824  real_type linfty_norm () const;
825 
844  TrilinosScalar add_and_dot (const TrilinosScalar a,
845  const Vector &V,
846  const Vector &W);
847 
853  bool all_zero () const;
854 
860  bool is_non_negative () const;
862 
863 
868 
876  reference
877  operator() (const size_type index);
878 
886  TrilinosScalar
887  operator() (const size_type index) const;
888 
894  reference
895  operator[] (const size_type index);
896 
902  TrilinosScalar
903  operator[] (const size_type index) const;
904 
920  void extract_subvector_to (const std::vector<size_type> &indices,
921  std::vector<TrilinosScalar> &values) const;
922 
950  template <typename ForwardIterator, typename OutputIterator>
951  void extract_subvector_to (ForwardIterator indices_begin,
952  const ForwardIterator indices_end,
953  OutputIterator values_begin) const;
954 
965  iterator begin ();
966 
971  const_iterator begin () const;
972 
977  iterator end ();
978 
983  const_iterator end () const;
984 
986 
987 
992 
999  void set (const std::vector<size_type> &indices,
1000  const std::vector<TrilinosScalar> &values);
1001 
1006  void set (const std::vector<size_type> &indices,
1007  const ::Vector<TrilinosScalar> &values);
1008 
1014  void set (const size_type n_elements,
1015  const size_type *indices,
1016  const TrilinosScalar *values);
1017 
1022  void add (const std::vector<size_type> &indices,
1023  const std::vector<TrilinosScalar> &values);
1024 
1029  void add (const std::vector<size_type> &indices,
1030  const ::Vector<TrilinosScalar> &values);
1031 
1037  void add (const size_type n_elements,
1038  const size_type *indices,
1039  const TrilinosScalar *values);
1040 
1044  Vector &operator*= (const TrilinosScalar factor);
1045 
1049  Vector &operator/= (const TrilinosScalar factor);
1050 
1054  Vector &operator+= (const Vector &V);
1055 
1059  Vector &operator-= (const Vector &V);
1060 
1065  void add (const TrilinosScalar s);
1066 
1079  void add (const Vector &V,
1080  const bool allow_different_maps = false);
1081 
1085  void add (const TrilinosScalar a,
1086  const Vector &V);
1087 
1091  void add (const TrilinosScalar a,
1092  const Vector &V,
1093  const TrilinosScalar b,
1094  const Vector &W);
1095 
1100  void sadd (const TrilinosScalar s,
1101  const Vector &V);
1102 
1106  void sadd (const TrilinosScalar s,
1107  const TrilinosScalar a,
1108  const Vector &V);
1109 
1115  void scale (const Vector &scaling_factors);
1116 
1120  void equ (const TrilinosScalar a,
1121  const Vector &V);
1123 
1128 
1133  const Epetra_MultiVector &trilinos_vector () const;
1134 
1139  Epetra_FEVector &trilinos_vector ();
1140 
1145  const Epetra_Map &vector_partitioner () const;
1146 
1154  void print (std::ostream &out,
1155  const unsigned int precision = 3,
1156  const bool scientific = true,
1157  const bool across = true) const;
1158 
1172  void swap (Vector &v);
1173 
1177  std::size_t memory_consumption () const;
1178 
1183  const MPI_Comm &get_mpi_communicator () const;
1185 
1190 
1195  int,
1196  << "An error with error number " << arg1
1197  << " occurred while calling a Trilinos function");
1198 
1203  size_type, size_type, size_type, size_type,
1204  << "You tried to access element " << arg1
1205  << " of a distributed vector, but this element is not stored "
1206  << "on the current processor. Note: There are "
1207  << arg2 << " elements stored "
1208  << "on the current processor from within the range "
1209  << arg3 << " through " << arg4
1210  << " but Trilinos vectors need not store contiguous "
1211  << "ranges on each processor, and not every element in "
1212  << "this range may in fact be stored locally.");
1213 
1214  private:
1226  Epetra_CombineMode last_action;
1227 
1233 
1239 
1245  std::unique_ptr<Epetra_FEVector> vector;
1246 
1252  std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1253 
1258 
1262  friend class internal::VectorReference;
1263  };
1264 
1265 
1266 
1267 
1268 // ------------------- inline and template functions --------------
1269 
1270 
1279  inline
1280  void swap (Vector &u, Vector &v)
1281  {
1282  u.swap (v);
1283  }
1284  }
1285 
1286 #ifndef DOXYGEN
1287 
1288  namespace internal
1289  {
1290  inline
1291  VectorReference::VectorReference (MPI::Vector &vector,
1292  const size_type index)
1293  :
1294  vector (vector),
1295  index (index)
1296  {}
1297 
1298 
1299  inline
1300  const VectorReference &
1301  VectorReference::operator= (const VectorReference &r) const
1302  {
1303  // as explained in the class
1304  // documentation, this is not the copy
1305  // operator. so simply pass on to the
1306  // "correct" assignment operator
1307  *this = static_cast<TrilinosScalar> (r);
1308 
1309  return *this;
1310  }
1311 
1312 
1313 
1314  inline
1315  VectorReference &
1316  VectorReference::operator= (const VectorReference &r)
1317  {
1318  // as above
1319  *this = static_cast<TrilinosScalar> (r);
1320 
1321  return *this;
1322  }
1323 
1324 
1325  inline
1326  const VectorReference &
1327  VectorReference::operator= (const TrilinosScalar &value) const
1328  {
1329  vector.set (1, &index, &value);
1330  return *this;
1331  }
1332 
1333 
1334 
1335  inline
1336  const VectorReference &
1337  VectorReference::operator+= (const TrilinosScalar &value) const
1338  {
1339  vector.add (1, &index, &value);
1340  return *this;
1341  }
1342 
1343 
1344 
1345  inline
1346  const VectorReference &
1347  VectorReference::operator-= (const TrilinosScalar &value) const
1348  {
1349  TrilinosScalar new_value = -value;
1350  vector.add (1, &index, &new_value);
1351  return *this;
1352  }
1353 
1354 
1355 
1356  inline
1357  const VectorReference &
1358  VectorReference::operator*= (const TrilinosScalar &value) const
1359  {
1360  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1361  vector.set (1, &index, &new_value);
1362  return *this;
1363  }
1364 
1365 
1366 
1367  inline
1368  const VectorReference &
1369  VectorReference::operator/= (const TrilinosScalar &value) const
1370  {
1371  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1372  vector.set (1, &index, &new_value);
1373  return *this;
1374  }
1375  }
1376 
1377  namespace MPI
1378  {
1379  inline
1380  bool
1381  Vector::in_local_range (const size_type index) const
1382  {
1383  std::pair<size_type, size_type> range = local_range();
1384 
1385  return ((index >= range.first) && (index < range.second));
1386  }
1387 
1388 
1389 
1390  inline
1391  IndexSet
1393  {
1394  Assert(owned_elements.size()==size(),
1395  ExcMessage("The locally owned elements have not been properly initialized!"
1396  " This happens for example if this object has been initialized"
1397  " with exactly one overlapping IndexSet."));
1398  return owned_elements;
1399  }
1400 
1401 
1402 
1403  inline
1404  bool
1405  Vector::has_ghost_elements() const
1406  {
1407  return has_ghosts;
1408  }
1409 
1410 
1411 
1412  inline
1413  void
1415  {}
1416 
1417 
1418 
1419  inline
1420  internal::VectorReference
1421  Vector::operator() (const size_type index)
1422  {
1423  return internal::VectorReference (*this, index);
1424  }
1425 
1426 
1427 
1428  inline
1429  internal::VectorReference
1430  Vector::operator[] (const size_type index)
1431  {
1432  return operator() (index);
1433  }
1434 
1435 
1436 
1437  inline
1438  TrilinosScalar
1439  Vector::operator[] (const size_type index) const
1440  {
1441  return operator() (index);
1442  }
1443 
1444 
1445 
1446  inline
1447  void Vector::extract_subvector_to (const std::vector<size_type> &indices,
1448  std::vector<TrilinosScalar> &values) const
1449  {
1450  for (size_type i = 0; i < indices.size(); ++i)
1451  values[i] = operator()(indices[i]);
1452  }
1453 
1454 
1455 
1456  template <typename ForwardIterator, typename OutputIterator>
1457  inline
1458  void Vector::extract_subvector_to (ForwardIterator indices_begin,
1459  const ForwardIterator indices_end,
1460  OutputIterator values_begin) const
1461  {
1462  while (indices_begin != indices_end)
1463  {
1464  *values_begin = operator()(*indices_begin);
1465  indices_begin++;
1466  values_begin++;
1467  }
1468  }
1469 
1470 
1471 
1472  inline
1473  Vector::iterator
1474  Vector::begin()
1475  {
1476  return (*vector)[0];
1477  }
1478 
1479 
1480 
1481  inline
1482  Vector::iterator
1483  Vector::end()
1484  {
1485  return (*vector)[0]+local_size();
1486  }
1487 
1488 
1489 
1490  inline
1491  Vector::const_iterator
1492  Vector::begin() const
1493  {
1494  return (*vector)[0];
1495  }
1496 
1497 
1498 
1499  inline
1500  Vector::const_iterator
1501  Vector::end() const
1502  {
1503  return (*vector)[0]+local_size();
1504  }
1505 
1506 
1507 
1508  inline
1509  void
1510  Vector::set (const std::vector<size_type> &indices,
1511  const std::vector<TrilinosScalar> &values)
1512  {
1513  // if we have ghost values, do not allow
1514  // writing to this vector at all.
1515  Assert (!has_ghost_elements(), ExcGhostsPresent());
1516 
1517  Assert (indices.size() == values.size(),
1518  ExcDimensionMismatch(indices.size(),values.size()));
1519 
1520  set (indices.size(), indices.data(), values.data());
1521  }
1522 
1523 
1524 
1525  inline
1526  void
1527  Vector::set (const std::vector<size_type> &indices,
1528  const ::Vector<TrilinosScalar> &values)
1529  {
1530  // if we have ghost values, do not allow
1531  // writing to this vector at all.
1532  Assert (!has_ghost_elements(), ExcGhostsPresent());
1533 
1534  Assert (indices.size() == values.size(),
1535  ExcDimensionMismatch(indices.size(),values.size()));
1536 
1537  set (indices.size(), indices.data(), values.begin());
1538  }
1539 
1540 
1541 
1542  inline
1543  void
1544  Vector::set (const size_type n_elements,
1545  const size_type *indices,
1546  const TrilinosScalar *values)
1547  {
1548  // if we have ghost values, do not allow
1549  // writing to this vector at all.
1550  Assert (!has_ghost_elements(), ExcGhostsPresent());
1551 
1552  if (last_action == Add)
1553  {
1554  const int ierr = vector->GlobalAssemble(Add);
1555  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1556  }
1557 
1558  if (last_action != Insert)
1559  last_action = Insert;
1560 
1561  for (size_type i=0; i<n_elements; ++i)
1562  {
1563  const size_type row = indices[i];
1564  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1565  if (local_row != -1)
1566  (*vector)[0][local_row] = values[i];
1567  else
1568  {
1569  const int ierr = vector->ReplaceGlobalValues (1,
1570  (const TrilinosWrappers::types::int_type *)(&row),
1571  &values[i]);
1572  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1573  compressed = false;
1574  }
1575  // in set operation, do not use the pre-allocated vector for nonlocal
1576  // entries even if it exists. This is to ensure that we really only
1577  // set the elements touched by the set() method and not all contained
1578  // in the nonlocal entries vector (there is no way to distinguish them
1579  // on the receiving processor)
1580  }
1581  }
1582 
1583 
1584 
1585  inline
1586  void
1587  Vector::add (const std::vector<size_type> &indices,
1588  const std::vector<TrilinosScalar> &values)
1589  {
1590  // if we have ghost values, do not allow
1591  // writing to this vector at all.
1592  Assert (!has_ghost_elements(), ExcGhostsPresent());
1593  Assert (indices.size() == values.size(),
1594  ExcDimensionMismatch(indices.size(),values.size()));
1595 
1596  add (indices.size(), indices.data(), values.data());
1597  }
1598 
1599 
1600 
1601  inline
1602  void
1603  Vector::add (const std::vector<size_type> &indices,
1604  const ::Vector<TrilinosScalar> &values)
1605  {
1606  // if we have ghost values, do not allow
1607  // writing to this vector at all.
1608  Assert (!has_ghost_elements(), ExcGhostsPresent());
1609  Assert (indices.size() == values.size(),
1610  ExcDimensionMismatch(indices.size(),values.size()));
1611 
1612  add (indices.size(), indices.data(), values.begin());
1613  }
1614 
1615 
1616 
1617  inline
1618  void
1619  Vector::add (const size_type n_elements,
1620  const size_type *indices,
1621  const TrilinosScalar *values)
1622  {
1623  // if we have ghost values, do not allow
1624  // writing to this vector at all.
1625  Assert (!has_ghost_elements(), ExcGhostsPresent());
1626 
1627  if (last_action != Add)
1628  {
1629  if (last_action == Insert)
1630  {
1631  const int ierr = vector->GlobalAssemble(Insert);
1632  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1633  }
1634  last_action = Add;
1635  }
1636 
1637  for (size_type i=0; i<n_elements; ++i)
1638  {
1639  const size_type row = indices[i];
1640  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1641  if (local_row != -1)
1642  (*vector)[0][local_row] += values[i];
1643  else if (nonlocal_vector.get() == nullptr)
1644  {
1645  const int ierr = vector->SumIntoGlobalValues (1,
1646  (const TrilinosWrappers::types::int_type *)(&row),
1647  &values[i]);
1648  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1649  compressed = false;
1650  }
1651  else
1652  {
1653  // use pre-allocated vector for non-local entries if it exists for
1654  // addition operation
1655  const TrilinosWrappers::types::int_type my_row = nonlocal_vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1656  Assert(my_row != -1,
1657  ExcMessage("Attempted to write into off-processor vector entry "
1658  "that has not be specified as being writable upon "
1659  "initialization"));
1660  (*nonlocal_vector)[0][my_row] += values[i];
1661  compressed = false;
1662  }
1663  }
1664  }
1665 
1666 
1667 
1668  inline
1669  Vector::size_type
1670  Vector::size () const
1671  {
1672 #ifndef DEAL_II_WITH_64BIT_INDICES
1673  return (size_type) (vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID());
1674 #else
1675  return (size_type) (vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64());
1676 #endif
1677  }
1678 
1679 
1680 
1681  inline
1682  Vector::size_type
1683  Vector::local_size () const
1684  {
1685  return (size_type) vector->Map().NumMyElements();
1686  }
1687 
1688 
1689 
1690  inline
1691  std::pair<Vector::size_type, Vector::size_type>
1692  Vector::local_range () const
1693  {
1694 #ifndef DEAL_II_WITH_64BIT_INDICES
1695  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1696  const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID()+1;
1697 #else
1698  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID64();
1699  const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID64()+1;
1700 #endif
1701 
1702  Assert (end-begin == vector->Map().NumMyElements(),
1703  ExcMessage ("This function only makes sense if the elements that this "
1704  "vector stores on the current processor form a contiguous range. "
1705  "This does not appear to be the case for the current vector."));
1706 
1707  return std::make_pair (begin, end);
1708  }
1709 
1710 
1711 
1712  inline
1713  TrilinosScalar
1714  Vector::operator* (const Vector &vec) const
1715  {
1716  Assert (vector->Map().SameAs(vec.vector->Map()),
1717  ExcDifferentParallelPartitioning());
1718  Assert (!has_ghost_elements(), ExcGhostsPresent());
1719 
1720  TrilinosScalar result;
1721 
1722  const int ierr = vector->Dot(*(vec.vector), &result);
1723  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1724 
1725  return result;
1726  }
1727 
1728 
1729 
1730  inline
1732  Vector::norm_sqr () const
1733  {
1734  const TrilinosScalar d = l2_norm();
1735  return d*d;
1736  }
1737 
1738 
1739 
1740  inline
1741  TrilinosScalar
1742  Vector::mean_value () const
1743  {
1744  Assert (!has_ghost_elements(), ExcGhostsPresent());
1745 
1746  TrilinosScalar mean;
1747  const int ierr = vector->MeanValue (&mean);
1748  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1749 
1750  return mean;
1751  }
1752 
1753 
1754 
1755  inline
1756  TrilinosScalar
1757  Vector::min () const
1758  {
1759  TrilinosScalar min_value;
1760  const int ierr = vector->MinValue (&min_value);
1761  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1762 
1763  return min_value;
1764  }
1765 
1766 
1767 
1768  inline
1769  TrilinosScalar
1770  Vector::max () const
1771  {
1772  TrilinosScalar max_value;
1773  const int ierr = vector->MaxValue (&max_value);
1774  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1775 
1776  return max_value;
1777  }
1778 
1779 
1780 
1781  inline
1783  Vector::l1_norm () const
1784  {
1785  Assert (!has_ghost_elements(), ExcGhostsPresent());
1786 
1787  TrilinosScalar d;
1788  const int ierr = vector->Norm1 (&d);
1789  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1790 
1791  return d;
1792  }
1793 
1794 
1795 
1796  inline
1798  Vector::l2_norm () const
1799  {
1800  Assert (!has_ghost_elements(), ExcGhostsPresent());
1801 
1802  TrilinosScalar d;
1803  const int ierr = vector->Norm2 (&d);
1804  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1805 
1806  return d;
1807  }
1808 
1809 
1810 
1811  inline
1813  Vector::lp_norm (const TrilinosScalar p) const
1814  {
1815  Assert (!has_ghost_elements(), ExcGhostsPresent());
1816 
1817  TrilinosScalar norm = 0;
1818  TrilinosScalar sum=0;
1819  const size_type n_local = local_size();
1820 
1821  // loop over all the elements because
1822  // Trilinos does not support lp norms
1823  for (size_type i=0; i<n_local; i++)
1824  sum += std::pow(std::fabs((*vector)[0][i]), p);
1825 
1826  norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
1827 
1828  return norm;
1829  }
1830 
1831 
1832 
1833  inline
1835  Vector::linfty_norm () const
1836  {
1837  // while we disallow the other
1838  // norm operations on ghosted
1839  // vectors, this particular norm
1840  // is safe to run even in the
1841  // presence of ghost elements
1842  TrilinosScalar d;
1843  const int ierr = vector->NormInf (&d);
1844  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1845 
1846  return d;
1847  }
1848 
1849 
1850 
1851  inline
1852  TrilinosScalar
1853  Vector::add_and_dot (const TrilinosScalar a,
1854  const Vector &V,
1855  const Vector &W)
1856  {
1857  this->add(a, V);
1858  return *this * W;
1859  }
1860 
1861 
1862 
1863  // inline also scalar products, vector
1864  // additions etc. since they are all
1865  // representable by a single Trilinos
1866  // call. This reduces the overhead of the
1867  // wrapper class.
1868  inline
1869  Vector &
1870  Vector::operator*= (const TrilinosScalar a)
1871  {
1872  AssertIsFinite(a);
1873 
1874  const int ierr = vector->Scale(a);
1875  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1876 
1877  return *this;
1878  }
1879 
1880 
1881 
1882  inline
1883  Vector &
1884  Vector::operator/= (const TrilinosScalar a)
1885  {
1886  AssertIsFinite(a);
1887 
1888  const TrilinosScalar factor = 1./a;
1889 
1890  AssertIsFinite(factor);
1891 
1892  const int ierr = vector->Scale(factor);
1893  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1894 
1895  return *this;
1896  }
1897 
1898 
1899 
1900  inline
1901  Vector &
1902  Vector::operator+= (const Vector &v)
1903  {
1904  Assert (size() == v.size(),
1905  ExcDimensionMismatch(size(), v.size()));
1906  Assert (vector->Map().SameAs(v.vector->Map()),
1907  ExcDifferentParallelPartitioning());
1908 
1909  const int ierr = vector->Update (1.0, *(v.vector), 1.0);
1910  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1911 
1912  return *this;
1913  }
1914 
1915 
1916 
1917  inline
1918  Vector &
1919  Vector::operator-= (const Vector &v)
1920  {
1921  Assert (size() == v.size(),
1922  ExcDimensionMismatch(size(), v.size()));
1923  Assert (vector->Map().SameAs(v.vector->Map()),
1924  ExcDifferentParallelPartitioning());
1925 
1926  const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
1927  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1928 
1929  return *this;
1930  }
1931 
1932 
1933 
1934  inline
1935  void
1936  Vector::add (const TrilinosScalar s)
1937  {
1938  // if we have ghost values, do not allow
1939  // writing to this vector at all.
1940  Assert (!has_ghost_elements(), ExcGhostsPresent());
1941  AssertIsFinite(s);
1942 
1943  size_type n_local = local_size();
1944  for (size_type i=0; i<n_local; i++)
1945  (*vector)[0][i] += s;
1946  }
1947 
1948 
1949 
1950  inline
1951  void
1952  Vector::add (const TrilinosScalar a,
1953  const Vector &v)
1954  {
1955  // if we have ghost values, do not allow
1956  // writing to this vector at all.
1957  Assert (!has_ghost_elements(), ExcGhostsPresent());
1958  Assert (local_size() == v.local_size(),
1959  ExcDimensionMismatch(local_size(), v.local_size()));
1960 
1961  AssertIsFinite(a);
1962 
1963  const int ierr = vector->Update(a, *(v.vector), 1.);
1964  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1965  }
1966 
1967 
1968 
1969  inline
1970  void
1971  Vector::add (const TrilinosScalar a,
1972  const Vector &v,
1973  const TrilinosScalar b,
1974  const Vector &w)
1975  {
1976  // if we have ghost values, do not allow
1977  // writing to this vector at all.
1978  Assert (!has_ghost_elements(), ExcGhostsPresent());
1979  Assert (local_size() == v.local_size(),
1980  ExcDimensionMismatch(local_size(), v.local_size()));
1981  Assert (local_size() == w.local_size(),
1982  ExcDimensionMismatch(local_size(), w.local_size()));
1983 
1984  AssertIsFinite(a);
1985  AssertIsFinite(b);
1986 
1987  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
1988 
1989  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1990  }
1991 
1992 
1993 
1994  inline
1995  void
1996  Vector::sadd (const TrilinosScalar s,
1997  const Vector &v)
1998  {
1999  // if we have ghost values, do not allow
2000  // writing to this vector at all.
2001  Assert (!has_ghost_elements(), ExcGhostsPresent());
2002  Assert (size() == v.size(),
2003  ExcDimensionMismatch (size(), v.size()));
2004 
2005  AssertIsFinite(s);
2006 
2007  // We assume that the vectors have the same Map
2008  // if the local size is the same and if the vectors are not ghosted
2009  if (local_size() == v.local_size() && !v.has_ghost_elements())
2010  {
2011  Assert (this->vector->Map().SameAs(v.vector->Map())==true,
2012  ExcDifferentParallelPartitioning());
2013  const int ierr = vector->Update(1., *(v.vector), s);
2014  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2015  }
2016  else
2017  {
2018  (*this) *= s;
2019  this->add(v, true);
2020  }
2021  }
2022 
2023 
2024 
2025  inline
2026  void
2027  Vector::sadd (const TrilinosScalar s,
2028  const TrilinosScalar a,
2029  const Vector &v)
2030  {
2031  // if we have ghost values, do not allow
2032  // writing to this vector at all.
2033  Assert (!has_ghost_elements(), ExcGhostsPresent());
2034  Assert (size() == v.size(),
2035  ExcDimensionMismatch (size(), v.size()));
2036  AssertIsFinite(s);
2037  AssertIsFinite(a);
2038 
2039  // We assume that the vectors have the same Map
2040  // if the local size is the same and if the vectors are not ghosted
2041  if (local_size() == v.local_size() && !v.has_ghost_elements())
2042  {
2043  Assert (this->vector->Map().SameAs(v.vector->Map())==true,
2044  ExcDifferentParallelPartitioning());
2045  const int ierr = vector->Update(a, *(v.vector), s);
2046  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2047  }
2048  else
2049  {
2050  (*this) *= s;
2051  Vector tmp = v;
2052  tmp *= a;
2053  this->add(tmp, true);
2054  }
2055  }
2056 
2057 
2058 
2059  inline
2060  void
2061  Vector::scale (const Vector &factors)
2062  {
2063  // if we have ghost values, do not allow
2064  // writing to this vector at all.
2065  Assert (!has_ghost_elements(), ExcGhostsPresent());
2066  Assert (local_size() == factors.local_size(),
2067  ExcDimensionMismatch(local_size(), factors.local_size()));
2068 
2069  const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
2070  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2071  }
2072 
2073 
2074 
2075  inline
2076  void
2077  Vector::equ (const TrilinosScalar a,
2078  const Vector &v)
2079  {
2080  // if we have ghost values, do not allow
2081  // writing to this vector at all.
2082  Assert (!has_ghost_elements(), ExcGhostsPresent());
2083  AssertIsFinite(a);
2084 
2085  // If we don't have the same map, copy.
2086  if (vector->Map().SameAs(v.vector->Map())==false)
2087  {
2088  this->sadd(0., a, v);
2089  }
2090  else
2091  {
2092  // Otherwise, just update
2093  int ierr = vector->Update(a, *v.vector, 0.0);
2094  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2095 
2096  last_action = Zero;
2097  }
2098 
2099  }
2100 
2101 
2102 
2103  inline
2104  const Epetra_MultiVector &
2105  Vector::trilinos_vector () const
2106  {
2107  return static_cast<const Epetra_MultiVector &>(*vector);
2108  }
2109 
2110 
2111 
2112  inline
2113  Epetra_FEVector &
2114  Vector::trilinos_vector ()
2115  {
2116  return *vector;
2117  }
2118 
2119 
2120 
2121  inline
2122  const Epetra_Map &
2123  Vector::vector_partitioner () const
2124  {
2125  return static_cast<const Epetra_Map &>(vector->Map());
2126  }
2127 
2128 
2129 
2130  inline
2131  const MPI_Comm &
2132  Vector::get_mpi_communicator () const
2133  {
2134  static MPI_Comm comm;
2135 
2136 #ifdef DEAL_II_WITH_MPI
2137 
2138  const Epetra_MpiComm *mpi_comm
2139  = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2140  comm = mpi_comm->Comm();
2141 
2142 #else
2143 
2144  comm = MPI_COMM_SELF;
2145 
2146 #endif
2147 
2148  return comm;
2149  }
2150 
2151  template <typename number>
2152  Vector::Vector (const IndexSet &parallel_partitioner,
2153  const ::Vector<number> &v,
2154  const MPI_Comm &communicator)
2155  {
2156  *this = Vector(parallel_partitioner.make_trilinos_map (communicator, true),
2157  v);
2158  owned_elements = parallel_partitioner;
2159  }
2160 
2161 
2162 
2163  inline
2164  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> class ReinitHelper;
2194 
2199  template <>
2200  class ReinitHelper<TrilinosWrappers::MPI::Vector>
2201  {
2202  public:
2203  template <typename Matrix>
2204  static
2205  void reinit_range_vector (const Matrix &matrix,
2207  bool omit_zeroing_entries)
2208  {
2209  v.reinit(matrix.locally_owned_range_indices(), matrix.get_mpi_communicator(), omit_zeroing_entries);
2210  }
2211 
2212  template <typename Matrix>
2213  static
2214  void reinit_domain_vector(const Matrix &matrix,
2216  bool omit_zeroing_entries)
2217  {
2218  v.reinit(matrix.locally_owned_domain_indices(), matrix.get_mpi_communicator(), omit_zeroing_entries);
2219  }
2220  };
2221 
2222  } /* namespace LinearOperator */
2223 } /* namespace internal */
2224 
2225 
2226 
2232 template <>
2233 struct is_serial_vector< TrilinosWrappers::MPI::Vector > : std::false_type
2234 {
2235 };
2236 
2237 
2238 DEAL_II_NAMESPACE_CLOSE
2239 
2240 #endif // DEAL_II_WITH_TRILINOS
2241 
2242 /*---------------------------- trilinos_vector.h ---------------------------*/
2243 
2244 #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:1193
TrilinosScalar min() const
void update_ghost_values() const
iterator end()
Vector< Number > & operator=(const Number s)
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:128
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:344
Vector< Number > & operator+=(const Vector< Number > &V)
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1114
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:321
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:549
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
double norm(const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Du)
Definition: divergence.h:534
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:380
size_type local_size() const
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:1271
const Epetra_MultiVector & trilinos_vector() const