dense_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 #ifndef LIBMESH_DENSE_VECTOR_H
21 #define LIBMESH_DENSE_VECTOR_H
22 
23 // Local Includes
24 #include "libmesh/libmesh_common.h"
25 #include "libmesh/compare_types.h"
27 #include "libmesh/tensor_tools.h"
28 
29 #ifdef LIBMESH_HAVE_EIGEN
31 #include <Eigen/Core>
33 #endif
34 
35 // C++ includes
36 #include <vector>
37 
38 namespace libMesh
39 {
40 
51 template<typename T>
52 class DenseVector : public DenseVectorBase<T>
53 {
54 public:
55 
59  explicit
60  DenseVector(const unsigned int n=0);
61 
66  explicit
67  DenseVector(const unsigned int n,
68  const T & val);
69 
73  template <typename T2>
74  DenseVector (const DenseVector<T2> & other_vector);
75 
79  template <typename T2>
80  DenseVector (const std::vector<T2> & other_vector);
81 
86  DenseVector (DenseVector &&) = default;
87  DenseVector (const DenseVector &) = default;
88  DenseVector & operator= (const DenseVector &) = default;
89  DenseVector & operator= (DenseVector &&) = default;
90  virtual ~DenseVector() = default;
91 
92  virtual unsigned int size() const override
93  {
94  return cast_int<unsigned int>(_val.size());
95  }
96 
97  virtual bool empty() const override
98  { return _val.empty(); }
99 
100  virtual void zero() override;
101 
105  const T & operator() (const unsigned int i) const;
106 
110  T & operator() (const unsigned int i);
111 
112  virtual T el(const unsigned int i) const override
113  { return (*this)(i); }
114 
115  virtual T & el(const unsigned int i) override
116  { return (*this)(i); }
117 
123  template <typename T2>
124  DenseVector<T> & operator = (const DenseVector<T2> & other_vector);
125 
129  void swap(DenseVector<T> & other_vector);
130 
134  void resize (const unsigned int n);
135 
140  template <typename T2>
141  void append (const DenseVector<T2> & other_vector);
142 
146  void scale (const T factor);
147 
153  DenseVector<T> & operator*= (const T factor);
154 
162  template <typename T2, typename T3>
163  typename boostcopy::enable_if_c<
164  ScalarTraits<T2>::value, void >::type
165  add (const T2 factor,
166  const DenseVector<T3> & vec);
167 
173  template <typename T2>
174  typename CompareTypes<T, T2>::supertype dot (const DenseVector<T2> & vec) const;
175 
182  template <typename T2>
184 
188  template <typename T2>
189  bool operator== (const DenseVector<T2> & vec) const;
190 
194  template <typename T2>
195  bool operator!= (const DenseVector<T2> & vec) const;
196 
202  template <typename T2>
204 
210  template <typename T2>
212 
217  Real min () const;
218 
223  Real max () const;
224 
229  Real l1_norm () const;
230 
235  Real l2_norm () const;
236 
241  Real linfty_norm () const;
242 
247  void get_principal_subvector (unsigned int sub_n, DenseVector<T> & dest) const;
248 
256  std::vector<T> & get_values() { return _val; }
257 
261  const std::vector<T> & get_values() const { return _val; }
262 
263 private:
264 
268  std::vector<T> _val;
269 };
270 
271 
272 
273 // ------------------------------------------------------------
274 // DenseVector member functions
275 template<typename T>
276 inline
277 DenseVector<T>::DenseVector(const unsigned int n) :
278  _val (n, T(0.))
279 {
280 }
281 
282 
283 template<typename T>
284 inline
285 DenseVector<T>::DenseVector(const unsigned int n,
286  const T & val) :
287  _val (n, val)
288 {
289 }
290 
291 
292 
293 template<typename T>
294 template<typename T2>
295 inline
297  DenseVectorBase<T>()
298 {
299  const std::vector<T2> & other_vals = other_vector.get_values();
300 
301  _val.clear();
302 
303  const int N = cast_int<int>(other_vals.size());
304  _val.reserve(N);
305 
306  for (int i=0; i<N; i++)
307  _val.push_back(other_vals[i]);
308 }
309 
310 
311 
312 template<typename T>
313 template<typename T2>
314 inline
315 DenseVector<T>::DenseVector (const std::vector<T2> & other_vector) :
316  _val(other_vector)
317 {
318 }
319 
320 
321 
322 
323 
324 template<typename T>
325 template<typename T2>
326 inline
328 {
329  const std::vector<T2> & other_vals = other_vector.get_values();
330 
331  _val.clear();
332 
333  const int N = cast_int<int>(other_vals.size());
334  _val.reserve(N);
335 
336  for (int i=0; i<N; i++)
337  _val.push_back(other_vals[i]);
338 
339  return *this;
340 }
341 
342 
343 
344 template<typename T>
345 inline
347 {
348  _val.swap(other_vector._val);
349 }
350 
351 
352 
353 template<typename T>
354 inline
355 void DenseVector<T>::resize(const unsigned int n)
356 {
357  _val.resize(n);
358 
359  zero();
360 }
361 
362 
363 
364 template<typename T>
365 template<typename T2>
366 inline
367 void DenseVector<T>::append (const DenseVector<T2> & other_vector)
368 {
369  const std::vector<T2> & other_vals = other_vector.get_values();
370 
371  _val.reserve(this->size() + other_vals.size());
372  _val.insert(_val.end(), other_vals.begin(), other_vals.end());
373 }
374 
375 
376 
377 template<typename T>
378 inline
380 {
381  std::fill (_val.begin(),
382  _val.end(),
383  T(0.));
384 }
385 
386 
387 
388 template<typename T>
389 inline
390 const T & DenseVector<T>::operator () (const unsigned int i) const
391 {
392  libmesh_assert_less (i, _val.size());
393 
394  return _val[i];
395 }
396 
397 
398 
399 template<typename T>
400 inline
401 T & DenseVector<T>::operator () (const unsigned int i)
402 {
403  libmesh_assert_less (i, _val.size());
404 
405  return _val[i];
406 }
407 
408 
409 
410 template<typename T>
411 inline
412 void DenseVector<T>::scale (const T factor)
413 {
414  const int N = cast_int<int>(_val.size());
415  for (int i=0; i<N; i++)
416  _val[i] *= factor;
417 }
418 
419 
420 
421 template<typename T>
422 inline
424 {
425  this->scale(factor);
426  return *this;
427 }
428 
429 
430 
431 template<typename T>
432 template<typename T2, typename T3>
433 inline
434 typename boostcopy::enable_if_c<
435  ScalarTraits<T2>::value, void >::type
436 DenseVector<T>::add (const T2 factor,
437  const DenseVector<T3> & vec)
438 {
439  libmesh_assert_equal_to (this->size(), vec.size());
440 
441  const int N = cast_int<int>(_val.size());
442  for (int i=0; i<N; i++)
443  (*this)(i) += static_cast<T>(factor)*vec(i);
444 }
445 
446 
447 
448 template<typename T>
449 template<typename T2>
450 inline
452 {
453  if (!_val.size())
454  return 0.;
455 
456  libmesh_assert_equal_to (this->size(), vec.size());
457 
458 #ifdef LIBMESH_HAVE_EIGEN
459  // We reverse the order of the arguments to dot() here since
460  // the convention in Eigen is to take the complex conjugate of the
461  // *first* argument, while ours is to take the complex conjugate of
462  // the second.
463  return Eigen::Map<const typename Eigen::Matrix<T2, Eigen::Dynamic, 1>>(vec.get_values().data(), vec.size())
464  .dot(Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()));
465 #else
466  typename CompareTypes<T, T2>::supertype val = 0.;
467 
468  const int N = cast_int<int>(_val.size());
469  // The following pragma tells clang's vectorizer that it is safe to
470  // reorder floating point operations for this loop.
471 #ifdef __clang__
472 #pragma clang loop vectorize(enable)
473 #endif
474  for (int i=0; i<N; i++)
475  val += (*this)(i)*libmesh_conj(vec(i));
476 
477  return val;
478 #endif
479 }
480 
481 template<typename T>
482 template<typename T2>
483 inline
485 {
486  libmesh_assert_equal_to (this->size(), vec.size());
487 
488  typename CompareTypes<T, T2>::supertype val = 0.;
489 
490  const int N = cast_int<int>(_val.size());
491  for (int i=0; i<N; i++)
492  val += (*this)(i)*(vec(i));
493 
494  return val;
495 }
496 
497 template<typename T>
498 template<typename T2>
499 inline
501 {
502  libmesh_assert_equal_to (this->size(), vec.size());
503 
504  const int N = cast_int<int>(_val.size());
505  for (int i=0; i<N; i++)
506  if ((*this)(i) != vec(i))
507  return false;
508 
509  return true;
510 }
511 
512 
513 
514 template<typename T>
515 template<typename T2>
516 inline
518 {
519  libmesh_assert_equal_to (this->size(), vec.size());
520 
521  const int N = cast_int<int>(_val.size());
522  for (int i=0; i<N; i++)
523  if ((*this)(i) != vec(i))
524  return true;
525 
526  return false;
527 }
528 
529 
530 
531 template<typename T>
532 template<typename T2>
533 inline
535 {
536  libmesh_assert_equal_to (this->size(), vec.size());
537 
538  const int N = cast_int<int>(_val.size());
539  for (int i=0; i<N; i++)
540  (*this)(i) += vec(i);
541 
542  return *this;
543 }
544 
545 
546 
547 template<typename T>
548 template<typename T2>
549 inline
551 {
552  libmesh_assert_equal_to (this->size(), vec.size());
553 
554  const int N = cast_int<int>(_val.size());
555  for (int i=0; i<N; i++)
556  (*this)(i) -= vec(i);
557 
558  return *this;
559 }
560 
561 
562 
563 template<typename T>
564 inline
566 {
567  libmesh_assert (this->size());
568  Real my_min = libmesh_real((*this)(0));
569 
570  const int N = cast_int<int>(_val.size());
571  for (int i=1; i!=N; i++)
572  {
573  Real current = libmesh_real((*this)(i));
574  my_min = (my_min < current? my_min : current);
575  }
576  return my_min;
577 }
578 
579 
580 
581 template<typename T>
582 inline
584 {
585  libmesh_assert (this->size());
586  Real my_max = libmesh_real((*this)(0));
587 
588  const int N = cast_int<int>(_val.size());
589  for (int i=1; i!=N; i++)
590  {
591  Real current = libmesh_real((*this)(i));
592  my_max = (my_max > current? my_max : current);
593  }
594  return my_max;
595 }
596 
597 
598 
599 template<typename T>
600 inline
602 {
603  if (!_val.size())
604  return 0.;
605 
606 #ifdef LIBMESH_HAVE_EIGEN
607  return Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()).template lpNorm<1>();
608 #else
609  Real my_norm = 0.;
610  const int N = cast_int<int>(_val.size());
611  for (int i=0; i!=N; i++)
612  my_norm += std::abs((*this)(i));
613 
614  return my_norm;
615 #endif
616 }
617 
618 
619 
620 template<typename T>
621 inline
623 {
624  if (!_val.size())
625  return 0.;
626 
627 #ifdef LIBMESH_HAVE_EIGEN
628  return Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()).norm();
629 #else
630  Real my_norm = 0.;
631  const int N = cast_int<int>(_val.size());
632  // The following pragma tells clang's vectorizer that it is safe to
633  // reorder floating point operations for this loop.
634 #ifdef __clang__
635 #pragma clang loop vectorize(enable)
636 #endif
637  for (int i=0; i!=N; i++)
638  my_norm += TensorTools::norm_sq((*this)(i));
639 
640  return sqrt(my_norm);
641 #endif
642 }
643 
644 
645 
646 template<typename T>
647 inline
649 {
650  if (!_val.size())
651  return 0.;
652 
653 #ifdef LIBMESH_HAVE_EIGEN
654  return Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()).template lpNorm<Eigen::Infinity>();
655 #else
656  Real my_norm = TensorTools::norm_sq((*this)(0));
657 
658  const int N = cast_int<int>(_val.size());
659  for (int i=1; i!=N; i++)
660  {
661  Real current = TensorTools::norm_sq((*this)(i));
662  my_norm = (my_norm > current? my_norm : current);
663  }
664  return sqrt(my_norm);
665 #endif
666 }
667 
668 
669 
670 template<typename T>
671 inline
673  DenseVector<T> & dest) const
674 {
675  libmesh_assert_less_equal ( sub_n, this->size() );
676 
677  dest.resize(sub_n);
678  const int N = cast_int<int>(sub_n);
679  for (int i=0; i<N; i++)
680  dest(i) = _val[i];
681 }
682 
683 } // namespace libMesh
684 
685 #endif // LIBMESH_DENSE_VECTOR_H
T libmesh_real(T a)
virtual unsigned int size() const override
Definition: dense_vector.h:92
DenseVector< T > & operator*=(const T factor)
Definition: dense_vector.h:423
double abs(double a)
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseVector< T3 > &vec)
Definition: dense_vector.h:436
T libmesh_conj(T a)
void scale(MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
void resize(const unsigned int n)
Definition: dense_vector.h:355
std::vector< T > & get_values()
Definition: dense_vector.h:256
bool operator==(const DenseVector< T2 > &vec) const
Definition: dense_vector.h:500
virtual T el(const unsigned int i) const override
Definition: dense_vector.h:112
const Number zero
Definition: libmesh.h:239
void scale(const T factor)
Definition: dense_vector.h:412
long double max(long double a, double b)
void get_principal_subvector(unsigned int sub_n, DenseVector< T > &dest) const
Definition: dense_vector.h:672
const std::vector< T > & get_values() const
Definition: dense_vector.h:261
virtual bool empty() const override
Definition: dense_vector.h:97
void swap(DenseVector< T > &other_vector)
Definition: dense_vector.h:346
DenseVector & operator=(const DenseVector &)=default
Real linfty_norm() const
Definition: dense_vector.h:648
bool operator!=(const variant_filter_iterator< Predicate, Type, ReferenceType, PointerType > &x, const variant_filter_iterator< Predicate, Type, ReferenceType, PointerType > &y)
Real l2_norm() const
Definition: dense_vector.h:622
virtual ~DenseVector()=default
virtual T & el(const unsigned int i) override
Definition: dense_vector.h:115
std::vector< T > _val
Definition: dense_vector.h:268
bool operator!=(const DenseVector< T2 > &vec) const
Definition: dense_vector.h:517
const T & operator()(const unsigned int i) const
Definition: dense_vector.h:390
CompareTypes< T, T2 >::supertype indefinite_dot(const DenseVector< T2 > &vec) const
Definition: dense_vector.h:484
CompareTypes< T, T2 >::supertype dot(const DenseVector< T2 > &vec) const
Definition: dense_vector.h:451
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void swap(Iterator &lhs, Iterator &rhs)
Real l1_norm() const
Definition: dense_vector.h:601
DenseVector(const unsigned int n=0)
Definition: dense_vector.h:277
bool operator==(const variant_filter_iterator< Predicate, Type, ReferenceType, PointerType > &x, const variant_filter_iterator< Predicate, Type, ReferenceType, PointerType > &y)
void append(const DenseVector< T2 > &other_vector)
Definition: dense_vector.h:367
Iterator & operator=(const Iterator &rhs)
DenseVector< T > & operator+=(const DenseVector< T2 > &vec)
Definition: dense_vector.h:534
long double min(long double a, double b)
DenseVector< T > & operator-=(const DenseVector< T2 > &vec)
Definition: dense_vector.h:550
virtual void zero() override
Definition: dense_vector.h:379