Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
petsc_vector_base.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_petsc_vector_base_h
17 # define dealii_petsc_vector_base_h
18 
19 
20 # include <deal.II/base/config.h>
21 
22 # ifdef DEAL_II_WITH_PETSC
23 
24 # include <deal.II/base/index_set.h>
26 
27 # include <deal.II/lac/exceptions.h>
28 # include <deal.II/lac/vector.h>
30 
31 # include <petscvec.h>
32 
33 # include <utility>
34 # include <vector>
35 
37 
38 // forward declaration
39 # ifndef DOXYGEN
40 template <typename number>
41 class Vector;
42 
43 namespace PETScWrappers
44 {
45  class VectorBase;
46 }
47 # endif
48 
56 namespace PETScWrappers
57 {
67  namespace internal
68  {
81  class VectorReference
82  {
83  public:
88 
89  private:
94  VectorReference(const VectorBase &vector, const size_type index);
95 
96 
97  public:
109  const VectorReference &
110  operator=(const VectorReference &r) const;
111 
117  VectorReference &
118  operator=(const VectorReference &r);
119 
123  const VectorReference &
124  operator=(const PetscScalar &s) const;
125 
129  const VectorReference &
130  operator+=(const PetscScalar &s) const;
131 
135  const VectorReference &
136  operator-=(const PetscScalar &s) const;
137 
141  const VectorReference &
142  operator*=(const PetscScalar &s) const;
143 
147  const VectorReference &
148  operator/=(const PetscScalar &s) const;
149 
153  PetscReal
154  real() const;
155 
162  PetscReal
163  imag() const;
164 
169  operator PetscScalar() const;
173  DeclException3(ExcAccessToNonlocalElement,
174  int,
175  int,
176  int,
177  << "You tried to access element " << arg1
178  << " of a distributed vector, but only elements " << arg2
179  << " through " << arg3
180  << " are stored locally and can be accessed.");
184  DeclException2(ExcWrongMode,
185  int,
186  int,
187  << "You tried to do a "
188  << (arg1 == 1 ? "'set'" : (arg1 == 2 ? "'add'" : "???"))
189  << " operation but the vector is currently in "
190  << (arg2 == 1 ? "'set'" : (arg2 == 2 ? "'add'" : "???"))
191  << " mode. You first have to call 'compress()'.");
192 
193  private:
197  const VectorBase &vector;
198 
202  const size_type index;
203 
204  // Make the vector class a friend, so that it can create objects of the
205  // present type.
206  friend class ::PETScWrappers::VectorBase;
207  };
208  } // namespace internal
240  class VectorBase : public Subscriptor
241  {
242  public:
248  using value_type = PetscScalar;
249  using real_type = PetscReal;
251  using reference = internal::VectorReference;
252  using const_reference = const internal::VectorReference;
253 
258  VectorBase();
259 
264  VectorBase(const VectorBase &v);
265 
271  explicit VectorBase(const Vec &v);
272 
277  VectorBase &
278  operator=(const VectorBase &) = delete;
279 
283  virtual ~VectorBase() override;
284 
289  virtual void
290  clear();
291 
302  void
303  compress(const VectorOperation::values operation);
304 
318  VectorBase &
319  operator=(const PetscScalar s);
320 
326  bool
327  operator==(const VectorBase &v) const;
328 
334  bool
335  operator!=(const VectorBase &v) const;
336 
340  size_type
341  size() const;
342 
351  size_type
352  local_size() const;
353 
362  std::pair<size_type, size_type>
363  local_range() const;
364 
369  bool
370  in_local_range(const size_type index) const;
371 
385  IndexSet
386  locally_owned_elements() const;
387 
394  bool
395  has_ghost_elements() const;
396 
403  void
404  update_ghost_values() const;
405 
409  reference
410  operator()(const size_type index);
411 
415  PetscScalar
416  operator()(const size_type index) const;
417 
423  reference operator[](const size_type index);
424 
430  PetscScalar operator[](const size_type index) const;
431 
438  void
439  set(const std::vector<size_type> & indices,
440  const std::vector<PetscScalar> &values);
441 
457  void
458  extract_subvector_to(const std::vector<size_type> &indices,
459  std::vector<PetscScalar> & values) const;
460 
488  template <typename ForwardIterator, typename OutputIterator>
489  void
490  extract_subvector_to(const ForwardIterator indices_begin,
491  const ForwardIterator indices_end,
492  OutputIterator values_begin) const;
493 
498  void
499  add(const std::vector<size_type> & indices,
500  const std::vector<PetscScalar> &values);
501 
506  void
507  add(const std::vector<size_type> & indices,
508  const ::Vector<PetscScalar> &values);
509 
515  void
516  add(const size_type n_elements,
517  const size_type * indices,
518  const PetscScalar *values);
519 
526  PetscScalar operator*(const VectorBase &vec) const;
527 
531  real_type
532  norm_sqr() const;
533 
537  PetscScalar
538  mean_value() const;
539 
547  real_type
548  l1_norm() const;
549 
554  real_type
555  l2_norm() const;
556 
561  real_type
562  lp_norm(const real_type p) const;
563 
568  real_type
569  linfty_norm() const;
570 
590  PetscScalar
591  add_and_dot(const PetscScalar a, const VectorBase &V, const VectorBase &W);
592 
596  real_type
597  min() const;
598 
602  real_type
603  max() const;
604 
610  bool
611  all_zero() const;
612 
618  bool
619  is_non_negative() const;
620 
624  VectorBase &
625  operator*=(const PetscScalar factor);
626 
630  VectorBase &
631  operator/=(const PetscScalar factor);
632 
636  VectorBase &
637  operator+=(const VectorBase &V);
638 
642  VectorBase &
643  operator-=(const VectorBase &V);
644 
649  void
650  add(const PetscScalar s);
651 
655  void
656  add(const PetscScalar a, const VectorBase &V);
657 
661  void
662  add(const PetscScalar a,
663  const VectorBase &V,
664  const PetscScalar b,
665  const VectorBase &W);
666 
670  void
671  sadd(const PetscScalar s, const VectorBase &V);
672 
676  void
677  sadd(const PetscScalar s, const PetscScalar a, const VectorBase &V);
678 
684  void
685  scale(const VectorBase &scaling_factors);
686 
690  void
691  equ(const PetscScalar a, const VectorBase &V);
692 
700  void
701  write_ascii(const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
702 
710  void
711  print(std::ostream & out,
712  const unsigned int precision = 3,
713  const bool scientific = true,
714  const bool across = true) const;
715 
728  void
729  swap(VectorBase &v);
730 
738  operator const Vec &() const;
739 
743  std::size_t
744  memory_consumption() const;
745 
750  virtual const MPI_Comm &
751  get_mpi_communicator() const;
752 
753  protected:
758  Vec vector;
759 
765  bool ghosted;
766 
773 
780 
781  // Make the reference class a friend.
782  friend class internal::VectorReference;
783 
790 
796  void
797  do_set_add_operation(const size_type n_elements,
798  const size_type * indices,
799  const PetscScalar *values,
800  const bool add_values);
801  };
802 
803 
804 
805  // ------------------- inline and template functions --------------
806 
815  inline void
817  {
818  u.swap(v);
819  }
820 
821 # ifndef DOXYGEN
822  namespace internal
823  {
824  inline VectorReference::VectorReference(const VectorBase &vector,
825  const size_type index)
826  : vector(vector)
827  , index(index)
828  {}
829 
830 
831  inline const VectorReference &
832  VectorReference::operator=(const VectorReference &r) const
833  {
834  // as explained in the class
835  // documentation, this is not the copy
836  // operator. so simply pass on to the
837  // "correct" assignment operator
838  *this = static_cast<PetscScalar>(r);
839 
840  return *this;
841  }
842 
843 
844 
845  inline VectorReference &
846  VectorReference::operator=(const VectorReference &r)
847  {
848  // as explained in the class
849  // documentation, this is not the copy
850  // operator. so simply pass on to the
851  // "correct" assignment operator
852  *this = static_cast<PetscScalar>(r);
853 
854  return *this;
855  }
856 
857 
858 
859  inline const VectorReference &
860  VectorReference::operator=(const PetscScalar &value) const
861  {
864  ExcWrongMode(VectorOperation::insert, vector.last_action));
865 
867 
868  const PetscInt petsc_i = index;
869 
870  const PetscErrorCode ierr =
871  VecSetValues(vector, 1, &petsc_i, &value, INSERT_VALUES);
872  AssertThrow(ierr == 0, ExcPETScError(ierr));
873 
875 
876  return *this;
877  }
878 
879 
880 
881  inline const VectorReference &
882  VectorReference::operator+=(const PetscScalar &value) const
883  {
886  ExcWrongMode(VectorOperation::add, vector.last_action));
887 
889 
891 
892  // we have to do above actions in any
893  // case to be consistent with the MPI
894  // communication model (see the
895  // comments in the documentation of
896  // PETScWrappers::MPI::Vector), but we
897  // can save some work if the addend is
898  // zero
899  if (value == PetscScalar())
900  return *this;
901 
902  // use the PETSc function to add something
903  const PetscInt petsc_i = index;
904  const PetscErrorCode ierr =
905  VecSetValues(vector, 1, &petsc_i, &value, ADD_VALUES);
906  AssertThrow(ierr == 0, ExcPETScError(ierr));
907 
908 
909  return *this;
910  }
911 
912 
913 
914  inline const VectorReference &
915  VectorReference::operator-=(const PetscScalar &value) const
916  {
919  ExcWrongMode(VectorOperation::add, vector.last_action));
920 
922 
924 
925  // we have to do above actions in any
926  // case to be consistent with the MPI
927  // communication model (see the
928  // comments in the documentation of
929  // PETScWrappers::MPI::Vector), but we
930  // can save some work if the addend is
931  // zero
932  if (value == PetscScalar())
933  return *this;
934 
935  // use the PETSc function to
936  // add something
937  const PetscInt petsc_i = index;
938  const PetscScalar subtractand = -value;
939  const PetscErrorCode ierr =
940  VecSetValues(vector, 1, &petsc_i, &subtractand, ADD_VALUES);
941  AssertThrow(ierr == 0, ExcPETScError(ierr));
942 
943  return *this;
944  }
945 
946 
947 
948  inline const VectorReference &
949  VectorReference::operator*=(const PetscScalar &value) const
950  {
953  ExcWrongMode(VectorOperation::insert, vector.last_action));
954 
956 
958 
959  // we have to do above actions in any
960  // case to be consistent with the MPI
961  // communication model (see the
962  // comments in the documentation of
963  // PETScWrappers::MPI::Vector), but we
964  // can save some work if the factor is
965  // one
966  if (value == 1.)
967  return *this;
968 
969  const PetscInt petsc_i = index;
970  const PetscScalar new_value = static_cast<PetscScalar>(*this) * value;
971 
972  const PetscErrorCode ierr =
973  VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
974  AssertThrow(ierr == 0, ExcPETScError(ierr));
975 
976  return *this;
977  }
978 
979 
980 
981  inline const VectorReference &
982  VectorReference::operator/=(const PetscScalar &value) const
983  {
986  ExcWrongMode(VectorOperation::insert, vector.last_action));
987 
989 
991 
992  // we have to do above actions in any
993  // case to be consistent with the MPI
994  // communication model (see the
995  // comments in the documentation of
996  // PETScWrappers::MPI::Vector), but we
997  // can save some work if the factor is
998  // one
999  if (value == 1.)
1000  return *this;
1001 
1002  const PetscInt petsc_i = index;
1003  const PetscScalar new_value = static_cast<PetscScalar>(*this) / value;
1004 
1005  const PetscErrorCode ierr =
1006  VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1007  AssertThrow(ierr == 0, ExcPETScError(ierr));
1008 
1009  return *this;
1010  }
1011 
1012 
1013 
1014  inline PetscReal
1015  VectorReference::real() const
1016  {
1017 # ifndef PETSC_USE_COMPLEX
1018  return static_cast<PetscScalar>(*this);
1019 # else
1020  return PetscRealPart(static_cast<PetscScalar>(*this));
1021 # endif
1022  }
1023 
1024 
1025 
1026  inline PetscReal
1027  VectorReference::imag() const
1028  {
1029 # ifndef PETSC_USE_COMPLEX
1030  return PetscReal(0);
1031 # else
1032  return PetscImaginaryPart(static_cast<PetscScalar>(*this));
1033 # endif
1034  }
1035 
1036  } // namespace internal
1037 
1038  inline bool
1039  VectorBase::in_local_range(const size_type index) const
1040  {
1041  PetscInt begin, end;
1042  const PetscErrorCode ierr =
1043  VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
1044  AssertThrow(ierr == 0, ExcPETScError(ierr));
1045 
1046  return ((index >= static_cast<size_type>(begin)) &&
1047  (index < static_cast<size_type>(end)));
1048  }
1049 
1050 
1051  inline IndexSet
1053  {
1054  IndexSet is(size());
1055 
1056  // PETSc only allows for contiguous local ranges, so this is simple
1057  const std::pair<size_type, size_type> x = local_range();
1058  is.add_range(x.first, x.second);
1059  return is;
1060  }
1061 
1062 
1063 
1064  inline bool
1066  {
1067  return ghosted;
1068  }
1069 
1070 
1071 
1072  inline void
1074  {}
1075 
1076 
1077 
1078  inline internal::VectorReference
1079  VectorBase::operator()(const size_type index)
1080  {
1081  return internal::VectorReference(*this, index);
1082  }
1083 
1084 
1085 
1086  inline PetscScalar
1087  VectorBase::operator()(const size_type index) const
1088  {
1089  return static_cast<PetscScalar>(internal::VectorReference(*this, index));
1090  }
1091 
1092 
1093 
1094  inline internal::VectorReference VectorBase::operator[](const size_type index)
1095  {
1096  return operator()(index);
1097  }
1098 
1099 
1100 
1101  inline PetscScalar VectorBase::operator[](const size_type index) const
1102  {
1103  return operator()(index);
1104  }
1105 
1106  inline const MPI_Comm &
1108  {
1109  static MPI_Comm comm;
1110  PetscObjectGetComm(reinterpret_cast<PetscObject>(vector), &comm);
1111  return comm;
1112  }
1113 
1114  inline void
1115  VectorBase::extract_subvector_to(const std::vector<size_type> &indices,
1116  std::vector<PetscScalar> & values) const
1117  {
1118  Assert(indices.size() <= values.size(),
1119  ExcDimensionMismatch(indices.size(), values.size()));
1120  extract_subvector_to(indices.begin(), indices.end(), values.begin());
1121  }
1122 
1123  template <typename ForwardIterator, typename OutputIterator>
1124  inline void
1125  VectorBase::extract_subvector_to(const ForwardIterator indices_begin,
1126  const ForwardIterator indices_end,
1127  OutputIterator values_begin) const
1128  {
1129  const PetscInt n_idx = static_cast<PetscInt>(indices_end - indices_begin);
1130  if (n_idx == 0)
1131  return;
1132 
1133  // if we are dealing
1134  // with a parallel vector
1135  if (ghosted)
1136  {
1137  // there is the possibility
1138  // that the vector has
1139  // ghost elements. in that
1140  // case, we first need to
1141  // figure out which
1142  // elements we own locally,
1143  // then get a pointer to
1144  // the elements that are
1145  // stored here (both the
1146  // ones we own as well as
1147  // the ghost elements). in
1148  // this array, the locally
1149  // owned elements come
1150  // first followed by the
1151  // ghost elements whose
1152  // position we can get from
1153  // an index set
1154  PetscInt begin, end;
1155  PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1156  AssertThrow(ierr == 0, ExcPETScError(ierr));
1157 
1158  Vec locally_stored_elements = nullptr;
1159  ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
1160  AssertThrow(ierr == 0, ExcPETScError(ierr));
1161 
1162  PetscInt lsize;
1163  ierr = VecGetSize(locally_stored_elements, &lsize);
1164  AssertThrow(ierr == 0, ExcPETScError(ierr));
1165 
1166  PetscScalar *ptr;
1167  ierr = VecGetArray(locally_stored_elements, &ptr);
1168  AssertThrow(ierr == 0, ExcPETScError(ierr));
1169 
1170  for (PetscInt i = 0; i < n_idx; ++i)
1171  {
1172  const unsigned int index = *(indices_begin + i);
1173  if (index >= static_cast<unsigned int>(begin) &&
1174  index < static_cast<unsigned int>(end))
1175  {
1176  // local entry
1177  *(values_begin + i) = *(ptr + index - begin);
1178  }
1179  else
1180  {
1181  // ghost entry
1182  const unsigned int ghostidx =
1183  ghost_indices.index_within_set(index);
1184 
1185  AssertIndexRange(ghostidx + end - begin, lsize);
1186  *(values_begin + i) = *(ptr + ghostidx + end - begin);
1187  }
1188  }
1189 
1190  ierr = VecRestoreArray(locally_stored_elements, &ptr);
1191  AssertThrow(ierr == 0, ExcPETScError(ierr));
1192 
1193  ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
1194  AssertThrow(ierr == 0, ExcPETScError(ierr));
1195  }
1196  // if the vector is local or the
1197  // caller, then simply access the
1198  // element we are interested in
1199  else
1200  {
1201  PetscInt begin, end;
1202  PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1203  AssertThrow(ierr == 0, ExcPETScError(ierr));
1204 
1205  PetscScalar *ptr;
1206  ierr = VecGetArray(vector, &ptr);
1207  AssertThrow(ierr == 0, ExcPETScError(ierr));
1208 
1209  for (PetscInt i = 0; i < n_idx; ++i)
1210  {
1211  const unsigned int index = *(indices_begin + i);
1212 
1213  Assert(index >= static_cast<unsigned int>(begin) &&
1214  index < static_cast<unsigned int>(end),
1215  ExcInternalError());
1216 
1217  *(values_begin + i) = *(ptr + index - begin);
1218  }
1219 
1220  ierr = VecRestoreArray(vector, &ptr);
1221  AssertThrow(ierr == 0, ExcPETScError(ierr));
1222  }
1223  }
1224 
1225 # endif // DOXYGEN
1226 } // namespace PETScWrappers
1227 
1229 
1230 # endif // DEAL_II_WITH_PETSC
1231 
1232 #endif
1233 /*---------------------------- petsc_vector_base.h --------------------------*/
reference operator[](const size_type index)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
types::global_dof_index size_type
Definition: cuda_kernels.h:45
void update_ghost_values() const
const internal::VectorReference const_reference
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
std::vector< value_type > l2_norm(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:837
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
void swap(VectorBase &u, VectorBase &v)
__global__ void add_and_dot(Number *res, Number *v1, const Number *v2, const Number *v3, const Number a, const size_type N)
static const char V
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool in_local_range(const size_type index) const
bool is_non_negative(const T &t)
std::string compress(const std::string &input)
Definition: utilities.cc:392
VectorOperation::values last_action
#define Assert(cond, exc)
Definition: exceptions.h:1419
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2732
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
VectorType::value_type * end(VectorType &V)
__global__ void equ(Number *val, const Number a, const Number *V_val, const size_type N)
__global__ void sadd(const Number s, Number *val, const Number a, const Number *V_val, const size_type N)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
static ::ExceptionBase & ExcGhostsPresent()
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1674
unsigned int global_dof_index
Definition: types.h:76
bool has_ghost_elements() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
VectorType::value_type * begin(VectorType &V)
reference operator()(const size_type index)
T min(const T &t, const MPI_Comm &mpi_communicator)
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2705
IndexSet locally_owned_elements() const
internal::VectorReference reference
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:564
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< PetscScalar > &values) const
virtual const MPI_Comm & get_mpi_communicator() const
static const bool value
T max(const T &t, const MPI_Comm &mpi_communicator)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()