Reference documentation for deal.II version GIT relicensing-416-gd521cc2ad4 2024-04-17 12:40: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
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, public ReadVector<TrilinosScalar>
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
702 Vector &
703 operator=(const Vector &v);
704
709 Vector &
710 operator=(Vector &&v) noexcept;
711
719 template <typename Number>
720 Vector &
721 operator=(const ::Vector<Number> &v);
722
740 void
742 const ::TrilinosWrappers::SparseMatrix &matrix,
743 const Vector &vector);
744
751 void
753 const VectorOperation::values operation);
754
759 void
761 const VectorOperation::values operation)
762 {
763 import_elements(rwv, operation);
764 }
765
766
772 bool
773 operator==(const Vector &v) const;
774
780 bool
781 operator!=(const Vector &v) const;
782
787 size() const override;
788
795
818 std::pair<size_type, size_type>
819 local_range() const;
820
828 bool
829 in_local_range(const size_type index) const;
830
846
854 bool
856
863 void
865
871 operator*(const Vector &vec) const;
872
877 norm_sqr() const;
878
883 mean_value() const;
884
889 min() const;
890
895 max() const;
896
901 l1_norm() const;
902
908 l2_norm() const;
909
915 lp_norm(const TrilinosScalar p) const;
916
921 linfty_norm() const;
922
943 add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
944
950 bool
951 all_zero() const;
952
958 bool
959 is_non_negative() const;
976 operator()(const size_type index);
977
986 operator()(const size_type index) const;
987
994 operator[](const size_type index);
995
1002 operator[](const size_type index) const;
1003
1019 void
1020 extract_subvector_to(const std::vector<size_type> &indices,
1021 std::vector<TrilinosScalar> &values) const;
1022
1026 virtual void
1028 ArrayView<TrilinosScalar> &elements) const override;
1029
1057 template <typename ForwardIterator, typename OutputIterator>
1058 void
1059 extract_subvector_to(ForwardIterator indices_begin,
1060 const ForwardIterator indices_end,
1061 OutputIterator values_begin) const;
1062
1071 iterator
1073
1079 begin() const;
1080
1085 iterator
1087
1093 end() const;
1094
1109 void
1110 set(const std::vector<size_type> &indices,
1111 const std::vector<TrilinosScalar> &values);
1112
1117 void
1118 set(const std::vector<size_type> &indices,
1119 const ::Vector<TrilinosScalar> &values);
1120
1126 void
1127 set(const size_type n_elements,
1128 const size_type *indices,
1129 const TrilinosScalar *values);
1130
1135 void
1136 add(const std::vector<size_type> &indices,
1137 const std::vector<TrilinosScalar> &values);
1138
1143 void
1144 add(const std::vector<size_type> &indices,
1145 const ::Vector<TrilinosScalar> &values);
1146
1152 void
1153 add(const size_type n_elements,
1154 const size_type *indices,
1155 const TrilinosScalar *values);
1156
1160 Vector &
1162
1166 Vector &
1168
1172 Vector &
1174
1178 Vector &
1180
1185 void
1187
1200 void
1201 add(const Vector &V, const bool allow_different_maps = false);
1202
1206 void
1207 add(const TrilinosScalar a, const Vector &V);
1208
1212 void
1214 const Vector &V,
1215 const TrilinosScalar b,
1216 const Vector &W);
1217
1222 void
1223 sadd(const TrilinosScalar s, const Vector &V);
1224
1228 void
1229 sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1230
1236 void
1237 scale(const Vector &scaling_factors);
1238
1242 void
1243 equ(const TrilinosScalar a, const Vector &V);
1255 const Epetra_MultiVector &
1257
1262 Epetra_FEVector &
1264
1269 const Epetra_BlockMap &
1271
1279 void
1280 print(std::ostream &out,
1281 const unsigned int precision = 3,
1282 const bool scientific = true,
1283 const bool across = true) const;
1284
1298 void
1299 swap(Vector &v);
1300
1304 std::size_t
1305 memory_consumption() const;
1306
1310 MPI_Comm
1318
1323 int,
1324 << "An error with error number " << arg1
1325 << " occurred while calling a Trilinos function");
1326
1332 size_type,
1333 size_type,
1334 size_type,
1335 size_type,
1336 << "You are trying to access element " << arg1
1337 << " of a distributed vector, but this element is not stored "
1338 << "on the current processor. Note: There are " << arg2
1339 << " elements stored "
1340 << "on the current processor from within the range [" << arg3 << ','
1341 << arg4 << "] but Trilinos vectors need not store contiguous "
1342 << "ranges on each processor, and not every element in "
1343 << "this range may in fact be stored locally."
1344 << "\n\n"
1345 << "A common source for this kind of problem is that you "
1346 << "are passing a 'fully distributed' vector into a function "
1347 << "that needs read access to vector elements that correspond "
1348 << "to degrees of freedom on ghost cells (or at least to "
1349 << "'locally active' degrees of freedom that are not also "
1350 << "'locally owned'). You need to pass a vector that has these "
1351 << "elements as ghost entries.");
1352
1353 private:
1365 Epetra_CombineMode last_action;
1366
1372
1378
1384 std::unique_ptr<Epetra_FEVector> vector;
1385
1391 std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1392
1397
1398 // Make the reference class a friend.
1400 };
1401
1402
1403
1404 // ------------------- inline and template functions --------------
1405
1406
1414 inline void
1416 {
1417 u.swap(v);
1418 }
1419 } // namespace MPI
1420
1421# ifndef DOXYGEN
1422
1423 namespace internal
1424 {
1425 inline VectorReference::VectorReference(MPI::Vector &vector,
1426 const size_type index)
1427 : vector(vector)
1428 , index(index)
1429 {}
1430
1431
1432 inline const VectorReference &
1433 VectorReference::operator=(const VectorReference &r) const
1434 {
1435 // as explained in the class
1436 // documentation, this is not the copy
1437 // operator. so simply pass on to the
1438 // "correct" assignment operator
1439 *this = static_cast<TrilinosScalar>(r);
1440
1441 return *this;
1442 }
1443
1444
1445
1446 inline VectorReference &
1447 VectorReference::operator=(const VectorReference &r)
1448 {
1449 // as above
1450 *this = static_cast<TrilinosScalar>(r);
1451
1452 return *this;
1453 }
1454
1455
1456 inline const VectorReference &
1457 VectorReference::operator=(const TrilinosScalar &value) const
1458 {
1459 vector.set(1, &index, &value);
1460 return *this;
1461 }
1462
1463
1464
1465 inline const VectorReference &
1466 VectorReference::operator+=(const TrilinosScalar &value) const
1467 {
1468 vector.add(1, &index, &value);
1469 return *this;
1470 }
1471
1472
1473
1474 inline const VectorReference &
1475 VectorReference::operator-=(const TrilinosScalar &value) const
1476 {
1477 TrilinosScalar new_value = -value;
1478 vector.add(1, &index, &new_value);
1479 return *this;
1480 }
1481
1482
1483
1484 inline const VectorReference &
1485 VectorReference::operator*=(const TrilinosScalar &value) const
1486 {
1487 TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1488 vector.set(1, &index, &new_value);
1489 return *this;
1490 }
1491
1492
1493
1494 inline const VectorReference &
1495 VectorReference::operator/=(const TrilinosScalar &value) const
1496 {
1497 TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1498 vector.set(1, &index, &new_value);
1499 return *this;
1500 }
1501 } // namespace internal
1502
1503 namespace MPI
1504 {
1505 inline bool
1506 Vector::in_local_range(const size_type index) const
1507 {
1508 std::pair<size_type, size_type> range = local_range();
1509
1510 return ((index >= range.first) && (index < range.second));
1511 }
1512
1513
1514
1515 inline IndexSet
1517 {
1519 ExcMessage(
1520 "The locally owned elements have not been properly initialized!"
1521 " This happens for example if this object has been initialized"
1522 " with exactly one overlapping IndexSet."));
1523 return owned_elements;
1524 }
1525
1526
1527
1528 inline bool
1530 {
1531 return has_ghosts;
1532 }
1533
1534
1535
1536 inline void
1538 {}
1539
1540
1541
1542 inline internal::VectorReference
1543 Vector::operator()(const size_type index)
1544 {
1545 return internal::VectorReference(*this, index);
1546 }
1547
1548
1549
1550 inline internal::VectorReference
1551 Vector::operator[](const size_type index)
1552 {
1553 return operator()(index);
1554 }
1555
1556
1557
1558 inline TrilinosScalar
1559 Vector::operator[](const size_type index) const
1560 {
1561 return operator()(index);
1562 }
1563
1564
1565
1566 inline void
1567 Vector::extract_subvector_to(const std::vector<size_type> &indices,
1568 std::vector<TrilinosScalar> &values) const
1569 {
1570 for (size_type i = 0; i < indices.size(); ++i)
1571 values[i] = operator()(indices[i]);
1572 }
1573
1574
1575
1576 inline void
1578 ArrayView<TrilinosScalar> &elements) const
1579 {
1580 AssertDimension(indices.size(), elements.size());
1581 for (unsigned int i = 0; i < indices.size(); ++i)
1582 {
1583 AssertIndexRange(indices[i], size());
1584 elements[i] = (*this)[indices[i]];
1585 }
1586 }
1587
1588
1589
1590 template <typename ForwardIterator, typename OutputIterator>
1591 inline void
1592 Vector::extract_subvector_to(ForwardIterator indices_begin,
1593 const ForwardIterator indices_end,
1594 OutputIterator values_begin) const
1595 {
1596 while (indices_begin != indices_end)
1597 {
1598 *values_begin = operator()(*indices_begin);
1599 ++indices_begin;
1600 ++values_begin;
1601 }
1602 }
1603
1604
1605
1606 inline Vector::iterator
1608 {
1609 return (*vector)[0];
1610 }
1611
1612
1613
1614 inline Vector::iterator
1615 Vector::end()
1616 {
1617 return (*vector)[0] + locally_owned_size();
1618 }
1619
1620
1621
1623 Vector::begin() const
1624 {
1625 return (*vector)[0];
1626 }
1627
1628
1629
1631 Vector::end() const
1632 {
1633 return (*vector)[0] + locally_owned_size();
1634 }
1635
1636
1637
1638 inline void
1639 Vector::set(const std::vector<size_type> &indices,
1640 const std::vector<TrilinosScalar> &values)
1641 {
1642 // if we have ghost values, do not allow
1643 // writing to this vector at all.
1645
1646 AssertDimension(indices.size(), values.size());
1647
1648 set(indices.size(), indices.data(), values.data());
1649 }
1650
1651
1652
1653 inline void
1654 Vector::set(const std::vector<size_type> &indices,
1655 const ::Vector<TrilinosScalar> &values)
1656 {
1657 // if we have ghost values, do not allow
1658 // writing to this vector at all.
1660
1661 AssertDimension(indices.size(), values.size());
1662
1663 set(indices.size(), indices.data(), values.begin());
1664 }
1665
1666
1667
1668 inline void
1669 Vector::set(const size_type n_elements,
1670 const size_type *indices,
1671 const TrilinosScalar *values)
1672 {
1673 // if we have ghost values, do not allow
1674 // writing to this vector at all.
1676
1677 if (last_action == Add)
1678 {
1679 const int ierr = vector->GlobalAssemble(Add);
1680 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1681 }
1682
1683 if (last_action != Insert)
1684 last_action = Insert;
1685
1686 for (size_type i = 0; i < n_elements; ++i)
1687 {
1688 const TrilinosWrappers::types::int_type row = indices[i];
1689 const TrilinosWrappers::types::int_type local_row =
1690 vector->Map().LID(row);
1691 if (local_row != -1)
1692 (*vector)[0][local_row] = values[i];
1693 else
1694 {
1695 const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1696 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1697 compressed = false;
1698 }
1699 // in set operation, do not use the pre-allocated vector for nonlocal
1700 // entries even if it exists. This is to ensure that we really only
1701 // set the elements touched by the set() method and not all contained
1702 // in the nonlocal entries vector (there is no way to distinguish them
1703 // on the receiving processor)
1704 }
1705 }
1706
1707
1708
1709 inline void
1710 Vector::add(const std::vector<size_type> &indices,
1711 const std::vector<TrilinosScalar> &values)
1712 {
1713 // if we have ghost values, do not allow
1714 // writing to this vector at all.
1716 AssertDimension(indices.size(), values.size());
1717
1718 add(indices.size(), indices.data(), values.data());
1719 }
1720
1721
1722
1723 inline void
1724 Vector::add(const std::vector<size_type> &indices,
1725 const ::Vector<TrilinosScalar> &values)
1726 {
1727 // if we have ghost values, do not allow
1728 // writing to this vector at all.
1730 AssertDimension(indices.size(), values.size());
1731
1732 add(indices.size(), indices.data(), values.begin());
1733 }
1734
1735
1736
1737 inline void
1738 Vector::add(const size_type n_elements,
1739 const size_type *indices,
1740 const TrilinosScalar *values)
1741 {
1742 // if we have ghost values, do not allow
1743 // writing to this vector at all.
1745
1746 if (last_action != Add)
1747 {
1748 if (last_action == Insert)
1749 {
1750 const int ierr = vector->GlobalAssemble(Insert);
1751 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1752 }
1753 last_action = Add;
1754 }
1755
1756 for (size_type i = 0; i < n_elements; ++i)
1757 {
1758 const size_type row = indices[i];
1759 const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1760 static_cast<TrilinosWrappers::types::int_type>(row));
1761 if (local_row != -1)
1762 (*vector)[0][local_row] += values[i];
1763 else if (nonlocal_vector.get() == nullptr)
1764 {
1765 const int ierr = vector->SumIntoGlobalValues(
1766 1,
1767 reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1768 &row),
1769 &values[i]);
1770 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1771 compressed = false;
1772 }
1773 else
1774 {
1775 // use pre-allocated vector for non-local entries if it exists for
1776 // addition operation
1778 nonlocal_vector->Map().LID(
1779 static_cast<TrilinosWrappers::types::int_type>(row));
1780 Assert(my_row != -1,
1781 ExcMessage(
1782 "Attempted to write into off-processor vector entry "
1783 "that has not be specified as being writable upon "
1784 "initialization"));
1785 (*nonlocal_vector)[0][my_row] += values[i];
1786 compressed = false;
1787 }
1788 }
1789 }
1790
1791
1792
1793 inline Vector::size_type
1794 Vector::size() const
1795 {
1796# ifndef DEAL_II_WITH_64BIT_INDICES
1797 return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1798# else
1799 return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1800# endif
1801 }
1802
1803
1804
1805 inline Vector::size_type
1807 {
1808 return owned_elements.n_elements();
1809 }
1810
1811
1812
1813 inline std::pair<Vector::size_type, Vector::size_type>
1814 Vector::local_range() const
1815 {
1816# ifndef DEAL_II_WITH_64BIT_INDICES
1817 const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1819 vector->Map().MaxMyGID() + 1;
1820# else
1822 vector->Map().MinMyGID64();
1824 vector->Map().MaxMyGID64() + 1;
1825# endif
1826
1827 Assert(
1828 end - begin == vector->Map().NumMyElements(),
1829 ExcMessage(
1830 "This function only makes sense if the elements that this "
1831 "vector stores on the current processor form a contiguous range. "
1832 "This does not appear to be the case for the current vector."));
1833
1834 return std::make_pair(begin, end);
1835 }
1836
1837
1838
1839 inline TrilinosScalar
1840 Vector::operator*(const Vector &vec) const
1841 {
1842 Assert(vector->Map().SameAs(vec.vector->Map()),
1845
1846 TrilinosScalar result;
1847
1848 const int ierr = vector->Dot(*(vec.vector), &result);
1849 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1850
1851 return result;
1852 }
1853
1854
1855
1856 inline Vector::real_type
1857 Vector::norm_sqr() const
1858 {
1859 const TrilinosScalar d = l2_norm();
1860 return d * d;
1861 }
1862
1863
1864
1865 inline TrilinosScalar
1866 Vector::mean_value() const
1867 {
1869
1871 const int ierr = vector->MeanValue(&mean);
1872 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1873
1874 return mean;
1875 }
1876
1877
1878
1879 inline TrilinosScalar
1880 Vector::min() const
1881 {
1882 TrilinosScalar min_value;
1883 const int ierr = vector->MinValue(&min_value);
1884 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1885
1886 return min_value;
1887 }
1888
1889
1890
1891 inline TrilinosScalar
1892 Vector::max() const
1893 {
1894 TrilinosScalar max_value;
1895 const int ierr = vector->MaxValue(&max_value);
1896 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1897
1898 return max_value;
1899 }
1900
1901
1902
1903 inline Vector::real_type
1904 Vector::l1_norm() const
1905 {
1907
1909 const int ierr = vector->Norm1(&d);
1910 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1911
1912 return d;
1913 }
1914
1915
1916
1917 inline Vector::real_type
1918 Vector::l2_norm() const
1919 {
1921
1923 const int ierr = vector->Norm2(&d);
1924 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1925
1926 return d;
1927 }
1928
1929
1930
1931 inline Vector::real_type
1932 Vector::lp_norm(const TrilinosScalar p) const
1933 {
1935
1936 TrilinosScalar norm = 0;
1937 TrilinosScalar sum = 0;
1938 const size_type n_local = locally_owned_size();
1939
1940 // loop over all the elements because
1941 // Trilinos does not support lp norms
1942 for (size_type i = 0; i < n_local; ++i)
1943 sum += std::pow(std::fabs((*vector)[0][i]), p);
1944
1945 norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1946
1947 return norm;
1948 }
1949
1950
1951
1952 inline Vector::real_type
1953 Vector::linfty_norm() const
1954 {
1955 // while we disallow the other
1956 // norm operations on ghosted
1957 // vectors, this particular norm
1958 // is safe to run even in the
1959 // presence of ghost elements
1961 const int ierr = vector->NormInf(&d);
1962 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1963
1964 return d;
1965 }
1966
1967
1968
1969 inline TrilinosScalar
1971 const Vector &V,
1972 const Vector &W)
1973 {
1974 this->add(a, V);
1975 return *this * W;
1976 }
1977
1978
1979
1980 // inline also scalar products, vector
1981 // additions etc. since they are all
1982 // representable by a single Trilinos
1983 // call. This reduces the overhead of the
1984 // wrapper class.
1985 inline Vector &
1987 {
1988 AssertIsFinite(a);
1989
1990 const int ierr = vector->Scale(a);
1991 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1992
1993 return *this;
1994 }
1995
1996
1997
1998 inline Vector &
2000 {
2001 AssertIsFinite(a);
2002
2003 const TrilinosScalar factor = 1. / a;
2004
2005 AssertIsFinite(factor);
2006
2007 const int ierr = vector->Scale(factor);
2008 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2009
2010 return *this;
2011 }
2012
2013
2014
2015 inline Vector &
2016 Vector::operator+=(const Vector &v)
2017 {
2018 AssertDimension(size(), v.size());
2019 Assert(vector->Map().SameAs(v.vector->Map()),
2021
2022 const int ierr = vector->Update(1.0, *(v.vector), 1.0);
2023 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2024
2025 return *this;
2026 }
2027
2028
2029
2030 inline Vector &
2031 Vector::operator-=(const Vector &v)
2032 {
2033 AssertDimension(size(), v.size());
2034 Assert(vector->Map().SameAs(v.vector->Map()),
2036
2037 const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
2038 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2039
2040 return *this;
2041 }
2042
2043
2044
2045 inline void
2047 {
2048 // if we have ghost values, do not allow
2049 // writing to this vector at all.
2051 AssertIsFinite(s);
2052
2053 size_type n_local = locally_owned_size();
2054 for (size_type i = 0; i < n_local; ++i)
2055 (*vector)[0][i] += s;
2056 }
2057
2058
2059
2060 inline void
2061 Vector::add(const TrilinosScalar a, const Vector &v)
2062 {
2063 // if we have ghost values, do not allow
2064 // writing to this vector at all.
2067
2068 AssertIsFinite(a);
2069
2070 const int ierr = vector->Update(a, *(v.vector), 1.);
2071 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2072 }
2073
2074
2075
2076 inline void
2078 const Vector &v,
2079 const TrilinosScalar b,
2080 const Vector &w)
2081 {
2082 // if we have ghost values, do not allow
2083 // writing to this vector at all.
2086 AssertDimension(locally_owned_size(), w.locally_owned_size());
2087
2088 AssertIsFinite(a);
2089 AssertIsFinite(b);
2090
2091 const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2092
2093 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2094 }
2095
2096
2097
2098 inline void
2099 Vector::sadd(const TrilinosScalar s, const Vector &v)
2100 {
2101 // if we have ghost values, do not allow
2102 // writing to this vector at all.
2104 AssertDimension(size(), v.size());
2105
2106 AssertIsFinite(s);
2107
2108 // We assume that the vectors have the same Map
2109 // if the local size is the same and if the vectors are not ghosted
2111 !v.has_ghost_elements())
2112 {
2113 Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2115 const int ierr = vector->Update(1., *(v.vector), s);
2116 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2117 }
2118 else
2119 {
2120 (*this) *= s;
2121 this->add(v, true);
2122 }
2123 }
2124
2125
2126
2127 inline void
2129 const TrilinosScalar a,
2130 const Vector &v)
2131 {
2132 // if we have ghost values, do not allow
2133 // writing to this vector at all.
2135 AssertDimension(size(), v.size());
2136 AssertIsFinite(s);
2137 AssertIsFinite(a);
2138
2139 // We assume that the vectors have the same Map
2140 // if the local size is the same and if the vectors are not ghosted
2142 !v.has_ghost_elements())
2143 {
2144 Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2146 const int ierr = vector->Update(a, *(v.vector), s);
2147 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2148 }
2149 else
2150 {
2151 (*this) *= s;
2152 Vector tmp = v;
2153 tmp *= a;
2154 this->add(tmp, true);
2155 }
2156 }
2157
2158
2159
2160 inline void
2161 Vector::scale(const Vector &factors)
2162 {
2163 // if we have ghost values, do not allow
2164 // writing to this vector at all.
2167
2168 const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2169 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2170 }
2171
2172
2173
2174 inline void
2175 Vector::equ(const TrilinosScalar a, const Vector &v)
2176 {
2177 // if we have ghost values, do not allow
2178 // writing to this vector at all.
2180 AssertIsFinite(a);
2181
2182 // If we don't have the same map, copy.
2183 if (vector->Map().SameAs(v.vector->Map()) == false)
2184 {
2185 this->sadd(0., a, v);
2186 }
2187 else
2188 {
2189 // Otherwise, just update
2190 int ierr = vector->Update(a, *v.vector, 0.0);
2191 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2192
2193 last_action = Zero;
2194 }
2195 }
2196
2197
2198
2199 inline const Epetra_MultiVector &
2201 {
2202 return static_cast<const Epetra_MultiVector &>(*vector);
2203 }
2204
2205
2206
2207 inline Epetra_FEVector &
2209 {
2210 return *vector;
2211 }
2212
2213
2214
2215 inline const Epetra_BlockMap &
2217 {
2218 return vector->Map();
2219 }
2220
2221
2222
2223 inline MPI_Comm
2225 {
2226 const Epetra_MpiComm *mpi_comm =
2227 dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2228 return mpi_comm->Comm();
2229 }
2230
2231
2232
2233 template <typename number>
2234 Vector::Vector(const IndexSet &parallel_partitioner,
2235 const ::Vector<number> &v,
2236 const MPI_Comm communicator)
2237 {
2238 *this =
2239 Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2240 owned_elements = parallel_partitioner;
2241 }
2242
2243
2244
2245 inline Vector &
2247 {
2248 AssertIsFinite(s);
2249
2250 int ierr = vector->PutScalar(s);
2251 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2252
2253 if (nonlocal_vector.get() != nullptr)
2254 {
2255 ierr = nonlocal_vector->PutScalar(0.);
2256 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2257 }
2258
2259 return *this;
2260 }
2261 } /* end of namespace MPI */
2262
2263# endif /* DOXYGEN */
2264
2265} /* end of namespace TrilinosWrappers */
2266
2270namespace internal
2271{
2272 namespace LinearOperatorImplementation
2273 {
2274 template <typename>
2275 class ReinitHelper;
2276
2281 template <>
2282 class ReinitHelper<TrilinosWrappers::MPI::Vector>
2283 {
2284 public:
2285 template <typename Matrix>
2286 static void
2287 reinit_range_vector(const Matrix &matrix,
2289 bool omit_zeroing_entries)
2290 {
2291 v.reinit(matrix.locally_owned_range_indices(),
2292 matrix.get_mpi_communicator(),
2293 omit_zeroing_entries);
2294 }
2295
2296 template <typename Matrix>
2297 static void
2298 reinit_domain_vector(const Matrix &matrix,
2300 bool omit_zeroing_entries)
2301 {
2302 v.reinit(matrix.locally_owned_domain_indices(),
2303 matrix.get_mpi_communicator(),
2304 omit_zeroing_entries);
2305 }
2306 };
2307
2308 } // namespace LinearOperatorImplementation
2309} /* namespace internal */
2310
2311
2312
2316template <>
2317struct is_serial_vector<TrilinosWrappers::MPI::Vector> : std::false_type
2318{};
2319
2320
2322
2323#endif
2324
2325#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)
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)
::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: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)
types::global_dof_index size_type
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