vector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2016 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, 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.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__vector_h
#define dealii__vector_h


#include <deal.II/base/config.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/index_set.h>
#include <boost/serialization/array.hpp>
#include <boost/serialization/split_member.hpp>

#include <cstdio>
#include <iostream>
#include <cstring>
#include <vector>

DEAL_II_NAMESPACE_OPEN


#ifdef DEAL_II_WITH_PETSC
namespace PETScWrappers
{
  class Vector;
  namespace MPI
  {
    class Vector;
  }
}
#endif

#ifdef DEAL_II_WITH_TRILINOS
namespace TrilinosWrappers
{
  namespace MPI
  {
    class Vector;
  }
  class Vector;
}
#endif

template<typename number> class LAPACKFullMatrix;

template <typename> class BlockVector;

template <typename> class VectorView;




struct VectorOperation
{
  enum values { unknown, insert, add };
};


template <typename Number>
class Vector : public Subscriptor
{
public:
  typedef Number                                            value_type;
  typedef value_type                                       *pointer;
  typedef const value_type                                 *const_pointer;
  typedef value_type                                       *iterator;
  typedef const value_type                                 *const_iterator;
  typedef value_type                                       &reference;
  typedef const value_type                                 &const_reference;
  typedef types::global_dof_index                           size_type;

  typedef typename numbers::NumberTraits<Number>::real_type real_type;

  static const bool supports_distributed_data = false;

public:

  Vector ();

  Vector (const Vector<Number> &v);

#ifdef DEAL_II_WITH_CXX11

  Vector (Vector<Number> &&v);
#endif


#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG

  template <typename OtherNumber>
  explicit
  Vector (const Vector<OtherNumber> &v);
#endif

#ifdef DEAL_II_WITH_PETSC

  Vector (const PETScWrappers::Vector &v);

  Vector (const PETScWrappers::MPI::Vector &v);
#endif

#ifdef DEAL_II_WITH_TRILINOS

  Vector (const TrilinosWrappers::MPI::Vector &v);

  Vector (const TrilinosWrappers::Vector &v);
#endif

  explicit Vector (const size_type n);

  template <typename InputIterator>
  Vector (const InputIterator first,
          const InputIterator last);

  virtual ~Vector ();

  void compress (::VectorOperation::values operation
                 =::VectorOperation::unknown) const;

  virtual void reinit (const size_type N,
                       const bool      omit_zeroing_entries=false);

  template <typename Number2>
  void reinit (const Vector<Number2> &V,
               const bool            omit_zeroing_entries=false);

  virtual void swap (Vector<Number> &v);

  Vector<Number> &operator= (const Number s);

  Vector<Number> &operator= (const Vector<Number> &v);

#ifdef DEAL_II_WITH_CXX11

  Vector<Number> &operator= (Vector<Number> &&v);
#endif

  template <typename Number2>
  Vector<Number> &operator= (const Vector<Number2> &v);

  Vector<Number> &operator= (const BlockVector<Number> &v);

#ifdef DEAL_II_WITH_PETSC

  Vector<Number> &
  operator= (const PETScWrappers::Vector &v);

  Vector<Number> &
  operator= (const PETScWrappers::MPI::Vector &v);
#endif


#ifdef DEAL_II_WITH_TRILINOS

  Vector<Number> &
  operator= (const TrilinosWrappers::MPI::Vector &v);

  Vector<Number> &
  operator= (const TrilinosWrappers::Vector &v);
#endif

  template <typename Number2>
  bool operator== (const Vector<Number2> &v) const;

  template <typename Number2>
  bool operator != (const Vector<Number2> &v) const;




  template <typename Number2>
  Number operator * (const Vector<Number2> &V) const;

  real_type norm_sqr () const;

  Number mean_value () const;

  real_type l1_norm () const;

  real_type l2_norm () const;

  real_type lp_norm (const real_type p) const;

  real_type linfty_norm () const;

