Reference documentation for deal.II version Git 54f37b890c 2020-02-28 19:52:28 +0100
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
petsc_vector_base.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2019 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>
25 # include <deal.II/base/subscriptor.h>
26 
27 # include <deal.II/lac/exceptions.h>
28 # include <deal.II/lac/vector.h>
29 # include <deal.II/lac/vector_operation.h>
30 
31 # include <petscvec.h>
32 
33 # include <utility>
34 # include <vector>
35 
36 DEAL_II_NAMESPACE_OPEN
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:
87  using size_type = types::global_dof_index;
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;
250  using size_type = types::global_dof_index;
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 
703  DEAL_II_DEPRECATED
704  void
705  ratio(const VectorBase &a, const VectorBase &b);
706 
714  void
715  write_ascii(const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
716 
724  void
725  print(std::ostream & out,
726  const unsigned int precision = 3,
727  const bool scientific = true,
728  const bool across = true) const;
729 
742  void
743  swap(VectorBase &v);
744 
752  operator const Vec &() const;
753 
757  std::size_t
758  memory_consumption() const;
759 
764  virtual const MPI_Comm &
765  get_mpi_communicator() const;
766 
767  protected:
772  Vec vector;
773 
779  bool ghosted;
780 
787 
794 
795  // Make the reference class a friend.
796  friend class internal::VectorReference;
797 
804 
810  void
811  do_set_add_operation(const size_type n_elements,
812  const size_type * indices,
813  const PetscScalar *values,
814  const bool add_values);
815  };
816 
817 
818 
819  // ------------------- inline and template functions --------------
820 
829  inline void
831  {
832  u.swap(v);
833  }
834 
835 # ifndef DOXYGEN
836  namespace internal
837  {
838  inline VectorReference::VectorReference(const VectorBase &vector,
839  const size_type index)
840  : vector(vector)
841  , index(index)
842  {}
843 
844 
845  inline const VectorReference &
846  VectorReference::operator=(const VectorReference &r) const
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 VectorReference &
860  VectorReference::operator=(const VectorReference &r)
861  {
862  // as explained in the class
863  // documentation, this is not the copy
864  // operator. so simply pass on to the
865  // "correct" assignment operator
866  *this = static_cast<PetscScalar>(r);
867 
868  return *this;
869  }
870 
871 
872 
873  inline const VectorReference &
874  VectorReference::operator=(const PetscScalar &value) const
875  {
878  ExcWrongMode(VectorOperation::insert, vector.last_action));
879 
881 
882  const PetscInt petsc_i = index;
883 
884  const PetscErrorCode ierr =
885  VecSetValues(vector, 1, &petsc_i, &value, INSERT_VALUES);
886  AssertThrow(ierr == 0, ExcPETScError(ierr));
887 
889 
890  return *this;
891  }
892 
893 
894 
895  inline const VectorReference &
896  VectorReference::operator+=(const PetscScalar &value) const
897  {
900  ExcWrongMode(VectorOperation::add, vector.last_action));
901 
903 
905 
906  // we have to do above actions in any
907  // case to be consistent with the MPI
908  // communication model (see the
909  // comments in the documentation of
910  // PETScWrappers::MPI::Vector), but we
911  // can save some work if the addend is
912  // zero
913  if (value == PetscScalar())
914  return *this;
915 
916  // use the PETSc function to add something
917  const PetscInt petsc_i = index;
918  const PetscErrorCode ierr =
919  VecSetValues(vector, 1, &petsc_i, &value, ADD_VALUES);
920  AssertThrow(ierr == 0, ExcPETScError(ierr));
921 
922 
923  return *this;
924  }
925 
926 
927 
928  inline const VectorReference &
929  VectorReference::operator-=(const PetscScalar &value) const
930  {
933  ExcWrongMode(VectorOperation::add, vector.last_action));
934 
936 
938 
939  // we have to do above actions in any
940  // case to be consistent with the MPI
941  // communication model (see the
942  // comments in the documentation of
943  // PETScWrappers::MPI::Vector), but we
944  // can save some work if the addend is
945  // zero
946  if (value == PetscScalar())
947  return *this;
948 
949  // use the PETSc function to
950  // add something
951  const PetscInt petsc_i = index;
952  const PetscScalar subtractand = -value;
953  const PetscErrorCode ierr =
954  VecSetValues(vector, 1, &petsc_i, &subtractand, ADD_VALUES);
955  AssertThrow(ierr == 0, ExcPETScError(ierr));
956 
957  return *this;
958  }
959 
960 
961 
962  inline const VectorReference &
963  VectorReference::operator*=(const PetscScalar &value) const
964  {
967  ExcWrongMode(VectorOperation::insert, vector.last_action));
968 
970 
972 
973  // we have to do above actions in any
974  // case to be consistent with the MPI
975  // communication model (see the
976  // comments in the documentation of
977  // PETScWrappers::MPI::Vector), but we
978  // can save some work if the factor is
979  // one
980  if (value == 1.)
981  return *this;
982 
983  const PetscInt petsc_i = index;
984  const PetscScalar new_value = static_cast<PetscScalar>(*this) * value;
985 
986  const PetscErrorCode ierr =
987  VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
988  AssertThrow(ierr == 0, ExcPETScError(ierr));
989 
990  return *this;
991  }
992 
993 
994 
995  inline const VectorReference &
996  VectorReference::operator/=(const PetscScalar &value) const
997  {
1000  ExcWrongMode(VectorOperation::insert, vector.last_action));
1001 
1003 
1005 
1006  // we have to do above actions in any
1007  // case to be consistent with the MPI
1008  // communication model (see the
1009  // comments in the documentation of
1010  // PETScWrappers::MPI::Vector), but we
1011  // can save some work if the factor is
1012  // one
1013  if (value == 1.)
1014  return *this;
1015 
1016  const PetscInt petsc_i = index;
1017  const PetscScalar new_value = static_cast<PetscScalar>(*this) / value;
1018 
1019  const PetscErrorCode ierr =
1020  VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1021  AssertThrow(ierr == 0, ExcPETScError(ierr));
1022 
1023  return *this;
1024  }
1025 
1026 
1027 
1028  inline PetscReal
1029  VectorReference::real() const
1030  {
1031 # ifndef PETSC_USE_COMPLEX
1032  return static_cast<PetscScalar>(*this);
1033 # else
1034  return PetscRealPart(static_cast<PetscScalar>(*this));
1035 # endif
1036  }
1037 
1038 
1039 
1040  inline PetscReal
1041  VectorReference::imag() const
1042  {
1043 # ifndef PETSC_USE_COMPLEX
1044  return PetscReal(0);
1045 # else
1046  return PetscImaginaryPart(static_cast<PetscScalar>(*this));
1047 # endif
1048  }
1049 
1050  } // namespace internal
1051 
1052  inline bool
1053  VectorBase::in_local_range(const size_type index) const
1054  {
1055  PetscInt begin, end;
1056  const PetscErrorCode ierr =
1057  VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
1058  AssertThrow(ierr == 0, ExcPETScError(ierr));
1059 
1060  return ((index >= static_cast<size_type>(begin)) &&
1061  (index < static_cast<size_type>(end)));
1062  }
1063 
1064 
1065  inline IndexSet
1067  {
1068  IndexSet is(size());
1069 
1070  // PETSc only allows for contiguous local ranges, so this is simple
1071  const std::pair<size_type, size_type> x = local_range();
1072  is.add_range(x.first, x.second);
1073  return is;
1074  }
1075 
1076 
1077 
1078  inline bool
1080  {
1081  return ghosted;
1082  }
1083 
1084 
1085 
1086  inline void
1088  {}
1089 
1090 
1091 
1092  inline internal::VectorReference
1093  VectorBase::operator()(const size_type index)
1094  {
1095  return internal::VectorReference(*this, index);
1096  }
1097 
1098 
1099 
1100  inline PetscScalar
1101  VectorBase::operator()(const size_type index) const
1102  {
1103  return static_cast<PetscScalar>(internal::VectorReference(*this, index));
1104  }
1105 
1106 
1107 
1108  inline internal::VectorReference VectorBase::operator[](const size_type index)
1109  {
1110  return operator()(index);
1111  }
1112 
1113 
1114 
1115  inline PetscScalar VectorBase::operator[](const size_type index) const
1116  {
1117  return operator()(index);
1118  }
1119 
1120  inline const MPI_Comm &
1122  {
1123  static MPI_Comm comm;
1124  PetscObjectGetComm(reinterpret_cast<PetscObject>(vector), &comm);
1125  return comm;
1126  }
1127 
1128  inline void
1129  VectorBase::extract_subvector_to(const std::vector<size_type> &indices,
1130  std::vector<PetscScalar> & values) const
1131  {
1132  Assert(indices.size() <= values.size(),
1133  ExcDimensionMismatch(indices.size(), values.size()));
1134  extract_subvector_to(indices.begin(), indices.end(), values.begin());
1135  }
1136 
1137  template <typename ForwardIterator, typename OutputIterator>
1138  inline void
1139  VectorBase::extract_subvector_to(const ForwardIterator indices_begin,
1140  const ForwardIterator indices_end,
1141  OutputIterator values_begin) const
1142  {
1143  const PetscInt n_idx = static_cast<PetscInt>(indices_end - indices_begin);
1144  if (n_idx == 0)
1145  return;
1146 
1147  // if we are dealing
1148  // with a parallel vector
1149  if (ghosted)
1150  {
1151  // there is the possibility
1152  // that the vector has
1153  // ghost elements. in that
1154  // case, we first need to
1155  // figure out which
1156  // elements we own locally,
1157  // then get a pointer to
1158  // the elements that are
1159  // stored here (both the
1160  // ones we own as well as
1161  // the ghost elements). in
1162  // this array, the locally
1163  // owned elements come
1164  // first followed by the
1165  // ghost elements whose
1166  // position we can get from
1167  // an index set
1168  PetscInt begin, end;
1169  PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1170  AssertThrow(ierr == 0, ExcPETScError(ierr));
1171 
1172  Vec locally_stored_elements = nullptr;
1173  ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
1174  AssertThrow(ierr == 0, ExcPETScError(ierr));
1175 
1176  PetscInt lsize;
1177  ierr = VecGetSize(locally_stored_elements, &lsize);
1178  AssertThrow(ierr == 0, ExcPETScError(ierr));
1179 
1180  PetscScalar *ptr;
1181  ierr = VecGetArray(locally_stored_elements, &ptr);
1182  AssertThrow(ierr == 0, ExcPETScError(ierr));
1183 
1184  for (PetscInt i = 0; i < n_idx; ++i)
1185  {
1186  const unsigned int index = *(indices_begin + i);
1187  if (index >= static_cast<unsigned int>(begin) &&
1188  index < static_cast<unsigned int>(end))
1189  {
1190  // local entry
1191  *(values_begin + i) = *(ptr + index - begin);
1192  }
1193  else
1194  {
1195  // ghost entry
1196  const unsigned int ghostidx =
1197  ghost_indices.index_within_set(index);
1198 
1199  AssertIndexRange(ghostidx + end - begin, lsize);
1200  *(values_begin + i) = *(ptr + ghostidx + end - begin);
1201  }
1202  }
1203 
1204  ierr = VecRestoreArray(locally_stored_elements, &ptr);
1205  AssertThrow(ierr == 0, ExcPETScError(ierr));
1206 
1207  ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
1208  AssertThrow(ierr == 0, ExcPETScError(ierr));
1209  }
1210  // if the vector is local or the
1211  // caller, then simply access the
1212  // element we are interested in
1213  else
1214  {
1215  PetscInt begin, end;
1216  PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1217  AssertThrow(ierr == 0, ExcPETScError(ierr));
1218 
1219  PetscScalar *ptr;
1220  ierr = VecGetArray(vector, &ptr);
1221  AssertThrow(ierr == 0, ExcPETScError(ierr));
1222 
1223  for (PetscInt i = 0; i < n_idx; ++i)
1224  {
1225  const unsigned int index = *(indices_begin + i);
1226 
1227  Assert(index >= static_cast<unsigned int>(begin) &&
1228  index < static_cast<unsigned int>(end),
1229  ExcInternalError());
1230 
1231  *(values_begin + i) = *(ptr + index - begin);
1232  }
1233 
1234  ierr = VecRestoreArray(vector, &ptr);
1235  AssertThrow(ierr == 0, ExcPETScError(ierr));
1236  }
1237  }
1238 
1239 # endif // DOXYGEN
1240 } // namespace PETScWrappers
1241 
1242 DEAL_II_NAMESPACE_CLOSE
1243 
1244 # endif // DEAL_II_WITH_PETSC
1245 
1246 #endif
1247 /*---------------------------- petsc_vector_base.h --------------------------*/
reference operator[](const size_type index)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
void update_ghost_values() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
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)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
bool in_local_range(const size_type index) const
LinearAlgebra::distributed::Vector< Number > Vector
__global__ void scale(Number *val, const Number *V_val, const size_type N)
VectorOperation::values last_action
#define Assert(cond, exc)
Definition: exceptions.h:1411
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
constexpr SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
__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)
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:89
bool has_ghost_elements() const
reference operator()(const size_type index)
IndexSet locally_owned_elements() const
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:564
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< PetscScalar > &values) const
virtual const MPI_Comm & get_mpi_communicator() const
static ::ExceptionBase & ExcInternalError()