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
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
134  inline operator double *() { return data; }
135
136  /// Conversion to const double *.
137  /** @note This conversion function makes it possible to use [] for indexing
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
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
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 {
307  for (int i = 0; i < n; i++)
308  {
309  if (!IsFinite(v[i])) { bad++; }
310  }
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 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
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:666
double & operator()(int i)
Access Vector entries using () for 0-based indexing.
Definition: vector.hpp:357
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Definition: vector.cpp:458
void StealData(double **p)
Changes the ownership of the data; after the call the Vector is empty.
Definition: vector.hpp:144
int Capacity() const
Return the size of the currently allocated data array.
Definition: vector.hpp:124
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&lt;&gt; for class Vector.
Definition: vector.hpp:381
double Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:703
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:647
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:129
double operator*(const double *) const
Dot product with a double * array.
Definition: vector.cpp:89
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
bool OwnsData() const
Definition: vector.hpp:141
double Normlp(double p) const
Returns the l_p norm of the vector.
Definition: vector.cpp:723
int dim
Definition: ex3.cpp:47
void Swap(Vector &other)
Swap the contents of two Vectors.
Definition: vector.hpp:373
void median(const Vector &lo, const Vector &hi)
v = median(v,lo,hi) entrywise. Implementation assumes lo &lt;= hi.
Definition: vector.cpp:441
double DistanceSquared(const double *x, const double *y, const int n)
Definition: vector.hpp:394
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
void SetData(double *d)
Definition: vector.hpp:94
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
double DistanceTo(const double *p) const
Compute the Euclidean distance to another vector.
Definition: vector.hpp:416
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
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
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:785
double Norml1() const
Returns the l_1 norm of the vector.
Definition: vector.cpp:713
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
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:605
void Print_HYPRE(std::ostream &out) const
Prints vector to stream out in HYPRE_Vector format.
Definition: vector.cpp:629
double DistanceSquaredTo(const double *p) const
Compute the square of the Euclidean distance to another vector.
Definition: vector.hpp:411
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:101
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
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:772
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&lt;double&gt;::infinity()
Definition: vector.hpp:42
int CheckFinite() const
Definition: vector.hpp:272
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
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:798
virtual N_Vector ToNVector()
Return a new wrapper SUNDIALS N_Vector of type SUNDIALS_NVEC_SERIAL.
Definition: vector.hpp:283
virtual ~Vector()
Destroys vector.
Definition: vector.hpp:386
void Neg()
(*this) = -(*this)
Definition: vector.cpp:252