  Number add_and_dot (const Number          a,
                      const Vector<Number> &V,
                      const Vector<Number> &W);




  iterator begin ();

  const_iterator begin () const;

  iterator end ();

  const_iterator end () const;

  Number operator() (const size_type i) const;

  Number &operator() (const size_type i);

  Number operator[] (const size_type i) const;

  Number &operator[] (const size_type i);

  template <typename OtherNumber>
  void extract_subvector_to (const std::vector<size_type> &indices,
                             std::vector<OtherNumber> &values) const;

  template <typename ForwardIterator, typename OutputIterator>
  void extract_subvector_to (ForwardIterator       indices_begin,
                             const ForwardIterator indices_end,
                             OutputIterator        values_begin) const;



  Vector<Number> &operator += (const Vector<Number> &V);

  Vector<Number> &operator -= (const Vector<Number> &V);

  template <typename OtherNumber>
  void add (const std::vector<size_type>   &indices,
            const std::vector<OtherNumber>  &values);

  template <typename OtherNumber>
  void add (const std::vector<size_type> &indices,
            const Vector<OtherNumber>    &values);

  template <typename OtherNumber>
  void add (const size_type    n_elements,
            const size_type   *indices,
            const OtherNumber  *values);

  void add (const Number s);

  void add (const Vector<Number> &V) DEAL_II_DEPRECATED;


  void add (const Number a, const Vector<Number> &V,
            const Number b, const Vector<Number> &W);

  void add (const Number a, const Vector<Number> &V);

  void sadd (const Number          s,
             const Vector<Number> &V);

  void sadd (const Number          s,
             const Number          a,
             const Vector<Number> &V);

  void sadd (const Number          s,
             const Number          a,
             const Vector<Number> &V,
             const Number          b,
             const Vector<Number> &W) DEAL_II_DEPRECATED;

  void sadd (const Number          s,
             const Number          a,
             const Vector<Number> &V,
             const Number          b,
             const Vector<Number> &W,
             const Number          c,
             const Vector<Number> &X) DEAL_II_DEPRECATED;

  Vector<Number> &operator *= (const Number factor);

  Vector<Number> &operator /= (const Number factor);

  void scale (const Vector<Number> &scaling_factors);

  template <typename Number2>
  void scale (const Vector<Number2> &scaling_factors);

  void equ (const Number a, const Vector<Number> &u);

  template <typename Number2>
  void equ (const Number a, const Vector<Number2> &u);

  void equ (const Number a, const Vector<Number> &u,
            const Number b, const Vector<Number> &v) DEAL_II_DEPRECATED;

  void equ (const Number a, const Vector<Number> &u,
            const Number b, const Vector<Number> &v,
            const Number c, const Vector<Number> &w) DEAL_II_DEPRECATED;

  void ratio (const Vector<Number> &a,
              const Vector<Number> &b) DEAL_II_DEPRECATED;

  void update_ghost_values () const;


  void print (const char *format = 0) const DEAL_II_DEPRECATED;

  void print (std::ostream &out,
              const unsigned int precision  = 3,
              const bool scientific = true,
              const bool across     = true) const;

  void print (LogStream &out,
              const unsigned int width = 6,
              const bool across = true) const DEAL_II_DEPRECATED;

  void block_write (std::ostream &out) const;

  void block_read (std::istream &in);

  template <class Archive>
  void save (Archive &ar, const unsigned int version) const;

  template <class Archive>
  void load (Archive &ar, const unsigned int version);

  BOOST_SERIALIZATION_SPLIT_MEMBER()


  bool in_local_range (const size_type global_index) const;

  IndexSet locally_owned_elements () const;

  std::size_t size () const;

  bool all_zero () const;

  bool is_non_negative () const;

  std::size_t memory_consumption () const;

protected:

  size_type vec_size;

  size_type max_vec_size;

  Number *val;

  template <typename Number2> friend class Vector;

