18 #ifdef DEAL_II_WITH_PETSC
35 VectorReference::operator PetscScalar()
const
49 VecGetOwnershipRange(vector.vector, &
begin, &
end);
52 Vec locally_stored_elements =
nullptr;
53 ierr = VecGhostGetLocalForm(vector.vector, &locally_stored_elements);
57 ierr = VecGetSize(locally_stored_elements, &lsize);
60 const PetscScalar *ptr;
61 ierr = VecGetArrayRead(locally_stored_elements, &ptr);
70 value = *(ptr + index -
begin);
75 Assert(vector.ghost_indices.is_element(index),
77 "You are trying to access an element of a vector "
78 "that is neither a locally owned element nor a "
79 "ghost element of the vector."));
81 vector.ghost_indices.index_within_set(index);
84 value = *(ptr + ghostidx +
end -
begin);
87 ierr = VecRestoreArrayRead(locally_stored_elements, &ptr);
91 VecGhostRestoreLocalForm(vector.vector, &locally_stored_elements);
103 PetscErrorCode ierr = VecGetOwnershipRange(vector.vector, &
begin, &
end);
108 ExcAccessToNonlocalElement(index,
begin,
end - 1));
110 const PetscScalar *ptr;
112 ierr = VecGetArrayRead(vector.vector, &ptr);
114 value = *(ptr + index -
begin);
115 ierr = VecRestoreArrayRead(vector.vector, &ptr);
134 , ghost_indices(v.ghost_indices)
152 const PetscErrorCode ierr =
153 PetscObjectReference(
reinterpret_cast<PetscObject
>(
vector));
163 const PetscErrorCode ierr = VecDestroy(&
vector);
175 PetscErrorCode ierr =
176 PetscObjectReference(
reinterpret_cast<PetscObject
>(v));
178 ierr = VecDestroy(&
vector);
188 template <
typename Iterator,
typename OutType>
189 class ConvertingIterator
194 using difference_type =
195 typename std::iterator_traits<Iterator>::difference_type;
196 using value_type = OutType;
197 using pointer = OutType *;
198 using reference = OutType &;
199 using iterator_category = std::forward_iterator_tag;
201 ConvertingIterator(
const Iterator &iterator)
208 return static_cast<OutType
>(std::real(*
m_iterator));
221 ConvertingIterator old = *
this;
227 operator==(
const ConvertingIterator &other)
const
229 return this->m_iterator == other.m_iterator;
233 operator!=(
const ConvertingIterator &other)
const
235 return this->m_iterator != other.m_iterator;
255 ierr = VecGhostGetLocalForm(
vector, &ghosted_vec);
257 if (ghosted_vec && ghosted_vec !=
vector)
261 PetscInt ghost_start_index, end_index, n_elements_stored_locally;
263 ierr = VecGhostRestoreLocalForm(
vector, &ghosted_vec);
266 ierr = VecGetOwnershipRange(
vector, &ghost_start_index, &end_index);
270 ierr = VecGetArray(tvector, &array);
278 for (PetscInt i = 0; i < end_index - ghost_start_index; i++)
280 Assert(
static_cast<PetscInt
>(std::real(
static_cast<PetscScalar
>(
281 ghost_start_index + i))) == (ghost_start_index + i),
283 array[i] = ghost_start_index + i;
286 ierr = VecRestoreArray(tvector, &array);
288 ierr = VecGhostUpdateBegin(tvector, INSERT_VALUES, SCATTER_FORWARD);
290 ierr = VecGhostUpdateEnd(tvector, INSERT_VALUES, SCATTER_FORWARD);
292 ierr = VecGhostGetLocalForm(tvector, &ghosted_vec);
294 ierr = VecGetLocalSize(ghosted_vec, &n_elements_stored_locally);
296 ierr = VecGetArrayRead(ghosted_vec, (
const PetscScalar **)&array);
308 ConvertingIterator<PetscScalar *, types::global_dof_index> begin_ghosts(
309 &array[end_index - ghost_start_index]);
310 ConvertingIterator<PetscScalar *, types::global_dof_index> end_ghosts(
311 &array[n_elements_stored_locally]);
312 if (std::is_sorted(&array[end_index - ghost_start_index],
313 &array[n_elements_stored_locally],
314 [](PetscScalar left, PetscScalar right) {
315 return static_cast<PetscInt
>(std::real(left)) <
316 static_cast<PetscInt
>(std::real(right));
323 std::vector<PetscInt> sorted_indices(begin_ghosts, end_ghosts);
324 std::sort(sorted_indices.begin(), sorted_indices.end());
326 sorted_indices.end());
330 ierr = VecRestoreArrayRead(ghosted_vec, (
const PetscScalar **)&array);
332 ierr = VecGhostRestoreLocalForm(tvector, &ghosted_vec);
334 ierr = VecDestroy(&tvector);
339 ierr = VecGhostRestoreLocalForm(
vector, &ghosted_vec);
348 const PetscErrorCode ierr = VecDestroy(&
vector);
363 PetscErrorCode ierr = VecCopy(v,
vector);
378 PetscErrorCode ierr = VecSet(
vector, s);
384 ierr = VecGhostGetLocalForm(
vector, &ghost);
387 ierr = VecSet(ghost, s);
390 ierr = VecGhostRestoreLocalForm(
vector, &ghost);
405 const PetscErrorCode ierr = VecEqual(
vector, v.
vector, &flag);
408 return (flag == PETSC_TRUE);
419 const PetscErrorCode ierr = VecEqual(
vector, v.
vector, &flag);
422 return (flag == PETSC_FALSE);
431 const PetscErrorCode ierr = VecGetSize(
vector, &sz);
443 const PetscErrorCode ierr = VecGetLocalSize(
vector, &sz);
451 std::pair<VectorBase::size_type, VectorBase::size_type>
455 const PetscErrorCode ierr =
456 VecGetOwnershipRange(
static_cast<const Vec &
>(
vector), &
begin, &
end);
466 const std::vector<PetscScalar> &
values)
469 ExcMessage(
"Function called with arguments of different sizes"));
477 const std::vector<PetscScalar> &
values)
480 ExcMessage(
"Function called with arguments of different sizes"));
488 const ::Vector<PetscScalar> &
values)
491 ExcMessage(
"Function called with arguments of different sizes"));
500 const PetscScalar *
values)
521 const PetscErrorCode ierr = VecDot(vec.
vector,
vector, &result);
548 int all_int_last_action;
550 const int ierr = MPI_Allreduce(&my_int_last_action,
551 &all_int_last_action,
560 ExcMessage(
"Error: not all processors agree on the last "
561 "VectorOperation before this compress() call."));
568 "Missing compress() or calling with wrong VectorOperation argument."));
584 PetscErrorCode ierr = VecAssemblyBegin(
vector);
586 ierr = VecAssemblyEnd(
vector);
614 const PetscErrorCode ierr = VecSum(
vector, &
sum);
616 return sum /
static_cast<PetscReal
>(
size());
621 const PetscScalar *start_ptr;
622 PetscErrorCode ierr = VecGetArrayRead(
vector, &start_ptr);
625 PetscScalar
mean = 0;
627 PetscScalar sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
632 const PetscScalar *ptr = start_ptr;
633 const PetscScalar *eptr = ptr + (
size() / 4) * 4;
642 while (ptr != start_ptr +
size())
645 mean = (sum0 + sum1 + sum2 + sum3) /
static_cast<PetscReal
>(
size());
650 ierr = VecRestoreArrayRead(
vector, &start_ptr);
662 const PetscErrorCode ierr = VecNorm(
vector, NORM_1, &
d);
675 const PetscErrorCode ierr = VecNorm(
vector, NORM_2, &
d);
688 const PetscScalar *start_ptr;
689 PetscErrorCode ierr = VecGetArrayRead(
vector, &start_ptr);
694 real_type sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
699 const PetscScalar *ptr = start_ptr;
700 const PetscScalar *eptr = ptr + (
size() / 4) * 4;
709 while (ptr != start_ptr +
size())
712 norm = std::pow(sum0 + sum1 + sum2 + sum3, 1. / p);
717 ierr = VecRestoreArrayRead(
vector, &start_ptr);
730 const PetscErrorCode ierr = VecNorm(
vector, NORM_INFINITY, &
d);
743 const PetscScalar *start_ptr;
744 PetscErrorCode ierr = VecGetArrayRead(
vector, &start_ptr);
747 const PetscScalar *ptr = start_ptr,
762 ierr = VecRestoreArrayRead(
vector, &start_ptr);
771 template <
typename T>
780 template <
typename T>
786 "whether it is non-negative."))
return true;
798 const PetscErrorCode ierr = VecScale(
vector, a);
812 const PetscScalar factor = 1. / a;
815 const PetscErrorCode ierr = VecScale(
vector, factor);
827 const PetscErrorCode ierr = VecAXPY(
vector, 1, v);
839 const PetscErrorCode ierr = VecAXPY(
vector, -1, v);
853 const PetscErrorCode ierr = VecShift(
vector, s);
865 const PetscErrorCode ierr = VecAXPY(
vector, a, v);
881 const PetscScalar weights[2] = {a,
b};
882 Vec addends[2] = {v.
vector,
w.vector};
884 const PetscErrorCode ierr = VecMAXPY(
vector, 2, weights, addends);
896 const PetscErrorCode ierr = VecAYPX(
vector, s, v);
924 const PetscErrorCode ierr = VecPointwiseMult(
vector, factors,
vector);
938 const PetscErrorCode ierr = VecAXPBY(
vector, a, 0.0, v.
vector);
948 MPI_Comm
comm = PetscObjectComm((PetscObject)
vector);
951 PetscErrorCode ierr =
952 PetscViewerSetFormat(PETSC_VIEWER_STDOUT_(
comm), format);
956 ierr = VecView(
vector, PETSC_VIEWER_STDOUT_(
comm));
964 const unsigned int precision,
965 const bool scientific,
966 const bool across)
const
972 const PetscScalar *val;
973 PetscErrorCode ierr = VecGetArrayRead(
vector, &val);
978 const std::ios::fmtflags old_flags = out.flags();
979 const unsigned int old_precision = out.precision(precision);
981 out.precision(precision);
983 out.setf(std::ios::scientific, std::ios::floatfield);
985 out.setf(std::ios::fixed, std::ios::floatfield);
989 out << val[i] <<
' ';
992 out << val[i] << std::endl;
996 out.flags(old_flags);
997 out.precision(old_precision);
1001 ierr = VecRestoreArrayRead(
vector, &val);
1022 VectorBase::operator
const Vec &()
const
1038 std::size_t mem =
sizeof(Vec) +
sizeof(
last_action) +
1060 const PetscScalar *
values,
1061 const bool add_values)
1066 internal::VectorReference::ExcWrongMode(action,
last_action));
1074 if (n_elements != 0)
1076 const PetscInt *petsc_indices =
1077 reinterpret_cast<const PetscInt *
>(indices);
1079 const InsertMode mode = (add_values ? ADD_VALUES : INSERT_VALUES);
1080 const PetscErrorCode ierr =
1081 VecSetValues(
vector, n_elements, petsc_indices,
values, mode);
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
size_type n_elements() const
void set_size(const size_type size)
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
real_type lp_norm(const real_type p) const
real_type l1_norm() const
VectorBase & operator+=(const VectorBase &V)
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
std::size_t memory_consumption() const
VectorBase & operator/=(const PetscScalar factor)
MPI_Comm get_mpi_communicator() const
void scale(const VectorBase &scaling_factors)
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
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
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
std::enable_if_t< std::is_floating_point_v< T > &&std::is_floating_point_v< U >, typename ProductType< std::complex< T >, std::complex< U > >::type > operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define AssertThrowMPI(error_code)
#define AssertNothrow(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcGhostsPresent()
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)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
void swap(MemorySpaceData< T, MemorySpace > &u, MemorySpaceData< T, MemorySpace > &v)
bool is_non_negative(const T &t)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)