Reference documentation for deal.II version Git 5b8b897cb2 2021-04-22 22:24:19 -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 - 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 
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 constrcutor.
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 
465  void
466  set(const std::vector<size_type> & indices,
467  const std::vector<PetscScalar> &values);
468 
484  void
485  extract_subvector_to(const std::vector<size_type> &indices,
486  std::vector<PetscScalar> & values) const;
487 
515  template <typename ForwardIterator, typename OutputIterator>
516  void
517  extract_subvector_to(const ForwardIterator indices_begin,
518  const ForwardIterator indices_end,
519  OutputIterator values_begin) const;
520 
525  void
526  add(const std::vector<size_type> & indices,
527  const std::vector<PetscScalar> &values);
528 
533  void
534  add(const std::vector<size_type> & indices,
535  const ::Vector<PetscScalar> &values);
536 
542  void
543  add(const size_type n_elements,
544  const size_type * indices,
545  const PetscScalar *values);
546 
553  PetscScalar operator*(const VectorBase &vec) const;
554 
558  real_type
559  norm_sqr() const;
560 
564  PetscScalar
565  mean_value() const;
566 
574  real_type
575  l1_norm() const;
576 
581  real_type
582  l2_norm() const;
583 
588  real_type
589  lp_norm(const real_type p) const;
590 
595  real_type
596  linfty_norm() const;
597 
617  PetscScalar
618  add_and_dot(const PetscScalar a, const VectorBase &V, const VectorBase &W);
619 
628  real_type
629  min() const;
630 
639  real_type
640  max() const;
641 
647  bool
648  all_zero() const;
649 
659  bool
660  is_non_negative() const;
661 
665  VectorBase &
666  operator*=(const PetscScalar factor);
667 
671  VectorBase &
672  operator/=(const PetscScalar factor);
673 
677  VectorBase &
678  operator+=(const VectorBase &V);
679 
683  VectorBase &
684  operator-=(const VectorBase &V);
685 
690  void
691  add(const PetscScalar s);
692 
696  void
697  add(const PetscScalar a, const VectorBase &V);
698 
702  void
703  add(const PetscScalar a,
704  const VectorBase &V,
705  const PetscScalar b,
706  const VectorBase &W);
707 
711  void
712  sadd(const PetscScalar s, const VectorBase &V);
713 
717  void
718  sadd(const PetscScalar s, const PetscScalar a, const VectorBase &V);
719 
725  void
726  scale(const VectorBase &scaling_factors);
727 
731  void
732  equ(const PetscScalar a, const VectorBase &V);
733 
741  void
742  write_ascii(const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
743 
751  void
752  print(std::ostream & out,
753  const unsigned int precision = 3,
754  const bool scientific = true,
755  const bool across = true) const;
756 
769  void
770  swap(VectorBase &v);
771 
779  operator const Vec &() const;
780 
784  std::size_t
785  memory_consumption() const;
786 
791  virtual const MPI_Comm &
792  get_mpi_communicator() const;
793 
794  protected:
799  Vec vector;
800 
806  bool ghosted;
807 
814 
821 
822  // Make the reference class a friend.
823  friend class internal::VectorReference;
824 
831 
837  void
838  do_set_add_operation(const size_type n_elements,
839  const size_type * indices,
840  const PetscScalar *values,
841  const bool add_values);
842  };
843 
844 
845 
846  // ------------------- inline and template functions --------------
847 
855  inline void
857  {
858  u.swap(v);
859  }
860 
861 # ifndef DOXYGEN
862  namespace internal
863  {
864  inline VectorReference::VectorReference(const VectorBase &vector,
865  const size_type index)
866  : vector(vector)
867  , index(index)
868  {}
869 
870 
871  inline const VectorReference &
872  VectorReference::operator=(const VectorReference &r) const
873  {
874  // as explained in the class
875  // documentation, this is not the copy
876  // operator. so simply pass on to the
877  // "correct" assignment operator
878  *this = static_cast<PetscScalar>(r);
879 
880  return *this;
881  }
882 
883 
884 
885  inline VectorReference &
886  VectorReference::operator=(const VectorReference &r)
887  {
888  // as explained in the class
889  // documentation, this is not the copy
890  // operator. so simply pass on to the
891  // "correct" assignment operator
892  *this = static_cast<PetscScalar>(r);
893 
894  return *this;
895  }
896 
897 
898 
899  inline const VectorReference &
900  VectorReference::operator=(const PetscScalar &value) const
901  {
904  ExcWrongMode(VectorOperation::insert, vector.last_action));
905 
907 
908  const PetscInt petsc_i = index;
909 
910  const PetscErrorCode ierr =
911  VecSetValues(vector, 1, &petsc_i, &value, INSERT_VALUES);
912  AssertThrow(ierr == 0, ExcPETScError(ierr));
913 
915 
916  return *this;
917  }
918 
919 
920 
921  inline const VectorReference &
922  VectorReference::operator+=(const PetscScalar &value) const
923  {
926  ExcWrongMode(VectorOperation::add, vector.last_action));
927 
929 
931 
932  // we have to do above actions in any
933  // case to be consistent with the MPI
934  // communication model (see the
935  // comments in the documentation of
936  // PETScWrappers::MPI::Vector), but we
937  // can save some work if the addend is
938  // zero
939  if (value == PetscScalar())
940  return *this;
941 
942  // use the PETSc function to add something
943  const PetscInt petsc_i = index;
944  const PetscErrorCode ierr =
945  VecSetValues(vector, 1, &petsc_i, &value, ADD_VALUES);
946  AssertThrow(ierr == 0, ExcPETScError(ierr));
947 
948 
949  return *this;
950  }
951 
952 
953 
954  inline const VectorReference &
955  VectorReference::operator-=(const PetscScalar &value) const
956  {
959  ExcWrongMode(VectorOperation::add, vector.last_action));
960 
962 
964 
965  // we have to do above actions in any
966  // case to be consistent with the MPI
967  // communication model (see the
968  // comments in the documentation of
969  // PETScWrappers::MPI::Vector), but we
970  // can save some work if the addend is
971  // zero
972  if (value == PetscScalar())
973  return *this;
974 
975  // use the PETSc function to
976  // add something
977  const PetscInt petsc_i = index;
978  const PetscScalar subtractand = -value;
979  const PetscErrorCode ierr =
980  VecSetValues(vector, 1, &petsc_i, &subtractand, ADD_VALUES);
981  AssertThrow(ierr == 0, ExcPETScError(ierr));
982 
983  return *this;
984  }
985 
986 
987 
988  inline const VectorReference &
989  VectorReference::operator*=(const PetscScalar &value) const
990  {
993  ExcWrongMode(VectorOperation::insert, vector.last_action));
994 
996 
998 
999  // we have to do above actions in any
1000  // case to be consistent with the MPI
1001  // communication model (see the
1002  // comments in the documentation of
1003  // PETScWrappers::MPI::Vector), but we
1004  // can save some work if the factor is
1005  // one
1006  if (value == 1.)
1007  return *this;
1008 
1009  const PetscInt petsc_i = index;
1010  const PetscScalar new_value = static_cast<PetscScalar>(*this) * value;
1011 
1012  const PetscErrorCode ierr =
1013  VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1014  AssertThrow(ierr == 0, ExcPETScError(ierr));
1015 
1016  return *this;
1017  }
1018 
1019 
1020 
1021  inline const VectorReference &
1022  VectorReference::operator/=(const PetscScalar &value) const
1023  {
1026  ExcWrongMode(VectorOperation::insert, vector.last_action));
1027 
1029 
1031 
1032  // we have to do above actions in any
1033  // case to be consistent with the MPI
1034  // communication model (see the
1035  // comments in the documentation of
1036  // PETScWrappers::MPI::Vector), but we
1037  // can save some work if the factor is
1038  // one
1039  if (value == 1.)
1040  return *this;
1041 
1042  const PetscInt petsc_i = index;
1043  const PetscScalar new_value = static_cast<PetscScalar>(*this) / value;
1044 
1045  const PetscErrorCode ierr =
1046  VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1047  AssertThrow(ierr == 0, ExcPETScError(ierr));
1048 
1049  return *this;
1050  }
1051 
1052 
1053 
1054  inline PetscReal
1055  VectorReference::real() const
1056  {
1057 # ifndef PETSC_USE_COMPLEX
1058  return static_cast<PetscScalar>(*this);
1059 # else
1060  return PetscRealPart(static_cast<PetscScalar>(*this));
1061 # endif
1062  }
1063 
1064 
1065 
1066  inline PetscReal
1067  VectorReference::imag() const
1068  {
1069 # ifndef PETSC_USE_COMPLEX
1070  return PetscReal(0);
1071 # else
1072  return PetscImaginaryPart(static_cast<PetscScalar>(*this));
1073 # endif
1074  }
1075 
1076  } // namespace internal
1077 
1078  inline bool
1079  VectorBase::in_local_range(const size_type index) const
1080  {
1081  PetscInt begin, end;
1082  const PetscErrorCode ierr =
1083  VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
1084  AssertThrow(ierr == 0, ExcPETScError(ierr));
1085 
1086  return ((index >= static_cast<size_type>(begin)) &&
1087  (index < static_cast<size_type>(end)));
1088  }
1089 
1090 
1091  inline IndexSet
1093  {
1094  IndexSet is(size());
1095 
1096  // PETSc only allows for contiguous local ranges, so this is simple
1097  const std::pair<size_type, size_type> x = local_range();
1098  is.add_range(x.first, x.second);
1099  return is;
1100  }
1101 
1102 
1103 
1104  inline bool
1106  {
1107  return ghosted;
1108  }
1109 
1110 
1111 
1112  inline void
1114  {}
1115 
1116 
1117 
1118  inline internal::VectorReference
1119  VectorBase::operator()(const size_type index)
1120  {
1121  return internal::VectorReference(*this, index);
1122  }
1123 
1124 
1125 
1126  inline PetscScalar
1127  VectorBase::operator()(const size_type index) const
1128  {
1129  return static_cast<PetscScalar>(internal::VectorReference(*this, index));
1130  }
1131 
1132 
1133 
1134  inline internal::VectorReference VectorBase::operator[](const size_type index)
1135  {
1136  return operator()(index);
1137  }
1138 
1139 
1140 
1141  inline PetscScalar VectorBase::operator[](const size_type index) const
1142  {
1143  return operator()(index);
1144  }
1145 
1146  inline const MPI_Comm &
1148  {
1149  static MPI_Comm comm;
1150  PetscObjectGetComm(reinterpret_cast<PetscObject>(vector), &comm);
1151  return comm;
1152  }
1153 
1154  inline void
1155  VectorBase::extract_subvector_to(const std::vector<size_type> &indices,
1156  std::vector<PetscScalar> & values) const
1157  {
1158  Assert(indices.size() <= values.size(),
1159  ExcDimensionMismatch(indices.size(), values.size()));
1160  extract_subvector_to(indices.begin(), indices.end(), values.begin());
1161  }
1162 
1163  template <typename ForwardIterator, typename OutputIterator>
1164  inline void
1165  VectorBase::extract_subvector_to(const ForwardIterator indices_begin,
1166  const ForwardIterator indices_end,
1167  OutputIterator values_begin) const
1168  {
1169  const PetscInt n_idx = static_cast<PetscInt>(indices_end - indices_begin);
1170  if (n_idx == 0)
1171  return;
1172 
1173  // if we are dealing
1174  // with a parallel vector
1175  if (ghosted)
1176  {
1177  // there is the possibility
1178  // that the vector has
1179  // ghost elements. in that
1180  // case, we first need to
1181  // figure out which
1182  // elements we own locally,
1183  // then get a pointer to
1184  // the elements that are
1185  // stored here (both the
1186  // ones we own as well as
1187  // the ghost elements). in
1188  // this array, the locally
1189  // owned elements come
1190  // first followed by the
1191  // ghost elements whose
1192  // position we can get from
1193  // an index set
1194  PetscInt begin, end;
1195  PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1196  AssertThrow(ierr == 0, ExcPETScError(ierr));
1197 
1198  Vec locally_stored_elements = nullptr;
1199  ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
1200  AssertThrow(ierr == 0, ExcPETScError(ierr));
1201 
1202  PetscInt lsize;
1203  ierr = VecGetSize(locally_stored_elements, &lsize);
1204  AssertThrow(ierr == 0, ExcPETScError(ierr));
1205 
1206  PetscScalar *ptr;
1207  ierr = VecGetArray(locally_stored_elements, &ptr);
1208  AssertThrow(ierr == 0, ExcPETScError(ierr));
1209 
1210  for (PetscInt i = 0; i < n_idx; ++i)
1211  {
1212  const unsigned int index = *(indices_begin + i);
1213  if (index >= static_cast<unsigned int>(begin) &&
1214  index < static_cast<unsigned int>(end))
1215  {
1216  // local entry
1217  *(values_begin + i) = *(ptr + index - begin);
1218  }
1219  else
1220  {
1221  // ghost entry
1222  const unsigned int ghostidx =
1223  ghost_indices.index_within_set(index);
1224 
1225  AssertIndexRange(ghostidx + end - begin, lsize);
1226  *(values_begin + i) = *(ptr + ghostidx + end - begin);
1227  }
1228  }
1229 
1230  ierr = VecRestoreArray(locally_stored_elements, &ptr);
1231  AssertThrow(ierr == 0, ExcPETScError(ierr));
1232 
1233  ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
1234  AssertThrow(ierr == 0, ExcPETScError(ierr));
1235  }
1236  // if the vector is local or the
1237  // caller, then simply access the
1238  // element we are interested in
1239  else
1240  {
1241  PetscInt begin, end;
1242  PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1243  AssertThrow(ierr == 0, ExcPETScError(ierr));
1244 
1245  PetscScalar *ptr;
1246  ierr = VecGetArray(vector, &ptr);
1247  AssertThrow(ierr == 0, ExcPETScError(ierr));
1248 
1249  for (PetscInt i = 0; i < n_idx; ++i)
1250  {
1251  const unsigned int index = *(indices_begin + i);
1252 
1253  Assert(index >= static_cast<unsigned int>(begin) &&
1254  index < static_cast<unsigned int>(end),
1255  ExcInternalError());
1256 
1257  *(values_begin + i) = *(ptr + index - begin);
1258  }
1259 
1260  ierr = VecRestoreArray(vector, &ptr);
1261  AssertThrow(ierr == 0, ExcPETScError(ierr));
1262  }
1263  }
1264 
1265 # endif // DOXYGEN
1266 } // namespace PETScWrappers
1267 
1269 
1270 # endif // DEAL_II_WITH_PETSC
1271 
1272 #endif
1273 /*---------------------------- petsc_vector_base.h --------------------------*/
reference operator[](const size_type index)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
#define DEAL_II_DEPRECATED_EARLY
Definition: config.h:166
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:2043
#define AssertIndexRange(index, range)
Definition: exceptions.h:1698
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:1583
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:1473
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2817
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
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:395
VectorType::value_type * begin(VectorType &V)
reference operator()(const size_type index)
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2791
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:158
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()