numeric_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_NUMERIC_VECTOR_H
#define LIBMESH_NUMERIC_VECTOR_H

// Local includes
#include "libmesh/libmesh_common.h"
#include "libmesh/enum_parallel_type.h" // AUTOMATIC
#include "libmesh/id_types.h"
#include "libmesh/int_range.h"
#include "libmesh/reference_counted_object.h"
#include "libmesh/libmesh.h"
#include "libmesh/parallel_object.h"
#include "libmesh/dense_subvector.h"
#include "libmesh/dense_vector.h"

// C++ includes
#include <cstddef>
#include <set>
#include <vector>
#include <memory>
#include <mutex>

namespace libMesh
{


// forward declarations
template <typename T> class NumericVector;
template <typename T> class DenseVector;
template <typename T> class DenseSubVector;
template <typename T> class SparseMatrix;
template <typename T> class ShellMatrix;
enum SolverPackage : int;

template <typename T>
class NumericVector : public ReferenceCountedObject<NumericVector<T>>,
                      public ParallelObject
{
public:

  explicit
  NumericVector (const Parallel::Communicator & comm_in,
                 const ParallelType ptype = AUTOMATIC);

  explicit
  NumericVector (const Parallel::Communicator & comm_in,
                 const numeric_index_type n,
                 const ParallelType ptype = AUTOMATIC);

  NumericVector (const Parallel::Communicator & comm_in,
                 const numeric_index_type n,
                 const numeric_index_type n_local,
                 const ParallelType ptype = AUTOMATIC);

  NumericVector (const Parallel::Communicator & comm_in,
                 const numeric_index_type N,
                 const numeric_index_type n_local,
                 const std::vector<numeric_index_type> & ghost,
                 const ParallelType ptype = AUTOMATIC);

  virtual NumericVector<T> & operator= (const NumericVector<T> & v) = 0;

  NumericVector (NumericVector &&) = default;
  NumericVector (const NumericVector &) = default;
  NumericVector & operator= (NumericVector &&) = default;

  virtual ~NumericVector() = default;

  static std::unique_ptr<NumericVector<T>>
  build(const Parallel::Communicator & comm,
        SolverPackage solver_package = libMesh::default_solver_package(),
        ParallelType parallel_type = AUTOMATIC);

  virtual bool initialized() const { return _is_initialized; }

  ParallelType type() const { return _type; }

  bool is_effectively_serial() const
  { return (this->n_processors()==1) || (_type == SERIAL); }

  bool is_effectively_ghosted() const
  { return (this->n_processors()!=1) && (_type == GHOSTED); }

#ifdef LIBMESH_ENABLE_DEPRECATED
  ParallelType & type() { return _type; }
#endif

  void set_type(ParallelType t);

  virtual bool closed() const { return _is_closed; }

  virtual void close () = 0;

  virtual void clear ();

  virtual void zero () = 0;

  virtual std::unique_ptr<NumericVector<T>> zero_clone () const = 0;

  virtual std::unique_ptr<NumericVector<T>> clone () const = 0;

  virtual void init (const numeric_index_type n,
                     const numeric_index_type n_local,
                     const bool fast = false,
                     const ParallelType ptype = AUTOMATIC) = 0;

  virtual void init (const numeric_index_type n,
                     const bool fast = false,
                     const ParallelType ptype = AUTOMATIC) = 0;

  virtual void init (const numeric_index_type n,
                     const numeric_index_type n_local,
                     const std::vector<numeric_index_type> & ghost,
                     const bool fast = false,
                     const ParallelType ptype = AUTOMATIC) = 0;

  virtual void init (const NumericVector<T> & other,
                     const bool fast = false) = 0;

  virtual NumericVector<T> & operator= (const T s) = 0;

  virtual NumericVector<T> & operator= (const std::vector<T> & v) = 0;

  virtual Real min () const = 0;

  virtual Real max () const = 0;

  virtual T sum() const = 0;

  virtual Real l1_norm () const = 0;

  virtual Real l2_norm () const = 0;

  virtual Real linfty_norm () const = 0;

  virtual Real subset_l1_norm (const std::set<numeric_index_type> & indices) const;

  virtual Real subset_l2_norm (const std::set<numeric_index_type> & indices) const;

  virtual Real subset_linfty_norm (const std::set<numeric_index_type> & indices) const;

  Real l2_norm_diff (const NumericVector<T> & other_vec) const;

  Real l1_norm_diff (const NumericVector<T> & other_vec) const;

  virtual numeric_index_type size () const = 0;

  virtual numeric_index_type local_size() const = 0;

  virtual numeric_index_type first_local_index() const = 0;

  virtual numeric_index_type last_local_index() const = 0;

  virtual T operator() (const numeric_index_type i) const = 0;

  virtual T el(const numeric_index_type i) const { return (*this)(i); }

  virtual void get(const std::vector<numeric_index_type> & index,
                   T * values) const;

  void get(const std::vector<numeric_index_type> & index,
           std::vector<T> & values) const;

  virtual NumericVector<T> & operator += (const NumericVector<T> & v) = 0;

  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) = 0;

