Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
trilinos_vector.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2008 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_trilinos_vector_h
17#define dealii_trilinos_vector_h
18
19
20#include <deal.II/base/config.h>
21
22#ifdef DEAL_II_WITH_TRILINOS
27
29# include <deal.II/lac/vector.h>
32
33# include <Epetra_ConfigDefs.h>
34# include <Epetra_FEVector.h>
35# include <Epetra_LocalMap.h>
36# include <Epetra_Map.h>
37# include <Epetra_MpiComm.h>
38
39# include <memory>
40# include <utility>
41# include <vector>
42
44
45// Forward declarations
46# ifndef DOXYGEN
47namespace LinearAlgebra
48{
49 // Forward declaration
50 template <typename Number>
51 class ReadWriteVector;
52} // namespace LinearAlgebra
53# endif
54
65namespace TrilinosWrappers
66{
67 class SparseMatrix;
68
79 namespace internal
80 {
85
95 class VectorReference
96 {
97 private:
102 VectorReference(MPI::Vector &vector, const size_type index);
103
104 public:
108 VectorReference(const VectorReference &) = default;
109
121 const VectorReference &
122 operator=(const VectorReference &r) const;
123
127 VectorReference &
128 operator=(const VectorReference &r);
129
133 const VectorReference &
134 operator=(const TrilinosScalar &s) const;
135
139 const VectorReference &
140 operator+=(const TrilinosScalar &s) const;
141
145 const VectorReference &
146 operator-=(const TrilinosScalar &s) const;
147
151 const VectorReference &
152 operator*=(const TrilinosScalar &s) const;
153
157 const VectorReference &
158 operator/=(const TrilinosScalar &s) const;
159
164 operator TrilinosScalar() const;
165
170 int,
171 << "An error with error number " << arg1
172 << " occurred while calling a Trilinos function");
173
174 private:
178 MPI::Vector &vector;
179
183 const size_type index;
184
185 // Make the vector class a friend, so that it can create objects of the
186 // present type.
187 friend class ::TrilinosWrappers::MPI::Vector;
188 };
189 } // namespace internal
194# ifndef DEAL_II_WITH_64BIT_INDICES
195 // define a helper function that queries the global ID of local ID of
196 // an Epetra_BlockMap object by calling either the 32- or 64-bit
197 // function necessary.
198 inline int
199 gid(const Epetra_BlockMap &map, int i)
200 {
201 return map.GID(i);
202 }
203# else
204 // define a helper function that queries the global ID of local ID of
205 // an Epetra_BlockMap object by calling either the 32- or 64-bit
206 // function necessary.
207 inline long long int
208 gid(const Epetra_BlockMap &map, int i)
209 {
210 return map.GID64(i);
211 }
212# endif
213
219 namespace MPI
220 {
221 class BlockVector;
222
397 class Vector : public Subscriptor
398 {
399 public:
409 using const_iterator = const value_type *;
412
422 Vector();
423
427 Vector(const Vector &v);
428
447 explicit Vector(const IndexSet &parallel_partitioning,
448 const MPI_Comm communicator = MPI_COMM_WORLD);
449
461 Vector(const IndexSet &local,
462 const IndexSet &ghost,
463 const MPI_Comm communicator = MPI_COMM_WORLD);
464
479 Vector(const IndexSet &parallel_partitioning,
480 const Vector & v,
481 const MPI_Comm communicator = MPI_COMM_WORLD);
482
495 template <typename Number>
496 Vector(const IndexSet & parallel_partitioning,
497 const ::Vector<Number> &v,
498 const MPI_Comm communicator = MPI_COMM_WORLD);
499
508 Vector(Vector &&v); // NOLINT
509
513 ~Vector() override = default;
514
519 void
520 clear();
521
545 void
546 reinit(const Vector &v,
547 const bool omit_zeroing_entries = false,
548 const bool allow_different_maps = false);
549
572 void
573 reinit(const IndexSet &parallel_partitioning,
574 const MPI_Comm communicator = MPI_COMM_WORLD,
575 const bool omit_zeroing_entries = false);
576
611 void
612 reinit(const IndexSet &locally_owned_entries,
613 const IndexSet &locally_relevant_or_ghost_entries,
614 const MPI_Comm communicator = MPI_COMM_WORLD,
615 const bool vector_writable = false);
616
627 void
628 reinit(
629 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
630 const bool make_ghosted = true,
631 const bool vector_writable = false);
632
636 void
637 reinit(const BlockVector &v, const bool import_data = false);
638
655 void
657
670 Vector &
672
678 Vector &
679 operator=(const Vector &v);
680
685 Vector &
686 operator=(Vector &&v) noexcept;
687
695 template <typename Number>
696 Vector &
697 operator=(const ::Vector<Number> &v);
698
716 void
718 const ::TrilinosWrappers::SparseMatrix &matrix,
719 const Vector & vector);
720
727 void
729 const VectorOperation::values operation);
730
734 DEAL_II_DEPRECATED_EARLY
735 void
737 const VectorOperation::values operation)
738 {
739 import_elements(rwv, operation);
740 }
741
742
748 bool
749 operator==(const Vector &v) const;
750
756 bool
757 operator!=(const Vector &v) const;
758
763 size() const;
764
780 local_size() const;
781
788
810 std::pair<size_type, size_type>
811 local_range() const;
812
820 bool
821 in_local_range(const size_type index) const;
822
838
846 bool
848
855 void
857
863 operator*(const Vector &vec) const;
864
869 norm_sqr() const;
870
875 mean_value() const;
876
881 min() const;
882
887 max() const;
888
893 l1_norm() const;
894
900 l2_norm() const;
901
907 lp_norm(const TrilinosScalar p) const;
908
913 linfty_norm() const;
914
935 add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
936
942 bool
943 all_zero() const;
944
950 bool
951 is_non_negative() const;
968 operator()(const size_type index);
969
978 operator()(const size_type index) const;
979
986 operator[](const size_type index);
987
994 operator[](const size_type index) const;
995
1011 void
1012 extract_subvector_to(const std::vector<size_type> &indices,
1013 std::vector<TrilinosScalar> & values) const;
1014
1042 template <typename ForwardIterator, typename OutputIterator>
1043 void
1044 extract_subvector_to(ForwardIterator indices_begin,
1045 const ForwardIterator indices_end,
1046 OutputIterator values_begin) const;
1047
1058 iterator
1060
1066 begin() const;
1067
1072 iterator
1074
1080 end() const;
1081
1096 void
1097 set(const std::vector<size_type> & indices,
1098 const std::vector<TrilinosScalar> &values);
1099
1104 void
1105 set(const std::vector<size_type> & indices,
1106 const ::Vector<TrilinosScalar> &values);
1107
1113 void
1114 set(const size_type n_elements,
1115 const size_type * indices,
1116 const TrilinosScalar *values);
1117
1122 void
1123 add(const std::vector<size_type> & indices,
1124 const std::vector<TrilinosScalar> &values);
1125
1130 void
1131 add(const std::vector<size_type> & indices,
1132 const ::Vector<TrilinosScalar> &values);
1133
1139 void
1140 add(const size_type n_elements,
1141 const size_type * indices,
1142 const TrilinosScalar *values);
1143
1147 Vector &
1149
1153 Vector &
1155
1159 Vector &
1161
1165 Vector &
1167
1172 void
1174
1187 void
1188 add(const Vector &V, const bool allow_different_maps = false);
1189
1193 void
1194 add(const TrilinosScalar a, const Vector &V);
1195
1199 void
1201 const Vector & V,
1202 const TrilinosScalar b,
1203 const Vector & W);
1204
1209 void
1210 sadd(const TrilinosScalar s, const Vector &V);
1211
1215 void
1216 sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1217
1223 void
1224 scale(const Vector &scaling_factors);
1225
1229 void
1230 equ(const TrilinosScalar a, const Vector &V);
1242 const Epetra_MultiVector &
1244
1249 Epetra_FEVector &
1251
1256 const Epetra_BlockMap &
1258
1266 void
1267 print(std::ostream & out,
1268 const unsigned int precision = 3,
1269 const bool scientific = true,
1270 const bool across = true) const;
1271
1285 void
1286 swap(Vector &v);
1287
1291 std::size_t
1292 memory_consumption() const;
1293
1297 MPI_Comm
1305
1310 int,
1311 << "An error with error number " << arg1
1312 << " occurred while calling a Trilinos function");
1313
1319 size_type,
1320 size_type,
1321 size_type,
1322 size_type,
1323 << "You are trying to access element " << arg1
1324 << " of a distributed vector, but this element is not stored "
1325 << "on the current processor. Note: There are " << arg2
1326 << " elements stored "
1327 << "on the current processor from within the range [" << arg3 << ','
1328 << arg4 << "] but Trilinos vectors need not store contiguous "
1329 << "ranges on each processor, and not every element in "
1330 << "this range may in fact be stored locally."
1331 << "\n\n"
1332 << "A common source for this kind of problem is that you "
1333 << "are passing a 'fully distributed' vector into a function "
1334 << "that needs read access to vector elements that correspond "
1335 << "to degrees of freedom on ghost cells (or at least to "
1336 << "'locally active' degrees of freedom that are not also "
1337 << "'locally owned'). You need to pass a vector that has these "
1338 << "elements as ghost entries.");
1339
1340 private:
1352 Epetra_CombineMode last_action;
1353
1359
1365
1371 std::unique_ptr<Epetra_FEVector> vector;
1372
1378 std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1379
1384
1385 // Make the reference class a friend.
1387 };
1388
1389
1390
1391 // ------------------- inline and template functions --------------
1392
1393
1401 inline void
1403 {
1404 u.swap(v);
1405 }
1406 } // namespace MPI
1407
1408# ifndef DOXYGEN
1409
1410 namespace internal
1411 {
1412 inline VectorReference::VectorReference(MPI::Vector & vector,
1413 const size_type index)
1414 : vector(vector)
1415 , index(index)
1416 {}
1417
1418
1419 inline const VectorReference &
1420 VectorReference::operator=(const VectorReference &r) const
1421 {
1422 // as explained in the class
1423 // documentation, this is not the copy
1424 // operator. so simply pass on to the
1425 // "correct" assignment operator
1426 *this = static_cast<TrilinosScalar>(r);
1427
1428 return *this;
1429 }
1430
1431
1432
1433 inline VectorReference &
1434 VectorReference::operator=(const VectorReference &r)
1435 {
1436 // as above
1437 *this = static_cast<TrilinosScalar>(r);
1438
1439 return *this;
1440 }
1441
1442
1443 inline const VectorReference &
1444 VectorReference::operator=(const TrilinosScalar &value) const
1445 {
1446 vector.set(1, &index, &value);
1447 return *this;
1448 }
1449
1450
1451
1452 inline const VectorReference &
1453 VectorReference::operator+=(const TrilinosScalar &value) const
1454 {
1455 vector.add(1, &index, &value);
1456 return *this;
1457 }
1458
1459
1460
1461 inline const VectorReference &
1462 VectorReference::operator-=(const TrilinosScalar &value) const
1463 {
1464 TrilinosScalar new_value = -value;
1465 vector.add(1, &index, &new_value);
1466 return *this;
1467 }
1468
1469
1470
1471 inline const VectorReference &
1472 VectorReference::operator*=(const TrilinosScalar &value) const
1473 {
1474 TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1475 vector.set(1, &index, &new_value);
1476 return *this;
1477 }
1478
1479
1480
1481 inline const VectorReference &
1482 VectorReference::operator/=(const TrilinosScalar &value) const
1483 {
1484 TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1485 vector.set(1, &index, &new_value);
1486 return *this;
1487 }
1488 } // namespace internal
1489
1490 namespace MPI
1491 {
1492 inline bool
1493 Vector::in_local_range(const size_type index) const
1494 {
1495 std::pair<size_type, size_type> range = local_range();
1496
1497 return ((index >= range.first) && (index < range.second));
1498 }
1499
1500
1501
1502 inline IndexSet
1504 {
1506 ExcMessage(
1507 "The locally owned elements have not been properly initialized!"
1508 " This happens for example if this object has been initialized"
1509 " with exactly one overlapping IndexSet."));
1510 return owned_elements;
1511 }
1512
1513
1514
1515 inline bool
1517 {
1518 return has_ghosts;
1519 }
1520
1521
1522
1523 inline void
1525 {}
1526
1527
1528
1529 inline internal::VectorReference
1530 Vector::operator()(const size_type index)
1531 {
1532 return internal::VectorReference(*this, index);
1533 }
1534
1535
1536
1537 inline internal::VectorReference
1538 Vector::operator[](const size_type index)
1539 {
1540 return operator()(index);
1541 }
1542
1543
1544
1545 inline TrilinosScalar
1546 Vector::operator[](const size_type index) const
1547 {
1548 return operator()(index);
1549 }
1550
1551
1552
1553 inline void
1554 Vector::extract_subvector_to(const std::vector<size_type> &indices,
1555 std::vector<TrilinosScalar> & values) const
1556 {
1557 for (size_type i = 0; i < indices.size(); ++i)
1558 values[i] = operator()(indices[i]);
1559 }
1560
1561
1562
1563 template <typename ForwardIterator, typename OutputIterator>
1564 inline void
1565 Vector::extract_subvector_to(ForwardIterator indices_begin,
1566 const ForwardIterator indices_end,
1567 OutputIterator values_begin) const
1568 {
1569 while (indices_begin != indices_end)
1570 {
1571 *values_begin = operator()(*indices_begin);
1572 indices_begin++;
1573 values_begin++;
1574 }
1575 }
1576
1577
1578
1579 inline Vector::iterator
1581 {
1582 return (*vector)[0];
1583 }
1584
1585
1586
1587 inline Vector::iterator
1588 Vector::end()
1589 {
1590 return (*vector)[0] + locally_owned_size();
1591 }
1592
1593
1594
1596 Vector::begin() const
1597 {
1598 return (*vector)[0];
1599 }
1600
1601
1602
1604 Vector::end() const
1605 {
1606 return (*vector)[0] + locally_owned_size();
1607 }
1608
1609
1610
1611 inline void
1612 Vector::set(const std::vector<size_type> & indices,
1613 const std::vector<TrilinosScalar> &values)
1614 {
1615 // if we have ghost values, do not allow
1616 // writing to this vector at all.
1618
1619 AssertDimension(indices.size(), values.size());
1620
1621 set(indices.size(), indices.data(), values.data());
1622 }
1623
1624
1625
1626 inline void
1627 Vector::set(const std::vector<size_type> & indices,
1628 const ::Vector<TrilinosScalar> &values)
1629 {
1630 // if we have ghost values, do not allow
1631 // writing to this vector at all.
1633
1634 AssertDimension(indices.size(), values.size());
1635
1636 set(indices.size(), indices.data(), values.begin());
1637 }
1638
1639
1640
1641 inline void
1642 Vector::set(const size_type n_elements,
1643 const size_type * indices,
1644 const TrilinosScalar *values)
1645 {
1646 // if we have ghost values, do not allow
1647 // writing to this vector at all.
1649
1650 if (last_action == Add)
1651 {
1652 const int ierr = vector->GlobalAssemble(Add);
1653 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1654 }
1655
1656 if (last_action != Insert)
1657 last_action = Insert;
1658
1659 for (size_type i = 0; i < n_elements; ++i)
1660 {
1661 const TrilinosWrappers::types::int_type row = indices[i];
1662 const TrilinosWrappers::types::int_type local_row =
1663 vector->Map().LID(row);
1664 if (local_row != -1)
1665 (*vector)[0][local_row] = values[i];
1666 else
1667 {
1668 const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1669 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1670 compressed = false;
1671 }
1672 // in set operation, do not use the pre-allocated vector for nonlocal
1673 // entries even if it exists. This is to ensure that we really only
1674 // set the elements touched by the set() method and not all contained
1675 // in the nonlocal entries vector (there is no way to distinguish them
1676 // on the receiving processor)
1677 }
1678 }
1679
1680
1681
1682 inline void
1683 Vector::add(const std::vector<size_type> & indices,
1684 const std::vector<TrilinosScalar> &values)
1685 {
1686 // if we have ghost values, do not allow
1687 // writing to this vector at all.
1689 AssertDimension(indices.size(), values.size());
1690
1691 add(indices.size(), indices.data(), values.data());
1692 }
1693
1694
1695
1696 inline void
1697 Vector::add(const std::vector<size_type> & indices,
1698 const ::Vector<TrilinosScalar> &values)
1699 {
1700 // if we have ghost values, do not allow
1701 // writing to this vector at all.
1703 AssertDimension(indices.size(), values.size());
1704
1705 add(indices.size(), indices.data(), values.begin());
1706 }
1707
1708
1709
1710 inline void
1711 Vector::add(const size_type n_elements,
1712 const size_type * indices,
1713 const TrilinosScalar *values)
1714 {
1715 // if we have ghost values, do not allow
1716 // writing to this vector at all.
1718
1719 if (last_action != Add)
1720 {
1721 if (last_action == Insert)
1722 {
1723 const int ierr = vector->GlobalAssemble(Insert);
1724 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1725 }
1726 last_action = Add;
1727 }
1728
1729 for (size_type i = 0; i < n_elements; ++i)
1730 {
1731 const size_type row = indices[i];
1732 const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1733 static_cast<TrilinosWrappers::types::int_type>(row));
1734 if (local_row != -1)
1735 (*vector)[0][local_row] += values[i];
1736 else if (nonlocal_vector.get() == nullptr)
1737 {
1738 const int ierr = vector->SumIntoGlobalValues(
1739 1,
1740 reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1741 &row),
1742 &values[i]);
1743 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1744 compressed = false;
1745 }
1746 else
1747 {
1748 // use pre-allocated vector for non-local entries if it exists for
1749 // addition operation
1751 nonlocal_vector->Map().LID(
1752 static_cast<TrilinosWrappers::types::int_type>(row));
1753 Assert(my_row != -1,
1754 ExcMessage(
1755 "Attempted to write into off-processor vector entry "
1756 "that has not be specified as being writable upon "
1757 "initialization"));
1758 (*nonlocal_vector)[0][my_row] += values[i];
1759 compressed = false;
1760 }
1761 }
1762 }
1763
1764
1765
1766 inline Vector::size_type
1767 Vector::size() const
1768 {
1769# ifndef DEAL_II_WITH_64BIT_INDICES
1770 return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1771# else
1772 return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1773# endif
1774 }
1775
1776
1777
1778 inline Vector::size_type
1779 Vector::local_size() const
1780 {
1781 return vector->Map().NumMyElements();
1782 }
1783
1784
1785
1786 inline Vector::size_type
1788 {
1789 return owned_elements.n_elements();
1790 }
1791
1792
1793
1794 inline std::pair<Vector::size_type, Vector::size_type>
1795 Vector::local_range() const
1796 {
1797# ifndef DEAL_II_WITH_64BIT_INDICES
1798 const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1800 vector->Map().MaxMyGID() + 1;
1801# else
1803 vector->Map().MinMyGID64();
1805 vector->Map().MaxMyGID64() + 1;
1806# endif
1807
1808 Assert(
1809 end - begin == vector->Map().NumMyElements(),
1810 ExcMessage(
1811 "This function only makes sense if the elements that this "
1812 "vector stores on the current processor form a contiguous range. "
1813 "This does not appear to be the case for the current vector."));
1814
1815 return std::make_pair(begin, end);
1816 }
1817
1818
1819
1820 inline TrilinosScalar
1821 Vector::operator*(const Vector &vec) const
1822 {
1823 Assert(vector->Map().SameAs(vec.vector->Map()),
1826
1827 TrilinosScalar result;
1828
1829 const int ierr = vector->Dot(*(vec.vector), &result);
1830 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1831
1832 return result;
1833 }
1834
1835
1836
1837 inline Vector::real_type
1838 Vector::norm_sqr() const
1839 {
1840 const TrilinosScalar d = l2_norm();
1841 return d * d;
1842 }
1843
1844
1845
1846 inline TrilinosScalar
1847 Vector::mean_value() const
1848 {
1850
1852 const int ierr = vector->MeanValue(&mean);
1853 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1854
1855 return mean;
1856 }
1857
1858
1859
1860 inline TrilinosScalar
1861 Vector::min() const
1862 {
1863 TrilinosScalar min_value;
1864 const int ierr = vector->MinValue(&min_value);
1865 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1866
1867 return min_value;
1868 }
1869
1870
1871
1872 inline TrilinosScalar
1873 Vector::max() const
1874 {
1875 TrilinosScalar max_value;
1876 const int ierr = vector->MaxValue(&max_value);
1877 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1878
1879 return max_value;
1880 }
1881
1882
1883
1884 inline Vector::real_type
1885 Vector::l1_norm() const
1886 {
1888
1890 const int ierr = vector->Norm1(&d);
1891 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1892
1893 return d;
1894 }
1895
1896
1897
1898 inline Vector::real_type
1899 Vector::l2_norm() const
1900 {
1902
1904 const int ierr = vector->Norm2(&d);
1905 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1906
1907 return d;
1908 }
1909
1910
1911
1912 inline Vector::real_type
1913 Vector::lp_norm(const TrilinosScalar p) const
1914 {
1916
1917 TrilinosScalar norm = 0;
1918 TrilinosScalar sum = 0;
1919 const size_type n_local = locally_owned_size();
1920
1921 // loop over all the elements because
1922 // Trilinos does not support lp norms
1923 for (size_type i = 0; i < n_local; ++i)
1924 sum += std::pow(std::fabs((*vector)[0][i]), p);
1925
1926 norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1927
1928 return norm;
1929 }
1930
1931
1932
1933 inline Vector::real_type
1934 Vector::linfty_norm() const
1935 {
1936 // while we disallow the other
1937 // norm operations on ghosted
1938 // vectors, this particular norm
1939 // is safe to run even in the
1940 // presence of ghost elements
1942 const int ierr = vector->NormInf(&d);
1943 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1944
1945 return d;
1946 }
1947
1948
1949
1950 inline TrilinosScalar
1952 const Vector & V,
1953 const Vector & W)
1954 {
1955 this->add(a, V);
1956 return *this * W;
1957 }
1958
1959
1960
1961 // inline also scalar products, vector
1962 // additions etc. since they are all
1963 // representable by a single Trilinos
1964 // call. This reduces the overhead of the
1965 // wrapper class.
1966 inline Vector &
1968 {
1969 AssertIsFinite(a);
1970
1971 const int ierr = vector->Scale(a);
1972 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1973
1974 return *this;
1975 }
1976
1977
1978
1979 inline Vector &
1981 {
1982 AssertIsFinite(a);
1983
1984 const TrilinosScalar factor = 1. / a;
1985
1986 AssertIsFinite(factor);
1987
1988 const int ierr = vector->Scale(factor);
1989 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1990
1991 return *this;
1992 }
1993
1994
1995
1996 inline Vector &
1997 Vector::operator+=(const Vector &v)
1998 {
1999 AssertDimension(size(), v.size());
2000 Assert(vector->Map().SameAs(v.vector->Map()),
2002
2003 const int ierr = vector->Update(1.0, *(v.vector), 1.0);
2004 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2005
2006 return *this;
2007 }
2008
2009
2010
2011 inline Vector &
2012 Vector::operator-=(const Vector &v)
2013 {
2014 AssertDimension(size(), v.size());
2015 Assert(vector->Map().SameAs(v.vector->Map()),
2017
2018 const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
2019 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2020
2021 return *this;
2022 }
2023
2024
2025
2026 inline void
2028 {
2029 // if we have ghost values, do not allow
2030 // writing to this vector at all.
2032 AssertIsFinite(s);
2033
2034 size_type n_local = locally_owned_size();
2035 for (size_type i = 0; i < n_local; ++i)
2036 (*vector)[0][i] += s;
2037 }
2038
2039
2040
2041 inline void
2042 Vector::add(const TrilinosScalar a, const Vector &v)
2043 {
2044 // if we have ghost values, do not allow
2045 // writing to this vector at all.
2048
2049 AssertIsFinite(a);
2050
2051 const int ierr = vector->Update(a, *(v.vector), 1.);
2052 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2053 }
2054
2055
2056
2057 inline void
2059 const Vector & v,
2060 const TrilinosScalar b,
2061 const Vector & w)
2062 {
2063 // if we have ghost values, do not allow
2064 // writing to this vector at all.
2067 AssertDimension(locally_owned_size(), w.locally_owned_size());
2068
2069 AssertIsFinite(a);
2070 AssertIsFinite(b);
2071
2072 const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2073
2074 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2075 }
2076
2077
2078
2079 inline void
2080 Vector::sadd(const TrilinosScalar s, const Vector &v)
2081 {
2082 // if we have ghost values, do not allow
2083 // writing to this vector at all.
2085 AssertDimension(size(), v.size());
2086
2087 AssertIsFinite(s);
2088
2089 // We assume that the vectors have the same Map
2090 // if the local size is the same and if the vectors are not ghosted
2092 !v.has_ghost_elements())
2093 {
2094 Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2096 const int ierr = vector->Update(1., *(v.vector), s);
2097 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2098 }
2099 else
2100 {
2101 (*this) *= s;
2102 this->add(v, true);
2103 }
2104 }
2105
2106
2107
2108 inline void
2110 const TrilinosScalar a,
2111 const Vector & v)
2112 {
2113 // if we have ghost values, do not allow
2114 // writing to this vector at all.
2116 AssertDimension(size(), v.size());
2117 AssertIsFinite(s);
2118 AssertIsFinite(a);
2119
2120 // We assume that the vectors have the same Map
2121 // if the local size is the same and if the vectors are not ghosted
2123 !v.has_ghost_elements())
2124 {
2125 Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2127 const int ierr = vector->Update(a, *(v.vector), s);
2128 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2129 }
2130 else
2131 {
2132 (*this) *= s;
2133 Vector tmp = v;
2134 tmp *= a;
2135 this->add(tmp, true);
2136 }
2137 }
2138
2139
2140
2141 inline void
2142 Vector::scale(const Vector &factors)
2143 {
2144 // if we have ghost values, do not allow
2145 // writing to this vector at all.
2148
2149 const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2150 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2151 }
2152
2153
2154
2155 inline void
2156 Vector::equ(const TrilinosScalar a, const Vector &v)
2157 {
2158 // if we have ghost values, do not allow
2159 // writing to this vector at all.
2161 AssertIsFinite(a);
2162
2163 // If we don't have the same map, copy.
2164 if (vector->Map().SameAs(v.vector->Map()) == false)
2165 {
2166 this->sadd(0., a, v);
2167 }
2168 else
2169 {
2170 // Otherwise, just update
2171 int ierr = vector->Update(a, *v.vector, 0.0);
2172 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2173
2174 last_action = Zero;
2175 }
2176 }
2177
2178
2179
2180 inline const Epetra_MultiVector &
2182 {
2183 return static_cast<const Epetra_MultiVector &>(*vector);
2184 }
2185
2186
2187
2188 inline Epetra_FEVector &
2190 {
2191 return *vector;
2192 }
2193
2194
2195
2196 inline const Epetra_BlockMap &
2198 {
2199 return vector->Map();
2200 }
2201
2202
2203
2204 inline MPI_Comm
2206 {
2207 const Epetra_MpiComm *mpi_comm =
2208 dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2209 return mpi_comm->Comm();
2210 }
2211
2212
2213
2214 template <typename number>
2215 Vector::Vector(const IndexSet & parallel_partitioner,
2216 const ::Vector<number> &v,
2217 const MPI_Comm communicator)
2218 {
2219 *this =
2220 Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2221 owned_elements = parallel_partitioner;
2222 }
2223
2224
2225
2226 inline Vector &
2228 {
2229 AssertIsFinite(s);
2230
2231 int ierr = vector->PutScalar(s);
2232 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2233
2234 if (nonlocal_vector.get() != nullptr)
2235 {
2236 ierr = nonlocal_vector->PutScalar(0.);
2237 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2238 }
2239
2240 return *this;
2241 }
2242 } /* end of namespace MPI */
2243
2244# endif /* DOXYGEN */
2245
2246} /* end of namespace TrilinosWrappers */
2247
2251namespace internal
2252{
2253 namespace LinearOperatorImplementation
2254 {
2255 template <typename>
2256 class ReinitHelper;
2257
2262 template <>
2263 class ReinitHelper<TrilinosWrappers::MPI::Vector>
2264 {
2265 public:
2266 template <typename Matrix>
2267 static void
2268 reinit_range_vector(const Matrix & matrix,
2270 bool omit_zeroing_entries)
2271 {
2272 v.reinit(matrix.locally_owned_range_indices(),
2273 matrix.get_mpi_communicator(),
2274 omit_zeroing_entries);
2275 }
2276
2277 template <typename Matrix>
2278 static void
2279 reinit_domain_vector(const Matrix & matrix,
2281 bool omit_zeroing_entries)
2282 {
2283 v.reinit(matrix.locally_owned_domain_indices(),
2284 matrix.get_mpi_communicator(),
2285 omit_zeroing_entries);
2286 }
2287 };
2288
2289 } // namespace LinearOperatorImplementation
2290} /* namespace internal */
2291
2292
2293
2297template <>
2298struct is_serial_vector<TrilinosWrappers::MPI::Vector> : std::false_type
2299{};
2300
2301
2303
2304#endif
2305
2306#endif
size_type size() const
Definition index_set.h:1661
size_type n_elements() const
Definition index_set.h:1816
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition index_set.cc:789
Vector & operator/=(const TrilinosScalar factor)
void compress(VectorOperation::values operation)
Vector(const IndexSet &parallel_partitioning, const ::Vector< Number > &v, const MPI_Comm communicator=MPI_COMM_WORLD)
void sadd(const TrilinosScalar s, const Vector &V)
TrilinosScalar mean_value() const
void add(const size_type n_elements, const size_type *indices, const TrilinosScalar *values)
void add(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
std::unique_ptr< Epetra_MultiVector > nonlocal_vector
void add(const TrilinosScalar s)
void import_elements(const LinearAlgebra::ReadWriteVector< double > &rwv, const VectorOperation::values operation)
MPI_Comm get_mpi_communicator() const
const Epetra_BlockMap & trilinos_partitioner() const
reference operator()(const size_type index)
void import_nonlocal_data_for_fe(const ::TrilinosWrappers::SparseMatrix &matrix, const Vector &vector)
std::unique_ptr< Epetra_FEVector > vector
bool in_local_range(const size_type index) const
size_type local_size() const
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
friend class internal::VectorReference
real_type lp_norm(const TrilinosScalar p) const
IndexSet locally_owned_elements() const
Vector & operator-=(const Vector &V)
Vector & operator+=(const Vector &V)
const Epetra_MultiVector & trilinos_vector() const
const internal::VectorReference const_reference
void reinit(const Vector &v, const bool omit_zeroing_entries=false, const bool allow_different_maps=false)
std::pair< size_type, size_type > local_range() const
real_type norm_sqr() const
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< TrilinosScalar > &values) const
bool operator!=(const Vector &v) const
~Vector() override=default
Epetra_FEVector & trilinos_vector()
const_iterator begin() const
real_type linfty_norm() const
internal::VectorReference reference
TrilinosScalar min() const
void set(const std::vector< size_type > &indices, const ::Vector< TrilinosScalar > &values)
void set(const size_type n_elements, const size_type *indices, const TrilinosScalar *values)
void scale(const Vector &scaling_factors)
void equ(const TrilinosScalar a, const Vector &V)
bool operator==(const Vector &v) const
void sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V)
reference operator[](const size_type index)
const value_type * const_iterator
size_type locally_owned_size() const
TrilinosScalar operator[](const size_type index) const
TrilinosScalar add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W)
Vector & operator=(const TrilinosScalar s)
Vector & operator*=(const TrilinosScalar factor)
::types::global_dof_index size_type
TrilinosScalar operator*(const Vector &vec) const
void add(const std::vector< size_type > &indices, const ::Vector< TrilinosScalar > &values)
void extract_subvector_to(ForwardIterator indices_begin, const ForwardIterator indices_end, OutputIterator values_begin) const
std::size_t memory_consumption() const
void add(const TrilinosScalar a, const Vector &V, const TrilinosScalar b, const Vector &W)
void add(const TrilinosScalar a, const Vector &V)
void set(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
void swap(Vector &u, Vector &v)
TrilinosScalar max() const
Vector & operator=(const ::Vector< Number > &v)
const_iterator end() const
typename numbers::NumberTraits< Number >::real_type real_type
Definition vector.h:146
bool has_ghost_elements() const
const value_type * const_iterator
Definition vector.h:132
size_type size() const
value_type * iterator
Definition vector.h:131
size_type locally_owned_size() const
static void reinit_domain_vector(const Matrix &matrix, TrilinosWrappers::MPI::Vector &v, bool omit_zeroing_entries)
static void reinit_range_vector(const Matrix &matrix, TrilinosWrappers::MPI::Vector &v, bool omit_zeroing_entries)
#define DEAL_II_DEPRECATED
Definition config.h:172
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define DeclException0(Exception0)
Definition exceptions.h:465
static ::ExceptionBase & ExcGhostsPresent()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:579
static ::ExceptionBase & ExcDifferentParallelPartitioning()
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:510
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
types::global_dof_index size_type
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:472
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
int gid(const Epetra_BlockMap &map, int i)
T sum(const T &t, const MPI_Comm mpi_communicator)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
unsigned int global_dof_index
Definition types.h:82
double TrilinosScalar
Definition types.h:175