  template <typename Number2> friend class LAPACKFullMatrix;

  friend class VectorView<Number>;

private:

  void allocate();

  void deallocate();
};

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


#ifndef DOXYGEN


template <typename Number>
inline
Vector<Number>::Vector ()
  :
  vec_size(0),
  max_vec_size(0),
  val(0)
{}



template <typename Number>
template <typename InputIterator>
Vector<Number>::Vector (const InputIterator first, const InputIterator last)
  :
  vec_size (0),
  max_vec_size (0),
  val (0)
{
  // allocate memory. do not initialize it, as we will copy over to it in a
  // second
  reinit (std::distance (first, last), true);
  std::copy (first, last, begin());
}



template <typename Number>
inline
Vector<Number>::Vector (const size_type n)
  :
  vec_size(0),
  max_vec_size(0),
  val(0)
{
  reinit (n, false);
}



template <typename Number>
inline
Vector<Number>::~Vector ()
{
  if (val)
    {
      deallocate();
      val=0;
    }
}



template <typename Number>
inline
void Vector<Number>::reinit (const size_type n,
                             const bool omit_zeroing_entries)
{
  if (n==0)
    {
      if (val) deallocate();
      val = 0;
      max_vec_size = vec_size = 0;
      return;
    };

  if (n>max_vec_size)
    {
      if (val) deallocate();
      max_vec_size = n;
      allocate();
    };
  vec_size = n;
  if (omit_zeroing_entries == false)
    *this = static_cast<Number>(0);
}



// declare function that is implemented in vector.templates.h
namespace internal
{
  namespace Vector
  {
    template <typename T, typename U>
    void copy_vector (const ::Vector<T> &src,
                      ::Vector<U>       &dst);
  }
}



template <typename Number>
inline
Vector<Number> &
Vector<Number>::operator= (const Vector<Number> &v)
{
  ::internal::Vector::copy_vector (v, *this);
  return *this;
}



#ifdef DEAL_II_WITH_CXX11
template <typename Number>
inline
Vector<Number> &
Vector<Number>::operator= (Vector<Number> &&v)
{
  swap(v);

  return *this;
}
#endif



template <typename Number>
template <typename Number2>
inline
Vector<Number> &
Vector<Number>::operator= (const Vector<Number2> &v)
{
  internal::Vector::copy_vector (v, *this);
  return *this;
}



template <typename Number>
inline
std::size_t Vector<Number>::size () const
{
  return vec_size;
}


template <typename Number>
inline
bool Vector<Number>::in_local_range
(const size_type) const
{
  return true;
}



template <typename Number>
inline
typename Vector<Number>::iterator
Vector<Number>::begin ()
{
  return &val[0];
}



template <typename Number>
inline
typename Vector<Number>::const_iterator
Vector<Number>::begin () const
{
  return &val[0];
}



template <typename Number>
inline
typename Vector<Number>::iterator
Vector<Number>::end ()
{
  return &val[vec_size];
}



template <typename Number>
inline
typename Vector<Number>::const_iterator
Vector<Number>::end () const
{
  return &val[vec_size];
}



template <typename Number>
inline
Number Vector<Number>::operator() (const size_type i) const
{
  Assert (i<vec_size, ExcIndexRange(i,0,vec_size));
  return val[i];
}



template <typename Number>
inline
Number &Vector<Number>::operator() (const size_type i)
{
  Assert (i<vec_size, ExcIndexRangeType<size_type>(i,0,vec_size));
  return val[i];
}



template <typename Number>
inline
Number Vector<Number>::operator[] (const size_type i) const
{
  return operator()(i);
}



template <typename Number>
inline
Number &Vector<Number>::operator[] (const size_type i)
{
  return operator()(i);
}



template <typename Number>
template <typename OtherNumber>
inline
void Vector<Number>::extract_subvector_to (const std::vector<size_type> &indices,
                                           std::vector<OtherNumber> &values) const
{
  for (size_type i = 0; i < indices.size(); ++i)
    values[i] = operator()(indices[i]);
}



