Reference documentation for deal.II version Git f264146 2017-10-18 16:34:34 -0400
trilinos_vector.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2017 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_type_traits.h>
30 
31 # include <vector>
32 # include <utility>
33 # include <memory>
34 
35 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
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 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
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,
99  const size_type index);
100 
101  public:
102 
114  const VectorReference &
115  operator= (const VectorReference &r) const;
116 
120  VectorReference &
121  operator= (const VectorReference &r);
122 
126  const VectorReference &
127  operator= (const TrilinosScalar &s) const;
128 
132  const VectorReference &
133  operator+= (const TrilinosScalar &s) const;
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 
157  operator TrilinosScalar () const;
158 
163  int,
164  << "An error with error number " << arg1
165  << " occurred while calling a Trilinos function");
166 
167  private:
171  MPI::Vector &vector;
172 
176  const size_type index;
177 
182  friend class ::TrilinosWrappers::MPI::Vector;
183  };
184  }
189  namespace
190  {
191 #ifndef DEAL_II_WITH_64BIT_INDICES
192  // define a helper function that queries the global ID of local ID of
193  // an Epetra_BlockMap object by calling either the 32- or 64-bit
194  // function necessary.
195  inline
196  int gid(const Epetra_BlockMap &map, int i)
197  {
198  return map.GID(i);
199  }
200 #else
201  // define a helper function that queries the global ID of local ID of
202  // an Epetra_BlockMap object by calling either the 32- or 64-bit
203  // function necessary.
204  inline
205  long long int gid(const Epetra_BlockMap &map, int i)
206  {
207  return map.GID64(i);
208  }
209 #endif
210  }
211 
218  namespace MPI
219  {
220  class BlockVector;
221 
398  class Vector : public Subscriptor
399  {
400  public:
406  typedef TrilinosScalar value_type;
407  typedef TrilinosScalar real_type;
408  typedef ::types::global_dof_index size_type;
409  typedef value_type *iterator;
410  typedef const value_type *const_iterator;
411  typedef internal::VectorReference reference;
412  typedef const internal::VectorReference const_reference;
413 
423  Vector ();
424 
428  Vector (const Vector &v);
429 
448  explicit Vector (const IndexSet &parallel_partitioning,
449  const MPI_Comm &communicator = MPI_COMM_WORLD);
450 
462  Vector (const IndexSet &local,
463  const IndexSet &ghost,
464  const MPI_Comm &communicator = MPI_COMM_WORLD);
465 
480  Vector (const IndexSet &parallel_partitioning,
481  const Vector &v,
482  const MPI_Comm &communicator = MPI_COMM_WORLD);
483 
496  template <typename Number>
497  Vector (const IndexSet &parallel_partitioning,
498  const ::Vector<Number> &v,
499  const MPI_Comm &communicator = MPI_COMM_WORLD);
500 
505  Vector (Vector &&v);
506 
510  ~Vector () = default;
511 
516  void clear ();
517 
541  void reinit (const Vector &v,
542  const bool omit_zeroing_entries = false,
543  const bool allow_different_maps = false);
544 
567  void reinit (const IndexSet &parallel_partitioning,
568  const MPI_Comm &communicator = MPI_COMM_WORLD,
569  const bool omit_zeroing_entries = false);
570 
596  void reinit (const IndexSet &locally_owned_entries,
597  const IndexSet &ghost_entries,
598  const MPI_Comm &communicator = MPI_COMM_WORLD,
599  const bool vector_writable = false);
600 
604  void reinit (const BlockVector &v,
605  const bool import_data = false);
606 
623  void compress (::VectorOperation::values operation);
624 
637  Vector &operator= (const TrilinosScalar s);
638 
644  Vector &operator= (const Vector &v);
645 
650  Vector &operator= (Vector &&v);
651 
659  template <typename Number>
660  Vector &operator= (const ::Vector<Number> &v);
661 
680  (const ::TrilinosWrappers::SparseMatrix &matrix,
681  const Vector &vector);
682 
688  bool operator== (const Vector &v) const;
689 
695  bool operator!= (const Vector &v) const;
696 
700  size_type size () const;
701 
713  size_type local_size () const;
714 
736  std::pair<size_type, size_type> local_range () const;
737 
745  bool in_local_range (const size_type index) const;
746 
761 
769  bool has_ghost_elements() const;
770 
777  void update_ghost_values () const;
778 
783  TrilinosScalar operator* (const Vector &vec) const;
784 
788  real_type norm_sqr () const;
789 
793  TrilinosScalar mean_value () const;
794 
798  TrilinosScalar min () const;
799 
803  TrilinosScalar max () const;
804 
808  real_type l1_norm () const;
809 
814  real_type l2_norm () const;
815 
820  real_type lp_norm (const TrilinosScalar p) const;
821 
825  real_type linfty_norm () const;
826 
842  TrilinosScalar add_and_dot (const TrilinosScalar a,
843  const Vector &V,
844  const Vector &W);
845 
851  bool all_zero () const;
852 
858  bool is_non_negative () const;
860 
861 
866 
876  reference
877  operator() (const size_type index);
878 
888  TrilinosScalar
889  operator() (const size_type index) const;
890 
896  reference
897  operator[] (const size_type index);
898 
904  TrilinosScalar
905  operator[] (const size_type index) const;
906 
922  void extract_subvector_to (const std::vector<size_type> &indices,
923  std::vector<TrilinosScalar> &values) const;
924 
952  template <typename ForwardIterator, typename OutputIterator>
953  void extract_subvector_to (ForwardIterator indices_begin,
954  const ForwardIterator indices_end,
955  OutputIterator values_begin) const;
956 
967  iterator begin ();
968 
973  const_iterator begin () const;
974 
979  iterator end ();
980 
985  const_iterator end () const;
986 
988 
989 
994 
1001  void set (const std::vector<size_type> &indices,
1002  const std::vector<TrilinosScalar> &values);
1003 
1008  void set (const std::vector<size_type> &indices,
1009  const ::Vector<TrilinosScalar> &values);
1010 
1016  void set (const size_type n_elements,
1017  const size_type *indices,
1018  const TrilinosScalar *values);
1019 
1024  void add (const std::vector<size_type> &indices,
1025  const std::vector<TrilinosScalar> &values);
1026 
1031  void add (const std::vector<size_type> &indices,
1032  const ::Vector<TrilinosScalar> &values);
1033 
1039  void add (const size_type n_elements,
1040  const size_type *indices,
1041  const TrilinosScalar *values);
1042 
1046  Vector &operator*= (const TrilinosScalar factor);
1047 
1051  Vector &operator/= (const TrilinosScalar factor);
1052 
1056  Vector &operator+= (const Vector &V);
1057 
1061  Vector &operator-= (const Vector &V);
1062 
1067  void add (const TrilinosScalar s);
1068 
1081  void add (const Vector &V,
1082  const bool allow_different_maps = false);
1083 
1087  void add (const TrilinosScalar a,
1088  const Vector &V);
1089 
1093  void add (const TrilinosScalar a,
1094  const Vector &V,
1095  const TrilinosScalar b,
1096  const Vector &W);
1097 
1102  void sadd (const TrilinosScalar s,
1103  const Vector &V);
1104 
1108  void sadd (const TrilinosScalar s,
1109  const TrilinosScalar a,
1110  const Vector &V);
1111 
1117  void scale (const Vector &scaling_factors);
1118 
1122  void equ (const TrilinosScalar a,
1123  const Vector &V);
1125 
1130 
1135  const Epetra_MultiVector &trilinos_vector () const;
1136 
1141  Epetra_FEVector &trilinos_vector ();
1142 
1147  const Epetra_Map &vector_partitioner () const;
1148 
1156  void print (std::ostream &out,
1157  const unsigned int precision = 3,
1158  const bool scientific = true,
1159  const bool across = true) const;
1160 
1174  void swap (Vector &v);
1175 
1179  std::size_t memory_consumption () const;
1180 
1185  const MPI_Comm &get_mpi_communicator () const;
1187 
1192 
1197  int,
1198  << "An error with error number " << arg1
1199  << " occurred while calling a Trilinos function");
1200 
1205  size_type, size_type, size_type, size_type,
1206  << "You tried to access element " << arg1
1207  << " of a distributed vector, but this element is not stored "
1208  << "on the current processor. Note: There are "
1209  << arg2 << " elements stored "
1210  << "on the current processor from within the range "
1211  << arg3 << " through " << arg4
1212  << " but Trilinos vectors need not store contiguous "
1213  << "ranges on each processor, and not every element in "
1214  << "this range may in fact be stored locally.");
1215 
1216  private:
1228  Epetra_CombineMode last_action;
1229 
1235 
1241 
1247  std::shared_ptr<Epetra_FEVector> vector;
1248 
1254  std::shared_ptr<Epetra_MultiVector> nonlocal_vector;
1255 
1260 
1265  };
1266 
1267 
1268 
1269 
1270 // ------------------- inline and template functions --------------
1271 
1272 
1281  inline
1282  void swap (Vector &u, Vector &v)
1283  {
1284  u.swap (v);
1285  }
1286  }
1287 
1288 #ifndef DOXYGEN
1289 
1290  namespace internal
1291  {
1292  inline
1293  VectorReference::VectorReference (MPI::Vector &vector,
1294  const size_type index)
1295  :
1296  vector (vector),
1297  index (index)
1298  {}
1299 
1300 
1301  inline
1302  const VectorReference &
1303  VectorReference::operator= (const VectorReference &r) const
1304  {
1305  // as explained in the class
1306  // documentation, this is not the copy
1307  // operator. so simply pass on to the
1308  // "correct" assignment operator
1309  *this = static_cast<TrilinosScalar> (r);
1310 
1311  return *this;
1312  }
1313 
1314 
1315 
1316  inline
1317  VectorReference &
1318  VectorReference::operator= (const VectorReference &r)
1319  {
1320  // as above
1321  *this = static_cast<TrilinosScalar> (r);
1322 
1323  return *this;
1324  }
1325 
1326 
1327  inline
1328  const VectorReference &
1329  VectorReference::operator= (const TrilinosScalar &value) const
1330  {
1331  vector.set (1, &index, &value);
1332  return *this;
1333  }
1334 
1335 
1336 
1337  inline
1338  const VectorReference &
1339  VectorReference::operator+= (const TrilinosScalar &value) const
1340  {
1341  vector.add (1, &index, &value);
1342  return *this;
1343  }
1344 
1345 
1346 
1347  inline
1348  const VectorReference &
1349  VectorReference::operator-= (const TrilinosScalar &value) const
1350  {
1351  TrilinosScalar new_value = -value;
1352  vector.add (1, &index, &new_value);
1353  return *this;
1354  }
1355 
1356 
1357 
1358  inline
1359  const VectorReference &
1360  VectorReference::operator*= (const TrilinosScalar &value) const
1361  {
1362  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1363  vector.set (1, &index, &new_value);
1364  return *this;
1365  }
1366 
1367 
1368 
1369  inline
1370  const VectorReference &
1371  VectorReference::operator/= (const TrilinosScalar &value) const
1372  {
1373  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1374  vector.set (1, &index, &new_value);
1375  return *this;
1376  }
1377  }
1378 
1379  namespace MPI
1380  {
1381  inline
1382  bool
1383  Vector::in_local_range (const size_type index) const
1384  {
1385  std::pair<size_type, size_type> range = local_range();
1386 
1387  return ((index >= range.first) && (index < range.second));
1388  }
1389 
1390 
1391 
1392  inline
1393  IndexSet
1395  {
1397  ExcMessage("The locally owned elements have not been properly initialized!"
1398  " This happens for example if this object has been initialized"
1399  " with exactly one overlapping IndexSet."));
1400  return owned_elements;
1401  }
1402 
1403 
1404 
1405  inline
1406  bool
1408  {
1409  return has_ghosts;
1410  }
1411 
1412 
1413 
1414  inline
1415  void
1417  {}
1418 
1419 
1420 
1421  inline
1422  internal::VectorReference
1423  Vector::operator() (const size_type index)
1424  {
1425  return internal::VectorReference (*this, index);
1426  }
1427 
1428 
1429 
1430  inline
1431  internal::VectorReference
1432  Vector::operator[] (const size_type index)
1433  {
1434  return operator() (index);
1435  }
1436 
1437 
1438 
1439  inline
1440  TrilinosScalar
1441  Vector::operator[] (const size_type index) const
1442  {
1443  return operator() (index);
1444  }
1445 
1446 
1447 
1448  inline
1449  void Vector::extract_subvector_to (const std::vector<size_type> &indices,
1450  std::vector<TrilinosScalar> &values) const
1451  {
1452  for (size_type i = 0; i < indices.size(); ++i)
1453  values[i] = operator()(indices[i]);
1454  }
1455 
1456 
1457 
1458  template <typename ForwardIterator, typename OutputIterator>
1459  inline
1460  void Vector::extract_subvector_to (ForwardIterator indices_begin,
1461  const ForwardIterator indices_end,
1462  OutputIterator values_begin) const
1463  {
1464  while (indices_begin != indices_end)
1465  {
1466  *values_begin = operator()(*indices_begin);
1467  indices_begin++;
1468  values_begin++;
1469  }
1470  }
1471 
1472 
1473 
1474  inline
1475  Vector::iterator
1476  Vector::begin()
1477  {
1478  return (*vector)[0];
1479  }
1480 
1481 
1482 
1483  inline
1484  Vector::iterator
1485  Vector::end()
1486  {
1487  return (*vector)[0]+local_size();
1488  }
1489 
1490 
1491 
1492  inline
1493  Vector::const_iterator
1494  Vector::begin() const
1495  {
1496  return (*vector)[0];
1497  }
1498 
1499 
1500 
1501  inline
1502  Vector::const_iterator
1503  Vector::end() const
1504  {
1505  return (*vector)[0]+local_size();
1506  }
1507 
1508 
1509 
1510  inline
1511  void
1512  Vector::set (const std::vector<size_type> &indices,
1513  const std::vector<TrilinosScalar> &values)
1514  {
1515  // if we have ghost values, do not allow
1516  // writing to this vector at all.
1518 
1519  Assert (indices.size() == values.size(),
1520  ExcDimensionMismatch(indices.size(),values.size()));
1521 
1522  set (indices.size(), indices.data(), values.data());
1523  }
1524 
1525 
1526 
1527  inline
1528  void
1529  Vector::set (const std::vector<size_type> &indices,
1530  const ::Vector<TrilinosScalar> &values)
1531  {
1532  // if we have ghost values, do not allow
1533  // writing to this vector at all.
1535 
1536  Assert (indices.size() == values.size(),
1537  ExcDimensionMismatch(indices.size(),values.size()));
1538 
1539  set (indices.size(), indices.data(), values.begin());
1540  }
1541 
1542 
1543 
1544  inline
1545  void
1546  Vector::set (const size_type n_elements,
1547  const size_type *indices,
1548  const TrilinosScalar *values)
1549  {
1550  // if we have ghost values, do not allow
1551  // writing to this vector at all.
1553 
1554  if (last_action == Add)
1555  {
1556  const int ierr = vector->GlobalAssemble(Add);
1557  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1558  }
1559 
1560  if (last_action != Insert)
1561  last_action = Insert;
1562 
1563  for (size_type i=0; i<n_elements; ++i)
1564  {
1565  const size_type row = indices[i];
1566  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1567  if (local_row != -1)
1568  (*vector)[0][local_row] = values[i];
1569  else
1570  {
1571  const int ierr = vector->ReplaceGlobalValues (1,
1572  (const TrilinosWrappers::types::int_type *)(&row),
1573  &values[i]);
1574  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1575  compressed = false;
1576  }
1577  // in set operation, do not use the pre-allocated vector for nonlocal
1578  // entries even if it exists. This is to ensure that we really only
1579  // set the elements touched by the set() method and not all contained
1580  // in the nonlocal entries vector (there is no way to distinguish them
1581  // on the receiving processor)
1582  }
1583  }
1584 
1585 
1586 
1587  inline
1588  void
1589  Vector::add (const std::vector<size_type> &indices,
1590  const std::vector<TrilinosScalar> &values)
1591  {
1592  // if we have ghost values, do not allow
1593  // writing to this vector at all.
1595  Assert (indices.size() == values.size(),
1596  ExcDimensionMismatch(indices.size(),values.size()));
1597 
1598  add (indices.size(), indices.data(), values.data());
1599  }
1600 
1601 
1602 
1603  inline
1604  void
1605  Vector::add (const std::vector<size_type> &indices,
1606  const ::Vector<TrilinosScalar> &values)
1607  {
1608  // if we have ghost values, do not allow
1609  // writing to this vector at all.
1611  Assert (indices.size() == values.size(),
1612  ExcDimensionMismatch(indices.size(),values.size()));
1613 
1614  add (indices.size(), indices.data(), values.begin());
1615  }
1616 
1617 
1618 
1619  inline
1620  void
1621  Vector::add (const size_type n_elements,
1622  const size_type *indices,
1623  const TrilinosScalar *values)
1624  {
1625  // if we have ghost values, do not allow
1626  // writing to this vector at all.
1628 
1629  if (last_action != Add)
1630  {
1631  if (last_action == Insert)
1632  {
1633  const int ierr = vector->GlobalAssemble(Insert);
1634  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1635  }
1636  last_action = Add;
1637  }
1638 
1639  for (size_type i=0; i<n_elements; ++i)
1640  {
1641  const size_type row = indices[i];
1642  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1643  if (local_row != -1)
1644  (*vector)[0][local_row] += values[i];
1645  else if (nonlocal_vector.get() == nullptr)
1646  {
1647  const int ierr = vector->SumIntoGlobalValues (1,
1648  (const TrilinosWrappers::types::int_type *)(&row),
1649  &values[i]);
1650  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1651  compressed = false;
1652  }
1653  else
1654  {
1655  // use pre-allocated vector for non-local entries if it exists for
1656  // addition operation
1657  const TrilinosWrappers::types::int_type my_row = nonlocal_vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1658  Assert(my_row != -1,
1659  ExcMessage("Attempted to write into off-processor vector entry "
1660  "that has not be specified as being writable upon "
1661  "initialization"));
1662  (*nonlocal_vector)[0][my_row] += values[i];
1663  compressed = false;
1664  }
1665  }
1666  }
1667 
1668 
1669 
1670  inline
1671  Vector::size_type
1672  Vector::size () const
1673  {
1674 #ifndef DEAL_II_WITH_64BIT_INDICES
1675  return (size_type) (vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID());
1676 #else
1677  return (size_type) (vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64());
1678 #endif
1679  }
1680 
1681 
1682 
1683  inline
1684  Vector::size_type
1685  Vector::local_size () const
1686  {
1687  return (size_type) vector->Map().NumMyElements();
1688  }
1689 
1690 
1691 
1692  inline
1693  std::pair<Vector::size_type, Vector::size_type>
1694  Vector::local_range () const
1695  {
1696 #ifndef DEAL_II_WITH_64BIT_INDICES
1697  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1698  const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID()+1;
1699 #else
1700  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID64();
1701  const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID64()+1;
1702 #endif
1703 
1704  Assert (end-begin == vector->Map().NumMyElements(),
1705  ExcMessage ("This function only makes sense if the elements that this "
1706  "vector stores on the current processor form a contiguous range. "
1707  "This does not appear to be the case for the current vector."));
1708 
1709  return std::make_pair (begin, end);
1710  }
1711 
1712 
1713 
1714  inline
1715  TrilinosScalar
1716  Vector::operator* (const Vector &vec) const
1717  {
1718  Assert (vector->Map().SameAs(vec.vector->Map()),
1721 
1722  TrilinosScalar result;
1723 
1724  const int ierr = vector->Dot(*(vec.vector), &result);
1725  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1726 
1727  return result;
1728  }
1729 
1730 
1731 
1732  inline
1734  Vector::norm_sqr () const
1735  {
1736  const TrilinosScalar d = l2_norm();
1737  return d*d;
1738  }
1739 
1740 
1741 
1742  inline
1743  TrilinosScalar
1744  Vector::mean_value () const
1745  {
1747 
1748  TrilinosScalar mean;
1749  const int ierr = vector->MeanValue (&mean);
1750  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1751 
1752  return mean;
1753  }
1754 
1755 
1756 
1757  inline
1758  TrilinosScalar
1759  Vector::min () const
1760  {
1761  TrilinosScalar min_value;
1762  const int ierr = vector->MinValue (&min_value);
1763  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1764 
1765  return min_value;
1766  }
1767 
1768 
1769 
1770  inline
1771  TrilinosScalar
1772  Vector::max () const
1773  {
1774  TrilinosScalar max_value;
1775  const int ierr = vector->MaxValue (&max_value);
1776  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1777 
1778  return max_value;
1779  }
1780 
1781 
1782 
1783  inline
1785  Vector::l1_norm () const
1786  {
1788 
1789  TrilinosScalar d;
1790  const int ierr = vector->Norm1 (&d);
1791  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1792 
1793  return d;
1794  }
1795 
1796 
1797 
1798  inline
1800  Vector::l2_norm () const
1801  {
1803 
1804  TrilinosScalar d;
1805  const int ierr = vector->Norm2 (&d);
1806  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1807 
1808  return d;
1809  }
1810 
1811 
1812 
1813  inline
1815  Vector::lp_norm (const TrilinosScalar p) const
1816  {
1818 
1819  TrilinosScalar norm = 0;
1820  TrilinosScalar sum=0;
1821  const size_type n_local = local_size();
1822 
1823  // loop over all the elements because
1824  // Trilinos does not support lp norms
1825  for (size_type i=0; i<n_local; i++)
1826  sum += std::pow(std::fabs((*vector)[0][i]), p);
1827 
1828  norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
1829 
1830  return norm;
1831  }
1832 
1833 
1834 
1835  inline
1837  Vector::linfty_norm () const
1838  {
1839  // while we disallow the other
1840  // norm operations on ghosted
1841  // vectors, this particular norm
1842  // is safe to run even in the
1843  // presence of ghost elements
1844  TrilinosScalar d;
1845  const int ierr = vector->NormInf (&d);
1846  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1847 
1848  return d;
1849  }
1850 
1851 
1852 
1853  inline
1854  TrilinosScalar
1855  Vector::add_and_dot (const TrilinosScalar a,
1856  const Vector &V,
1857  const Vector &W)
1858  {
1859  this->add(a, V);
1860  return *this * W;
1861  }
1862 
1863 
1864 
1865  // inline also scalar products, vector
1866  // additions etc. since they are all
1867  // representable by a single Trilinos
1868  // call. This reduces the overhead of the
1869  // wrapper class.
1870  inline
1871  Vector &
1872  Vector::operator*= (const TrilinosScalar a)
1873  {
1874  AssertIsFinite(a);
1875 
1876  const int ierr = vector->Scale(a);
1877  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1878 
1879  return *this;
1880  }
1881 
1882 
1883 
1884  inline
1885  Vector &
1886  Vector::operator/= (const TrilinosScalar a)
1887  {
1888  AssertIsFinite(a);
1889 
1890  const TrilinosScalar factor = 1./a;
1891 
1892  AssertIsFinite(factor);
1893 
1894  const int ierr = vector->Scale(factor);
1895  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1896 
1897  return *this;
1898  }
1899 
1900 
1901 
1902  inline
1903  Vector &
1904  Vector::operator+= (const Vector &v)
1905  {
1906  Assert (size() == v.size(),
1907  ExcDimensionMismatch(size(), v.size()));
1908  Assert (vector->Map().SameAs(v.vector->Map()),
1910 
1911  const int ierr = vector->Update (1.0, *(v.vector), 1.0);
1912  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1913 
1914  return *this;
1915  }
1916 
1917 
1918 
1919  inline
1920  Vector &
1921  Vector::operator-= (const Vector &v)
1922  {
1923  Assert (size() == v.size(),
1924  ExcDimensionMismatch(size(), v.size()));
1925  Assert (vector->Map().SameAs(v.vector->Map()),
1927 
1928  const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
1929  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1930 
1931  return *this;
1932  }
1933 
1934 
1935 
1936  inline
1937  void
1938  Vector::add (const TrilinosScalar s)
1939  {
1940  // if we have ghost values, do not allow
1941  // writing to this vector at all.
1943  AssertIsFinite(s);
1944 
1945  size_type n_local = local_size();
1946  for (size_type i=0; i<n_local; i++)
1947  (*vector)[0][i] += s;
1948  }
1949 
1950 
1951 
1952  inline
1953  void
1954  Vector::add (const TrilinosScalar a,
1955  const Vector &v)
1956  {
1957  // if we have ghost values, do not allow
1958  // writing to this vector at all.
1960  Assert (local_size() == v.local_size(),
1961  ExcDimensionMismatch(local_size(), v.local_size()));
1962 
1963  AssertIsFinite(a);
1964 
1965  const int ierr = vector->Update(a, *(v.vector), 1.);
1966  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1967  }
1968 
1969 
1970 
1971  inline
1972  void
1973  Vector::add (const TrilinosScalar a,
1974  const Vector &v,
1975  const TrilinosScalar b,
1976  const Vector &w)
1977  {
1978  // if we have ghost values, do not allow
1979  // writing to this vector at all.
1981  Assert (local_size() == v.local_size(),
1982  ExcDimensionMismatch(local_size(), v.local_size()));
1983  Assert (local_size() == w.local_size(),
1984  ExcDimensionMismatch(local_size(), w.local_size()));
1985 
1986  AssertIsFinite(a);
1987  AssertIsFinite(b);
1988 
1989  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
1990 
1991  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1992  }
1993 
1994 
1995 
1996  inline
1997  void
1998  Vector::sadd (const TrilinosScalar s,
1999  const Vector &v)
2000  {
2001  // if we have ghost values, do not allow
2002  // writing to this vector at all.
2004  Assert (size() == v.size(),
2005  ExcDimensionMismatch (size(), v.size()));
2006 
2007  AssertIsFinite(s);
2008 
2009  // We assume that the vectors have the same Map
2010  // if the local size is the same and if the vectors are not ghosted
2011  if (local_size() == v.local_size() && !v.has_ghost_elements())
2012  {
2013  Assert (this->vector->Map().SameAs(v.vector->Map())==true,
2015  const int ierr = vector->Update(1., *(v.vector), s);
2016  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2017  }
2018  else
2019  {
2020  (*this) *= s;
2021  this->add(v, true);
2022  }
2023  }
2024 
2025 
2026 
2027  inline
2028  void
2029  Vector::sadd (const TrilinosScalar s,
2030  const TrilinosScalar a,
2031  const Vector &v)
2032  {
2033  // if we have ghost values, do not allow
2034  // writing to this vector at all.
2036  Assert (size() == v.size(),
2037  ExcDimensionMismatch (size(), v.size()));
2038  AssertIsFinite(s);
2039  AssertIsFinite(a);
2040 
2041  // We assume that the vectors have the same Map
2042  // if the local size is the same and if the vectors are not ghosted
2043  if (local_size() == v.local_size() && !v.has_ghost_elements())
2044  {
2045  Assert (this->vector->Map().SameAs(v.vector->Map())==true,
2047  const int ierr = vector->Update(a, *(v.vector), s);
2048  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2049  }
2050  else
2051  {
2052  (*this) *= s;
2053  Vector tmp = v;
2054  tmp *= a;
2055  this->add(tmp, true);
2056  }
2057  }
2058 
2059 
2060 
2061  inline
2062  void
2063  Vector::scale (const Vector &factors)
2064  {
2065  // if we have ghost values, do not allow
2066  // writing to this vector at all.
2068  Assert (local_size() == factors.local_size(),
2069  ExcDimensionMismatch(local_size(), factors.local_size()));
2070 
2071  const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
2072  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2073  }
2074 
2075 
2076 
2077  inline
2078  void
2079  Vector::equ (const TrilinosScalar a,
2080  const Vector &v)
2081  {
2082  // if we have ghost values, do not allow
2083  // writing to this vector at all.
2085  AssertIsFinite(a);
2086 
2087  // If we don't have the same map, copy.
2088  if (vector->Map().SameAs(v.vector->Map())==false)
2089  {
2090  this->sadd(0., a, v);
2091  }
2092  else
2093  {
2094  // Otherwise, just update
2095  int ierr = vector->Update(a, *v.vector, 0.0);
2096  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2097 
2098  last_action = Zero;
2099  }
2100 
2101  }
2102 
2103 
2104 
2105  inline
2106  const Epetra_MultiVector &
2107  Vector::trilinos_vector () const
2108  {
2109  return static_cast<const Epetra_MultiVector &>(*vector);
2110  }
2111 
2112 
2113 
2114  inline
2115  Epetra_FEVector &
2117  {
2118  return *vector;
2119  }
2120 
2121 
2122 
2123  inline
2124  const Epetra_Map &
2126  {
2127  return static_cast<const Epetra_Map &>(vector->Map());
2128  }
2129 
2130 
2131 
2132  inline
2133  const MPI_Comm &
2135  {
2136  static MPI_Comm comm;
2137 
2138 #ifdef DEAL_II_WITH_MPI
2139 
2140  const Epetra_MpiComm *mpi_comm
2141  = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2142  comm = mpi_comm->Comm();
2143 
2144 #else
2145 
2146  comm = MPI_COMM_SELF;
2147 
2148 #endif
2149 
2150  return comm;
2151  }
2152 
2153  template <typename number>
2154  Vector::Vector (const IndexSet &parallel_partitioner,
2155  const ::Vector<number> &v,
2156  const MPI_Comm &communicator)
2157  {
2158  *this = Vector(parallel_partitioner.make_trilinos_map (communicator, true),
2159  v);
2160  owned_elements = parallel_partitioner;
2161  }
2162 
2163 
2164 
2165  inline
2166  Vector &
2167  Vector::operator= (const TrilinosScalar s)
2168  {
2169  AssertIsFinite(s);
2170 
2171  int ierr = vector->PutScalar(s);
2172  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2173 
2174  if (nonlocal_vector.get() != nullptr)
2175  {
2176  ierr = nonlocal_vector->PutScalar(0.);
2177  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2178  }
2179 
2180  return *this;
2181  }
2182  } /* end of namespace MPI */
2183 
2184 #endif /* DOXYGEN */
2185 
2186 } /* end of namespace TrilinosWrappers */
2187 
2191 namespace internal
2192 {
2193  namespace LinearOperator
2194  {
2195  template <typename> class ReinitHelper;
2196 
2201  template <>
2202  class ReinitHelper<TrilinosWrappers::MPI::Vector>
2203  {
2204  public:
2205  template <typename Matrix>
2206  static
2207  void reinit_range_vector (const Matrix &matrix,
2209  bool omit_zeroing_entries)
2210  {
2211  v.reinit(matrix.locally_owned_range_indices(), matrix.get_mpi_communicator(), omit_zeroing_entries);
2212  }
2213 
2214  template <typename Matrix>
2215  static
2216  void reinit_domain_vector(const Matrix &matrix,
2218  bool omit_zeroing_entries)
2219  {
2220  v.reinit(matrix.locally_owned_domain_indices(), matrix.get_mpi_communicator(), omit_zeroing_entries);
2221  }
2222  };
2223 
2224  } /* namespace LinearOperator */
2225 } /* namespace internal */
2226 
2227 
2228 
2234 template <>
2235 struct is_serial_vector< TrilinosWrappers::MPI::Vector > : std::false_type
2236 {
2237 };
2238 
2239 
2240 DEAL_II_NAMESPACE_CLOSE
2241 
2242 #endif // DEAL_II_WITH_TRILINOS
2243 
2244 /*---------------------------- trilinos_vector.h ---------------------------*/
2245 
2246 #endif // dealii_trilinos_vector_h
Vector & operator=(const TrilinosScalar s)
static::ExceptionBase & ExcTrilinosError(int arg1)
real_type linfty_norm() const
void sadd(const TrilinosScalar s, const Vector &V)
TrilinosScalar max() const
reference operator[](const size_type index)
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)
real_type l1_norm() const
void equ(const TrilinosScalar a, const Vector &V)
std::pair< size_type, size_type > local_range() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:398
TrilinosScalar min() const
size_type size() const
Definition: index_set.h:1553
static void reinit_domain_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
std::shared_ptr< Epetra_MultiVector > nonlocal_vector
std::shared_ptr< Epetra_FEVector > vector
static::ExceptionBase & ExcMessage(std::string arg1)
numbers::NumberTraits< Number >::real_type real_type
Definition: vector.h:158
bool operator!=(const Vector &v) const
reference operator()(const size_type index)
void set(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:593
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:337
real_type l2_norm() const
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void compress(::VectorOperation::values operation)
#define DeclException0(Exception0)
Definition: exceptions.h:570
IndexSet locally_owned_elements() const
std::size_t size() const
Vector & operator/=(const TrilinosScalar factor)
real_type norm_sqr() const
const Epetra_Map & vector_partitioner() const
void swap(Vector &u, Vector &v)
const MPI_Comm & get_mpi_communicator() const
friend class internal::VectorReference
void import_nonlocal_data_for_fe(const ::TrilinosWrappers::SparseMatrix &matrix, const Vector &vector)
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 void reinit_range_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
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:532
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< TrilinosScalar > &values) const
real_type lp_norm(const TrilinosScalar p) const
bool operator==(const Vector &v) const
TrilinosScalar operator*(const Vector &vec) const
bool in_local_range(const size_type index) const
void scale(const Vector &scaling_factors)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:629
Vector & operator-=(const Vector &V)
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)
Vector & operator+=(const Vector &V)
static::ExceptionBase & ExcDifferentParallelPartitioning()
Vector & operator*=(const TrilinosScalar factor)
static::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertIsFinite(number)
Definition: exceptions.h:1219
const Epetra_MultiVector & trilinos_vector() const