17 #ifndef __deal2__petsc_vector_base_h
18 #define __deal2__petsc_vector_base_h
21 #include <deal.II/base/config.h>
23 #ifdef DEAL_II_WITH_PETSC
25 # include <deal.II/base/subscriptor.h>
26 # include <deal.II/lac/exceptions.h>
27 # include <deal.II/lac/vector.h>
32 # include <petscvec.h>
33 # include <deal.II/base/index_set.h>
35 DEAL_II_NAMESPACE_OPEN
38 template <
typename number>
class Vector;
48 namespace PETScWrappers
97 VectorReference (
const VectorBase &vector,
123 const VectorReference &operator = (
const VectorReference &r)
const;
133 VectorReference &operator = (
const VectorReference &r);
139 const VectorReference &operator = (
const PetscScalar &s)
const;
145 const VectorReference &operator += (
const PetscScalar &s)
const;
151 const VectorReference &operator -= (
const PetscScalar &s)
const;
157 const VectorReference &operator *= (
const PetscScalar &s)
const;
163 const VectorReference &operator /= (
const PetscScalar &s)
const;
171 operator PetscScalar ()
const;
178 <<
"An error with error number " << arg1
179 <<
" occurred while calling a PETSc function");
185 <<
"You tried to access element " << arg1
186 <<
" of a distributed vector, but only elements "
187 << arg2 <<
" through " << arg3
188 <<
" are stored locally and can be accessed.");
194 <<
"You tried to do a "
199 <<
" operation but the vector is currently in "
204 <<
" mode. You first have to call 'compress()'.");
211 const VectorBase &vector;
224 friend class ::PETScWrappers::VectorBase;
267 typedef PetscReal real_type;
312 void compress (::VectorOperation::values operation);
365 size_type
size () const;
398 std::pair<size_type, size_type>
435 operator () (const size_type index);
442 operator () (const size_type index) const;
451 operator [] (const size_type index);
461 operator [] (const size_type index) const;
473 void set (const std::
vector<size_type> &indices,
474 const std::
vector<PetscScalar> &values);
487 std::
vector<PetscScalar> &values) const;
493 template <typename ForwardIterator, typename OutputIterator>
495 const ForwardIterator indices_end,
496 OutputIterator values_begin) const;
504 void add (const std::
vector<size_type> &indices,
505 const std::
vector<PetscScalar> &values);
514 void add (const std::
vector<size_type> &indices,
515 const ::
Vector<PetscScalar> &values);
526 void add (const size_type n_elements,
527 const size_type *indices,
528 const PetscScalar *values);
535 PetscScalar operator * (const
VectorBase &vec) const;
566 real_type
lp_norm (const real_type p) const;
583 real_type
min () const;
588 real_type
max () const;
652 VectorBase &operator *= (const PetscScalar factor);
658 VectorBase &operator /= (const PetscScalar factor);
677 void add (const PetscScalar s);
703 void sadd (const PetscScalar s,
710 void sadd (const PetscScalar s,
717 void sadd (const PetscScalar s,
727 void sadd (const PetscScalar s,
795 void write_ascii (const PetscViewerFormat format = PETSC_VIEWER_DEFAULT) ;
811 void print (std::ostream &out,
812 const
unsigned int precision = 3,
813 const
bool scientific = true,
814 const
bool across = true) const;
852 operator const Vec &() const;
907 friend class internal::VectorReference;
927 const size_type *indices,
928 const PetscScalar *values,
929 const
bool add_values);
965 const VectorReference &
966 VectorReference::operator = (
const VectorReference &r)
const
972 *
this =
static_cast<PetscScalar
> (r);
981 VectorReference::operator = (
const VectorReference &r)
987 *
this =
static_cast<PetscScalar
> (r);
995 const VectorReference &
996 VectorReference::operator = (
const PetscScalar &value)
const
998 Assert ((vector.last_action == VectorOperation::insert)
1000 (vector.last_action == VectorOperation::unknown),
1001 ExcWrongMode (VectorOperation::insert,
1002 vector.last_action));
1006 const PetscInt petsc_i = index;
1009 = VecSetValues (vector, 1, &petsc_i, &value, INSERT_VALUES);
1012 vector.last_action = VectorOperation::insert;
1020 const VectorReference &
1021 VectorReference::operator += (
const PetscScalar &value)
const
1023 Assert ((vector.last_action == VectorOperation::add)
1025 (vector.last_action == VectorOperation::unknown),
1026 ExcWrongMode (VectorOperation::add,
1027 vector.last_action));
1031 vector.last_action = VectorOperation::add;
1040 if (value == PetscScalar())
1044 const PetscInt petsc_i = index;
1046 = VecSetValues (vector, 1, &petsc_i, &value, ADD_VALUES);
1056 const VectorReference &
1057 VectorReference::operator -= (
const PetscScalar &value)
const
1059 Assert ((vector.last_action == VectorOperation::add)
1061 (vector.last_action == VectorOperation::unknown),
1062 ExcWrongMode (VectorOperation::add,
1063 vector.last_action));
1067 vector.last_action = VectorOperation::add;
1076 if (value == PetscScalar())
1081 const PetscInt petsc_i = index;
1082 const PetscScalar subtractand = -value;
1084 = VecSetValues (vector, 1, &petsc_i, &subtractand, ADD_VALUES);
1093 const VectorReference &
1094 VectorReference::operator *= (
const PetscScalar &value)
const
1096 Assert ((vector.last_action == VectorOperation::insert)
1098 (vector.last_action == VectorOperation::unknown),
1099 ExcWrongMode (VectorOperation::insert,
1100 vector.last_action));
1104 vector.last_action = VectorOperation::insert;
1116 const PetscInt petsc_i = index;
1117 const PetscScalar new_value
1118 =
static_cast<PetscScalar
>(*this) * value;
1121 = VecSetValues (vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1130 const VectorReference &
1131 VectorReference::operator /= (
const PetscScalar &value)
const
1133 Assert ((vector.last_action == VectorOperation::insert)
1135 (vector.last_action == VectorOperation::unknown),
1136 ExcWrongMode (VectorOperation::insert,
1137 vector.last_action));
1141 vector.last_action = VectorOperation::insert;
1153 const PetscInt petsc_i = index;
1154 const PetscScalar new_value
1155 =
static_cast<PetscScalar
>(*this) / value;
1158 = VecSetValues (vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1169 VectorBase::in_local_range (
const size_type index)
const
1171 PetscInt begin, end;
1172 const int ierr = VecGetOwnershipRange (static_cast<const Vec &>(vector),
1176 return ((index >= static_cast<size_type>(begin)) &&
1177 (index < static_cast<size_type>(end)));
1183 VectorBase::locally_owned_elements()
const
1188 const std::pair<size_type, size_type> x = local_range();
1189 is.add_range (x.first, x.second);
1197 VectorBase::has_ghost_elements()
const
1205 internal::VectorReference
1206 VectorBase::operator () (
const size_type index)
1208 return internal::VectorReference (*
this, index);
1215 VectorBase::operator () (
const size_type index)
const
1217 return static_cast<PetscScalar
>(internal::VectorReference (*
this, index));
1223 internal::VectorReference
1224 VectorBase::operator [] (
const size_type index)
1226 return operator()(index);
1233 VectorBase::operator [] (
const size_type index)
const
1235 return operator()(index);
1240 VectorBase::get_mpi_communicator ()
const
1242 static MPI_Comm comm;
1243 PetscObjectGetComm((PetscObject)vector, &comm);
1248 void VectorBase::extract_subvector_to (
const std::vector<size_type> &indices,
1249 std::vector<PetscScalar> &values)
const
1251 extract_subvector_to(&(indices[0]), &(indices[0]) + indices.size(), &(values[0]));
1254 template <
typename ForwardIterator,
typename OutputIterator>
1256 void VectorBase::extract_subvector_to (
const ForwardIterator indices_begin,
1257 const ForwardIterator indices_end,
1258 OutputIterator values_begin)
const
1260 const PetscInt n_idx =
static_cast<PetscInt
>(indices_end - indices_begin);
1288 PetscInt begin, end, i;
1289 ierr = VecGetOwnershipRange (vector, &begin, &end);
1292 Vec locally_stored_elements = PETSC_NULL;
1293 ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
1297 ierr = VecGetSize(locally_stored_elements, &lsize);
1301 ierr = VecGetArray(locally_stored_elements, &ptr);
1304 for (i = 0; i < n_idx; i++)
1306 const unsigned int index = *(indices_begin+i);
1307 if ( index>=static_cast<unsigned int>(begin)
1308 && index<static_cast<unsigned int>(end) )
1311 *(values_begin+i) = *(ptr+index-begin);
1316 const unsigned int ghostidx
1317 = ghost_indices.index_within_set(index);
1320 *(values_begin+i) = *(ptr+ghostidx+end-begin);
1324 ierr = VecRestoreArray(locally_stored_elements, &ptr);
1327 ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
1338 PetscInt begin, end;
1339 ierr = VecGetOwnershipRange (vector, &begin, &end);
1343 ierr = VecGetArray(vector, &ptr);
1346 for (PetscInt i = 0; i < n_idx; i++)
1348 const unsigned int index = *(indices_begin+i);
1350 Assert(index>=static_cast<unsigned int>(begin)
1353 *(values_begin+i) = *(ptr+index-begin);
1356 ierr = VecRestoreArray(vector, &ptr);
1365 DEAL_II_NAMESPACE_CLOSE
1367 #endif // DEAL_II_WITH_PETSC
real_type l2_norm() const
real_type linfty_norm() const
#define DeclException2(Exception2, type1, type2, outsequence)
real_type l1_norm() const
void sadd(const PetscScalar s, const VectorBase &V)
void ratio(const VectorBase &a, const VectorBase &b)
bool has_ghost_elements() const
void update_ghost_values() const DEAL_II_DEPRECATED
real_type norm_sqr() const
#define AssertThrow(cond, exc)
void scale(const VectorBase &scaling_factors)
real_type lp_norm(const real_type p) const
void compress() DEAL_II_DEPRECATED
virtual const MPI_Comm & get_mpi_communicator() const
::ExceptionBase & ExcGhostsPresent()
#define DeclException1(Exception1, type1, outsequence)
mutable::VectorOperation::values last_action
unsigned int global_dof_index
std::pair< size_type, size_type > local_range() const
void do_set_add_operation(const size_type n_elements, const size_type *indices, const PetscScalar *values, const bool add_values)
#define Assert(cond, exc)
PetscScalar mean_value() const
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< PetscScalar > &values) const
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
void equ(const PetscScalar a, const VectorBase &V)
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
friend class internal::VectorReference
size_type local_size() const
bool is_non_negative() const
std::size_t memory_consumption() const
void add(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)
void write_ascii(const PetscViewerFormat format=PETSC_VIEWER_DEFAULT)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
::ExceptionBase & ExcInternalError()
IndexSet locally_owned_elements() const
real_type normalize() const
bool in_local_range(const size_type index) const
void set(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)