16 #ifndef dealii_petsc_vector_base_h
17 #define dealii_petsc_vector_base_h
22 #ifdef DEAL_II_WITH_PETSC
31 # include <petscvec.h>
40 template <
typename number>
93 VectorReference(
const VectorBase &vector,
const size_type index);
99 VectorReference(
const VectorReference &vector) =
default;
112 const VectorReference &
113 operator=(
const VectorReference &r)
const;
121 operator=(
const VectorReference &r);
126 const VectorReference &
127 operator=(
const PetscScalar &s)
const;
132 const VectorReference &
133 operator+=(
const PetscScalar &s)
const;
138 const VectorReference &
139 operator-=(
const PetscScalar &s)
const;
144 const VectorReference &
145 operator*=(
const PetscScalar &s)
const;
150 const VectorReference &
151 operator/=(
const PetscScalar &s)
const;
172 operator PetscScalar()
const;
177 ExcAccessToNonlocalElement,
181 <<
"You tried to access element " << arg1
182 <<
" of a distributed vector, but only elements in range [" << arg2
183 <<
',' << arg3 <<
"] are stored locally and can be accessed."
185 <<
"A common source for this kind of problem is that you "
186 <<
"are passing a 'fully distributed' vector into a function "
187 <<
"that needs read access to vector elements that correspond "
188 <<
"to degrees of freedom on ghost cells (or at least to "
189 <<
"'locally active' degrees of freedom that are not also "
190 <<
"'locally owned'). You need to pass a vector that has these "
191 <<
"elements as ghost entries.");
198 <<
"You tried to do a "
199 << (arg1 == 1 ?
"'set'" : (arg1 == 2 ?
"'add'" :
"???"))
200 <<
" operation but the vector is currently in "
201 << (arg2 == 1 ?
"'set'" : (arg2 == 2 ?
"'add'" :
"???"))
202 <<
" mode. You first have to call 'compress()'.");
208 const VectorBase &vector;
217 friend class ::PETScWrappers::VectorBase;
358 size()
const override;
379 std::pair<size_type, size_type>
461 set(
const std::vector<size_type> &indices,
462 const std::vector<PetscScalar> &
values);
481 std::vector<PetscScalar> &
values)
const;
518 template <
typename ForwardIterator,
typename OutputIterator>
521 const ForwardIterator indices_end,
522 OutputIterator values_begin)
const;
529 add(
const std::vector<size_type> &indices,
530 const std::vector<PetscScalar> &
values);
537 add(
const std::vector<size_type> &indices,
538 const ::Vector<PetscScalar> &
values);
548 const PetscScalar *
values);
661 add(
const PetscScalar s);
673 add(
const PetscScalar a,
712 write_ascii(
const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
722 print(std::ostream &out,
723 const unsigned int precision = 3,
724 const bool scientific =
true,
725 const bool across =
true)
const;
749 operator const Vec &()
const;
810 const PetscScalar *
values,
811 const bool add_values);
840 inline VectorReference::VectorReference(
const VectorBase &vector,
847 inline const VectorReference &
848 VectorReference::operator=(
const VectorReference &r)
const
854 *
this =
static_cast<PetscScalar
>(r);
861 inline VectorReference &
862 VectorReference::operator=(
const VectorReference &r)
868 *
this =
static_cast<PetscScalar
>(r);
875 inline const VectorReference &
876 VectorReference::operator=(
const PetscScalar &value)
const
884 const PetscInt petsc_i =
index;
886 const PetscErrorCode ierr =
887 VecSetValues(vector, 1, &petsc_i, &value, INSERT_VALUES);
897 inline const VectorReference &
898 VectorReference::operator+=(
const PetscScalar &value)
const
915 if (value == PetscScalar())
919 const PetscInt petsc_i =
index;
920 const PetscErrorCode ierr =
921 VecSetValues(vector, 1, &petsc_i, &value, ADD_VALUES);
930 inline const VectorReference &
931 VectorReference::operator-=(
const PetscScalar &value)
const
948 if (value == PetscScalar())
953 const PetscInt petsc_i =
index;
954 const PetscScalar subtractand = -
value;
955 const PetscErrorCode ierr =
956 VecSetValues(vector, 1, &petsc_i, &subtractand, ADD_VALUES);
964 inline const VectorReference &
965 VectorReference::operator*=(
const PetscScalar &value)
const
985 const PetscInt petsc_i =
index;
986 const PetscScalar new_value =
static_cast<PetscScalar
>(*this) *
value;
988 const PetscErrorCode ierr =
989 VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
997 inline const VectorReference &
998 VectorReference::operator/=(
const PetscScalar &value)
const
1018 const PetscInt petsc_i =
index;
1019 const PetscScalar new_value =
static_cast<PetscScalar
>(*this) /
value;
1021 const PetscErrorCode ierr =
1022 VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1031 VectorReference::real()
const
1033 # ifndef PETSC_USE_COMPLEX
1034 return static_cast<PetscScalar
>(*this);
1036 return PetscRealPart(
static_cast<PetscScalar
>(*
this));
1043 VectorReference::imag()
const
1045 # ifndef PETSC_USE_COMPLEX
1046 return PetscReal(0);
1048 return PetscImaginaryPart(
static_cast<PetscScalar
>(*
this));
1055 VectorBase::in_local_range(
const size_type index)
const
1058 const PetscErrorCode ierr =
1059 VecGetOwnershipRange(
static_cast<const Vec &
>(vector), &
begin, &
end);
1068 VectorBase::locally_owned_elements()
const
1073 const std::pair<size_type, size_type> x = local_range();
1074 is.add_range(x.first, x.second);
1081 VectorBase::has_ghost_elements()
const
1088 VectorBase::ghost_elements()
const
1090 return ghost_indices;
1095 VectorBase::update_ghost_values()
const
1099 PetscErrorCode ierr;
1101 ierr = VecGhostUpdateBegin(vector, INSERT_VALUES, SCATTER_FORWARD);
1103 ierr = VecGhostUpdateEnd(vector, INSERT_VALUES, SCATTER_FORWARD);
1110 inline internal::VectorReference
1111 VectorBase::operator()(
const size_type index)
1113 return internal::VectorReference(*
this, index);
1119 VectorBase::operator()(
const size_type index)
const
1121 return static_cast<PetscScalar
>(internal::VectorReference(*
this, index));
1126 inline internal::VectorReference
1127 VectorBase::operator[](
const size_type index)
1129 return operator()(index);
1135 VectorBase::operator[](
const size_type index)
const
1137 return operator()(index);
1141 VectorBase::get_mpi_communicator()
const
1143 return PetscObjectComm(
reinterpret_cast<PetscObject
>(vector));
1147 VectorBase::extract_subvector_to(
const std::vector<size_type> &indices,
1148 std::vector<PetscScalar> &
values)
const
1152 extract_subvector_to(indices.begin(), indices.end(),
values.begin());
1156 VectorBase::extract_subvector_to(
1161 extract_subvector_to(indices.
begin(), indices.
end(), elements.
begin());
1165 template <
typename ForwardIterator,
typename OutputIterator>
1167 VectorBase::extract_subvector_to(
const ForwardIterator indices_begin,
1168 const ForwardIterator indices_end,
1169 OutputIterator values_begin)
const
1171 const PetscInt n_idx =
static_cast<PetscInt
>(indices_end - indices_begin);
1197 PetscErrorCode ierr = VecGetOwnershipRange(vector, &
begin, &
end);
1200 Vec locally_stored_elements =
nullptr;
1201 ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
1205 ierr = VecGetSize(locally_stored_elements, &lsize);
1208 const PetscScalar *ptr;
1209 ierr = VecGetArrayRead(locally_stored_elements, &ptr);
1212 for (PetscInt i = 0; i < n_idx; ++i)
1214 const unsigned int index = *(indices_begin + i);
1215 if (index >=
static_cast<unsigned int>(
begin) &&
1216 index < static_cast<unsigned int>(
end))
1219 *(values_begin + i) = *(ptr + index -
begin);
1224 const unsigned int ghostidx =
1225 ghost_indices.index_within_set(index);
1228 *(values_begin + i) = *(ptr + ghostidx +
end -
begin);
1232 ierr = VecRestoreArrayRead(locally_stored_elements, &ptr);
1235 ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
1244 PetscErrorCode ierr = VecGetOwnershipRange(vector, &
begin, &
end);
1247 const PetscScalar *ptr;
1248 ierr = VecGetArrayRead(vector, &ptr);
1251 for (PetscInt i = 0; i < n_idx; ++i)
1253 const unsigned int index = *(indices_begin + i);
1255 Assert(index >=
static_cast<unsigned int>(
begin) &&
1256 index <
static_cast<unsigned int>(
end),
1257 ExcMessage(
"You are accessing elements of a vector without "
1258 "ghost elements that are not actually owned by "
1259 "this vector. A typical case where this may "
1260 "happen is if you are passing a non-ghosted "
1261 "(completely distributed) vector to a function "
1262 "that expects a vector that stores ghost "
1263 "elements for all locally relevant or locally "
1264 "active vector entries."));
1266 *(values_begin + i) = *(ptr + index -
begin);
1269 ierr = VecRestoreArrayRead(vector, &ptr);
real_type lp_norm(const real_type p) const
real_type l1_norm() const
VectorBase & operator+=(const VectorBase &V)
const IndexSet & ghost_elements() const
PetscScalar mean_value() const
void determine_ghost_indices()
VectorOperation::values last_action
PetscScalar operator*(const VectorBase &vec) const
bool operator==(const VectorBase &v) const
VectorBase & operator*=(const PetscScalar factor)
VectorBase & operator-=(const VectorBase &V)
bool operator!=(const VectorBase &v) const
IndexSet locally_owned_elements() const
bool in_local_range(const size_type index) const
std::size_t memory_consumption() const
internal::VectorReference reference
VectorBase & operator/=(const PetscScalar factor)
friend class internal::VectorReference
PetscScalar operator[](const size_type index) const
MPI_Comm get_mpi_communicator() const
void scale(const VectorBase &scaling_factors)
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< PetscScalar > &values) const
void update_ghost_values() const
std::pair< size_type, size_type > local_range() const
void sadd(const PetscScalar s, const VectorBase &V)
real_type norm_sqr() const
void compress(const VectorOperation::values operation)
void set(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
PetscScalar add_and_dot(const PetscScalar a, const VectorBase &V, const VectorBase &W)
virtual ~VectorBase() override
real_type linfty_norm() const
reference operator()(const size_type index)
const internal::VectorReference const_reference
void add(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)
void write_ascii(const PetscViewerFormat format=PETSC_VIEWER_DEFAULT)
bool has_ghost_elements() const
reference operator[](const size_type index)
virtual void extract_subvector_to(const ArrayView< const types::global_dof_index > &indices, ArrayView< PetscScalar > &elements) const override
size_type locally_owned_size() const
void equ(const PetscScalar a, const VectorBase &V)
void do_set_add_operation(const size_type n_elements, const size_type *indices, const PetscScalar *values, const bool add_values)
VectorBase & operator=(const VectorBase &)
real_type l2_norm() const
types::global_dof_index size_type
size_type size() const override
void swap(VectorBase &u, VectorBase &v)
PetscScalar operator()(const size_type index) const
void extract_subvector_to(const ForwardIterator indices_begin, const ForwardIterator indices_end, OutputIterator values_begin) const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
static ::ExceptionBase & ExcGhostsPresent()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
types::global_dof_index size_type
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
unsigned int global_dof_index