Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
petsc_vector_base.h
1 // ---------------------------------------------------------------------
2 // @f$Id: petsc_vector_base.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2004 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__petsc_vector_base_h
18 #define __deal2__petsc_vector_base_h
19 
20 
21 #include <deal.II/base/config.h>
22 
23 #ifdef DEAL_II_WITH_PETSC
24 
25 # include <deal.II/base/subscriptor.h>
26 # include <deal.II/lac/exceptions.h>
27 # include <deal.II/lac/vector.h>
28 
29 # include <vector>
30 # include <utility>
31 
32 # include <petscvec.h>
33 # include <deal.II/base/index_set.h>
34 
36 
37 // forward declaration
38 template <typename number> class Vector;
39 
40 
48 namespace PETScWrappers
49 {
50  // forward declaration
51  class VectorBase;
52 
62  namespace internal
63  {
83  class VectorReference
84  {
85  public:
90 
91  private:
97  VectorReference (const VectorBase &vector,
98  const size_type index);
99 
100 
101  public:
102 
123  const VectorReference &operator = (const VectorReference &r) const;
124 
133  VectorReference &operator = (const VectorReference &r);
134 
139  const VectorReference &operator = (const PetscScalar &s) const;
140 
145  const VectorReference &operator += (const PetscScalar &s) const;
146 
151  const VectorReference &operator -= (const PetscScalar &s) const;
152 
157  const VectorReference &operator *= (const PetscScalar &s) const;
158 
163  const VectorReference &operator /= (const PetscScalar &s) const;
164 
171  operator PetscScalar () const;
172 
176  DeclException1 (ExcPETScError,
177  int,
178  << "An error with error number " << arg1
179  << " occurred while calling a PETSc function");
183  DeclException3 (ExcAccessToNonlocalElement,
184  int, int, int,
185  << "You tried to access element " << arg1
186  << " of a distributed vector, but only elements "
187  << arg2 << " through " << arg3
188  << " are stored locally and can be accessed.");
192  DeclException2 (ExcWrongMode,
193  int, int,
194  << "You tried to do a "
195  << (arg1 == 1 ?
196  "'set'" :
197  (arg1 == 2 ?
198  "'add'" : "???"))
199  << " operation but the vector is currently in "
200  << (arg2 == 1 ?
201  "'set'" :
202  (arg2 == 2 ?
203  "'add'" : "???"))
204  << " mode. You first have to call 'compress()'.");
205 
206  private:
211  const VectorBase &vector;
212 
217  const size_type index;
218 
224  friend class ::PETScWrappers::VectorBase;
225  };
226  }
256  class VectorBase : public Subscriptor
257  {
258  public:
266  typedef PetscScalar value_type;
267  typedef PetscReal real_type;
268  typedef types::global_dof_index size_type;
269  typedef internal::VectorReference reference;
270  typedef const internal::VectorReference const_reference;
271 
277  VectorBase ();
278 
284  VectorBase (const VectorBase &v);
285 
293  explicit VectorBase (const Vec &v);
294 
298  virtual ~VectorBase ();
299 
312  void compress (::VectorOperation::values operation);
313 
318 
339  VectorBase &operator = (const PetscScalar s);
340 
349  bool operator == (const VectorBase &v) const;
350 
359  bool operator != (const VectorBase &v) const;
360 
365  size_type size () const;
366 
380  size_type local_size () const;
381 
398  std::pair<size_type, size_type>
399  local_range () const;
400 
406  bool in_local_range (const size_type index) const;
407 
423 
428  bool has_ghost_elements() const;
429 
434  reference
435  operator () (const size_type index);
436 
441  PetscScalar
442  operator () (const size_type index) const;
443 
450  reference
451  operator [] (const size_type index);
452 
460  PetscScalar
461  operator [] (const size_type index) const;
462 
473  void set (const std::vector<size_type> &indices,
474  const std::vector<PetscScalar> &values);
475 
486  void extract_subvector_to (const std::vector<size_type> &indices,
487  std::vector<PetscScalar> &values) const;
488 
493  template <typename ForwardIterator, typename OutputIterator>
494  void extract_subvector_to (const ForwardIterator indices_begin,
495  const ForwardIterator indices_end,
496  OutputIterator values_begin) const;
497 
504  void add (const std::vector<size_type> &indices,
505  const std::vector<PetscScalar> &values);
506 
514  void add (const std::vector<size_type> &indices,
515  const ::Vector<PetscScalar> &values);
516 
526  void add (const size_type n_elements,
527  const size_type *indices,
528  const PetscScalar *values);
529 
535  PetscScalar operator * (const VectorBase &vec) const;
536 
540  real_type norm_sqr () const;
541 
545  PetscScalar mean_value () const;
546 
551  real_type l1_norm () const;
552 
558  real_type l2_norm () const;
559 
566  real_type lp_norm (const real_type p) const;
567 
572  real_type linfty_norm () const;
573 
578  real_type normalize () const;
579 
583  real_type min () const;
584 
588  real_type max () const;
589 
593  VectorBase &abs ();
594 
598  VectorBase &conjugate ();
599 
608  VectorBase &mult ();
609 
617  VectorBase &mult (const VectorBase &v);
618 
625  VectorBase &mult (const VectorBase &u,
626  const VectorBase &v);
627 
636  bool all_zero () const;
637 
646  bool is_non_negative () const;
647 
652  VectorBase &operator *= (const PetscScalar factor);
653 
658  VectorBase &operator /= (const PetscScalar factor);
659 
664  VectorBase &operator += (const VectorBase &V);
665 
670  VectorBase &operator -= (const VectorBase &V);
671 
677  void add (const PetscScalar s);
678 
683  void add (const VectorBase &V);
684 
689  void add (const PetscScalar a, const VectorBase &V);
690 
695  void add (const PetscScalar a, const VectorBase &V,
696  const PetscScalar b, const VectorBase &W);
697 
703  void sadd (const PetscScalar s,
704  const VectorBase &V);
705 
710  void sadd (const PetscScalar s,
711  const PetscScalar a,
712  const VectorBase &V);
713 
717  void sadd (const PetscScalar s,
718  const PetscScalar a,
719  const VectorBase &V,
720  const PetscScalar b,
721  const VectorBase &W);
722 
727  void sadd (const PetscScalar s,
728  const PetscScalar a,
729  const VectorBase &V,
730  const PetscScalar b,
731  const VectorBase &W,
732  const PetscScalar c,
733  const VectorBase &X);
734 
744  void scale (const VectorBase &scaling_factors);
745 
749  void equ (const PetscScalar a, const VectorBase &V);
750 
754  void equ (const PetscScalar a, const VectorBase &V,
755  const PetscScalar b, const VectorBase &W);
756 
773  void ratio (const VectorBase &a,
774  const VectorBase &b);
775 
784 
785 
795  void write_ascii (const PetscViewerFormat format = PETSC_VIEWER_DEFAULT) ;
796 
811  void print (std::ostream &out,
812  const unsigned int precision = 3,
813  const bool scientific = true,
814  const bool across = true) const;
815 
839  void swap (VectorBase &v);
840 
852  operator const Vec &() const;
853 
859  std::size_t memory_consumption () const;
860 
866  virtual const MPI_Comm &get_mpi_communicator () const;
867 
868  protected:
874  Vec vector;
875 
883  bool ghosted;
884 
893 
902  mutable ::VectorOperation::values last_action;
903 
907  friend class internal::VectorReference;
908 
917 
926  void do_set_add_operation (const size_type n_elements,
927  const size_type *indices,
928  const PetscScalar *values,
929  const bool add_values);
930 
931 
932  };
933 
934 
935 
936 // ------------------- inline and template functions --------------
937 
946  inline
947  void swap (VectorBase &u, VectorBase &v)
948  {
949  u.swap (v);
950  }
951 
952 #ifndef DOXYGEN
953  namespace internal
954  {
955  inline
956  VectorReference::VectorReference (const VectorBase &vector,
957  const size_type index)
958  :
959  vector (vector),
960  index (index)
961  {}
962 
963 
964  inline
965  const VectorReference &
966  VectorReference::operator = (const VectorReference &r) const
967  {
968  // as explained in the class
969  // documentation, this is not the copy
970  // operator. so simply pass on to the
971  // "correct" assignment operator
972  *this = static_cast<PetscScalar> (r);
973 
974  return *this;
975  }
976 
977 
978 
979  inline
980  VectorReference &
981  VectorReference::operator = (const VectorReference &r)
982  {
983  // as explained in the class
984  // documentation, this is not the copy
985  // operator. so simply pass on to the
986  // "correct" assignment operator
987  *this = static_cast<PetscScalar> (r);
988 
989  return *this;
990  }
991 
992 
993 
994  inline
995  const VectorReference &
996  VectorReference::operator = (const PetscScalar &value) const
997  {
998  Assert ((vector.last_action == VectorOperation::insert)
999  ||
1000  (vector.last_action == VectorOperation::unknown),
1001  ExcWrongMode (VectorOperation::insert,
1002  vector.last_action));
1003 
1004  Assert (!vector.has_ghost_elements(), ExcGhostsPresent());
1005 
1006  const PetscInt petsc_i = index;
1007 
1008  const int ierr
1009  = VecSetValues (vector, 1, &petsc_i, &value, INSERT_VALUES);
1010  AssertThrow (ierr == 0, ExcPETScError(ierr));
1011 
1012  vector.last_action = VectorOperation::insert;
1013 
1014  return *this;
1015  }
1016 
1017 
1018 
1019  inline
1020  const VectorReference &
1021  VectorReference::operator += (const PetscScalar &value) const
1022  {
1023  Assert ((vector.last_action == VectorOperation::add)
1024  ||
1025  (vector.last_action == VectorOperation::unknown),
1026  ExcWrongMode (VectorOperation::add,
1027  vector.last_action));
1028 
1029  Assert (!vector.has_ghost_elements(), ExcGhostsPresent());
1030 
1031  vector.last_action = VectorOperation::add;
1032 
1033  // we have to do above actions in any
1034  // case to be consistent with the MPI
1035  // communication model (see the
1036  // comments in the documentation of
1037  // PETScWrappers::MPI::Vector), but we
1038  // can save some work if the addend is
1039  // zero
1040  if (value == PetscScalar())
1041  return *this;
1042 
1043  // use the PETSc function to add something
1044  const PetscInt petsc_i = index;
1045  const int ierr
1046  = VecSetValues (vector, 1, &petsc_i, &value, ADD_VALUES);
1047  AssertThrow (ierr == 0, ExcPETScError(ierr));
1048 
1049 
1050  return *this;
1051  }
1052 
1053 
1054 
1055  inline
1056  const VectorReference &
1057  VectorReference::operator -= (const PetscScalar &value) const
1058  {
1059  Assert ((vector.last_action == VectorOperation::add)
1060  ||
1061  (vector.last_action == VectorOperation::unknown),
1062  ExcWrongMode (VectorOperation::add,
1063  vector.last_action));
1064 
1065  Assert (!vector.has_ghost_elements(), ExcGhostsPresent());
1066 
1067  vector.last_action = VectorOperation::add;
1068 
1069  // we have to do above actions in any
1070  // case to be consistent with the MPI
1071  // communication model (see the
1072  // comments in the documentation of
1073  // PETScWrappers::MPI::Vector), but we
1074  // can save some work if the addend is
1075  // zero
1076  if (value == PetscScalar())
1077  return *this;
1078 
1079  // use the PETSc function to
1080  // add something
1081  const PetscInt petsc_i = index;
1082  const PetscScalar subtractand = -value;
1083  const int ierr
1084  = VecSetValues (vector, 1, &petsc_i, &subtractand, ADD_VALUES);
1085  AssertThrow (ierr == 0, ExcPETScError(ierr));
1086 
1087  return *this;
1088  }
1089 
1090 
1091 
1092  inline
1093  const VectorReference &
1094  VectorReference::operator *= (const PetscScalar &value) const
1095  {
1096  Assert ((vector.last_action == VectorOperation::insert)
1097  ||
1098  (vector.last_action == VectorOperation::unknown),
1099  ExcWrongMode (VectorOperation::insert,
1100  vector.last_action));
1101 
1102  Assert (!vector.has_ghost_elements(), ExcGhostsPresent());
1103 
1104  vector.last_action = VectorOperation::insert;
1105 
1106  // we have to do above actions in any
1107  // case to be consistent with the MPI
1108  // communication model (see the
1109  // comments in the documentation of
1110  // PETScWrappers::MPI::Vector), but we
1111  // can save some work if the factor is
1112  // one
1113  if (value == 1.)
1114  return *this;
1115 
1116  const PetscInt petsc_i = index;
1117  const PetscScalar new_value
1118  = static_cast<PetscScalar>(*this) * value;
1119 
1120  const int ierr
1121  = VecSetValues (vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1122  AssertThrow (ierr == 0, ExcPETScError(ierr));
1123 
1124  return *this;
1125  }
1126 
1127 
1128 
1129  inline
1130  const VectorReference &
1131  VectorReference::operator /= (const PetscScalar &value) const
1132  {
1133  Assert ((vector.last_action == VectorOperation::insert)
1134  ||
1135  (vector.last_action == VectorOperation::unknown),
1136  ExcWrongMode (VectorOperation::insert,
1137  vector.last_action));
1138 
1139  Assert (!vector.has_ghost_elements(), ExcGhostsPresent());
1140 
1141  vector.last_action = VectorOperation::insert;
1142 
1143  // we have to do above actions in any
1144  // case to be consistent with the MPI
1145  // communication model (see the
1146  // comments in the documentation of
1147  // PETScWrappers::MPI::Vector), but we
1148  // can save some work if the factor is
1149  // one
1150  if (value == 1.)
1151  return *this;
1152 
1153  const PetscInt petsc_i = index;
1154  const PetscScalar new_value
1155  = static_cast<PetscScalar>(*this) / value;
1156 
1157  const int ierr
1158  = VecSetValues (vector, 1, &petsc_i, &new_value, INSERT_VALUES);
1159  AssertThrow (ierr == 0, ExcPETScError(ierr));
1160 
1161  return *this;
1162  }
1163  }
1164 
1165 
1166 
1167  inline
1168  bool
1169  VectorBase::in_local_range (const size_type index) const
1170  {
1171  PetscInt begin, end;
1172  const int ierr = VecGetOwnershipRange (static_cast<const Vec &>(vector),
1173  &begin, &end);
1174  AssertThrow (ierr == 0, ExcPETScError(ierr));
1175 
1176  return ((index >= static_cast<size_type>(begin)) &&
1177  (index < static_cast<size_type>(end)));
1178  }
1179 
1180 
1181  inline
1182  IndexSet
1183  VectorBase::locally_owned_elements() const
1184  {
1185  IndexSet is (size());
1186 
1187  // PETSc only allows for contiguous local ranges, so this is simple
1188  const std::pair<size_type, size_type> x = local_range();
1189  is.add_range (x.first, x.second);
1190  return is;
1191  }
1192 
1193 
1194 
1195  inline
1196  bool
1197  VectorBase::has_ghost_elements() const
1198  {
1199  return ghosted;
1200  }
1201 
1202 
1203 
1204  inline
1205  internal::VectorReference
1206  VectorBase::operator () (const size_type index)
1207  {
1208  return internal::VectorReference (*this, index);
1209  }
1210 
1211 
1212 
1213  inline
1214  PetscScalar
1215  VectorBase::operator () (const size_type index) const
1216  {
1217  return static_cast<PetscScalar>(internal::VectorReference (*this, index));
1218  }
1219 
1220 
1221 
1222  inline
1223  internal::VectorReference
1224  VectorBase::operator [] (const size_type index)
1225  {
1226  return operator()(index);
1227  }
1228 
1229 
1230 
1231  inline
1232  PetscScalar
1233  VectorBase::operator [] (const size_type index) const
1234  {
1235  return operator()(index);
1236  }
1237 
1238  inline
1239  const MPI_Comm &
1240  VectorBase::get_mpi_communicator () const
1241  {
1242  static MPI_Comm comm;
1243  PetscObjectGetComm((PetscObject)vector, &comm);
1244  return comm;
1245  }
1246 
1247  inline
1248  void VectorBase::extract_subvector_to (const std::vector<size_type> &indices,
1249  std::vector<PetscScalar> &values) const
1250  {
1251  extract_subvector_to(&(indices[0]), &(indices[0]) + indices.size(), &(values[0]));
1252  }
1253 
1254  template <typename ForwardIterator, typename OutputIterator>
1255  inline
1256  void VectorBase::extract_subvector_to (const ForwardIterator indices_begin,
1257  const ForwardIterator indices_end,
1258  OutputIterator values_begin) const
1259  {
1260  const PetscInt n_idx = static_cast<PetscInt>(indices_end - indices_begin);
1261  if (n_idx == 0)
1262  return;
1263 
1264  // if we are dealing
1265  // with a parallel vector
1266  if (ghosted )
1267  {
1268 
1269  int ierr;
1270 
1271  // there is the possibility
1272  // that the vector has
1273  // ghost elements. in that
1274  // case, we first need to
1275  // figure out which
1276  // elements we own locally,
1277  // then get a pointer to
1278  // the elements that are
1279  // stored here (both the
1280  // ones we own as well as
1281  // the ghost elements). in
1282  // this array, the locally
1283  // owned elements come
1284  // first followed by the
1285  // ghost elements whose
1286  // position we can get from
1287  // an index set
1288  PetscInt begin, end, i;
1289  ierr = VecGetOwnershipRange (vector, &begin, &end);
1290  AssertThrow (ierr == 0, ExcPETScError(ierr));
1291 
1292  Vec locally_stored_elements = PETSC_NULL;
1293  ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
1294  AssertThrow (ierr == 0, ExcPETScError(ierr));
1295 
1296  PetscInt lsize;
1297  ierr = VecGetSize(locally_stored_elements, &lsize);
1298  AssertThrow (ierr == 0, ExcPETScError(ierr));
1299 
1300  PetscScalar *ptr;
1301  ierr = VecGetArray(locally_stored_elements, &ptr);
1302  AssertThrow (ierr == 0, ExcPETScError(ierr));
1303 
1304  for (i = 0; i < n_idx; i++)
1305  {
1306  const unsigned int index = *(indices_begin+i);
1307  if ( index>=static_cast<unsigned int>(begin)
1308  && index<static_cast<unsigned int>(end) )
1309  {
1310  //local entry
1311  *(values_begin+i) = *(ptr+index-begin);
1312  }
1313  else
1314  {
1315  //ghost entry
1316  const unsigned int ghostidx
1317  = ghost_indices.index_within_set(index);
1318 
1319  Assert(ghostidx+end-begin<(unsigned int)lsize, ExcInternalError());
1320  *(values_begin+i) = *(ptr+ghostidx+end-begin);
1321  }
1322  }
1323 
1324  ierr = VecRestoreArray(locally_stored_elements, &ptr);
1325  AssertThrow (ierr == 0, ExcPETScError(ierr));
1326 
1327  ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
1328  AssertThrow (ierr == 0, ExcPETScError(ierr));
1329 
1330  }
1331  // if the vector is local or the
1332  // caller, then simply access the
1333  // element we are interested in
1334  else
1335  {
1336  int ierr;
1337 
1338  PetscInt begin, end;
1339  ierr = VecGetOwnershipRange (vector, &begin, &end);
1340  AssertThrow (ierr == 0, ExcPETScError(ierr));
1341 
1342  PetscScalar *ptr;
1343  ierr = VecGetArray(vector, &ptr);
1344  AssertThrow (ierr == 0, ExcPETScError(ierr));
1345 
1346  for (PetscInt i = 0; i < n_idx; i++)
1347  {
1348  const unsigned int index = *(indices_begin+i);
1349 
1350  Assert(index>=static_cast<unsigned int>(begin)
1351  && index<static_cast<unsigned int>(end), ExcInternalError());
1352 
1353  *(values_begin+i) = *(ptr+index-begin);
1354  }
1355 
1356  ierr = VecRestoreArray(vector, &ptr);
1357  AssertThrow (ierr == 0, ExcPETScError(ierr));
1358 
1359  }
1360  }
1361 
1362 #endif // DOXYGEN
1363 }
1364 
1365 DEAL_II_NAMESPACE_CLOSE
1366 
1367 #endif // DEAL_II_WITH_PETSC
1368 
1369 /*---------------------------- petsc_vector_base.h ---------------------------*/
1370 
1371 #endif
1372 /*---------------------------- petsc_vector_base.h ---------------------------*/
real_type l2_norm() const
real_type linfty_norm() const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:525
real_type l1_norm() const
void sadd(const PetscScalar s, const VectorBase &V)
void ratio(const VectorBase &a, const VectorBase &b)
bool has_ghost_elements() const
void update_ghost_values() const DEAL_II_DEPRECATED
STL namespace.
real_type norm_sqr() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
void scale(const VectorBase &scaling_factors)
real_type min() const
real_type lp_norm(const real_type p) const
void compress() DEAL_II_DEPRECATED
virtual const MPI_Comm & get_mpi_communicator() const
::ExceptionBase & ExcGhostsPresent()
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
mutable::VectorOperation::values last_action
size_type size() const
unsigned int global_dof_index
Definition: types.h:100
std::pair< size_type, size_type > local_range() const
void do_set_add_operation(const size_type n_elements, const size_type *indices, const PetscScalar *values, const bool add_values)
#define Assert(cond, exc)
Definition: exceptions.h:299
PetscScalar mean_value() const
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< PetscScalar > &values) const
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
void equ(const PetscScalar a, const VectorBase &V)
void swap(VectorBase &v)
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
friend class internal::VectorReference
size_type local_size() const
bool is_non_negative() const
std::size_t memory_consumption() const
void add(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)
void write_ascii(const PetscViewerFormat format=PETSC_VIEWER_DEFAULT)
real_type max() const
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:535
::ExceptionBase & ExcInternalError()
IndexSet locally_owned_elements() const
VectorBase & conjugate()
real_type normalize() const
bool in_local_range(const size_type index) const
void set(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)