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