template <typename Number>
template <typename ForwardIterator, typename OutputIterator>
inline
void Vector<Number>::extract_subvector_to (ForwardIterator          indices_begin,
                                           const ForwardIterator    indices_end,
                                           OutputIterator           values_begin) const
{
  while (indices_begin != indices_end)
    {
      *values_begin = operator()(*indices_begin);
      indices_begin++;
      values_begin++;
    }
}



template <typename Number>
inline
Vector<Number> &
Vector<Number>::operator /= (const Number factor)
{
  AssertIsFinite(factor);
  Assert (factor != Number(0.), ExcZero() );

  this->operator *= (Number(1.)/factor);
  return *this;
}



template <typename Number>
template <typename OtherNumber>
inline
void
Vector<Number>::add (const std::vector<size_type> &indices,
                     const std::vector<OtherNumber>  &values)
{
  Assert (indices.size() == values.size(),
          ExcDimensionMismatch(indices.size(), values.size()));
  add (indices.size(), &indices[0], &values[0]);
}



template <typename Number>
template <typename OtherNumber>
inline
void
Vector<Number>::add (const std::vector<size_type> &indices,
                     const Vector<OtherNumber>    &values)
{
  Assert (indices.size() == values.size(),
          ExcDimensionMismatch(indices.size(), values.size()));
  add (indices.size(), &indices[0], values.val);
}



template <typename Number>
template <typename OtherNumber>
inline
void
Vector<Number>::add (const size_type  n_indices,
                     const size_type *indices,
                     const OtherNumber  *values)
{
  for (size_type i=0; i<n_indices; ++i)
    {
      Assert (indices[i] < vec_size, ExcIndexRange(indices[i],0,vec_size));
      Assert (numbers::is_finite(values[i]),
              ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));

      val[indices[i]] += values[i];
    }
}



template <typename Number>
template <typename Number2>
inline
bool
Vector<Number>::operator != (const Vector<Number2> &v) const
{
  return ! (*this == v);
}



template <typename Number>
inline
void
Vector<Number>::compress (::VectorOperation::values) const
{}



template <typename Number>
inline
void
Vector<Number>::update_ghost_values () const
{}



// Moved from vector.templates.h as an inline function by Luca Heltai
// on 2009/04/12 to prevent strange compiling errors, after making
// swap virtual.
template <typename Number>
inline
void
Vector<Number>::swap (Vector<Number> &v)
{
  std::swap (vec_size,     v.vec_size);
  std::swap (max_vec_size, v.max_vec_size);
  std::swap (val,          v.val);
}



template <typename Number>
template <class Archive>
inline
void
Vector<Number>::save (Archive &ar, const unsigned int) const
{
  // forward to serialization function in the base class.
  ar &static_cast<const Subscriptor &>(*this);

  ar &vec_size &max_vec_size ;
  ar &boost::serialization::make_array(val, max_vec_size);
}



template <typename Number>
template <class Archive>
inline
void
Vector<Number>::load (Archive &ar, const unsigned int)
{
  // get rid of previous content
  deallocate();

  // the load stuff again from the archive
  ar &static_cast<Subscriptor &>(*this);
  ar &vec_size &max_vec_size ;

  allocate();
  ar &boost::serialization::make_array(val, max_vec_size);
}

#endif


template <typename Number>
inline
void swap (Vector<Number> &u, Vector<Number> &v)
{
  u.swap (v);
}


template <typename number>
inline
std::ostream &
operator << (std::ostream &os, const Vector<number> &v)
{
  v.print(os);
  return os;
}

template <typename number>
inline
LogStream &
operator << (LogStream &os, const Vector<number> &v)
{
  v.print(os);
  return os;
}


DEAL_II_NAMESPACE_CLOSE

#endif