Reference documentation for deal.II version GIT relicensing-878-g25c11e0287 2024-06-21 14:50:02+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2008 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_trilinos_vector_h
16#define dealii_trilinos_vector_h
17
18
19#include <deal.II/base/config.h>
20
21#ifdef DEAL_II_WITH_TRILINOS
26
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
74 {
75 public:
77 };
78
89 namespace internal
90 {
100 class VectorReference
101 {
102 private:
103 using size_type = VectorTraits::size_type;
104
109 VectorReference(MPI::Vector &vector, const size_type index);
110
111 public:
115 VectorReference(const VectorReference &) = default;
116
128 const VectorReference &
129 operator=(const VectorReference &r) const;
130
134 VectorReference &
135 operator=(const VectorReference &r);
136
140 const VectorReference &
141 operator=(const TrilinosScalar &s) const;
142
146 const VectorReference &
147 operator+=(const TrilinosScalar &s) const;
148
152 const VectorReference &
153 operator-=(const TrilinosScalar &s) const;
154
158 const VectorReference &
159 operator*=(const TrilinosScalar &s) const;
160
164 const VectorReference &
165 operator/=(const TrilinosScalar &s) const;
166
171 operator TrilinosScalar() const;
172
177 int,
178 << "An error with error number " << arg1
179 << " occurred while calling a Trilinos function");
180
181 private:
185 MPI::Vector &vector;
186
190 const size_type index;
191
192 // Make the vector class a friend, so that it can create objects of the
193 // present type.
194 friend class ::TrilinosWrappers::MPI::Vector;
195 };
196 } // namespace internal
201# ifndef DEAL_II_WITH_64BIT_INDICES
202 // define a helper function that queries the global ID of local ID of
203 // an Epetra_BlockMap object by calling either the 32- or 64-bit
204 // function necessary.
205 inline int
206 gid(const Epetra_BlockMap &map, int i)
207 {
208 return map.GID(i);
209 }
210# else
211 // define a helper function that queries the global ID of local ID of
212 // an Epetra_BlockMap object by calling either the 32- or 64-bit
213 // function necessary.
214 inline long long int
215 gid(const Epetra_BlockMap &map, int i)
216 {
217 return map.GID64(i);
218 }
219# endif
220
226 namespace MPI
227 {
228 class BlockVector;
229
404 class Vector : public Subscriptor, public ReadVector<TrilinosScalar>
405 {
406 public:
416 using const_iterator = const value_type *;
419
429 Vector();
430
434 Vector(const Vector &v);
435
454 explicit Vector(const IndexSet &parallel_partitioning,
455 const MPI_Comm communicator = MPI_COMM_WORLD);
456
468 Vector(const IndexSet &local,
469 const IndexSet &ghost,
470 const MPI_Comm communicator = MPI_COMM_WORLD);
471
486 Vector(const IndexSet &parallel_partitioning,
487 const Vector &v,
488 const MPI_Comm communicator = MPI_COMM_WORLD);
489
502 template <typename Number>
503 Vector(const IndexSet &parallel_partitioning,
504 const ::Vector<Number> &v,
505 const MPI_Comm communicator = MPI_COMM_WORLD);
506
515 Vector(Vector &&v); // NOLINT
516
520 ~Vector() override = default;
521
526 void
527 clear();
528
552 void
553 reinit(const Vector &v,
554 const bool omit_zeroing_entries = false,
555 const bool allow_different_maps = false);
556
579 void
580 reinit(const IndexSet &parallel_partitioning,
581 const MPI_Comm communicator = MPI_COMM_WORLD,
582 const bool omit_zeroing_entries = false);
583
618 void
619 reinit(const IndexSet &locally_owned_entries,
620 const IndexSet &locally_relevant_or_ghost_entries,
621 const MPI_Comm communicator = MPI_COMM_WORLD,
622 const bool vector_writable = false);
623
634 void
635 reinit(
636 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
637 const bool make_ghosted = true,
638 const bool vector_writable = false);
639
643 void
644 reinit(const BlockVector &v, const bool import_data = false);
645
662 void
664
677 Vector &
679
709 Vector &
710 operator=(const Vector &v);
711
716 Vector &
717 operator=(Vector &&v) noexcept;
718
726 template <typename Number>
727 Vector &
728 operator=(const ::Vector<Number> &v);
729
747 void
749 const ::TrilinosWrappers::SparseMatrix &matrix,
750 const Vector &vector);
751
758 void
760 const VectorOperation::values operation);
761
766 void
768 const VectorOperation::values operation)
769 {
770 import_elements(rwv, operation);
771 }
772
773
779 bool
780 operator==(const Vector &v) const;
781
787 bool
788 operator!=(const Vector &v) const;
789
794 size() const override;
795
802
825 std::pair<size_type, size_type>
826 local_range() const;
827
835 bool
836 in_local_range(const size_type index) const;
837
853
861 bool
863
870 void
872
878 operator*(const Vector &vec) const;
879
884 norm_sqr() const;
885
890 mean_value() const;
891
896 min() const;
897
902 max() const;
903
908 l1_norm() const;
909
915 l2_norm() const;
916
922 lp_norm(const TrilinosScalar p) const;
923
928 linfty_norm() const;
929
950 add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
951
957 bool
958 all_zero() const;
959
965 bool
966 is_non_negative() const;
983 operator()(const size_type index);
984
993 operator()(const size_type index) const;
994
1000 reference
1001 operator[](const size_type index);
1002
1009 operator[](const size_type index) const;
1010
1026 void
1027 extract_subvector_to(const std::vector<size_type> &indices,
1028 std::vector<TrilinosScalar> &values) const;
1029
1033 virtual void
1035 ArrayView<TrilinosScalar> &elements) const override;
1036
1064 template <typename ForwardIterator, typename OutputIterator>
1065 void
1066 extract_subvector_to(ForwardIterator indices_begin,
1067 const ForwardIterator indices_end,
1068 OutputIterator values_begin) const;
1069
1078 iterator
1080
1086 begin() const;
1087
1092 iterator
1094
1100 end() const;
1101
1116 void
1117 set(const std::vector<size_type> &indices,
1118 const std::vector<TrilinosScalar> &values);
1119
1124 void
1125 set(const std::vector<size_type> &indices,
1126 const ::Vector<TrilinosScalar> &values);
1127
1133 void
1134 set(const size_type n_elements,
1135 const size_type *indices,
1136 const TrilinosScalar *values);
1137
1142 void
1143 add(const std::vector<size_type> &indices,
1144 const std::vector<TrilinosScalar> &values);
1145
1150 void
1151 add(const std::vector<size_type> &indices,
1152 const ::Vector<TrilinosScalar> &values);
1153
1159 void
1160 add(const size_type n_elements,
1161 const size_type *indices,
1162 const TrilinosScalar *values);
1163
1167 Vector &
1169
1173 Vector &
1175
1179 Vector &
1181
1185 Vector &
1187
1192 void
1194
1207 void
1208 add(const Vector &V, const bool allow_different_maps = false);
1209
1213 void
1214 add(const TrilinosScalar a, const Vector &V);
1215
1219 void
1221 const Vector &V,
1222 const TrilinosScalar b,
1223 const Vector &W);
1224
1229 void
1230 sadd(const TrilinosScalar s, const Vector &V);
1231
1235 void
1236 sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1237
1243 void
1244 scale(const Vector &scaling_factors);
1245
1249 void
1250 equ(const TrilinosScalar a, const Vector &V);
1262 const Epetra_MultiVector &
1264
1269 Epetra_FEVector &
1271
1276 const Epetra_BlockMap &
1278
1286 void
1287 print(std::ostream &out,
1288 const unsigned int precision = 3,
1289 const bool scientific = true,
1290 const bool across = true) const;
1291
1305 void
1306 swap(Vector &v);
1307
1311 std::size_t
1312 memory_consumption() const;
1313
1317 MPI_Comm
1325
1330 int,
1331 << "An error with error number " << arg1
1332 << " occurred while calling a Trilinos function");
1333
1339 size_type,
1340 size_type,
1341 size_type,
1342 size_type,
1343 << "You are trying to access element " << arg1
1344 << " of a distributed vector, but this element is not stored "
1345 << "on the current processor. Note: There are " << arg2
1346 << " elements stored "
1347 << "on the current processor from within the range [" << arg3 << ','
1348 << arg4 << "] but Trilinos vectors need not store contiguous "
1349 << "ranges on each processor, and not every element in "
1350 << "this range may in fact be stored locally."
1351 << "\n\n"
1352 << "A common source for this kind of problem is that you "
1353 << "are passing a 'fully distributed' vector into a function "
1354 << "that needs read access to vector elements that correspond "
1355 << "to degrees of freedom on ghost cells (or at least to "
1356 << "'locally active' degrees of freedom that are not also "
1357 << "'locally owned'). You need to pass a vector that has these "
1358 << "elements as ghost entries.");
1359
1360 private:
1372 Epetra_CombineMode last_action;
1373
1379
1385
1391 std::unique_ptr<Epetra_FEVector> vector;
1392
1398 std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1399
1404
1405 // Make the reference class a friend.
1407 };
1408
1409
1410
1411 // ------------------- inline and template functions --------------
1412
1413
1421 inline void
1423 {
1424 u.swap(v);
1425 }
1426 } // namespace MPI
1427
1428# ifndef DOXYGEN
1429
1430 namespace internal
1431 {
1432 inline VectorReference::VectorReference(MPI::Vector &vector,
1433 const size_type index)
1434 : vector(vector)
1435 , index(index)
1436 {}
1437
1438
1439 inline const VectorReference &
1440 VectorReference::operator=(const VectorReference &r) const
1441 {
1442 // as explained in the class
1443 // documentation, this is not the copy
1444 // operator. so simply pass on to the
1445 // "correct" assignment operator
1446 *this = static_cast<TrilinosScalar>(r);
1447
1448 return *this;
1449 }
1450
1451
1452
1453 inline VectorReference &
1454 VectorReference::operator=(const VectorReference &r)
1455 {
1456 // as above
1457 *this = static_cast<TrilinosScalar>(r);
1458
1459 return *this;
1460 }
1461
1462
1463 inline const VectorReference &
1464 VectorReference::operator=(const TrilinosScalar &value) const
1465 {
1466 vector.set(1, &index, &value);
1467 return *this;
1468 }
1469
1470
1471
1472 inline const VectorReference &
1473 VectorReference::operator+=(const TrilinosScalar &value) const
1474 {
1475 vector.add(1, &index, &value);
1476 return *this;
1477 }
1478
1479
1480
1481 inline const VectorReference &
1482 VectorReference::operator-=(const TrilinosScalar &value) const
1483 {
1484 TrilinosScalar new_value = -value;
1485 vector.add(1, &index, &new_value);
1486 return *this;
1487 }
1488
1489
1490
1491 inline const VectorReference &
1492 VectorReference::operator*=(const TrilinosScalar &value) const
1493 {
1494 TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1495 vector.set(1, &index, &new_value);
1496 return *this;
1497 }
1498
1499
1500
1501 inline const VectorReference &
1502 VectorReference::operator/=(const TrilinosScalar &value) const
1503 {
1504 TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1505 vector.set(1, &index, &new_value);
1506 return *this;
1507 }
1508 } // namespace internal
1509
1510 namespace MPI
1511 {
1512 inline bool
1513 Vector::in_local_range(const size_type index) const
1514 {
1515 std::pair<size_type, size_type> range = local_range();
1516
1517 return ((index >= range.first) && (index < range.second));
1518 }
1519
1520
1521
1522 inline IndexSet
1524 {
1526 ExcMessage(
1527 "The locally owned elements have not been properly initialized!"
1528 " This happens for example if this object has been initialized"
1529 " with exactly one overlapping IndexSet."));
1530 return owned_elements;
1531 }
1532
1533
1534
1535 inline bool
1537 {
1538 return has_ghosts;
1539 }
1540
1541
1542
1543 inline void
1545 {}
1546
1547
1548
1549 inline internal::VectorReference
1550 Vector::operator()(const size_type index)
1551 {
1552 return internal::VectorReference(*this, index);
1553 }
1554
1555
1556
1557 inline internal::VectorReference
1558 Vector::operator[](const size_type index)
1559 {
1560 return operator()(index);
1561 }
1562
1563
1564
1565 inline TrilinosScalar
1566 Vector::operator[](const size_type index) const
1567 {
1568 return operator()(index);
1569 }
1570
1571
1572
1573 inline void
1574 Vector::extract_subvector_to(const std::vector<size_type> &indices,
1575 std::vector<TrilinosScalar> &values) const
1576 {
1577 for (size_type i = 0; i < indices.size(); ++i)
1578 values[i] = operator()(indices[i]);
1579 }
1580
1581
1582
1583 inline void
1585 ArrayView<TrilinosScalar> &elements) const
1586 {
1587 AssertDimension(indices.size(), elements.size());
1588 for (unsigned int i = 0; i < indices.size(); ++i)
1589 {
1590 AssertIndexRange(indices[i], size());
1591 elements[i] = (*this)[indices[i]];
1592 }
1593 }
1594
1595
1596
1597 template <typename ForwardIterator, typename OutputIterator>
1598 inline void
1599 Vector::extract_subvector_to(ForwardIterator indices_begin,
1600 const ForwardIterator indices_end,
1601 OutputIterator values_begin) const
1602 {
1603 while (indices_begin != indices_end)
1604 {
1605 *values_begin = operator()(*indices_begin);
1606 ++indices_begin;
1607 ++values_begin;
1608 }
1609 }
1610
1611
1612
1613 inline Vector::iterator
1615 {
1616 return (*vector)[0];
1617 }
1618
1619
1620
1621 inline Vector::iterator
1622 Vector::end()
1623 {
1624 return (*vector)[0] + locally_owned_size();
1625 }
1626
1627
1628
1630 Vector::begin() const
1631 {
1632 return (*vector)[0];
1633 }
1634
1635
1636
1638 Vector::end() const
1639 {
1640 return (*vector)[0] + locally_owned_size();
1641 }
1642
1643
1644
1645 inline void
1646 Vector::set(const std::vector<size_type> &indices,
1647 const std::vector<TrilinosScalar> &values)
1648 {
1649 // if we have ghost values, do not allow
1650 // writing to this vector at all.
1652
1653 AssertDimension(indices.size(), values.size());
1654
1655 set(indices.size(), indices.data(), values.data());
1656 }
1657
1658
1659
1660 inline void
1661 Vector::set(const std::vector<size_type> &indices,
1662 const ::Vector<TrilinosScalar> &values)
1663 {
1664 // if we have ghost values, do not allow
1665 // writing to this vector at all.
1667
1668 AssertDimension(indices.size(), values.size());
1669
1670 set(indices.size(), indices.data(), values.begin());
1671 }
1672
1673
1674
1675 inline void
1676 Vector::set(const size_type n_elements,
1677 const size_type *indices,
1678 const TrilinosScalar *values)
1679 {
1680 // if we have ghost values, do not allow
1681 // writing to this vector at all.
1683
1684 if (last_action == Add)
1685 {
1686 const int ierr = vector->GlobalAssemble(Add);
1687 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1688 }
1689
1690 if (last_action != Insert)
1691 last_action = Insert;
1692
1693 for (size_type i = 0; i < n_elements; ++i)
1694 {
1695 const TrilinosWrappers::types::int_type row = indices[i];
1696 const TrilinosWrappers::types::int_type local_row =
1697 vector->Map().LID(row);
1698 if (local_row != -1)
1699 (*vector)[0][local_row] = values[i];
1700 else
1701 {
1702 const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1703 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1704 compressed = false;
1705 }
1706 // in set operation, do not use the pre-allocated vector for nonlocal
1707 // entries even if it exists. This is to ensure that we really only
1708 // set the elements touched by the set() method and not all contained
1709 // in the nonlocal entries vector (there is no way to distinguish them
1710 // on the receiving processor)
1711 }
1712 }
1713
1714
1715
1716 inline void
1717 Vector::add(const std::vector<size_type> &indices,
1718 const std::vector<TrilinosScalar> &values)
1719 {
1720 // if we have ghost values, do not allow
1721 // writing to this vector at all.
1723 AssertDimension(indices.size(), values.size());
1724
1725 add(indices.size(), indices.data(), values.data());
1726 }
1727
1728
1729
1730 inline void
1731 Vector::add(const std::vector<size_type> &indices,
1732 const ::Vector<TrilinosScalar> &values)
1733 {
1734 // if we have ghost values, do not allow
1735 // writing to this vector at all.
1737 AssertDimension(indices.size(), values.size());
1738
1739 add(indices.size(), indices.data(), values.begin());
1740 }
1741
1742
1743
1744 inline void
1745 Vector::add(const size_type n_elements,
1746 const size_type *indices,
1747 const TrilinosScalar *values)
1748 {
1749 // if we have ghost values, do not allow
1750 // writing to this vector at all.
1752
1753 if (last_action != Add)
1754 {
1755 if (last_action == Insert)
1756 {
1757 const int ierr = vector->GlobalAssemble(Insert);
1758 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1759 }
1760 last_action = Add;
1761 }
1762
1763 for (size_type i = 0; i < n_elements; ++i)
1764 {
1765 const size_type row = indices[i];
1766 const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1767 static_cast<TrilinosWrappers::types::int_type>(row));
1768 if (local_row != -1)
1769 (*vector)[0][local_row] += values[i];
1770 else if (nonlocal_vector.get() == nullptr)
1771 {
1772 const int ierr = vector->SumIntoGlobalValues(
1773 1,
1774 reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1775 &row),
1776 &values[i]);
1777 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1778 compressed = false;
1779 }
1780 else
1781 {
1782 // use pre-allocated vector for non-local entries if it exists for
1783 // addition operation
1785 nonlocal_vector->Map().LID(
1786 static_cast<TrilinosWrappers::types::int_type>(row));
1787 Assert(my_row != -1,
1788 ExcMessage(
1789 "Attempted to write into off-processor vector entry "
1790 "that has not be specified as being writable upon "
1791 "initialization"));
1792 (*nonlocal_vector)[0][my_row] += values[i];
1793 compressed = false;
1794 }
1795 }
1796 }
1797
1798
1799
1800 inline Vector::size_type
1801 Vector::size() const
1802 {
1803# ifndef DEAL_II_WITH_64BIT_INDICES
1804 return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1805# else
1806 return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1807# endif
1808 }
1809
1810
1811
1812 inline Vector::size_type
1814 {
1815 return owned_elements.n_elements();
1816 }
1817
1818
1819
1820 inline std::pair<Vector::size_type, Vector::size_type>
1821 Vector::local_range() const
1822 {
1823# ifndef DEAL_II_WITH_64BIT_INDICES
1824 const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1826 vector->Map().MaxMyGID() + 1;
1827# else
1829 vector->Map().MinMyGID64();
1831 vector->Map().MaxMyGID64() + 1;
1832# endif
1833
1834 Assert(
1835 end - begin == vector->Map().NumMyElements(),
1836 ExcMessage(
1837 "This function only makes sense if the elements that this "
1838 "vector stores on the current processor form a contiguous range. "
1839 "This does not appear to be the case for the current vector."));
1840
1841 return std::make_pair(begin, end);
1842 }
1843
1844
1845
1846 inline TrilinosScalar
1847 Vector::operator*(const Vector &vec) const
1848 {
1849 Assert(vector->Map().SameAs(vec.vector->Map()),
1852
1853 TrilinosScalar result;
1854
1855 const int ierr = vector->Dot(*(vec.vector), &result);
1856 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1857
1858 return result;
1859 }
1860
1861
1862
1863 inline Vector::real_type
1864 Vector::norm_sqr() const
1865 {
1866 const TrilinosScalar d = l2_norm();
1867 return d * d;
1868 }
1869
1870
1871
1872 inline TrilinosScalar
1873 Vector::mean_value() const
1874 {
1876
1878 const int ierr = vector->MeanValue(&mean);
1879 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1880
1881 return mean;
1882 }
1883
1884
1885
1886 inline TrilinosScalar
1887 Vector::min() const
1888 {
1889 TrilinosScalar min_value;
1890 const int ierr = vector->MinValue(&min_value);
1891 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1892
1893 return min_value;
1894 }
1895
1896
1897
1898 inline TrilinosScalar
1899 Vector::max() const
1900 {
1901 TrilinosScalar max_value;
1902 const int ierr = vector->MaxValue(&max_value);
1903 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1904
1905 return max_value;
1906 }
1907
1908
1909
1910 inline Vector::real_type
1911 Vector::l1_norm() const
1912 {
1914
1916 const int ierr = vector->Norm1(&d);
1917 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1918
1919 return d;
1920 }
1921
1922
1923
1924 inline Vector::real_type
1925 Vector::l2_norm() const
1926 {
1928
1930 const int ierr = vector->Norm2(&d);
1931 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1932
1933 return d;
1934 }
1935
1936
1937
1938 inline Vector::real_type
1939 Vector::lp_norm(const TrilinosScalar p) const
1940 {
1942
1943 TrilinosScalar norm = 0;
1944 TrilinosScalar sum = 0;
1945 const size_type n_local = locally_owned_size();
1946
1947 // loop over all the elements because
1948 // Trilinos does not support lp norms
1949 for (size_type i = 0; i < n_local; ++i)
1950 sum += std::pow(std::fabs((*vector)[0][i]), p);
1951
1952 norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1953
1954 return norm;
1955 }
1956
1957
1958
1959 inline Vector::real_type
1960 Vector::linfty_norm() const
1961 {
1962 // while we disallow the other
1963 // norm operations on ghosted
1964 // vectors, this particular norm
1965 // is safe to run even in the
1966 // presence of ghost elements
1968 const int ierr = vector->NormInf(&d);
1969 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1970
1971 return d;
1972 }
1973
1974
1975
1976 inline TrilinosScalar
1978 const Vector &V,
1979 const Vector &W)
1980 {
1981 this->add(a, V);
1982 return *this * W;
1983 }
1984
1985
1986
1987 // inline also scalar products, vector
1988 // additions etc. since they are all
1989 // representable by a single Trilinos
1990 // call. This reduces the overhead of the
1991 // wrapper class.
1992 inline Vector &
1994 {
1995 AssertIsFinite(a);
1996
1997 const int ierr = vector->Scale(a);
1998 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1999
2000 return *this;
2001 }
2002
2003
2004
2005 inline Vector &
2007 {
2008 AssertIsFinite(a);
2009
2010 const TrilinosScalar factor = 1. / a;
2011
2012 AssertIsFinite(factor);
2013
2014 const int ierr = vector->Scale(factor);
2015 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2016
2017 return *this;
2018 }
2019
2020
2021
2022 inline Vector &
2023 Vector::operator+=(const Vector &v)
2024 {
2025 AssertDimension(size(), v.size());
2026 Assert(vector->Map().SameAs(v.vector->Map()),
2028
2029 const int ierr = vector->Update(1.0, *(v.vector), 1.0);
2030 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2031
2032 return *this;
2033 }
2034
2035
2036
2037 inline Vector &
2038 Vector::operator-=(const Vector &v)
2039 {
2040 AssertDimension(size(), v.size());
2041 Assert(vector->Map().SameAs(v.vector->Map()),
2043
2044 const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
2045 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2046
2047 return *this;
2048 }
2049
2050
2051
2052 inline void
2054 {
2055 // if we have ghost values, do not allow
2056 // writing to this vector at all.
2058 AssertIsFinite(s);
2059
2060 size_type n_local = locally_owned_size();
2061 for (size_type i = 0; i < n_local; ++i)
2062 (*vector)[0][i] += s;
2063 }
2064
2065
2066
2067 inline void
2068 Vector::add(const TrilinosScalar a, const Vector &v)
2069 {
2070 // if we have ghost values, do not allow
2071 // writing to this vector at all.
2074
2075 AssertIsFinite(a);
2076
2077 const int ierr = vector->Update(a, *(v.vector), 1.);
2078 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2079 }
2080
2081
2082
2083 inline void
2085 const Vector &v,
2086 const TrilinosScalar b,
2087 const Vector &w)
2088 {
2089 // if we have ghost values, do not allow
2090 // writing to this vector at all.
2093 AssertDimension(locally_owned_size(), w.locally_owned_size());
2094
2095 AssertIsFinite(a);
2096 AssertIsFinite(b);
2097
2098 const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2099
2100 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2101 }
2102
2103
2104
2105 inline void
2106 Vector::sadd(const TrilinosScalar s, const Vector &v)
2107 {
2108 // if we have ghost values, do not allow
2109 // writing to this vector at all.
2111 AssertDimension(size(), v.size());
2112
2113 AssertIsFinite(s);
2114
2115 // We assume that the vectors have the same Map
2116 // if the local size is the same and if the vectors are not ghosted
2118 !v.has_ghost_elements())
2119 {
2120 Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2122 const int ierr = vector->Update(1., *(v.vector), s);
2123 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2124 }
2125 else
2126 {
2127 (*this) *= s;
2128 this->add(v, true);
2129 }
2130 }
2131
2132
2133
2134 inline void
2136 const TrilinosScalar a,
2137 const Vector &v)
2138 {
2139 // if we have ghost values, do not allow
2140 // writing to this vector at all.
2142 AssertDimension(size(), v.size());
2143 AssertIsFinite(s);
2144 AssertIsFinite(a);
2145
2146 // We assume that the vectors have the same Map
2147 // if the local size is the same and if the vectors are not ghosted
2149 !v.has_ghost_elements())
2150 {
2151 Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2153 const int ierr = vector->Update(a, *(v.vector), s);
2154 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2155 }
2156 else
2157 {
2158 (*this) *= s;
2159 Vector tmp = v;
2160 tmp *= a;
2161 this->add(tmp, true);
2162 }
2163 }
2164
2165
2166
2167 inline void
2168 Vector::scale(const Vector &factors)
2169 {
2170 // if we have ghost values, do not allow
2171 // writing to this vector at all.
2174
2175 const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2176 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2177 }
2178
2179
2180
2181 inline void
2182 Vector::equ(const TrilinosScalar a, const Vector &v)
2183 {
2184 // if we have ghost values, do not allow
2185 // writing to this vector at all.
2187 AssertIsFinite(a);
2188
2189 // If we don't have the same map, copy.
2190 if (vector->Map().SameAs(v.vector->Map()) == false)
2191 {
2192 this->sadd(0., a, v);
2193 }
2194 else
2195 {
2196 // Otherwise, just update
2197 int ierr = vector->Update(a, *v.vector, 0.0);
2198 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2199
2200 last_action = Zero;
2201 }
2202 }
2203
2204
2205
2206 inline const Epetra_MultiVector &
2208 {
2209 return static_cast<const Epetra_MultiVector &>(*vector);
2210 }
2211
2212
2213
2214 inline Epetra_FEVector &
2216 {
2217 return *vector;
2218 }
2219
2220
2221
2222 inline const Epetra_BlockMap &
2224 {
2225 return vector->Map();
2226 }
2227
2228
2229
2230 inline MPI_Comm
2232 {
2233 const Epetra_MpiComm *mpi_comm =
2234 dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2235 return mpi_comm->Comm();
2236 }
2237
2238
2239
2240 template <typename number>
2241 Vector::Vector(const IndexSet &parallel_partitioner,
2242 const ::Vector<number> &v,
2243 const MPI_Comm communicator)
2244 {
2245 *this =
2246 Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2247 owned_elements = parallel_partitioner;
2248 }
2249
2250
2251
2252 inline Vector &
2254 {
2255 AssertIsFinite(s);
2256
2257 int ierr = vector->PutScalar(s);
2258 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2259
2260 if (nonlocal_vector.get() != nullptr)
2261 {
2262 ierr = nonlocal_vector->PutScalar(0.);
2263 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2264 }
2265
2266 return *this;
2267 }
2268 } /* end of namespace MPI */
2269
2270# endif /* DOXYGEN */
2271
2272} /* end of namespace TrilinosWrappers */
2273
2277namespace internal
2278{
2279 namespace LinearOperatorImplementation
2280 {
2281 template <typename>
2282 class ReinitHelper;
2283
2288 template <>
2289 class ReinitHelper<TrilinosWrappers::MPI::Vector>
2290 {
2291 public:
2292 template <typename Matrix>
2293 static void
2294 reinit_range_vector(const Matrix &matrix,
2296 bool omit_zeroing_entries)
2297 {
2298 v.reinit(matrix.locally_owned_range_indices(),
2299 matrix.get_mpi_communicator(),
2300 omit_zeroing_entries);
2301 }
2302
2303 template <typename Matrix>
2304 static void
2305 reinit_domain_vector(const Matrix &matrix,
2307 bool omit_zeroing_entries)
2308 {
2309 v.reinit(matrix.locally_owned_domain_indices(),
2310 matrix.get_mpi_communicator(),
2311 omit_zeroing_entries);
2312 }
2313 };
2314
2315 } // namespace LinearOperatorImplementation
2316} /* namespace internal */
2317
2318
2319
2323template <>
2324struct is_serial_vector<TrilinosWrappers::MPI::Vector> : std::false_type
2325{};
2326
2327
2329
2330#endif
2331
2332#endif
std::size_t size() const
Definition array_view.h:684
size_type size() const
Definition index_set.h:1765
size_type n_elements() const
Definition index_set.h:1923
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
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)
VectorTraits::size_type size_type
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
size_type size() const override
bool in_local_range(const size_type index) 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
virtual void extract_subvector_to(const ArrayView< const size_type > &indices, ArrayView< TrilinosScalar > &elements) const override
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)
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
::types::global_dof_index size_type
typename numbers::NumberTraits< Number >::real_type real_type
Definition vector.h:157
bool has_ghost_elements() const
const value_type * const_iterator
Definition vector.h:143
virtual size_type size() const override
value_type * iterator
Definition vector.h:142
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:207
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define DeclException0(Exception0)
Definition exceptions.h:471
static ::ExceptionBase & ExcGhostsPresent()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:585
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)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
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:81
double TrilinosScalar
Definition types.h:178