dense_vector.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2025 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 boostcopy::enable_if_c<
    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 boostcopy::enable_if_c<
  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];
}

} // 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