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