MFEM v4.8.0
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-2025, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#ifndef MFEM_VECTOR
13#define MFEM_VECTOR
14
15#include "../general/array.hpp"
16#ifdef MFEM_USE_ADIOS2
18#endif
21#include "../general/device.hpp"
22
23#include <cmath>
24#include <cstdlib>
25#include <iostream>
26#include <limits>
27#include <type_traits>
28#include <initializer_list>
29#if defined(_MSC_VER) && (_MSC_VER < 1800)
30#include <float.h>
31#define isfinite _finite
32#endif
33
34#ifdef MFEM_USE_MPI
35#include <mpi.h>
36#endif
37
38namespace mfem
39{
40
41/** Count the number of entries in an array of doubles for which isfinite
42 is false, i.e. the entry is a NaN or +/-Inf. */
43inline int CheckFinite(const real_t *v, const int n);
44
45/// Define a shortcut for std::numeric_limits<double>::infinity()
46#ifndef __CYGWIN__
48{
49 return std::numeric_limits<real_t>::infinity();
50}
51#else
52// On Cygwin math.h defines a function 'infinity()' which will conflict with the
53// above definition if we have 'using namespace mfem;' and try to use something
54// like 'double a = infinity();'. This 'infinity()' function is non-standard and
55// is defined by the Newlib C standard library implementation used by Cygwin,
56// see https://en.wikipedia.org/wiki/Newlib, http://www.sourceware.org/newlib.
57using ::infinity;
58#endif
59
60/// Generate a random `real_t` number in the interval [0,1) using rand().
62{
63 constexpr real_t max = (real_t)(RAND_MAX) + 1_r;
64#if defined(MFEM_USE_SINGLE)
65 // Note: For RAND_MAX = 2^31-1, float(RAND_MAX) = 2^31, and
66 // max = float(RAND_MAX)+1.0f is equal to 2^31 too; this is the actual value
67 // that we want for max. However, this rounding behavior means that for a
68 // range of values of rand() close to RAND_MAX, the expression
69 // float(rand())/max will give 1.0f.
70 //
71 // Therefore, to ensure we return a number less than 1, we take the minimum
72 // of float(rand())/max and (1 - 2^(-24)), where the latter number is the
73 // largest float less than 1.
74 return std::fmin(real_t(rand())/max, 0.9999999403953552_r);
75#else
76 return real_t(rand())/max;
77#endif
78}
79
80/// Vector data type.
81class Vector
82{
83protected:
84
86 int size;
87
88public:
89
90 /** Default constructor for Vector. Sets size = 0, and calls Memory::Reset on
91 data through Memory<double>'s default constructor. */
92 Vector(): size(0) { }
93
94 /// Copy constructor. Allocates a new data array and copies the data.
95 Vector(const Vector &);
96
97 /// Move constructor. "Steals" data from its argument.
98 Vector(Vector&& v);
99
100 /// @brief Creates vector of size s.
101 /// @warning Entries are not initialized to zero!
102 explicit Vector(int s);
103
104 /// Creates a vector referencing an array of doubles, owned by someone else.
105 /** The pointer @a data_ can be NULL. The data array can be replaced later
106 with SetData(). */
107 Vector(real_t *data_, int size_)
108 { data.Wrap(data_, size_, false); size = size_; }
109
110 /** @brief Create a Vector referencing a sub-vector of the Vector @a base
111 starting at the given offset, @a base_offset, and size @a size_. */
112 Vector(Vector &base, int base_offset, int size_)
113 : data(base.data, base_offset, size_), size(size_) { }
114
115 /// Create a Vector of size @a size_ using MemoryType @a mt.
116 Vector(int size_, MemoryType mt)
117 : data(size_, mt), size(size_) { }
118
119 /** @brief Create a Vector of size @a size_ using host MemoryType @a h_mt and
120 device MemoryType @a d_mt. */
121 Vector(int size_, MemoryType h_mt, MemoryType d_mt)
122 : data(size_, h_mt, d_mt), size(size_) { }
123
124 /// Create a vector from a statically sized C-style array of convertible type
125 template <typename CT, int N>
126 explicit Vector(const CT (&values)[N]) : Vector(N)
127 { std::copy(values, values + N, begin()); }
128
129 /// Create a vector using a braced initializer list
130 template <typename CT, typename std::enable_if<
131 std::is_convertible<CT,real_t>::value,bool>::type = true>
132 explicit Vector(std::initializer_list<CT> values) : Vector(values.size())
133 { std::copy(values.begin(), values.end(), begin()); }
134
135 /// Enable execution of Vector operations using the mfem::Device.
136 /** The default is to use Backend::CPU (serial execution on each MPI rank),
137 regardless of the mfem::Device configuration.
138
139 When appropriate, MFEM functions and class methods will enable the use
140 of the mfem::Device for their Vector parameters.
141
142 Some derived classes, e.g. GridFunction, enable the use of the
143 mfem::Device by default. */
144 virtual void UseDevice(bool use_dev) const { data.UseDevice(use_dev); }
145
146 /// Return the device flag of the Memory object used by the Vector
147 virtual bool UseDevice() const { return data.UseDevice(); }
148
149 /// Reads a vector from multiple files
150 void Load(std::istream ** in, int np, int * dim);
151
152 /// Load a vector from an input stream.
153 void Load(std::istream &in, int Size);
154
155 /// Load a vector from an input stream, reading the size from the stream.
156 void Load(std::istream &in) { int s; in >> s; Load(in, s); }
157
158 /// @brief Resize the vector to size @a s.
159 /** If the new size is less than or equal to Capacity() then the internal
160 data array remains the same. Otherwise, the old array is deleted, if
161 owned, and a new array of size @a s is allocated without copying the
162 previous content of the Vector.
163 @warning In the second case above (new size greater than current one),
164 the vector will allocate new data array, even if it did not own the
165 original data! Also, new entries are not initialized! */
166 void SetSize(int s);
167
168 /// Resize the vector to size @a s using MemoryType @a mt.
169 void SetSize(int s, MemoryType mt);
170
171 /// Resize the vector to size @a s using the MemoryType of @a v.
172 void SetSize(int s, const Vector &v) { SetSize(s, v.GetMemory().GetMemoryType()); }
173
174 /// Set the Vector data.
175 /// @warning This method should be called only when OwnsData() is false.
176 void SetData(real_t *d) { data.Wrap(d, data.Capacity(), false); }
177
178 /// Set the Vector data and size.
179 /** The Vector does not assume ownership of the new data. The new size is
180 also used as the new Capacity().
181 @warning This method should be called only when OwnsData() is false.
182 @sa NewDataAndSize(). */
183 void SetDataAndSize(real_t *d, int s) { data.Wrap(d, s, false); size = s; }
184
185 /// Set the Vector data and size, deleting the old data, if owned.
186 /** The Vector does not assume ownership of the new data. The new size is
187 also used as the new Capacity().
188 @sa SetDataAndSize(). */
189 void NewDataAndSize(real_t *d, int s)
190 {
191 data.Delete();
192 SetDataAndSize(d, s);
193 }
194
195 /// Reset the Vector to use the given external Memory @a mem and size @a s.
196 /** If @a own_mem is false, the Vector will not own any of the pointers of
197 @a mem.
198
199 Note that when @a own_mem is true, the @a mem object can be destroyed
200 immediately by the caller but `mem.Delete()` should NOT be called since
201 the Vector object takes ownership of all pointers owned by @a mem.
202
203 @sa NewDataAndSize(). */
204 inline void NewMemoryAndSize(const Memory<real_t> &mem, int s, bool own_mem);
205
206 /// Reset the Vector to be a reference to a sub-vector of @a base.
207 inline void MakeRef(Vector &base, int offset, int size);
208
209 /** @brief Reset the Vector to be a reference to a sub-vector of @a base
210 without changing its current size. */
211 inline void MakeRef(Vector &base, int offset);
212
213 /// Set the Vector data (host pointer) ownership flag.
214 void MakeDataOwner() const { data.SetHostPtrOwner(true); }
215
216 /// Destroy a vector
217 void Destroy();
218
219 /** @brief Delete the device pointer, if owned. If @a copy_to_host is true
220 and the data is valid only on device, move it to host before deleting.
221 Invalidates the device memory. */
222 void DeleteDevice(bool copy_to_host = true)
223 { data.DeleteDevice(copy_to_host); }
224
225 /// Returns the size of the vector.
226 inline int Size() const { return size; }
227
228 /// Return the size of the currently allocated data array.
229 /** It is always true that Capacity() >= Size(). */
230 inline int Capacity() const { return data.Capacity(); }
231
232 /// Return a pointer to the beginning of the Vector data.
233 /** @warning This method should be used with caution as it gives write access
234 to the data of const-qualified Vector%s. */
235 inline real_t *GetData() const
236 { return const_cast<real_t*>((const real_t*)data); }
237
238 /// Conversion to `double *`. Deprecated.
239 MFEM_DEPRECATED inline operator real_t *() { return data; }
240
241 /// Conversion to `const double *`. Deprecated.
242 MFEM_DEPRECATED inline operator const real_t *() const { return data; }
243
244 /// STL-like begin.
245 inline real_t *begin() { return data; }
246
247 /// STL-like end.
248 inline real_t *end() { return data + size; }
249
250 /// STL-like begin (const version).
251 inline const real_t *begin() const { return data; }
252
253 /// STL-like end (const version).
254 inline const real_t *end() const { return data + size; }
255
256 /// Return a reference to the Memory object used by the Vector.
258
259 /** @brief Return a reference to the Memory object used by the Vector, const
260 version. */
261 const Memory<real_t> &GetMemory() const { return data; }
262
263 /// Update the memory location of the vector to match @a v.
264 void SyncMemory(const Vector &v) const { GetMemory().Sync(v.GetMemory()); }
265
266 /// Update the alias memory location of the vector to match @a v.
267 void SyncAliasMemory(const Vector &v) const
268 { GetMemory().SyncAlias(v.GetMemory(),Size()); }
269
270 /// Read the Vector data (host pointer) ownership flag.
271 inline bool OwnsData() const { return data.OwnsHostPtr(); }
272
273 /// Changes the ownership of the data; after the call the Vector is empty
274 inline void StealData(real_t **p)
275 { *p = data; data.Reset(); size = 0; }
276
277 /// Changes the ownership of the data; after the call the Vector is empty
278 inline real_t *StealData() { real_t *p; StealData(&p); return p; }
279
280 /// Access Vector entries. Index i = 0 .. size-1.
281 real_t &Elem(int i);
282
283 /// Read only access to Vector entries. Index i = 0 .. size-1.
284 const real_t &Elem(int i) const;
285
286 /// Access Vector entries using () for 0-based indexing.
287 /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
288 inline real_t &operator()(int i);
289
290 /// Read only access to Vector entries using () for 0-based indexing.
291 /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
292 inline const real_t &operator()(int i) const;
293
294 /// Access Vector entries using [] for 0-based indexing.
295 /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
296 inline real_t &operator[](int i) { return (*this)(i); }
297
298 /// Read only access to Vector entries using [] for 0-based indexing.
299 /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
300 inline const real_t &operator[](int i) const { return (*this)(i); }
301
302 /// Dot product with a `double *` array.
303 /// This function always executes on the CPU. A HostRead() will be called if
304 /// required.
305 /// To optionally execute on the device:
306 /// Vector tmp(v, Size());
307 /// res = (*this) * tmp;
308 real_t operator*(const real_t *v) const;
309
310 /// Return the inner-product.
311 real_t operator*(const Vector &v) const;
312
313 /// Copy Size() entries from @a v.
314 Vector &operator=(const real_t *v);
315
316 /// Copy assignment.
317 /** @note Defining this method overwrites the implicitly defined copy
318 assignment operator. */
319 Vector &operator=(const Vector &v);
320
321 /// Move assignment
322 Vector &operator=(Vector&& v);
323
324 /// Redefine '=' for vector = constant.
325 Vector &operator=(real_t value);
326
328
329 /// Component-wise scaling: (*this)(i) *= v(i)
330 Vector &operator*=(const Vector &v);
331
333
334 /// Component-wise division: (*this)(i) /= v(i)
335 Vector &operator/=(const Vector &v);
336
338
339 Vector &operator-=(const Vector &v);
340
342
343 Vector &operator+=(const Vector &v);
344
345 /// (*this) += a * Va
346 Vector &Add(const real_t a, const Vector &Va);
347
348 /// (*this) = a * x
349 Vector &Set(const real_t a, const Vector &x);
350
351 void SetVector(const Vector &v, int offset);
352
353 void AddSubVector(const Vector &v, int offset);
354
355 /// (*this) = -(*this)
356 void Neg();
357
358 /// (*this)(i) = 1.0 / (*this)(i)
359 void Reciprocal();
360
361 /// Swap the contents of two Vectors
362 inline void Swap(Vector &other);
363
364 /// Set v = v1 + v2.
365 friend void add(const Vector &v1, const Vector &v2, Vector &v);
366
367 /// Set v = v1 + alpha * v2.
368 friend void add(const Vector &v1, real_t alpha, const Vector &v2, Vector &v);
369
370 /// z = a * (x + y)
371 friend void add(const real_t a, const Vector &x, const Vector &y, Vector &z);
372
373 /// z = a * x + b * y
374 friend void add(const real_t a, const Vector &x,
375 const real_t b, const Vector &y, Vector &z);
376
377 /// Set v = v1 - v2.
378 friend void subtract(const Vector &v1, const Vector &v2, Vector &v);
379
380 /// z = a * (x - y)
381 friend void subtract(const real_t a, const Vector &x,
382 const Vector &y, Vector &z);
383
384 /// Computes cross product of this vector with another 3D vector.
385 /// vout = this x vin.
386 void cross3D(const Vector &vin, Vector &vout) const;
387
388 /// v = median(v,lo,hi) entrywise. Implementation assumes lo <= hi.
389 void median(const Vector &lo, const Vector &hi);
390
391 /// Extract entries listed in @a dofs to the output Vector @a elemvect.
392 /** Negative dof values cause the -dof-1 position in @a elemvect to receive
393 the -val in from this Vector. */
394 void GetSubVector(const Array<int> &dofs, Vector &elemvect) const;
395
396 /// Extract entries listed in @a dofs to the output array @a elem_data.
397 /** Negative dof values cause the -dof-1 position in @a elem_data to receive
398 the -val in from this Vector. */
399 void GetSubVector(const Array<int> &dofs, real_t *elem_data) const;
400
401 /// Set the entries listed in @a dofs to the given @a value.
402 /** Negative dof values cause the -dof-1 position in this Vector to receive
403 the -value. */
404 void SetSubVector(const Array<int> &dofs, const real_t value);
405
406 /** @brief Set the entries listed in @a dofs to the values given in the @a
407 elemvect Vector. Negative dof values cause the -dof-1 position in this
408 Vector to receive the -val from @a elemvect. */
409 void SetSubVector(const Array<int> &dofs, const Vector &elemvect);
410
411 /** @brief Set the entries listed in @a dofs to the values given the @a ,
412 elem_data array. Negative dof values cause the -dof-1 position in this
413 Vector to receive the -val from @a elem_data. */
414 void SetSubVector(const Array<int> &dofs, real_t *elem_data);
415
416 /** @brief Add elements of the @a elemvect Vector to the entries listed in @a
417 dofs. Negative dof values cause the -dof-1 position in this Vector to add
418 the -val from @a elemvect. */
419 void AddElementVector(const Array<int> & dofs, const Vector & elemvect);
420
421 /** @brief Add elements of the @a elem_data array to the entries listed in @a
422 dofs. Negative dof values cause the -dof-1 position in this Vector to add
423 the -val from @a elem_data. */
424 void AddElementVector(const Array<int> & dofs, real_t *elem_data);
425
426 /** @brief Add @a times the elements of the @a elemvect Vector to the entries
427 listed in @a dofs. Negative dof values cause the -dof-1 position in this
428 Vector to add the -a*val from @a elemvect. */
429 void AddElementVector(const Array<int> & dofs, const real_t a,
430 const Vector & elemvect);
431
432 /// Set all vector entries NOT in the @a dofs Array to the given @a val.
433 void SetSubVectorComplement(const Array<int> &dofs, const real_t val);
434
435 /// Prints vector to stream out.
436 void Print(std::ostream &out = mfem::out, int width = 8) const;
437
438#ifdef MFEM_USE_ADIOS2
439 /// Prints vector to stream out.
440 /// @param out adios2stream output
441 /// @param variable_name variable name associated with current Vector
442 void Print(adios2stream & out, const std::string& variable_name) const;
443#endif
444
445 /// Prints vector to stream out in HYPRE_Vector format.
446 void Print_HYPRE(std::ostream &out) const;
447
448 /// Prints vector as a List for importing into Mathematica.
449 /** The resulting file can be read into Mathematica using an expression such
450 as: myVec = Get["output_file_name"]
451 The Mathematica variable "myVec" will then be assigned to a new
452 List object containing the data from this MFEM Vector. */
453 void PrintMathematica(std::ostream &out = mfem::out) const;
454
455 /// Print the Vector size and hash of its data.
456 /** This is a compact text representation of the Vector contents that can be
457 used to compare vectors from different runs without the need to save the
458 whole vector. */
459 void PrintHash(std::ostream &out) const;
460
461 /// Set random values in the vector.
462 void Randomize(int seed = 0);
463 /// Returns the l2 norm of the vector.
464 real_t Norml2() const;
465 /// Returns the l_infinity norm of the vector.
466 real_t Normlinf() const;
467 /// Returns the l_1 norm of the vector.
468 real_t Norml1() const;
469 /// Returns the l_p norm of the vector.
470 real_t Normlp(real_t p) const;
471 /// Returns the maximal element of the vector.
472 real_t Max() const;
473 /// Returns the minimal element of the vector.
474 real_t Min() const;
475 /// Return the sum of the vector entries
476 real_t Sum() const;
477 /// Compute the square of the Euclidean distance to another vector.
478 inline real_t DistanceSquaredTo(const real_t *p) const;
479 /// Compute the square of the Euclidean distance to another vector.
480 inline real_t DistanceSquaredTo(const Vector &p) const;
481 /// Compute the Euclidean distance to another vector.
482 inline real_t DistanceTo(const real_t *p) const;
483 /// Compute the Euclidean distance to another vector.
484 inline real_t DistanceTo(const Vector &p) const;
485
486 /** @brief Count the number of entries in the Vector for which isfinite
487 is false, i.e. the entry is a NaN or +/-Inf. */
488 int CheckFinite() const { return mfem::CheckFinite(HostRead(), size); }
489
490 /// Destroys vector.
491 virtual ~Vector();
492
493 /// Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
494 virtual const real_t *Read(bool on_dev = true) const
495 { return mfem::Read(data, size, on_dev); }
496
497 /// Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
498 virtual const real_t *HostRead() const
499 { return mfem::Read(data, size, false); }
500
501 /// Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
502 virtual real_t *Write(bool on_dev = true)
503 { return mfem::Write(data, size, on_dev); }
504
505 /// Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
506 virtual real_t *HostWrite()
507 { return mfem::Write(data, size, false); }
508
509 /// Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
510 virtual real_t *ReadWrite(bool on_dev = true)
511 { return mfem::ReadWrite(data, size, on_dev); }
512
513 /// Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
515 { return mfem::ReadWrite(data, size, false); }
516
517};
518
519// Inline methods
520
521template <typename T>
522inline T ZeroSubnormal(T val)
523{
524 return (std::fpclassify(val) == FP_SUBNORMAL) ? 0.0 : val;
525}
526
527inline bool IsFinite(const real_t &val)
528{
529 // isfinite didn't appear in a standard until C99, and later C++11. It wasn't
530 // standard in C89 or C++98. PGI as of 14.7 still defines it as a macro.
531#ifdef isfinite
532 return isfinite(val);
533#else
534 return std::isfinite(val);
535#endif
536}
537
538inline int CheckFinite(const real_t *v, const int n)
539{
540 int bad = 0;
541 for (int i = 0; i < n; i++)
542 {
543 if (!IsFinite(v[i])) { bad++; }
544 }
545 return bad;
546}
547
548inline Vector::Vector(int s)
549{
550 MFEM_ASSERT(s>=0,"Unexpected negative size.");
551 size = s;
552 if (s > 0)
553 {
554 data.New(s);
555 }
556}
557
558inline void Vector::SetSize(int s)
559{
560 if (s == size)
561 {
562 return;
563 }
564 if (s <= data.Capacity())
565 {
566 size = s;
567 return;
568 }
569 // preserve a valid MemoryType and device flag
570 const MemoryType mt = data.GetMemoryType();
571 const bool use_dev = data.UseDevice();
572 data.Delete();
573 size = s;
574 data.New(s, mt);
575 data.UseDevice(use_dev);
576}
577
578inline void Vector::SetSize(int s, MemoryType mt)
579{
580 if (mt == data.GetMemoryType())
581 {
582 if (s == size)
583 {
584 return;
585 }
586 if (s <= data.Capacity())
587 {
588 size = s;
589 return;
590 }
591 }
592 const bool use_dev = data.UseDevice();
593 data.Delete();
594 if (s > 0)
595 {
596 data.New(s, mt);
597 size = s;
598 }
599 else
600 {
601 data.Reset();
602 size = 0;
603 }
604 data.UseDevice(use_dev);
605}
606
607inline void Vector::NewMemoryAndSize(const Memory<real_t> &mem, int s,
608 bool own_mem)
609{
610 data.Delete();
611 size = s;
612 if (own_mem)
613 {
614 data = mem;
615 }
616 else
617 {
618 data.MakeAlias(mem, 0, s);
619 }
620}
621
622inline void Vector::MakeRef(Vector &base, int offset, int s)
623{
624 data.Delete();
625 size = s;
626 data.MakeAlias(base.GetMemory(), offset, s);
627}
628
629inline void Vector::MakeRef(Vector &base, int offset)
630{
631 data.Delete();
632 data.MakeAlias(base.GetMemory(), offset, size);
633}
634
635inline void Vector::Destroy()
636{
637 const bool use_dev = data.UseDevice();
638 data.Delete();
639 size = 0;
640 data.Reset();
641 data.UseDevice(use_dev);
642}
643
645{
646 MFEM_ASSERT(data && i >= 0 && i < size,
647 "index [" << i << "] is out of range [0," << size << ")");
648
649 return data[i];
650}
651
652inline const real_t &Vector::operator()(int i) const
653{
654 MFEM_ASSERT(data && i >= 0 && i < size,
655 "index [" << i << "] is out of range [0," << size << ")");
656
657 return data[i];
658}
659
660inline void Vector::Swap(Vector &other)
661{
662 mfem::Swap(data, other.data);
663 mfem::Swap(size, other.size);
664}
665
666/// Specialization of the template function Swap<> for class Vector
667template<> inline void Swap<Vector>(Vector &a, Vector &b)
668{
669 a.Swap(b);
670}
671
673{
674 data.Delete();
675}
676
677inline real_t DistanceSquared(const real_t *x, const real_t *y, const int n)
678{
679 real_t d = 0.0;
680
681 for (int i = 0; i < n; i++)
682 {
683 d += (x[i]-y[i])*(x[i]-y[i]);
684 }
685
686 return d;
687}
688
689inline real_t Distance(const real_t *x, const real_t *y, const int n)
690{
691 return std::sqrt(DistanceSquared(x, y, n));
692}
693
694inline real_t Distance(const Vector &x, const Vector &y)
695{
696 return x.DistanceTo(y);
697}
698
700{
701 return DistanceSquared(data, p, size);
702}
703
705{
706 MFEM_ASSERT(p.Size() == Size(), "Incompatible vector sizes.");
707 return DistanceSquared(data, p.data, size);
708}
709
710inline real_t Vector::DistanceTo(const real_t *p) const
711{
712 return Distance(data, p, size);
713}
714
715inline real_t Vector::DistanceTo(const Vector &p) const
716{
717 MFEM_ASSERT(p.Size() == Size(), "Incompatible vector sizes.");
718 return Distance(data, p.data, size);
719}
720
721/// Returns the inner product of x and y
722/** In parallel this computes the inner product of the local vectors,
723 producing different results on each MPI rank.
724*/
725inline real_t InnerProduct(const Vector &x, const Vector &y)
726{
727 return x * y;
728}
729
730#ifdef MFEM_USE_MPI
731/// Returns the inner product of x and y in parallel
732/** In parallel this computes the inner product of the global vectors,
733 producing identical results on each MPI rank.
734*/
735inline real_t InnerProduct(MPI_Comm comm, const Vector &x, const Vector &y)
736{
737 real_t loc_prod = x * y;
738 real_t glb_prod;
739 MPI_Allreduce(&loc_prod, &glb_prod, 1, MFEM_MPI_REAL_T, MPI_SUM, comm);
740 return glb_prod;
741}
742#endif
743
744} // namespace mfem
745
746#endif
Class used by MFEM to store pointers to host and/or device memory.
void SetHostPtrOwner(bool own) const
Set/clear the ownership flag for the host pointer. Ownership indicates whether the pointer will be de...
int Capacity() const
Return the size of the allocated memory.
void MakeAlias(const Memory &base, int offset, int size)
Create a memory object that points inside the memory object base.
bool UseDevice() const
Read the internal device flag.
void SyncAlias(const Memory &base, int alias_size) const
Update the alias Memory *this to match the memory location (all valid locations) of its base Memory,...
bool OwnsHostPtr() const
Return true if the host pointer is owned. Ownership indicates whether the pointer will be deleted by ...
MemoryType GetMemoryType() const
Return a MemoryType that is currently valid. If both the host and the device pointers are currently v...
void Sync(const Memory &other) const
Copy the host/device pointer validity flags from other to *this.
void DeleteDevice(bool copy_to_host=true)
Delete the device pointer, if owned. If copy_to_host is true and the data is valid only on device,...
void Reset()
Reset the memory to be empty, ensuring that Delete() will be a no-op.
void Wrap(T *ptr, int size, bool own)
Wrap an externally allocated host pointer, ptr with the current host memory type returned by MemoryMa...
void Delete()
Delete the owned pointers and reset the Memory object.
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
Vector data type.
Definition vector.hpp:82
void Randomize(int seed=0)
Set random values in the vector.
Definition vector.cpp:915
void PrintHash(std::ostream &out) const
Print the Vector size and hash of its data.
Definition vector.cpp:907
real_t & Elem(int i)
Access Vector entries. Index i = 0 .. size-1.
Definition vector.cpp:172
Vector(int size_, MemoryType mt)
Create a Vector of size size_ using MemoryType mt.
Definition vector.hpp:116
void StealData(real_t **p)
Changes the ownership of the data; after the call the Vector is empty.
Definition vector.hpp:274
void SetSize(int s, const Vector &v)
Resize the vector to size s using the MemoryType of v.
Definition vector.hpp:172
void SetVector(const Vector &v, int offset)
Definition vector.cpp:349
virtual const real_t * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition vector.hpp:498
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:494
Vector(const CT(&values)[N])
Create a vector from a statically sized C-style array of convertible type.
Definition vector.hpp:126
const real_t & operator[](int i) const
Read only access to Vector entries using [] for 0-based indexing.
Definition vector.hpp:300
void median(const Vector &lo, const Vector &hi)
v = median(v,lo,hi) entrywise. Implementation assumes lo <= hi.
Definition vector.cpp:629
void Load(std::istream &in)
Load a vector from an input stream, reading the size from the stream.
Definition vector.hpp:156
void Neg()
(*this) = -(*this)
Definition vector.cpp:375
real_t * begin()
STL-like begin.
Definition vector.hpp:245
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition vector.cpp:830
virtual real_t * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:510
void SetDataAndSize(real_t *d, int s)
Set the Vector data and size.
Definition vector.hpp:183
real_t operator*(const real_t *v) const
Definition vector.cpp:182
Vector(Vector &base, int base_offset, int size_)
Create a Vector referencing a sub-vector of the Vector base starting at the given offset,...
Definition vector.hpp:112
real_t Normlinf() const
Returns the l_infinity norm of the vector.
Definition vector.cpp:972
void MakeDataOwner() const
Set the Vector data (host pointer) ownership flag.
Definition vector.hpp:214
real_t Norml1() const
Returns the l_1 norm of the vector.
Definition vector.cpp:985
bool OwnsData() const
Read the Vector data (host pointer) ownership flag.
Definition vector.hpp:271
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:679
friend void subtract(const Vector &v1, const Vector &v2, Vector &v)
Set v = v1 - v2.
Definition vector.cpp:547
Memory< real_t > & GetMemory()
Return a reference to the Memory object used by the Vector.
Definition vector.hpp:257
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add elements of the elemvect Vector to the entries listed in dofs. Negative dof values cause the -dof...
Definition vector.cpp:745
Memory< real_t > data
Definition vector.hpp:85
real_t & operator[](int i)
Access Vector entries using [] for 0-based indexing.
Definition vector.hpp:296
int CheckFinite() const
Count the number of entries in the Vector for which isfinite is false, i.e. the entry is a NaN or +/-...
Definition vector.hpp:488
void AddSubVector(const Vector &v, int offset)
Definition vector.cpp:362
const real_t * end() const
STL-like end (const version).
Definition vector.hpp:254
void Swap(Vector &other)
Swap the contents of two Vectors.
Definition vector.hpp:660
Vector & operator*=(real_t c)
Definition vector.cpp:237
void SyncAliasMemory(const Vector &v) const
Update the alias memory location of the vector to match v.
Definition vector.hpp:267
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:931
void SyncMemory(const Vector &v) const
Update the memory location of the vector to match v.
Definition vector.hpp:264
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Definition vector.cpp:126
real_t DistanceSquaredTo(const real_t *p) const
Compute the square of the Euclidean distance to another vector.
Definition vector.hpp:699
Vector & Set(const real_t a, const Vector &x)
(*this) = a * x
Definition vector.cpp:337
real_t * end()
STL-like end.
Definition vector.hpp:248
const Memory< real_t > & GetMemory() const
Return a reference to the Memory object used by the Vector, const version.
Definition vector.hpp:261
void NewMemoryAndSize(const Memory< real_t > &mem, int s, bool own_mem)
Reset the Vector to use the given external Memory mem and size s.
Definition vector.hpp:607
real_t Max() const
Returns the maximal element of the vector.
Definition vector.cpp:1164
virtual bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
Definition vector.hpp:147
Vector(int size_, MemoryType h_mt, MemoryType d_mt)
Create a Vector of size size_ using host MemoryType h_mt and device MemoryType d_mt.
Definition vector.hpp:121
void Destroy()
Destroy a vector.
Definition vector.hpp:635
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition vector.hpp:144
real_t Sum() const
Return the sum of the vector entries.
Definition vector.cpp:1204
void PrintMathematica(std::ostream &out=mfem::out) const
Prints vector as a List for importing into Mathematica.
Definition vector.cpp:883
void Print_HYPRE(std::ostream &out) const
Prints vector to stream out in HYPRE_Vector format.
Definition vector.cpp:864
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
void Reciprocal()
(*this)(i) = 1.0 / (*this)(i)
Definition vector.cpp:383
virtual real_t * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition vector.hpp:506
void NewDataAndSize(real_t *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition vector.hpp:189
real_t Normlp(real_t p) const
Returns the l_p norm of the vector.
Definition vector.cpp:998
Vector & operator-=(real_t c)
Definition vector.cpp:280
const real_t * begin() const
STL-like begin (const version).
Definition vector.hpp:251
int Capacity() const
Return the size of the currently allocated data array.
Definition vector.hpp:230
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:235
Vector & operator=(const real_t *v)
Copy Size() entries from v.
Definition vector.cpp:196
void SetSubVectorComplement(const Array< int > &dofs, const real_t val)
Set all vector entries NOT in the dofs Array to the given val.
Definition vector.cpp:814
void SetData(real_t *d)
Definition vector.hpp:176
virtual real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition vector.hpp:514
virtual ~Vector()
Destroys vector.
Definition vector.hpp:672
real_t Min() const
Returns the minimal element of the vector.
Definition vector.cpp:1123
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition vector.cpp:653
Vector & operator+=(real_t c)
Definition vector.cpp:301
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:502
real_t * StealData()
Changes the ownership of the data; after the call the Vector is empty.
Definition vector.hpp:278
Vector & Add(const real_t a, const Vector &Va)
(*this) += a * Va
Definition vector.cpp:322
real_t & operator()(int i)
Access Vector entries using () for 0-based indexing.
Definition vector.hpp:644
void DeleteDevice(bool copy_to_host=true)
Delete the device pointer, if owned. If copy_to_host is true and the data is valid only on device,...
Definition vector.hpp:222
Vector(real_t *data_, int size_)
Creates a vector referencing an array of doubles, owned by someone else.
Definition vector.hpp:107
void cross3D(const Vector &vin, Vector &vout) const
Definition vector.cpp:616
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition vector.hpp:622
Vector(std::initializer_list< CT > values)
Create a vector using a braced initializer list.
Definition vector.hpp:132
Vector & operator/=(real_t c)
Definition vector.cpp:258
friend void add(const Vector &v1, const Vector &v2, Vector &v)
Set v = v1 + v2.
Definition vector.cpp:391
real_t DistanceTo(const real_t *p) const
Compute the Euclidean distance to another vector.
Definition vector.hpp:710
const real_t alpha
Definition ex15.cpp:369
int dim
Definition ex24.cpp:53
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const T * Read(const Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read access to mem with the mfem::Device's DeviceMemoryClass, if on_dev = true,...
Definition device.hpp:341
void Swap(Array< T > &, Array< T > &)
Definition array.hpp:672
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition vector.hpp:47
T * Write(Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for write access to mem with the mfem::Device's DeviceMemoryClass, if on_dev = true,...
Definition device.hpp:358
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:66
real_t rand_real()
Generate a random real_t number in the interval [0,1) using rand().
Definition vector.hpp:61
bool IsFinite(const real_t &val)
Definition vector.hpp:527
real_t Distance(const real_t *x, const real_t *y, const int n)
Definition vector.hpp:689
real_t InnerProduct(HypreParVector *x, HypreParVector *y)
Definition hypre.cpp:468
real_t DistanceSquared(const real_t *x, const real_t *y, const int n)
Definition vector.hpp:677
T * ReadWrite(Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read+write access to mem with the mfem::Device's DeviceMemoryClass,...
Definition device.hpp:375
int CheckFinite(const real_t *v, const int n)
Definition vector.hpp:538
T ZeroSubnormal(T val)
Definition vector.hpp:522
float real_t
Definition config.hpp:43
MemoryType
Memory types supported by MFEM.
void Swap< Vector >(Vector &a, Vector &b)
Specialization of the template function Swap<> for class Vector.
Definition vector.hpp:667
real_t p(const Vector &x, real_t t)