dense_vector.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
 
 
#ifndef LIBMESH_DENSE_VECTOR_H
#define LIBMESH_DENSE_VECTOR_H
 
// Local Includes
#include "libmesh/libmesh_common.h"
#include "libmesh/dense_vector_base.h"
#include "libmesh/compare_types.h"
#include "libmesh/int_range.h"
#include "libmesh/tensor_tools.h"
 
#ifdef LIBMESH_HAVE_EIGEN
#include "libmesh/ignore_warnings.h"
#include <Eigen/Core>
#include "libmesh/restore_warnings.h"
#endif
 
#ifdef LIBMESH_HAVE_METAPHYSICL
#include "metaphysicl/raw_type.h"
#endif
 
// C++ includes
#include <vector>
#include <initializer_list>
 
namespace libMesh
{
 
template<typename T>
class DenseVector : public DenseVectorBase<T>
{
public:
 
  explicit
  DenseVector(const unsigned int n=0);
 
  explicit
  DenseVector(const unsigned int n,
              const T & val);
 
  template <typename T2>
  DenseVector (const DenseVector<T2> & other_vector);
 
  template <typename T2>
  DenseVector (const std::vector<T2> & other_vector);
 
  template <typename T2>
  DenseVector (std::initializer_list<T2> init_list);
 
  DenseVector (DenseVector &&) = default;
  DenseVector (const DenseVector &) = default;
  DenseVector & operator= (const DenseVector &) = default;
  DenseVector & operator= (DenseVector &&) = default;
  virtual ~DenseVector() = default;
 
  virtual unsigned int size() const override final
  {
    return cast_int<unsigned int>(_val.size());
  }
 
  virtual bool empty() const override final
  { return _val.empty(); }
 
  virtual void zero() override final;
 
  const T & operator() (const unsigned int i) const;
 
  T & operator() (const unsigned int i);
 
  const T & operator[] (const unsigned int i) const { return (*this)(i); }
 
  T & operator[] (const unsigned int i) { return (*this)(i); }
 
  virtual T el(const unsigned int i) const override final
  { return (*this)(i); }
 
  virtual T & el(const unsigned int i) override final
  { return (*this)(i); }
 
  template <typename T2>
  DenseVector<T> & operator = (const DenseVector<T2> & other_vector);
 
  void swap(DenseVector<T> & other_vector);
 
  void resize (const unsigned int n);
 
  template <typename T2>
  void append (const DenseVector<T2> & other_vector);
 
  void scale (const T factor);
 
  DenseVector<T> & operator*= (const T factor);
 
  template <typename T2, typename T3>
  typename std::enable_if<
    ScalarTraits<T2>::value, void >::type
  add (const T2 factor,
       const DenseVector<T3> & vec);
 
  template <typename T2>
  typename CompareTypes<T, T2>::supertype dot (const DenseVector<T2> & vec) const;
 
  template <typename T2>
  typename CompareTypes<T, T2>::supertype indefinite_dot (const DenseVector<T2> & vec) const;
 
  template <typename T2>
  bool operator== (const DenseVector<T2> & vec) const;
 
  template <typename T2>
  bool operator!= (const DenseVector<T2> & vec) const;
 
  template <typename T2>
  DenseVector<T> & operator+= (const DenseVector<T2> & vec);
 
  template <typename T2>
  DenseVector<T> & operator-= (const DenseVector<T2> & vec);
 
  Real min () const;
 
  Real max () const;
 
  Real l1_norm () const;
 
  Real l2_norm () const;
 
  Real linfty_norm () const;
 
  void get_principal_subvector (unsigned int sub_n, DenseVector<T> & dest) const;
 
  std::vector<T> & get_values() { return _val; }
 
  const std::vector<T> & get_values() const { return _val; }
 
  typename std::vector<T>::const_iterator begin() const { return _val.begin(); }
  typename std::vector<T>::iterator begin() { return _val.begin(); }
 
  typename std::vector<T>::const_iterator end() const { return _val.end(); }
  typename std::vector<T>::iterator end() { return _val.end(); }
 
private:
 
