MFEM  v3.3.2
Finite element discretization library
vector.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #ifndef MFEM_VECTOR
13 #define MFEM_VECTOR
14 
15 // Data type vector
16 
17 #include "../general/array.hpp"
18 #include "../general/globals.hpp"
19 #ifdef MFEM_USE_SUNDIALS
20 #include <nvector/nvector_serial.h>
21 #endif
22 #include <cmath>
23 #include <iostream>
24 #if defined(_MSC_VER) && (_MSC_VER < 1800)
25 #include <float.h>
26 #define isfinite _finite
27 #endif
28 
29 #ifdef MFEM_USE_MPI
30 #include <mpi.h>
31 #endif
32 
33 namespace mfem
34 {
35 
36 /** Count the number of entries in an array of doubles for which isfinite
37  is false, i.e. the entry is a NaN or +/-Inf. */
38 inline int CheckFinite(const double *v, const int n);
39 
40 /// Vector data type.
41 class Vector
42 {
43 protected:
44 
46  double * data;
47 
48 public:
49 
50  /// Default constructor for Vector. Sets size = 0 and data = NULL.
51  Vector () { allocsize = size = 0; data = 0; }
52 
53  /// Copy constructor. Allocates a new data array and copies the data.
54  Vector(const Vector &);
55 
56  /// @brief Creates vector of size s.
57  /// @warning Entries are not initialized to zero!
58  explicit Vector (int s);
59 
60  /// Creates a vector referencing an array of doubles, owned by someone else.
61  /** The pointer @a _data can be NULL. The data array can be replaced later
62  with SetData(). */
63  Vector (double *_data, int _size)
64  { data = _data; size = _size; allocsize = -size; }
65 
66  /// Reads a vector from multiple files
67  void Load (std::istream ** in, int np, int * dim);
68 
69  /// Load a vector from an input stream.
70  void Load(std::istream &in, int Size);
71 
72  /// Load a vector from an input stream, reading the size from the stream.
73  void Load(std::istream &in) { int s; in >> s; Load (in, s); }
74 
75  /// @brief Resize the vector to size @a s.
76  /** If the new size is less than or equal to Capacity() then the internal
77  data array remains the same. Otherwise, the old array is deleted, if
78  owned, and a new array of size @a s is allocated without copying the
79  previous content of the Vector.
80  @warning In the second case above (new size greater than current one),
81  the vector will allocate new data array, even if it did not own the
82  original data! Also, new entries are not initialized! */
83  void SetSize(int s);
84 
85  /// Set the Vector data.
86  /// @warning This method should be called only when OwnsData() is false.
87  void SetData(double *d) { data = d; }
88 
89  /// Set the Vector data and size.
90  /** The Vector does not assume ownership of the new data. The new size is
91  also used as the new Capacity().
92  @warning This method should be called only when OwnsData() is false.
93  @sa NewDataAndSize(). */
94  void SetDataAndSize(double *d, int s)
95  { data = d; size = s; allocsize = -s; }
96 
97  /// Set the Vector data and size, deleting the old data, if owned.
98  /** The Vector does not assume ownership of the new data. The new size is
99  also used as the new Capacity().
100  @sa SetDataAndSize(). */
101  void NewDataAndSize(double *d, int s)
102  {
103  if (allocsize > 0) { delete [] data; }
104  SetDataAndSize(d, s);
105  }
106 
107  void MakeDataOwner() { allocsize = abs(allocsize); }
108 
109  /// Destroy a vector
110  void Destroy();
111 
112  /// Returns the size of the vector.
113  inline int Size() const { return size; }
114 
115  /// Return the size of the currently allocated data array.
116  /** It is always true that Capacity() >= Size(). */
117  inline int Capacity() const { return abs(allocsize); }
118 
119  // double *GetData() { return data; }
120 
121  inline double *GetData() const { return data; }
122 
123  /// Conversion to `double *`.
124  /** @note This conversion function makes it possible to use [] for indexing
125  in addition to the overloaded operator()(int). */
126  inline operator double *() { return data; }
127 
128  /// Conversion to `const double *`.
129  /** @note This conversion function makes it possible to use [] for indexing
130  in addition to the overloaded operator()(int). */
131  inline operator const double *() const { return data; }
132 
133  inline bool OwnsData() const { return (allocsize > 0); }
134 
135  /// Changes the ownership of the data; after the call the Vector is empty
136  inline void StealData(double **p)
137  { *p = data; data = 0; size = allocsize = 0; }
138 
139  /// Changes the ownership of the data; after the call the Vector is empty
140  inline double *StealData() { double *p; StealData(&p); return p; }
141 
142  /// Access Vector entries. Index i = 0 .. size-1.
143  double & Elem (int i);
144 
145  /// Read only access to Vector entries. Index i = 0 .. size-1.
146  const double & Elem (int i) const;
147 
148  /// Access Vector entries using () for 0-based indexing.
149  /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
150  inline double & operator() (int i);
151 
152  /// Read only access to Vector entries using () for 0-based indexing.
153  /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
154  inline const double & operator() (int i) const;
155 
156  /// Dot product with a `double *` array.
157  double operator*(const double *) const;
158 
159  /// Return the inner-product.
160  double operator*(const Vector &v) const;
161 
162  /// Copy Size() entries from @a v.
163  Vector & operator=(const double *v);
164 
165  /// Redefine '=' for vector = vector.
166  Vector & operator=(const Vector &v);
167 
168  /// Redefine '=' for vector = constant.
169  Vector & operator=(double value);
170 
171  Vector & operator*=(double c);
172 
173  Vector & operator/=(double c);
174 
175  Vector & operator-=(double c);
176 
177  Vector & operator-=(const Vector &v);
178 
179  Vector & operator+=(const Vector &v);
180 
181  /// (*this) += a * Va
182  Vector & Add(const double a, const Vector &Va);
183 
184  /// (*this) = a * x
185  Vector & Set(const double a, const Vector &x);
186 
187  void SetVector (const Vector &v, int offset);
188 
189  /// (*this) = -(*this)
190  void Neg();
191 
192  /// Swap the contents of two Vectors
193  inline void Swap(Vector &other);
194 
195  /// Set v = v1 + v2.
196  friend void add(const Vector &v1, const Vector &v2, Vector &v);
197 
198  /// Set v = v1 + alpha * v2.
199  friend void add(const Vector &v1, double alpha, const Vector &v2, Vector &v);
200 
201  /// z = a * (x + y)
202  friend void add(const double a, const Vector &x, const Vector &y, Vector &z);
203 
204  /// z = a * x + b * y
205  friend void add (const double a, const Vector &x,
206  const double b, const Vector &y, Vector &z);
207 
208  /// Set v = v1 - v2.
209  friend void subtract(const Vector &v1, const Vector &v2, Vector &v);
210 
211  /// z = a * (x - y)
212  friend void subtract(const double a, const Vector &x,
213  const Vector &y, Vector &z);
214 
215  /// v = median(v,lo,hi) entrywise. Implementation assumes lo <= hi.
216  void median(const Vector &lo, const Vector &hi);
217 
218  void GetSubVector(const Array<int> &dofs, Vector &elemvect) const;
219  void GetSubVector(const Array<int> &dofs, double *elem_data) const;
220 
221  /// Set the entries listed in `dofs` to the given `value`.
222  void SetSubVector(const Array<int> &dofs, const double value);
223  void SetSubVector(const Array<int> &dofs, const Vector &elemvect);
224  void SetSubVector(const Array<int> &dofs, double *elem_data);
225 
226  /// Add (element) subvector to the vector.
227  void AddElementVector(const Array<int> & dofs, const Vector & elemvect);
228  void AddElementVector(const Array<int> & dofs, double *elem_data);
229  void AddElementVector(const Array<int> & dofs, const double a,
230  const Vector & elemvect);
231 
232  /// Set all vector entries NOT in the 'dofs' array to the given 'val'.
233  void SetSubVectorComplement(const Array<int> &dofs, const double val);
234 
235  /// Prints vector to stream out.
236  void Print(std::ostream & out = mfem::out, int width = 8) const;
237 
238  /// Prints vector to stream out in HYPRE_Vector format.
239  void Print_HYPRE(std::ostream &out) const;
240 
241  /// Set random values in the vector.
242  void Randomize(int seed = 0);
243  /// Returns the l2 norm of the vector.
244  double Norml2() const;
245  /// Returns the l_infinity norm of the vector.
246  double Normlinf() const;
247  /// Returns the l_1 norm of the vector.
248  double Norml1() const;
249  /// Returns the l_p norm of the vector.
250  double Normlp(double p) const;
251  /// Returns the maximal element of the vector.
252  double Max() const;
253  /// Returns the minimal element of the vector.
254  double Min() const;
255  /// Return the sum of the vector entries
256  double Sum() const;
257  /// Compute the square of the Euclidean distance to another vector.
258  inline double DistanceSquaredTo(const double *p) const;
259  /// Compute the Euclidean distance to another vector.
260  inline double DistanceTo(const double *p) const;
261 
262  /** Count the number of entries in the Vector for which isfinite
263  is false, i.e. the entry is a NaN or +/-Inf. */
264  int CheckFinite() const { return mfem::CheckFinite(data, size); }
265 
266  /// Destroys vector.
267  virtual ~Vector ();
268 
269 #ifdef MFEM_USE_SUNDIALS
270  /// Construct a wrapper Vector from SUNDIALS N_Vector.
271  explicit Vector(N_Vector nv);
272 
273  /// Return a new wrapper SUNDIALS N_Vector of type SUNDIALS_NVEC_SERIAL.
274  /** The returned N_Vector must be destroyed by the caller. */
275  virtual N_Vector ToNVector() { return N_VMake_Serial(Size(), GetData()); }
276 
277  /** @brief Update an existing wrapper SUNDIALS N_Vector to point to this
278  Vector. */
279  virtual void ToNVector(N_Vector &nv);
280 #endif
281 };
282 
283 // Inline methods
284 
285 inline bool IsFinite(const double &val)
286 {
287  // isfinite didn't appear in a standard until C99, and later C++11. It wasn't
288  // standard in C89 or C++98. PGI as of 14.7 still defines it as a macro.
289 #ifdef isfinite
290  return isfinite(val);
291 #else
292  return std::isfinite(val);
293 #endif
294 }
295 
296 inline int CheckFinite(const double *v, const int n)
297 {
298  int bad = 0;
299  for (int i = 0; i < n; i++)
300  {
301  if (!IsFinite(v[i])) { bad++; }
302  }
303  return bad;
304 }
305 
306 inline Vector::Vector (int s)
307 {
308  if (s > 0)
309  {
310  allocsize = size = s;
311  data = new double[s];
312  }
313  else
314  {
315  allocsize = size = 0;
316  data = NULL;
317  }
318 }
319 
320 inline void Vector::SetSize(int s)
321 {
322  if (s == size)
323  {
324  return;
325  }
326  if (s <= abs(allocsize))
327  {
328  size = s;
329  return;
330  }
331  if (allocsize > 0)
332  {
333  delete [] data;
334  }
335  allocsize = size = s;
336  data = new double[s];
337 }
338 
339 inline void Vector::Destroy()
340 {
341  if (allocsize > 0)
342  {
343  delete [] data;
344  }
345  allocsize = size = 0;
346  data = NULL;
347 }
348 
349 inline double & Vector::operator() (int i)
350 {
351  MFEM_ASSERT(data && i >= 0 && i < size,
352  "index [" << i << "] is out of range [0," << size << ")");
353 
354  return data[i];
355 }
356 
357 inline const double & Vector::operator() (int i) const
358 {
359  MFEM_ASSERT(data && i >= 0 && i < size,
360  "index [" << i << "] is out of range [0," << size << ")");
361 
362  return data[i];
363 }
364 
365 inline void Vector::Swap(Vector &other)
366 {
367  mfem::Swap(size, other.size);
369  mfem::Swap(data, other.data);
370 }
371 
372 /// Specialization of the template function Swap<> for class Vector
373 template<> inline void Swap<Vector>(Vector &a, Vector &b)
374 {
375  a.Swap(b);
376 }
377 
379 {
380  if (allocsize > 0)
381  {
382  delete [] data;
383  }
384 }
385 
386 inline double DistanceSquared(const double *x, const double *y, const int n)
387 {
388  double d = 0.0;
389 
390  for (int i = 0; i < n; i++)
391  {
392  d += (x[i]-y[i])*(x[i]-y[i]);
393  }
394 
395  return d;
396 }
397 
398 inline double Distance(const double *x, const double *y, const int n)
399 {
400  return std::sqrt(DistanceSquared(x, y, n));
401 }
402 
403 inline double Vector::DistanceSquaredTo(const double *p) const
404 {
405  return DistanceSquared(data, p, size);
406 }
407 
408 inline double Vector::DistanceTo(const double *p) const
409 {
410  return Distance(data, p, size);
411 }
412 
413 /// Returns the inner product of x and y
414 /** In parallel this computes the inner product of the local vectors,
415  producing different results on each MPI rank.
416 */
417 inline double InnerProduct(const Vector &x, const Vector &y)
418 {
419  return x * y;
420 }
421 
422 #ifdef MFEM_USE_MPI
423 /// Returns the inner product of x and y in parallel
424 /** In parallel this computes the inner product of the global vectors,
425  producing identical results on each MPI rank.
426 */
427 inline double InnerProduct(MPI_Comm comm, const Vector &x, const Vector &y)
428 {
429  double loc_prod = x * y;
430  double glb_prod;
431  MPI_Allreduce(&loc_prod, &glb_prod, 1, MPI_DOUBLE, MPI_SUM, comm);
432  return glb_prod;
433 }
434 #endif
435 
436 }
437 
438 #endif
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:498
void SetVector(const Vector &v, int offset)
Definition: vector.cpp:238
int CheckFinite(const double *v, const int n)
Definition: vector.hpp:296
Vector()
Default constructor for Vector. Sets size = 0 and data = NULL.
Definition: vector.hpp:51
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:101
double & Elem(int i)
Access Vector entries. Index i = 0 .. size-1.
Definition: vector.cpp:80
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:320
void MakeDataOwner()
Definition: vector.hpp:107
friend void subtract(const Vector &v1, const Vector &v2, Vector &v)
Set v = v1 - v2.
Definition: vector.cpp:390
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:670
double & operator()(int i)
Access Vector entries using () for 0-based indexing.
Definition: vector.hpp:349
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Definition: vector.cpp:462
void StealData(double **p)
Changes the ownership of the data; after the call the Vector is empty.
Definition: vector.hpp:136
int Capacity() const
Return the size of the currently allocated data array.
Definition: vector.hpp:117
int Size() const
Returns the size of the vector.
Definition: vector.hpp:113
void Swap< Vector >(Vector &a, Vector &b)
Specialization of the template function Swap<> for class Vector.
Definition: vector.hpp:373
double Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:707
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:651
double * GetData() const
Definition: vector.hpp:121
double operator*(const double *) const
Dot product with a double * array.
Definition: vector.cpp:90
Vector & operator=(const double *v)
Copy Size() entries from v.
Definition: vector.cpp:117
bool IsFinite(const double &val)
Definition: vector.hpp:285
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Definition: vector.cpp:50
bool OwnsData() const
Definition: vector.hpp:133
double Normlp(double p) const
Returns the l_p norm of the vector.
Definition: vector.cpp:727
int dim
Definition: ex3.cpp:47
void Swap(Vector &other)
Swap the contents of two Vectors.
Definition: vector.hpp:365
void median(const Vector &lo, const Vector &hi)
v = median(v,lo,hi) entrywise. Implementation assumes lo <= hi.
Definition: vector.cpp:445
double DistanceSquared(const double *x, const double *y, const int n)
Definition: vector.hpp:386
void Load(std::istream &in)
Load a vector from an input stream, reading the size from the stream.
Definition: vector.hpp:73
double * StealData()
Changes the ownership of the data; after the call the Vector is empty.
Definition: vector.hpp:140
void SetData(double *d)
Definition: vector.hpp:87
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add (element) subvector to the vector.
Definition: vector.cpp:550
Vector & operator/=(double c)
Definition: vector.cpp:156
double DistanceTo(const double *p) const
Compute the Euclidean distance to another vector.
Definition: vector.hpp:408
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the 'dofs' array to the given 'val'.
Definition: vector.cpp:601
Vector(double *_data, int _size)
Creates a vector referencing an array of doubles, owned by someone else.
Definition: vector.hpp:63
Vector & operator*=(double c)
Definition: vector.cpp:147
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:789
double Norml1() const
Returns the l_1 norm of the vector.
Definition: vector.cpp:717
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:494
double Distance(const double *x, const double *y, const int n)
Definition: vector.hpp:398
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:609
void Print_HYPRE(std::ostream &out) const
Prints vector to stream out in HYPRE_Vector format.
Definition: vector.cpp:633
double DistanceSquaredTo(const double *p) const
Compute the square of the Euclidean distance to another vector.
Definition: vector.hpp:403
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:94
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:223
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:203
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:205
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:776
void Destroy()
Destroy a vector.
Definition: vector.hpp:339
Vector & operator-=(double c)
Definition: vector.cpp:166
int CheckFinite() const
Definition: vector.hpp:264
const double alpha
Definition: ex15.cpp:337
Vector data type.
Definition: vector.hpp:41
friend void add(const Vector &v1, const Vector &v2, Vector &v)
Set v = v1 + v2.
Definition: vector.cpp:264
int allocsize
Definition: vector.hpp:45
Vector & operator+=(const Vector &v)
Definition: vector.cpp:190
double * data
Definition: vector.hpp:46
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:64
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:802
virtual N_Vector ToNVector()
Return a new wrapper SUNDIALS N_Vector of type SUNDIALS_NVEC_SERIAL.
Definition: vector.hpp:275
virtual ~Vector()
Destroys vector.
Definition: vector.hpp:378
void Neg()
(*this) = -(*this)
Definition: vector.cpp:256