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;
393 std::pair<size_type, size_type>
472 set(
const std::vector<size_type> & indices,
473 const std::vector<PetscScalar> &
values);
492 std::vector<PetscScalar> &
values)
const;
521 template <
typename ForwardIterator,
typename OutputIterator>
524 const ForwardIterator indices_end,
525 OutputIterator values_begin)
const;
532 add(
const std::vector<size_type> & indices,
533 const std::vector<PetscScalar> &
values);
540 add(
const std::vector<size_type> & indices,
541 const ::Vector<PetscScalar> &
values);
551 const PetscScalar *
values);
698 add(
const PetscScalar s);
710 add(
const PetscScalar a,
749 write_ascii(
const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
759 print(std::ostream & out,
760 const unsigned int precision = 3,
761 const bool scientific =
true,
762 const bool across =
true)
const;
786 operator const Vec &()
const;
848 const PetscScalar *
values,
849 const bool add_values);
872 inline VectorReference::VectorReference(
const VectorBase &vector,
879 inline const VectorReference &
880 VectorReference::operator=(
const VectorReference &r)
const
886 *
this =
static_cast<PetscScalar
>(r);
893 inline VectorReference &
894 VectorReference::operator=(
const VectorReference &r)
900 *
this =
static_cast<PetscScalar
>(r);
907 inline const VectorReference &
908 VectorReference::operator=(
const PetscScalar &value)
const
916 const PetscInt petsc_i =
index;
918 const PetscErrorCode ierr =
919 VecSetValues(vector, 1, &petsc_i, &value, INSERT_VALUES);
929 inline const VectorReference &
930 VectorReference::operator+=(
const PetscScalar &value)
const
947 if (value == PetscScalar())
951 const PetscInt petsc_i =
index;
952 const PetscErrorCode ierr =
953 VecSetValues(vector, 1, &petsc_i, &value, ADD_VALUES);
962 inline const VectorReference &
963 VectorReference::operator-=(
const PetscScalar &value)
const
980 if (value == PetscScalar())
985 const PetscInt petsc_i =
index;
986 const PetscScalar subtractand = -
value;
987 const PetscErrorCode ierr =
988 VecSetValues(vector, 1, &petsc_i, &subtractand, ADD_VALUES);
996 inline const VectorReference &
997 VectorReference::operator*=(
const PetscScalar &value)
const
1017 const PetscInt petsc_i =
index;
1018 const PetscScalar new_value =
static_cast<PetscScalar
>(*this) *
value;
1020 const PetscErrorCode ierr =
1021 VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1029 inline const VectorReference &
1030 VectorReference::operator/=(
const PetscScalar &value)
const
1050 const PetscInt petsc_i =
index;
1051 const PetscScalar new_value =
static_cast<PetscScalar
>(*this) /
value;
1053 const PetscErrorCode ierr =
1054 VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1063 VectorReference::real()
const
1065 # ifndef PETSC_USE_COMPLEX
1066 return static_cast<PetscScalar
>(*this);
1068 return PetscRealPart(
static_cast<PetscScalar
>(*
this));
1075 VectorReference::imag()
const
1077 # ifndef PETSC_USE_COMPLEX
1078 return PetscReal(0);
1080 return PetscImaginaryPart(
static_cast<PetscScalar
>(*
this));
1087 VectorBase::in_local_range(
const size_type index)
const
1090 const PetscErrorCode ierr =
1091 VecGetOwnershipRange(
static_cast<const Vec &
>(vector), &
begin, &
end);
1100 VectorBase::locally_owned_elements()
const
1105 const std::pair<size_type, size_type> x = local_range();
1106 is.add_range(x.first, x.second);
1113 VectorBase::has_ghost_elements()
const
1121 VectorBase::update_ghost_values()
const
1126 inline internal::VectorReference
1127 VectorBase::operator()(
const size_type index)
1129 return internal::VectorReference(*
this, index);
1135 VectorBase::operator()(
const size_type index)
const
1137 return static_cast<PetscScalar
>(internal::VectorReference(*
this, index));
1142 inline internal::VectorReference
1143 VectorBase::operator[](
const size_type index)
1145 return operator()(index);
1151 VectorBase::operator[](
const size_type index)
const
1153 return operator()(index);
1156 inline const MPI_Comm &
1157 VectorBase::get_mpi_communicator()
const
1159 static MPI_Comm
comm = PETSC_COMM_SELF;
1160 MPI_Comm pcomm = PetscObjectComm(
reinterpret_cast<PetscObject
>(vector));
1161 if (pcomm != MPI_COMM_NULL)
1167 VectorBase::extract_subvector_to(
const std::vector<size_type> &indices,
1168 std::vector<PetscScalar> &
values)
const
1172 extract_subvector_to(indices.begin(), indices.end(),
values.begin());
1175 template <
typename ForwardIterator,
typename OutputIterator>
1177 VectorBase::extract_subvector_to(
const ForwardIterator indices_begin,
1178 const ForwardIterator indices_end,
1179 OutputIterator values_begin)
const
1181 const PetscInt n_idx =
static_cast<PetscInt
>(indices_end - indices_begin);
1207 PetscErrorCode ierr = VecGetOwnershipRange(vector, &
begin, &
end);
1210 Vec locally_stored_elements =
nullptr;
1211 ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
1215 ierr = VecGetSize(locally_stored_elements, &lsize);
1218 const PetscScalar *ptr;
1219 ierr = VecGetArrayRead(locally_stored_elements, &ptr);
1222 for (PetscInt i = 0; i < n_idx; ++i)
1224 const unsigned int index = *(indices_begin + i);
1225 if (index >=
static_cast<unsigned int>(
begin) &&
1226 index < static_cast<unsigned int>(
end))
1229 *(values_begin + i) = *(ptr + index -
begin);
1234 const unsigned int ghostidx =
1235 ghost_indices.index_within_set(index);
1238 *(values_begin + i) = *(ptr + ghostidx +
end -
begin);
1242 ierr = VecRestoreArrayRead(locally_stored_elements, &ptr);
1245 ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
1254 PetscErrorCode ierr = VecGetOwnershipRange(vector, &
begin, &
end);
1257 const PetscScalar *ptr;
1258 ierr = VecGetArrayRead(vector, &ptr);
1261 for (PetscInt i = 0; i < n_idx; ++i)
1263 const unsigned int index = *(indices_begin + i);
1265 Assert(index >=
static_cast<unsigned int>(
begin) &&
1266 index <
static_cast<unsigned int>(
end),
1267 ExcMessage(
"You are accessing elements of a vector without "
1268 "ghost elements that are not actually owned by "
1269 "this vector. A typical case where this may "
1270 "happen is if you are passing a non-ghosted "
1271 "(completely distributed) vector to a function "
1272 "that expects a vector that stores ghost "
1273 "elements for all locally relevant or locally "
1274 "active vector entries."));
1276 *(values_begin + i) = *(ptr + index -
begin);
1279 ierr = VecRestoreArrayRead(vector, &ptr);
real_type lp_norm(const real_type p) const
real_type l1_norm() const
VectorBase & operator+=(const VectorBase &V)
PetscScalar mean_value() const
size_type local_size() const
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)
bool is_non_negative() const
friend class internal::VectorReference
PetscScalar operator[](const size_type index) 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)
const MPI_Comm & get_mpi_communicator() const
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)
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
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_DEPRECATED
#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 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