Reference documentation for deal.II version Git 035d1ad 2018-08-17 19:42:27 -0400
vector.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2018 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_vector_h
17 #define dealii_vector_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/index_set.h>
24 #include <deal.II/base/logstream.h>
25 #include <deal.II/base/subscriptor.h>
26 
27 #include <deal.II/differentiation/ad/ad_number_traits.h>
28 
29 #include <deal.II/lac/vector_operation.h>
30 #include <deal.II/lac/vector_type_traits.h>
31 
32 // boost::serialization::make_array used to be in array.hpp, but was
33 // moved to a different file in BOOST 1.64
34 #include <boost/version.hpp>
35 #if BOOST_VERSION >= 106400
36 # include <boost/serialization/array_wrapper.hpp>
37 #else
38 # include <boost/serialization/array.hpp>
39 #endif
40 #include <boost/serialization/split_member.hpp>
41 
42 #include <cstdio>
43 #include <cstring>
44 #include <iostream>
45 #include <vector>
46 
47 DEAL_II_NAMESPACE_OPEN
48 
49 
50 #ifdef DEAL_II_WITH_PETSC
51 namespace PETScWrappers
52 {
53  class VectorBase;
54 }
55 #endif
56 
57 #ifdef DEAL_II_WITH_TRILINOS
58 namespace TrilinosWrappers
59 {
60  namespace MPI
61  {
62  class Vector;
63  }
64 } // namespace TrilinosWrappers
65 #endif
66 
67 template <typename number>
68 class LAPACKFullMatrix;
69 
70 template <typename>
71 class BlockVector;
72 
73 template <typename>
74 class VectorView;
75 
76 namespace parallel
77 {
78  namespace internal
79  {
80  class TBBPartitioner;
81  }
82 } // namespace parallel
83 
84 
85 
109 template <typename Number>
110 class Vector : public Subscriptor
111 {
112 public:
113  // The assertion in vector.templates.h for whether or not a number is
114  // finite is not compatible for AD number types.
115  static_assert(
117  "The Vector class does not support auto-differentiable numbers.");
118 
123  using value_type = Number;
124  using pointer = value_type *;
125  using const_pointer = const value_type *;
126  using iterator = value_type *;
127  using const_iterator = const value_type *;
128  using reference = value_type &;
129  using const_reference = const value_type &;
130  using size_type = types::global_dof_index;
131 
142 
143 public:
151  Vector();
152 
160  Vector(const Vector<Number> &v);
161 
166  Vector(Vector<Number> &&v) noexcept;
167 
180  template <typename OtherNumber>
181  explicit Vector(const Vector<OtherNumber> &v);
182 
183 #ifdef DEAL_II_WITH_PETSC
184 
195  explicit Vector(const PETScWrappers::VectorBase &v);
196 #endif
197 
198 #ifdef DEAL_II_WITH_TRILINOS
199 
213  explicit Vector(const TrilinosWrappers::MPI::Vector &v);
214 #endif
215 
225  explicit Vector(const size_type n);
226 
231  template <typename InputIterator>
232  Vector(const InputIterator first, const InputIterator last);
233 
238  virtual ~Vector() override = default;
239 
253  void
254  compress(::VectorOperation::values operation =
255  ::VectorOperation::unknown) const;
256 
272  virtual void
273  reinit(const size_type N, const bool omit_zeroing_entries = false);
274 
298  void
299  grow_or_shrink(const size_type N);
300 
308  template <typename Number2>
309  void
310  reinit(const Vector<Number2> &V, const bool omit_zeroing_entries = false);
311 
327  virtual void
328  swap(Vector<Number> &v);
329 
343  Vector<Number> &
344  operator=(const Number s);
345 
351  Vector<Number> &
352  operator=(const Vector<Number> &v);
353 
359  Vector<Number> &
360  operator=(Vector<Number> &&v) noexcept;
361 
367  template <typename Number2>
368  Vector<Number> &
369  operator=(const Vector<Number2> &v);
370 
374  Vector<Number> &
375  operator=(const BlockVector<Number> &v);
376 
377 #ifdef DEAL_II_WITH_PETSC
378 
389  Vector<Number> &
391 #endif
392 
393 
394 #ifdef DEAL_II_WITH_TRILINOS
395 
410  Vector<Number> &
412 #endif
413 
419  template <typename Number2>
420  bool
421  operator==(const Vector<Number2> &v) const;
422 
428  template <typename Number2>
429  bool
430  operator!=(const Vector<Number2> &v) const;
431 
433 
434 
439 
453  template <typename Number2>
454  Number operator*(const Vector<Number2> &V) const;
455 
463  real_type
464  norm_sqr() const;
465 
473  Number
474  mean_value() const;
475 
483  real_type
484  l1_norm() const;
485 
494  real_type
495  l2_norm() const;
496 
505  real_type
506  lp_norm(const real_type p) const;
507 
511  real_type
512  linfty_norm() const;
513 
538  Number
539  add_and_dot(const Number a, const Vector<Number> &V, const Vector<Number> &W);
540 
542 
543 
548 
554  iterator
555  begin();
556 
560  const_iterator
561  begin() const;
562 
566  iterator
567  end();
568 
573  const_iterator
574  end() const;
575 
579  Number
580  operator()(const size_type i) const;
581 
585  Number &
586  operator()(const size_type i);
587 
593  Number operator[](const size_type i) const;
594 
600  Number &operator[](const size_type i);
601 
617  template <typename OtherNumber>
618  void
619  extract_subvector_to(const std::vector<size_type> &indices,
620  std::vector<OtherNumber> & values) const;
621 
649  template <typename ForwardIterator, typename OutputIterator>
650  void
651  extract_subvector_to(ForwardIterator indices_begin,
652  const ForwardIterator indices_end,
653  OutputIterator values_begin) const;
655 
656 
661 
667  Vector<Number> &
668  operator+=(const Vector<Number> &V);
669 
675  Vector<Number> &
676  operator-=(const Vector<Number> &V);
677 
682  template <typename OtherNumber>
683  void
684  add(const std::vector<size_type> & indices,
685  const std::vector<OtherNumber> &values);
686 
691  template <typename OtherNumber>
692  void
693  add(const std::vector<size_type> &indices, const Vector<OtherNumber> &values);
694 
700  template <typename OtherNumber>
701  void
702  add(const size_type n_elements,
703  const size_type * indices,
704  const OtherNumber *values);
705 
712  void
713  add(const Number s);
714 
720  void
721  add(const Number a,
722  const Vector<Number> &V,
723  const Number b,
724  const Vector<Number> &W);
725 
731  void
732  add(const Number a, const Vector<Number> &V);
733 
739  void
740  sadd(const Number s, const Vector<Number> &V);
741 
747  void
748  sadd(const Number s, const Number a, const Vector<Number> &V);
749 
755  Vector<Number> &
756  operator*=(const Number factor);
757 
763  Vector<Number> &
764  operator/=(const Number factor);
765 
773  void
774  scale(const Vector<Number> &scaling_factors);
775 
781  template <typename Number2>
782  void
783  scale(const Vector<Number2> &scaling_factors);
784 
790  void
791  equ(const Number a, const Vector<Number> &u);
792 
796  template <typename Number2>
797  void
798  equ(const Number a, const Vector<Number2> &u);
799 
812  DEAL_II_DEPRECATED
813  void
814  ratio(const Vector<Number> &a, const Vector<Number> &b);
815 
820  void
821  update_ghost_values() const;
823 
824 
835  DEAL_II_DEPRECATED
836  void
837  print(const char *format = nullptr) const;
838 
845  void
846  print(std::ostream & out,
847  const unsigned int precision = 3,
848  const bool scientific = true,
849  const bool across = true) const;
850 
859  DEAL_II_DEPRECATED
860  void
861  print(LogStream & out,
862  const unsigned int width = 6,
863  const bool across = true) const;
864 
870  void
871  block_write(std::ostream &out) const;
872 
884  void
885  block_read(std::istream &in);
886 
891  template <class Archive>
892  void
893  save(Archive &ar, const unsigned int version) const;
894 
899  template <class Archive>
900  void
901  load(Archive &ar, const unsigned int version);
902 
903  BOOST_SERIALIZATION_SPLIT_MEMBER()
904 
905 
913 
919  bool
920  in_local_range(const size_type global_index) const;
921 
937  IndexSet
938  locally_owned_elements() const;
939 
943  size_type
944  size() const;
945 
951  bool
952  all_zero() const;
953 
963  bool
964  is_non_negative() const;
965 
970  std::size_t
971  memory_consumption() const;
973 
974 protected:
979  size_type vec_size;
980 
987  size_type max_vec_size;
988 
996  std::unique_ptr<Number[], decltype(&free)> values;
997 
1002  mutable std::shared_ptr<parallel::internal::TBBPartitioner>
1004 
1008  template <typename Number2>
1009  friend class Vector;
1010 
1014  template <typename Number2>
1015  friend class LAPACKFullMatrix;
1016 
1020  friend class VectorView<Number>;
1021 
1022 private:
1028  void
1029  allocate(const size_type copy_n_el = 0);
1030 };
1031 
1033 /*----------------------- Inline functions ----------------------------------*/
1034 
1035 
1036 #ifndef DOXYGEN
1037 
1038 
1039 //------------------------ declarations for explicit specializations
1040 template <>
1042 Vector<int>::lp_norm(const real_type) const;
1043 
1044 
1045 //------------------------ inline functions
1046 
1047 template <typename Number>
1048 inline Vector<Number>::Vector()
1049  : vec_size(0)
1050  , max_vec_size(0)
1051  , values(nullptr, &free)
1052 {
1053  // virtual functions called in constructors and destructors never use the
1054  // override in a derived class
1055  // for clarity be explicit on which function is called
1057 }
1058 
1059 
1060 
1061 template <typename Number>
1062 template <typename InputIterator>
1063 Vector<Number>::Vector(const InputIterator first, const InputIterator last)
1064  : vec_size(0)
1065  , max_vec_size(0)
1066  , values(nullptr, &free)
1067 {
1068  // allocate memory. do not initialize it, as we will copy over to it in a
1069  // second
1070  reinit(std::distance(first, last), true);
1071  std::copy(first, last, begin());
1072 }
1073 
1074 
1075 
1076 template <typename Number>
1077 inline Vector<Number>::Vector(const size_type n)
1078  : vec_size(0)
1079  , max_vec_size(0)
1080  , values(nullptr, &free)
1081 {
1082  // virtual functions called in constructors and destructors never use the
1083  // override in a derived class
1084  // for clarity be explicit on which function is called
1085  Vector<Number>::reinit(n, false);
1086 }
1087 
1088 
1089 
1090 template <typename Number>
1091 inline typename Vector<Number>::size_type
1092 Vector<Number>::size() const
1093 {
1094  return vec_size;
1095 }
1096 
1097 
1098 template <typename Number>
1099 inline bool
1100 Vector<Number>::in_local_range(const size_type) const
1101 {
1102  return true;
1103 }
1104 
1105 
1106 
1107 template <typename Number>
1108 inline typename Vector<Number>::iterator
1110 {
1111  return values.get();
1112 }
1113 
1114 
1115 
1116 template <typename Number>
1117 inline typename Vector<Number>::const_iterator
1118 Vector<Number>::begin() const
1119 {
1120  return values.get();
1121 }
1122 
1123 
1124 
1125 template <typename Number>
1126 inline typename Vector<Number>::iterator
1128 {
1129  return values.get() + vec_size;
1130 }
1131 
1132 
1133 
1134 template <typename Number>
1135 inline typename Vector<Number>::const_iterator
1136 Vector<Number>::end() const
1137 {
1138  return values.get() + vec_size;
1139 }
1140 
1141 
1142 
1143 template <typename Number>
1144 inline Number
1145 Vector<Number>::operator()(const size_type i) const
1146 {
1147  Assert(i < vec_size, ExcIndexRange(i, 0, vec_size));
1148  return values[i];
1149 }
1150 
1151 
1152 
1153 template <typename Number>
1154 inline Number &
1155 Vector<Number>::operator()(const size_type i)
1156 {
1157  Assert(i < vec_size, ExcIndexRangeType<size_type>(i, 0, vec_size));
1158  return values[i];
1159 }
1160 
1161 
1162 
1163 template <typename Number>
1164 inline Number Vector<Number>::operator[](const size_type i) const
1165 {
1166  return operator()(i);
1167 }
1168 
1169 
1170 
1171 template <typename Number>
1172 inline Number &Vector<Number>::operator[](const size_type i)
1173 {
1174  return operator()(i);
1175 }
1176 
1177 
1178 
1179 template <typename Number>
1180 template <typename OtherNumber>
1181 inline void
1182 Vector<Number>::extract_subvector_to(const std::vector<size_type> &indices,
1183  std::vector<OtherNumber> & values) const
1184 {
1185  for (size_type i = 0; i < indices.size(); ++i)
1186  values[i] = operator()(indices[i]);
1187 }
1188 
1189 
1190 
1191 template <typename Number>
1192 template <typename ForwardIterator, typename OutputIterator>
1193 inline void
1194 Vector<Number>::extract_subvector_to(ForwardIterator indices_begin,
1195  const ForwardIterator indices_end,
1196  OutputIterator values_begin) const
1197 {
1198  while (indices_begin != indices_end)
1199  {
1200  *values_begin = operator()(*indices_begin);
1201  indices_begin++;
1202  values_begin++;
1203  }
1204 }
1205 
1206 
1207 
1208 template <typename Number>
1209 inline Vector<Number> &
1210 Vector<Number>::operator/=(const Number factor)
1211 {
1212  AssertIsFinite(factor);
1213  Assert(factor != Number(0.), ExcZero());
1214 
1215  this->operator*=(Number(1.) / factor);
1216  return *this;
1217 }
1218 
1219 
1220 
1221 template <typename Number>
1222 template <typename OtherNumber>
1223 inline void
1224 Vector<Number>::add(const std::vector<size_type> & indices,
1225  const std::vector<OtherNumber> &values)
1226 {
1227  Assert(indices.size() == values.size(),
1228  ExcDimensionMismatch(indices.size(), values.size()));
1229  add(indices.size(), indices.data(), values.data());
1230 }
1231 
1232 
1233 
1234 template <typename Number>
1235 template <typename OtherNumber>
1236 inline void
1237 Vector<Number>::add(const std::vector<size_type> &indices,
1238  const Vector<OtherNumber> & values)
1239 {
1240  Assert(indices.size() == values.size(),
1241  ExcDimensionMismatch(indices.size(), values.size()));
1242  add(indices.size(), indices.data(), values.values.get());
1243 }
1244 
1245 
1246 
1247 template <typename Number>
1248 template <typename OtherNumber>
1249 inline void
1250 Vector<Number>::add(const size_type n_indices,
1251  const size_type * indices,
1252  const OtherNumber *values)
1253 {
1254  for (size_type i = 0; i < n_indices; ++i)
1255  {
1256  Assert(indices[i] < vec_size, ExcIndexRange(indices[i], 0, vec_size));
1257  Assert(
1258  numbers::is_finite(values[i]),
1259  ExcMessage(
1260  "The given value is not finite but either infinite or Not A Number (NaN)"));
1261 
1262  this->values[indices[i]] += values[i];
1263  }
1264 }
1265 
1266 
1267 
1268 template <typename Number>
1269 template <typename Number2>
1270 inline bool
1272 {
1273  return !(*this == v);
1274 }
1275 
1276 
1277 
1278 template <typename Number>
1280 {}
1281 
1282 
1283 
1284 template <typename Number>
1285 inline void
1287 {}
1288 
1289 
1290 
1291 // Moved from vector.templates.h as an inline function by Luca Heltai
1292 // on 2009/04/12 to prevent strange compiling errors, after making
1293 // swap virtual.
1294 template <typename Number>
1295 inline void
1297 {
1298  std::swap(vec_size, v.vec_size);
1299  std::swap(max_vec_size, v.max_vec_size);
1300  std::swap(values, v.values);
1301 }
1302 
1303 
1304 
1305 template <typename Number>
1306 template <class Archive>
1307 inline void
1308 Vector<Number>::save(Archive &ar, const unsigned int) const
1309 {
1310  // forward to serialization function in the base class.
1311  ar &static_cast<const Subscriptor &>(*this);
1312 
1313  ar &vec_size &max_vec_size;
1314  ar & boost::serialization::make_array(values.get(), max_vec_size);
1315 }
1316 
1317 
1318 
1319 template <typename Number>
1320 template <class Archive>
1321 inline void
1322 Vector<Number>::load(Archive &ar, const unsigned int)
1323 {
1324  // get rid of previous content
1325  values.reset();
1326 
1327  // the load stuff again from the archive
1328  ar &static_cast<Subscriptor &>(*this);
1329  ar &vec_size &max_vec_size;
1330 
1331  allocate();
1332  ar &boost::serialization::make_array(values.get(), max_vec_size);
1333 }
1334 
1335 #endif
1336 
1337 
1351 template <typename Number>
1352 inline void
1354 {
1355  u.swap(v);
1356 }
1357 
1358 
1362 template <typename number>
1363 inline std::ostream &
1364 operator<<(std::ostream &os, const Vector<number> &v)
1365 {
1366  v.print(os);
1367  return os;
1368 }
1369 
1373 template <typename number>
1374 inline LogStream &
1375 operator<<(LogStream &os, const Vector<number> &v)
1376 {
1377  v.print(os);
1378  return os;
1379 }
1380 
1389 template <typename Number>
1390 struct is_serial_vector<Vector<Number>> : std::true_type
1391 {};
1392 
1393 
1394 DEAL_II_NAMESPACE_CLOSE
1395 
1396 #endif
void block_read(std::istream &in)
real_type lp_norm(const real_type p) const
real_type norm_sqr() const
IndexSet locally_owned_elements() const
real_type l2_norm() const
Number add_and_dot(const Number a, const Vector< Number > &V, const Vector< Number > &W)
void ratio(const Vector< Number > &a, const Vector< Number > &b)
void load(Archive &ar, const unsigned int version)
Number mean_value() const
void sadd(const Number s, const Vector< Number > &V)
typename numbers::NumberTraits< typename VectorType::value_type >::real_type real_type
Definition: vector.h:141
bool is_non_negative() const
size_type size() const
void equ(const Number a, const Vector< Number > &u)
bool in_local_range(const size_type global_index) const
STL namespace.
void compress(::VectorOperation::values operation=::VectorOperation::unknown) const
void update_ghost_values() const
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
iterator end()
Vector< Number > & operator=(const Number s)
bool is_finite(const double x)
Definition: numbers.h:429
void allocate(const size_type copy_n_el=0)
Number operator*(const Vector< Number2 > &V) const
size_type max_vec_size
Definition: vector.h:987
bool all_zero() const
real_type linfty_norm() const
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
bool operator==(const Vector< Number2 > &v) const
static::ExceptionBase & ExcMessage(std::string arg1)
void swap(BlockIndices &u, BlockIndices &v)
typename VectorType::value_type value_type
Definition: vector.h:123
void scale(const Vector< Number > &scaling_factors)
Vector< Number > & operator+=(const Vector< Number > &V)
#define Assert(cond, exc)
Definition: exceptions.h:1239
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Number operator[](const size_type i) const
void grow_or_shrink(const size_type N)
bool operator!=(const Vector< Number2 > &v) const
Vector< Number > & operator-=(const Vector< Number > &V)
void block_write(std::ostream &out) const
Vector< Number > & operator*=(const Number factor)
iterator begin()
real_type l1_norm() const
std::unique_ptr< Number[], decltype(&free)> values
Definition: vector.h:996
virtual ~Vector() override=default
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1353
unsigned int global_dof_index
Definition: types.h:89
virtual void swap(Vector< Number > &v)
Vector< Number > & operator/=(const Number factor)
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< OtherNumber > &values) const
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
size_type vec_size
Definition: vector.h:979
std::shared_ptr< parallel::internal::TBBPartitioner > thread_loop_partitioner
Definition: vector.h:1003
void save(Archive &ar, const unsigned int version) const
void free(T *&pointer)
Definition: cuda.h:95
Number operator()(const size_type i) const
static::ExceptionBase & ExcZero()
void print(const char *format=nullptr) const
#define AssertIsFinite(number)
Definition: exceptions.h:1428
std::size_t memory_consumption() const