16 #ifndef dealii_trilinos_vector_h 17 #define dealii_trilinos_vector_h 22 #ifdef DEAL_II_WITH_TRILINOS 33 # include <Epetra_ConfigDefs.h> 38 # ifdef DEAL_II_WITH_MPI // only if MPI is installed 39 # include <Epetra_MpiComm.h> 42 # include <Epetra_SerialComm.h> 44 # include <Epetra_FEVector.h> 45 # include <Epetra_LocalMap.h> 46 # include <Epetra_Map.h> 55 template <
typename Number>
56 class ReadWriteVector;
100 class VectorReference
107 VectorReference(
MPI::Vector &vector,
const size_type index);
113 VectorReference(
const VectorReference &) =
default;
126 const VectorReference &
127 operator=(
const VectorReference &r)
const;
133 operator=(
const VectorReference &r);
138 const VectorReference &
144 const VectorReference &
150 const VectorReference &
156 const VectorReference &
162 const VectorReference &
176 <<
"An error with error number " << arg1
177 <<
" occurred while calling a Trilinos function");
188 const size_type index;
199 # ifndef DEAL_II_WITH_64BIT_INDICES 204 gid(
const Epetra_BlockMap &map,
int i)
213 gid(
const Epetra_BlockMap &map,
int i)
454 const MPI_Comm &communicator = MPI_COMM_WORLD);
469 const MPI_Comm &communicator = MPI_COMM_WORLD);
487 const MPI_Comm &communicator = MPI_COMM_WORLD);
501 template <
typename Number>
503 const ::Vector<Number> &v,
504 const MPI_Comm & communicator = MPI_COMM_WORLD);
515 ~
Vector()
override =
default;
549 const bool omit_zeroing_entries =
false,
550 const bool allow_different_maps =
false);
576 const MPI_Comm &communicator = MPI_COMM_WORLD,
577 const bool omit_zeroing_entries =
false);
607 const MPI_Comm &communicator = MPI_COMM_WORLD,
608 const bool vector_writable =
false);
656 operator=(
const Vector &v);
663 operator=(
Vector &&v) noexcept;
672 template <
typename Number>
674 operator=(const ::Vector<Number> &v);
694 import_nonlocal_data_for_fe(
753 locally_owned_size()
const;
776 std::pair<size_type, size_type>
787 in_local_range(
const size_type index)
const;
803 locally_owned_elements()
const;
813 has_ghost_elements()
const;
822 update_ghost_values()
const;
975 extract_subvector_to(
const std::vector<size_type> &indices,
976 std::vector<TrilinosScalar> &
values)
const;
1005 template <
typename ForwardIterator,
typename OutputIterator>
1007 extract_subvector_to(ForwardIterator indices_begin,
1008 const ForwardIterator indices_end,
1009 OutputIterator values_begin)
const;
1060 set(
const std::vector<size_type> & indices,
1061 const std::vector<TrilinosScalar> &
values);
1068 set(
const std::vector<size_type> & indices,
1069 const ::Vector<TrilinosScalar> &
values);
1086 add(
const std::vector<size_type> & indices,
1087 const std::vector<TrilinosScalar> &
values);
1094 add(
const std::vector<size_type> & indices,
1095 const ::Vector<TrilinosScalar> &
values);
1151 add(
const Vector &
V,
const bool allow_different_maps =
false);
1205 const Epetra_MultiVector &
1206 trilinos_vector()
const;
1219 const Epetra_BlockMap &
1220 trilinos_partitioner()
const;
1230 print(std::ostream & out,
1231 const unsigned int precision = 3,
1232 const bool scientific =
true,
1233 const bool across =
true)
const;
1262 get_mpi_communicator()
const;
1275 <<
"An error with error number " << arg1
1276 <<
" occurred while calling a Trilinos function");
1282 ExcAccessToNonLocalElement,
1287 <<
"You are trying to access element " << arg1
1288 <<
" of a distributed vector, but this element is not stored " 1289 <<
"on the current processor. Note: There are " << arg2
1290 <<
" elements stored " 1291 <<
"on the current processor from within the range [" << arg3 <<
"," 1292 << arg4 <<
"] but Trilinos vectors need not store contiguous " 1293 <<
"ranges on each processor, and not every element in " 1294 <<
"this range may in fact be stored locally." 1296 <<
"A common source for this kind of problem is that you " 1297 <<
"are passing a 'fully distributed' vector into a function " 1298 <<
"that needs read access to vector elements that correspond " 1299 <<
"to degrees of freedom on ghost cells (or at least to " 1300 <<
"'locally active' degrees of freedom that are not also " 1301 <<
"'locally owned'). You need to pass a vector that has these " 1302 <<
"elements as ghost entries.");
1350 friend class internal::VectorReference;
1376 inline VectorReference::VectorReference(
MPI::Vector & vector,
1383 inline const VectorReference &
1384 VectorReference::operator=(
const VectorReference &r)
const 1397 inline VectorReference &
1398 VectorReference::operator=(
const VectorReference &r)
1407 inline const VectorReference &
1410 vector.
set(1, &index, &value);
1416 inline const VectorReference &
1419 vector.
add(1, &index, &value);
1425 inline const VectorReference &
1429 vector.
add(1, &index, &new_value);
1435 inline const VectorReference &
1439 vector.
set(1, &index, &new_value);
1445 inline const VectorReference &
1449 vector.
set(1, &index, &new_value);
1459 std::pair<size_type, size_type> range = local_range();
1461 return ((index >= range.first) && (index < range.second));
1469 Assert(owned_elements.size() == size(),
1471 "The locally owned elements have not been properly initialized!" 1472 " This happens for example if this object has been initialized" 1473 " with exactly one overlapping IndexSet."));
1474 return owned_elements;
1493 inline internal::VectorReference
1496 return internal::VectorReference(*
this, index);
1503 return operator()(index);
1510 return operator()(index);
1517 std::vector<TrilinosScalar> &
values)
const 1519 for (
size_type i = 0; i < indices.size(); ++i)
1520 values[i] =
operator()(indices[i]);
1525 template <
typename ForwardIterator,
typename OutputIterator>
1528 const ForwardIterator indices_end,
1529 OutputIterator values_begin)
const 1531 while (indices_begin != indices_end)
1533 *values_begin = operator()(*indices_begin);
1544 return (*vector)[0];
1552 return (*vector)[0] + local_size();
1560 return (*vector)[0];
1568 return (*vector)[0] + local_size();
1574 Vector::set(
const std::vector<size_type> & indices,
1575 const std::vector<TrilinosScalar> &values)
1581 Assert(indices.size() == values.size(),
1584 set(indices.size(), indices.data(), values.data());
1590 Vector::set(
const std::vector<size_type> & indices,
1591 const ::Vector<TrilinosScalar> &values)
1597 Assert(indices.size() == values.size(),
1600 set(indices.size(), indices.data(), values.begin());
1614 if (last_action == Add)
1616 const int ierr = vector->GlobalAssemble(Add);
1620 if (last_action != Insert)
1621 last_action = Insert;
1623 for (
size_type i = 0; i < n_elements; ++i)
1627 vector->Map().LID(row);
1628 if (local_row != -1)
1629 (*vector)[0][local_row] = values[i];
1632 const int ierr = vector->ReplaceGlobalValues(1, &row, &values[i]);
1647 Vector::add(
const std::vector<size_type> & indices,
1648 const std::vector<TrilinosScalar> &values)
1653 Assert(indices.size() == values.size(),
1656 add(indices.size(), indices.data(), values.data());
1662 Vector::add(
const std::vector<size_type> & indices,
1663 const ::Vector<TrilinosScalar> &values)
1668 Assert(indices.size() == values.size(),
1671 add(indices.size(), indices.data(), values.begin());
1685 if (last_action != Add)
1687 if (last_action == Insert)
1689 const int ierr = vector->GlobalAssemble(Insert);
1695 for (
size_type i = 0; i < n_elements; ++i)
1699 static_cast<TrilinosWrappers::types::int_type>(row));
1700 if (local_row != -1)
1701 (*vector)[0][local_row] += values[i];
1702 else if (nonlocal_vector.get() ==
nullptr)
1704 const int ierr = vector->SumIntoGlobalValues(
1706 reinterpret_cast<const TrilinosWrappers::types::int_type *>(
1717 nonlocal_vector->Map().LID(
1718 static_cast<TrilinosWrappers::types::int_type>(row));
1721 "Attempted to write into off-processor vector entry " 1722 "that has not be specified as being writable upon " 1724 (*nonlocal_vector)[0][my_row] += values[i];
1735 # ifndef DEAL_II_WITH_64BIT_INDICES 1736 return vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
1738 return vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64();
1745 Vector::local_size()
const 1747 return vector->Map().NumMyElements();
1755 return owned_elements.n_elements();
1760 inline std::pair<Vector::size_type, Vector::size_type>
1761 Vector::local_range()
const 1763 # ifndef DEAL_II_WITH_64BIT_INDICES 1766 vector->Map().MaxMyGID() + 1;
1769 vector->Map().MinMyGID64();
1771 vector->Map().MaxMyGID64() + 1;
1775 end - begin == vector->Map().NumMyElements(),
1777 "This function only makes sense if the elements that this " 1778 "vector stores on the current processor form a contiguous range. " 1779 "This does not appear to be the case for the current vector."));
1781 return std::make_pair(begin, end);
1789 ExcDifferentParallelPartitioning());
1794 const int ierr = vector->Dot(*(vec.
vector), &result);
1817 const int ierr = vector->MeanValue(&mean);
1829 const int ierr = vector->MinValue(&min_value);
1841 const int ierr = vector->MaxValue(&max_value);
1855 const int ierr = vector->Norm1(&d);
1869 const int ierr = vector->Norm2(&d);
1889 sum += std::pow(
std::fabs((*vector)[0][i]), p);
1891 norm = std::pow(sum, static_cast<TrilinosScalar>(1. / p));
1907 const int ierr = vector->NormInf(&d);
1936 const int ierr = vector->Scale(a);
1953 const int ierr = vector->Scale(factor);
1966 ExcDifferentParallelPartitioning());
1968 const int ierr = vector->Update(1.0, *(v.
vector), 1.0);
1981 ExcDifferentParallelPartitioning());
1983 const int ierr = vector->Update(-1.0, *(v.
vector), 1.0);
2001 (*vector)[0][i] += s;
2017 const int ierr = vector->Update(a, *(v.
vector), 1.);
2040 const int ierr = vector->Update(a, *(v.
vector), b, *(w.
vector), 1.);
2061 Assert(this->vector->Map().SameAs(v.
vector->Map()) ==
true,
2062 ExcDifferentParallelPartitioning());
2063 const int ierr = vector->Update(1., *(v.
vector), s);
2091 Assert(this->vector->Map().SameAs(v.
vector->Map()) ==
true,
2092 ExcDifferentParallelPartitioning());
2093 const int ierr = vector->Update(a, *(v.
vector), s);
2101 this->add(tmp,
true);
2116 const int ierr = vector->Multiply(1.0, *(factors.
vector), *vector, 0.0);
2131 if (vector->Map().SameAs(v.
vector->Map()) ==
false)
2133 this->
sadd(0., a, v);
2138 int ierr = vector->Update(a, *v.
vector, 0.0);
2147 inline const Epetra_MultiVector &
2148 Vector::trilinos_vector()
const 2150 return static_cast<const Epetra_MultiVector &
>(*vector);
2155 inline Epetra_FEVector &
2156 Vector::trilinos_vector()
2163 inline const Epetra_BlockMap &
2164 Vector::trilinos_partitioner()
const 2166 return vector->Map();
2172 Vector::get_mpi_communicator()
const 2176 # ifdef DEAL_II_WITH_MPI 2178 const Epetra_MpiComm *mpi_comm =
2179 dynamic_cast<const Epetra_MpiComm *
>(&vector->Map().Comm());
2180 comm = mpi_comm->Comm();
2184 comm = MPI_COMM_SELF;
2191 template <
typename number>
2193 const ::Vector<number> &v,
2198 owned_elements = parallel_partitioner;
2208 int ierr = vector->PutScalar(s);
2211 if (nonlocal_vector.get() !=
nullptr)
2213 ierr = nonlocal_vector->PutScalar(0.);
2230 namespace LinearOperatorImplementation
2243 template <
typename Matrix>
2247 bool omit_zeroing_entries)
2249 v.
reinit(matrix.locally_owned_range_indices(),
2250 matrix.get_mpi_communicator(),
2251 omit_zeroing_entries);
2254 template <
typename Matrix>
2258 bool omit_zeroing_entries)
2260 v.
reinit(matrix.locally_owned_domain_indices(),
2261 matrix.get_mpi_communicator(),
2262 omit_zeroing_entries);
2281 #endif // DEAL_II_WITH_TRILINOS 2285 #endif // dealii_trilinos_vector_h
std::unique_ptr< Epetra_MultiVector > nonlocal_vector
static ::ExceptionBase & ExcTrilinosError(int arg1)
void update_ghost_values() const
int gid(const Epetra_BlockMap &map, int i)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Number operator[](const size_type global_index) const
virtual VectorSpaceVector< Number >::real_type linfty_norm() const override
::types::global_dof_index size_type
Contents is actually a matrix.
#define DEAL_II_DEPRECATED_EARLY
types::global_dof_index size_type
static void reinit_range_vector(const Matrix &matrix, TrilinosWrappers::MPI::Vector &v, bool omit_zeroing_entries)
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
std::vector< value_type > l2_norm(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
void reinit(const Vector &v, const bool omit_zeroing_entries=false, const bool allow_different_maps=false)
virtual void add(const Number a) override
void add(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
__global__ void add_and_dot(Number *res, Number *v1, const Number *v2, const Number *v3, const Number a, const size_type N)
const internal::VectorReference const_reference
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Vector< Number > & operator=(const Vector< Number > &in_vector)
virtual Vector< Number > & operator*=(const Number factor) override
#define AssertThrow(cond, exc)
virtual Number operator*(const VectorSpaceVector< Number > &V) const override
virtual Vector< Number > & operator+=(const VectorSpaceVector< Number > &V) override
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
virtual void equ(const Number a, const VectorSpaceVector< Number > &V) override
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
virtual ::IndexSet locally_owned_elements() const override
static void reinit_domain_vector(const Matrix &matrix, TrilinosWrappers::MPI::Vector &v, bool omit_zeroing_entries)
virtual Vector< Number > & operator-=(const VectorSpaceVector< Number > &V) override
virtual void sadd(const Number s, const Number a, const VectorSpaceVector< Number > &V) override
static ::ExceptionBase & ExcMessage(std::string arg1)
std::string compress(const std::string &input)
bool in_local_range(const size_type global_index) const
void set(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
#define DeclException1(Exception1, type1, outsequence)
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::unique_ptr< Epetra_FEVector > vector
real_type lp_norm(const real_type p) const
size_type locally_owned_size() const
Number operator()(const size_type global_index) const
void swap(Vector &u, Vector &v)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
#define DeclException0(Exception0)
#define DEAL_II_NAMESPACE_CLOSE
virtual Vector< Number > & operator/=(const Number factor) override
virtual void scale(const VectorSpaceVector< Number > &scaling_factors) override
bool is_non_negative(const T &t)
VectorType::value_type * end(VectorType &V)
virtual size_type size() const override
Expression fabs(const Expression &x)
bool has_ghost_elements() const
__global__ void equ(Number *val, const Number a, const Number *V_val, const size_type N)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SparseMatrix< double > SparseMatrix
__global__ void sadd(const Number s, Number *val, const Number a, const Number *V_val, const size_type N)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
static ::ExceptionBase & ExcGhostsPresent()
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int global_dof_index
__global__ void set(Number *val, const Number s, const size_type N)
virtual VectorSpaceVector< Number >::real_type l1_norm() const override
size_type local_size() const
virtual VectorSpaceVector< Number >::real_type l2_norm() const override
*braid_SplitCommworld & comm
internal::VectorReference reference
#define DEAL_II_NAMESPACE_OPEN
VectorType::value_type * begin(VectorType &V)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Number l1_norm(const Tensor< 2, dim, Number > &t)
virtual value_type mean_value() const override
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< Number2 > &values) const
real_type norm_sqr() const
Epetra_CombineMode last_action
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
virtual Number add_and_dot(const Number a, const VectorSpaceVector< Number > &V, const VectorSpaceVector< Number > &W) override
bool has_ghost_elements() const
#define AssertIsFinite(number)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)