  NumericVector<T> & operator *= (const T a) { this->scale(a); return *this; }

  virtual NumericVector<T> & operator *= (const NumericVector<T> & v) = 0;

  NumericVector<T> & operator /= (const T a) { this->scale(1./a); return *this; }

  virtual NumericVector<T> & operator /= (const NumericVector<T> & v) = 0;

  virtual void reciprocal() = 0;

  virtual void conjugate() = 0;

  virtual void set (const numeric_index_type i, const T value) = 0;

  virtual void add (const numeric_index_type i, const T value) = 0;

  virtual void add (const T s) = 0;

  virtual void add (const NumericVector<T> & v) = 0;

  virtual void add (const T a, const NumericVector<T> & v) = 0;

  virtual void add_vector (const T * v,
                           const std::vector<numeric_index_type> & dof_indices);

  void add_vector (const std::vector<T> & v,
                   const std::vector<numeric_index_type> & dof_indices);

  void add_vector (const NumericVector<T> & v,
                   const std::vector<numeric_index_type> & dof_indices);

  void add_vector (const DenseVector<T> & v,
                   const std::vector<numeric_index_type> & dof_indices);

  virtual void add_vector (const NumericVector<T> & v,
                           const SparseMatrix<T> & A) = 0;

  void add_vector (const NumericVector<T> & v,
                   const ShellMatrix<T> & A);

  virtual void add_vector_transpose (const NumericVector<T> & v,
                                     const SparseMatrix<T> & A) = 0;

  virtual void insert (const T * v,
                       const std::vector<numeric_index_type> & dof_indices);

  void insert (const std::vector<T> & v,
               const std::vector<numeric_index_type> & dof_indices);

  void insert (const NumericVector<T> & v,
               const std::vector<numeric_index_type> & dof_indices);

  void insert (const DenseVector<T> & v,
               const std::vector<numeric_index_type> & dof_indices);

  void insert (const DenseSubVector<T> & v,
               const std::vector<numeric_index_type> & dof_indices);

  virtual void scale (const T factor) = 0;

  virtual void abs() = 0;

  virtual T dot(const NumericVector<T> & v) const = 0;

  virtual void localize (std::vector<T> & v_local) const = 0;

  virtual void localize (NumericVector<T> & v_local) const = 0;

  virtual void localize (NumericVector<T> & v_local,
                         const std::vector<numeric_index_type> & send_list) const = 0;

  virtual void localize (std::vector<T> & v_local,
                         const std::vector<numeric_index_type> & indices) const = 0;

  virtual void localize (const numeric_index_type first_local_idx,
                         const numeric_index_type last_local_idx,
                         const std::vector<numeric_index_type> & send_list) = 0;

  virtual void localize_to_one (std::vector<T> & v_local,
                                const processor_id_type proc_id=0) const = 0;

  virtual int compare (const NumericVector<T> & other_vector,
                       const Real threshold = TOLERANCE) const;

  virtual int local_relative_compare (const NumericVector<T> & other_vector,
                                      const Real threshold = TOLERANCE) const;

  virtual int global_relative_compare (const NumericVector<T> & other_vector,
                                       const Real threshold = TOLERANCE) const;

  virtual void pointwise_mult (const NumericVector<T> & vec1,
                               const NumericVector<T> & vec2) = 0;

