Reference documentation for deal.II version GIT relicensing-12-g89ba5bf9db 2024-02-24 00:00:03+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
678 Vector &
679 operator=(const Vector &v);
680
685 Vector &
686 operator=(Vector &&v) noexcept;
687
695 template <typename Number>
696 Vector &
697 operator=(const ::Vector<Number> &v);
698
716 void
718 const ::TrilinosWrappers::SparseMatrix &matrix,
719 const Vector &vector);
720
727 void
729 const VectorOperation::values operation);
730
735 void
737 const VectorOperation::values operation)
738 {
739 import_elements(rwv, operation);
740 }
741
742
748 bool
749 operator==(const Vector &v) const;
750
756 bool
757 operator!=(const Vector &v) const;
758
763 size() const override;
764
771
794 std::pair<size_type, size_type>
795 local_range() const;
796
804 bool
805 in_local_range(const size_type index) const;
806
822
830 bool
832
839 void
841
847 operator*(const Vector &vec) const;
848
853 norm_sqr() const;
854
859 mean_value() const;
860
865 min() const;
866
871 max() const;
872
877 l1_norm() const;
878
884 l2_norm() const;
885
891 lp_norm(const TrilinosScalar p) const;
892
897 linfty_norm() const;
898
919 add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
920
926 bool
927 all_zero() const;
928
934 bool
935 is_non_negative() const;
952 operator()(const size_type index);
953
962 operator()(const size_type index) const;
963
970 operator[](const size_type index);
971
978 operator[](const size_type index) const;
979
995 void
996 extract_subvector_to(const std::vector<size_type> &indices,
997 std::vector<TrilinosScalar> &values) const;
998
1002 virtual void
1004 ArrayView<TrilinosScalar> &elements) const override;
1005
1033 template <typename ForwardIterator, typename OutputIterator>
1034 void
1035 extract_subvector_to(ForwardIterator indices_begin,
1036 const ForwardIterator indices_end,
1037 OutputIterator values_begin) const;
1038
1047 iterator
1049
1055 begin() const;
1056
1061 iterator
1063
1069 end() const;
1070
1085 void
1086 set(const std::vector<size_type> &indices,
1087 const std::vector<TrilinosScalar> &values);
1088
1093 void
1094 set(const std::vector<size_type> &indices,
1095 const ::Vector<TrilinosScalar> &values);
1096
1102 void
1103 set(const size_type n_elements,
1104 const size_type *indices,
1105 const TrilinosScalar *values);
1106
1111 void
1112 add(const std::vector<size_type> &indices,
1113 const std::vector<TrilinosScalar> &values);
1114
1119 void
1120 add(const std::vector<size_type> &indices,
1121 const ::Vector<TrilinosScalar> &values);
1122
1128 void
1129 add(const size_type n_elements,
1130 const size_type *indices,
1131 const TrilinosScalar *values);
1132
1136 Vector &
1138
1142 Vector &
1144
1148 Vector &
1150
1154 Vector &
1156
1161 void
1163
1176 void
1177 add(const Vector &V, const bool allow_different_maps = false);
1178
1182 void
1183 add(const TrilinosScalar a, const Vector &V);
1184
1188 void
1190 const Vector &V,
1191 const TrilinosScalar b,
1192 const Vector &W);
1193
1198 void
1199 sadd(const TrilinosScalar s, const Vector &V);
1200
1204 void
1205 sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1206
1212 void
1213 scale(const Vector &scaling_factors);
1214
1218 void
1219 equ(const TrilinosScalar a, const Vector &V);
1231 const Epetra_MultiVector &
1233
1238 Epetra_FEVector &
1240
1245 const Epetra_BlockMap &
1247
1255 void
1256 print(std::ostream &out,
1257 const unsigned int precision = 3,
1258 const bool scientific = true,
1259 const bool across = true) const;
1260
1274 void
1275 swap(Vector &v);
1276
1280 std::size_t
1281 memory_consumption() const;
1282
1286 MPI_Comm
1294
1299 int,
1300 << "An error with error number " << arg1
1301 << " occurred while calling a Trilinos function");
1302
1308 size_type,
1309 size_type,
1310 size_type,
1311 size_type,
1312 << "You are trying to access element " << arg1
1313 << " of a distributed vector, but this element is not stored "
1314 << "on the current processor. Note: There are " << arg2
1315 << " elements stored "
1316 << "on the current processor from within the range [" << arg3 << ','
1317 << arg4 << "] but Trilinos vectors need not store contiguous "
1318 << "ranges on each processor, and not every element in "
1319 << "this range may in fact be stored locally."
1320 << "\n\n"
1321 << "A common source for this kind of problem is that you "
1322 << "are passing a 'fully distributed' vector into a function "
1323 << "that needs read access to vector elements that correspond "
1324 << "to degrees of freedom on ghost cells (or at least to "
1325 << "'locally active' degrees of freedom that are not also "
1326 << "'locally owned'). You need to pass a vector that has these "
1327 << "elements as ghost entries.");
1328
1329 private:
1341 Epetra_CombineMode last_action;
1342
1348
1354
1360 std::unique_ptr<Epetra_FEVector> vector;
1361
1367 std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1368
1373
1374 // Make the reference class a friend.
1376 };
1377
1378
1379
1380 // ------------------- inline and template functions --------------
1381
1382
1390 inline void
1392 {
1393 u.swap(v);
1394 }
1395 } // namespace MPI
1396
1397# ifndef DOXYGEN
1398
1399 namespace internal
1400 {
1401 inline VectorReference::VectorReference(MPI::Vector &vector,
1402 const size_type index)
1403 : vector(vector)
1404 , index(index)
1405 {}
1406
1407
1408 inline const VectorReference &
1409 VectorReference::operator=(const VectorReference &r) const
1410 {
1411 // as explained in the class
1412 // documentation, this is not the copy
1413 // operator. so simply pass on to the
1414 // "correct" assignment operator
1415 *this = static_cast<TrilinosScalar>(r);
1416
1417 return *this;
1418 }
1419
1420
1421
1422 inline VectorReference &
1423 VectorReference::operator=(const VectorReference &r)
1424 {
1425 // as above
1426 *this = static_cast<TrilinosScalar>(r);
1427
1428 return *this;
1429 }
1430
1431
1432 inline const VectorReference &
1433 VectorReference::operator=(const TrilinosScalar &value) const
1434 {
1435 vector.set(1, &index, &value);
1436 return *this;
1437 }
1438
1439
1440
1441 inline const VectorReference &
1442 VectorReference::operator+=(const TrilinosScalar &value) const
1443 {
1444 vector.add(1, &index, &value);
1445 return *this;
1446 }
1447
1448
1449
1450 inline const VectorReference &
1451 VectorReference::operator-=(const TrilinosScalar &value) const
1452 {
1453 TrilinosScalar new_value = -value;
1454 vector.add(1, &index, &new_value);
1455 return *this;
1456 }
1457
1458
1459
1460 inline const VectorReference &
1461 VectorReference::operator*=(const TrilinosScalar &value) const
1462 {
1463 TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1464 vector.set(1, &index, &new_value);
1465 return *this;
1466 }
1467
1468
1469
1470 inline const VectorReference &
1471 VectorReference::operator/=(const TrilinosScalar &value) const
1472 {
1473 TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1474 vector.set(1, &index, &new_value);
1475 return *this;
1476 }
1477 } // namespace internal
1478
1479 namespace MPI
1480 {
1481 inline bool
1482 Vector::in_local_range(const size_type index) const
1483 {
1484 std::pair<size_type, size_type> range = local_range();
1485
1486 return ((index >= range.first) && (index < range.second));
1487 }
1488
1489
1490
1491 inline IndexSet
1493 {
1495 ExcMessage(
1496 "The locally owned elements have not been properly initialized!"
1497 " This happens for example if this object has been initialized"
1498 " with exactly one overlapping IndexSet."));
1499 return owned_elements;
1500 }
1501
1502
1503
1504 inline bool
1506 {
1507 return has_ghosts;
1508 }
1509
1510
1511
1512 inline void
1514 {}
1515
1516
1517
1518 inline internal::VectorReference
1519 Vector::operator()(const size_type index)
1520 {
1521 return internal::VectorReference(*this, index);
1522 }
1523
1524
1525
1526 inline internal::VectorReference
1527 Vector::operator[](const size_type index)
1528 {
1529 return operator()(index);
1530 }
1531
1532
1533
1534 inline TrilinosScalar
1535 Vector::operator[](const size_type index) const
1536 {
1537 return operator()(index);
1538 }
1539
1540
1541
1542 inline void
1543 Vector::extract_subvector_to(const std::vector<size_type> &indices,
1544 std::vector<TrilinosScalar> &values) const
1545 {
1546 for (size_type i = 0; i < indices.size(); ++i)
1547 values[i] = operator()(indices[i]);
1548 }
1549
1550
1551
1552 inline void
1554 ArrayView<TrilinosScalar> &elements) const
1555 {
1556 AssertDimension(indices.size(), elements.size());
1557 for (unsigned int i = 0; i < indices.size(); ++i)
1558 {
1559 AssertIndexRange(indices[i], size());
1560 elements[i] = (*this)[indices[i]];
1561 }
1562 }
1563
1564
1565
1566 template <typename ForwardIterator, typename OutputIterator>
1567 inline void
1568 Vector::extract_subvector_to(ForwardIterator indices_begin,
1569 const ForwardIterator indices_end,
1570 OutputIterator values_begin) const
1571 {
1572 while (indices_begin != indices_end)
1573 {
1574 *values_begin = operator()(*indices_begin);
1575 ++indices_begin;
1576 ++values_begin;
1577 }
1578 }
1579
1580
1581
1582 inline Vector::iterator
1584 {
1585 return (*vector)[0];
1586 }
1587
1588
1589
1590 inline Vector::iterator
1591 Vector::end()
1592 {
1593 return (*vector)[0] + locally_owned_size();
1594 }
1595
1596
1597
1599 Vector::begin() const
1600 {
1601 return (*vector)[0];
1602 }
1603
1604
1605
1607 Vector::end() const
1608 {
1609 return (*vector)[0] + locally_owned_size();
1610 }
1611
1612
1613
1614 inline void
1615 Vector::set(const std::vector<size_type> &indices,
1616 const std::vector<TrilinosScalar> &values)
1617 {
1618 // if we have ghost values, do not allow
1619 // writing to this vector at all.
1621
1622 AssertDimension(indices.size(), values.size());
1623
1624 set(indices.size(), indices.data(), values.data());
1625 }
1626
1627
1628
1629 inline void
1630 Vector::set(const std::vector<size_type> &indices,
1631 const ::Vector<TrilinosScalar> &values)
1632 {
1633 // if we have ghost values, do not allow
1634 // writing to this vector at all.
1636
1637 AssertDimension(indices.size(), values.size());
1638
1639 set(indices.size(), indices.data(), values.begin());
1640 }
1641
1642
1643
1644 inline void
1645 Vector::set(const size_type n_elements,
1646 const size_type *indices,
1647 const TrilinosScalar *values)
1648 {
1649 // if we have ghost values, do not allow
1650 // writing to this vector at all.
1652
1653 if (last_action == Add)
1654 {
1655 const int ierr = vector->GlobalAssemble(Add);
1656 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1657 }
1658
1659 if (last_action != Insert)
1660 last_action = Insert;
1661
1662 for (size_type i = 0; i < n_elements; ++i)
1663 {
1664 const TrilinosWrappers::types::int_type row = indices[i];
1665 const TrilinosWrappers::types::int_type local_row =
1666 vector->Map().LID(row);
1667 if (local_row != -1)
1668 (*vector)[0][local_row] = values[i];
1669 else
1670 {
1671 const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1672 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1673 compressed = false;
1674 }
1675 // in set operation, do not use the pre-allocated vector for nonlocal
1676 // entries even if it exists. This is to ensure that we really only
1677 // set the elements touched by the set() method and not all contained
1678 // in the nonlocal entries vector (there is no way to distinguish them
1679 // on the receiving processor)
1680 }
1681 }
1682
1683
1684
1685 inline void
1686 Vector::add(const std::vector<size_type> &indices,
1687 const std::vector<TrilinosScalar> &values)
1688 {
1689 // if we have ghost values, do not allow
1690 // writing to this vector at all.
1692 AssertDimension(indices.size(), values.size());
1693
1694 add(indices.size(), indices.data(), values.data());
1695 }
1696
1697
1698
1699 inline void
1700 Vector::add(const std::vector<size_type> &indices,
1701 const ::Vector<TrilinosScalar> &values)
1702 {
1703 // if we have ghost values, do not allow
1704 // writing to this vector at all.
1706 AssertDimension(indices.size(), values.size());
1707
1708 add(indices.size(), indices.data(), values.begin());
1709 }
1710
1711
1712
1713 inline void
1714 Vector::add(const size_type n_elements,
1715 const size_type *indices,
1716 const TrilinosScalar *values)
1717 {
1718 // if we have ghost values, do not allow
1719 // writing to this vector at all.
1721
1722 if (last_action != Add)
1723 {
1724 if (last_action == Insert)
1725 {
1726 const int ierr = vector->GlobalAssemble(Insert);
1727 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1728 }
1729 last_action = Add;
1730 }
1731
1732 for (size_type i = 0; i < n_elements; ++i)
1733 {
1734 const size_type row = indices[i];
1735 const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1736 static_cast<TrilinosWrappers::types::int_type>(row));
1737 if (local_row != -1)
1738 (*vector)[0][local_row] += values[i];
1739 else if (nonlocal_vector.get() == nullptr)
1740 {
1741 const int ierr = vector->SumIntoGlobalValues(
1742 1,
1743 reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1744 &row),
1745 &values[i]);
1746 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1747 compressed = false;
1748 }
1749 else
1750 {
1751 // use pre-allocated vector for non-local entries if it exists for
1752 // addition operation
1754 nonlocal_vector->Map().LID(
1755 static_cast<TrilinosWrappers::types::int_type>(row));
1756 Assert(my_row != -1,
1757 ExcMessage(
1758 "Attempted to write into off-processor vector entry "
1759 "that has not be specified as being writable upon "
1760 "initialization"));
1761 (*nonlocal_vector)[0][my_row] += values[i];
1762 compressed = false;
1763 }
1764 }
1765 }
1766
1767
1768
1769 inline Vector::size_type
1770 Vector::size() const
1771 {
1772# ifndef DEAL_II_WITH_64BIT_INDICES
1773 return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1774# else
1775 return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1776# endif
1777 }
1778
1779
1780
1781 inline Vector::size_type
1783 {
1784 return owned_elements.n_elements();
1785 }
1786
1787
1788
1789 inline std::pair<Vector::size_type, Vector::size_type>
1790 Vector::local_range() const
1791 {
1792# ifndef DEAL_II_WITH_64BIT_INDICES
1793 const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1795 vector->Map().MaxMyGID() + 1;
1796# else
1798 vector->Map().MinMyGID64();
1800 vector->Map().MaxMyGID64() + 1;
1801# endif
1802
1803 Assert(
1804 end - begin == vector->Map().NumMyElements(),
1805 ExcMessage(
1806 "This function only makes sense if the elements that this "
1807 "vector stores on the current processor form a contiguous range. "
1808 "This does not appear to be the case for the current vector."));
1809
1810 return std::make_pair(begin, end);
1811 }
1812
1813
1814
1815 inline TrilinosScalar
1816 Vector::operator*(const Vector &vec) const
1817 {
1818 Assert(vector->Map().SameAs(vec.vector->Map()),
1821
1822 TrilinosScalar result;
1823
1824 const int ierr = vector->Dot(*(vec.vector), &result);
1825 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1826
1827 return result;
1828 }
1829
1830
1831
1832 inline Vector::real_type
1833 Vector::norm_sqr() const
1834 {
1835 const TrilinosScalar d = l2_norm();
1836 return d * d;
1837 }
1838
1839
1840
1841 inline TrilinosScalar
1842 Vector::mean_value() const
1843 {
1845
1847 const int ierr = vector->MeanValue(&mean);
1848 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1849
1850 return mean;
1851 }
1852
1853
1854
1855 inline TrilinosScalar
1856 Vector::min() const
1857 {
1858 TrilinosScalar min_value;
1859 const int ierr = vector->MinValue(&min_value);
1860 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1861
1862 return min_value;
1863 }
1864
1865
1866
1867 inline TrilinosScalar
1868 Vector::max() const
1869 {
1870 TrilinosScalar max_value;
1871 const int ierr = vector->MaxValue(&max_value);
1872 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1873
1874 return max_value;
1875 }
1876
1877
1878
1879 inline Vector::real_type
1880 Vector::l1_norm() const
1881 {
1883
1885 const int ierr = vector->Norm1(&d);
1886 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1887
1888 return d;
1889 }
1890
1891
1892
1893 inline Vector::real_type
1894 Vector::l2_norm() const
1895 {
1897
1899 const int ierr = vector->Norm2(&d);
1900 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1901
1902 return d;
1903 }
1904
1905
1906
1907 inline Vector::real_type
1908 Vector::lp_norm(const TrilinosScalar p) const
1909 {
1911
1912 TrilinosScalar norm = 0;
1913 TrilinosScalar sum = 0;
1914 const size_type n_local = locally_owned_size();
1915
1916 // loop over all the elements because
1917 // Trilinos does not support lp norms
1918 for (size_type i = 0; i < n_local; ++i)
1919 sum += std::pow(std::fabs((*vector)[0][i]), p);
1920
1921 norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1922
1923 return norm;
1924 }
1925
1926
1927
1928 inline Vector::real_type
1929 Vector::linfty_norm() const
1930 {
1931 // while we disallow the other
1932 // norm operations on ghosted
1933 // vectors, this particular norm
1934 // is safe to run even in the
1935 // presence of ghost elements
1937 const int ierr = vector->NormInf(&d);
1938 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1939
1940 return d;
1941 }
1942
1943
1944
1945 inline TrilinosScalar
1947 const Vector &V,
1948 const Vector &W)
1949 {
1950 this->add(a, V);
1951 return *this * W;
1952 }
1953
1954
1955
1956 // inline also scalar products, vector
1957 // additions etc. since they are all
1958 // representable by a single Trilinos
1959 // call. This reduces the overhead of the
1960 // wrapper class.
1961 inline Vector &
1963 {
1964 AssertIsFinite(a);
1965
1966 const int ierr = vector->Scale(a);
1967 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1968
1969 return *this;
1970 }
1971
1972
1973
1974 inline Vector &
1976 {
1977 AssertIsFinite(a);
1978
1979 const TrilinosScalar factor = 1. / a;
1980
1981 AssertIsFinite(factor);
1982
1983 const int ierr = vector->Scale(factor);
1984 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1985
1986 return *this;
1987 }
1988
1989
1990
1991 inline Vector &
1992 Vector::operator+=(const Vector &v)
1993 {
1994 AssertDimension(size(), v.size());
1995 Assert(vector->Map().SameAs(v.vector->Map()),
1997
1998 const int ierr = vector->Update(1.0, *(v.vector), 1.0);
1999 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2000
2001 return *this;
2002 }
2003
2004
2005
2006 inline Vector &
2007 Vector::operator-=(const Vector &v)
2008 {
2009 AssertDimension(size(), v.size());
2010 Assert(vector->Map().SameAs(v.vector->Map()),
2012
2013 const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
2014 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2015
2016 return *this;
2017 }
2018
2019
2020
2021 inline void
2023 {
2024 // if we have ghost values, do not allow
2025 // writing to this vector at all.
2027 AssertIsFinite(s);
2028
2029 size_type n_local = locally_owned_size();
2030 for (size_type i = 0; i < n_local; ++i)
2031 (*vector)[0][i] += s;
2032 }
2033
2034
2035
2036 inline void
2037 Vector::add(const TrilinosScalar a, const Vector &v)
2038 {
2039 // if we have ghost values, do not allow
2040 // writing to this vector at all.
2043
2044 AssertIsFinite(a);
2045
2046 const int ierr = vector->Update(a, *(v.vector), 1.);
2047 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2048 }
2049
2050
2051
2052 inline void
2054 const Vector &v,
2055 const TrilinosScalar b,
2056 const Vector &w)
2057 {
2058 // if we have ghost values, do not allow
2059 // writing to this vector at all.
2062 AssertDimension(locally_owned_size(), w.locally_owned_size());
2063
2064 AssertIsFinite(a);
2065 AssertIsFinite(b);
2066
2067 const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2068
2069 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2070 }
2071
2072
2073
2074 inline void
2075 Vector::sadd(const TrilinosScalar s, const Vector &v)
2076 {
2077 // if we have ghost values, do not allow
2078 // writing to this vector at all.
2080 AssertDimension(size(), v.size());
2081
2082 AssertIsFinite(s);
2083
2084 // We assume that the vectors have the same Map
2085 // if the local size is the same and if the vectors are not ghosted
2087 !v.has_ghost_elements())
2088 {
2089 Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2091 const int ierr = vector->Update(1., *(v.vector), s);
2092 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2093 }
2094 else
2095 {
2096 (*this) *= s;
2097 this->add(v, true);
2098 }
2099 }
2100
2101
2102
2103 inline void
2105 const TrilinosScalar a,
2106 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 AssertIsFinite(s);
2113 AssertIsFinite(a);
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(a, *(v.vector), s);
2123 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2124 }
2125 else
2126 {
2127 (*this) *= s;
2128 Vector tmp = v;
2129 tmp *= a;
2130 this->add(tmp, true);
2131 }
2132 }
2133
2134
2135
2136 inline void
2137 Vector::scale(const Vector &factors)
2138 {
2139 // if we have ghost values, do not allow
2140 // writing to this vector at all.
2143
2144 const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2145 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2146 }
2147
2148
2149
2150 inline void
2151 Vector::equ(const TrilinosScalar a, const Vector &v)
2152 {
2153 // if we have ghost values, do not allow
2154 // writing to this vector at all.
2156 AssertIsFinite(a);
2157
2158 // If we don't have the same map, copy.
2159 if (vector->Map().SameAs(v.vector->Map()) == false)
2160 {
2161 this->sadd(0., a, v);
2162 }
2163 else
2164 {
2165 // Otherwise, just update
2166 int ierr = vector->Update(a, *v.vector, 0.0);
2167 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2168
2169 last_action = Zero;
2170 }
2171 }
2172
2173
2174
2175 inline const Epetra_MultiVector &
2177 {
2178 return static_cast<const Epetra_MultiVector &>(*vector);
2179 }
2180
2181
2182
2183 inline Epetra_FEVector &
2185 {
2186 return *vector;
2187 }
2188
2189
2190
2191 inline const Epetra_BlockMap &
2193 {
2194 return vector->Map();
2195 }
2196
2197
2198
2199 inline MPI_Comm
2201 {
2202 const Epetra_MpiComm *mpi_comm =
2203 dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2204 return mpi_comm->Comm();
2205 }
2206
2207
2208
2209 template <typename number>
2210 Vector::Vector(const IndexSet &parallel_partitioner,
2211 const ::Vector<number> &v,
2212 const MPI_Comm communicator)
2213 {
2214 *this =
2215 Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2216 owned_elements = parallel_partitioner;
2217 }
2218
2219
2220
2221 inline Vector &
2223 {
2224 AssertIsFinite(s);
2225
2226 int ierr = vector->PutScalar(s);
2227 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2228
2229 if (nonlocal_vector.get() != nullptr)
2230 {
2231 ierr = nonlocal_vector->PutScalar(0.);
2232 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2233 }
2234
2235 return *this;
2236 }
2237 } /* end of namespace MPI */
2238
2239# endif /* DOXYGEN */
2240
2241} /* end of namespace TrilinosWrappers */
2242
2246namespace internal
2247{
2248 namespace LinearOperatorImplementation
2249 {
2250 template <typename>
2251 class ReinitHelper;
2252
2257 template <>
2258 class ReinitHelper<TrilinosWrappers::MPI::Vector>
2259 {
2260 public:
2261 template <typename Matrix>
2262 static void
2263 reinit_range_vector(const Matrix &matrix,
2265 bool omit_zeroing_entries)
2266 {
2267 v.reinit(matrix.locally_owned_range_indices(),
2268 matrix.get_mpi_communicator(),
2269 omit_zeroing_entries);
2270 }
2271
2272 template <typename Matrix>
2273 static void
2274 reinit_domain_vector(const Matrix &matrix,
2276 bool omit_zeroing_entries)
2277 {
2278 v.reinit(matrix.locally_owned_domain_indices(),
2279 matrix.get_mpi_communicator(),
2280 omit_zeroing_entries);
2281 }
2282 };
2283
2284 } // namespace LinearOperatorImplementation
2285} /* namespace internal */
2286
2287
2288
2292template <>
2293struct is_serial_vector<TrilinosWrappers::MPI::Vector> : std::false_type
2294{};
2295
2296
2298
2299#endif
2300
2301#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