  std::vector<T> _val;
};
 
 
 
// ------------------------------------------------------------
// DenseVector member functions
template<typename T>
inline
DenseVector<T>::DenseVector(const unsigned int n) :
  _val (n, T{})
{
}
 
 
template<typename T>
inline
DenseVector<T>::DenseVector(const unsigned int n,
                            const T & val) :
  _val (n, val)
{
}
 
 
 
template<typename T>
template<typename T2>
inline
DenseVector<T>::DenseVector (const DenseVector<T2> & other_vector) :
  DenseVectorBase<T>()
{
  const std::vector<T2> & other_vals = other_vector.get_values();
 
  _val.clear();
 
  const int N = cast_int<int>(other_vals.size());
  _val.reserve(N);
 
  for (int i=0; i<N; i++)
    _val.push_back(other_vals[i]);
}
 
 
 
template<typename T>
template<typename T2>
inline
DenseVector<T>::DenseVector (const std::vector<T2> & other_vector) :
  _val(other_vector)
{
}
 
 
template<typename T>
template <typename T2>
inline
DenseVector<T>::DenseVector (std::initializer_list<T2> init_list) :
  _val(init_list.begin(), init_list.end())
{
}
 
 
 
template<typename T>
template<typename T2>
inline
DenseVector<T> & DenseVector<T>::operator = (const DenseVector<T2> & other_vector)
{
  const std::vector<T2> & other_vals = other_vector.get_values();
 
  _val.clear();
 
  const int N = cast_int<int>(other_vals.size());
  _val.reserve(N);
 
  for (int i=0; i<N; i++)
    _val.push_back(other_vals[i]);
 
  return *this;
}
 
 
 
template<typename T>
inline
void DenseVector<T>::swap(DenseVector<T> & other_vector)
{
  _val.swap(other_vector._val);
}
 
 
 
template<typename T>
inline
void DenseVector<T>::resize(const unsigned int n)
{
  _val.resize(n);
 
  zero();
}
 
 
 
template<typename T>
template<typename T2>
inline
void DenseVector<T>::append (const DenseVector<T2> & other_vector)
{
  const std::vector<T2> & other_vals = other_vector.get_values();
 
  _val.reserve(this->size() + other_vals.size());
  _val.insert(_val.end(), other_vals.begin(), other_vals.end());
}
 
 
 
template<typename T>
inline
void DenseVector<T>::zero()
{
  std::fill (_val.begin(),
             _val.end(),
             T{});
}
 
 
 
template<typename T>
inline
const T & DenseVector<T>::operator () (const unsigned int i) const
{
  libmesh_assert_less (i, _val.size());
 
  return _val[i];
}
 
 
 
template<typename T>
inline
T & DenseVector<T>::operator () (const unsigned int i)
{
  libmesh_assert_less (i, _val.size());
 
  return _val[i];
}
 
 
 
template<typename T>
inline
void DenseVector<T>::scale (const T factor)
{
  const int N = cast_int<int>(_val.size());
  for (int i=0; i<N; i++)
    _val[i] *= factor;
}
 
 
 
template<typename T>
inline
DenseVector<T> & DenseVector<T>::operator*= (const T factor)
{
  this->scale(factor);
  return *this;
}
 
 
 
template<typename T>
template<typename T2, typename T3>
inline
typename std::enable_if<
  ScalarTraits<T2>::value, void >::type
DenseVector<T>::add (const T2 factor,
                     const DenseVector<T3> & vec)
{
  libmesh_assert_equal_to (this->size(), vec.size());
 
  const int N = cast_int<int>(_val.size());
  for (int i=0; i<N; i++)
    (*this)(i) += static_cast<T>(factor)*vec(i);
}
 
 
 
template<typename T>
template<typename T2>
inline
typename CompareTypes<T, T2>::supertype DenseVector<T>::dot (const DenseVector<T2> & vec) const
{
  if (!_val.size())
    return 0.;
 
  libmesh_assert_equal_to (this->size(), vec.size());
 
#ifdef LIBMESH_HAVE_EIGEN
  // We reverse the order of the arguments to dot() here since
  // the convention in Eigen is to take the complex conjugate of the
  // *first* argument, while ours is to take the complex conjugate of
  // the second.
  return Eigen::Map<const typename Eigen::Matrix<T2, Eigen::Dynamic, 1>>(vec.get_values().data(), vec.size())
    .dot(Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()));
#else
  typename CompareTypes<T, T2>::supertype val = 0.;
 
