17#ifdef DEAL_II_WITH_PETSC
34 VectorReference::operator PetscScalar()
const
48 VecGetOwnershipRange(vector.vector, &begin, &end);
51 Vec locally_stored_elements =
nullptr;
52 ierr = VecGhostGetLocalForm(vector.vector, &locally_stored_elements);
56 ierr = VecGetSize(locally_stored_elements, &lsize);
59 const PetscScalar *ptr;
60 ierr = VecGetArrayRead(locally_stored_elements, &ptr);
65 if (index >=
static_cast<size_type
>(begin) &&
66 index <
static_cast<size_type
>(end))
69 value = *(ptr + index - begin);
74 Assert(vector.ghost_indices.is_element(index),
76 "You are trying to access an element of a vector "
77 "that is neither a locally owned element nor a "
78 "ghost element of the vector."));
79 const size_type ghostidx =
80 vector.ghost_indices.index_within_set(index);
83 value = *(ptr + ghostidx + end - begin);
86 ierr = VecRestoreArrayRead(locally_stored_elements, &ptr);
90 VecGhostRestoreLocalForm(vector.vector, &locally_stored_elements);
102 PetscErrorCode ierr = VecGetOwnershipRange(vector.vector, &begin, &end);
105 AssertThrow((index >=
static_cast<size_type
>(begin)) &&
106 (index <
static_cast<size_type
>(end)),
107 ExcAccessToNonlocalElement(index, begin, end - 1));
109 const PetscScalar *ptr;
111 ierr = VecGetArrayRead(vector.vector, &ptr);
113 value = *(ptr + index - begin);
114 ierr = VecRestoreArrayRead(vector.vector, &ptr);
132 , ghost_indices(v.ghost_indices)
149 const PetscErrorCode ierr =
150 PetscObjectReference(
reinterpret_cast<PetscObject
>(
vector));
159 const PetscErrorCode ierr = VecDestroy(&
vector);
170 PetscErrorCode ierr =
171 PetscObjectReference(
reinterpret_cast<PetscObject
>(v));
173 ierr = VecDestroy(&
vector);
183 template <
typename Iterator,
typename OutType>
184 class ConvertingIterator
189 using difference_type =
190 typename std::iterator_traits<Iterator>::difference_type;
191 using value_type = OutType;
192 using pointer = OutType *;
193 using reference = OutType &;
194 using iterator_category = std::forward_iterator_tag;
196 ConvertingIterator(
const Iterator &iterator)
203 return static_cast<OutType
>(std::real(*m_iterator));
216 ConvertingIterator old = *
this;
222 operator==(
const ConvertingIterator &other)
const
224 return this->m_iterator == other.m_iterator;
228 operator!=(
const ConvertingIterator &other)
const
230 return this->m_iterator != other.m_iterator;
250 ierr = VecGhostGetLocalForm(
vector, &ghosted_vec);
252 if (ghosted_vec && ghosted_vec !=
vector)
256 PetscInt ghost_start_index, end_index, n_elements_stored_locally;
258 ierr = VecGhostRestoreLocalForm(
vector, &ghosted_vec);
261 ierr = VecGetOwnershipRange(
vector, &ghost_start_index, &end_index);
263 ierr = VecDuplicate(
vector, &tvector);
265 ierr = VecGetArray(tvector, &array);
273 for (PetscInt i = 0; i < end_index - ghost_start_index; i++)
275 Assert(
static_cast<PetscInt
>(std::real(
static_cast<PetscScalar
>(
276 ghost_start_index + i))) == (ghost_start_index + i),
278 array[i] = ghost_start_index + i;
281 ierr = VecRestoreArray(tvector, &array);
283 ierr = VecGhostUpdateBegin(tvector, INSERT_VALUES, SCATTER_FORWARD);
285 ierr = VecGhostUpdateEnd(tvector, INSERT_VALUES, SCATTER_FORWARD);
287 ierr = VecGhostGetLocalForm(tvector, &ghosted_vec);
289 ierr = VecGetLocalSize(ghosted_vec, &n_elements_stored_locally);
291 ierr = VecGetArrayRead(ghosted_vec, (
const PetscScalar **)&array);
303 ConvertingIterator<PetscScalar *, types::global_dof_index> begin_ghosts(
304 &array[end_index - ghost_start_index]);
305 ConvertingIterator<PetscScalar *, types::global_dof_index> end_ghosts(
306 &array[n_elements_stored_locally]);
307 if (std::is_sorted(&array[end_index - ghost_start_index],
308 &array[n_elements_stored_locally],
309 [](PetscScalar left, PetscScalar right) {
310 return static_cast<PetscInt
>(std::real(left)) <
311 static_cast<PetscInt
>(std::real(right));
318 std::vector<PetscInt> sorted_indices(begin_ghosts, end_ghosts);
319 std::sort(sorted_indices.begin(), sorted_indices.end());
321 sorted_indices.end());
325 ierr = VecRestoreArrayRead(ghosted_vec, (
const PetscScalar **)&array);
327 ierr = VecGhostRestoreLocalForm(tvector, &ghosted_vec);
329 ierr = VecDestroy(&tvector);
334 ierr = VecGhostRestoreLocalForm(
vector, &ghosted_vec);
343 const PetscErrorCode ierr = VecDestroy(&
vector);
358 PetscErrorCode ierr = VecCopy(v,
vector);
373 PetscErrorCode ierr = VecSet(
vector, s);
379 ierr = VecGhostGetLocalForm(
vector, &ghost);
382 ierr = VecSet(ghost, s);
385 ierr = VecGhostRestoreLocalForm(
vector, &ghost);
400 const PetscErrorCode ierr = VecEqual(
vector, v.
vector, &flag);
403 return (flag == PETSC_TRUE);
414 const PetscErrorCode ierr = VecEqual(
vector, v.
vector, &flag);
417 return (flag == PETSC_FALSE);
426 const PetscErrorCode ierr = VecGetSize(
vector, &sz);
438 const PetscErrorCode ierr = VecGetLocalSize(
vector, &sz);
446 std::pair<VectorBase::size_type, VectorBase::size_type>
450 const PetscErrorCode ierr =
451 VecGetOwnershipRange(
static_cast<const Vec &
>(
vector), &begin, &end);
454 return std::make_pair(begin, end);
461 const std::vector<PetscScalar> &values)
463 Assert(indices.size() == values.size(),
464 ExcMessage(
"Function called with arguments of different sizes"));
472 const std::vector<PetscScalar> &values)
474 Assert(indices.size() == values.size(),
475 ExcMessage(
"Function called with arguments of different sizes"));
483 const ::Vector<PetscScalar> &values)
485 Assert(indices.size() == values.size(),
486 ExcMessage(
"Function called with arguments of different sizes"));
495 const PetscScalar *values)
516 const PetscErrorCode ierr = VecDot(vec.
vector,
vector, &result);
539 ExcMessage(
"Calling compress() is only useful if a vector "
540 "has been written into, but this is a vector with ghost "
541 "elements and consequently is read-only. It does "
542 "not make sense to call compress() for such "
552 int all_int_last_action;
554 const int ierr = MPI_Allreduce(&my_int_last_action,
555 &all_int_last_action,
565 "Error: not all processors agree on the last "
566 "VectorOperation before this compress() call."));
573 "Missing compress() or calling with wrong VectorOperation argument."));
589 PetscErrorCode ierr = VecAssemblyBegin(
vector);
591 ierr = VecAssemblyEnd(
vector);
619 const PetscErrorCode ierr = VecSum(
vector, &sum);
621 return sum /
static_cast<PetscReal
>(
size());
626 const PetscScalar *start_ptr;
627 PetscErrorCode ierr = VecGetArrayRead(
vector, &start_ptr);
630 PetscScalar mean = 0;
632 PetscScalar sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
637 const PetscScalar *ptr = start_ptr;
652 static_cast<PetscReal
>(
size());
657 ierr = VecRestoreArrayRead(
vector, &start_ptr);
669 const PetscErrorCode ierr = VecNorm(
vector, NORM_1, &d);
682 const PetscErrorCode ierr = VecNorm(
vector, NORM_2, &d);
695 const PetscScalar *start_ptr;
696 PetscErrorCode ierr = VecGetArrayRead(
vector, &start_ptr);
701 real_type sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0;
706 const PetscScalar *ptr = start_ptr;
726 ierr = VecRestoreArrayRead(
vector, &start_ptr);
739 const PetscErrorCode ierr = VecNorm(
vector, NORM_INFINITY, &d);
752 const PetscScalar *start_ptr;
753 PetscErrorCode ierr = VecGetArrayRead(
vector, &start_ptr);
756 const PetscScalar *ptr = start_ptr,
771 ierr = VecRestoreArrayRead(
vector, &start_ptr);
780 template <
typename T>
789 template <
typename T>
795 "whether it is non-negative."));
808 const PetscErrorCode ierr = VecScale(
vector, a);
822 const PetscScalar factor = 1. / a;
825 const PetscErrorCode ierr = VecScale(
vector, factor);
837 const PetscErrorCode ierr = VecAXPY(
vector, 1, v);
849 const PetscErrorCode ierr = VecAXPY(
vector, -1, v);
863 const PetscErrorCode ierr = VecShift(
vector, s);
875 const PetscErrorCode ierr = VecAXPY(
vector, a, v);
891 const PetscScalar weights[2] = {a, b};
892 Vec addends[2] = {v.
vector, w.vector};
894 const PetscErrorCode ierr = VecMAXPY(
vector, 2, weights, addends);
906 const PetscErrorCode ierr = VecAYPX(
vector, s, v);
934 const PetscErrorCode ierr = VecPointwiseMult(
vector, factors,
vector);
948 const PetscErrorCode ierr = VecAXPBY(
vector, a, 0.0, v.
vector);
961 PetscErrorCode ierr =
962 PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(
comm), format);
966 ierr = VecView(
vector, PETSC_VIEWER_STDOUT_(
comm));
968 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(
comm));
976 const unsigned int precision,
977 const bool scientific,
978 const bool across)
const
984 const PetscScalar *val;
985 PetscErrorCode ierr = VecGetArrayRead(
vector, &val);
990 const std::ios::fmtflags old_flags = out.flags();
991 const unsigned int old_precision = out.precision(precision);
993 out.precision(precision);
995 out.setf(std::ios::scientific, std::ios::floatfield);
997 out.setf(std::ios::fixed, std::ios::floatfield);
1001 out << val[i] <<
' ';
1004 out << val[i] << std::endl;
1008 out.flags(old_flags);
1009 out.precision(old_precision);
1013 ierr = VecRestoreArrayRead(
vector, &val);
1024 std::swap(this->vector, v.vector);
1025 std::swap(this->ghosted, v.ghosted);
1026 std::swap(this->last_action, v.last_action);
1029 this->ghost_indices = v.ghost_indices;
1030 v.ghost_indices = t;
1034 VectorBase::operator
const Vec &()
const
1050 std::size_t mem =
sizeof(Vec) +
sizeof(
last_action) +
1072 const PetscScalar *values,
1073 const bool add_values)
1078 internal::VectorReference::ExcWrongMode(action,
last_action));
1081 std::vector<PetscInt> petsc_indices(n_elements);
1082 for (
size_type i = 0; i < n_elements; ++i)
1084 const auto petsc_index =
static_cast<PetscInt
>(indices[i]);
1086 petsc_indices[i] = petsc_index;
1089 const InsertMode mode = (add_values ? ADD_VALUES : INSERT_VALUES);
1090 const PetscErrorCode ierr = VecSetValues(
1091 vector, petsc_indices.size(), petsc_indices.data(), 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 swap(VectorBase &v) noexcept
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
size_type size() const override
std::enable_if_t< std::is_floating_point_v< T > &&std::is_floating_point_v< U > &&!std::is_same_v< T, 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
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define AssertIntegerConversion(index1, index2)
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcGhostsPresent()
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define AssertThrowMPI(error_code)
#define AssertNothrow(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
bool is_non_negative(const T &t)
T sum(const T &t, const MPI_Comm mpi_communicator)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
SynchronousIterators< Iterators > & operator++(SynchronousIterators< Iterators > &a)