Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
petsc_vector_base.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2004 - 2023 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
26
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
40template <typename number>
41class Vector;
42
43namespace PETScWrappers
44{
45 class VectorBase;
46}
47# endif
48
55namespace 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;
263
268 VectorBase();
269
274 VectorBase(const VectorBase &v);
275
280 explicit VectorBase(const Vec &v);
281
285 virtual ~VectorBase() override;
286
291 virtual void
292 clear();
293
304 void
305 compress(const VectorOperation::values operation);
306
310 VectorBase &
311 operator=(const VectorBase &);
312
326 VectorBase &
327 operator=(const PetscScalar s);
328
335 void
336 reinit(Vec v);
337
343 bool
344 operator==(const VectorBase &v) const;
345
351 bool
352 operator!=(const VectorBase &v) const;
353
358 size() const;
359
372 local_size() const;
373
383 locally_owned_size() const;
384
393 std::pair<size_type, size_type>
394 local_range() const;
395
400 bool
401 in_local_range(const size_type index) const;
402
418
425 bool
427
431 const IndexSet &
433
437 void
439
444 operator()(const size_type index);
445
449 PetscScalar
450 operator()(const size_type index) const;
451
458 operator[](const size_type index);
459
465 PetscScalar
466 operator[](const size_type index) const;
467
474 void
475 set(const std::vector<size_type> & indices,
476 const std::vector<PetscScalar> &values);
477
493 void
494 extract_subvector_to(const std::vector<size_type> &indices,
495 std::vector<PetscScalar> & values) const;
496
524 template <typename ForwardIterator, typename OutputIterator>
525 void
526 extract_subvector_to(const ForwardIterator indices_begin,
527 const ForwardIterator indices_end,
528 OutputIterator values_begin) const;
529
534 void
535 add(const std::vector<size_type> & indices,
536 const std::vector<PetscScalar> &values);
537
542 void
543 add(const std::vector<size_type> & indices,
544 const ::Vector<PetscScalar> &values);
545
551 void
552 add(const size_type n_elements,
553 const size_type * indices,
554 const PetscScalar *values);
555
562 PetscScalar
563 operator*(const VectorBase &vec) const;
564
569 norm_sqr() const;
570
574 PetscScalar
575 mean_value() const;
576
585 l1_norm() const;
586
592 l2_norm() const;
593
599 lp_norm(const real_type p) const;
600
606 linfty_norm() const;
607
627 PetscScalar
628 add_and_dot(const PetscScalar a, const VectorBase &V, const VectorBase &W);
629
639 min() const;
640
650 max() const;
651
657 bool
658 all_zero() const;
659
669 bool
670 is_non_negative() const;
671
675 VectorBase &
676 operator*=(const PetscScalar factor);
677
681 VectorBase &
682 operator/=(const PetscScalar factor);
683
687 VectorBase &
688 operator+=(const VectorBase &V);
689
693 VectorBase &
694 operator-=(const VectorBase &V);
695
700 void
701 add(const PetscScalar s);
702
706 void
707 add(const PetscScalar a, const VectorBase &V);
708
712 void
713 add(const PetscScalar a,
714 const VectorBase &V,
715 const PetscScalar b,
716 const VectorBase &W);
717
721 void
722 sadd(const PetscScalar s, const VectorBase &V);
723
727 void
728 sadd(const PetscScalar s, const PetscScalar a, const VectorBase &V);
729
735 void
736 scale(const VectorBase &scaling_factors);
737
741 void
742 equ(const PetscScalar a, const VectorBase &V);
743
751 void
752 write_ascii(const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
753
761 void
762 print(std::ostream & out,
763 const unsigned int precision = 3,
764 const bool scientific = true,
765 const bool across = true) const;
766
779 void
780 swap(VectorBase &v);
781
789 operator const Vec &() const;
790
796 Vec &
797 petsc_vector();
798
802 std::size_t
803 memory_consumption() const;
804
810
811 protected:
817
824
831
838
839 // Make the reference class a friend.
841
847 void
848 do_set_add_operation(const size_type n_elements,
849 const size_type * indices,
850 const PetscScalar *values,
851 const bool add_values);
852
856 void
858 };
859
860
861
862 // ------------------- inline and template functions --------------
863
871 inline void
873 {
874 u.swap(v);
875 }
876
877# ifndef DOXYGEN
878 namespace internal
879 {
880 inline VectorReference::VectorReference(const VectorBase &vector,
881 const size_type index)
882 : vector(vector)
883 , index(index)
884 {}
885
886
887 inline const VectorReference &
888 VectorReference::operator=(const VectorReference &r) const
889 {
890 // as explained in the class
891 // documentation, this is not the copy
892 // operator. so simply pass on to the
893 // "correct" assignment operator
894 *this = static_cast<PetscScalar>(r);
895
896 return *this;
897 }
898
899
900
901 inline VectorReference &
902 VectorReference::operator=(const VectorReference &r)
903 {
904 // as explained in the class
905 // documentation, this is not the copy
906 // operator. so simply pass on to the
907 // "correct" assignment operator
908 *this = static_cast<PetscScalar>(r);
909
910 return *this;
911 }
912
913
914
915 inline const VectorReference &
916 VectorReference::operator=(const PetscScalar &value) const
917 {
918 Assert((vector.last_action == VectorOperation::insert) ||
919 (vector.last_action == VectorOperation::unknown),
920 ExcWrongMode(VectorOperation::insert, vector.last_action));
921
922 Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
923
924 const PetscInt petsc_i = index;
925
926 const PetscErrorCode ierr =
927 VecSetValues(vector, 1, &petsc_i, &value, INSERT_VALUES);
928 AssertThrow(ierr == 0, ExcPETScError(ierr));
929
930 vector.last_action = VectorOperation::insert;
931
932 return *this;
933 }
934
935
936
937 inline const VectorReference &
938 VectorReference::operator+=(const PetscScalar &value) const
939 {
940 Assert((vector.last_action == VectorOperation::add) ||
941 (vector.last_action == VectorOperation::unknown),
942 ExcWrongMode(VectorOperation::add, vector.last_action));
943
944 Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
945
946 vector.last_action = VectorOperation::add;
947
948 // we have to do above actions in any
949 // case to be consistent with the MPI
950 // communication model (see the
951 // comments in the documentation of
952 // PETScWrappers::MPI::Vector), but we
953 // can save some work if the addend is
954 // zero
955 if (value == PetscScalar())
956 return *this;
957
958 // use the PETSc function to add something
959 const PetscInt petsc_i = index;
960 const PetscErrorCode ierr =
961 VecSetValues(vector, 1, &petsc_i, &value, ADD_VALUES);
962 AssertThrow(ierr == 0, ExcPETScError(ierr));
963
964
965 return *this;
966 }
967
968
969
970 inline const VectorReference &
971 VectorReference::operator-=(const PetscScalar &value) const
972 {
973 Assert((vector.last_action == VectorOperation::add) ||
974 (vector.last_action == VectorOperation::unknown),
975 ExcWrongMode(VectorOperation::add, vector.last_action));
976
977 Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
978
979 vector.last_action = VectorOperation::add;
980
981 // we have to do above actions in any
982 // case to be consistent with the MPI
983 // communication model (see the
984 // comments in the documentation of
985 // PETScWrappers::MPI::Vector), but we
986 // can save some work if the addend is
987 // zero
988 if (value == PetscScalar())
989 return *this;
990
991 // use the PETSc function to
992 // add something
993 const PetscInt petsc_i = index;
994 const PetscScalar subtractand = -value;
995 const PetscErrorCode ierr =
996 VecSetValues(vector, 1, &petsc_i, &subtractand, ADD_VALUES);
997 AssertThrow(ierr == 0, ExcPETScError(ierr));
998
999 return *this;
1000 }
1001
1002
1003
1004 inline const VectorReference &
1005 VectorReference::operator*=(const PetscScalar &value) const
1006 {
1007 Assert((vector.last_action == VectorOperation::insert) ||
1008 (vector.last_action == VectorOperation::unknown),
1009 ExcWrongMode(VectorOperation::insert, vector.last_action));
1010
1011 Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
1012
1013 vector.last_action = VectorOperation::insert;
1014
1015 // we have to do above actions in any
1016 // case to be consistent with the MPI
1017 // communication model (see the
1018 // comments in the documentation of
1019 // PETScWrappers::MPI::Vector), but we
1020 // can save some work if the factor is
1021 // one
1022 if (value == 1.)
1023 return *this;
1024
1025 const PetscInt petsc_i = index;
1026 const PetscScalar new_value = static_cast<PetscScalar>(*this) * value;
1027
1028 const PetscErrorCode ierr =
1029 VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1030 AssertThrow(ierr == 0, ExcPETScError(ierr));
1031
1032 return *this;
1033 }
1034
1035
1036
1037 inline const VectorReference &
1038 VectorReference::operator/=(const PetscScalar &value) const
1039 {
1040 Assert((vector.last_action == VectorOperation::insert) ||
1041 (vector.last_action == VectorOperation::unknown),
1042 ExcWrongMode(VectorOperation::insert, vector.last_action));
1043
1044 Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
1045
1046 vector.last_action = VectorOperation::insert;
1047
1048 // we have to do above actions in any
1049 // case to be consistent with the MPI
1050 // communication model (see the
1051 // comments in the documentation of
1052 // PETScWrappers::MPI::Vector), but we
1053 // can save some work if the factor is
1054 // one
1055 if (value == 1.)
1056 return *this;
1057
1058 const PetscInt petsc_i = index;
1059 const PetscScalar new_value = static_cast<PetscScalar>(*this) / value;
1060
1061 const PetscErrorCode ierr =
1062 VecSetValues(vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1063 AssertThrow(ierr == 0, ExcPETScError(ierr));
1064
1065 return *this;
1066 }
1067
1068
1069
1070 inline PetscReal
1071 VectorReference::real() const
1072 {
1073# ifndef PETSC_USE_COMPLEX
1074 return static_cast<PetscScalar>(*this);
1075# else
1076 return PetscRealPart(static_cast<PetscScalar>(*this));
1077# endif
1078 }
1079
1080
1081
1082 inline PetscReal
1083 VectorReference::imag() const
1084 {
1085# ifndef PETSC_USE_COMPLEX
1086 return PetscReal(0);
1087# else
1088 return PetscImaginaryPart(static_cast<PetscScalar>(*this));
1089# endif
1090 }
1091
1092 } // namespace internal
1093
1094 inline bool
1095 VectorBase::in_local_range(const size_type index) const
1096 {
1097 PetscInt begin, end;
1098 const PetscErrorCode ierr =
1099 VecGetOwnershipRange(static_cast<const Vec &>(vector), &begin, &end);
1100 AssertThrow(ierr == 0, ExcPETScError(ierr));
1101
1102 return ((index >= static_cast<size_type>(begin)) &&
1103 (index < static_cast<size_type>(end)));
1104 }
1105
1106
1107 inline IndexSet
1108 VectorBase::locally_owned_elements() const
1109 {
1110 IndexSet is(size());
1111
1112 // PETSc only allows for contiguous local ranges, so this is simple
1113 const std::pair<size_type, size_type> x = local_range();
1114 is.add_range(x.first, x.second);
1115 return is;
1116 }
1117
1118
1119
1120 inline bool
1121 VectorBase::has_ghost_elements() const
1122 {
1123 return ghosted;
1124 }
1125
1126
1127 inline const IndexSet &
1128 VectorBase::ghost_elements() const
1129 {
1130 return ghost_indices;
1131 }
1132
1133
1134 inline void
1135 VectorBase::update_ghost_values() const
1136 {
1137 if (ghosted)
1138 {
1139 PetscErrorCode ierr;
1140
1141 ierr = VecGhostUpdateBegin(vector, INSERT_VALUES, SCATTER_FORWARD);
1142 AssertThrow(ierr == 0, ExcPETScError(ierr));
1143 ierr = VecGhostUpdateEnd(vector, INSERT_VALUES, SCATTER_FORWARD);
1144 AssertThrow(ierr == 0, ExcPETScError(ierr));
1145 }
1146 }
1147
1148
1149
1150 inline internal::VectorReference
1151 VectorBase::operator()(const size_type index)
1152 {
1153 return internal::VectorReference(*this, index);
1154 }
1155
1156
1157
1158 inline PetscScalar
1159 VectorBase::operator()(const size_type index) const
1160 {
1161 return static_cast<PetscScalar>(internal::VectorReference(*this, index));
1162 }
1163
1164
1165
1166 inline internal::VectorReference
1167 VectorBase::operator[](const size_type index)
1168 {
1169 return operator()(index);
1170 }
1171
1172
1173
1174 inline PetscScalar
1175 VectorBase::operator[](const size_type index) const
1176 {
1177 return operator()(index);
1178 }
1179
1180 inline MPI_Comm
1181 VectorBase::get_mpi_communicator() const
1182 {
1183 return PetscObjectComm(reinterpret_cast<PetscObject>(vector));
1184 }
1185
1186 inline void
1187 VectorBase::extract_subvector_to(const std::vector<size_type> &indices,
1188 std::vector<PetscScalar> & values) const
1189 {
1190 Assert(indices.size() <= values.size(),
1191 ExcDimensionMismatch(indices.size(), values.size()));
1192 extract_subvector_to(indices.begin(), indices.end(), values.begin());
1193 }
1194
1195 template <typename ForwardIterator, typename OutputIterator>
1196 inline void
1197 VectorBase::extract_subvector_to(const ForwardIterator indices_begin,
1198 const ForwardIterator indices_end,
1199 OutputIterator values_begin) const
1200 {
1201 const PetscInt n_idx = static_cast<PetscInt>(indices_end - indices_begin);
1202 if (n_idx == 0)
1203 return;
1204
1205 // if we are dealing
1206 // with a parallel vector
1207 if (ghosted)
1208 {
1209 // there is the possibility
1210 // that the vector has
1211 // ghost elements. in that
1212 // case, we first need to
1213 // figure out which
1214 // elements we own locally,
1215 // then get a pointer to
1216 // the elements that are
1217 // stored here (both the
1218 // ones we own as well as
1219 // the ghost elements). in
1220 // this array, the locally
1221 // owned elements come
1222 // first followed by the
1223 // ghost elements whose
1224 // position we can get from
1225 // an index set
1226 PetscInt begin, end;
1227 PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1228 AssertThrow(ierr == 0, ExcPETScError(ierr));
1229
1230 Vec locally_stored_elements = nullptr;
1231 ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
1232 AssertThrow(ierr == 0, ExcPETScError(ierr));
1233
1234 PetscInt lsize;
1235 ierr = VecGetSize(locally_stored_elements, &lsize);
1236 AssertThrow(ierr == 0, ExcPETScError(ierr));
1237
1238 const PetscScalar *ptr;
1239 ierr = VecGetArrayRead(locally_stored_elements, &ptr);
1240 AssertThrow(ierr == 0, ExcPETScError(ierr));
1241
1242 for (PetscInt i = 0; i < n_idx; ++i)
1243 {
1244 const unsigned int index = *(indices_begin + i);
1245 if (index >= static_cast<unsigned int>(begin) &&
1246 index < static_cast<unsigned int>(end))
1247 {
1248 // local entry
1249 *(values_begin + i) = *(ptr + index - begin);
1250 }
1251 else
1252 {
1253 // ghost entry
1254 const unsigned int ghostidx =
1255 ghost_indices.index_within_set(index);
1256
1257 AssertIndexRange(ghostidx + end - begin, lsize);
1258 *(values_begin + i) = *(ptr + ghostidx + end - begin);
1259 }
1260 }
1261
1262 ierr = VecRestoreArrayRead(locally_stored_elements, &ptr);
1263 AssertThrow(ierr == 0, ExcPETScError(ierr));
1264
1265 ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
1266 AssertThrow(ierr == 0, ExcPETScError(ierr));
1267 }
1268 // if the vector is local or the
1269 // caller, then simply access the
1270 // element we are interested in
1271 else
1272 {
1273 PetscInt begin, end;
1274 PetscErrorCode ierr = VecGetOwnershipRange(vector, &begin, &end);
1275 AssertThrow(ierr == 0, ExcPETScError(ierr));
1276
1277 const PetscScalar *ptr;
1278 ierr = VecGetArrayRead(vector, &ptr);
1279 AssertThrow(ierr == 0, ExcPETScError(ierr));
1280
1281 for (PetscInt i = 0; i < n_idx; ++i)
1282 {
1283 const unsigned int index = *(indices_begin + i);
1284
1285 Assert(index >= static_cast<unsigned int>(begin) &&
1286 index < static_cast<unsigned int>(end),
1287 ExcMessage("You are accessing elements of a vector without "
1288 "ghost elements that are not actually owned by "
1289 "this vector. A typical case where this may "
1290 "happen is if you are passing a non-ghosted "
1291 "(completely distributed) vector to a function "
1292 "that expects a vector that stores ghost "
1293 "elements for all locally relevant or locally "
1294 "active vector entries."));
1295
1296 *(values_begin + i) = *(ptr + index - begin);
1297 }
1298
1299 ierr = VecRestoreArrayRead(vector, &ptr);
1300 AssertThrow(ierr == 0, ExcPETScError(ierr));
1301 }
1302 }
1303
1304# endif // DOXYGEN
1305} // namespace PETScWrappers
1306
1308
1309#endif // DEAL_II_WITH_PETSC
1310
1311#endif
real_type lp_norm(const real_type p) const
VectorBase & operator+=(const VectorBase &V)
PetscScalar mean_value() const
VectorOperation::values last_action
PetscScalar operator*(const VectorBase &vec) const
bool operator==(const VectorBase &v) const
VectorBase & operator*=(const PetscScalar factor)
VectorBase & operator-=(const VectorBase &V)
bool operator!=(const VectorBase &v) const
IndexSet locally_owned_elements() const
bool in_local_range(const size_type index) const
std::size_t memory_consumption() const
internal::VectorReference reference
VectorBase & operator/=(const PetscScalar factor)
friend class internal::VectorReference
PetscScalar operator[](const size_type index) const
MPI_Comm get_mpi_communicator() const
void scale(const VectorBase &scaling_factors)
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< PetscScalar > &values) const
void update_ghost_values() const
std::pair< size_type, size_type > local_range() const
void sadd(const PetscScalar s, const VectorBase &V)
void compress(const VectorOperation::values operation)
void set(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
PetscScalar add_and_dot(const PetscScalar a, const VectorBase &V, const VectorBase &W)
const IndexSet & ghost_elements() const
reference operator()(const size_type index)
const internal::VectorReference const_reference
void add(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)
void write_ascii(const PetscViewerFormat format=PETSC_VIEWER_DEFAULT)
bool has_ghost_elements() const
reference operator[](const size_type index)
size_type locally_owned_size() const
void equ(const PetscScalar a, const VectorBase &V)
void do_set_add_operation(const size_type n_elements, const size_type *indices, const PetscScalar *values, const bool add_values)
VectorBase & operator=(const VectorBase &)
void swap(VectorBase &u, VectorBase &v)
PetscScalar operator()(const size_type index) const
void extract_subvector_to(const ForwardIterator indices_begin, const ForwardIterator indices_end, OutputIterator values_begin) const
#define DEAL_II_DEPRECATED
Definition config.h:172
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcGhostsPresent()
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:533
#define AssertIndexRange(index, range)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition exceptions.h:556
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
types::global_dof_index size_type
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
unsigned int global_dof_index
Definition types.h:82