deal.II version GIT relicensing-2017-g7ae82fad8b 2024-10-22 19:20:00+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
25
28# include <deal.II/lac/vector.h>
31
32# include <Epetra_ConfigDefs.h>
33# include <Epetra_FEVector.h>
34# include <Epetra_LocalMap.h>
35# include <Epetra_Map.h>
36# include <Epetra_MpiComm.h>
37
38# include <memory>
39# include <utility>
40# include <vector>
41
43
44// Forward declarations
45# ifndef DOXYGEN
46namespace LinearAlgebra
47{
48 // Forward declaration
49 template <typename Number>
50 class ReadWriteVector;
51} // namespace LinearAlgebra
52# endif
53
64namespace TrilinosWrappers
65{
66 class SparseMatrix;
67
73 {
74 public:
76 };
77
88 namespace internal
89 {
99 class VectorReference
100 {
101 private:
102 using size_type = VectorTraits::size_type;
103
108 VectorReference(MPI::Vector &vector, const size_type index);
109
110 public:
114 VectorReference(const VectorReference &) = default;
115
127 const VectorReference &
128 operator=(const VectorReference &r) const;
129
133 VectorReference &
134 operator=(const VectorReference &r);
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
163 const VectorReference &
164 operator/=(const TrilinosScalar &s) const;
165
170 operator TrilinosScalar() const;
171
176 int,
177 << "An error with error number " << arg1
178 << " occurred while calling a Trilinos function");
179
180 private:
184 MPI::Vector &vector;
185
189 const size_type index;
190
191 // Make the vector class a friend, so that it can create objects of the
192 // present type.
193 friend class ::TrilinosWrappers::MPI::Vector;
194 };
195 } // namespace internal
200# ifndef DEAL_II_WITH_64BIT_INDICES
201 // define a helper function that queries the global ID of local ID of
202 // an Epetra_BlockMap object by calling either the 32- or 64-bit
203 // function necessary.
204 inline int
205 gid(const Epetra_BlockMap &map, int i)
206 {
207 return map.GID(i);
208 }
209# else
210 // define a helper function that queries the global ID of local ID of
211 // an Epetra_BlockMap object by calling either the 32- or 64-bit
212 // function necessary.
213 inline long long int
214 gid(const Epetra_BlockMap &map, int i)
215 {
216 return map.GID64(i);
217 }
218# endif
219
225 namespace MPI
226 {
227 class BlockVector;
228
403 class Vector : public ReadVector<TrilinosScalar>
404 {
405 public:
415 using const_iterator = const value_type *;
418
428 Vector();
429
433 Vector(const Vector &v);
434
453 explicit Vector(const IndexSet &parallel_partitioning,
454 const MPI_Comm communicator = MPI_COMM_WORLD);
455
467 Vector(const IndexSet &local,
468 const IndexSet &ghost,
469 const MPI_Comm communicator = MPI_COMM_WORLD);
470
485 Vector(const IndexSet &parallel_partitioning,
486 const Vector &v,
487 const MPI_Comm communicator = MPI_COMM_WORLD);
488
501 template <typename Number>
502 Vector(const IndexSet &parallel_partitioning,
503 const ::Vector<Number> &v,
504 const MPI_Comm communicator = MPI_COMM_WORLD);
505
514 Vector(Vector &&v); // NOLINT
515
519 ~Vector() override = default;
520
525 void
526 clear();
527
551 void
552 reinit(const Vector &v,
553 const bool omit_zeroing_entries = false,
554 const bool allow_different_maps = false);
555
578 void
579 reinit(const IndexSet &parallel_partitioning,
580 const MPI_Comm communicator = MPI_COMM_WORLD,
581 const bool omit_zeroing_entries = false);
582
617 void
618 reinit(const IndexSet &locally_owned_entries,
619 const IndexSet &locally_relevant_or_ghost_entries,
620 const MPI_Comm communicator = MPI_COMM_WORLD,
621 const bool vector_writable = false);
622
633 void
634 reinit(
635 const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner,
636 const bool make_ghosted = true,
637 const bool vector_writable = false);
638
642 void
643 reinit(const BlockVector &v, const bool import_data = false);
644
661 void
663
676 Vector &
678
708 Vector &
709 operator=(const Vector &v);
710
715 Vector &
716 operator=(Vector &&v) noexcept;
717
725 template <typename Number>
726 Vector &
727 operator=(const ::Vector<Number> &v);
728
746 void
748 const ::TrilinosWrappers::SparseMatrix &matrix,
749 const Vector &vector);
750
757 void
759 const VectorOperation::values operation);
760
765 void
767 const VectorOperation::values operation)
768 {
769 import_elements(rwv, operation);
770 }
771
772
778 bool
779 operator==(const Vector &v) const;
780
786 bool
787 operator!=(const Vector &v) const;
788
793 size() const override;
794
801
824 std::pair<size_type, size_type>
825 local_range() const;
826
834 bool
835 in_local_range(const size_type index) const;
836
852
860 bool
862
869 void
871
877 operator*(const Vector &vec) const;
878
883 norm_sqr() const;
884
889 mean_value() const;
890
895 min() const;
896
901 max() const;
902
907 l1_norm() const;
908
914 l2_norm() const;
915
921 lp_norm(const TrilinosScalar p) const;
922
927 linfty_norm() const;
928
949 add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
950
956 bool
957 all_zero() const;
958
964 bool
965 is_non_negative() const;
982 operator()(const size_type index);
983
992 operator()(const size_type index) const;
993
1000 operator[](const size_type index);
1001
1008 operator[](const size_type index) const;
1009
1025 void
1026 extract_subvector_to(const std::vector<size_type> &indices,
1027 std::vector<TrilinosScalar> &values) const;
1028
1032 virtual void
1034 const ArrayView<const size_type> &indices,
1035 const 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) noexcept;
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
1422 swap(Vector &u, Vector &v) noexcept
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 const ArrayView<const size_type> &indices,
1586 const ArrayView<TrilinosScalar> &elements) const
1587 {
1588 AssertDimension(indices.size(), elements.size());
1589 for (unsigned int i = 0; i < indices.size(); ++i)
1590 {
1591 AssertIndexRange(indices[i], size());
1592 elements[i] = (*this)[indices[i]];
1593 }
1594 }
1595
1596
1597
1598 template <typename ForwardIterator, typename OutputIterator>
1599 inline void
1600 Vector::extract_subvector_to(ForwardIterator indices_begin,
1601 const ForwardIterator indices_end,
1602 OutputIterator values_begin) const
1603 {
1604 while (indices_begin != indices_end)
1605 {
1606 *values_begin = operator()(*indices_begin);
1607 ++indices_begin;
1608 ++values_begin;
1609 }
1610 }
1611
1612
1613
1614 inline Vector::iterator
1616 {
1617 return (*vector)[0];
1618 }
1619
1620
1621
1622 inline Vector::iterator
1623 Vector::end()
1624 {
1625 return (*vector)[0] + locally_owned_size();
1626 }
1627
1628
1629
1631 Vector::begin() const
1632 {
1633 return (*vector)[0];
1634 }
1635
1636
1637
1639 Vector::end() const
1640 {
1641 return (*vector)[0] + locally_owned_size();
1642 }
1643
1644
1645
1646 inline void
1647 Vector::set(const std::vector<size_type> &indices,
1648 const std::vector<TrilinosScalar> &values)
1649 {
1650 // if we have ghost values, do not allow
1651 // writing to this vector at all.
1653
1654 AssertDimension(indices.size(), values.size());
1655
1656 set(indices.size(), indices.data(), values.data());
1657 }
1658
1659
1660
1661 inline void
1662 Vector::set(const std::vector<size_type> &indices,
1663 const ::Vector<TrilinosScalar> &values)
1664 {
1665 // if we have ghost values, do not allow
1666 // writing to this vector at all.
1668
1669 AssertDimension(indices.size(), values.size());
1670
1671 set(indices.size(), indices.data(), values.begin());
1672 }
1673
1674
1675
1676 inline void
1677 Vector::set(const size_type n_elements,
1678 const size_type *indices,
1679 const TrilinosScalar *values)
1680 {
1681 // if we have ghost values, do not allow
1682 // writing to this vector at all.
1684
1685 if (last_action == Add)
1686 {
1687 const int ierr = vector->GlobalAssemble(Add);
1688 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1689 }
1690
1691 if (last_action != Insert)
1692 last_action = Insert;
1693
1694 for (size_type i = 0; i < n_elements; ++i)
1695 {
1696 const TrilinosWrappers::types::int_type row = indices[i];
1697 const TrilinosWrappers::types::int_type local_row =
1698 vector->Map().LID(row);
1699 if (local_row != -1)
1700 (*vector)[0][local_row] = values[i];
1701 else
1702 {
1703 const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1704 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1705 compressed = false;
1706 }
1707 // in set operation, do not use the pre-allocated vector for nonlocal
1708 // entries even if it exists. This is to ensure that we really only
1709 // set the elements touched by the set() method and not all contained
1710 // in the nonlocal entries vector (there is no way to distinguish them
1711 // on the receiving processor)
1712 }
1713 }
1714
1715
1716
1717 inline void
1718 Vector::add(const std::vector<size_type> &indices,
1719 const std::vector<TrilinosScalar> &values)
1720 {
1721 // if we have ghost values, do not allow
1722 // writing to this vector at all.
1724 AssertDimension(indices.size(), values.size());
1725
1726 add(indices.size(), indices.data(), values.data());
1727 }
1728
1729
1730
1731 inline void
1732 Vector::add(const std::vector<size_type> &indices,
1733 const ::Vector<TrilinosScalar> &values)
1734 {
1735 // if we have ghost values, do not allow
1736 // writing to this vector at all.
1738 AssertDimension(indices.size(), values.size());
1739
1740 add(indices.size(), indices.data(), values.begin());
1741 }
1742
1743
1744
1745 inline void
1746 Vector::add(const size_type n_elements,
1747 const size_type *indices,
1748 const TrilinosScalar *values)
1749 {
1750 // if we have ghost values, do not allow
1751 // writing to this vector at all.
1753
1754 if (last_action != Add)
1755 {
1756 if (last_action == Insert)
1757 {
1758 const int ierr = vector->GlobalAssemble(Insert);
1759 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1760 }
1761 last_action = Add;
1762 }
1763
1764 for (size_type i = 0; i < n_elements; ++i)
1765 {
1766 const size_type row = indices[i];
1767 const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1768 static_cast<TrilinosWrappers::types::int_type>(row));
1769 if (local_row != -1)
1770 (*vector)[0][local_row] += values[i];
1771 else if (nonlocal_vector.get() == nullptr)
1772 {
1773 const int ierr = vector->SumIntoGlobalValues(
1774 1,
1775 reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1776 &row),
1777 &values[i]);
1778 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1779 compressed = false;
1780 }
1781 else
1782 {
1783 // use pre-allocated vector for non-local entries if it exists for
1784 // addition operation
1786 nonlocal_vector->Map().LID(
1787 static_cast<TrilinosWrappers::types::int_type>(row));
1788 Assert(my_row != -1,
1789 ExcMessage(
1790 "Attempted to write into off-processor vector entry "
1791 "that has not be specified as being writable upon "
1792 "initialization"));
1793 (*nonlocal_vector)[0][my_row] += values[i];
1794 compressed = false;
1795 }
1796 }
1797 }
1798
1799
1800
1801 inline Vector::size_type
1802 Vector::size() const
1803 {
1804# ifndef DEAL_II_WITH_64BIT_INDICES
1805 return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1806# else
1807 return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1808# endif
1809 }
1810
1811
1812
1813 inline Vector::size_type
1815 {
1816 return owned_elements.n_elements();
1817 }
1818
1819
1820
1821 inline std::pair<Vector::size_type, Vector::size_type>
1822 Vector::local_range() const
1823 {
1824# ifndef DEAL_II_WITH_64BIT_INDICES
1825 const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1827 vector->Map().MaxMyGID() + 1;
1828# else
1830 vector->Map().MinMyGID64();
1832 vector->Map().MaxMyGID64() + 1;
1833# endif
1834
1835 Assert(
1836 end - begin == vector->Map().NumMyElements(),
1837 ExcMessage(
1838 "This function only makes sense if the elements that this "
1839 "vector stores on the current processor form a contiguous range. "
1840 "This does not appear to be the case for the current vector."));
1841
1842 return std::make_pair(begin, end);
1843 }
1844
1845
1846
1847 inline TrilinosScalar
1848 Vector::operator*(const Vector &vec) const
1849 {
1850 Assert(vector->Map().SameAs(vec.vector->Map()),
1853
1854 TrilinosScalar result;
1855
1856 const int ierr = vector->Dot(*(vec.vector), &result);
1857 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1858
1859 return result;
1860 }
1861
1862
1863
1864 inline Vector::real_type
1865 Vector::norm_sqr() const
1866 {
1867 const TrilinosScalar d = l2_norm();
1868 return d * d;
1869 }
1870
1871
1872
1873 inline TrilinosScalar
1874 Vector::mean_value() const
1875 {
1877
1879 const int ierr = vector->MeanValue(&mean);
1880 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1881
1882 return mean;
1883 }
1884
1885
1886
1887 inline TrilinosScalar
1888 Vector::min() const
1889 {
1890 TrilinosScalar min_value;
1891 const int ierr = vector->MinValue(&min_value);
1892 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1893
1894 return min_value;
1895 }
1896
1897
1898
1899 inline TrilinosScalar
1900 Vector::max() const
1901 {
1902 TrilinosScalar max_value;
1903 const int ierr = vector->MaxValue(&max_value);
1904 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1905
1906 return max_value;
1907 }
1908
1909
1910
1911 inline Vector::real_type
1912 Vector::l1_norm() const
1913 {
1915
1917 const int ierr = vector->Norm1(&d);
1918 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1919
1920 return d;
1921 }
1922
1923
1924
1925 inline Vector::real_type
1926 Vector::l2_norm() const
1927 {
1929
1931 const int ierr = vector->Norm2(&d);
1932 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1933
1934 return d;
1935 }
1936
1937
1938
1939 inline Vector::real_type
1940 Vector::lp_norm(const TrilinosScalar p) const
1941 {
1943
1944 TrilinosScalar norm = 0;
1945 TrilinosScalar sum = 0;
1946 const size_type n_local = locally_owned_size();
1947
1948 // loop over all the elements because
1949 // Trilinos does not support lp norms
1950 for (size_type i = 0; i < n_local; ++i)
1951 sum += std::pow(std::fabs((*vector)[0][i]), p);
1952
1953 norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1954
1955 return norm;
1956 }
1957
1958
1959
1960 inline Vector::real_type
1961 Vector::linfty_norm() const
1962 {
1963 // while we disallow the other
1964 // norm operations on ghosted
1965 // vectors, this particular norm
1966 // is safe to run even in the
1967 // presence of ghost elements
1969 const int ierr = vector->NormInf(&d);
1970 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1971
1972 return d;
1973 }
1974
1975
1976
1977 inline TrilinosScalar
1979 const Vector &V,
1980 const Vector &W)
1981 {
1982 this->add(a, V);
1983 return *this * W;
1984 }
1985
1986
1987
1988 // inline also scalar products, vector
1989 // additions etc. since they are all
1990 // representable by a single Trilinos
1991 // call. This reduces the overhead of the
1992 // wrapper class.
1993 inline Vector &
1995 {
1996 AssertIsFinite(a);
1997
1998 const int ierr = vector->Scale(a);
1999 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2000
2001 return *this;
2002 }
2003
2004
2005
2006 inline Vector &
2008 {
2009 AssertIsFinite(a);
2010
2011 const TrilinosScalar factor = 1. / a;
2012
2013 AssertIsFinite(factor);
2014
2015 const int ierr = vector->Scale(factor);
2016 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2017
2018 return *this;
2019 }
2020
2021
2022
2023 inline Vector &
2024 Vector::operator+=(const Vector &v)
2025 {
2026 AssertDimension(size(), v.size());
2027 Assert(vector->Map().SameAs(v.vector->Map()),
2029
2030 const int ierr = vector->Update(1.0, *(v.vector), 1.0);
2031 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2032
2033 return *this;
2034 }
2035
2036
2037
2038 inline Vector &
2039 Vector::operator-=(const Vector &v)
2040 {
2041 AssertDimension(size(), v.size());
2042 Assert(vector->Map().SameAs(v.vector->Map()),
2044
2045 const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
2046 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2047
2048 return *this;
2049 }
2050
2051
2052
2053 inline void
2055 {
2056 // if we have ghost values, do not allow
2057 // writing to this vector at all.
2059 AssertIsFinite(s);
2060
2061 size_type n_local = locally_owned_size();
2062 for (size_type i = 0; i < n_local; ++i)
2063 (*vector)[0][i] += s;
2064 }
2065
2066
2067
2068 inline void
2069 Vector::add(const TrilinosScalar a, const Vector &v)
2070 {
2071 // if we have ghost values, do not allow
2072 // writing to this vector at all.
2075
2076 AssertIsFinite(a);
2077
2078 const int ierr = vector->Update(a, *(v.vector), 1.);
2079 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2080 }
2081
2082
2083
2084 inline void
2086 const Vector &v,
2087 const TrilinosScalar b,
2088 const Vector &w)
2089 {
2090 // if we have ghost values, do not allow
2091 // writing to this vector at all.
2094 AssertDimension(locally_owned_size(), w.locally_owned_size());
2095
2096 AssertIsFinite(a);
2097 AssertIsFinite(b);
2098
2099 const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2100
2101 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2102 }
2103
2104
2105
2106 inline void
2107 Vector::sadd(const TrilinosScalar s, const Vector &v)
2108 {
2109 // if we have ghost values, do not allow
2110 // writing to this vector at all.
2112 AssertDimension(size(), v.size());
2113
2114 AssertIsFinite(s);
2115
2116 // We assume that the vectors have the same Map
2117 // if the local size is the same and if the vectors are not ghosted
2119 !v.has_ghost_elements())
2120 {
2121 Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2123 const int ierr = vector->Update(1., *(v.vector), s);
2124 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2125 }
2126 else
2127 {
2128 (*this) *= s;
2129 this->add(v, true);
2130 }
2131 }
2132
2133
2134
2135 inline void
2137 const TrilinosScalar a,
2138 const Vector &v)
2139 {
2140 // if we have ghost values, do not allow
2141 // writing to this vector at all.
2143 AssertDimension(size(), v.size());
2144 AssertIsFinite(s);
2145 AssertIsFinite(a);
2146
2147 // We assume that the vectors have the same Map
2148 // if the local size is the same and if the vectors are not ghosted
2150 !v.has_ghost_elements())
2151 {
2152 Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2154 const int ierr = vector->Update(a, *(v.vector), s);
2155 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2156 }
2157 else
2158 {
2159 (*this) *= s;
2160 Vector tmp = v;
2161 tmp *= a;
2162 this->add(tmp, true);
2163 }
2164 }
2165
2166
2167
2168 inline void
2169 Vector::scale(const Vector &factors)
2170 {
2171 // if we have ghost values, do not allow
2172 // writing to this vector at all.
2175
2176 const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2177 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2178 }
2179
2180
2181
2182 inline void
2183 Vector::equ(const TrilinosScalar a, const Vector &v)
2184 {
2185 // if we have ghost values, do not allow
2186 // writing to this vector at all.
2188 AssertIsFinite(a);
2189
2190 // If we don't have the same map, copy.
2191 if (vector->Map().SameAs(v.vector->Map()) == false)
2192 {
2193 this->sadd(0., a, v);
2194 }
2195 else
2196 {
2197 // Otherwise, just update
2198 int ierr = vector->Update(a, *v.vector, 0.0);
2199 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2200
2201 last_action = Zero;
2202 }
2203 }
2204
2205
2206
2207 inline const Epetra_MultiVector &
2209 {
2210 return static_cast<const Epetra_MultiVector &>(*vector);
2211 }
2212
2213
2214
2215 inline Epetra_FEVector &
2217 {
2218 return *vector;
2219 }
2220
2221
2222
2223 inline const Epetra_BlockMap &
2225 {
2226 return vector->Map();
2227 }
2228
2229
2230
2231 inline MPI_Comm
2233 {
2234 const Epetra_MpiComm *mpi_comm =
2235 dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2236 return mpi_comm->Comm();
2237 }
2238
2239
2240
2241 template <typename number>
2242 Vector::Vector(const IndexSet &parallel_partitioner,
2243 const ::Vector<number> &v,
2244 const MPI_Comm communicator)
2245 {
2246 *this =
2247 Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2248 owned_elements = parallel_partitioner;
2249 }
2250
2251
2252
2253 inline Vector &
2255 {
2256 AssertIsFinite(s);
2257
2258 int ierr = vector->PutScalar(s);
2259 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2260
2261 if (nonlocal_vector.get() != nullptr)
2262 {
2263 ierr = nonlocal_vector->PutScalar(0.);
2264 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2265 }
2266
2267 return *this;
2268 }
2269 } /* end of namespace MPI */
2270
2271# endif /* DOXYGEN */
2272
2273} /* end of namespace TrilinosWrappers */
2274
2278namespace internal
2279{
2280 namespace LinearOperatorImplementation
2281 {
2282 template <typename>
2283 class ReinitHelper;
2284
2289 template <>
2290 class ReinitHelper<TrilinosWrappers::MPI::Vector>
2291 {
2292 public:
2293 template <typename Matrix>
2294 static void
2295 reinit_range_vector(const Matrix &matrix,
2297 bool omit_zeroing_entries)
2298 {
2299 v.reinit(matrix.locally_owned_range_indices(),
2300 matrix.get_mpi_communicator(),
2301 omit_zeroing_entries);
2302 }
2303
2304 template <typename Matrix>
2305 static void
2306 reinit_domain_vector(const Matrix &matrix,
2308 bool omit_zeroing_entries)
2309 {
2310 v.reinit(matrix.locally_owned_domain_indices(),
2311 matrix.get_mpi_communicator(),
2312 omit_zeroing_entries);
2313 }
2314 };
2315
2316 } // namespace LinearOperatorImplementation
2317} /* namespace internal */
2318
2319
2320
2324template <>
2325struct is_serial_vector<TrilinosWrappers::MPI::Vector> : std::false_type
2326{};
2327
2328
2330
2331#endif
2332
2333#endif
std::size_t size() const
Definition array_view.h:684
size_type size() const
Definition index_set.h:1776
size_type n_elements() const
Definition index_set.h:1934
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
void swap(Vector &v) noexcept
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
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
virtual void extract_subvector_to(const ArrayView< const size_type > &indices, const ArrayView< TrilinosScalar > &elements) const override
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
void swap(Vector &u, Vector &v) noexcept
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)
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:155
bool has_ghost_elements() const
const value_type * const_iterator
Definition vector.h:141
virtual size_type size() const override
value_type * iterator
Definition vector.h:140
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:205
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#define DeclException0(Exception0)
Definition exceptions.h:466
static ::ExceptionBase & ExcGhostsPresent()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:580
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:511
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