  const int N = cast_int<int>(_val.size());
  // The following pragma tells clang's vectorizer that it is safe to
  // reorder floating point operations for this loop.
#ifdef __clang__
#pragma clang loop vectorize(enable)
#endif
  for (int i=0; i<N; i++)
    val += (*this)(i)*libmesh_conj(vec(i));
 
  return val;
#endif
}
 
template<typename T>
template<typename T2>
inline
typename CompareTypes<T, T2>::supertype DenseVector<T>::indefinite_dot (const DenseVector<T2> & vec) const
{
  libmesh_assert_equal_to (this->size(), vec.size());
 
  typename CompareTypes<T, T2>::supertype val = 0.;
 
  const int N = cast_int<int>(_val.size());
  for (int i=0; i<N; i++)
    val += (*this)(i)*(vec(i));
 
  return val;
}
 
template<typename T>
template<typename T2>
inline
bool DenseVector<T>::operator== (const DenseVector<T2> & vec) const
{
  libmesh_assert_equal_to (this->size(), vec.size());
 
  const int N = cast_int<int>(_val.size());
  for (int i=0; i<N; i++)
    if ((*this)(i) != vec(i))
      return false;
 
  return true;
}
 
 
 
template<typename T>
template<typename T2>
inline
bool DenseVector<T>::operator!= (const DenseVector<T2> & vec) const
{
  libmesh_assert_equal_to (this->size(), vec.size());
 
  const int N = cast_int<int>(_val.size());
  for (int i=0; i<N; i++)
    if ((*this)(i) != vec(i))
      return true;
 
  return false;
}
 
 
 
template<typename T>
template<typename T2>
inline
DenseVector<T> & DenseVector<T>::operator+= (const DenseVector<T2> & vec)
{
  libmesh_assert_equal_to (this->size(), vec.size());
 
  const int N = cast_int<int>(_val.size());
  for (int i=0; i<N; i++)
    (*this)(i) += vec(i);
 
  return *this;
}
 
 
 
template<typename T>
template<typename T2>
inline
DenseVector<T> & DenseVector<T>::operator-= (const DenseVector<T2> & vec)
{
  libmesh_assert_equal_to (this->size(), vec.size());
 
  const int N = cast_int<int>(_val.size());
  for (int i=0; i<N; i++)
    (*this)(i) -= vec(i);
 
  return *this;
}
 
 
 
template<typename T>
inline
Real DenseVector<T>::min () const
{
  libmesh_assert (this->size());
  Real my_min = libmesh_real((*this)(0));
 
  const int N = cast_int<int>(_val.size());
  for (int i=1; i!=N; i++)
    {
      Real current = libmesh_real((*this)(i));
      my_min = (my_min < current? my_min : current);
    }
  return my_min;
}
 
 
 
template<typename T>
inline
Real DenseVector<T>::max () const
{
  libmesh_assert (this->size());
  Real my_max = libmesh_real((*this)(0));
 
  const int N = cast_int<int>(_val.size());
  for (int i=1; i!=N; i++)
    {
      Real current = libmesh_real((*this)(i));
      my_max = (my_max > current? my_max : current);
    }
  return my_max;
}
 
 
 
template<typename T>
inline
Real DenseVector<T>::l1_norm () const
{
  if (!_val.size())
    return 0.;
 
#ifdef LIBMESH_HAVE_EIGEN
  return Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()).template lpNorm<1>();
