MFEM  v3.1
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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 // Implementation of data type vector
13 
14 #include <iostream>
15 #include <iomanip>
16 #include <cmath>
17 #include <cstdlib>
18 #include <ctime>
19 #include <limits>
20 
21 #include "vector.hpp"
22 
23 namespace mfem
24 {
25 
27 {
28  int s = v.Size();
29  if (s > 0)
30  {
31  allocsize = size = s;
32  data = new double[s];
33  for (int i = 0; i < s; i++)
34  {
35  data[i] = v(i);
36  }
37  }
38  else
39  {
40  allocsize = size = 0;
41  data = NULL;
42  }
43 }
44 
45 void Vector::Load(std::istream **in, int np, int *dim)
46 {
47  int i, j, s;
48 
49  s = 0;
50  for (i = 0; i < np; i++)
51  {
52  s += dim[i];
53  }
54 
55  SetSize(s);
56 
57  int p = 0;
58  for (i = 0; i < np; i++)
59  for (j = 0; j < dim[i]; j++)
60  {
61  *in[i] >> data[p++];
62  }
63 }
64 
65 void Vector::Load(std::istream &in, int Size)
66 {
67  SetSize(Size);
68 
69  for (int i = 0; i < size; i++)
70  {
71  in >> data[i];
72  }
73 }
74 
75 double &Vector::Elem(int i)
76 {
77  return operator()(i);
78 }
79 
80 const double &Vector::Elem(int i) const
81 {
82  return operator()(i);
83 }
84 
85 double Vector::operator*(const double *v) const
86 {
87  int s = size;
88  const double *d = data;
89  double prod = 0.0;
90 #ifdef MFEM_USE_OPENMP
91  #pragma omp parallel for reduction(+:prod)
92 #endif
93  for (int i = 0; i < s; i++)
94  {
95  prod += d[i] * v[i];
96  }
97  return prod;
98 }
99 
100 double Vector::operator*(const Vector &v) const
101 {
102 #ifdef MFEM_DEBUG
103  if (v.size != size)
104  {
105  mfem_error("Vector::operator*(const Vector &) const");
106  }
107 #endif
108 
109  return operator*(v.data);
110 }
111 
112 Vector &Vector::operator=(const double *v)
113 {
114  for (int i = 0; i < size; i++)
115  {
116  data[i] = v[i];
117  }
118  return *this;
119 }
120 
122 {
123  SetSize(v.Size());
124  for (int i = 0; i < size; i++)
125  {
126  data[i] = v.data[i];
127  }
128  return *this;
129 }
130 
131 Vector &Vector::operator=(double value)
132 {
133  int i, s = size;
134  double *p = data, v = value;
135  for (i = 0; i < s; i++)
136  {
137  *(p++) = v;
138  }
139  return *this;
140 }
141 
143 {
144  for (int i = 0; i < size; i++)
145  {
146  data[i] *= c;
147  }
148  return *this;
149 }
150 
152 {
153  double m = 1.0/c;
154  for (int i = 0; i < size; i++)
155  {
156  data[i] *= m;
157  }
158  return *this;
159 }
160 
162 {
163  for (int i = 0; i < size; i++)
164  {
165  data[i] -= c;
166  }
167  return *this;
168 }
169 
171 {
172 #ifdef MFEM_DEBUG
173  if (size != v.size)
174  {
175  mfem_error("Vector::operator-=(const Vector &)");
176  }
177 #endif
178  for (int i = 0; i < size; i++)
179  {
180  data[i] -= v(i);
181  }
182  return *this;
183 }
184 
186 {
187 #ifdef MFEM_DEBUG
188  if (size != v.size)
189  {
190  mfem_error("Vector::operator+=(const Vector &)");
191  }
192 #endif
193  for (int i = 0; i < size; i++)
194  {
195  data[i] += v(i);
196  }
197  return *this;
198 }
199 
200 Vector &Vector::Add(const double a, const Vector &Va)
201 {
202 #ifdef MFEM_DEBUG
203  if (size != Va.size)
204  {
205  mfem_error("Vector::Add(const double, const Vector &)");
206  }
207 #endif
208  if (a != 0.0)
209  {
210  for (int i = 0; i < size; i++)
211  {
212  data[i] += a * Va(i);
213  }
214  }
215  return *this;
216 }
217 
218 Vector &Vector::Set(const double a, const Vector &Va)
219 {
220 #ifdef MFEM_DEBUG
221  if (size != Va.size)
222  {
223  mfem_error("Vector::Set(const double, const Vector &)");
224  }
225 #endif
226  for (int i = 0; i < size; i++)
227  {
228  data[i] = a * Va(i);
229  }
230  return *this;
231 }
232 
233 void Vector::SetVector(const Vector &v, int offset)
234 {
235  int vs = v.Size();
236  double *vp = v.data, *p = data + offset;
237 
238 #ifdef MFEM_DEBUG
239  if (offset+vs > size)
240  {
241  mfem_error("Vector::SetVector(const Vector &, int)");
242  }
243 #endif
244 
245  for (int i = 0; i < vs; i++)
246  {
247  p[i] = vp[i];
248  }
249 }
250 
252 {
253  for (int i = 0; i < size; i++)
254  {
255  data[i] = -data[i];
256  }
257 }
258 
259 void add(const Vector &v1, const Vector &v2, Vector &v)
260 {
261 #ifdef MFEM_DEBUG
262  if (v.size != v1.size || v.size != v2.size)
263  {
264  mfem_error("add(Vector &v1, Vector &v2, Vector &v)");
265  }
266 #endif
267 
268 #ifdef MFEM_USE_OPENMP
269  #pragma omp parallel for
270 #endif
271  for (int i = 0; i < v.size; i++)
272  {
273  v.data[i] = v1.data[i] + v2.data[i];
274  }
275 }
276 
277 void add(const Vector &v1, double alpha, const Vector &v2, Vector &v)
278 {
279 #ifdef MFEM_DEBUG
280  if (v.size != v1.size || v.size != v2.size)
281  {
282  mfem_error ("add(Vector &v1, double alpha, Vector &v2, Vector &v)");
283  }
284 #endif
285  if (alpha == 0.0)
286  {
287  v = v1;
288  }
289  else if (alpha == 1.0)
290  {
291  add(v1, v2, v);
292  }
293  else
294  {
295  const double *v1p = v1.data, *v2p = v2.data;
296  double *vp = v.data;
297  int s = v.size;
298 #ifdef MFEM_USE_OPENMP
299  #pragma omp parallel for
300 #endif
301  for (int i = 0; i < s; i++)
302  {
303  vp[i] = v1p[i] + alpha*v2p[i];
304  }
305  }
306 }
307 
308 void add(const double a, const Vector &x, const Vector &y, Vector &z)
309 {
310 #ifdef MFEM_DEBUG
311  if (x.size != y.size || x.size != z.size)
312  mfem_error ("add(const double a, const Vector &x, const Vector &y,"
313  " Vector &z)");
314 #endif
315  if (a == 0.0)
316  {
317  z = 0.0;
318  }
319  else if (a == 1.0)
320  {
321  add(x, y, z);
322  }
323  else
324  {
325  const double *xp = x.data;
326  const double *yp = y.data;
327  double *zp = z.data;
328  int s = x.size;
329 
330 #ifdef MFEM_USE_OPENMP
331  #pragma omp parallel for
332 #endif
333  for (int i = 0; i < s; i++)
334  {
335  zp[i] = a * (xp[i] + yp[i]);
336  }
337  }
338 }
339 
340 void add(const double a, const Vector &x,
341  const double b, const Vector &y, Vector &z)
342 {
343 #ifdef MFEM_DEBUG
344  if (x.size != y.size || x.size != z.size)
345  mfem_error("add(const double a, const Vector &x,\n"
346  " const double b, const Vector &y, Vector &z)");
347 #endif
348  if (a == 0.0)
349  {
350  z.Set(b, y);
351  }
352  else if (b == 0.0)
353  {
354  z.Set(a, x);
355  }
356  else if (a == 1.0)
357  {
358  add(x, b, y, z);
359  }
360  else if (b == 1.0)
361  {
362  add(y, a, x, z);
363  }
364  else if (a == b)
365  {
366  add(a, x, y, z);
367  }
368  else
369  {
370  const double *xp = x.data;
371  const double *yp = y.data;
372  double *zp = z.data;
373  int s = x.size;
374 
375 #ifdef MFEM_USE_OPENMP
376  #pragma omp parallel for
377 #endif
378  for (int i = 0; i < s; i++)
379  {
380  zp[i] = a * xp[i] + b * yp[i];
381  }
382  }
383 }
384 
385 void subtract(const Vector &x, const Vector &y, Vector &z)
386 {
387 #ifdef MFEM_DEBUG
388  if (x.size != y.size || x.size != z.size)
389  {
390  mfem_error ("subtract(const Vector &, const Vector &, Vector &)");
391  }
392 #endif
393  const double *xp = x.data;
394  const double *yp = y.data;
395  double *zp = z.data;
396  int s = x.size;
397 
398 #ifdef MFEM_USE_OPENMP
399  #pragma omp parallel for
400 #endif
401  for (int i = 0; i < s; i++)
402  {
403  zp[i] = xp[i] - yp[i];
404  }
405 }
406 
407 void subtract(const double a, const Vector &x, const Vector &y, Vector &z)
408 {
409 #ifdef MFEM_DEBUG
410  if (x.size != y.size || x.size != z.size)
411  mfem_error("subtract(const double a, const Vector &x,"
412  " const Vector &y, Vector &z)");
413 #endif
414 
415  if (a == 0.)
416  {
417  z = 0.;
418  }
419  else if (a == 1.)
420  {
421  subtract(x, y, z);
422  }
423  else
424  {
425  const double *xp = x.data;
426  const double *yp = y.data;
427  double *zp = z.data;
428  int s = x.size;
429 
430 #ifdef MFEM_USE_OPENMP
431  #pragma omp parallel for
432 #endif
433  for (int i = 0; i < s; i++)
434  {
435  zp[i] = a * (xp[i] - yp[i]);
436  }
437  }
438 }
439 
440 void Vector::median(const Vector &lo, const Vector &hi)
441 {
442  double *v = data;
443 
444  for (int i = 0; i < size; i++)
445  {
446  if (v[i] < lo[i])
447  {
448  v[i] = lo[i];
449  }
450  else if (v[i] > hi[i])
451  {
452  v[i] = hi[i];
453  }
454  }
455 }
456 
457 void Vector::GetSubVector(const Array<int> &dofs, Vector &elemvect) const
458 {
459  int i, j, n = dofs.Size();
460 
461  elemvect.SetSize (n);
462 
463  for (i = 0; i < n; i++)
464  {
465  if ((j=dofs[i]) >= 0)
466  {
467  elemvect(i) = data[j];
468  }
469  else
470  {
471  elemvect(i) = -data[-1-j];
472  }
473  }
474 }
475 
476 void Vector::GetSubVector(const Array<int> &dofs, double *elem_data) const
477 {
478  int i, j, n = dofs.Size();
479 
480  for (i = 0; i < n; i++)
481  {
482  if ((j=dofs[i]) >= 0)
483  {
484  elem_data[i] = data[j];
485  }
486  else
487  {
488  elem_data[i] = -data[-1-j];
489  }
490  }
491 }
492 
493 void Vector::SetSubVector(const Array<int> &dofs, const Vector &elemvect)
494 {
495  int i, j, n = dofs.Size();
496 
497  for (i = 0; i < n; i++)
498  {
499  if ((j=dofs[i]) >= 0)
500  {
501  data[j] = elemvect(i);
502  }
503  else
504  {
505  data[-1-j] = -elemvect(i);
506  }
507  }
508 }
509 
510 void Vector::SetSubVector(const Array<int> &dofs, double *elem_data)
511 {
512  int i, j, n = dofs.Size();
513 
514  for (i = 0; i < n; i++)
515  {
516  if ((j=dofs[i]) >= 0)
517  {
518  data[j] = elem_data[i];
519  }
520  else
521  {
522  data[-1-j] = -elem_data[i];
523  }
524  }
525 }
526 
527 void Vector::AddElementVector(const Array<int> &dofs, const Vector &elemvect)
528 {
529  int i, j, n = dofs.Size();
530 
531  for (i = 0; i < n; i++)
532  if ((j=dofs[i]) >= 0)
533  {
534  data[j] += elemvect(i);
535  }
536  else
537  {
538  data[-1-j] -= elemvect(i);
539  }
540 }
541 
542 void Vector::AddElementVector(const Array<int> &dofs, double *elem_data)
543 {
544  int i, j, n = dofs.Size();
545 
546  for (i = 0; i < n; i++)
547  {
548  if ((j = dofs[i]) >= 0)
549  {
550  data[j] += elem_data[i];
551  }
552  else
553  {
554  data[-1-j] -= elem_data[i];
555  }
556  }
557 }
558 
559 void Vector::AddElementVector(const Array<int> &dofs, const double a,
560  const Vector &elemvect)
561 {
562  int i, j, n = dofs.Size();
563 
564  for (i = 0; i < n; i++)
565  if ((j=dofs[i]) >= 0)
566  {
567  data[j] += a * elemvect(i);
568  }
569  else
570  {
571  data[-1-j] -= a * elemvect(i);
572  }
573 }
574 
575 void Vector::SetSubVectorComplement(const Array<int> &dofs, const double val)
576 {
577  Vector dofs_vals;
578  GetSubVector(dofs, dofs_vals);
579  operator=(val);
580  SetSubVector(dofs, dofs_vals);
581 }
582 
583 void Vector::Print(std::ostream &out, int width) const
584 {
585  if (!size) { return; }
586 
587  for (int i = 0; 1; )
588  {
589  out << data[i];
590  i++;
591  if (i == size)
592  {
593  break;
594  }
595  if ( i % width == 0 )
596  {
597  out << '\n';
598  }
599  else
600  {
601  out << ' ';
602  }
603  }
604  out << '\n';
605 }
606 
607 void Vector::Print_HYPRE(std::ostream &out) const
608 {
609  int i;
610  std::ios::fmtflags old_fmt = out.flags();
611  out.setf(std::ios::scientific);
612  std::streamsize old_prec = out.precision(14);
613 
614  out << size << '\n'; // number of rows
615 
616  for (i = 0; i < size; i++)
617  {
618  out << data[i] << '\n';
619  }
620 
621  out.precision(old_prec);
622  out.flags(old_fmt);
623 }
624 
625 void Vector::Randomize(int seed)
626 {
627  // static unsigned int seed = time(0);
628  const double max = (double)(RAND_MAX) + 1.;
629 
630  if (seed == 0)
631  {
632  seed = (int)time(0);
633  }
634 
635  // srand(seed++);
636  srand((unsigned)seed);
637 
638  for (int i = 0; i < size; i++)
639  {
640  data[i] = fabs(rand()/max);
641  }
642 }
643 
644 double Vector::Norml2() const
645 {
646  return sqrt((*this)*(*this));
647 }
648 
649 double Vector::Normlinf() const
650 {
651  double max = fabs(data[0]);
652 
653  for (int i = 1; i < size; i++)
654  if (fabs(data[i]) > max)
655  {
656  max = fabs(data[i]);
657  }
658 
659  return max;
660 }
661 
662 double Vector::Norml1() const
663 {
664  double sum = 0.0;
665 
666  for (int i = 0; i < size; i++)
667  {
668  sum += fabs (data[i]);
669  }
670 
671  return sum;
672 }
673 
674 double Vector::Normlp(double p) const
675 {
676  MFEM_ASSERT(p > 0.0, "Vector::Normlp");
677  if (p == 1.0)
678  {
679  return Norml1();
680  }
681  if (p == 2.0)
682  {
683  return Norml2();
684  }
685  if (p < std::numeric_limits<double>::infinity())
686  {
687  double sum = 0.0;
688  for (int i = 0; i < size; i++)
689  {
690  sum += pow(fabs(data[i]), p);
691  }
692  return pow(sum, 1.0/p);
693  }
694  else
695  {
696  return Normlinf();
697  }
698 }
699 
700 double Vector::Max() const
701 {
702  double max = data[0];
703 
704  for (int i = 1; i < size; i++)
705  if (data[i] > max)
706  {
707  max = data[i];
708  }
709 
710  return max;
711 }
712 
713 double Vector::Min() const
714 {
715  double min = data[0];
716 
717  for (int i = 1; i < size; i++)
718  if (data[i] < min)
719  {
720  min = data[i];
721  }
722 
723  return min;
724 }
725 
726 double Vector::Sum() const
727 {
728  double sum = 0.0;
729 
730  for (int i = 0; i < size; i++)
731  {
732  sum += data[i];
733  }
734 
735  return sum;
736 }
737 
738 double Vector::DistanceTo(const double *p) const
739 {
740  return Distance(data, p, size);
741 }
742 
743 }
int Size() const
Logical size of the array.
Definition: array.hpp:109
void SetVector(const Vector &v, int offset)
Definition: vector.cpp:233
Vector()
Default constructor for Vector. Sets size = 0 and data = NULL.
Definition: vector.hpp:43
double & Elem(int i)
Sets value in vector. Index i = 0 .. size-1.
Definition: vector.cpp:75
void Print(std::ostream &out=std::cout, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:583
void SetSize(int s)
Resizes the vector if the new size is different.
Definition: vector.hpp:259
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:644
double & operator()(int i)
Sets value in vector. Index i = 0 .. size-1.
Definition: vector.hpp:288
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Definition: vector.cpp:457
int Size() const
Returns the size of the vector.
Definition: vector.hpp:84
double Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:649
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:625
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:259
double operator*(const double *) const
Definition: vector.cpp:85
Vector & operator=(const double *v)
Definition: vector.cpp:112
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Definition: vector.cpp:45
double Normlp(double p) const
Returns the l_p norm of the vector.
Definition: vector.cpp:674
int dim
Definition: ex3.cpp:48
void median(const Vector &lo, const Vector &hi)
v = median(v,lo,hi) entrywise. Implementation assumes lo &lt;= hi.
Definition: vector.cpp:440
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add (element) subvector to the vector.
Definition: vector.cpp:527
Vector & operator/=(double c)
Definition: vector.cpp:151
double DistanceTo(const double *p) const
Compute the Euclidean distance to another vector.
Definition: vector.cpp:738
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the &#39;dofs&#39; array to the given &#39;val&#39;.
Definition: vector.cpp:575
Vector & operator*=(double c)
Definition: vector.cpp:142
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:713
double Norml1() const
Returns the l_1 norm of the vector.
Definition: vector.cpp:662
double Distance(const double *x, const double *y, const int n)
Definition: vector.hpp:325
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:385
void mfem_error(const char *msg)
Definition: error.cpp:23
void SetSubVector(const Array< int > &dofs, const Vector &elemvect)
Definition: vector.cpp:493
void Print_HYPRE(std::ostream &out) const
Prints vector to stream out in HYPRE_Vector format.
Definition: vector.cpp:607
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:218
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:200
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:700
Vector & operator-=(double c)
Definition: vector.cpp:161
Vector data type.
Definition: vector.hpp:33
int allocsize
Definition: vector.hpp:37
Vector & operator+=(const Vector &v)
Definition: vector.cpp:185
double * data
Definition: vector.hpp:38
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:726
void Neg()
(*this) = -(*this)
Definition: vector.cpp:251