Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
vector.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 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_vector_h
17 #define dealii_vector_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/aligned_vector.h>
23 #include <deal.II/base/exceptions.h>
24 #include <deal.II/base/index_set.h>
25 #include <deal.II/base/logstream.h>
26 #include <deal.II/base/subscriptor.h>
27 
28 #include <deal.II/differentiation/ad/ad_number_traits.h>
29 
30 #include <deal.II/lac/vector_operation.h>
31 #include <deal.II/lac/vector_type_traits.h>
32 
33 #include <boost/serialization/split_member.hpp>
34 
35 #include <algorithm>
36 #include <initializer_list>
37 #include <iosfwd>
38 #include <iterator>
39 #include <vector>
40 
41 DEAL_II_NAMESPACE_OPEN
42 
43 
44 #ifdef DEAL_II_WITH_PETSC
45 namespace PETScWrappers
46 {
47  class VectorBase;
48 }
49 #endif
50 
51 #ifdef DEAL_II_WITH_TRILINOS
52 namespace TrilinosWrappers
53 {
54  namespace MPI
55  {
56  class Vector;
57  }
58 } // namespace TrilinosWrappers
59 #endif
60 
61 template <typename number>
62 class LAPACKFullMatrix;
63 
64 template <typename>
65 class BlockVector;
66 
67 namespace parallel
68 {
69  namespace internal
70  {
71  class TBBPartitioner;
72  }
73 } // namespace parallel
74 
75 
76 
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 
150  Vector();
151 
161  Vector(const Vector<Number> &v);
162 
167  Vector(Vector<Number> &&v) noexcept = default;
168 
177  template <typename OtherNumber>
178  explicit Vector(const Vector<OtherNumber> &v);
179 
195  template <typename OtherNumber>
196  explicit Vector(const std::initializer_list<OtherNumber> &v);
197 
198 #ifdef DEAL_II_WITH_PETSC
199 
210  explicit Vector(const PETScWrappers::VectorBase &v);
211 #endif
212 
213 #ifdef DEAL_II_WITH_TRILINOS
214 
228  explicit Vector(const TrilinosWrappers::MPI::Vector &v);
229 #endif
230 
240  explicit Vector(const size_type n);
241 
246  template <typename InputIterator>
247  Vector(const InputIterator first, const InputIterator last);
248 
253  virtual ~Vector() override = default;
254 
268  void
269  compress(::VectorOperation::values operation =
270  ::VectorOperation::unknown) const;
271 
287  virtual void
288  reinit(const size_type N, const bool omit_zeroing_entries = false);
289 
313  void
314  grow_or_shrink(const size_type N);
315 
323  void
324  apply_givens_rotation(const std::array<Number, 3> &csr,
325  const size_type i,
326  const size_type k);
327 
335  template <typename Number2>
336  void
337  reinit(const Vector<Number2> &V, const bool omit_zeroing_entries = false);
338 
354  virtual void
355  swap(Vector<Number> &v);
356 
367  Vector<Number> &
368  operator=(const Number s);
369 
375  Vector<Number> &
376  operator=(const Vector<Number> &v);
377 
383  Vector<Number> &
384  operator=(Vector<Number> &&v) noexcept = default;
385 
391  template <typename Number2>
392  Vector<Number> &
393  operator=(const Vector<Number2> &v);
394 
398  Vector<Number> &
399  operator=(const BlockVector<Number> &v);
400 
401 #ifdef DEAL_II_WITH_PETSC
402 
413  Vector<Number> &
414  operator=(const PETScWrappers::VectorBase &v);
415 #endif
416 
417 
418 #ifdef DEAL_II_WITH_TRILINOS
419 
434  Vector<Number> &
435  operator=(const TrilinosWrappers::MPI::Vector &v);
436 #endif
437 
443  template <typename Number2>
444  bool
445  operator==(const Vector<Number2> &v) const;
446 
452  template <typename Number2>
453  bool
454  operator!=(const Vector<Number2> &v) const;
455 
457 
458 
463 
477  template <typename Number2>
478  Number operator*(const Vector<Number2> &V) const;
479 
487  real_type
488  norm_sqr() const;
489 
497  Number
498  mean_value() const;
499 
507  real_type
508  l1_norm() const;
509 
518  real_type
519  l2_norm() const;
520 
529  real_type
530  lp_norm(const real_type p) const;
531 
535  real_type
536  linfty_norm() const;
537 
562  Number
563  add_and_dot(const Number a, const Vector<Number> &V, const Vector<Number> &W);
564 
566 
567 
572 
576  pointer
577  data();
578 
582  const_pointer
583  data() const;
584 
590  iterator
591  begin();
592 
596  const_iterator
597  begin() const;
598 
602  iterator
603  end();
604 
609  const_iterator
610  end() const;
611 
615  Number
616  operator()(const size_type i) const;
617 
621  Number &
622  operator()(const size_type i);
623 
629  Number operator[](const size_type i) const;
630 
636  Number &operator[](const size_type i);
637 
653  template <typename OtherNumber>
654  void
655  extract_subvector_to(const std::vector<size_type> &indices,
656  std::vector<OtherNumber> & values) const;
657 
685  template <typename ForwardIterator, typename OutputIterator>
686  void
687  extract_subvector_to(ForwardIterator indices_begin,
688  const ForwardIterator indices_end,
689  OutputIterator values_begin) const;
691 
692 
697 
703  Vector<Number> &
704  operator+=(const Vector<Number> &V);
705 
711  Vector<Number> &
712  operator-=(const Vector<Number> &V);
713 
718  template <typename OtherNumber>
719  void
720  add(const std::vector<size_type> & indices,
721  const std::vector<OtherNumber> &values);
722 
727  template <typename OtherNumber>
728  void
729  add(const std::vector<size_type> &indices, const Vector<OtherNumber> &values);
730 
736  template <typename OtherNumber>
737  void
738  add(const size_type n_elements,
739  const size_type * indices,
740  const OtherNumber *values);
741 
748  void
749  add(const Number s);
750 
756  void
757  add(const Number a,
758  const Vector<Number> &V,
759  const Number b,
760  const Vector<Number> &W);
761 
767  void
768  add(const Number a, const Vector<Number> &V);
769 
775  void
776  sadd(const Number s, const Vector<Number> &V);
777 
783  void
784  sadd(const Number s, const Number a, const Vector<Number> &V);
785 
791  Vector<Number> &
792  operator*=(const Number factor);
793 
799  Vector<Number> &
800  operator/=(const Number factor);
801 
809  void
810  scale(const Vector<Number> &scaling_factors);
811 
817  template <typename Number2>
818  void
819  scale(const Vector<Number2> &scaling_factors);
820 
826  void
827  equ(const Number a, const Vector<Number> &u);
828 
832  template <typename Number2>
833  void
834  equ(const Number a, const Vector<Number2> &u);
835 
848  DEAL_II_DEPRECATED
849  void
850  ratio(const Vector<Number> &a, const Vector<Number> &b);
851 
856  void
857  update_ghost_values() const;
859 
860 
871  DEAL_II_DEPRECATED
872  void
873  print(const char *format = nullptr) const;
874 
881  void
882  print(std::ostream & out,
883  const unsigned int precision = 3,
884  const bool scientific = true,
885  const bool across = true) const;
886 
895  DEAL_II_DEPRECATED
896  void
897  print(LogStream & out,
898  const unsigned int width = 6,
899  const bool across = true) const;
900 
906  void
907  block_write(std::ostream &out) const;
908 
920  void
921  block_read(std::istream &in);
922 
927  template <class Archive>
928  void
929  save(Archive &ar, const unsigned int version) const;
930 
935  template <class Archive>
936  void
937  load(Archive &ar, const unsigned int version);
938 
939  BOOST_SERIALIZATION_SPLIT_MEMBER()
940 
941 
949 
955  bool
956  in_local_range(const size_type global_index) const;
957 
973  IndexSet
974  locally_owned_elements() const;
975 
979  size_type
980  size() const;
981 
987  bool
988  all_zero() const;
989 
999  bool
1000  is_non_negative() const;
1001 
1006  std::size_t
1007  memory_consumption() const;
1009 
1010 private:
1014  AlignedVector<Number> values;
1015 
1021  void
1022  maybe_reset_thread_partitioner();
1023 
1027  void
1028  do_reinit(const size_type new_size,
1029  const bool omit_zeroing_entries,
1030  const bool reset_partitioner);
1031 
1036  mutable std::shared_ptr<parallel::internal::TBBPartitioner>
1037  thread_loop_partitioner;
1038 
1042  template <typename Number2>
1043  friend class Vector;
1044 };
1045 
1047 /*----------------------- Inline functions ----------------------------------*/
1048 
1049 
1050 #ifndef DOXYGEN
1051 
1052 
1053 //------------------------ declarations for explicit specializations
1054 template <>
1056 Vector<int>::lp_norm(const real_type) const;
1057 
1058 
1059 //------------------------ inline functions
1060 
1061 template <typename Number>
1062 inline Vector<Number>::Vector()
1063 {
1064  // virtual functions called in constructors and destructors never use the
1065  // override in a derived class
1066  // for clarity be explicit on which function is called
1068 }
1069 
1070 
1071 
1072 template <typename Number>
1073 template <typename OtherNumber>
1074 Vector<Number>::Vector(const std::initializer_list<OtherNumber> &v)
1075  : Vector(v.begin(), v.end())
1076 {}
1077 
1078 
1079 
1080 template <typename Number>
1081 template <typename InputIterator>
1082 Vector<Number>::Vector(const InputIterator first, const InputIterator last)
1083 {
1084  // allocate memory. do not initialize it, as we will copy over to it in a
1085  // second
1086  reinit(std::distance(first, last), true);
1087  std::copy(first, last, begin());
1088 }
1089 
1090 
1091 
1092 template <typename Number>
1093 inline Vector<Number>::Vector(const size_type n)
1094 {
1095  // virtual functions called in constructors and destructors never use the
1096  // override in a derived class
1097  // for clarity be explicit on which function is called
1098  Vector<Number>::reinit(n, false);
1099 }
1100 
1101 
1102 
1103 template <typename Number>
1104 inline typename Vector<Number>::size_type
1105 Vector<Number>::size() const
1106 {
1107  return values.size();
1108 }
1109 
1110 
1111 template <typename Number>
1112 inline bool
1113 Vector<Number>::in_local_range(const size_type) const
1114 {
1115  return true;
1116 }
1117 
1118 
1119 
1120 template <typename Number>
1121 inline typename Vector<Number>::pointer
1123 {
1124  return values.data();
1125 }
1126 
1127 
1128 
1129 template <typename Number>
1130 inline typename Vector<Number>::const_pointer
1131 Vector<Number>::data() const
1132 {
1133  return values.data();
1134 }
1135 
1136 
1137 
1138 template <typename Number>
1139 inline typename Vector<Number>::iterator
1141 {
1142  return values.begin();
1143 }
1144 
1145 
1146 
1147 template <typename Number>
1148 inline typename Vector<Number>::const_iterator
1149 Vector<Number>::begin() const
1150 {
1151  return values.begin();
1152 }
1153 
1154 
1155 
1156 template <typename Number>
1157 inline typename Vector<Number>::iterator
1159 {
1160  return values.end();
1161 }
1162 
1163 
1164 
1165 template <typename Number>
1166 inline typename Vector<Number>::const_iterator
1167 Vector<Number>::end() const
1168 {
1169  return values.end();
1170 }
1171 
1172 
1173 
1174 template <typename Number>
1175 inline Number
1176 Vector<Number>::operator()(const size_type i) const
1177 {
1178  Assert(i < size(), ExcIndexRange(i, 0, size()));
1179  return values[i];
1180 }
1181 
1182 
1183 
1184 template <typename Number>
1185 inline Number &
1186 Vector<Number>::operator()(const size_type i)
1187 {
1188  Assert(i < size(), ExcIndexRangeType<size_type>(i, 0, size()));
1189  return values[i];
1190 }
1191 
1192 
1193 
1194 template <typename Number>
1195 inline Number Vector<Number>::operator[](const size_type i) const
1196 {
1197  return operator()(i);
1198 }
1199 
1200 
1201 
1202 template <typename Number>
1203 inline Number &Vector<Number>::operator[](const size_type i)
1204 {
1205  return operator()(i);
1206 }
1207 
1208 
1209 
1210 template <typename Number>
1211 template <typename OtherNumber>
1212 inline void
1213 Vector<Number>::extract_subvector_to(const std::vector<size_type> &indices,
1214  std::vector<OtherNumber> & values) const
1215 {
1216  for (size_type i = 0; i < indices.size(); ++i)
1217  values[i] = operator()(indices[i]);
1218 }
1219 
1220 
1221 
1222 template <typename Number>
1223 template <typename ForwardIterator, typename OutputIterator>
1224 inline void
1225 Vector<Number>::extract_subvector_to(ForwardIterator indices_begin,
1226  const ForwardIterator indices_end,
1227  OutputIterator values_begin) const
1228 {
1229  while (indices_begin != indices_end)
1230  {
1231  *values_begin = operator()(*indices_begin);
1232  indices_begin++;
1233  values_begin++;
1234  }
1235 }
1236 
1237 
1238 
1239 template <typename Number>
1240 inline Vector<Number> &
1241 Vector<Number>::operator/=(const Number factor)
1242 {
1243  AssertIsFinite(factor);
1244  Assert(factor != Number(0.), ExcZero());
1245 
1246  this->operator*=(Number(1.) / factor);
1247  return *this;
1248 }
1249 
1250 
1251 
1252 template <typename Number>
1253 template <typename OtherNumber>
1254 inline void
1255 Vector<Number>::add(const std::vector<size_type> & indices,
1256  const std::vector<OtherNumber> &values)
1257 {
1258  Assert(indices.size() == values.size(),
1259  ExcDimensionMismatch(indices.size(), values.size()));
1260  add(indices.size(), indices.data(), values.data());
1261 }
1262 
1263 
1264 
1265 template <typename Number>
1266 template <typename OtherNumber>
1267 inline void
1268 Vector<Number>::add(const std::vector<size_type> &indices,
1269  const Vector<OtherNumber> & values)
1270 {
1271  Assert(indices.size() == values.size(),
1272  ExcDimensionMismatch(indices.size(), values.size()));
1273  add(indices.size(), indices.data(), values.values.begin());
1274 }
1275 
1276 
1277 
1278 template <typename Number>
1279 template <typename OtherNumber>
1280 inline void
1281 Vector<Number>::add(const size_type n_indices,
1282  const size_type * indices,
1283  const OtherNumber *values)
1284 {
1285  for (size_type i = 0; i < n_indices; ++i)
1286  {
1287  Assert(indices[i] < size(), ExcIndexRange(indices[i], 0, size()));
1288  Assert(
1289  numbers::is_finite(values[i]),
1290  ExcMessage(
1291  "The given value is not finite but either infinite or Not A Number (NaN)"));
1292 
1293  this->values[indices[i]] += values[i];
1294  }
1295 }
1296 
1297 
1298 
1299 template <typename Number>
1300 template <typename Number2>
1301 inline bool
1302 Vector<Number>::operator!=(const Vector<Number2> &v) const
1303 {
1304  return !(*this == v);
1305 }
1306 
1307 
1308 
1309 template <typename Number>
1311 {}
1312 
1313 
1314 
1315 template <typename Number>
1316 inline void
1318 {}
1319 
1320 
1321 
1322 // Moved from vector.templates.h as an inline function by Luca Heltai
1323 // on 2009/04/12 to prevent strange compiling errors, after making
1324 // swap virtual.
1325 template <typename Number>
1326 inline void
1327 Vector<Number>::swap(Vector<Number> &v)
1328 {
1329  values.swap(v.values);
1330  std::swap(thread_loop_partitioner, v.thread_loop_partitioner);
1331 }
1332 
1333 
1334 
1335 template <typename Number>
1336 template <class Archive>
1337 inline void
1338 Vector<Number>::save(Archive &ar, const unsigned int) const
1339 {
1340  // forward to serialization function in the base class.
1341  ar &static_cast<const Subscriptor &>(*this);
1342  ar &values;
1343 }
1344 
1345 
1346 
1347 template <typename Number>
1348 template <class Archive>
1349 inline void
1350 Vector<Number>::load(Archive &ar, const unsigned int)
1351 {
1352  // the load stuff again from the archive
1353  ar &static_cast<Subscriptor &>(*this);
1354  ar &values;
1356 }
1357 
1358 #endif
1359 
1360 
1374 template <typename Number>
1375 inline void
1376 swap(Vector<Number> &u, Vector<Number> &v)
1377 {
1378  u.swap(v);
1379 }
1380 
1381 
1385 template <typename number>
1386 inline std::ostream &
1387 operator<<(std::ostream &os, const Vector<number> &v)
1388 {
1389  v.print(os);
1390  return os;
1391 }
1392 
1396 template <typename number>
1397 inline LogStream &
1398 operator<<(LogStream &os, const Vector<number> &v)
1399 {
1400  v.print(os);
1401  return os;
1402 }
1403 
1412 template <typename Number>
1413 struct is_serial_vector<Vector<Number>> : std::true_type
1414 {};
1415 
1416 
1417 DEAL_II_NAMESPACE_CLOSE
1418 
1419 #endif
void update_ghost_values() const
Number operator[](const size_type i) const
bool operator!=(const Vector< Number2 > &v) const
void load(Archive &ar, const unsigned int version)
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< OtherNumber > &values) const
Number operator()(const size_type i) const
typename numbers::NumberTraits< typename VectorType::value_type >::real_type real_type
Definition: vector.h:141
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)
STL namespace.
pointer data()
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
iterator end()
LinearAlgebra::distributed::Vector< Number > Vector
bool is_finite(const double x)
Definition: numbers.h:463
__global__ void scale(Number *val, const Number *V_val, const size_type N)
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
static ::ExceptionBase & ExcMessage(std::string arg1)
typename VectorType::value_type value_type
Definition: vector.h:123
bool in_local_range(const size_type global_index) const
void maybe_reset_thread_partitioner()
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
real_type lp_norm(const real_type p) const
Vector< Number > & operator*=(const Number factor)
iterator begin()
__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)
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
unsigned int global_dof_index
Definition: types.h:89
virtual void swap(Vector< Number > &v)
Vector< Number > & operator/=(const Number factor)
void save(Archive &ar, const unsigned int version) const
void compress(::VectorOperation::values operation=::VectorOperation::unknown) const
size_type size() const
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
AlignedVector< Number > values
Definition: vector.h:1014
std::shared_ptr< parallel::internal::TBBPartitioner > thread_loop_partitioner
Definition: vector.h:1037
static ::ExceptionBase & ExcZero()
#define AssertIsFinite(number)
Definition: exceptions.h:1669