Reference documentation for deal.II version Git b3b9f019ec 2021-07-26 17:39:47 -0400
\(\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 - 2021 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 
55 namespace PETScWrappers
56 {
66  namespace internal
67  {
80  class VectorReference
81  {
82  public:
87 
88  private:
93  VectorReference(const VectorBase &vector, const size_type index);
94 
95  public:
96  /*
97  * Copy constructor.
98  */
99  VectorReference(const VectorReference &vector) = default;
100 
112  const VectorReference &
113  operator=(const VectorReference &r) const;
114 
120  VectorReference &
121  operator=(const VectorReference &r);
122 
126  const VectorReference &
127  operator=(const PetscScalar &s) const;
128 
132  const VectorReference &
133  operator+=(const PetscScalar &s) const;
134 
138  const VectorReference &
139  operator-=(const PetscScalar &s) const;
140 
144  const VectorReference &
145  operator*=(const PetscScalar &s) const;
146 
150  const VectorReference &
151  operator/=(const PetscScalar &s) const;
152 
156  PetscReal
157  real() const;
158 
165  PetscReal
166  imag() const;
167 
172  operator PetscScalar() const;
177  ExcAccessToNonlocalElement,
178  int,
179  int,
180  int,
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."
184  << "\n\n"
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.");
195  DeclException2(ExcWrongMode,
196  int,
197  int,
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()'.");
203 
204  private:
208  const VectorBase &vector;
209 
213  const size_type index;
214 
215  // Make the vector class a friend, so that it can create objects of the
216  // present type.
217  friend class ::PETScWrappers::VectorBase;
218  };
219  } // namespace internal
250  class VectorBase : public Subscriptor
251  {
252  public:
258  using value_type = PetscScalar;
259  using real_type = PetscReal;
261  using reference = internal::VectorReference;
262  using const_reference = const internal::VectorReference;
263 
268  VectorBase();
269 
274  VectorBase(const VectorBase &v);
275 
281  explicit VectorBase(const Vec &v);
282 
287  VectorBase &
288  operator=(const VectorBase &) = delete;
289 
293  virtual ~VectorBase() override;
294 
299  virtual void
300  clear();
301 
312  void
313  compress(const VectorOperation::values operation);
314 
328  VectorBase &
329  operator=(const PetscScalar s);
330 
336  bool
337  operator==(const VectorBase &v) const;
338 
344  bool
345  operator!=(const VectorBase &v) const;
346 
350  size_type
351  size() const;
352 
364  size_type
365  local_size() const;
366 
375  size_type
376  locally_owned_size() const;
377 
386  std::pair<size_type, size_type>
387  local_range() const;
388 
393  bool
394  in_local_range(const size_type index) const;
395 
409  IndexSet
410  locally_owned_elements() const;
411 
418  bool
419  has_ghost_elements() const;
420 
427  void
428  update_ghost_values() const;
429 
433  reference
434  operator()(const size_type index);
435 
439  PetscScalar
440  operator()(const size_type index) const;
441 
447  reference operator[](const size_type index);
448 
454  PetscScalar operator[](const size_type index) const;
455 
462  void
463  set(const std::vector<size_type> & indices,
464  const std::vector<PetscScalar> &values);
465 
481  void
482  extract_subvector_to(const std::vector<size_type> &indices,
483  std::vector<PetscScalar> & values) const;
484 
512  template <typename ForwardIterator, typename OutputIterator>
513  void
514  extract_subvector_to(const ForwardIterator indices_begin,
515  const ForwardIterator indices_end,
516  OutputIterator values_begin) const;
517 
522  void
523  add(const std::vector<size_type> & indices,
524  const std::vector<PetscScalar> &values);
525 
530  void
531  add(const std::vector<size_type> & indices,
532  const ::Vector<PetscScalar> &values);
533 
539  void
540  add(const size_type n_elements,
541  const size_type * indices,
542  const PetscScalar *values);
543 
550  PetscScalar operator*(const VectorBase &vec) const;
551 
555  real_type
556  norm_sqr() const;
557 
561  PetscScalar
562  mean_value() const;
563 
571  real_type
572  l1_norm() const;
573 
578  real_type
579  l2_norm() const;
580 
585  real_type
586  lp_norm(const real_type p) const;
587 
592  real_type
593  linfty_norm() const;
594 
614  PetscScalar
615  add_and_dot(const PetscScalar a, const VectorBase &V, const VectorBase &W);
616 
625  real_type
626  min() const;
627 
636  real_type
637  max() const;
638 
644  bool
645  all_zero() const;
646 
656  bool
657  is_non_negative() const;
658 
662  VectorBase &
663  operator*=(const PetscScalar factor);
664 
668  VectorBase &
669  operator/=(const PetscScalar factor);
670 
674  VectorBase &
675  operator+=(const VectorBase &V);
676 
680  VectorBase &
681  operator-=(const VectorBase &V);
682 
687  void
688  add(const PetscScalar s);
689 
693  void
694  add(const PetscScalar a, const VectorBase &V);
695 
699  void
700  add(const PetscScalar a,
701  const VectorBase &V,
702  const PetscScalar b,
703  const VectorBase &W);
704 
708  void
709  sadd(const PetscScalar s, const VectorBase &V);
710 
714  void
715  sadd(const PetscScalar s, const PetscScalar a, const VectorBase &V);
716 
722  void
723  scale(const VectorBase &scaling_factors);
724 
728  void
729  equ(const PetscScalar a, const VectorBase &V);
730 
738  void
739  write_ascii(const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
740 
748  void
749  print(std::ostream & out,
750  const unsigned int precision = 3,
751  const bool scientific = true,
752  const bool across = true) const;
753 
766  void
767  swap(VectorBase &v);
768 
776  operator const Vec &() const;
777 
781  std::size_t
782  memory_consumption() const;
783 
788  virtual const MPI_Comm &
789  get_mpi_communicator() const;
790 
791  protected:
796  Vec vector;
797 
803  bool ghosted;
804 
811 
818 
819  // Make the reference class a friend.
820  friend class internal::VectorReference;
821 
828 
834  void
835  do_set_add_operation(const size_type n_elements,
836  const size_type * indices,
837  const PetscScalar *values,
838  const bool add_values);
839  };
840 
841 
842 
843  // ------------------- inline and template functions --------------
844 
852  inline void
854  {
855  u.swap(v);
856  }
857 
858 # ifndef DOXYGEN
859  namespace internal
860  {
861  inline VectorReference::VectorReference(const VectorBase &vector,
862  const size_type index)
863  : vector(vector)
864  , index(index)
865  {}
866 
867 
868  inline const VectorReference &
869  VectorReference::operator=(const VectorReference &r) const
870  {
871  // as explained in the class
872  // documentation, this is not the copy
873  // operator. so simply pass on to the
874  // "correct" assignment operator
875  *this = static_cast<PetscScalar>(r);
876 
877  return *this;
878  }
879 
880 
881 
882  inline VectorReference &
883  VectorReference::operator=(const VectorReference &r)
884  {
885  // as explained in the class
886  // documentation, this is not the copy
887  // operator. so simply pass on to the
888  // "correct" assignment operator
889  *this = static_cast<PetscScalar>(r);
890 
891  return *this;
892  }
893 
894 
895 
896  inline const VectorReference &
897  VectorReference::operator=(const PetscScalar &value) const
898  {
901  ExcWrongMode(VectorOperation::insert, vector.last_action));
902 
904 
905  const PetscInt petsc_i = index;
906 
907  const PetscErrorCode ierr =
908  VecSetValues(vector, 1, &petsc_i, &value, INSERT_VALUES);
909  AssertThrow(ierr == 0, ExcPETScError(ierr));
910 
912 
913  return *this;
914  }
915 
916 
917 
918  inline const VectorReference &
919  VectorReference::operator+=(const PetscScalar &value) const
920  {
923  ExcWrongMode(VectorOperation::add, vector.last_action));
924 
926 
928 
929  // we have to do above actions in any
930  // case to be consistent with the MPI
931  // communication model (see the
932  // comments in the documentation of
933  // PETScWrappers::MPI::Vector), but we
934  // can save some work if the addend is
935  // zero
936  if (value == PetscScalar())
937  return *this;
938 
939  // use the PETSc function to add something
940  const PetscInt petsc_i = index;
941  const PetscErrorCode ierr =
942  VecSetValues(vector, 1, &petsc_i, &value, ADD_VALUES);
943  AssertThrow(ierr == 0, ExcPETScError(ierr));
944 
945 
946  return *this;
947  }
948 
949 
950 
951  inline const VectorReference &
952  VectorReference::operator-=(const PetscScalar &value) const
953  {
956  ExcWrongMode(VectorOperation::add, vector.last_action));
957 
959 
961 
962  // we have to do above actions in any
963  // case to be consistent with the MPI
964  // communication model (see the
965  // comments in the documentation of
966  // PETScWrappers::MPI::Vector), but we
967  // can save some work if the addend is
968  // zero
969  if (value == PetscScalar())
970  return *this;
971 
972  // use the PETSc function to
973  // add something
974  const PetscInt petsc_i = index;
975  const PetscScalar subtractand = -value;
976  const PetscErrorCode ierr =
977  VecSetValues(vector, 1, &petsc_i, &subtractand, ADD_VALUES);
978  AssertThrow(ierr == 0, ExcPETScError(ierr));
979 
980  return *this;
981  }
982 
983 
984 
985  inline const VectorReference &
986  VectorReference::operator*=(const PetscScalar &value) const
987  {
990  ExcWrongMode(VectorOperation::insert, vector.last_action));
991 
993 
995 
996  // we have to do above actions in any
997  // case to be consistent with the MPI
998  // communication model (see the
999  // comments in the documentation of
1000  // PETScWrappers::MPI::Vector), but we
1001  // can save some work if the factor is
1002  // one
1003  if (value == 1.)
1004  return *this;
1005 
1006  const PetscInt petsc_i = index;
1007  const PetscScalar new_value = static_cast<PetscScalar>(*this) * value;
1008 
1009  const PetscErrorCode ierr =
1010  VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1011  AssertThrow(ierr == 0, ExcPETScError(ierr));
1012 
1013  return *this;
1014  }
1015 
1016 
1017 
1018  inline const VectorReference &
1019  VectorReference::operator/=(const PetscScalar &value) const
1020  {
1023  ExcWrongMode(VectorOperation::insert, vector.last_action));
1024 
1026 
1028 
1029  // we have to do above actions in any
1030  // case to be consistent with the MPI
1031  // communication model (see the
1032  // comments in the documentation of
1033  // PETScWrappers::MPI::Vector), but we
1034  // can save some work if the factor is
1035  // one
1036  if (value == 1.)
1037  return *this;
1038 
1039  const PetscInt petsc_i = index;
1040  const PetscScalar new_value = static_cast<PetscScalar>(*this) / value;
1041 
1042  const PetscErrorCode ierr =
1043  VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1044  AssertThrow(ierr == 0, ExcPETScError(ierr));
1045 
1046  return *this;
1047  }
1048 
1049 
1050 
1051  inline PetscReal
1052  VectorReference::real() const
1053  {
1054 # ifndef PETSC_USE_COMPLEX
1055  return static_cast<PetscScalar>(*this);
1056 # else
1057  return PetscRealPart(static_cast<PetscScalar>(*this));
1058 # endif
1059  }
1060 
1061 
1062 
1063  inline PetscReal
1064  VectorReference::imag() const
1065  {
1066 # ifndef PETSC_USE_COMPLEX
1067  return PetscReal(0);
1068 # else
1069  return PetscImaginaryPart(static_cast<PetscScalar>(*this));
1070 # endif
1071  }
1072 
1073  } // namespace internal
1074 
1075  inline bool
1076  VectorBase::in_local_range(const size_type index) const
1077  {
1078  PetscInt begin, end;
1079  const PetscErrorCode ierr =
1080  VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
1081  AssertThrow(ierr == 0, ExcPETScError(ierr));
1082 
1083  return ((index >= static_cast<size_type>(begin)) &&
1084  (index < static_cast<size_type>(end)));
1085  }
1086 
1087 
1088  inline IndexSet
1090  {
1091  IndexSet is(size());
1092 
1093  // PETSc only allows for contiguous local ranges, so this is simple
1094  const std::pair<size_type, size_type> x = local_range();
1095  is.add_range(x.first, x.second);
1096  return is;
1097  }
1098 
1099 
1100 
1101  inline bool
1103  {
1104  return ghosted;
1105  }
1106 
1107 
1108 
1109  inline void
1111  {}
1112 
1113 
1114 
1115  inline internal::VectorReference
1116  VectorBase::operator()(const size_type index)
1117  {
1118  return internal::VectorReference(*this, index);
1119  }
1120 
1121 
1122 
1123  inline PetscScalar
1124  VectorBase::operator()(const size_type index) const
1125  {
1126  return static_cast<PetscScalar>(internal::VectorReference(*this, index));
1127  }
1128 
1129 
1130 
1131  inline internal::VectorReference VectorBase::operator[](const size_type index)
1132  {
1133  return operator()(index);
1134  }
1135 
1136 
1137 
1138  inline PetscScalar VectorBase::operator[](const size_type index) const
1139  {
1140  return operator()(index);
1141  }
1142 
1143  inline const MPI_Comm &
1145  {
1146  static MPI_Comm comm;
1147  PetscObjectGetComm(reinterpret_cast<PetscObject>(vector), &comm);
1148  return comm;
1149  }
1150 
1151  inline void
1152  VectorBase::extract_subvector_to(const std::vector<size_type> &indices,
1153  std::vector<PetscScalar> & values) const
1154  {
1155  Assert(indices.size() <= values.size(),
1156  ExcDimensionMismatch(indices.size(), values.size()));
1157  extract_subvector_to(indices.begin(), indices.end(), values.begin());
1158  }
1159 
1160  template <typename ForwardIterator, typename OutputIterator>
1161  inline void
1162  VectorBase::extract_subvector_to(const ForwardIterator indices_begin,
1163  const ForwardIterator indices_end,
1164  OutputIterator values_begin) const
1165  {
1166  const PetscInt n_idx = static_cast<PetscInt>(indices_end - indices_begin);
1167  if (n_idx == 0)
1168  return;
1169 
1170  // if we are dealing
1171  // with a parallel vector
1172  if (ghosted)
1173  {
1174  // there is the possibility
1175  // that the vector has
1176  // ghost elements. in that
1177  // case, we first need to
1178  // figure out which
1179  // elements we own locally,
1180  // then get a pointer to
1181  // the elements that are
1182  // stored here (both the
1183  // ones we own as well as
1184  // the ghost elements). in
1185  // this array, the locally
1186  // owned elements come
1187  // first followed by the
1188  // ghost elements whose
1189  // position we can get from
1190  // an index set
1191  PetscInt begin, end;
1192  PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1193  AssertThrow(ierr == 0, ExcPETScError(ierr));
1194 
1195  Vec locally_stored_elements = nullptr;
1196  ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
1197  AssertThrow(ierr == 0, ExcPETScError(ierr));
1198 
1199  PetscInt lsize;
1200  ierr = VecGetSize(locally_stored_elements, &lsize);
1201  AssertThrow(ierr == 0, ExcPETScError(ierr));
1202 
1203  PetscScalar *ptr;
1204  ierr = VecGetArray(locally_stored_elements, &ptr);
1205  AssertThrow(ierr == 0, ExcPETScError(ierr));
1206 
1207  for (PetscInt i = 0; i < n_idx; ++i)
1208  {
1209  const unsigned int index = *(indices_begin + i);
1210  if (index >= static_cast<unsigned int>(begin) &&
1211  index < static_cast<unsigned int>(end))
1212  {
1213  // local entry
1214  *(values_begin + i) = *(ptr + index - begin);
1215  }
1216  else
1217  {
1218  // ghost entry
1219  const unsigned int ghostidx =
1220  ghost_indices.index_within_set(index);
1221 
1222  AssertIndexRange(ghostidx + end - begin, lsize);
1223  *(values_begin + i) = *(ptr + ghostidx + end - begin);
1224  }
1225  }
1226 
1227  ierr = VecRestoreArray(locally_stored_elements, &ptr);
1228  AssertThrow(ierr == 0, ExcPETScError(ierr));
1229 
1230  ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
1231  AssertThrow(ierr == 0, ExcPETScError(ierr));
1232  }
1233  // if the vector is local or the
1234  // caller, then simply access the
1235  // element we are interested in
1236  else
1237  {
1238  PetscInt begin, end;
1239  PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1240  AssertThrow(ierr == 0, ExcPETScError(ierr));
1241 
1242  PetscScalar *ptr;
1243  ierr = VecGetArray(vector, &ptr);
1244  AssertThrow(ierr == 0, ExcPETScError(ierr));
1245 
1246  for (PetscInt i = 0; i < n_idx; ++i)
1247  {
1248  const unsigned int index = *(indices_begin + i);
1249 
1250  Assert(index >= static_cast<unsigned int>(begin) &&
1251  index < static_cast<unsigned int>(end),
1252  ExcInternalError());
1253 
1254  *(values_begin + i) = *(ptr + index - begin);
1255  }
1256 
1257  ierr = VecRestoreArray(vector, &ptr);
1258  AssertThrow(ierr == 0, ExcPETScError(ierr));
1259  }
1260  }
1261 
1262 # endif // DOXYGEN
1263 } // namespace PETScWrappers
1264 
1266 
1267 # endif // DEAL_II_WITH_PETSC
1268 
1269 #endif
1270 /*---------------------------- petsc_vector_base.h --------------------------*/
reference operator[](const size_type index)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
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:2040
#define AssertIndexRange(index, range)
Definition: exceptions.h:1722
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:1575
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:1465
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2939
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
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:1673
unsigned int global_dof_index
Definition: types.h:76
*braid_SplitCommworld & comm
bool has_ghost_elements() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
VectorType::value_type * begin(VectorType &V)
reference operator()(const size_type index)
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2913
IndexSet locally_owned_elements() const
internal::VectorReference reference
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:561
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
#define DEAL_II_DEPRECATED
Definition: config.h:160
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()