#else
  Real my_norm = 0.;
  const int N = cast_int<int>(_val.size());
  for (int i=0; i!=N; i++)
    my_norm += std::abs((*this)(i));
 
  return my_norm;
#endif
}
 
 
 
template<typename T>
inline
Real DenseVector<T>::l2_norm () const
{
  if (!_val.size())
    return 0.;
 
#ifdef LIBMESH_HAVE_EIGEN
  return Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()).norm();
#else
  Real my_norm = 0.;
  const int N = cast_int<int>(_val.size());
  // The following pragma tells clang's vectorizer that it is safe to
  // reorder floating point operations for this loop.
#ifdef __clang__
#pragma clang loop vectorize(enable)
#endif
  for (int i=0; i!=N; i++)
    my_norm += TensorTools::norm_sq((*this)(i));
 
  return sqrt(my_norm);
#endif
}
 
 
 
template<typename T>
inline
Real DenseVector<T>::linfty_norm () const
{
  if (!_val.size())
    return 0.;
 
#ifdef LIBMESH_HAVE_EIGEN
  return Eigen::Map<const typename Eigen::Matrix<T, Eigen::Dynamic, 1>>(_val.data(), _val.size()).template lpNorm<Eigen::Infinity>();
#else
  Real my_norm = TensorTools::norm_sq((*this)(0));
 
  const int N = cast_int<int>(_val.size());
  for (int i=1; i!=N; i++)
    {
      Real current = TensorTools::norm_sq((*this)(i));
      my_norm = (my_norm > current? my_norm : current);
    }
  return sqrt(my_norm);
#endif
}
 
 
 
template<typename T>
inline
void DenseVector<T>::get_principal_subvector (unsigned int sub_n,
                                              DenseVector<T> & dest) const
{
  libmesh_assert_less_equal ( sub_n, this->size() );
 
  dest.resize(sub_n);
  const int N = cast_int<int>(sub_n);
  for (int i=0; i<N; i++)
    dest(i) = _val[i];
}
 
 
// A vector is finite iff every component is
template <typename T>
bool isfinite (const DenseVector<T> & var)
{
  using std::isfinite;
  using libMesh::isfinite; // for T==complex
  for (auto i : index_range(var))
    if (!isfinite(var(i)))
      return false;
  return true;
}
 
 
// A vector is infinite iff some component is infinite but no
// component is NaN.
//
// This is arguably inconsistent with our std::complex overload (and
// the C99 Annex G recommendations for _Complex, and C++ std::complex
// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but
// this is probably safer for users.
template <typename T>
bool isinf (const DenseVector<T> & var)
{
  using std::isinf;
  using libMesh::isinf; // for T==complex
  using std::isnan;
  using libMesh::isnan;
  bool has_inf = false;
  for (auto i : index_range(var))
    {
      // NaN anywhere makes us NaN, not inf
      if (isnan(var(i)))
        return false;
      has_inf = has_inf || isinf(var(i));
    }
  return has_inf;
}
 
 
// A vector is NaN iff some component is NaN
//
// This is arguably inconsistent with our std::complex overload (and
// the C99 Annex G recommendations for _Complex, and C++ std::complex
// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but
// this is probably safer for users.
template <typename T>
bool isnan (const DenseVector<T> & var)
{
  using std::isnan;
  using libMesh::isnan; // for T==complex
  for (auto i : index_range(var))
    if (isnan(var(i)))
      return true;
  return false;
}
 
 
} // namespace libMesh
 
#ifdef LIBMESH_HAVE_METAPHYSICL
namespace MetaPhysicL
{
template <typename T>
struct RawType<libMesh::DenseVector<T>>
{
  typedef libMesh::DenseVector<typename RawType<T>::value_type> value_type;
 
  static value_type value (const libMesh::DenseVector<T> & in)
    {
      const auto s = in.size();
      value_type ret(s);
      for (unsigned int i = 0; i < s; ++i)
          ret(i) = raw_value(in(i));
 
      return ret;
    }
};
}
#endif
 
#endif // LIBMESH_DENSE_VECTOR_H