Reference documentation for deal.II version Git ce1d1e1 2017-06-28 06:17:25 -0500
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 ();
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 
913  void extract_subvector_to (const std::vector<size_type> &indices,
914  std::vector<TrilinosScalar> &values) const;
915 
920  template <typename ForwardIterator, typename OutputIterator>
921  void extract_subvector_to (ForwardIterator indices_begin,
922  const ForwardIterator indices_end,
923  OutputIterator values_begin) const;
924 
935  iterator begin ();
936 
941  const_iterator begin () const;
942 
947  iterator end ();
948 
953  const_iterator end () const;
954 
956 
957 
962 
969  void set (const std::vector<size_type> &indices,
970  const std::vector<TrilinosScalar> &values);
971 
976  void set (const std::vector<size_type> &indices,
977  const ::Vector<TrilinosScalar> &values);
978 
984  void set (const size_type n_elements,
985  const size_type *indices,
986  const TrilinosScalar *values);
987 
992  void add (const std::vector<size_type> &indices,
993  const std::vector<TrilinosScalar> &values);
994 
999  void add (const std::vector<size_type> &indices,
1000  const ::Vector<TrilinosScalar> &values);
1001 
1007  void add (const size_type n_elements,
1008  const size_type *indices,
1009  const TrilinosScalar *values);
1010 
1014  Vector &operator*= (const TrilinosScalar factor);
1015 
1019  Vector &operator/= (const TrilinosScalar factor);
1020 
1024  Vector &operator+= (const Vector &V);
1025 
1029  Vector &operator-= (const Vector &V);
1030 
1035  void add (const TrilinosScalar s);
1036 
1049  void add (const Vector &V,
1050  const bool allow_different_maps = false);
1051 
1055  void add (const TrilinosScalar a,
1056  const Vector &V);
1057 
1061  void add (const TrilinosScalar a,
1062  const Vector &V,
1063  const TrilinosScalar b,
1064  const Vector &W);
1065 
1070  void sadd (const TrilinosScalar s,
1071  const Vector &V);
1072 
1076  void sadd (const TrilinosScalar s,
1077  const TrilinosScalar a,
1078  const Vector &V);
1079 
1085  void scale (const Vector &scaling_factors);
1086 
1090  void equ (const TrilinosScalar a,
1091  const Vector &V);
1093 
1098 
1103  const Epetra_MultiVector &trilinos_vector () const;
1104 
1109  Epetra_FEVector &trilinos_vector ();
1110 
1115  const Epetra_Map &vector_partitioner () const;
1116 
1124  void print (std::ostream &out,
1125  const unsigned int precision = 3,
1126  const bool scientific = true,
1127  const bool across = true) const;
1128 
1142  void swap (Vector &v);
1143 
1147  std::size_t memory_consumption () const;
1148 
1153  const MPI_Comm &get_mpi_communicator () const;
1155 
1160 
1165  int,
1166  << "An error with error number " << arg1
1167  << " occurred while calling a Trilinos function");
1168 
1173  size_type, size_type, size_type, size_type,
1174  << "You tried to access element " << arg1
1175  << " of a distributed vector, but this element is not stored "
1176  << "on the current processor. Note: There are "
1177  << arg2 << " elements stored "
1178  << "on the current processor from within the range "
1179  << arg3 << " through " << arg4
1180  << " but Trilinos vectors need not store contiguous "
1181  << "ranges on each processor, and not every element in "
1182  << "this range may in fact be stored locally.");
1183 
1184  private:
1196  Epetra_CombineMode last_action;
1197 
1203 
1209 
1215  std::shared_ptr<Epetra_FEVector> vector;
1216 
1222  std::shared_ptr<Epetra_MultiVector> nonlocal_vector;
1223 
1228 
1233  };
1234 
1235 
1236 
1237 
1238 // ------------------- inline and template functions --------------
1239 
1240 
1249  inline
1250  void swap (Vector &u, Vector &v)
1251  {
1252  u.swap (v);
1253  }
1254  }
1255 
1256 #ifndef DOXYGEN
1257 
1258  namespace internal
1259  {
1260  inline
1261  VectorReference::VectorReference (MPI::Vector &vector,
1262  const size_type index)
1263  :
1264  vector (vector),
1265  index (index)
1266  {}
1267 
1268 
1269  inline
1270  const VectorReference &
1271  VectorReference::operator= (const VectorReference &r) const
1272  {
1273  // as explained in the class
1274  // documentation, this is not the copy
1275  // operator. so simply pass on to the
1276  // "correct" assignment operator
1277  *this = static_cast<TrilinosScalar> (r);
1278 
1279  return *this;
1280  }
1281 
1282 
1283 
1284  inline
1285  VectorReference &
1286  VectorReference::operator= (const VectorReference &r)
1287  {
1288  // as above
1289  *this = static_cast<TrilinosScalar> (r);
1290 
1291  return *this;
1292  }
1293 
1294 
1295  inline
1296  const VectorReference &
1297  VectorReference::operator= (const TrilinosScalar &value) const
1298  {
1299  vector.set (1, &index, &value);
1300  return *this;
1301  }
1302 
1303 
1304 
1305  inline
1306  const VectorReference &
1307  VectorReference::operator+= (const TrilinosScalar &value) const
1308  {
1309  vector.add (1, &index, &value);
1310  return *this;
1311  }
1312 
1313 
1314 
1315  inline
1316  const VectorReference &
1317  VectorReference::operator-= (const TrilinosScalar &value) const
1318  {
1319  TrilinosScalar new_value = -value;
1320  vector.add (1, &index, &new_value);
1321  return *this;
1322  }
1323 
1324 
1325 
1326  inline
1327  const VectorReference &
1328  VectorReference::operator*= (const TrilinosScalar &value) const
1329  {
1330  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1331  vector.set (1, &index, &new_value);
1332  return *this;
1333  }
1334 
1335 
1336 
1337  inline
1338  const VectorReference &
1339  VectorReference::operator/= (const TrilinosScalar &value) const
1340  {
1341  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1342  vector.set (1, &index, &new_value);
1343  return *this;
1344  }
1345  }
1346 
1347  namespace MPI
1348  {
1349  inline
1350  bool
1351  Vector::in_local_range (const size_type index) const
1352  {
1353  std::pair<size_type, size_type> range = local_range();
1354 
1355  return ((index >= range.first) && (index < range.second));
1356  }
1357 
1358 
1359 
1360  inline
1361  IndexSet
1363  {
1365  ExcMessage("The locally owned elements have not been properly initialized!"
1366  " This happens for example if this object has been initialized"
1367  " with exactly one overlapping IndexSet."));
1368  return owned_elements;
1369  }
1370 
1371 
1372 
1373  inline
1374  bool
1376  {
1377  return has_ghosts;
1378  }
1379 
1380 
1381 
1382  inline
1383  void
1385  {}
1386 
1387 
1388 
1389  inline
1390  internal::VectorReference
1391  Vector::operator() (const size_type index)
1392  {
1393  return internal::VectorReference (*this, index);
1394  }
1395 
1396 
1397 
1398  inline
1399  internal::VectorReference
1400  Vector::operator[] (const size_type index)
1401  {
1402  return operator() (index);
1403  }
1404 
1405 
1406 
1407  inline
1408  TrilinosScalar
1409  Vector::operator[] (const size_type index) const
1410  {
1411  return operator() (index);
1412  }
1413 
1414 
1415 
1416  inline
1417  void Vector::extract_subvector_to (const std::vector<size_type> &indices,
1418  std::vector<TrilinosScalar> &values) const
1419  {
1420  for (size_type i = 0; i < indices.size(); ++i)
1421  values[i] = operator()(indices[i]);
1422  }
1423 
1424 
1425 
1426  template <typename ForwardIterator, typename OutputIterator>
1427  inline
1428  void Vector::extract_subvector_to (ForwardIterator indices_begin,
1429  const ForwardIterator indices_end,
1430  OutputIterator values_begin) const
1431  {
1432  while (indices_begin != indices_end)
1433  {
1434  *values_begin = operator()(*indices_begin);
1435  indices_begin++;
1436  values_begin++;
1437  }
1438  }
1439 
1440 
1441 
1442  inline
1443  Vector::iterator
1444  Vector::begin()
1445  {
1446  return (*vector)[0];
1447  }
1448 
1449 
1450 
1451  inline
1452  Vector::iterator
1453  Vector::end()
1454  {
1455  return (*vector)[0]+local_size();
1456  }
1457 
1458 
1459 
1460  inline
1461  Vector::const_iterator
1462  Vector::begin() const
1463  {
1464  return (*vector)[0];
1465  }
1466 
1467 
1468 
1469  inline
1470  Vector::const_iterator
1471  Vector::end() const
1472  {
1473  return (*vector)[0]+local_size();
1474  }
1475 
1476 
1477 
1478  inline
1479  void
1480  Vector::set (const std::vector<size_type> &indices,
1481  const std::vector<TrilinosScalar> &values)
1482  {
1483  // if we have ghost values, do not allow
1484  // writing to this vector at all.
1486 
1487  Assert (indices.size() == values.size(),
1488  ExcDimensionMismatch(indices.size(),values.size()));
1489 
1490  set (indices.size(), &indices[0], &values[0]);
1491  }
1492 
1493 
1494 
1495  inline
1496  void
1497  Vector::set (const std::vector<size_type> &indices,
1498  const ::Vector<TrilinosScalar> &values)
1499  {
1500  // if we have ghost values, do not allow
1501  // writing to this vector at all.
1503 
1504  Assert (indices.size() == values.size(),
1505  ExcDimensionMismatch(indices.size(),values.size()));
1506 
1507  set (indices.size(), &indices[0], values.begin());
1508  }
1509 
1510 
1511 
1512  inline
1513  void
1514  Vector::set (const size_type n_elements,
1515  const size_type *indices,
1516  const TrilinosScalar *values)
1517  {
1518  // if we have ghost values, do not allow
1519  // writing to this vector at all.
1521 
1522  if (last_action == Add)
1523  {
1524  const int ierr = vector->GlobalAssemble(Add);
1525  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1526  }
1527 
1528  if (last_action != Insert)
1529  last_action = Insert;
1530 
1531  for (size_type i=0; i<n_elements; ++i)
1532  {
1533  const size_type row = indices[i];
1534  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1535  if (local_row != -1)
1536  (*vector)[0][local_row] = values[i];
1537  else
1538  {
1539  const int ierr = vector->ReplaceGlobalValues (1,
1540  (const TrilinosWrappers::types::int_type *)(&row),
1541  &values[i]);
1542  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1543  compressed = false;
1544  }
1545  // in set operation, do not use the pre-allocated vector for nonlocal
1546  // entries even if it exists. This is to ensure that we really only
1547  // set the elements touched by the set() method and not all contained
1548  // in the nonlocal entries vector (there is no way to distinguish them
1549  // on the receiving processor)
1550  }
1551  }
1552 
1553 
1554 
1555  inline
1556  void
1557  Vector::add (const std::vector<size_type> &indices,
1558  const std::vector<TrilinosScalar> &values)
1559  {
1560  // if we have ghost values, do not allow
1561  // writing to this vector at all.
1563  Assert (indices.size() == values.size(),
1564  ExcDimensionMismatch(indices.size(),values.size()));
1565 
1566  add (indices.size(), &indices[0], &values[0]);
1567  }
1568 
1569 
1570 
1571  inline
1572  void
1573  Vector::add (const std::vector<size_type> &indices,
1574  const ::Vector<TrilinosScalar> &values)
1575  {
1576  // if we have ghost values, do not allow
1577  // writing to this vector at all.
1579  Assert (indices.size() == values.size(),
1580  ExcDimensionMismatch(indices.size(),values.size()));
1581 
1582  add (indices.size(), &indices[0], values.begin());
1583  }
1584 
1585 
1586 
1587  inline
1588  void
1589  Vector::add (const size_type n_elements,
1590  const size_type *indices,
1591  const TrilinosScalar *values)
1592  {
1593  // if we have ghost values, do not allow
1594  // writing to this vector at all.
1596 
1597  if (last_action != Add)
1598  {
1599  if (last_action == Insert)
1600  {
1601  const int ierr = vector->GlobalAssemble(Insert);
1602  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1603  }
1604  last_action = Add;
1605  }
1606 
1607  for (size_type i=0; i<n_elements; ++i)
1608  {
1609  const size_type row = indices[i];
1610  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1611  if (local_row != -1)
1612  (*vector)[0][local_row] += values[i];
1613  else if (nonlocal_vector.get() == nullptr)
1614  {
1615  const int ierr = vector->SumIntoGlobalValues (1,
1616  (const TrilinosWrappers::types::int_type *)(&row),
1617  &values[i]);
1618  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1619  compressed = false;
1620  }
1621  else
1622  {
1623  // use pre-allocated vector for non-local entries if it exists for
1624  // addition operation
1625  const TrilinosWrappers::types::int_type my_row = nonlocal_vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1626  Assert(my_row != -1,
1627  ExcMessage("Attempted to write into off-processor vector entry "
1628  "that has not be specified as being writable upon "
1629  "initialization"));
1630  (*nonlocal_vector)[0][my_row] += values[i];
1631  compressed = false;
1632  }
1633  }
1634  }
1635 
1636 
1637 
1638  inline
1639  Vector::size_type
1640  Vector::size () const
1641  {
1642 #ifndef DEAL_II_WITH_64BIT_INDICES
1643  return (size_type) (vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID());
1644 #else
1645  return (size_type) (vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64());
1646 #endif
1647  }
1648 
1649 
1650 
1651  inline
1652  Vector::size_type
1653  Vector::local_size () const
1654  {
1655  return (size_type) vector->Map().NumMyElements();
1656  }
1657 
1658 
1659 
1660  inline
1661  std::pair<Vector::size_type, Vector::size_type>
1662  Vector::local_range () const
1663  {
1664 #ifndef DEAL_II_WITH_64BIT_INDICES
1665  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1666  const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID()+1;
1667 #else
1668  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID64();
1669  const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID64()+1;
1670 #endif
1671 
1672  Assert (end-begin == vector->Map().NumMyElements(),
1673  ExcMessage ("This function only makes sense if the elements that this "
1674  "vector stores on the current processor form a contiguous range. "
1675  "This does not appear to be the case for the current vector."));
1676 
1677  return std::make_pair (begin, end);
1678  }
1679 
1680 
1681 
1682  inline
1683  TrilinosScalar
1684  Vector::operator* (const Vector &vec) const
1685  {
1686  Assert (vector->Map().SameAs(vec.vector->Map()),
1689 
1690  TrilinosScalar result;
1691 
1692  const int ierr = vector->Dot(*(vec.vector), &result);
1693  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1694 
1695  return result;
1696  }
1697 
1698 
1699 
1700  inline
1702  Vector::norm_sqr () const
1703  {
1704  const TrilinosScalar d = l2_norm();
1705  return d*d;
1706  }
1707 
1708 
1709 
1710  inline
1711  TrilinosScalar
1712  Vector::mean_value () const
1713  {
1715 
1716  TrilinosScalar mean;
1717  const int ierr = vector->MeanValue (&mean);
1718  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1719 
1720  return mean;
1721  }
1722 
1723 
1724 
1725  inline
1726  TrilinosScalar
1727  Vector::min () const
1728  {
1729  TrilinosScalar min_value;
1730  const int ierr = vector->MinValue (&min_value);
1731  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1732 
1733  return min_value;
1734  }
1735 
1736 
1737 
1738  inline
1739  TrilinosScalar
1740  Vector::max () const
1741  {
1742  TrilinosScalar max_value;
1743  const int ierr = vector->MaxValue (&max_value);
1744  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1745 
1746  return max_value;
1747  }
1748 
1749 
1750 
1751  inline
1753  Vector::l1_norm () const
1754  {
1756 
1757  TrilinosScalar d;
1758  const int ierr = vector->Norm1 (&d);
1759  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1760 
1761  return d;
1762  }
1763 
1764 
1765 
1766  inline
1768  Vector::l2_norm () const
1769  {
1771 
1772  TrilinosScalar d;
1773  const int ierr = vector->Norm2 (&d);
1774  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1775 
1776  return d;
1777  }
1778 
1779 
1780 
1781  inline
1783  Vector::lp_norm (const TrilinosScalar p) const
1784  {
1786 
1787  TrilinosScalar norm = 0;
1788  TrilinosScalar sum=0;
1789  const size_type n_local = local_size();
1790 
1791  // loop over all the elements because
1792  // Trilinos does not support lp norms
1793  for (size_type i=0; i<n_local; i++)
1794  sum += std::pow(std::fabs((*vector)[0][i]), p);
1795 
1796  norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
1797 
1798  return norm;
1799  }
1800 
1801 
1802 
1803  inline
1805  Vector::linfty_norm () const
1806  {
1807  // while we disallow the other
1808  // norm operations on ghosted
1809  // vectors, this particular norm
1810  // is safe to run even in the
1811  // presence of ghost elements
1812  TrilinosScalar d;
1813  const int ierr = vector->NormInf (&d);
1814  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1815 
1816  return d;
1817  }
1818 
1819 
1820 
1821  inline
1822  TrilinosScalar
1823  Vector::add_and_dot (const TrilinosScalar a,
1824  const Vector &V,
1825  const Vector &W)
1826  {
1827  this->add(a, V);
1828  return *this * W;
1829  }
1830 
1831 
1832 
1833  // inline also scalar products, vector
1834  // additions etc. since they are all
1835  // representable by a single Trilinos
1836  // call. This reduces the overhead of the
1837  // wrapper class.
1838  inline
1839  Vector &
1840  Vector::operator*= (const TrilinosScalar a)
1841  {
1842  AssertIsFinite(a);
1843 
1844  const int ierr = vector->Scale(a);
1845  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1846 
1847  return *this;
1848  }
1849 
1850 
1851 
1852  inline
1853  Vector &
1854  Vector::operator/= (const TrilinosScalar a)
1855  {
1856  AssertIsFinite(a);
1857 
1858  const TrilinosScalar factor = 1./a;
1859 
1860  AssertIsFinite(factor);
1861 
1862  const int ierr = vector->Scale(factor);
1863  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1864 
1865  return *this;
1866  }
1867 
1868 
1869 
1870  inline
1871  Vector &
1872  Vector::operator+= (const Vector &v)
1873  {
1874  Assert (size() == v.size(),
1875  ExcDimensionMismatch(size(), v.size()));
1876  Assert (vector->Map().SameAs(v.vector->Map()),
1878 
1879  const int ierr = vector->Update (1.0, *(v.vector), 1.0);
1880  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1881 
1882  return *this;
1883  }
1884 
1885 
1886 
1887  inline
1888  Vector &
1889  Vector::operator-= (const Vector &v)
1890  {
1891  Assert (size() == v.size(),
1892  ExcDimensionMismatch(size(), v.size()));
1893  Assert (vector->Map().SameAs(v.vector->Map()),
1895 
1896  const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
1897  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1898 
1899  return *this;
1900  }
1901 
1902 
1903 
1904  inline
1905  void
1906  Vector::add (const TrilinosScalar s)
1907  {
1908  // if we have ghost values, do not allow
1909  // writing to this vector at all.
1911  AssertIsFinite(s);
1912 
1913  size_type n_local = local_size();
1914  for (size_type i=0; i<n_local; i++)
1915  (*vector)[0][i] += s;
1916  }
1917 
1918 
1919 
1920  inline
1921  void
1922  Vector::add (const TrilinosScalar a,
1923  const Vector &v)
1924  {
1925  // if we have ghost values, do not allow
1926  // writing to this vector at all.
1928  Assert (local_size() == v.local_size(),
1929  ExcDimensionMismatch(local_size(), v.local_size()));
1930 
1931  AssertIsFinite(a);
1932 
1933  const int ierr = vector->Update(a, *(v.vector), 1.);
1934  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1935  }
1936 
1937 
1938 
1939  inline
1940  void
1941  Vector::add (const TrilinosScalar a,
1942  const Vector &v,
1943  const TrilinosScalar b,
1944  const Vector &w)
1945  {
1946  // if we have ghost values, do not allow
1947  // writing to this vector at all.
1949  Assert (local_size() == v.local_size(),
1950  ExcDimensionMismatch(local_size(), v.local_size()));
1951  Assert (local_size() == w.local_size(),
1952  ExcDimensionMismatch(local_size(), w.local_size()));
1953 
1954  AssertIsFinite(a);
1955  AssertIsFinite(b);
1956 
1957  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
1958 
1959  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1960  }
1961 
1962 
1963 
1964  inline
1965  void
1966  Vector::sadd (const TrilinosScalar s,
1967  const Vector &v)
1968  {
1969  // if we have ghost values, do not allow
1970  // writing to this vector at all.
1972  Assert (size() == v.size(),
1973  ExcDimensionMismatch (size(), v.size()));
1974 
1975  AssertIsFinite(s);
1976 
1977  // We assume that the vectors have the same Map
1978  // if the local size is the same and if the vectors are not ghosted
1979  if (local_size() == v.local_size() && !v.has_ghost_elements())
1980  {
1981  Assert (this->vector->Map().SameAs(v.vector->Map())==true,
1983  const int ierr = vector->Update(1., *(v.vector), s);
1984  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1985  }
1986  else
1987  {
1988  (*this) *= s;
1989  this->add(v, true);
1990  }
1991  }
1992 
1993 
1994 
1995  inline
1996  void
1997  Vector::sadd (const TrilinosScalar s,
1998  const TrilinosScalar a,
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  AssertIsFinite(s);
2007  AssertIsFinite(a);
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(a, *(v.vector), s);
2016  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2017  }
2018  else
2019  {
2020  (*this) *= s;
2021  Vector tmp = v;
2022  tmp *= a;
2023  this->add(tmp, true);
2024  }
2025  }
2026 
2027 
2028 
2029  inline
2030  void
2031  Vector::scale (const Vector &factors)
2032  {
2033  // if we have ghost values, do not allow
2034  // writing to this vector at all.
2036  Assert (local_size() == factors.local_size(),
2037  ExcDimensionMismatch(local_size(), factors.local_size()));
2038 
2039  const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
2040  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2041  }
2042 
2043 
2044 
2045  inline
2046  void
2047  Vector::equ (const TrilinosScalar a,
2048  const Vector &v)
2049  {
2050  // if we have ghost values, do not allow
2051  // writing to this vector at all.
2053  AssertIsFinite(a);
2054 
2055  // If we don't have the same map, copy.
2056  if (vector->Map().SameAs(v.vector->Map())==false)
2057  {
2058  this->sadd(0., a, v);
2059  }
2060  else
2061  {
2062  // Otherwise, just update
2063  int ierr = vector->Update(a, *v.vector, 0.0);
2064  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2065 
2066  last_action = Zero;
2067  }
2068 
2069  }
2070 
2071 
2072 
2073  inline
2074  const Epetra_MultiVector &
2075  Vector::trilinos_vector () const
2076  {
2077  return static_cast<const Epetra_MultiVector &>(*vector);
2078  }
2079 
2080 
2081 
2082  inline
2083  Epetra_FEVector &
2085  {
2086  return *vector;
2087  }
2088 
2089 
2090 
2091  inline
2092  const Epetra_Map &
2094  {
2095  return static_cast<const Epetra_Map &>(vector->Map());
2096  }
2097 
2098 
2099 
2100  inline
2101  const MPI_Comm &
2103  {
2104  static MPI_Comm comm;
2105 
2106 #ifdef DEAL_II_WITH_MPI
2107 
2108  const Epetra_MpiComm *mpi_comm
2109  = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2110  comm = mpi_comm->Comm();
2111 
2112 #else
2113 
2114  comm = MPI_COMM_SELF;
2115 
2116 #endif
2117 
2118  return comm;
2119  }
2120 
2121  template <typename number>
2122  Vector::Vector (const IndexSet &parallel_partitioner,
2123  const ::Vector<number> &v,
2124  const MPI_Comm &communicator)
2125  {
2126  *this = Vector(parallel_partitioner.make_trilinos_map (communicator, true),
2127  v);
2128  owned_elements = parallel_partitioner;
2129  }
2130 
2131 
2132 
2133  inline
2134  Vector &
2135  Vector::operator= (const TrilinosScalar s)
2136  {
2137  AssertIsFinite(s);
2138 
2139  int ierr = vector->PutScalar(s);
2140  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2141 
2142  if (nonlocal_vector.get() != nullptr)
2143  {
2144  ierr = nonlocal_vector->PutScalar(0.);
2145  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2146  }
2147 
2148  return *this;
2149  }
2150  } /* end of namespace MPI */
2151 
2152 #endif /* DOXYGEN */
2153 
2154 } /* end of namespace TrilinosWrappers */
2155 
2159 namespace internal
2160 {
2161  namespace LinearOperator
2162  {
2163  template <typename> class ReinitHelper;
2164 
2169  template <>
2170  class ReinitHelper<TrilinosWrappers::MPI::Vector>
2171  {
2172  public:
2173  template <typename Matrix>
2174  static
2175  void reinit_range_vector (const Matrix &matrix,
2177  bool omit_zeroing_entries)
2178  {
2179  v.reinit(matrix.locally_owned_range_indices(), matrix.get_mpi_communicator(), omit_zeroing_entries);
2180  }
2181 
2182  template <typename Matrix>
2183  static
2184  void reinit_domain_vector(const Matrix &matrix,
2186  bool omit_zeroing_entries)
2187  {
2188  v.reinit(matrix.locally_owned_domain_indices(), matrix.get_mpi_communicator(), omit_zeroing_entries);
2189  }
2190  };
2191 
2192  } /* namespace LinearOperator */
2193 } /* namespace internal */
2194 
2195 
2196 
2202 template <>
2203 struct is_serial_vector< TrilinosWrappers::MPI::Vector > : std::false_type
2204 {
2205 };
2206 
2207 
2208 DEAL_II_NAMESPACE_CLOSE
2209 
2210 #endif // DEAL_II_WITH_TRILINOS
2211 
2212 /*---------------------------- trilinos_vector.h ---------------------------*/
2213 
2214 #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:388
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:583
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:327
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:560
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:539
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:619
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:1201
const Epetra_MultiVector & trilinos_vector() const