  virtual void pointwise_divide (const NumericVector<T> & vec1,
                                 const NumericVector<T> & vec2) = 0;

  virtual void print(std::ostream & os=libMesh::out) const;

  virtual void print_global(std::ostream & os=libMesh::out) const;

  friend std::ostream & operator << (std::ostream & os, const NumericVector<T> & v)
  {
    v.print_global(os);
    return os;
  }

  virtual void print_matlab(const std::string & filename = "") const;

  virtual void read_matlab(const std::string & filename);

  virtual void create_subvector(NumericVector<T> & /* subvector */,
                                const std::vector<numeric_index_type> & /* rows */,
                                bool /* supplying_global_rows */ = true) const
  {
    libmesh_not_implemented();
  }

  virtual std::unique_ptr<NumericVector<T>>
  get_subvector(const std::vector<numeric_index_type> & /* rows */)
  {
    libmesh_not_implemented();
  }

  virtual void restore_subvector(std::unique_ptr<NumericVector<T>> /* subvector */,
                                 const std::vector<numeric_index_type> & /* rows */)
  {
    libmesh_not_implemented();
  }

  virtual void swap (NumericVector<T> & v);

  virtual std::size_t max_allowed_id() const = 0;

  bool readable() const;

  bool compatible(const NumericVector<T> & v) const;

protected:

  bool _is_closed;

  bool _is_initialized;

  ParallelType _type;

  std::mutex _numeric_vector_mutex;
};


/*----------------------- Inline functions ----------------------------------*/



template <typename T>
inline
NumericVector<T>::NumericVector (const Parallel::Communicator & comm_in,
                                 const ParallelType ptype) :
  ParallelObject(comm_in),
  _is_closed(false),
  _is_initialized(false),
  _type(ptype)
{
}



template <typename T>
inline
NumericVector<T>::NumericVector (const Parallel::Communicator & comm_in,
                                 const numeric_index_type /*n*/,
                                 const ParallelType ptype) :
  ParallelObject(comm_in),
  _is_closed(false),
  _is_initialized(false),
  _type(ptype)
{
  libmesh_not_implemented(); // Abstract base class!
  // init(n, n, false, ptype);
}



template <typename T>
inline
NumericVector<T>::NumericVector (const Parallel::Communicator & comm_in,
                                 const numeric_index_type /*n*/,
                                 const numeric_index_type /*n_local*/,
                                 const ParallelType ptype) :
  ParallelObject(comm_in),
  _is_closed(false),
  _is_initialized(false),
  _type(ptype)
{
  libmesh_not_implemented(); // Abstract base class!
  // init(n, n_local, false, ptype);
}



template <typename T>
inline
NumericVector<T>::NumericVector (const Parallel::Communicator & comm_in,
                                 const numeric_index_type /*n*/,
                                 const numeric_index_type /*n_local*/,
                                 const std::vector<numeric_index_type> & /*ghost*/,
                                 const ParallelType ptype) :
  ParallelObject(comm_in),
  _is_closed(false),
  _is_initialized(false),
  _type(ptype)
{
  libmesh_not_implemented(); // Abstract base class!
  // init(n, n_local, ghost, false, ptype);
}



template <typename T>
inline
void NumericVector<T>::clear ()
{
  _is_closed      = false;
  _is_initialized = false;
}



template <typename T>
inline
void NumericVector<T>::get(const std::vector<numeric_index_type> & index,
                           T * values) const
{
  const std::size_t num = index.size();
  for (std::size_t i=0; i<num; i++)
    {
      values[i] = (*this)(index[i]);
    }
}



template <typename T>
inline
void NumericVector<T>::get(const std::vector<numeric_index_type> & index,
                           std::vector<T> & values) const
{
  const std::size_t num = index.size();
  values.resize(num);
  if (!num)
    return;

  this->get(index, values.data());
}



template <typename T>
inline
void NumericVector<T>::add_vector(const std::vector<T> & v,
                                  const std::vector<numeric_index_type> & dof_indices)
{
  libmesh_assert(v.size() == dof_indices.size());
  if (!v.empty())
    this->add_vector(v.data(), dof_indices);
}



