MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
vector.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 // Implementation of data type vector
13 
14 #include "kernels.hpp"
15 #include "vector.hpp"
16 #include "../general/forall.hpp"
17 
18 #if defined(MFEM_USE_SUNDIALS)
19 #include "sundials.hpp"
20 #if defined(MFEM_USE_MPI)
21 #include <nvector/nvector_parallel.h>
22 #endif
23 #endif
24 
25 #ifdef MFEM_USE_OPENMP
26 #include <omp.h>
27 #endif
28 
29 #include <iostream>
30 #include <iomanip>
31 #include <cmath>
32 #include <cstdlib>
33 #include <ctime>
34 #include <limits>
35 
36 namespace mfem
37 {
38 
40 {
41  const int s = v.Size();
42  size = s;
43  if (s > 0)
44  {
45  MFEM_ASSERT(!v.data.Empty(), "invalid source vector");
46  data.New(s, v.data.GetMemoryType());
47  data.CopyFrom(v.data, s);
48  }
49  UseDevice(v.UseDevice());
50 }
51 
53 {
54  *this = std::move(v);
55 }
56 
57 void Vector::Load(std::istream **in, int np, int *dim)
58 {
59  int i, j, s;
60 
61  s = 0;
62  for (i = 0; i < np; i++)
63  {
64  s += dim[i];
65  }
66 
67  SetSize(s);
68 
69  int p = 0;
70  for (i = 0; i < np; i++)
71  {
72  for (j = 0; j < dim[i]; j++)
73  {
74  *in[i] >> data[p++];
75  // Clang's libc++ sets the failbit when (correctly) parsing subnormals,
76  // so we reset the failbit here.
77  if (!*in[i] && errno == ERANGE)
78  {
79  in[i]->clear();
80  }
81  }
82  }
83 }
84 
85 void Vector::Load(std::istream &in, int Size)
86 {
87  SetSize(Size);
88 
89  for (int i = 0; i < size; i++)
90  {
91  in >> data[i];
92  // Clang's libc++ sets the failbit when (correctly) parsing subnormals,
93  // so we reset the failbit here.
94  if (!in && errno == ERANGE)
95  {
96  in.clear();
97  }
98  }
99 }
100 
101 double &Vector::Elem(int i)
102 {
103  return operator()(i);
104 }
105 
106 const double &Vector::Elem(int i) const
107 {
108  return operator()(i);
109 }
110 
111 double Vector::operator*(const double *v) const
112 {
113  double dot = 0.0;
114 #ifdef MFEM_USE_LEGACY_OPENMP
115  #pragma omp parallel for reduction(+:dot)
116 #endif
117  for (int i = 0; i < size; i++)
118  {
119  dot += data[i] * v[i];
120  }
121  return dot;
122 }
123 
124 Vector &Vector::operator=(const double *v)
125 {
126  data.CopyFromHost(v, size);
127  return *this;
128 }
129 
131 {
132 #if 0
133  SetSize(v.Size(), v.data.GetMemoryType());
134  data.CopyFrom(v.data, v.Size());
135  UseDevice(v.UseDevice());
136 #else
137  SetSize(v.Size());
138  bool vuse = v.UseDevice();
139  const bool use_dev = UseDevice() || vuse;
140  v.UseDevice(use_dev);
141  // keep 'data' where it is, unless 'use_dev' is true
142  if (use_dev) { Write(); }
143  data.CopyFrom(v.data, v.Size());
144  v.UseDevice(vuse);
145 #endif
146  return *this;
147 }
148 
150 {
151  data = std::move(v.data);
152  size = v.size;
153  v.data.Reset();
154  v.size = 0;
155  return *this;
156 }
157 
158 Vector &Vector::operator=(double value)
159 {
160  const bool use_dev = UseDevice();
161  const int N = size;
162  auto y = Write(use_dev);
163  MFEM_FORALL_SWITCH(use_dev, i, N, y[i] = value;);
164  return *this;
165 }
166 
168 {
169  const bool use_dev = UseDevice();
170  const int N = size;
171  auto y = ReadWrite(use_dev);
172  MFEM_FORALL_SWITCH(use_dev, i, N, y[i] *= c;);
173  return *this;
174 }
175 
177 {
178  MFEM_ASSERT(size == v.size, "incompatible Vectors!");
179 
180  const bool use_dev = UseDevice() || v.UseDevice();
181  const int N = size;
182  auto y = ReadWrite(use_dev);
183  auto x = v.Read(use_dev);
184  MFEM_FORALL_SWITCH(use_dev, i, N, y[i] *= x[i];);
185  return *this;
186 }
187 
189 {
190  const bool use_dev = UseDevice();
191  const int N = size;
192  const double m = 1.0/c;
193  auto y = ReadWrite(use_dev);
194  MFEM_FORALL_SWITCH(use_dev, i, N, y[i] *= m;);
195  return *this;
196 }
197 
199 {
200  MFEM_ASSERT(size == v.size, "incompatible Vectors!");
201 
202  const bool use_dev = UseDevice() || v.UseDevice();
203  const int N = size;
204  auto y = ReadWrite(use_dev);
205  auto x = v.Read(use_dev);
206  MFEM_FORALL_SWITCH(use_dev, i, N, y[i] /= x[i];);
207  return *this;
208 }
209 
211 {
212  const bool use_dev = UseDevice();
213  const int N = size;
214  auto y = ReadWrite(use_dev);
215  MFEM_FORALL_SWITCH(use_dev, i, N, y[i] -= c;);
216  return *this;
217 }
218 
220 {
221  MFEM_ASSERT(size == v.size, "incompatible Vectors!");
222 
223  const bool use_dev = UseDevice() || v.UseDevice();
224  const int N = size;
225  auto y = ReadWrite(use_dev);
226  auto x = v.Read(use_dev);
227  MFEM_FORALL_SWITCH(use_dev, i, N, y[i] -= x[i];);
228  return *this;
229 }
230 
232 {
233  const bool use_dev = UseDevice();
234  const int N = size;
235  auto y = ReadWrite(use_dev);
236  MFEM_FORALL_SWITCH(use_dev, i, N, y[i] += c;);
237  return *this;
238 }
239 
241 {
242  MFEM_ASSERT(size == v.size, "incompatible Vectors!");
243 
244  const bool use_dev = UseDevice() || v.UseDevice();
245  const int N = size;
246  auto y = ReadWrite(use_dev);
247  auto x = v.Read(use_dev);
248  MFEM_FORALL_SWITCH(use_dev, i, N, y[i] += x[i];);
249  return *this;
250 }
251 
252 Vector &Vector::Add(const double a, const Vector &Va)
253 {
254  MFEM_ASSERT(size == Va.size, "incompatible Vectors!");
255 
256  if (a != 0.0)
257  {
258  const int N = size;
259  const bool use_dev = UseDevice() || Va.UseDevice();
260  auto y = ReadWrite(use_dev);
261  auto x = Va.Read(use_dev);
262  MFEM_FORALL_SWITCH(use_dev, i, N, y[i] += a * x[i];);
263  }
264  return *this;
265 }
266 
267 Vector &Vector::Set(const double a, const Vector &Va)
268 {
269  MFEM_ASSERT(size == Va.size, "incompatible Vectors!");
270 
271  const bool use_dev = UseDevice() || Va.UseDevice();
272  const int N = size;
273  auto x = Va.Read(use_dev);
274  auto y = Write(use_dev);
275  MFEM_FORALL_SWITCH(use_dev, i, N, y[i] = a * x[i];);
276  return *this;
277 }
278 
279 void Vector::SetVector(const Vector &v, int offset)
280 {
281  MFEM_ASSERT(v.Size() + offset <= size, "invalid sub-vector");
282 
283  const int vs = v.Size();
284  const double *vp = v.data;
285  double *p = data + offset;
286  for (int i = 0; i < vs; i++)
287  {
288  p[i] = vp[i];
289  }
290 }
291 
293 {
294  const bool use_dev = UseDevice();
295  const int N = size;
296  auto y = ReadWrite(use_dev);
297  MFEM_FORALL_SWITCH(use_dev, i, N, y[i] = -y[i];);
298 }
299 
300 void add(const Vector &v1, const Vector &v2, Vector &v)
301 {
302  MFEM_ASSERT(v.size == v1.size && v.size == v2.size,
303  "incompatible Vectors!");
304 
305 #if !defined(MFEM_USE_LEGACY_OPENMP)
306  const bool use_dev = v1.UseDevice() || v2.UseDevice() || v.UseDevice();
307  const int N = v.size;
308  // Note: get read access first, in case v is the same as v1/v2.
309  auto x1 = v1.Read(use_dev);
310  auto x2 = v2.Read(use_dev);
311  auto y = v.Write(use_dev);
312  MFEM_FORALL_SWITCH(use_dev, i, N, y[i] = x1[i] + x2[i];);
313 #else
314  #pragma omp parallel for
315  for (int i = 0; i < v.size; i++)
316  {
317  v.data[i] = v1.data[i] + v2.data[i];
318  }
319 #endif
320 }
321 
322 void add(const Vector &v1, double alpha, const Vector &v2, Vector &v)
323 {
324  MFEM_ASSERT(v.size == v1.size && v.size == v2.size,
325  "incompatible Vectors!");
326 
327  if (alpha == 0.0)
328  {
329  v = v1;
330  }
331  else if (alpha == 1.0)
332  {
333  add(v1, v2, v);
334  }
335  else
336  {
337 #if !defined(MFEM_USE_LEGACY_OPENMP)
338  const bool use_dev = v1.UseDevice() || v2.UseDevice() || v.UseDevice();
339  const int N = v.size;
340  // Note: get read access first, in case v is the same as v1/v2.
341  auto d_x = v1.Read(use_dev);
342  auto d_y = v2.Read(use_dev);
343  auto d_z = v.Write(use_dev);
344  MFEM_FORALL_SWITCH(use_dev, i, N, d_z[i] = d_x[i] + alpha * d_y[i];);
345 #else
346  const double *v1p = v1.data, *v2p = v2.data;
347  double *vp = v.data;
348  const int s = v.size;
349  #pragma omp parallel for
350  for (int i = 0; i < s; i++)
351  {
352  vp[i] = v1p[i] + alpha*v2p[i];
353  }
354 #endif
355  }
356 }
357 
358 void add(const double a, const Vector &x, const Vector &y, Vector &z)
359 {
360  MFEM_ASSERT(x.size == y.size && x.size == z.size,
361  "incompatible Vectors!");
362 
363  if (a == 0.0)
364  {
365  z = 0.0;
366  }
367  else if (a == 1.0)
368  {
369  add(x, y, z);
370  }
371  else
372  {
373 #if !defined(MFEM_USE_LEGACY_OPENMP)
374  const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
375  const int N = x.size;
376  // Note: get read access first, in case z is the same as x/y.
377  auto xd = x.Read(use_dev);
378  auto yd = y.Read(use_dev);
379  auto zd = z.Write(use_dev);
380  MFEM_FORALL_SWITCH(use_dev, i, N, zd[i] = a * (xd[i] + yd[i]););
381 #else
382  const double *xp = x.data;
383  const double *yp = y.data;
384  double *zp = z.data;
385  const int s = x.size;
386  #pragma omp parallel for
387  for (int i = 0; i < s; i++)
388  {
389  zp[i] = a * (xp[i] + yp[i]);
390  }
391 #endif
392  }
393 }
394 
395 void add(const double a, const Vector &x,
396  const double b, const Vector &y, Vector &z)
397 {
398  MFEM_ASSERT(x.size == y.size && x.size == z.size,
399  "incompatible Vectors!");
400 
401  if (a == 0.0)
402  {
403  z.Set(b, y);
404  }
405  else if (b == 0.0)
406  {
407  z.Set(a, x);
408  }
409 #if 0
410  else if (a == 1.0)
411  {
412  add(x, b, y, z);
413  }
414  else if (b == 1.0)
415  {
416  add(y, a, x, z);
417  }
418  else if (a == b)
419  {
420  add(a, x, y, z);
421  }
422 #endif
423  else
424  {
425 #if !defined(MFEM_USE_LEGACY_OPENMP)
426  const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
427  const int N = x.size;
428  // Note: get read access first, in case z is the same as x/y.
429  auto xd = x.Read(use_dev);
430  auto yd = y.Read(use_dev);
431  auto zd = z.Write(use_dev);
432  MFEM_FORALL_SWITCH(use_dev, i, N, zd[i] = a * xd[i] + b * yd[i];);
433 #else
434  const double *xp = x.data;
435  const double *yp = y.data;
436  double *zp = z.data;
437  const int s = x.size;
438  #pragma omp parallel for
439  for (int i = 0; i < s; i++)
440  {
441  zp[i] = a * xp[i] + b * yp[i];
442  }
443 #endif
444  }
445 }
446 
447 void subtract(const Vector &x, const Vector &y, Vector &z)
448 {
449  MFEM_ASSERT(x.size == y.size && x.size == z.size,
450  "incompatible Vectors!");
451 
452 #if !defined(MFEM_USE_LEGACY_OPENMP)
453  const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
454  const int N = x.size;
455  // Note: get read access first, in case z is the same as x/y.
456  auto xd = x.Read(use_dev);
457  auto yd = y.Read(use_dev);
458  auto zd = z.Write(use_dev);
459  MFEM_FORALL_SWITCH(use_dev, i, N, zd[i] = xd[i] - yd[i];);
460 #else
461  const double *xp = x.data;
462  const double *yp = y.data;
463  double *zp = z.data;
464  const int s = x.size;
465  #pragma omp parallel for
466  for (int i = 0; i < s; i++)
467  {
468  zp[i] = xp[i] - yp[i];
469  }
470 #endif
471 }
472 
473 void subtract(const double a, const Vector &x, const Vector &y, Vector &z)
474 {
475  MFEM_ASSERT(x.size == y.size && x.size == z.size,
476  "incompatible Vectors!");
477 
478  if (a == 0.)
479  {
480  z = 0.;
481  }
482  else if (a == 1.)
483  {
484  subtract(x, y, z);
485  }
486  else
487  {
488 #if !defined(MFEM_USE_LEGACY_OPENMP)
489  const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
490  const int N = x.size;
491  // Note: get read access first, in case z is the same as x/y.
492  auto xd = x.Read(use_dev);
493  auto yd = y.Read(use_dev);
494  auto zd = z.Write(use_dev);
495  MFEM_FORALL_SWITCH(use_dev, i, N, zd[i] = a * (xd[i] - yd[i]););
496 #else
497  const double *xp = x.data;
498  const double *yp = y.data;
499  double *zp = z.data;
500  const int s = x.size;
501  #pragma omp parallel for
502  for (int i = 0; i < s; i++)
503  {
504  zp[i] = a * (xp[i] - yp[i]);
505  }
506 #endif
507  }
508 }
509 
510 void Vector::median(const Vector &lo, const Vector &hi)
511 {
512  MFEM_ASSERT(size == lo.size && size == hi.size,
513  "incompatible Vectors!");
514 
515  const bool use_dev = UseDevice() || lo.UseDevice() || hi.UseDevice();
516  const int N = size;
517  // Note: get read access first, in case *this is the same as lo/hi.
518  auto l = lo.Read(use_dev);
519  auto h = hi.Read(use_dev);
520  auto m = Write(use_dev);
521  MFEM_FORALL_SWITCH(use_dev, i, N,
522  {
523  if (m[i] < l[i])
524  {
525  m[i] = l[i];
526  }
527  else if (m[i] > h[i])
528  {
529  m[i] = h[i];
530  }
531  });
532 }
533 
534 void Vector::GetSubVector(const Array<int> &dofs, Vector &elemvect) const
535 {
536  const int n = dofs.Size();
537  elemvect.SetSize(n);
538  const bool use_dev = dofs.UseDevice() || elemvect.UseDevice();
539  auto d_y = elemvect.Write(use_dev);
540  auto d_X = Read(use_dev);
541  auto d_dofs = dofs.Read(use_dev);
542  MFEM_FORALL_SWITCH(use_dev, i, n,
543  {
544  const int dof_i = d_dofs[i];
545  d_y[i] = dof_i >= 0 ? d_X[dof_i] : -d_X[-dof_i-1];
546  });
547 }
548 
549 void Vector::GetSubVector(const Array<int> &dofs, double *elem_data) const
550 {
552  const int n = dofs.Size();
553  for (int i = 0; i < n; i++)
554  {
555  const int j = dofs[i];
556  elem_data[i] = (j >= 0) ? data[j] : -data[-1-j];
557  }
558 }
559 
560 void Vector::SetSubVector(const Array<int> &dofs, const double value)
561 {
562  const bool use_dev = dofs.UseDevice();
563  const int n = dofs.Size();
564  // Use read+write access for *this - we only modify some of its entries
565  auto d_X = ReadWrite(use_dev);
566  auto d_dofs = dofs.Read(use_dev);
567  MFEM_FORALL_SWITCH(use_dev, i, n,
568  {
569  const int j = d_dofs[i];
570  if (j >= 0)
571  {
572  d_X[j] = value;
573  }
574  else
575  {
576  d_X[-1-j] = -value;
577  }
578  });
579 }
580 
581 void Vector::SetSubVector(const Array<int> &dofs, const Vector &elemvect)
582 {
583  MFEM_ASSERT(dofs.Size() == elemvect.Size(),
584  "Size mismatch: length of dofs is " << dofs.Size()
585  << ", length of elemvect is " << elemvect.Size());
586 
587  const bool use_dev = dofs.UseDevice() || elemvect.UseDevice();
588  const int n = dofs.Size();
589  // Use read+write access for X - we only modify some of its entries
590  auto d_X = ReadWrite(use_dev);
591  auto d_y = elemvect.Read(use_dev);
592  auto d_dofs = dofs.Read(use_dev);
593  MFEM_FORALL_SWITCH(use_dev, i, n,
594  {
595  const int dof_i = d_dofs[i];
596  if (dof_i >= 0)
597  {
598  d_X[dof_i] = d_y[i];
599  }
600  else
601  {
602  d_X[-1-dof_i] = -d_y[i];
603  }
604  });
605 }
606 
607 void Vector::SetSubVector(const Array<int> &dofs, double *elem_data)
608 {
609  // Use read+write access because we overwrite only part of the data.
611  const int n = dofs.Size();
612  for (int i = 0; i < n; i++)
613  {
614  const int j= dofs[i];
615  if (j >= 0)
616  {
617  operator()(j) = elem_data[i];
618  }
619  else
620  {
621  operator()(-1-j) = -elem_data[i];
622  }
623  }
624 }
625 
626 void Vector::AddElementVector(const Array<int> &dofs, const Vector &elemvect)
627 {
628  MFEM_ASSERT(dofs.Size() == elemvect.Size(), "Size mismatch: "
629  "length of dofs is " << dofs.Size() <<
630  ", length of elemvect is " << elemvect.Size());
631 
632  const bool use_dev = dofs.UseDevice() || elemvect.UseDevice();
633  const int n = dofs.Size();
634  auto d_y = elemvect.Read(use_dev);
635  auto d_X = ReadWrite(use_dev);
636  auto d_dofs = dofs.Read(use_dev);
637  MFEM_FORALL_SWITCH(use_dev, i, n,
638  {
639  const int j = d_dofs[i];
640  if (j >= 0)
641  {
642  d_X[j] += d_y[i];
643  }
644  else
645  {
646  d_X[-1-j] -= d_y[i];
647  }
648  });
649 }
650 
651 void Vector::AddElementVector(const Array<int> &dofs, double *elem_data)
652 {
654  const int n = dofs.Size();
655  for (int i = 0; i < n; i++)
656  {
657  const int j = dofs[i];
658  if (j >= 0)
659  {
660  operator()(j) += elem_data[i];
661  }
662  else
663  {
664  operator()(-1-j) -= elem_data[i];
665  }
666  }
667 }
668 
669 void Vector::AddElementVector(const Array<int> &dofs, const double a,
670  const Vector &elemvect)
671 {
672  MFEM_ASSERT(dofs.Size() == elemvect.Size(), "Size mismatch: "
673  "length of dofs is " << dofs.Size() <<
674  ", length of elemvect is " << elemvect.Size());
675 
676  const bool use_dev = dofs.UseDevice() || elemvect.UseDevice();
677  const int n = dofs.Size();
678  auto d_y = ReadWrite(use_dev);
679  auto d_x = elemvect.Read(use_dev);
680  auto d_dofs = dofs.Read(use_dev);
681  MFEM_FORALL_SWITCH(use_dev, i, n,
682  {
683  const int j = d_dofs[i];
684  if (j >= 0)
685  {
686  d_y[j] += a * d_x[i];
687  }
688  else
689  {
690  d_y[-1-j] -= a * d_x[i];
691  }
692  });
693 }
694 
695 void Vector::SetSubVectorComplement(const Array<int> &dofs, const double val)
696 {
697  const bool use_dev = UseDevice() || dofs.UseDevice();
698  const int n = dofs.Size();
699  const int N = size;
700  Vector dofs_vals(n, use_dev ?
703  auto d_data = ReadWrite(use_dev);
704  auto d_dofs_vals = dofs_vals.Write(use_dev);
705  auto d_dofs = dofs.Read(use_dev);
706  MFEM_FORALL_SWITCH(use_dev, i, n, d_dofs_vals[i] = d_data[d_dofs[i]];);
707  MFEM_FORALL_SWITCH(use_dev, i, N, d_data[i] = val;);
708  MFEM_FORALL_SWITCH(use_dev, i, n, d_data[d_dofs[i]] = d_dofs_vals[i];);
709 }
710 
711 void Vector::Print(std::ostream &os, int width) const
712 {
713  if (!size) { return; }
715  for (int i = 0; 1; )
716  {
717  os << ZeroSubnormal(data[i]);
718  i++;
719  if (i == size)
720  {
721  break;
722  }
723  if ( i % width == 0 )
724  {
725  os << '\n';
726  }
727  else
728  {
729  os << ' ';
730  }
731  }
732  os << '\n';
733 }
734 
735 #ifdef MFEM_USE_ADIOS2
737  const std::string& variable_name) const
738 {
739  if (!size) { return; }
741  os.engine.Put(variable_name, &data[0] );
742 }
743 #endif
744 
745 void Vector::Print_HYPRE(std::ostream &os) const
746 {
747  int i;
748  std::ios::fmtflags old_fmt = os.flags();
749  os.setf(std::ios::scientific);
750  std::streamsize old_prec = os.precision(14);
751 
752  os << size << '\n'; // number of rows
753 
754  data.Read(MemoryClass::HOST, size);
755  for (i = 0; i < size; i++)
756  {
757  os << ZeroSubnormal(data[i]) << '\n';
758  }
759 
760  os.precision(old_prec);
761  os.flags(old_fmt);
762 }
763 
764 void Vector::PrintHash(std::ostream &os) const
765 {
766  os << "size: " << size << '\n';
767  HashFunction hf;
768  hf.AppendDoubles(HostRead(), size);
769  os << "hash: " << hf.GetHash() << '\n';
770 }
771 
772 void Vector::Randomize(int seed)
773 {
774  // static unsigned int seed = time(0);
775  const double max = (double)(RAND_MAX) + 1.;
776 
777  if (seed == 0)
778  {
779  seed = (int)time(0);
780  }
781 
782  // srand(seed++);
783  srand((unsigned)seed);
784 
785  HostWrite();
786  for (int i = 0; i < size; i++)
787  {
788  data[i] = std::abs(rand()/max);
789  }
790 }
791 
792 double Vector::Norml2() const
793 {
794  // Scale entries of Vector on the fly, using algorithms from
795  // std::hypot() and LAPACK's drm2. This scaling ensures that the
796  // argument of each call to std::pow is <= 1 to avoid overflow.
797  if (0 == size)
798  {
799  return 0.0;
800  } // end if 0 == size
801 
803  if (1 == size)
804  {
805  return std::abs(data[0]);
806  } // end if 1 == size
807  return kernels::Norml2(size, (const double*) data);
808 }
809 
810 double Vector::Normlinf() const
811 {
812  HostRead();
813  double max = 0.0;
814  for (int i = 0; i < size; i++)
815  {
816  max = std::max(std::abs(data[i]), max);
817  }
818  return max;
819 }
820 
821 double Vector::Norml1() const
822 {
823  HostRead();
824  double sum = 0.0;
825  for (int i = 0; i < size; i++)
826  {
827  sum += std::abs(data[i]);
828  }
829  return sum;
830 }
831 
832 double Vector::Normlp(double p) const
833 {
834  MFEM_ASSERT(p > 0.0, "Vector::Normlp");
835 
836  if (p == 1.0)
837  {
838  return Norml1();
839  }
840  if (p == 2.0)
841  {
842  return Norml2();
843  }
844  if (p < infinity())
845  {
846  // Scale entries of Vector on the fly, using algorithms from
847  // std::hypot() and LAPACK's drm2. This scaling ensures that the
848  // argument of each call to std::pow is <= 1 to avoid overflow.
849  if (0 == size)
850  {
851  return 0.0;
852  } // end if 0 == size
853 
854  if (1 == size)
855  {
856  return std::abs(data[0]);
857  } // end if 1 == size
858 
859  double scale = 0.0;
860  double sum = 0.0;
861 
862  for (int i = 0; i < size; i++)
863  {
864  if (data[i] != 0.0)
865  {
866  const double absdata = std::abs(data[i]);
867  if (scale <= absdata)
868  {
869  sum = 1.0 + sum * std::pow(scale / absdata, p);
870  scale = absdata;
871  continue;
872  } // end if scale <= absdata
873  sum += std::pow(absdata / scale, p); // else scale > absdata
874  } // end if data[i] != 0
875  }
876  return scale * std::pow(sum, 1.0/p);
877  } // end if p < infinity()
878 
879  return Normlinf(); // else p >= infinity()
880 }
881 
882 double Vector::Max() const
883 {
884  if (size == 0) { return -infinity(); }
885 
886  HostRead();
887  double max = data[0];
888 
889  for (int i = 1; i < size; i++)
890  {
891  if (data[i] > max)
892  {
893  max = data[i];
894  }
895  }
896 
897  return max;
898 }
899 
900 double Vector::Sum() const
901 {
902  double sum = 0.0;
903 
904  const double *h_data = this->HostRead();
905  for (int i = 0; i < size; i++)
906  {
907  sum += h_data[i];
908  }
909 
910  return sum;
911 }
912 
913 #ifdef MFEM_USE_CUDA
914 static __global__ void cuKernelMin(const int N, double *gdsr, const double *x)
915 {
916  __shared__ double s_min[MFEM_CUDA_BLOCKS];
917  const int n = blockDim.x*blockIdx.x + threadIdx.x;
918  if (n>=N) { return; }
919  const int bid = blockIdx.x;
920  const int tid = threadIdx.x;
921  const int bbd = bid*blockDim.x;
922  const int rid = bbd+tid;
923  s_min[tid] = x[n];
924  for (int workers=blockDim.x>>1; workers>0; workers>>=1)
925  {
926  __syncthreads();
927  if (tid >= workers) { continue; }
928  if (rid >= N) { continue; }
929  const int dualTid = tid + workers;
930  if (dualTid >= N) { continue; }
931  const int rdd = bbd+dualTid;
932  if (rdd >= N) { continue; }
933  if (dualTid >= blockDim.x) { continue; }
934  s_min[tid] = fmin(s_min[tid], s_min[dualTid]);
935  }
936  if (tid==0) { gdsr[bid] = s_min[0]; }
937 }
938 
939 static Array<double> cuda_reduce_buf;
940 
941 static double cuVectorMin(const int N, const double *X)
942 {
943  const int tpb = MFEM_CUDA_BLOCKS;
944  const int blockSize = MFEM_CUDA_BLOCKS;
945  const int gridSize = (N+blockSize-1)/blockSize;
946  const int min_sz = (N%tpb)==0? (N/tpb) : (1+N/tpb);
947  cuda_reduce_buf.SetSize(min_sz);
948  Memory<double> &buf = cuda_reduce_buf.GetMemory();
949  double *d_min = buf.Write(MemoryClass::DEVICE, min_sz);
950  cuKernelMin<<<gridSize,blockSize>>>(N, d_min, X);
951  MFEM_GPU_CHECK(cudaGetLastError());
952  const double *h_min = buf.Read(MemoryClass::HOST, min_sz);
954  for (int i = 0; i < min_sz; i++) { min = fmin(min, h_min[i]); }
955  return min;
956 }
957 
958 static __global__ void cuKernelDot(const int N, double *gdsr,
959  const double *x, const double *y)
960 {
961  __shared__ double s_dot[MFEM_CUDA_BLOCKS];
962  const int n = blockDim.x*blockIdx.x + threadIdx.x;
963  if (n>=N) { return; }
964  const int bid = blockIdx.x;
965  const int tid = threadIdx.x;
966  const int bbd = bid*blockDim.x;
967  const int rid = bbd+tid;
968  s_dot[tid] = x[n] * y[n];
969  for (int workers=blockDim.x>>1; workers>0; workers>>=1)
970  {
971  __syncthreads();
972  if (tid >= workers) { continue; }
973  if (rid >= N) { continue; }
974  const int dualTid = tid + workers;
975  if (dualTid >= N) { continue; }
976  const int rdd = bbd+dualTid;
977  if (rdd >= N) { continue; }
978  if (dualTid >= blockDim.x) { continue; }
979  s_dot[tid] += s_dot[dualTid];
980  }
981  if (tid==0) { gdsr[bid] = s_dot[0]; }
982 }
983 
984 static double cuVectorDot(const int N, const double *X, const double *Y)
985 {
986  const int tpb = MFEM_CUDA_BLOCKS;
987  const int blockSize = MFEM_CUDA_BLOCKS;
988  const int gridSize = (N+blockSize-1)/blockSize;
989  const int dot_sz = (N%tpb)==0? (N/tpb) : (1+N/tpb);
990  cuda_reduce_buf.SetSize(dot_sz, Device::GetDeviceMemoryType());
991  Memory<double> &buf = cuda_reduce_buf.GetMemory();
992  double *d_dot = buf.Write(MemoryClass::DEVICE, dot_sz);
993  cuKernelDot<<<gridSize,blockSize>>>(N, d_dot, X, Y);
994  MFEM_GPU_CHECK(cudaGetLastError());
995  const double *h_dot = buf.Read(MemoryClass::HOST, dot_sz);
996  double dot = 0.0;
997  for (int i = 0; i < dot_sz; i++) { dot += h_dot[i]; }
998  return dot;
999 }
1000 #endif // MFEM_USE_CUDA
1001 
1002 #ifdef MFEM_USE_HIP
1003 static __global__ void hipKernelMin(const int N, double *gdsr, const double *x)
1004 {
1005  __shared__ double s_min[MFEM_HIP_BLOCKS];
1006  const int n = hipBlockDim_x*hipBlockIdx_x + hipThreadIdx_x;
1007  if (n>=N) { return; }
1008  const int bid = hipBlockIdx_x;
1009  const int tid = hipThreadIdx_x;
1010  const int bbd = bid*hipBlockDim_x;
1011  const int rid = bbd+tid;
1012  s_min[tid] = x[n];
1013  for (int workers=hipBlockDim_x>>1; workers>0; workers>>=1)
1014  {
1015  __syncthreads();
1016  if (tid >= workers) { continue; }
1017  if (rid >= N) { continue; }
1018  const int dualTid = tid + workers;
1019  if (dualTid >= N) { continue; }
1020  const int rdd = bbd+dualTid;
1021  if (rdd >= N) { continue; }
1022  if (dualTid >= hipBlockDim_x) { continue; }
1023  s_min[tid] = fmin(s_min[tid], s_min[dualTid]);
1024  }
1025  if (tid==0) { gdsr[bid] = s_min[0]; }
1026 }
1027 
1028 static Array<double> cuda_reduce_buf;
1029 
1030 static double hipVectorMin(const int N, const double *X)
1031 {
1032  const int tpb = MFEM_HIP_BLOCKS;
1033  const int blockSize = MFEM_HIP_BLOCKS;
1034  const int gridSize = (N+blockSize-1)/blockSize;
1035  const int min_sz = (N%tpb)==0 ? (N/tpb) : (1+N/tpb);
1036  cuda_reduce_buf.SetSize(min_sz);
1037  Memory<double> &buf = cuda_reduce_buf.GetMemory();
1038  double *d_min = buf.Write(MemoryClass::DEVICE, min_sz);
1039  hipLaunchKernelGGL(hipKernelMin,gridSize,blockSize,0,0,N,d_min,X);
1040  MFEM_GPU_CHECK(hipGetLastError());
1041  const double *h_min = buf.Read(MemoryClass::HOST, min_sz);
1043  for (int i = 0; i < min_sz; i++) { min = fmin(min, h_min[i]); }
1044  return min;
1045 }
1046 
1047 static __global__ void hipKernelDot(const int N, double *gdsr,
1048  const double *x, const double *y)
1049 {
1050  __shared__ double s_dot[MFEM_HIP_BLOCKS];
1051  const int n = hipBlockDim_x*hipBlockIdx_x + hipThreadIdx_x;
1052  if (n>=N) { return; }
1053  const int bid = hipBlockIdx_x;
1054  const int tid = hipThreadIdx_x;
1055  const int bbd = bid*hipBlockDim_x;
1056  const int rid = bbd+tid;
1057  s_dot[tid] = x[n] * y[n];
1058  for (int workers=hipBlockDim_x>>1; workers>0; workers>>=1)
1059  {
1060  __syncthreads();
1061  if (tid >= workers) { continue; }
1062  if (rid >= N) { continue; }
1063  const int dualTid = tid + workers;
1064  if (dualTid >= N) { continue; }
1065  const int rdd = bbd+dualTid;
1066  if (rdd >= N) { continue; }
1067  if (dualTid >= hipBlockDim_x) { continue; }
1068  s_dot[tid] += s_dot[dualTid];
1069  }
1070  if (tid==0) { gdsr[bid] = s_dot[0]; }
1071 }
1072 
1073 static double hipVectorDot(const int N, const double *X, const double *Y)
1074 {
1075  const int tpb = MFEM_HIP_BLOCKS;
1076  const int blockSize = MFEM_HIP_BLOCKS;
1077  const int gridSize = (N+blockSize-1)/blockSize;
1078  const int dot_sz = (N%tpb)==0 ? (N/tpb) : (1+N/tpb);
1079  cuda_reduce_buf.SetSize(dot_sz);
1080  Memory<double> &buf = cuda_reduce_buf.GetMemory();
1081  double *d_dot = buf.Write(MemoryClass::DEVICE, dot_sz);
1082  hipLaunchKernelGGL(hipKernelDot,gridSize,blockSize,0,0,N,d_dot,X,Y);
1083  MFEM_GPU_CHECK(hipGetLastError());
1084  const double *h_dot = buf.Read(MemoryClass::HOST, dot_sz);
1085  double dot = 0.0;
1086  for (int i = 0; i < dot_sz; i++) { dot += h_dot[i]; }
1087  return dot;
1088 }
1089 #endif // MFEM_USE_HIP
1090 
1091 double Vector::operator*(const Vector &v) const
1092 {
1093  MFEM_ASSERT(size == v.size, "incompatible Vectors!");
1094  if (size == 0) { return 0.0; }
1095 
1096  const bool use_dev = UseDevice() || v.UseDevice();
1097 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP) || defined(MFEM_USE_OPENMP)
1098  auto m_data = Read(use_dev);
1099 #else
1100  Read(use_dev);
1101 #endif
1102  auto v_data = v.Read(use_dev);
1103 
1104  if (!use_dev) { goto vector_dot_cpu; }
1105 
1106 #ifdef MFEM_USE_OCCA
1107  if (DeviceCanUseOcca())
1108  {
1109  return occa::linalg::dot<double,double,double>(
1111  }
1112 #endif
1113 
1114 #ifdef MFEM_USE_CUDA
1116  {
1117  return cuVectorDot(size, m_data, v_data);
1118  }
1119 #endif
1120 
1121 #ifdef MFEM_USE_HIP
1123  {
1124  return hipVectorDot(size, m_data, v_data);
1125  }
1126 #endif
1127 
1128 #ifdef MFEM_USE_OPENMP
1130  {
1131 #define MFEM_USE_OPENMP_DETERMINISTIC_DOT
1132 #ifdef MFEM_USE_OPENMP_DETERMINISTIC_DOT
1133  // By default, use a deterministic way of computing the dot product
1134  static Vector th_dot;
1135  #pragma omp parallel
1136  {
1137  const int nt = omp_get_num_threads();
1138  #pragma omp master
1139  th_dot.SetSize(nt);
1140  const int tid = omp_get_thread_num();
1141  const int stride = (size + nt - 1)/nt;
1142  const int start = tid*stride;
1143  const int stop = std::min(start + stride, size);
1144  double my_dot = 0.0;
1145  for (int i = start; i < stop; i++)
1146  {
1147  my_dot += m_data[i] * v_data[i];
1148  }
1149  #pragma omp barrier
1150  th_dot(tid) = my_dot;
1151  }
1152  return th_dot.Sum();
1153 #else
1154  // The standard way of computing the dot product is non-deterministic
1155  double prod = 0.0;
1156  #pragma omp parallel for reduction(+:prod)
1157  for (int i = 0; i < size; i++)
1158  {
1159  prod += m_data[i] * v_data[i];
1160  }
1161  return prod;
1162 #endif // MFEM_USE_OPENMP_DETERMINISTIC_DOT
1163  }
1164 #endif // MFEM_USE_OPENMP
1166  {
1167  const int N = size;
1168  auto v_data_ = v.Read();
1169  auto m_data_ = Read();
1170  Vector dot(1);
1171  dot.UseDevice(true);
1172  auto d_dot = dot.Write();
1173  dot = 0.0;
1174  MFEM_FORALL(i, N, d_dot[0] += m_data_[i] * v_data_[i];);
1175  dot.HostReadWrite();
1176  return dot[0];
1177  }
1178 vector_dot_cpu:
1179  return operator*(v_data);
1180 }
1181 
1182 double Vector::Min() const
1183 {
1184  if (size == 0) { return infinity(); }
1185 
1186  const bool use_dev = UseDevice();
1187  auto m_data = Read(use_dev);
1188 
1189  if (!use_dev) { goto vector_min_cpu; }
1190 
1191 #ifdef MFEM_USE_OCCA
1192  if (DeviceCanUseOcca())
1193  {
1194  return occa::linalg::min<double,double>(OccaMemoryRead(data, size));
1195  }
1196 #endif
1197 
1198 #ifdef MFEM_USE_CUDA
1200  {
1201  return cuVectorMin(size, m_data);
1202  }
1203 #endif
1204 
1205 #ifdef MFEM_USE_HIP
1207  {
1208  return hipVectorMin(size, m_data);
1209  }
1210 #endif
1211 
1212 #ifdef MFEM_USE_OPENMP
1214  {
1215  double minimum = m_data[0];
1216  #pragma omp parallel for reduction(min:minimum)
1217  for (int i = 0; i < size; i++)
1218  {
1219  minimum = std::min(minimum, m_data[i]);
1220  }
1221  return minimum;
1222  }
1223 #endif
1224 
1226  {
1227  const int N = size;
1228  auto m_data_ = Read();
1229  Vector min(1);
1230  min = infinity();
1231  min.UseDevice(true);
1232  auto d_min = min.ReadWrite();
1233  MFEM_FORALL(i, N, d_min[0] = (d_min[0]<m_data_[i])?d_min[0]:m_data_[i];);
1234  min.HostReadWrite();
1235  return min[0];
1236  }
1237 
1238 vector_min_cpu:
1239  double minimum = data[0];
1240  for (int i = 1; i < size; i++)
1241  {
1242  if (m_data[i] < minimum)
1243  {
1244  minimum = m_data[i];
1245  }
1246  }
1247  return minimum;
1248 }
1249 
1250 
1251 #ifdef MFEM_USE_SUNDIALS
1252 
1253 Vector::Vector(N_Vector nv)
1254 {
1255  N_Vector_ID nvid = N_VGetVectorID(nv);
1256 
1257  switch (nvid)
1258  {
1259  case SUNDIALS_NVEC_SERIAL:
1260  SetDataAndSize(NV_DATA_S(nv), NV_LENGTH_S(nv));
1261  break;
1262 #ifdef MFEM_USE_MPI
1263  case SUNDIALS_NVEC_PARALLEL:
1264  SetDataAndSize(NV_DATA_P(nv), NV_LOCLENGTH_P(nv));
1265  break;
1266 #endif
1267  default:
1268  MFEM_ABORT("N_Vector type " << nvid << " is not supported");
1269  }
1270 }
1271 
1272 void Vector::ToNVector(N_Vector &nv, long global_length)
1273 {
1274  MFEM_ASSERT(nv, "N_Vector handle is NULL");
1275  N_Vector_ID nvid = N_VGetVectorID(nv);
1276 
1277  switch (nvid)
1278  {
1279  case SUNDIALS_NVEC_SERIAL:
1280  MFEM_ASSERT(NV_OWN_DATA_S(nv) == SUNFALSE, "invalid serial N_Vector");
1281  NV_DATA_S(nv) = data;
1282  NV_LENGTH_S(nv) = size;
1283  break;
1284 #ifdef MFEM_USE_MPI
1285  case SUNDIALS_NVEC_PARALLEL:
1286  MFEM_ASSERT(NV_OWN_DATA_P(nv) == SUNFALSE, "invalid parallel N_Vector");
1287  NV_DATA_P(nv) = data;
1288  NV_LOCLENGTH_P(nv) = size;
1289  if (global_length == 0)
1290  {
1291  global_length = NV_GLOBLENGTH_P(nv);
1292 
1293  if (global_length == 0 && global_length != size)
1294  {
1295  MPI_Comm sundials_comm = NV_COMM_P(nv);
1296  long local_size = size;
1297  MPI_Allreduce(&local_size, &global_length, 1, MPI_LONG,
1298  MPI_SUM,sundials_comm);
1299  }
1300  }
1301  NV_GLOBLENGTH_P(nv) = global_length;
1302  break;
1303 #endif
1304  default:
1305  MFEM_ABORT("N_Vector type " << nvid << " is not supported");
1306  }
1307 }
1308 
1309 #endif // MFEM_USE_SUNDIALS
1310 
1311 }
Hash function for data sequences.
Definition: hash.hpp:455
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:560
void SetVector(const Vector &v, int offset)
Definition: vector.cpp:279
Memory< T > & GetMemory()
Return a reference to the Memory object used by the Array.
Definition: array.hpp:117
Memory< double > data
Definition: vector.hpp:64
double & Elem(int i)
Access Vector entries. Index i = 0 .. size-1.
Definition: vector.cpp:101
Device memory; using CUDA or HIP *Malloc and *Free.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:792
double & operator()(int i)
Access Vector entries using () for 0-based indexing.
Definition: vector.hpp:607
Biwise-OR of all HIP backends.
Definition: device.hpp:90
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:534
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:450
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
T * Write(MemoryClass mc, int size)
Get write-only access to the memory with the given MemoryClass.
bool UseDevice() const
Return the device flag of the Memory object used by the Array.
Definition: array.hpp:123
std::string GetHash() const
Return the hash string for the current sequence and reset (clear) the sequence.
Definition: hash.cpp:60
double Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:810
void CopyFrom(const Memory &src, int size)
Copy size entries from src to *this.
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:772
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:117
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:300
HashFunction & AppendDoubles(const double *doubles, size_t num_doubles)
Add a sequence of doubles for hashing, given as a c-array.
Definition: hash.hpp:508
double operator*(const double *) const
Dot product with a double * array.
Definition: vector.cpp:111
MemoryType GetMemoryType() const
Return a MemoryType that is currently valid. If both the host and the device pointers are currently v...
T ZeroSubnormal(T val)
Definition: vector.hpp:485
Vector & operator=(const double *v)
Copy Size() entries from v.
Definition: vector.cpp:124
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Definition: vector.cpp:57
double Normlp(double p) const
Returns the l_p norm of the vector.
Definition: vector.cpp:832
void CopyFromHost(const T *src, int size)
Copy size entries from the host pointer src to *this.
bool DeviceCanUseOcca()
Function that determines if an OCCA kernel should be used, based on the current mfem::Device configur...
Definition: occa.hpp:69
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:446
double b
Definition: lissajous.cpp:42
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition: device.hpp:273
const occa::memory OccaMemoryRead(const Memory< T > &mem, size_t size)
Wrap a Memory object as occa::memory for read only access with the mfem::Device MemoryClass. The returned occa::memory is associated with the default occa::device used by MFEM.
Definition: occa.hpp:37
void median(const Vector &lo, const Vector &hi)
v = median(v,lo,hi) entrywise. Implementation assumes lo &lt;= hi.
Definition: vector.cpp:510
Biwise-OR of all OpenMP backends.
Definition: device.hpp:92
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:304
MFEM_HOST_DEVICE double Norml2(const int size, const T *data)
Returns the l2 norm of the Vector with given size and data.
Definition: kernels.hpp:133
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:626
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:442
Vector & operator/=(double c)
Definition: vector.cpp:188
Biwise-OR of all CUDA backends.
Definition: device.hpp:88
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the dofs Array to the given val.
Definition: vector.cpp:695
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
Vector & operator*=(double c)
Definition: vector.cpp:167
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1182
Vector & operator+=(double c)
Definition: vector.cpp:231
double Norml1() const
Returns the l_1 norm of the vector.
Definition: vector.cpp:821
static MemoryType GetHostMemoryType()
Get the current Host MemoryType. This is the MemoryType used by most MFEM classes when allocating mem...
Definition: device.hpp:264
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:711
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:447
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
void Print_HYPRE(std::ostream &out) const
Prints vector to stream out in HYPRE_Vector format.
Definition: vector.cpp:745
static bool Allows(unsigned long b_mask)
Return true if any of the backends in the backend mask, b_mask, are allowed.
Definition: device.hpp:258
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:156
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:267
double a
Definition: lissajous.cpp:41
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:454
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:252
Host memory; using new[] and delete[].
virtual MFEM_DEPRECATED N_Vector ToNVector()
(DEPRECATED) Return a new wrapper SUNDIALS N_Vector of type SUNDIALS_NVEC_SERIAL. ...
Definition: vector.hpp:467
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:882
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
bool Empty() const
Return true if the Memory object is empty, see Reset().
T * ReadWrite(MemoryClass mc, int size)
Get read-write access to the memory with the given MemoryClass.
int dim
Definition: ex24.cpp:53
Vector & operator-=(double c)
Definition: vector.cpp:210
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:46
const double alpha
Definition: ex15.cpp:369
Vector data type.
Definition: vector.hpp:60
virtual bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
Definition: vector.hpp:120
void PrintHash(std::ostream &out) const
Print the Vector size and hash of its data.
Definition: vector.cpp:764
RefCoord s[3]
const T * Read(MemoryClass mc, int size) const
Get read-only access to the memory with the given MemoryClass.
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:900
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:438
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:458
[device] Debug backend: host memory is READ/WRITE protected while a device is in use. It allows to test the &quot;device&quot; code-path (using separate host/device memory pools and host &lt;-&gt; device transfers) without any GPU hardware. As &#39;DEBUG&#39; is sometimes used as a macro, _DEVICE has been added to avoid conflicts.
Definition: device.hpp:75
void Neg()
(*this) = -(*this)
Definition: vector.cpp:292