Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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 template <typename number>
40 class Vector;
41 
42 
50 namespace PETScWrappers
51 {
52  // forward declaration
53  class VectorBase;
54 
64  namespace internal
65  {
78  class VectorReference
79  {
80  public:
84  using size_type = types::global_dof_index;
85 
86  private:
91  VectorReference(const VectorBase &vector, const size_type index);
92 
93 
94  public:
106  const VectorReference &
107  operator=(const VectorReference &r) const;
108 
114  VectorReference &
115  operator=(const VectorReference &r);
116 
120  const VectorReference &
121  operator=(const PetscScalar &s) const;
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  PetscReal
151  real() const;
152 
159  PetscReal
160  imag() const;
161 
166  operator PetscScalar() const;
170  DeclException3(ExcAccessToNonlocalElement,
171  int,
172  int,
173  int,
174  << "You tried to access element " << arg1
175  << " of a distributed vector, but only elements " << arg2
176  << " through " << arg3
177  << " are stored locally and can be accessed.");
181  DeclException2(ExcWrongMode,
182  int,
183  int,
184  << "You tried to do a "
185  << (arg1 == 1 ? "'set'" : (arg1 == 2 ? "'add'" : "???"))
186  << " operation but the vector is currently in "
187  << (arg2 == 1 ? "'set'" : (arg2 == 2 ? "'add'" : "???"))
188  << " mode. You first have to call 'compress()'.");
189 
190  private:
194  const VectorBase &vector;
195 
199  const size_type index;
200 
205  friend class ::PETScWrappers::VectorBase;
206  };
207  } // namespace internal
239  class VectorBase : public Subscriptor
240  {
241  public:
247  using value_type = PetscScalar;
248  using real_type = PetscReal;
249  using size_type = types::global_dof_index;
250  using reference = internal::VectorReference;
251  using const_reference = const internal::VectorReference;
252 
257  VectorBase();
258 
263  VectorBase(const VectorBase &v);
264 
270  explicit VectorBase(const Vec &v);
271 
276  VectorBase &
277  operator=(const VectorBase &) = delete;
278 
282  virtual ~VectorBase() override;
283 
288  virtual void
289  clear();
290 
301  void
302  compress(const VectorOperation::values operation);
303 
317  VectorBase &
318  operator=(const PetscScalar s);
319 
325  bool
326  operator==(const VectorBase &v) const;
327 
333  bool
334  operator!=(const VectorBase &v) const;
335 
339  size_type
340  size() const;
341 
350  size_type
351  local_size() const;
352 
361  std::pair<size_type, size_type>
362  local_range() const;
363 
368  bool
369  in_local_range(const size_type index) const;
370 
384  IndexSet
385  locally_owned_elements() const;
386 
393  bool
394  has_ghost_elements() const;
395 
402  void
403  update_ghost_values() const;
404 
408  reference
409  operator()(const size_type index);
410 
414  PetscScalar
415  operator()(const size_type index) const;
416 
422  reference operator[](const size_type index);
423 
429  PetscScalar operator[](const size_type index) const;
430 
437  void
438  set(const std::vector<size_type> & indices,
439  const std::vector<PetscScalar> &values);
440 
456  void
457  extract_subvector_to(const std::vector<size_type> &indices,
458  std::vector<PetscScalar> & values) const;
459 
487  template <typename ForwardIterator, typename OutputIterator>
488  void
489  extract_subvector_to(const ForwardIterator indices_begin,
490  const ForwardIterator indices_end,
491  OutputIterator values_begin) const;
492 
497  void
498  add(const std::vector<size_type> & indices,
499  const std::vector<PetscScalar> &values);
500 
505  void
506  add(const std::vector<size_type> & indices,
507  const ::Vector<PetscScalar> &values);
508 
514  void
515  add(const size_type n_elements,
516  const size_type * indices,
517  const PetscScalar *values);
518 
525  PetscScalar operator*(const VectorBase &vec) const;
526 
530  real_type
531  norm_sqr() const;
532 
536  PetscScalar
537  mean_value() const;
538 
546  real_type
547  l1_norm() const;
548 
553  real_type
554  l2_norm() const;
555 
560  real_type
561  lp_norm(const real_type p) const;
562 
567  real_type
568  linfty_norm() const;
569 
589  PetscScalar
590  add_and_dot(const PetscScalar a, const VectorBase &V, const VectorBase &W);
591 
595  real_type
596  min() const;
597 
601  real_type
602  max() const;
603 
609  bool
610  all_zero() const;
611 
617  bool
618  is_non_negative() const;
619 
623  VectorBase &
624  operator*=(const PetscScalar factor);
625 
629  VectorBase &
630  operator/=(const PetscScalar factor);
631 
635  VectorBase &
636  operator+=(const VectorBase &V);
637 
641  VectorBase &
642  operator-=(const VectorBase &V);
643 
648  void
649  add(const PetscScalar s);
650 
654  void
655  add(const PetscScalar a, const VectorBase &V);
656 
660  void
661  add(const PetscScalar a,
662  const VectorBase &V,
663  const PetscScalar b,
664  const VectorBase &W);
665 
669  void
670  sadd(const PetscScalar s, const VectorBase &V);
671 
675  void
676  sadd(const PetscScalar s, const PetscScalar a, const VectorBase &V);
677 
683  void
684  scale(const VectorBase &scaling_factors);
685 
689  void
690  equ(const PetscScalar a, const VectorBase &V);
691 
702  DEAL_II_DEPRECATED
703  void
704  ratio(const VectorBase &a, const VectorBase &b);
705 
713  void
714  write_ascii(const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
715 
723  void
724  print(std::ostream & out,
725  const unsigned int precision = 3,
726  const bool scientific = true,
727  const bool across = true) const;
728 
741  void
742  swap(VectorBase &v);
743 
751  operator const Vec &() const;
752 
756  std::size_t
757  memory_consumption() const;
758 
763  virtual const MPI_Comm &
764  get_mpi_communicator() const;
765 
766  protected:
771  Vec vector;
772 
778  bool ghosted;
779 
786 
793 
797  friend class internal::VectorReference;
798 
805 
811  void
812  do_set_add_operation(const size_type n_elements,
813  const size_type * indices,
814  const PetscScalar *values,
815  const bool add_values);
816  };
817 
818 
819 
820  // ------------------- inline and template functions --------------
821 
830  inline void
832  {
833  u.swap(v);
834  }
835 
836 # ifndef DOXYGEN
837  namespace internal
838  {
839  inline VectorReference::VectorReference(const VectorBase &vector,
840  const size_type index)
841  : vector(vector)
842  , index(index)
843  {}
844 
845 
846  inline const VectorReference &
847  VectorReference::operator=(const VectorReference &r) const
848  {
849  // as explained in the class
850  // documentation, this is not the copy
851  // operator. so simply pass on to the
852  // "correct" assignment operator
853  *this = static_cast<PetscScalar>(r);
854 
855  return *this;
856  }
857 
858 
859 
860  inline VectorReference &
861  VectorReference::operator=(const VectorReference &r)
862  {
863  // as explained in the class
864  // documentation, this is not the copy
865  // operator. so simply pass on to the
866  // "correct" assignment operator
867  *this = static_cast<PetscScalar>(r);
868 
869  return *this;
870  }
871 
872 
873 
874  inline const VectorReference &
875  VectorReference::operator=(const PetscScalar &value) const
876  {
879  ExcWrongMode(VectorOperation::insert, vector.last_action));
880 
882 
883  const PetscInt petsc_i = index;
884 
885  const PetscErrorCode ierr =
886  VecSetValues(vector, 1, &petsc_i, &value, INSERT_VALUES);
887  AssertThrow(ierr == 0, ExcPETScError(ierr));
888 
890 
891  return *this;
892  }
893 
894 
895 
896  inline const VectorReference &
897  VectorReference::operator+=(const PetscScalar &value) const
898  {
901  ExcWrongMode(VectorOperation::add, vector.last_action));
902 
904 
906 
907  // we have to do above actions in any
908  // case to be consistent with the MPI
909  // communication model (see the
910  // comments in the documentation of
911  // PETScWrappers::MPI::Vector), but we
912  // can save some work if the addend is
913  // zero
914  if (value == PetscScalar())
915  return *this;
916 
917  // use the PETSc function to add something
918  const PetscInt petsc_i = index;
919  const PetscErrorCode ierr =
920  VecSetValues(vector, 1, &petsc_i, &value, ADD_VALUES);
921  AssertThrow(ierr == 0, ExcPETScError(ierr));
922 
923 
924  return *this;
925  }
926 
927 
928 
929  inline const VectorReference &
930  VectorReference::operator-=(const PetscScalar &value) const
931  {
934  ExcWrongMode(VectorOperation::add, vector.last_action));
935 
937 
939 
940  // we have to do above actions in any
941  // case to be consistent with the MPI
942  // communication model (see the
943  // comments in the documentation of
944  // PETScWrappers::MPI::Vector), but we
945  // can save some work if the addend is
946  // zero
947  if (value == PetscScalar())
948  return *this;
949 
950  // use the PETSc function to
951  // add something
952  const PetscInt petsc_i = index;
953  const PetscScalar subtractand = -value;
954  const PetscErrorCode ierr =
955  VecSetValues(vector, 1, &petsc_i, &subtractand, ADD_VALUES);
956  AssertThrow(ierr == 0, ExcPETScError(ierr));
957 
958  return *this;
959  }
960 
961 
962 
963  inline const VectorReference &
964  VectorReference::operator*=(const PetscScalar &value) const
965  {
968  ExcWrongMode(VectorOperation::insert, vector.last_action));
969 
971 
973 
974  // we have to do above actions in any
975  // case to be consistent with the MPI
976  // communication model (see the
977  // comments in the documentation of
978  // PETScWrappers::MPI::Vector), but we
979  // can save some work if the factor is
980  // one
981  if (value == 1.)
982  return *this;
983 
984  const PetscInt petsc_i = index;
985  const PetscScalar new_value = static_cast<PetscScalar>(*this) * value;
986 
987  const PetscErrorCode ierr =
988  VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
989  AssertThrow(ierr == 0, ExcPETScError(ierr));
990 
991  return *this;
992  }
993 
994 
995 
996  inline const VectorReference &
997  VectorReference::operator/=(const PetscScalar &value) const
998  {
1001  ExcWrongMode(VectorOperation::insert, vector.last_action));
1002 
1004 
1006 
1007  // we have to do above actions in any
1008  // case to be consistent with the MPI
1009  // communication model (see the
1010  // comments in the documentation of
1011  // PETScWrappers::MPI::Vector), but we
1012  // can save some work if the factor is
1013  // one
1014  if (value == 1.)
1015  return *this;
1016 
1017  const PetscInt petsc_i = index;
1018  const PetscScalar new_value = static_cast<PetscScalar>(*this) / value;
1019 
1020  const PetscErrorCode ierr =
1021  VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1022  AssertThrow(ierr == 0, ExcPETScError(ierr));
1023 
1024  return *this;
1025  }
1026 
1027 
1028 
1029  inline PetscReal
1030  VectorReference::real() const
1031  {
1032 # ifndef PETSC_USE_COMPLEX
1033  return static_cast<PetscScalar>(*this);
1034 # else
1035  return PetscRealPart(static_cast<PetscScalar>(*this));
1036 # endif
1037  }
1038 
1039 
1040 
1041  inline PetscReal
1042  VectorReference::imag() const
1043  {
1044 # ifndef PETSC_USE_COMPLEX
1045  return PetscReal(0);
1046 # else
1047  return PetscImaginaryPart(static_cast<PetscScalar>(*this));
1048 # endif
1049  }
1050 
1051  } // namespace internal
1052 
1053  inline bool
1054  VectorBase::in_local_range(const size_type index) const
1055  {
1056  PetscInt begin, end;
1057  const PetscErrorCode ierr =
1058  VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
1059  AssertThrow(ierr == 0, ExcPETScError(ierr));
1060 
1061  return ((index >= static_cast<size_type>(begin)) &&
1062  (index < static_cast<size_type>(end)));
1063  }
1064 
1065 
1066  inline IndexSet
1068  {
1069  IndexSet is(size());
1070 
1071  // PETSc only allows for contiguous local ranges, so this is simple
1072  const std::pair<size_type, size_type> x = local_range();
1073  is.add_range(x.first, x.second);
1074  return is;
1075  }
1076 
1077 
1078 
1079  inline bool
1081  {
1082  return ghosted;
1083  }
1084 
1085 
1086 
1087  inline void
1089  {}
1090 
1091 
1092 
1093  inline internal::VectorReference
1094  VectorBase::operator()(const size_type index)
1095  {
1096  return internal::VectorReference(*this, index);
1097  }
1098 
1099 
1100 
1101  inline PetscScalar
1102  VectorBase::operator()(const size_type index) const
1103  {
1104  return static_cast<PetscScalar>(internal::VectorReference(*this, index));
1105  }
1106 
1107 
1108 
1109  inline internal::VectorReference VectorBase::operator[](const size_type index)
1110  {
1111  return operator()(index);
1112  }
1113 
1114 
1115 
1116  inline PetscScalar VectorBase::operator[](const size_type index) const
1117  {
1118  return operator()(index);
1119  }
1120 
1121  inline const MPI_Comm &
1123  {
1124  static MPI_Comm comm;
1125  PetscObjectGetComm(reinterpret_cast<PetscObject>(vector), &comm);
1126  return comm;
1127  }
1128 
1129  inline void
1130  VectorBase::extract_subvector_to(const std::vector<size_type> &indices,
1131  std::vector<PetscScalar> & values) const
1132  {
1133  Assert(indices.size() <= values.size(),
1134  ExcDimensionMismatch(indices.size(), values.size()));
1135  extract_subvector_to(indices.begin(), indices.end(), values.begin());
1136  }
1137 
1138  template <typename ForwardIterator, typename OutputIterator>
1139  inline void
1140  VectorBase::extract_subvector_to(const ForwardIterator indices_begin,
1141  const ForwardIterator indices_end,
1142  OutputIterator values_begin) const
1143  {
1144  const PetscInt n_idx = static_cast<PetscInt>(indices_end - indices_begin);
1145  if (n_idx == 0)
1146  return;
1147 
1148  // if we are dealing
1149  // with a parallel vector
1150  if (ghosted)
1151  {
1152  // there is the possibility
1153  // that the vector has
1154  // ghost elements. in that
1155  // case, we first need to
1156  // figure out which
1157  // elements we own locally,
1158  // then get a pointer to
1159  // the elements that are
1160  // stored here (both the
1161  // ones we own as well as
1162  // the ghost elements). in
1163  // this array, the locally
1164  // owned elements come
1165  // first followed by the
1166  // ghost elements whose
1167  // position we can get from
1168  // an index set
1169  PetscInt begin, end;
1170  PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1171  AssertThrow(ierr == 0, ExcPETScError(ierr));
1172 
1173  Vec locally_stored_elements = nullptr;
1174  ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
1175  AssertThrow(ierr == 0, ExcPETScError(ierr));
1176 
1177  PetscInt lsize;
1178  ierr = VecGetSize(locally_stored_elements, &lsize);
1179  AssertThrow(ierr == 0, ExcPETScError(ierr));
1180 
1181  PetscScalar *ptr;
1182  ierr = VecGetArray(locally_stored_elements, &ptr);
1183  AssertThrow(ierr == 0, ExcPETScError(ierr));
1184 
1185  for (PetscInt i = 0; i < n_idx; ++i)
1186  {
1187  const unsigned int index = *(indices_begin + i);
1188  if (index >= static_cast<unsigned int>(begin) &&
1189  index < static_cast<unsigned int>(end))
1190  {
1191  // local entry
1192  *(values_begin + i) = *(ptr + index - begin);
1193  }
1194  else
1195  {
1196  // ghost entry
1197  const unsigned int ghostidx =
1198  ghost_indices.index_within_set(index);
1199 
1200  AssertIndexRange(ghostidx + end - begin, lsize);
1201  *(values_begin + i) = *(ptr + ghostidx + end - begin);
1202  }
1203  }
1204 
1205  ierr = VecRestoreArray(locally_stored_elements, &ptr);
1206  AssertThrow(ierr == 0, ExcPETScError(ierr));
1207 
1208  ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
1209  AssertThrow(ierr == 0, ExcPETScError(ierr));
1210  }
1211  // if the vector is local or the
1212  // caller, then simply access the
1213  // element we are interested in
1214  else
1215  {
1216  PetscInt begin, end;
1217  PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1218  AssertThrow(ierr == 0, ExcPETScError(ierr));
1219 
1220  PetscScalar *ptr;
1221  ierr = VecGetArray(vector, &ptr);
1222  AssertThrow(ierr == 0, ExcPETScError(ierr));
1223 
1224  for (PetscInt i = 0; i < n_idx; ++i)
1225  {
1226  const unsigned int index = *(indices_begin + i);
1227 
1228  Assert(index >= static_cast<unsigned int>(begin) &&
1229  index < static_cast<unsigned int>(end),
1230  ExcInternalError());
1231 
1232  *(values_begin + i) = *(ptr + index - begin);
1233  }
1234 
1235  ierr = VecRestoreArray(vector, &ptr);
1236  AssertThrow(ierr == 0, ExcPETScError(ierr));
1237  }
1238  }
1239 
1240 # endif // DOXYGEN
1241 } // namespace PETScWrappers
1242 
1243 DEAL_II_NAMESPACE_CLOSE
1244 
1245 # endif // DEAL_II_WITH_PETSC
1246 
1247 #endif
1248 /*---------------------------- 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:1637
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:1519
SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
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:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
__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.cc:98
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()