template <typename T>
inline
void NumericVector<T>::add_vector(const DenseVector<T> & v,
                                  const std::vector<numeric_index_type> & dof_indices)
{
  libmesh_assert(v.size() == dof_indices.size());
  if (!v.empty())
    this->add_vector(&v(0), dof_indices);
}



template <typename T>
inline
void NumericVector<T>::insert(const std::vector<T> & v,
                              const std::vector<numeric_index_type> & dof_indices)
{
  libmesh_assert(v.size() == dof_indices.size());
  if (!v.empty())
    this->insert(v.data(), dof_indices);
}



template <typename T>
inline
void NumericVector<T>::insert(const DenseVector<T> & v,
                              const std::vector<numeric_index_type> & dof_indices)
{
  libmesh_assert(v.size() == dof_indices.size());
  if (!v.empty())
    this->insert(&v(0), dof_indices);
}



template <typename T>
inline
void NumericVector<T>::insert(const DenseSubVector<T> & v,
                              const std::vector<numeric_index_type> & dof_indices)
{
  libmesh_assert(v.size() == dof_indices.size());
  if (!v.empty())
    this->insert(&v(0), dof_indices);
}



// Full specialization of the print() member for complex
// variables.  This must precede the non-specialized
// version, at least according to icc v7.1
template <>
inline
void NumericVector<Complex>::print(std::ostream & os) const
{
  libmesh_assert (this->initialized());
  os << "Size\tglobal =  " << this->size()
     << "\t\tlocal =  " << this->local_size() << std::endl;

  // std::complex<>::operator<<() is defined, but use this form
  os << "#\tReal part\t\tImaginary part" << std::endl;
  for (auto i : index_range(*this))
    os << i << "\t"
       << (*this)(i).real() << "\t\t"
       << (*this)(i).imag() << std::endl;
}



template <typename T>
inline
void NumericVector<T>::print(std::ostream & os) const
{
  libmesh_assert (this->initialized());
  os << "Size\tglobal =  " << this->size()
     << "\t\tlocal =  " << this->local_size() << std::endl;

  os << "#\tValue" << std::endl;
  for (auto i : index_range(*this))
    os << i << "\t" << (*this)(i) << std::endl;
}



template <>
inline
void NumericVector<Complex>::print_global(std::ostream & os) const
{
  libmesh_assert (this->initialized());

  std::vector<Complex> v(this->size());
  this->localize(v);

  // Right now we only want one copy of the output
  if (this->processor_id())
    return;

  os << "Size\tglobal =  " << this->size() << std::endl;
  os << "#\tReal part\t\tImaginary part" << std::endl;
  for (auto i : make_range(v.size()))
    os << i << "\t"
       << v[i].real() << "\t\t"
       << v[i].imag() << std::endl;
}


template <typename T>
inline
void NumericVector<T>::print_global(std::ostream & os) const
{
  libmesh_assert (this->initialized());

  std::vector<T> v(this->size());
  this->localize(v);

  // Right now we only want one copy of the output
  if (this->processor_id())
    return;

  os << "Size\tglobal =  " << this->size() << std::endl;
  os << "#\tValue" << std::endl;
  for (auto i : make_range(v.size()))
    os << i << "\t" << v[i] << std::endl;
}



template <typename T>
inline
void  NumericVector<T>::swap (NumericVector<T> & v)
{
  std::swap(_is_closed, v._is_closed);
  std::swap(_is_initialized, v._is_initialized);
  std::swap(_type, v._type);
}

template <typename T>
auto
l1_norm(const NumericVector<T> & vec)
{
  return vec.l1_norm();
}

template <typename T>
auto
l1_norm_diff(const NumericVector<T> & vec1, const NumericVector<T> & vec2)
{
  return vec1.l1_norm_diff(vec2);
}

} // namespace libMesh


// Workaround for weird boost/NumericVector interaction bug
#ifdef LIBMESH_DEFAULT_QUADRUPLE_PRECISION
#include <boost/mpl/bool.hpp>

namespace boost { namespace multiprecision { namespace detail {
template <typename T, typename To>
struct is_lossy_conversion<libMesh::NumericVector<T>, To> {
  typedef boost::mpl::true_ type;
  static const bool value = type::value;
};
}}}
#endif


#endif  // LIBMESH_NUMERIC_VECTOR_H