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