Reference documentation for deal.II version Git 4abc4a1666 2020-07-04 19:58:34 +0200
\(\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;
176  DeclException3(ExcAccessToNonlocalElement,
177  int,
178  int,
179  int,
180  << "You tried to access element " << arg1
181  << " of a distributed vector, but only elements " << arg2
182  << " through " << arg3
183  << " are stored locally and can be accessed.");
187  DeclException2(ExcWrongMode,
188  int,
189  int,
190  << "You tried to do a "
191  << (arg1 == 1 ? "'set'" : (arg1 == 2 ? "'add'" : "???"))
192  << " operation but the vector is currently in "
193  << (arg2 == 1 ? "'set'" : (arg2 == 2 ? "'add'" : "???"))
194  << " mode. You first have to call 'compress()'.");
195 
196  private:
200  const VectorBase &vector;
201 
205  const size_type index;
206 
207  // Make the vector class a friend, so that it can create objects of the
208  // present type.
209  friend class ::PETScWrappers::VectorBase;
210  };
211  } // namespace internal
242  class VectorBase : public Subscriptor
243  {
244  public:
250  using value_type = PetscScalar;
251  using real_type = PetscReal;
253  using reference = internal::VectorReference;
254  using const_reference = const internal::VectorReference;
255 
260  VectorBase();
261 
266  VectorBase(const VectorBase &v);
267 
273  explicit VectorBase(const Vec &v);
274 
279  VectorBase &
280  operator=(const VectorBase &) = delete;
281 
285  virtual ~VectorBase() override;
286 
291  virtual void
292  clear();
293 
304  void
305  compress(const VectorOperation::values operation);
306 
320  VectorBase &
321  operator=(const PetscScalar s);
322 
328  bool
329  operator==(const VectorBase &v) const;
330 
336  bool
337  operator!=(const VectorBase &v) const;
338 
342  size_type
343  size() const;
344 
353  size_type
354  local_size() const;
355 
364  std::pair<size_type, size_type>
365  local_range() const;
366 
371  bool
372  in_local_range(const size_type index) const;
373 
387  IndexSet
388  locally_owned_elements() const;
389 
396  bool
397  has_ghost_elements() const;
398 
405  void
406  update_ghost_values() const;
407 
411  reference
412  operator()(const size_type index);
413 
417  PetscScalar
418  operator()(const size_type index) const;
419 
425  reference operator[](const size_type index);
426 
432  PetscScalar operator[](const size_type index) const;
433 
443  void
444  set(const std::vector<size_type> & indices,
445  const std::vector<PetscScalar> &values);
446 
462  void
463  extract_subvector_to(const std::vector<size_type> &indices,
464  std::vector<PetscScalar> & values) const;
465 
493  template <typename ForwardIterator, typename OutputIterator>
494  void
495  extract_subvector_to(const ForwardIterator indices_begin,
496  const ForwardIterator indices_end,
497  OutputIterator values_begin) const;
498 
503  void
504  add(const std::vector<size_type> & indices,
505  const std::vector<PetscScalar> &values);
506 
511  void
512  add(const std::vector<size_type> & indices,
513  const ::Vector<PetscScalar> &values);
514 
520  void
521  add(const size_type n_elements,
522  const size_type * indices,
523  const PetscScalar *values);
524 
531  PetscScalar operator*(const VectorBase &vec) const;
532 
536  real_type
537  norm_sqr() const;
538 
542  PetscScalar
543  mean_value() const;
544 
552  real_type
553  l1_norm() const;
554 
559  real_type
560  l2_norm() const;
561 
566  real_type
567  lp_norm(const real_type p) const;
568 
573  real_type
574  linfty_norm() const;
575 
595  PetscScalar
596  add_and_dot(const PetscScalar a, const VectorBase &V, const VectorBase &W);
597 
606  real_type
607  min() const;
608 
617  real_type
618  max() const;
619 
625  bool
626  all_zero() const;
627 
637  bool
638  is_non_negative() const;
639 
643  VectorBase &
644  operator*=(const PetscScalar factor);
645 
649  VectorBase &
650  operator/=(const PetscScalar factor);
651 
655  VectorBase &
656  operator+=(const VectorBase &V);
657 
661  VectorBase &
662  operator-=(const VectorBase &V);
663 
668  void
669  add(const PetscScalar s);
670 
674  void
675  add(const PetscScalar a, const VectorBase &V);
676 
680  void
681  add(const PetscScalar a,
682  const VectorBase &V,
683  const PetscScalar b,
684  const VectorBase &W);
685 
689  void
690  sadd(const PetscScalar s, const VectorBase &V);
691 
695  void
696  sadd(const PetscScalar s, const PetscScalar a, const VectorBase &V);
697 
703  void
704  scale(const VectorBase &scaling_factors);
705 
709  void
710  equ(const PetscScalar a, const VectorBase &V);
711 
719  void
720  write_ascii(const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
721 
729  void
730  print(std::ostream & out,
731  const unsigned int precision = 3,
732  const bool scientific = true,
733  const bool across = true) const;
734 
747  void
748  swap(VectorBase &v);
749 
757  operator const Vec &() const;
758 
762  std::size_t
763  memory_consumption() const;
764 
769  virtual const MPI_Comm &
770  get_mpi_communicator() const;
771 
772  protected:
777  Vec vector;
778 
784  bool ghosted;
785 
792 
799 
800  // Make the reference class a friend.
801  friend class internal::VectorReference;
802 
809 
815  void
816  do_set_add_operation(const size_type n_elements,
817  const size_type * indices,
818  const PetscScalar *values,
819  const bool add_values);
820  };
821 
822 
823 
824  // ------------------- inline and template functions --------------
825 
833  inline void
835  {
836  u.swap(v);
837  }
838 
839 # ifndef DOXYGEN
840  namespace internal
841  {
842  inline VectorReference::VectorReference(const VectorBase &vector,
843  const size_type index)
844  : vector(vector)
845  , index(index)
846  {}
847 
848 
849  inline const VectorReference &
850  VectorReference::operator=(const VectorReference &r) const
851  {
852  // as explained in the class
853  // documentation, this is not the copy
854  // operator. so simply pass on to the
855  // "correct" assignment operator
856  *this = static_cast<PetscScalar>(r);
857 
858  return *this;
859  }
860 
861 
862 
863  inline VectorReference &
864  VectorReference::operator=(const VectorReference &r)
865  {
866  // as explained in the class
867  // documentation, this is not the copy
868  // operator. so simply pass on to the
869  // "correct" assignment operator
870  *this = static_cast<PetscScalar>(r);
871 
872  return *this;
873  }
874 
875 
876 
877  inline const VectorReference &
878  VectorReference::operator=(const PetscScalar &value) const
879  {
882  ExcWrongMode(VectorOperation::insert, vector.last_action));
883 
885 
886  const PetscInt petsc_i = index;
887 
888  const PetscErrorCode ierr =
889  VecSetValues(vector, 1, &petsc_i, &value, INSERT_VALUES);
890  AssertThrow(ierr == 0, ExcPETScError(ierr));
891 
893 
894  return *this;
895  }
896 
897 
898 
899  inline const VectorReference &
900  VectorReference::operator+=(const PetscScalar &value) const
901  {
904  ExcWrongMode(VectorOperation::add, vector.last_action));
905 
907 
909 
910  // we have to do above actions in any
911  // case to be consistent with the MPI
912  // communication model (see the
913  // comments in the documentation of
914  // PETScWrappers::MPI::Vector), but we
915  // can save some work if the addend is
916  // zero
917  if (value == PetscScalar())
918  return *this;
919 
920  // use the PETSc function to add something
921  const PetscInt petsc_i = index;
922  const PetscErrorCode ierr =
923  VecSetValues(vector, 1, &petsc_i, &value, ADD_VALUES);
924  AssertThrow(ierr == 0, ExcPETScError(ierr));
925 
926 
927  return *this;
928  }
929 
930 
931 
932  inline const VectorReference &
933  VectorReference::operator-=(const PetscScalar &value) const
934  {
937  ExcWrongMode(VectorOperation::add, vector.last_action));
938 
940 
942 
943  // we have to do above actions in any
944  // case to be consistent with the MPI
945  // communication model (see the
946  // comments in the documentation of
947  // PETScWrappers::MPI::Vector), but we
948  // can save some work if the addend is
949  // zero
950  if (value == PetscScalar())
951  return *this;
952 
953  // use the PETSc function to
954  // add something
955  const PetscInt petsc_i = index;
956  const PetscScalar subtractand = -value;
957  const PetscErrorCode ierr =
958  VecSetValues(vector, 1, &petsc_i, &subtractand, ADD_VALUES);
959  AssertThrow(ierr == 0, ExcPETScError(ierr));
960 
961  return *this;
962  }
963 
964 
965 
966  inline const VectorReference &
967  VectorReference::operator*=(const PetscScalar &value) const
968  {
971  ExcWrongMode(VectorOperation::insert, vector.last_action));
972 
974 
976 
977  // we have to do above actions in any
978  // case to be consistent with the MPI
979  // communication model (see the
980  // comments in the documentation of
981  // PETScWrappers::MPI::Vector), but we
982  // can save some work if the factor is
983  // one
984  if (value == 1.)
985  return *this;
986 
987  const PetscInt petsc_i = index;
988  const PetscScalar new_value = static_cast<PetscScalar>(*this) * value;
989 
990  const PetscErrorCode ierr =
991  VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
992  AssertThrow(ierr == 0, ExcPETScError(ierr));
993 
994  return *this;
995  }
996 
997 
998 
999  inline const VectorReference &
1000  VectorReference::operator/=(const PetscScalar &value) const
1001  {
1004  ExcWrongMode(VectorOperation::insert, vector.last_action));
1005 
1007 
1009 
1010  // we have to do above actions in any
1011  // case to be consistent with the MPI
1012  // communication model (see the
1013  // comments in the documentation of
1014  // PETScWrappers::MPI::Vector), but we
1015  // can save some work if the factor is
1016  // one
1017  if (value == 1.)
1018  return *this;
1019 
1020  const PetscInt petsc_i = index;
1021  const PetscScalar new_value = static_cast<PetscScalar>(*this) / value;
1022 
1023  const PetscErrorCode ierr =
1024  VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1025  AssertThrow(ierr == 0, ExcPETScError(ierr));
1026 
1027  return *this;
1028  }
1029 
1030 
1031 
1032  inline PetscReal
1033  VectorReference::real() const
1034  {
1035 # ifndef PETSC_USE_COMPLEX
1036  return static_cast<PetscScalar>(*this);
1037 # else
1038  return PetscRealPart(static_cast<PetscScalar>(*this));
1039 # endif
1040  }
1041 
1042 
1043 
1044  inline PetscReal
1045  VectorReference::imag() const
1046  {
1047 # ifndef PETSC_USE_COMPLEX
1048  return PetscReal(0);
1049 # else
1050  return PetscImaginaryPart(static_cast<PetscScalar>(*this));
1051 # endif
1052  }
1053 
1054  } // namespace internal
1055 
1056  inline bool
1057  VectorBase::in_local_range(const size_type index) const
1058  {
1059  PetscInt begin, end;
1060  const PetscErrorCode ierr =
1061  VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
1062  AssertThrow(ierr == 0, ExcPETScError(ierr));
1063 
1064  return ((index >= static_cast<size_type>(begin)) &&
1065  (index < static_cast<size_type>(end)));
1066  }
1067 
1068 
1069  inline IndexSet
1071  {
1072  IndexSet is(size());
1073 
1074  // PETSc only allows for contiguous local ranges, so this is simple
1075  const std::pair<size_type, size_type> x = local_range();
1076  is.add_range(x.first, x.second);
1077  return is;
1078  }
1079 
1080 
1081 
1082  inline bool
1084  {
1085  return ghosted;
1086  }
1087 
1088 
1089 
1090  inline void
1092  {}
1093 
1094 
1095 
1096  inline internal::VectorReference
1097  VectorBase::operator()(const size_type index)
1098  {
1099  return internal::VectorReference(*this, index);
1100  }
1101 
1102 
1103 
1104  inline PetscScalar
1105  VectorBase::operator()(const size_type index) const
1106  {
1107  return static_cast<PetscScalar>(internal::VectorReference(*this, index));
1108  }
1109 
1110 
1111 
1112  inline internal::VectorReference VectorBase::operator[](const size_type index)
1113  {
1114  return operator()(index);
1115  }
1116 
1117 
1118 
1119  inline PetscScalar VectorBase::operator[](const size_type index) const
1120  {
1121  return operator()(index);
1122  }
1123 
1124  inline const MPI_Comm &
1126  {
1127  static MPI_Comm comm;
1128  PetscObjectGetComm(reinterpret_cast<PetscObject>(vector), &comm);
1129  return comm;
1130  }
1131 
1132  inline void
1133  VectorBase::extract_subvector_to(const std::vector<size_type> &indices,
1134  std::vector<PetscScalar> & values) const
1135  {
1136  Assert(indices.size() <= values.size(),
1137  ExcDimensionMismatch(indices.size(), values.size()));
1138  extract_subvector_to(indices.begin(), indices.end(), values.begin());
1139  }
1140 
1141  template <typename ForwardIterator, typename OutputIterator>
1142  inline void
1143  VectorBase::extract_subvector_to(const ForwardIterator indices_begin,
1144  const ForwardIterator indices_end,
1145  OutputIterator values_begin) const
1146  {
1147  const PetscInt n_idx = static_cast<PetscInt>(indices_end - indices_begin);
1148  if (n_idx == 0)
1149  return;
1150 
1151  // if we are dealing
1152  // with a parallel vector
1153  if (ghosted)
1154  {
1155  // there is the possibility
1156  // that the vector has
1157  // ghost elements. in that
1158  // case, we first need to
1159  // figure out which
1160  // elements we own locally,
1161  // then get a pointer to
1162  // the elements that are
1163  // stored here (both the
1164  // ones we own as well as
1165  // the ghost elements). in
1166  // this array, the locally
1167  // owned elements come
1168  // first followed by the
1169  // ghost elements whose
1170  // position we can get from
1171  // an index set
1172  PetscInt begin, end;
1173  PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1174  AssertThrow(ierr == 0, ExcPETScError(ierr));
1175 
1176  Vec locally_stored_elements = nullptr;
1177  ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
1178  AssertThrow(ierr == 0, ExcPETScError(ierr));
1179 
1180  PetscInt lsize;
1181  ierr = VecGetSize(locally_stored_elements, &lsize);
1182  AssertThrow(ierr == 0, ExcPETScError(ierr));
1183 
1184  PetscScalar *ptr;
1185  ierr = VecGetArray(locally_stored_elements, &ptr);
1186  AssertThrow(ierr == 0, ExcPETScError(ierr));
1187 
1188  for (PetscInt i = 0; i < n_idx; ++i)
1189  {
1190  const unsigned int index = *(indices_begin + i);
1191  if (index >= static_cast<unsigned int>(begin) &&
1192  index < static_cast<unsigned int>(end))
1193  {
1194  // local entry
1195  *(values_begin + i) = *(ptr + index - begin);
1196  }
1197  else
1198  {
1199  // ghost entry
1200  const unsigned int ghostidx =
1201  ghost_indices.index_within_set(index);
1202 
1203  AssertIndexRange(ghostidx + end - begin, lsize);
1204  *(values_begin + i) = *(ptr + ghostidx + end - begin);
1205  }
1206  }
1207 
1208  ierr = VecRestoreArray(locally_stored_elements, &ptr);
1209  AssertThrow(ierr == 0, ExcPETScError(ierr));
1210 
1211  ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
1212  AssertThrow(ierr == 0, ExcPETScError(ierr));
1213  }
1214  // if the vector is local or the
1215  // caller, then simply access the
1216  // element we are interested in
1217  else
1218  {
1219  PetscInt begin, end;
1220  PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1221  AssertThrow(ierr == 0, ExcPETScError(ierr));
1222 
1223  PetscScalar *ptr;
1224  ierr = VecGetArray(vector, &ptr);
1225  AssertThrow(ierr == 0, ExcPETScError(ierr));
1226 
1227  for (PetscInt i = 0; i < n_idx; ++i)
1228  {
1229  const unsigned int index = *(indices_begin + i);
1230 
1231  Assert(index >= static_cast<unsigned int>(begin) &&
1232  index < static_cast<unsigned int>(end),
1233  ExcInternalError());
1234 
1235  *(values_begin + i) = *(ptr + index - begin);
1236  }
1237 
1238  ierr = VecRestoreArray(vector, &ptr);
1239  AssertThrow(ierr == 0, ExcPETScError(ierr));
1240  }
1241  }
1242 
1243 # endif // DOXYGEN
1244 } // namespace PETScWrappers
1245 
1247 
1248 # endif // DEAL_II_WITH_PETSC
1249 
1250 #endif
1251 /*---------------------------- 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:839
#define AssertIndexRange(index, range)
Definition: exceptions.h:1628
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:1513
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:1403
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2779
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
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:1671
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:362
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:2753
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:152
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()