Reference documentation for deal.II version 9.4.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_vector.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2008 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_trilinos_vector_h
17#define dealii_trilinos_vector_h
18
19
20#include <deal.II/base/config.h>
21
22#ifdef DEAL_II_WITH_TRILINOS
24# include <deal.II/base/mpi.h>
27
29# include <deal.II/lac/vector.h>
32
33# include <Epetra_ConfigDefs.h>
34# include <Epetra_FEVector.h>
35# include <Epetra_LocalMap.h>
36# include <Epetra_Map.h>
37# include <Epetra_MpiComm.h>
38# include <mpi.h>
39
40# include <memory>
41# include <utility>
42# include <vector>
43
45
46// Forward declarations
47# ifndef DOXYGEN
48namespace LinearAlgebra
49{
50 // Forward declaration
51 template <typename Number>
52 class ReadWriteVector;
53} // namespace LinearAlgebra
54# endif
55
66namespace TrilinosWrappers
67{
68 class SparseMatrix;
69
80 namespace internal
81 {
86
96 class VectorReference
97 {
98 private:
103 VectorReference(MPI::Vector &vector, const size_type index);
104
105 public:
109 VectorReference(const VectorReference &) = default;
110
122 const VectorReference &
123 operator=(const VectorReference &r) const;
124
128 VectorReference &
129 operator=(const VectorReference &r);
130
134 const VectorReference &
135 operator=(const TrilinosScalar &s) const;
136
140 const VectorReference &
141 operator+=(const TrilinosScalar &s) const;
142
146 const VectorReference &
147 operator-=(const TrilinosScalar &s) const;
148
152 const VectorReference &
153 operator*=(const TrilinosScalar &s) const;
154
158 const VectorReference &
159 operator/=(const TrilinosScalar &s) const;
160
165 operator TrilinosScalar() const;
166
171 int,
172 << "An error with error number " << arg1
173 << " occurred while calling a Trilinos function");
174
175 private:
179 MPI::Vector &vector;
180
184 const size_type index;
185
186 // Make the vector class a friend, so that it can create objects of the
187 // present type.
189 };
190 } // namespace internal
195# ifndef DEAL_II_WITH_64BIT_INDICES
196 // define a helper function that queries the global ID of local ID of
197 // an Epetra_BlockMap object by calling either the 32- or 64-bit
198 // function necessary.
199 inline int
200 gid(const Epetra_BlockMap &map, int i)
201 {
202 return map.GID(i);
203 }
204# else
205 // define a helper function that queries the global ID of local ID of
206 // an Epetra_BlockMap object by calling either the 32- or 64-bit
207 // function necessary.
208 inline long long int
209 gid(const Epetra_BlockMap &map, int i)
210 {
211 return map.GID64(i);
212 }
213# endif
214
220 namespace MPI
221 {
222 class BlockVector;
223
398 class Vector : public Subscriptor
399 {
400 public:
410 using const_iterator = const value_type *;
413
423 Vector();
424
428 Vector(const Vector &v);
429
448 explicit Vector(const IndexSet &parallel_partitioning,
449 const MPI_Comm &communicator = MPI_COMM_WORLD);
450
462 Vector(const IndexSet &local,
463 const IndexSet &ghost,
464 const MPI_Comm &communicator = MPI_COMM_WORLD);
465
480 Vector(const IndexSet &parallel_partitioning,
481 const Vector & v,
482 const MPI_Comm &communicator = MPI_COMM_WORLD);
483
496 template <typename Number>
497 Vector(const IndexSet & parallel_partitioning,
498 const ::Vector<Number> &v,
499 const MPI_Comm & communicator = MPI_COMM_WORLD);
500
505 Vector(Vector &&v) noexcept;
506
510 ~Vector() override = default;
511
516 void
517 clear();
518
542 void
543 reinit(const Vector &v,
544 const bool omit_zeroing_entries = false,
545 const bool allow_different_maps = false);
546
569 void
570 reinit(const IndexSet &parallel_partitioning,
571 const MPI_Comm &communicator = MPI_COMM_WORLD,
572 const bool omit_zeroing_entries = false);
573
599 void
600 reinit(const IndexSet &locally_owned_entries,
601 const IndexSet &ghost_entries,
602 const MPI_Comm &communicator = MPI_COMM_WORLD,
603 const bool vector_writable = false);
604
608 void
609 reinit(const BlockVector &v, const bool import_data = false);
610
627 void
629
642 Vector &
644
650 Vector &
651 operator=(const Vector &v);
652
657 Vector &
658 operator=(Vector &&v) noexcept;
659
667 template <typename Number>
668 Vector &
669 operator=(const ::Vector<Number> &v);
670
688 void
691 const Vector & vector);
692
699 void
701 const VectorOperation::values operation);
702
703
709 bool
710 operator==(const Vector &v) const;
711
717 bool
718 operator!=(const Vector &v) const;
719
724 size() const;
725
741 local_size() const;
742
749
771 std::pair<size_type, size_type>
772 local_range() const;
773
781 bool
782 in_local_range(const size_type index) const;
783
799
807 bool
809
816 void
818
824 operator*(const Vector &vec) const;
825
830 norm_sqr() const;
831
836 mean_value() const;
837
842 min() const;
843
848 max() const;
849
854 l1_norm() const;
855
861 l2_norm() const;
862
868 lp_norm(const TrilinosScalar p) const;
869
874 linfty_norm() const;
875
896 add_and_dot(const TrilinosScalar a, const Vector &V, const Vector &W);
897
903 bool
904 all_zero() const;
905
911 bool
912 is_non_negative() const;
914
915
920
929 operator()(const size_type index);
930
939 operator()(const size_type index) const;
940
947 operator[](const size_type index);
948
955 operator[](const size_type index) const;
956
972 void
973 extract_subvector_to(const std::vector<size_type> &indices,
974 std::vector<TrilinosScalar> & values) const;
975
1003 template <typename ForwardIterator, typename OutputIterator>
1004 void
1005 extract_subvector_to(ForwardIterator indices_begin,
1006 const ForwardIterator indices_end,
1007 OutputIterator values_begin) const;
1008
1019 iterator
1021
1027 begin() const;
1028
1033 iterator
1035
1041 end() const;
1042
1044
1045
1050
1057 void
1058 set(const std::vector<size_type> & indices,
1059 const std::vector<TrilinosScalar> &values);
1060
1065 void
1066 set(const std::vector<size_type> & indices,
1067 const ::Vector<TrilinosScalar> &values);
1068
1074 void
1075 set(const size_type n_elements,
1076 const size_type * indices,
1077 const TrilinosScalar *values);
1078
1083 void
1084 add(const std::vector<size_type> & indices,
1085 const std::vector<TrilinosScalar> &values);
1086
1091 void
1092 add(const std::vector<size_type> & indices,
1093 const ::Vector<TrilinosScalar> &values);
1094
1100 void
1101 add(const size_type n_elements,
1102 const size_type * indices,
1103 const TrilinosScalar *values);
1104
1108 Vector &
1110
1114 Vector &
1116
1120 Vector &
1122
1126 Vector &
1128
1133 void
1135
1148 void
1149 add(const Vector &V, const bool allow_different_maps = false);
1150
1154 void
1155 add(const TrilinosScalar a, const Vector &V);
1156
1160 void
1162 const Vector & V,
1163 const TrilinosScalar b,
1164 const Vector & W);
1165
1170 void
1171 sadd(const TrilinosScalar s, const Vector &V);
1172
1176 void
1177 sadd(const TrilinosScalar s, const TrilinosScalar a, const Vector &V);
1178
1184 void
1185 scale(const Vector &scaling_factors);
1186
1190 void
1191 equ(const TrilinosScalar a, const Vector &V);
1193
1198
1203 const Epetra_MultiVector &
1205
1210 Epetra_FEVector &
1212
1217 const Epetra_BlockMap &
1219
1227 void
1228 print(std::ostream & out,
1229 const unsigned int precision = 3,
1230 const bool scientific = true,
1231 const bool across = true) const;
1232
1246 void
1247 swap(Vector &v);
1248
1252 std::size_t
1253 memory_consumption() const;
1254
1259 const MPI_Comm &
1262
1267
1272 int,
1273 << "An error with error number " << arg1
1274 << " occurred while calling a Trilinos function");
1275
1281 size_type,
1282 size_type,
1283 size_type,
1284 size_type,
1285 << "You are trying to access element " << arg1
1286 << " of a distributed vector, but this element is not stored "
1287 << "on the current processor. Note: There are " << arg2
1288 << " elements stored "
1289 << "on the current processor from within the range [" << arg3 << ','
1290 << arg4 << "] but Trilinos vectors need not store contiguous "
1291 << "ranges on each processor, and not every element in "
1292 << "this range may in fact be stored locally."
1293 << "\n\n"
1294 << "A common source for this kind of problem is that you "
1295 << "are passing a 'fully distributed' vector into a function "
1296 << "that needs read access to vector elements that correspond "
1297 << "to degrees of freedom on ghost cells (or at least to "
1298 << "'locally active' degrees of freedom that are not also "
1299 << "'locally owned'). You need to pass a vector that has these "
1300 << "elements as ghost entries.");
1301
1302 private:
1314 Epetra_CombineMode last_action;
1315
1321
1327
1333 std::unique_ptr<Epetra_FEVector> vector;
1334
1340 std::unique_ptr<Epetra_MultiVector> nonlocal_vector;
1341
1346
1347 // Make the reference class a friend.
1349 };
1350
1351
1352
1353 // ------------------- inline and template functions --------------
1354
1355
1363 inline void
1365 {
1366 u.swap(v);
1367 }
1368 } // namespace MPI
1369
1370# ifndef DOXYGEN
1371
1372 namespace internal
1373 {
1374 inline VectorReference::VectorReference(MPI::Vector & vector,
1375 const size_type index)
1376 : vector(vector)
1377 , index(index)
1378 {}
1379
1380
1381 inline const VectorReference &
1382 VectorReference::operator=(const VectorReference &r) const
1383 {
1384 // as explained in the class
1385 // documentation, this is not the copy
1386 // operator. so simply pass on to the
1387 // "correct" assignment operator
1388 *this = static_cast<TrilinosScalar>(r);
1389
1390 return *this;
1391 }
1392
1393
1394
1395 inline VectorReference &
1396 VectorReference::operator=(const VectorReference &r)
1397 {
1398 // as above
1399 *this = static_cast<TrilinosScalar>(r);
1400
1401 return *this;
1402 }
1403
1404
1405 inline const VectorReference &
1406 VectorReference::operator=(const TrilinosScalar &value) const
1407 {
1408 vector.set(1, &index, &value);
1409 return *this;
1410 }
1411
1412
1413
1414 inline const VectorReference &
1415 VectorReference::operator+=(const TrilinosScalar &value) const
1416 {
1417 vector.add(1, &index, &value);
1418 return *this;
1419 }
1420
1421
1422
1423 inline const VectorReference &
1424 VectorReference::operator-=(const TrilinosScalar &value) const
1425 {
1426 TrilinosScalar new_value = -value;
1427 vector.add(1, &index, &new_value);
1428 return *this;
1429 }
1430
1431
1432
1433 inline const VectorReference &
1434 VectorReference::operator*=(const TrilinosScalar &value) const
1435 {
1436 TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
1437 vector.set(1, &index, &new_value);
1438 return *this;
1439 }
1440
1441
1442
1443 inline const VectorReference &
1444 VectorReference::operator/=(const TrilinosScalar &value) const
1445 {
1446 TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
1447 vector.set(1, &index, &new_value);
1448 return *this;
1449 }
1450 } // namespace internal
1451
1452 namespace MPI
1453 {
1454 inline bool
1455 Vector::in_local_range(const size_type index) const
1456 {
1457 std::pair<size_type, size_type> range = local_range();
1458
1459 return ((index >= range.first) && (index < range.second));
1460 }
1461
1462
1463
1464 inline IndexSet
1466 {
1468 ExcMessage(
1469 "The locally owned elements have not been properly initialized!"
1470 " This happens for example if this object has been initialized"
1471 " with exactly one overlapping IndexSet."));
1472 return owned_elements;
1473 }
1474
1475
1476
1477 inline bool
1479 {
1480 return has_ghosts;
1481 }
1482
1483
1484
1485 inline void
1487 {}
1488
1489
1490
1491 inline internal::VectorReference
1492 Vector::operator()(const size_type index)
1493 {
1494 return internal::VectorReference(*this, index);
1495 }
1496
1497
1498
1499 inline internal::VectorReference
1500 Vector::operator[](const size_type index)
1501 {
1502 return operator()(index);
1503 }
1504
1505
1506
1507 inline TrilinosScalar
1508 Vector::operator[](const size_type index) const
1509 {
1510 return operator()(index);
1511 }
1512
1513
1514
1515 inline void
1516 Vector::extract_subvector_to(const std::vector<size_type> &indices,
1517 std::vector<TrilinosScalar> & values) const
1518 {
1519 for (size_type i = 0; i < indices.size(); ++i)
1520 values[i] = operator()(indices[i]);
1521 }
1522
1523
1524
1525 template <typename ForwardIterator, typename OutputIterator>
1526 inline void
1527 Vector::extract_subvector_to(ForwardIterator indices_begin,
1528 const ForwardIterator indices_end,
1529 OutputIterator values_begin) const
1530 {
1531 while (indices_begin != indices_end)
1532 {
1533 *values_begin = operator()(*indices_begin);
1534 indices_begin++;
1535 values_begin++;
1536 }
1537 }
1538
1539
1540
1541 inline Vector::iterator
1543 {
1544 return (*vector)[0];
1545 }
1546
1547
1548
1549 inline Vector::iterator
1550 Vector::end()
1551 {
1552 return (*vector)[0] + locally_owned_size();
1553 }
1554
1555
1556
1558 Vector::begin() const
1559 {
1560 return (*vector)[0];
1561 }
1562
1563
1564
1566 Vector::end() const
1567 {
1568 return (*vector)[0] + locally_owned_size();
1569 }
1570
1571
1572
1573 inline void
1574 Vector::set(const std::vector<size_type> & indices,
1575 const std::vector<TrilinosScalar> &values)
1576 {
1577 // if we have ghost values, do not allow
1578 // writing to this vector at all.
1580
1581 AssertDimension(indices.size(), values.size());
1582
1583 set(indices.size(), indices.data(), values.data());
1584 }
1585
1586
1587
1588 inline void
1589 Vector::set(const std::vector<size_type> & indices,
1590 const ::Vector<TrilinosScalar> &values)
1591 {
1592 // if we have ghost values, do not allow
1593 // writing to this vector at all.
1595
1596 AssertDimension(indices.size(), values.size());
1597
1598 set(indices.size(), indices.data(), values.begin());
1599 }
1600
1601
1602
1603 inline void
1604 Vector::set(const size_type n_elements,
1605 const size_type * indices,
1606 const TrilinosScalar *values)
1607 {
1608 // if we have ghost values, do not allow
1609 // writing to this vector at all.
1611
1612 if (last_action == Add)
1613 {
1614 const int ierr = vector->GlobalAssemble(Add);
1615 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1616 }
1617
1618 if (last_action != Insert)
1619 last_action = Insert;
1620
1621 for (size_type i = 0; i < n_elements; ++i)
1622 {
1623 const TrilinosWrappers::types::int_type row = indices[i];
1624 const TrilinosWrappers::types::int_type local_row =
1625 vector->Map().LID(row);
1626 if (local_row != -1)
1627 (*vector)[0][local_row] = values[i];
1628 else
1629 {
1630 const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1631 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1632 compressed = false;
1633 }
1634 // in set operation, do not use the pre-allocated vector for nonlocal
1635 // entries even if it exists. This is to ensure that we really only
1636 // set the elements touched by the set() method and not all contained
1637 // in the nonlocal entries vector (there is no way to distinguish them
1638 // on the receiving processor)
1639 }
1640 }
1641
1642
1643
1644 inline void
1645 Vector::add(const std::vector<size_type> & indices,
1646 const std::vector<TrilinosScalar> &values)
1647 {
1648 // if we have ghost values, do not allow
1649 // writing to this vector at all.
1651 AssertDimension(indices.size(), values.size());
1652
1653 add(indices.size(), indices.data(), values.data());
1654 }
1655
1656
1657
1658 inline void
1659 Vector::add(const std::vector<size_type> & indices,
1660 const ::Vector<TrilinosScalar> &values)
1661 {
1662 // if we have ghost values, do not allow
1663 // writing to this vector at all.
1665 AssertDimension(indices.size(), values.size());
1666
1667 add(indices.size(), indices.data(), values.begin());
1668 }
1669
1670
1671
1672 inline void
1673 Vector::add(const size_type n_elements,
1674 const size_type * indices,
1675 const TrilinosScalar *values)
1676 {
1677 // if we have ghost values, do not allow
1678 // writing to this vector at all.
1680
1681 if (last_action != Add)
1682 {
1683 if (last_action == Insert)
1684 {
1685 const int ierr = vector->GlobalAssemble(Insert);
1686 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1687 }
1688 last_action = Add;
1689 }
1690
1691 for (size_type i = 0; i < n_elements; ++i)
1692 {
1693 const size_type row = indices[i];
1694 const TrilinosWrappers::types::int_type local_row = vector->Map().LID(
1695 static_cast<TrilinosWrappers::types::int_type>(row));
1696 if (local_row != -1)
1697 (*vector)[0][local_row] += values[i];
1698 else if (nonlocal_vector.get() == nullptr)
1699 {
1700 const int ierr = vector->SumIntoGlobalValues(
1701 1,
1702 reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1703 &row),
1704 &values[i]);
1705 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1706 compressed = false;
1707 }
1708 else
1709 {
1710 // use pre-allocated vector for non-local entries if it exists for
1711 // addition operation
1713 nonlocal_vector->Map().LID(
1714 static_cast<TrilinosWrappers::types::int_type>(row));
1715 Assert(my_row != -1,
1716 ExcMessage(
1717 "Attempted to write into off-processor vector entry "
1718 "that has not be specified as being writable upon "
1719 "initialization"));
1720 (*nonlocal_vector)[0][my_row] += values[i];
1721 compressed = false;
1722 }
1723 }
1724 }
1725
1726
1727
1728 inline Vector::size_type
1729 Vector::size() const
1730 {
1731# ifndef DEAL_II_WITH_64BIT_INDICES
1732 return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1733# else
1734 return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1735# endif
1736 }
1737
1738
1739
1740 inline Vector::size_type
1741 Vector::local_size() const
1742 {
1743 return vector->Map().NumMyElements();
1744 }
1745
1746
1747
1748 inline Vector::size_type
1750 {
1751 return owned_elements.n_elements();
1752 }
1753
1754
1755
1756 inline std::pair<Vector::size_type, Vector::size_type>
1757 Vector::local_range() const
1758 {
1759# ifndef DEAL_II_WITH_64BIT_INDICES
1760 const TrilinosWrappers::types::int_type begin = vector->Map().MinMyGID();
1762 vector->Map().MaxMyGID() + 1;
1763# else
1765 vector->Map().MinMyGID64();
1767 vector->Map().MaxMyGID64() + 1;
1768# endif
1769
1770 Assert(
1771 end - begin == vector->Map().NumMyElements(),
1772 ExcMessage(
1773 "This function only makes sense if the elements that this "
1774 "vector stores on the current processor form a contiguous range. "
1775 "This does not appear to be the case for the current vector."));
1776
1777 return std::make_pair(begin, end);
1778 }
1779
1780
1781
1782 inline TrilinosScalar
1783 Vector::operator*(const Vector &vec) const
1784 {
1785 Assert(vector->Map().SameAs(vec.vector->Map()),
1788
1789 TrilinosScalar result;
1790
1791 const int ierr = vector->Dot(*(vec.vector), &result);
1792 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1793
1794 return result;
1795 }
1796
1797
1798
1799 inline Vector::real_type
1800 Vector::norm_sqr() const
1801 {
1802 const TrilinosScalar d = l2_norm();
1803 return d * d;
1804 }
1805
1806
1807
1808 inline TrilinosScalar
1809 Vector::mean_value() const
1810 {
1812
1814 const int ierr = vector->MeanValue(&mean);
1815 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1816
1817 return mean;
1818 }
1819
1820
1821
1822 inline TrilinosScalar
1823 Vector::min() const
1824 {
1825 TrilinosScalar min_value;
1826 const int ierr = vector->MinValue(&min_value);
1827 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1828
1829 return min_value;
1830 }
1831
1832
1833
1834 inline TrilinosScalar
1835 Vector::max() const
1836 {
1837 TrilinosScalar max_value;
1838 const int ierr = vector->MaxValue(&max_value);
1839 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1840
1841 return max_value;
1842 }
1843
1844
1845
1846 inline Vector::real_type
1847 Vector::l1_norm() const
1848 {
1850
1852 const int ierr = vector->Norm1(&d);
1853 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1854
1855 return d;
1856 }
1857
1858
1859
1860 inline Vector::real_type
1861 Vector::l2_norm() const
1862 {
1864
1866 const int ierr = vector->Norm2(&d);
1867 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1868
1869 return d;
1870 }
1871
1872
1873
1874 inline Vector::real_type
1875 Vector::lp_norm(const TrilinosScalar p) const
1876 {
1878
1879 TrilinosScalar norm = 0;
1880 TrilinosScalar sum = 0;
1881 const size_type n_local = locally_owned_size();
1882
1883 // loop over all the elements because
1884 // Trilinos does not support lp norms
1885 for (size_type i = 0; i < n_local; ++i)
1886 sum += std::pow(std::fabs((*vector)[0][i]), p);
1887
1888 norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1889
1890 return norm;
1891 }
1892
1893
1894
1895 inline Vector::real_type
1896 Vector::linfty_norm() const
1897 {
1898 // while we disallow the other
1899 // norm operations on ghosted
1900 // vectors, this particular norm
1901 // is safe to run even in the
1902 // presence of ghost elements
1904 const int ierr = vector->NormInf(&d);
1905 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1906
1907 return d;
1908 }
1909
1910
1911
1912 inline TrilinosScalar
1914 const Vector & V,
1915 const Vector & W)
1916 {
1917 this->add(a, V);
1918 return *this * W;
1919 }
1920
1921
1922
1923 // inline also scalar products, vector
1924 // additions etc. since they are all
1925 // representable by a single Trilinos
1926 // call. This reduces the overhead of the
1927 // wrapper class.
1928 inline Vector &
1930 {
1931 AssertIsFinite(a);
1932
1933 const int ierr = vector->Scale(a);
1934 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1935
1936 return *this;
1937 }
1938
1939
1940
1941 inline Vector &
1943 {
1944 AssertIsFinite(a);
1945
1946 const TrilinosScalar factor = 1. / a;
1947
1948 AssertIsFinite(factor);
1949
1950 const int ierr = vector->Scale(factor);
1951 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1952
1953 return *this;
1954 }
1955
1956
1957
1958 inline Vector &
1959 Vector::operator+=(const Vector &v)
1960 {
1961 AssertDimension(size(), v.size());
1962 Assert(vector->Map().SameAs(v.vector->Map()),
1964
1965 const int ierr = vector->Update(1.0, *(v.vector), 1.0);
1966 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1967
1968 return *this;
1969 }
1970
1971
1972
1973 inline Vector &
1974 Vector::operator-=(const Vector &v)
1975 {
1976 AssertDimension(size(), v.size());
1977 Assert(vector->Map().SameAs(v.vector->Map()),
1979
1980 const int ierr = vector->Update(-1.0, *(v.vector), 1.0);
1981 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
1982
1983 return *this;
1984 }
1985
1986
1987
1988 inline void
1990 {
1991 // if we have ghost values, do not allow
1992 // writing to this vector at all.
1994 AssertIsFinite(s);
1995
1996 size_type n_local = locally_owned_size();
1997 for (size_type i = 0; i < n_local; ++i)
1998 (*vector)[0][i] += s;
1999 }
2000
2001
2002
2003 inline void
2004 Vector::add(const TrilinosScalar a, const Vector &v)
2005 {
2006 // if we have ghost values, do not allow
2007 // writing to this vector at all.
2010
2011 AssertIsFinite(a);
2012
2013 const int ierr = vector->Update(a, *(v.vector), 1.);
2014 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2015 }
2016
2017
2018
2019 inline void
2021 const Vector & v,
2022 const TrilinosScalar b,
2023 const Vector & w)
2024 {
2025 // if we have ghost values, do not allow
2026 // writing to this vector at all.
2029 AssertDimension(locally_owned_size(), w.locally_owned_size());
2030
2031 AssertIsFinite(a);
2033
2034 const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
2035
2036 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2037 }
2038
2039
2040
2041 inline void
2042 Vector::sadd(const TrilinosScalar s, const Vector &v)
2043 {
2044 // if we have ghost values, do not allow
2045 // writing to this vector at all.
2047 AssertDimension(size(), v.size());
2048
2049 AssertIsFinite(s);
2050
2051 // We assume that the vectors have the same Map
2052 // if the local size is the same and if the vectors are not ghosted
2054 !v.has_ghost_elements())
2055 {
2056 Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2058 const int ierr = vector->Update(1., *(v.vector), s);
2059 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2060 }
2061 else
2062 {
2063 (*this) *= s;
2064 this->add(v, true);
2065 }
2066 }
2067
2068
2069
2070 inline void
2072 const TrilinosScalar a,
2073 const Vector & v)
2074 {
2075 // if we have ghost values, do not allow
2076 // writing to this vector at all.
2078 AssertDimension(size(), v.size());
2079 AssertIsFinite(s);
2080 AssertIsFinite(a);
2081
2082 // We assume that the vectors have the same Map
2083 // if the local size is the same and if the vectors are not ghosted
2085 !v.has_ghost_elements())
2086 {
2087 Assert(this->vector->Map().SameAs(v.vector->Map()) == true,
2089 const int ierr = vector->Update(a, *(v.vector), s);
2090 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2091 }
2092 else
2093 {
2094 (*this) *= s;
2095 Vector tmp = v;
2096 tmp *= a;
2097 this->add(tmp, true);
2098 }
2099 }
2100
2101
2102
2103 inline void
2104 Vector::scale(const Vector &factors)
2105 {
2106 // if we have ghost values, do not allow
2107 // writing to this vector at all.
2110
2111 const int ierr = vector->Multiply(1.0, *(factors.vector), *vector, 0.0);
2112 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2113 }
2114
2115
2116
2117 inline void
2118 Vector::equ(const TrilinosScalar a, const Vector &v)
2119 {
2120 // if we have ghost values, do not allow
2121 // writing to this vector at all.
2123 AssertIsFinite(a);
2124
2125 // If we don't have the same map, copy.
2126 if (vector->Map().SameAs(v.vector->Map()) == false)
2127 {
2128 this->sadd(0., a, v);
2129 }
2130 else
2131 {
2132 // Otherwise, just update
2133 int ierr = vector->Update(a, *v.vector, 0.0);
2134 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2135
2136 last_action = Zero;
2137 }
2138 }
2139
2140
2141
2142 inline const Epetra_MultiVector &
2144 {
2145 return static_cast<const Epetra_MultiVector &>(*vector);
2146 }
2147
2148
2149
2150 inline Epetra_FEVector &
2152 {
2153 return *vector;
2154 }
2155
2156
2157
2158 inline const Epetra_BlockMap &
2160 {
2161 return vector->Map();
2162 }
2163
2164
2165
2166 inline const MPI_Comm &
2168 {
2169 static MPI_Comm comm;
2170
2171 const Epetra_MpiComm *mpi_comm =
2172 dynamic_cast<const Epetra_MpiComm *>(&vector->Map().Comm());
2173 comm = mpi_comm->Comm();
2174
2175 return comm;
2176 }
2177
2178 template <typename number>
2179 Vector::Vector(const IndexSet & parallel_partitioner,
2180 const ::Vector<number> &v,
2181 const MPI_Comm & communicator)
2182 {
2183 *this =
2184 Vector(parallel_partitioner.make_trilinos_map(communicator, true), v);
2185 owned_elements = parallel_partitioner;
2186 }
2187
2188
2189
2190 inline Vector &
2192 {
2193 AssertIsFinite(s);
2194
2195 int ierr = vector->PutScalar(s);
2196 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2197
2198 if (nonlocal_vector.get() != nullptr)
2199 {
2200 ierr = nonlocal_vector->PutScalar(0.);
2201 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2202 }
2203
2204 return *this;
2205 }
2206 } /* end of namespace MPI */
2207
2208# endif /* DOXYGEN */
2209
2210} /* end of namespace TrilinosWrappers */
2211
2215namespace internal
2216{
2217 namespace LinearOperatorImplementation
2218 {
2219 template <typename>
2220 class ReinitHelper;
2221
2226 template <>
2228 {
2229 public:
2230 template <typename Matrix>
2231 static void
2234 bool omit_zeroing_entries)
2235 {
2236 v.reinit(matrix.locally_owned_range_indices(),
2237 matrix.get_mpi_communicator(),
2238 omit_zeroing_entries);
2239 }
2240
2241 template <typename Matrix>
2242 static void
2245 bool omit_zeroing_entries)
2246 {
2247 v.reinit(matrix.locally_owned_domain_indices(),
2248 matrix.get_mpi_communicator(),
2249 omit_zeroing_entries);
2250 }
2251 };
2252
2253 } // namespace LinearOperatorImplementation
2254} /* namespace internal */
2255
2256
2257
2261template <>
2263{};
2264
2265
2267
2268#endif // DEAL_II_WITH_TRILINOS
2269
2270/*---------------------------- trilinos_vector.h ---------------------------*/
2271
2272#endif // dealii_trilinos_vector_h
size_type size() const
Definition: index_set.h:1636
size_type n_elements() const
Definition: index_set.h:1834
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:601
Definition: vector.h:109
#define DEAL_II_DEPRECATED
Definition: config.h:164
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DeclException0(Exception0)
Definition: exceptions.h:464
static ::ExceptionBase & ExcGhostsPresent()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:578
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertIsFinite(number)
Definition: exceptions.h:1758
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
Vector & operator/=(const TrilinosScalar factor)
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)
const MPI_Comm & get_mpi_communicator() const
static void reinit_domain_vector(const Matrix &matrix, TrilinosWrappers::MPI::Vector &v, bool omit_zeroing_entries)
std::unique_ptr< Epetra_MultiVector > nonlocal_vector
void add(const TrilinosScalar s)
static ::ExceptionBase & ExcDifferentParallelPartitioning()
const Epetra_BlockMap & trilinos_partitioner() const
int gid(const Epetra_BlockMap &map, int i)
reference operator()(const size_type index)
void import_nonlocal_data_for_fe(const ::TrilinosWrappers::SparseMatrix &matrix, const Vector &vector)
std::unique_ptr< Epetra_FEVector > vector
bool in_local_range(const size_type index) const
size_type local_size() const
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
friend class internal::VectorReference
real_type lp_norm(const TrilinosScalar p) const
IndexSet locally_owned_elements() const
Vector & operator-=(const Vector &V)
Vector & operator+=(const Vector &V)
const Epetra_MultiVector & trilinos_vector() const
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
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
void compress(::VectorOperation::values operation)
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)
static void reinit_range_vector(const Matrix &matrix, TrilinosWrappers::MPI::Vector &v, bool omit_zeroing_entries)
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)
static ::ExceptionBase & ExcTrilinosError(int arg1)
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
Vector(const IndexSet &parallel_partitioning, const ::Vector< Number > &v, const MPI_Comm &communicator=MPI_COMM_WORLD)
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:139
bool has_ghost_elements() const
const value_type * const_iterator
Definition: vector.h:125
size_type size() const
value_type * iterator
Definition: vector.h:124
size_type locally_owned_size() const
Expression fabs(const Expression &x)
@ matrix
Contents is actually a matrix.
static const char V
SparseMatrix< double > SparseMatrix
types::global_dof_index size_type
Definition: cuda_kernels.h:45
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition: divergence.h:472
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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:76
const MPI_Comm & comm
double TrilinosScalar
Definition: types.h:163