Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
trilinos_vector_base.h
1 // ---------------------------------------------------------------------
2 // @f$Id: trilinos_vector_base.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2008 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__trilinos_vector_base_h
18 #define __deal2__trilinos_vector_base_h
19 
20 
21 #include <deal.II/base/config.h>
22 
23 #ifdef DEAL_II_WITH_TRILINOS
24 
25 #include <deal.II/base/utilities.h>
26 # include <deal.II/base/std_cxx1x/shared_ptr.h>
27 # include <deal.II/base/subscriptor.h>
28 # include <deal.II/lac/exceptions.h>
29 # include <deal.II/lac/vector.h>
30 
31 # include <vector>
32 # include <utility>
33 # include <memory>
34 
35 # define TrilinosScalar double
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 
46 
47 // forward declaration
48 template <typename number> class Vector;
49 
50 
55 namespace TrilinosWrappers
56 {
57  // forward declaration
58  class VectorBase;
59 
70  namespace internal
71  {
75  typedef ::types::global_dof_index size_type;
76 
93  class VectorReference
94  {
95  private:
102  VectorReference (VectorBase &vector,
103  const size_type index);
104 
105  public:
106 
131  const VectorReference &
132  operator = (const VectorReference &r) const;
133 
138  const VectorReference &
139  operator = (const VectorReference &r);
140 
145  const VectorReference &
146  operator = (const TrilinosScalar &s) const;
147 
153  const VectorReference &
154  operator += (const TrilinosScalar &s) const;
155 
161  const VectorReference &
162  operator -= (const TrilinosScalar &s) const;
163 
169  const VectorReference &
170  operator *= (const TrilinosScalar &s) const;
171 
177  const VectorReference &
178  operator /= (const TrilinosScalar &s) const;
179 
186  operator TrilinosScalar () const;
187 
191  DeclException1 (ExcTrilinosError,
192  int,
193  << "An error with error number " << arg1
194  << " occurred while calling a Trilinos function");
195 
199  DeclException3 (ExcAccessToNonLocalElement,
200  size_type, size_type, size_type,
201  << "You tried to access element " << arg1
202  << " of a distributed vector, but it is not stored on "
203  << "the current processor. Note: the elements stored "
204  << "on the current processor are within the range "
205  << arg2 << " through " << arg3
206  << " but Trilinos vectors need not store contiguous "
207  << "ranges on each processor, and not every element in "
208  << "this range may in fact be stored locally.");
209 
210  private:
215  VectorBase &vector;
216 
221  const size_type index;
222 
229  friend class ::TrilinosWrappers::VectorBase;
230  };
231  }
269  class VectorBase : public Subscriptor
270  {
271  public:
280  typedef TrilinosScalar value_type;
281  typedef TrilinosScalar real_type;
282  typedef ::types::global_dof_index size_type;
283  typedef value_type *iterator;
284  typedef const value_type *const_iterator;
285  typedef internal::VectorReference reference;
286  typedef const internal::VectorReference const_reference;
287 
292 
303  VectorBase ();
304 
311  VectorBase (const VectorBase &v);
312 
316  virtual ~VectorBase ();
317 
324  void clear ();
325 
333  void reinit (const VectorBase &v,
334  const bool fast = false);
335 
364  void compress (::VectorOperation::values operation);
365 
371 
375  void compress (const Epetra_CombineMode last_action) DEAL_II_DEPRECATED;
376 
384  bool is_compressed () const;
385 
409  VectorBase &
410  operator = (const TrilinosScalar s);
411 
420  VectorBase &
421  operator = (const VectorBase &v);
422 
440  template <typename Number>
441  VectorBase &
442  operator = (const ::Vector<Number> &v);
443 
453  bool operator == (const VectorBase &v) const;
454 
464  bool operator != (const VectorBase &v) const;
465 
470  size_type size () const;
471 
490  size_type local_size () const;
491 
522  std::pair<size_type, size_type> local_range () const;
523 
532  bool in_local_range (const size_type index) const;
533 
549 
556  bool has_ghost_elements() const;
557 
564  TrilinosScalar operator * (const VectorBase &vec) const;
565 
570  real_type norm_sqr () const;
571 
576  TrilinosScalar mean_value () const;
577 
582  TrilinosScalar minimal_value () const;
583 
588  real_type l1_norm () const;
589 
595  real_type l2_norm () const;
596 
604  real_type lp_norm (const TrilinosScalar p) const;
605 
610  real_type linfty_norm () const;
611 
622  bool all_zero () const;
623 
634  bool is_non_negative () const;
636 
637 
642 
647  reference
648  operator () (const size_type index);
649 
655  TrilinosScalar
656  operator () (const size_type index) const;
657 
664  reference
665  operator [] (const size_type index);
666 
674  TrilinosScalar
675  operator [] (const size_type index) const;
676 
687  void extract_subvector_to (const std::vector<size_type> &indices,
688  std::vector<TrilinosScalar> &values) const;
689 
694  template <typename ForwardIterator, typename OutputIterator>
695  void extract_subvector_to (ForwardIterator indices_begin,
696  const ForwardIterator indices_end,
697  OutputIterator values_begin) const;
698 
710  TrilinosScalar el (const size_type index) const;
711 
719  iterator begin ();
720 
725  const_iterator begin () const;
726 
731  iterator end ();
732 
737  const_iterator end () const;
738 
750  void set (const std::vector<size_type> &indices,
751  const std::vector<TrilinosScalar> &values);
752 
760  void set (const std::vector<size_type> &indices,
761  const ::Vector<TrilinosScalar> &values);
763 
764 
769 
780  void set (const size_type n_elements,
781  const size_type *indices,
782  const TrilinosScalar *values);
783 
792  void add (const std::vector<size_type> &indices,
793  const std::vector<TrilinosScalar> &values);
794 
802  void add (const std::vector<size_type> &indices,
803  const ::Vector<TrilinosScalar> &values);
804 
814  void add (const size_type n_elements,
815  const size_type *indices,
816  const TrilinosScalar *values);
817 
822  VectorBase &operator *= (const TrilinosScalar factor);
823 
828  VectorBase &operator /= (const TrilinosScalar factor);
829 
834  VectorBase &operator += (const VectorBase &V);
835 
840  VectorBase &operator -= (const VectorBase &V);
841 
847  void add (const TrilinosScalar s);
848 
859  void add (const VectorBase &V,
860  const bool allow_different_maps = false);
861 
867  void add (const TrilinosScalar a,
868  const VectorBase &V);
869 
875  void add (const TrilinosScalar a,
876  const VectorBase &V,
877  const TrilinosScalar b,
878  const VectorBase &W);
879 
885  void sadd (const TrilinosScalar s,
886  const VectorBase &V);
887 
893  void sadd (const TrilinosScalar s,
894  const TrilinosScalar a,
895  const VectorBase &V);
896 
900  void sadd (const TrilinosScalar s,
901  const TrilinosScalar a,
902  const VectorBase &V,
903  const TrilinosScalar b,
904  const VectorBase &W);
905 
911  void sadd (const TrilinosScalar s,
912  const TrilinosScalar a,
913  const VectorBase &V,
914  const TrilinosScalar b,
915  const VectorBase &W,
916  const TrilinosScalar c,
917  const VectorBase &X);
918 
928  void scale (const VectorBase &scaling_factors);
929 
934  void equ (const TrilinosScalar a,
935  const VectorBase &V);
936 
941  void equ (const TrilinosScalar a,
942  const VectorBase &V,
943  const TrilinosScalar b,
944  const VectorBase &W);
945 
963  void ratio (const VectorBase &a,
964  const VectorBase &b);
966 
967 
972 
978  const Epetra_MultiVector &trilinos_vector () const;
979 
985  Epetra_FEVector &trilinos_vector ();
986 
993  const Epetra_Map &vector_partitioner () const;
994 
1001  void print (const char *format = 0) const;
1002 
1016  void print (std::ostream &out,
1017  const unsigned int precision = 3,
1018  const bool scientific = true,
1019  const bool across = true) const;
1020 
1047  void swap (VectorBase &v);
1048 
1053  std::size_t memory_consumption () const;
1054 
1060  const MPI_Comm &get_mpi_communicator () const;
1062 
1067 
1071  DeclException0 (ExcDifferentParallelPartitioning);
1072 
1076  DeclException1 (ExcTrilinosError,
1077  int,
1078  << "An error with error number " << arg1
1079  << " occurred while calling a Trilinos function");
1080 
1084  DeclException3 (ExcAccessToNonlocalElement,
1085  size_type, size_type, size_type,
1086  << "You tried to access element " << arg1
1087  << " of a distributed vector, but only entries "
1088  << arg2 << " through " << arg3
1089  << " are stored locally and can be accessed.");
1090 
1091 
1092  private:
1118  Epetra_CombineMode last_action;
1119 
1126 
1133 
1139  std_cxx1x::shared_ptr<Epetra_FEVector> vector;
1140 
1141 
1146  friend class internal::VectorReference;
1147  friend class Vector;
1148  friend class MPI::Vector;
1149  };
1150 
1151 
1152 
1153 
1154 // ------------------- inline and template functions --------------
1155 
1164  inline
1166  {
1167  u.swap (v);
1168  }
1169 
1170 
1171 #ifndef DOXYGEN
1172 
1173  namespace internal
1174  {
1175  inline
1176  VectorReference::VectorReference (VectorBase &vector,
1177  const size_type index)
1178  :
1179  vector (vector),
1180  index (index)
1181  {}
1182 
1183 
1184  inline
1185  const VectorReference &
1186  VectorReference::operator = (const VectorReference &r) const
1187  {
1188  // as explained in the class
1189  // documentation, this is not the copy
1190  // operator. so simply pass on to the
1191  // "correct" assignment operator
1192  *this = static_cast<TrilinosScalar> (r);
1193 
1194  return *this;
1195  }
1196 
1197 
1198 
1199  inline
1200  const VectorReference &
1201  VectorReference::operator = (const VectorReference &r)
1202  {
1203  // as above
1204  *this = static_cast<TrilinosScalar> (r);
1205 
1206  return *this;
1207  }
1208 
1209 
1210  inline
1211  const VectorReference &
1212  VectorReference::operator = (const TrilinosScalar &value) const
1213  {
1214  vector.set (1, &index, &value);
1215  return *this;
1216  }
1217 
1218 
1219 
1220  inline
1221  const VectorReference &
1222  VectorReference::operator += (const TrilinosScalar &value) const
1223  {
1224  vector.add (1, &index, &value);
1225  return *this;
1226  }
1227 
1228 
1229 
1230  inline
1231  const VectorReference &
1232  VectorReference::operator -= (const TrilinosScalar &value) const
1233  {
1234  TrilinosScalar new_value = -value;
1235  vector.add (1, &index, &new_value);
1236  return *this;
1237  }
1238 
1239 
1240 
1241  inline
1242  const VectorReference &
1243  VectorReference::operator *= (const TrilinosScalar &value) const
1244  {
1245  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1246  vector.set (1, &index, &new_value);
1247  return *this;
1248  }
1249 
1250 
1251 
1252  inline
1253  const VectorReference &
1254  VectorReference::operator /= (const TrilinosScalar &value) const
1255  {
1256  TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1257  vector.set (1, &index, &new_value);
1258  return *this;
1259  }
1260  }
1261 
1262 
1263 
1264  inline
1265  bool
1266  VectorBase::is_compressed () const
1267  {
1268  return compressed;
1269  }
1270 
1271 
1272 
1273  inline
1274  bool
1275  VectorBase::in_local_range (const size_type index) const
1276  {
1277  std::pair<size_type, size_type> range = local_range();
1278 
1279  return ((index >= range.first) && (index < range.second));
1280  }
1281 
1282 
1283 
1284  inline
1285  IndexSet
1286  VectorBase::locally_owned_elements() const
1287  {
1288  IndexSet is (size());
1289 
1290  // easy case: local range is contiguous
1291  if (vector->Map().LinearMap())
1292  {
1293  const std::pair<size_type, size_type> x = local_range();
1294  is.add_range (x.first, x.second);
1295  }
1296  else if (vector->Map().NumMyElements() > 0)
1297  {
1298  const size_type n_indices = vector->Map().NumMyElements();
1299 #ifndef DEAL_II_USE_LARGE_INDEX_TYPE
1300  unsigned int *vector_indices = (unsigned int *)vector->Map().MyGlobalElements();
1301 #else
1302  size_type *vector_indices = (size_type *)vector->Map().MyGlobalElements64();
1303 #endif
1304  is.add_indices(vector_indices, vector_indices+n_indices);
1305  is.compress();
1306  }
1307 
1308  return is;
1309  }
1310 
1311 
1312 
1313  inline
1314  bool
1315  VectorBase::has_ghost_elements() const
1316  {
1317  return has_ghosts;
1318  }
1319 
1320 
1321 
1322  inline
1323  internal::VectorReference
1324  VectorBase::operator () (const size_type index)
1325  {
1326  return internal::VectorReference (*this, index);
1327  }
1328 
1329 
1330 
1331  inline
1332  internal::VectorReference
1333  VectorBase::operator [] (const size_type index)
1334  {
1335  return operator() (index);
1336  }
1337 
1338 
1339  inline
1340  TrilinosScalar
1341  VectorBase::operator [] (const size_type index) const
1342  {
1343  return operator() (index);
1344  }
1345 
1346 
1347 
1348  inline
1349  void VectorBase::extract_subvector_to (const std::vector<size_type> &indices,
1350  std::vector<TrilinosScalar> &values) const
1351  {
1352  for (size_type i = 0; i < indices.size(); ++i)
1353  values[i] = operator()(indices[i]);
1354  }
1355 
1356 
1357 
1358  template <typename ForwardIterator, typename OutputIterator>
1359  inline
1360  void VectorBase::extract_subvector_to (ForwardIterator indices_begin,
1361  const ForwardIterator indices_end,
1362  OutputIterator values_begin) const
1363  {
1364  while (indices_begin != indices_end)
1365  {
1366  *values_begin = operator()(*indices_begin);
1367  indices_begin++;
1368  values_begin++;
1369  }
1370  }
1371 
1372 
1373 
1374  inline
1375  VectorBase::iterator
1376  VectorBase::begin()
1377  {
1378  return (*vector)[0];
1379  }
1380 
1381 
1382 
1383  inline
1384  VectorBase::iterator
1385  VectorBase::end()
1386  {
1387  return (*vector)[0]+local_size();
1388  }
1389 
1390 
1391 
1392  inline
1393  VectorBase::const_iterator
1394  VectorBase::begin() const
1395  {
1396  return (*vector)[0];
1397  }
1398 
1399 
1400 
1401  inline
1402  VectorBase::const_iterator
1403  VectorBase::end() const
1404  {
1405  return (*vector)[0]+local_size();
1406  }
1407 
1408 
1409 
1410  inline
1411  void
1412  VectorBase::reinit (const VectorBase &v,
1413  const bool fast)
1414  {
1415  Assert (vector.get() != 0,
1416  ExcMessage("Vector has not been constructed properly."));
1417 
1418  if (fast == false ||
1419  vector_partitioner().SameAs(v.vector_partitioner())==false)
1420  vector.reset (new Epetra_FEVector(*v.vector));
1421  }
1422 
1423 
1424 
1425  inline
1426  void
1427  VectorBase::compress (const Epetra_CombineMode last_action)
1428  {
1429  ::VectorOperation::values last_action_ =
1430  ::VectorOperation::unknown;
1431  if (last_action == Add)
1432  last_action_ = ::VectorOperation::add;
1433  else if (last_action == Insert)
1434  last_action_ = ::VectorOperation::insert;
1435  else
1436  AssertThrow(false, ExcNotImplemented());
1437 
1438  compress(last_action_);
1439  }
1440 
1441 
1442 
1443  inline
1444  void
1445  VectorBase::compress (::VectorOperation::values given_last_action)
1446  {
1447  //Select which mode to send to
1448  //Trilinos. Note that we use last_action
1449  //if available and ignore what the user
1450  //tells us to detect wrongly mixed
1451  //operations. Typically given_last_action
1452  //is only used on machines that do not
1453  //execute an operation (because they have
1454  //no own cells for example).
1455  Epetra_CombineMode mode = last_action;
1456  if (last_action == Zero)
1457  {
1458  if (given_last_action==::VectorOperation::add)
1459  mode = Add;
1460  else if (given_last_action==::VectorOperation::insert)
1461  mode = Insert;
1462  }
1463 
1464 #ifdef DEBUG
1465 # ifdef DEAL_II_WITH_MPI
1466  // check that every process has decided
1467  // to use the same mode. This will
1468  // otherwise result in undefined
1469  // behaviour in the call to
1470  // GlobalAssemble().
1471  double double_mode = mode;
1473  = Utilities::MPI::min_max_avg (double_mode,
1474  dynamic_cast<const Epetra_MpiComm *>
1475  (&vector_partitioner().Comm())->GetMpiComm());
1476  Assert(result.max-result.min<1e-5,
1477  ExcMessage ("Not all processors agree whether the last operation on "
1478  "this vector was an addition or a set operation. This will "
1479  "prevent the compress() operation from succeeding."));
1480 
1481 # endif
1482 #endif
1483 
1484  // Now pass over the information about
1485  // what we did last to the vector.
1486  const int ierr = vector->GlobalAssemble(mode);
1487  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1488  last_action = Zero;
1489 
1490  compressed = true;
1491  }
1492 
1493 
1494 
1495  inline
1496  void
1497  VectorBase::compress ()
1498  {
1499  compress(VectorOperation::unknown);
1500  }
1501 
1502 
1503 
1504  inline
1505  VectorBase &
1506  VectorBase::operator = (const TrilinosScalar s)
1507  {
1508  // if we have ghost values, do not allow
1509  // writing to this vector at all.
1510  Assert (!has_ghost_elements(), ExcGhostsPresent());
1511 
1513 
1514  const int ierr = vector->PutScalar(s);
1515 
1516  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1517 
1518  return *this;
1519  }
1520 
1521 
1522 
1523  inline
1524  void
1525  VectorBase::set (const std::vector<size_type> &indices,
1526  const std::vector<TrilinosScalar> &values)
1527  {
1528  // if we have ghost values, do not allow
1529  // writing to this vector at all.
1530  Assert (!has_ghost_elements(), ExcGhostsPresent());
1531 
1532  Assert (indices.size() == values.size(),
1533  ExcDimensionMismatch(indices.size(),values.size()));
1534 
1535  set (indices.size(), &indices[0], &values[0]);
1536  }
1537 
1538 
1539 
1540  inline
1541  void
1542  VectorBase::set (const std::vector<size_type> &indices,
1543  const ::Vector<TrilinosScalar> &values)
1544  {
1545  // if we have ghost values, do not allow
1546  // writing to this vector at all.
1547  Assert (!has_ghost_elements(), ExcGhostsPresent());
1548 
1549  Assert (indices.size() == values.size(),
1550  ExcDimensionMismatch(indices.size(),values.size()));
1551 
1552  set (indices.size(), &indices[0], values.begin());
1553  }
1554 
1555 
1556 
1557  inline
1558  void
1559  VectorBase::set (const size_type n_elements,
1560  const size_type *indices,
1561  const TrilinosScalar *values)
1562  {
1563  // if we have ghost values, do not allow
1564  // writing to this vector at all.
1565  Assert (!has_ghost_elements(), ExcGhostsPresent());
1566 
1567  if (last_action == Add)
1568  vector->GlobalAssemble(Add);
1569 
1570  if (last_action != Insert)
1571  last_action = Insert;
1572 
1573  for (size_type i=0; i<n_elements; ++i)
1574  {
1575  const size_type row = indices[i];
1576  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1577  if (local_row == -1)
1578  {
1579  const int ierr = vector->ReplaceGlobalValues (1,
1580  (const TrilinosWrappers::types::int_type *)(&row),
1581  &values[i]);
1582  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1583  compressed = false;
1584  }
1585  else
1586  (*vector)[0][local_row] = values[i];
1587  }
1588  }
1589 
1590 
1591 
1592  inline
1593  void
1594  VectorBase::add (const std::vector<size_type> &indices,
1595  const std::vector<TrilinosScalar> &values)
1596  {
1597  // if we have ghost values, do not allow
1598  // writing to this vector at all.
1599  Assert (!has_ghost_elements(), ExcGhostsPresent());
1600  Assert (indices.size() == values.size(),
1601  ExcDimensionMismatch(indices.size(),values.size()));
1602 
1603  add (indices.size(), &indices[0], &values[0]);
1604  }
1605 
1606 
1607 
1608  inline
1609  void
1610  VectorBase::add (const std::vector<size_type> &indices,
1611  const ::Vector<TrilinosScalar> &values)
1612  {
1613  // if we have ghost values, do not allow
1614  // writing to this vector at all.
1615  Assert (!has_ghost_elements(), ExcGhostsPresent());
1616  Assert (indices.size() == values.size(),
1617  ExcDimensionMismatch(indices.size(),values.size()));
1618 
1619  add (indices.size(), &indices[0], values.begin());
1620  }
1621 
1622 
1623 
1624  inline
1625  void
1626  VectorBase::add (const size_type n_elements,
1627  const size_type *indices,
1628  const TrilinosScalar *values)
1629  {
1630  // if we have ghost values, do not allow
1631  // writing to this vector at all.
1632  Assert (!has_ghost_elements(), ExcGhostsPresent());
1633 
1634  if (last_action != Add)
1635  {
1636  if (last_action == Insert)
1637  vector->GlobalAssemble(Insert);
1638  last_action = Add;
1639  }
1640 
1641  for (size_type i=0; i<n_elements; ++i)
1642  {
1643  const size_type row = indices[i];
1644  const TrilinosWrappers::types::int_type local_row = vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(row));
1645  if (local_row == -1)
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  (*vector)[0][local_row] += values[i];
1655  }
1656  }
1657 
1658 
1659 
1660  inline
1661  VectorBase::size_type
1662  VectorBase::size () const
1663  {
1664 #ifndef DEAL_II_USE_LARGE_INDEX_TYPE
1665  return (size_type) (vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID());
1666 #else
1667  return (size_type) (vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64());
1668 #endif
1669  }
1670 
1671 
1672 
1673  inline
1674  VectorBase::size_type
1675  VectorBase::local_size () const
1676  {
1677  return (size_type) vector->Map().NumMyElements();
1678  }
1679 
1680 
1681 
1682  inline
1683  std::pair<VectorBase::size_type, VectorBase::size_type>
1684  VectorBase::local_range () const
1685  {
1686 #ifndef DEAL_II_USE_LARGE_INDEX_TYPE
1687  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1688  const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID()+1;
1689 #else
1690  const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID64();
1691  const TrilinosWrappers::types::int_type end = vector->Map().MaxMyGID64()+1;
1692 #endif
1693 
1694  Assert (end-begin == vector->Map().NumMyElements(),
1695  ExcMessage ("This function only makes sense if the elements that this "
1696  "vector stores on the current processor form a contiguous range. "
1697  "This does not appear to be the case for the current vector."));
1698 
1699  return std::make_pair (begin, end);
1700  }
1701 
1702 
1703 
1704  inline
1705  TrilinosScalar
1706  VectorBase::operator * (const VectorBase &vec) const
1707  {
1708  Assert (vector->Map().SameAs(vec.vector->Map()),
1709  ExcDifferentParallelPartitioning());
1710  Assert (!has_ghost_elements(), ExcGhostsPresent());
1711 
1712  TrilinosScalar result;
1713 
1714  const int ierr = vector->Dot(*(vec.vector), &result);
1715  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1716 
1717  return result;
1718  }
1719 
1720 
1721 
1722  inline
1723  VectorBase::real_type
1724  VectorBase::norm_sqr () const
1725  {
1726  const TrilinosScalar d = l2_norm();
1727  return d*d;
1728  }
1729 
1730 
1731 
1732  inline
1733  TrilinosScalar
1734  VectorBase::mean_value () const
1735  {
1736  Assert (!has_ghost_elements(), ExcGhostsPresent());
1737 
1738  TrilinosScalar mean;
1739  const int ierr = vector->MeanValue (&mean);
1740  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1741 
1742  return mean;
1743  }
1744 
1745 
1746 
1747  inline
1748  TrilinosScalar
1749  VectorBase::minimal_value () const
1750  {
1751  TrilinosScalar min_value;
1752  const int ierr = vector->MinValue (&min_value);
1753  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1754 
1755  return min_value;
1756  }
1757 
1758 
1759 
1760  inline
1761  VectorBase::real_type
1762  VectorBase::l1_norm () const
1763  {
1764  Assert (!has_ghost_elements(), ExcGhostsPresent());
1765 
1766  TrilinosScalar d;
1767  const int ierr = vector->Norm1 (&d);
1768  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1769 
1770  return d;
1771  }
1772 
1773 
1774 
1775  inline
1776  VectorBase::real_type
1777  VectorBase::l2_norm () const
1778  {
1779  Assert (!has_ghost_elements(), ExcGhostsPresent());
1780 
1781  TrilinosScalar d;
1782  const int ierr = vector->Norm2 (&d);
1783  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1784 
1785  return d;
1786  }
1787 
1788 
1789 
1790  inline
1791  VectorBase::real_type
1792  VectorBase::lp_norm (const TrilinosScalar p) const
1793  {
1794  Assert (!has_ghost_elements(), ExcGhostsPresent());
1795 
1796  TrilinosScalar norm = 0;
1797  TrilinosScalar sum=0;
1798  const size_type n_local = local_size();
1799 
1800  // loop over all the elements because
1801  // Trilinos does not support lp norms
1802  for (size_type i=0; i<n_local; i++)
1803  sum += std::pow(std::fabs((*vector)[0][i]), p);
1804 
1805  norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
1806 
1807  return norm;
1808  }
1809 
1810 
1811 
1812  inline
1813  VectorBase::real_type
1814  VectorBase::linfty_norm () const
1815  {
1816  // while we disallow the other
1817  // norm operations on ghosted
1818  // vectors, this particular norm
1819  // is safe to run even in the
1820  // presence of ghost elements
1821  TrilinosScalar d;
1822  const int ierr = vector->NormInf (&d);
1823  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1824 
1825  return d;
1826  }
1827 
1828 
1829 
1830  // inline also scalar products, vector
1831  // additions etc. since they are all
1832  // representable by a single Trilinos
1833  // call. This reduces the overhead of the
1834  // wrapper class.
1835  inline
1836  VectorBase &
1837  VectorBase::operator *= (const TrilinosScalar a)
1838  {
1840 
1841  const int ierr = vector->Scale(a);
1842  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1843 
1844  return *this;
1845  }
1846 
1847 
1848 
1849  inline
1850  VectorBase &
1851  VectorBase::operator /= (const TrilinosScalar a)
1852  {
1854 
1855  const TrilinosScalar factor = 1./a;
1856 
1858 
1859  const int ierr = vector->Scale(factor);
1860  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1861 
1862  return *this;
1863  }
1864 
1865 
1866 
1867  inline
1868  VectorBase &
1869  VectorBase::operator += (const VectorBase &v)
1870  {
1871  Assert (size() == v.size(),
1872  ExcDimensionMismatch(size(), v.size()));
1873  Assert (vector->Map().SameAs(v.vector->Map()),
1874  ExcDifferentParallelPartitioning());
1875 
1876  const int ierr = vector->Update (1.0, *(v.vector), 1.0);
1877  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1878 
1879  return *this;
1880  }
1881 
1882 
1883 
1884  inline
1885  VectorBase &
1886  VectorBase::operator -= (const VectorBase &v)
1887  {
1888  Assert (size() == v.size(),
1889  ExcDimensionMismatch(size(), v.size()));
1890  Assert (vector->Map().SameAs(v.vector->Map()),
1891  ExcDifferentParallelPartitioning());
1892 
1893  const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
1894  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1895 
1896  return *this;
1897  }
1898 
1899 
1900 
1901  inline
1902  void
1903  VectorBase::add (const TrilinosScalar s)
1904  {
1905  // if we have ghost values, do not allow
1906  // writing to this vector at all.
1907  Assert (!has_ghost_elements(), ExcGhostsPresent());
1909 
1910  size_type n_local = local_size();
1911  for (size_type i=0; i<n_local; i++)
1912  (*vector)[0][i] += s;
1913  }
1914 
1915 
1916 
1917  inline
1918  void
1919  VectorBase::add (const TrilinosScalar a,
1920  const VectorBase &v)
1921  {
1922  // if we have ghost values, do not allow
1923  // writing to this vector at all.
1924  Assert (!has_ghost_elements(), ExcGhostsPresent());
1925  Assert (local_size() == v.local_size(),
1926  ExcDimensionMismatch(local_size(), v.local_size()));
1927 
1929 
1930  const int ierr = vector->Update(a, *(v.vector), 1.);
1931  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1932  }
1933 
1934 
1935 
1936  inline
1937  void
1938  VectorBase::add (const TrilinosScalar a,
1939  const VectorBase &v,
1940  const TrilinosScalar b,
1941  const VectorBase &w)
1942  {
1943  // if we have ghost values, do not allow
1944  // writing to this vector at all.
1945  Assert (!has_ghost_elements(), ExcGhostsPresent());
1946  Assert (local_size() == v.local_size(),
1947  ExcDimensionMismatch(local_size(), v.local_size()));
1948  Assert (local_size() == w.local_size(),
1949  ExcDimensionMismatch(local_size(), w.local_size()));
1950 
1953 
1954  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
1955 
1956  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1957  }
1958 
1959 
1960 
1961  inline
1962  void
1963  VectorBase::sadd (const TrilinosScalar s,
1964  const VectorBase &v)
1965  {
1966  // if we have ghost values, do not allow
1967  // writing to this vector at all.
1968  Assert (!has_ghost_elements(), ExcGhostsPresent());
1969  Assert (local_size() == v.local_size(),
1970  ExcDimensionMismatch(local_size(), v.local_size()));
1971 
1973 
1974  const int ierr = vector->Update(1., *(v.vector), s);
1975 
1976  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1977  }
1978 
1979 
1980 
1981  inline
1982  void
1983  VectorBase::sadd (const TrilinosScalar s,
1984  const TrilinosScalar a,
1985  const VectorBase &v)
1986  {
1987  // if we have ghost values, do not allow
1988  // writing to this vector at all.
1989  Assert (!has_ghost_elements(), ExcGhostsPresent());
1990  Assert (local_size() == v.local_size(),
1991  ExcDimensionMismatch(local_size(), v.local_size()));
1992 
1995 
1996  const int ierr = vector->Update(a, *(v.vector), s);
1997 
1998  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1999  }
2000 
2001 
2002 
2003  inline
2004  void
2005  VectorBase::sadd (const TrilinosScalar s,
2006  const TrilinosScalar a,
2007  const VectorBase &v,
2008  const TrilinosScalar b,
2009  const VectorBase &w)
2010  {
2011  // if we have ghost values, do not allow
2012  // writing to this vector at all.
2013  Assert (!has_ghost_elements(), ExcGhostsPresent());
2014  Assert (local_size() == v.local_size(),
2015  ExcDimensionMismatch(local_size(), v.local_size()));
2016  Assert (local_size() == w.local_size(),
2017  ExcDimensionMismatch(local_size(), w.local_size()));
2018 
2022 
2023  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
2024 
2025  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2026  }
2027 
2028 
2029 
2030  inline
2031  void
2032  VectorBase::sadd (const TrilinosScalar s,
2033  const TrilinosScalar a,
2034  const VectorBase &v,
2035  const TrilinosScalar b,
2036  const VectorBase &w,
2037  const TrilinosScalar c,
2038  const VectorBase &x)
2039  {
2040  // if we have ghost values, do not allow
2041  // writing to this vector at all.
2042  Assert (!has_ghost_elements(), ExcGhostsPresent());
2043  Assert (local_size() == v.local_size(),
2044  ExcDimensionMismatch(local_size(), v.local_size()));
2045  Assert (local_size() == w.local_size(),
2046  ExcDimensionMismatch(local_size(), w.local_size()));
2047  Assert (local_size() == x.local_size(),
2048  ExcDimensionMismatch(local_size(), x.local_size()));
2049 
2054 
2055  // Update member can only
2056  // input two other vectors so
2057  // do it in two steps
2058  const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
2059  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2060 
2061  const int jerr = vector->Update(c, *(x.vector), 1.);
2062  Assert (jerr == 0, ExcTrilinosError(jerr));
2063  (void)jerr; // removes -Wunused-parameter warning in optimized mode
2064  }
2065 
2066 
2067 
2068  inline
2069  void
2070  VectorBase::scale (const VectorBase &factors)
2071  {
2072  // if we have ghost values, do not allow
2073  // writing to this vector at all.
2074  Assert (!has_ghost_elements(), ExcGhostsPresent());
2075  Assert (local_size() == factors.local_size(),
2076  ExcDimensionMismatch(local_size(), factors.local_size()));
2077 
2078  const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
2079  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2080  }
2081 
2082 
2083 
2084  inline
2085  void
2086  VectorBase::equ (const TrilinosScalar a,
2087  const VectorBase &v)
2088  {
2089  // if we have ghost values, do not allow
2090  // writing to this vector at all.
2091  Assert (!has_ghost_elements(), ExcGhostsPresent());
2093 
2094  // If we don't have the same map, copy.
2095  if (vector->Map().SameAs(v.vector->Map())==false)
2096  {
2097  *vector = *v.vector;
2098  *this *= a;
2099  }
2100  else
2101  {
2102  // Otherwise, just update
2103  int ierr = vector->Update(a, *v.vector, 0.0);
2104  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2105 
2106  last_action = Zero;
2107  }
2108 
2109  }
2110 
2111 
2112 
2113  inline
2114  void
2115  VectorBase::equ (const TrilinosScalar a,
2116  const VectorBase &v,
2117  const TrilinosScalar b,
2118  const VectorBase &w)
2119  {
2120  // if we have ghost values, do not allow
2121  // writing to this vector at all.
2122  Assert (!has_ghost_elements(), ExcGhostsPresent());
2123  Assert (v.local_size() == w.local_size(),
2124  ExcDimensionMismatch (v.local_size(), w.local_size()));
2125 
2128 
2129  // If we don't have the same map, copy.
2130  if (vector->Map().SameAs(v.vector->Map())==false)
2131  {
2132  *vector = *v.vector;
2133  sadd(a, b, w);
2134  }
2135  else
2136  {
2137  // Otherwise, just update. verify
2138  // that *this does not only have
2139  // the same map as v (the
2140  // if-condition above) but also as
2141  // w
2142  Assert (vector->Map().SameAs(w.vector->Map()),
2143  ExcDifferentParallelPartitioning());
2144  int ierr = vector->Update(a, *v.vector, b, *w.vector, 0.0);
2145  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2146 
2147  last_action = Zero;
2148  }
2149  }
2150 
2151 
2152 
2153  inline
2154  void
2155  VectorBase::ratio (const VectorBase &v,
2156  const VectorBase &w)
2157  {
2158  Assert (v.local_size() == w.local_size(),
2159  ExcDimensionMismatch (v.local_size(), w.local_size()));
2160  Assert (local_size() == w.local_size(),
2161  ExcDimensionMismatch (local_size(), w.local_size()));
2162 
2163  const int ierr = vector->ReciprocalMultiply(1.0, *(w.vector),
2164  *(v.vector), 0.0);
2165 
2166  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
2167  }
2168 
2169 
2170 
2171  inline
2172  const Epetra_MultiVector &
2173  VectorBase::trilinos_vector () const
2174  {
2175  return static_cast<const Epetra_MultiVector &>(*vector);
2176  }
2177 
2178 
2179 
2180  inline
2181  Epetra_FEVector &
2182  VectorBase::trilinos_vector ()
2183  {
2184  return *vector;
2185  }
2186 
2187 
2188 
2189  inline
2190  const Epetra_Map &
2191  VectorBase::vector_partitioner () const
2192  {
2193  return static_cast<const Epetra_Map &>(vector->Map());
2194  }
2195 
2196 
2197 
2198  inline
2199  const MPI_Comm &
2200  VectorBase::get_mpi_communicator () const
2201  {
2202  static MPI_Comm comm;
2203 
2204 #ifdef DEAL_II_WITH_MPI
2205 
2206  const Epetra_MpiComm *mpi_comm
2207  = dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2208  comm = mpi_comm->Comm();
2209 
2210 #else
2211 
2212  comm = MPI_COMM_SELF;
2213 
2214 #endif
2215 
2216  return comm;
2217  }
2218 
2219 
2220 
2221 #endif // DOXYGEN
2222 
2223 }
2224 
2227 DEAL_II_NAMESPACE_CLOSE
2228 
2229 #endif // DEAL_II_WITH_TRILINOS
2230 
2231 /*---------------------------- trilinos_vector_base.h ---------------------------*/
2232 
2233 #endif
2234 /*---------------------------- trilinos_vector_base.h ---------------------------*/
DeclException0(ExcGhostsPresent)
size_type local_size() const
void reinit(const VectorBase &v, const bool fast=false)
const MPI_Comm & get_mpi_communicator() const
real_type l2_norm() const
IndexSet locally_owned_elements() const
void sadd(const TrilinosScalar s, const VectorBase &V)
::ExceptionBase & ExcMessage(std::string arg1)
void scale(const VectorBase &scaling_factors)
STL namespace.
void swap(VectorBase &v)
real_type l1_norm() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
TrilinosScalar minimal_value() const
bool is_finite(const double x)
void set(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
std::size_t memory_consumption() const
void add(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
bool has_ghost_elements() const
::ExceptionBase & ExcGhostsPresent()
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
void print(const char *format=0) const
T sum(const T &t, const MPI_Comm &mpi_communicator)
Definition: mpi.h:447
void ratio(const VectorBase &a, const VectorBase &b)
#define Assert(cond, exc)
Definition: exceptions.h:299
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
const Epetra_Map & vector_partitioner() const
void compress() DEAL_II_DEPRECATED
double norm(const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Du)
Definition: divergence.h:591
std_cxx1x::shared_ptr< Epetra_FEVector > vector
TrilinosScalar mean_value() const
::ExceptionBase & ExcNumberNotFinite()
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< TrilinosScalar > &values) const
::ExceptionBase & ExcNotImplemented()
real_type linfty_norm() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:535
const Epetra_MultiVector & trilinos_vector() const
bool in_local_range(const size_type index) const
DeclException3(ExcAccessToNonlocalElement, size_type, size_type, size_type,<< "You tried to access element "<< arg1<< " of a distributed vector, but only entries "<< arg2<< " through "<< arg3<< " are stored locally and can be accessed.")
friend class internal::VectorReference
MinMaxAvg min_max_avg(const double my_value, const MPI_Comm &mpi_communicator)
real_type lp_norm(const TrilinosScalar p) const
DeclException1(ExcTrilinosError, int,<< "An error with error number "<< arg1<< " occurred while calling a Trilinos function")
std::pair< size_type, size_type > local_range() const
size_type size() const
TrilinosScalar el(const size_type index) const
real_type norm_sqr() const
void equ(const TrilinosScalar a, const VectorBase &V)