bvector.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_ISTL_BVECTOR_HH
#define DUNE_ISTL_BVECTOR_HH
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <initializer_list>
#include <limits>
#include <memory>
#include <utility>
#include <vector>
 
#include <dune/common/dotproduct.hh>
#include <dune/common/ftraits.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/promotiontraits.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/scalarvectorview.hh>
 
#include <dune/istl/blocklevel.hh>
 
#include "basearray.hh"
#include "istlexception.hh"
 
namespace Dune {
 
namespace Imp {
 
  template <class B, bool isNumber>
  class BlockTraitsImp;
 
  template <class B>
  class BlockTraitsImp<B,true>
  {
  public:
    using field_type = B;
  };
 
  template <class B>
  class BlockTraitsImp<B,false>
  {
  public:
    using field_type = typename B::field_type;
  };
 
  template <class B>
  using BlockTraits = BlockTraitsImp<B,IsNumber<B>::value>;
 
  template<class B, class A=std::allocator<B> >
  class block_vector_unmanaged : public base_array_unmanaged<B,A>
  {
  public:
 
    //===== type definitions and constants
    using field_type = typename Imp::BlockTraits<B>::field_type;
 
    typedef B block_type;
 
    typedef A allocator_type;
 
    typedef typename A::size_type size_type;
 
    typedef typename base_array_unmanaged<B,A>::iterator Iterator;
 
    typedef typename base_array_unmanaged<B,A>::const_iterator ConstIterator;
 
    typedef B value_type;
 
    typedef B& reference;
 
    typedef const B& const_reference;
 
    //===== assignment from scalar
 
    block_vector_unmanaged& operator= (const field_type& k)
    {
      for (size_type i=0; i<this->n; i++)
        (*this)[i] = k;
      return *this;
    }
 
    //===== vector space arithmetic
    block_vector_unmanaged& operator+= (const block_vector_unmanaged& y)
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
#endif
      for (size_type i=0; i<this->n; ++i) (*this)[i] += y[i];
      return *this;
    }
 
    block_vector_unmanaged& operator-= (const block_vector_unmanaged& y)
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
#endif
      for (size_type i=0; i<this->n; ++i) (*this)[i] -= y[i];
      return *this;
    }
 
    block_vector_unmanaged& operator*= (const field_type& k)
    {
      for (size_type i=0; i<this->n; ++i) (*this)[i] *= k;
      return *this;
    }
 
    block_vector_unmanaged& operator/= (const field_type& k)
    {
      for (size_type i=0; i<this->n; ++i) (*this)[i] /= k;
      return *this;
    }
 
    block_vector_unmanaged& axpy (const field_type& a, const block_vector_unmanaged& y)
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
#endif
      for (size_type i=0; i<this->n; ++i)
        Impl::asVector((*this)[i]).axpy(a,Impl::asVector(y[i]));
 
      return *this;
    }
 
 
    template<class OtherB, class OtherA>
    auto operator* (const block_vector_unmanaged<OtherB,OtherA>& y) const
    {
      typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
      PromotedType sum(0);
#ifdef DUNE_ISTL_WITH_CHECKING
      if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
#endif
      for (size_type i=0; i<this->n; ++i) {
        sum += PromotedType(((*this)[i])*y[i]);
      }
      return sum;
    }
 
    template<class OtherB, class OtherA>
    auto dot(const block_vector_unmanaged<OtherB,OtherA>& y) const
    {
      typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
      PromotedType sum(0);
#ifdef DUNE_ISTL_WITH_CHECKING
      if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
#endif
 
      for (size_type i=0; i<this->n; ++i)
        sum += Impl::asVector((*this)[i]).dot(Impl::asVector(y[i]));
 
      return sum;
    }
 
    //===== norms
 
    typename FieldTraits<field_type>::real_type one_norm () const
    {
      typename FieldTraits<field_type>::real_type sum=0;
      for (size_type i=0; i<this->n; ++i)
        sum += Impl::asVector((*this)[i]).one_norm();
      return sum;
    }
 
    typename FieldTraits<field_type>::real_type one_norm_real () const
    {
      typename FieldTraits<field_type>::real_type sum=0;
      for (size_type i=0; i<this->n; ++i)
        sum += Impl::asVector((*this)[i]).one_norm_real();
      return sum;
    }
 
    typename FieldTraits<field_type>::real_type two_norm () const
    {
      using std::sqrt;
      return sqrt(two_norm2());
    }
 
    typename FieldTraits<field_type>::real_type two_norm2 () const
    {
      typename FieldTraits<field_type>::real_type sum=0;
      for (size_type i=0; i<this->n; ++i)
        sum += Impl::asVector((*this)[i]).two_norm2();
      return sum;
    }
 
    template <typename ft = field_type,
              typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
    typename FieldTraits<ft>::real_type infinity_norm() const {
      using real_type = typename FieldTraits<ft>::real_type;
      using std::max;
 
      real_type norm = 0;
      for (auto const &xi : *this) {
        real_type const a = Impl::asVector(xi).infinity_norm();
        norm = max(a, norm);
      }
      return norm;
    }
 
    template <typename ft = field_type,
              typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
    typename FieldTraits<ft>::real_type infinity_norm_real() const {
      using real_type = typename FieldTraits<ft>::real_type;
      using std::max;
 
      real_type norm = 0;
      for (auto const &xi : *this) {
        real_type const a = Impl::asVector(xi).infinity_norm_real();
        norm = max(a, norm);
      }
      return norm;
    }
 
    template <typename ft = field_type,
              typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
    typename FieldTraits<ft>::real_type infinity_norm() const {
      using real_type = typename FieldTraits<ft>::real_type;
      using std::max;
      using std::abs;
 
      real_type norm = 0;
      real_type isNaN = 1;
 
      for (auto const &xi : *this) {
        real_type const a = Impl::asVector(xi).infinity_norm();
        norm = max(a, norm);
        isNaN += a;
      }
      return norm * (isNaN / isNaN);
    }
 
    template <typename ft = field_type,
              typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
    typename FieldTraits<ft>::real_type infinity_norm_real() const {
      using real_type = typename FieldTraits<ft>::real_type;
      using std::max;
 
      real_type norm = 0;
      real_type isNaN = 1;
 
      for (auto const &xi : *this) {
        real_type const a = Impl::asVector(xi).infinity_norm_real();
        norm = max(a, norm);
        isNaN += a;
      }
 
      return norm * (isNaN / isNaN);
    }
 
    //===== sizes
 
    size_type N () const
    {
      return this->n;
    }
 
    size_type dim () const
    {
      size_type d=0;
 
      for (size_type i=0; i<this->n; i++)
        d += Impl::asVector((*this)[i]).dim();
 
      return d;
    }
 
  protected:
    block_vector_unmanaged () : base_array_unmanaged<B,A>()
    {       }
  };
 
 
  template<class F>
  class ScopeGuard {
    F cleanupFunc_;
  public:
    ScopeGuard(F cleanupFunc) : cleanupFunc_(std::move(cleanupFunc)) {}
    ScopeGuard(const ScopeGuard &) = delete;
    ScopeGuard(ScopeGuard &&) = delete;
    ScopeGuard &operator=(ScopeGuard) = delete;
    ~ScopeGuard() { cleanupFunc_(); }
  };
 
 
  template<class F>
  ScopeGuard<F> makeScopeGuard(F cleanupFunc)
  {
    return { std::move(cleanupFunc) };
  }
 
} // end namespace Imp
  template<class B, class A=std::allocator<B> >
  class BlockVector : public Imp::block_vector_unmanaged<B,A>
  {
  public:
 
    //===== type definitions and constants
 
    using field_type = typename Imp::BlockTraits<B>::field_type;
 
    typedef B block_type;
 
    typedef A allocator_type;
 
    typedef typename A::size_type size_type;
 
    [[deprecated("Use free function blockLevel(). Will be removed after 2.8.")]]
    static constexpr unsigned int blocklevel = blockLevel<B>()+1;
 
    typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
 
    typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
 
    //===== constructors and such
 
    BlockVector ()
    {
      syncBaseArray();
    }
 
    explicit BlockVector (size_type _n) : storage_(_n)
    {
      syncBaseArray();
    }
 
    BlockVector (std::initializer_list<B> const &l) : storage_(l)
    {
      syncBaseArray();
    }
 
 
    template<typename S>
    BlockVector (size_type _n, S _capacity)
    {
      static_assert(std::numeric_limits<S>::is_integer,
        "capacity must be an unsigned integral type (be aware, that this constructor does not set the default value!)" );
      if((size_type)_capacity > _n)
        storage_.reserve(_capacity);
      storage_.resize(_n);
      syncBaseArray();
    }
 
 
    void reserve(size_type capacity)
    {
      [[maybe_unused]] const auto &guard =
        Imp::makeScopeGuard([this]{ syncBaseArray(); });
      storage_.reserve(capacity);
    }
 
    size_type capacity() const
    {
      return storage_.capacity();
    }
 
    void resize(size_type size)
    {
      [[maybe_unused]] const auto &guard =
        Imp::makeScopeGuard([this]{ syncBaseArray(); });
      storage_.resize(size);
    }
 
    BlockVector(const BlockVector &a)
      noexcept(noexcept(std::declval<BlockVector>().storage_ = a.storage_))
    {
      storage_ = a.storage_;
      syncBaseArray();
    }
 
    BlockVector(BlockVector &&a)
      noexcept(noexcept(std::declval<BlockVector>().swap(a)))
    {
      swap(a);
    }
 
    BlockVector& operator= (const BlockVector& a)
      noexcept(noexcept(std::declval<BlockVector>().storage_ = a.storage_))
    {
      [[maybe_unused]] const auto &guard =
        Imp::makeScopeGuard([this]{ syncBaseArray(); });
      storage_ = a.storage_;
      return *this;
    }
 
    BlockVector& operator= (BlockVector&& a)
      noexcept(noexcept(std::declval<BlockVector>().swap(a)))
    {
      swap(a);
      return *this;
    }
 
    void swap(BlockVector &other)
      noexcept(noexcept(
            std::declval<BlockVector&>().storage_.swap(other.storage_)))
    {
      [[maybe_unused]] const auto &guard = Imp::makeScopeGuard([&]{
          syncBaseArray();
          other.syncBaseArray();
        });
      storage_.swap(other.storage_);
    }
 
    BlockVector& operator= (const field_type& k)
    {
      // forward to operator= in base class
      (static_cast<Imp::block_vector_unmanaged<B,A>&>(*this)) = k;
      return *this;
    }
 
  private:
    void syncBaseArray() noexcept
    {
      this->p = storage_.data();
      this->n = storage_.size();
    }
 
    std::vector<B, A> storage_;
  };
 
  template<class B, class A>
  struct FieldTraits< BlockVector<B, A> >
  {
    typedef typename FieldTraits<B>::field_type field_type;
    typedef typename FieldTraits<B>::real_type real_type;
  };
  template<class K, class A>
  std::ostream& operator<< (std::ostream& s, const BlockVector<K, A>& v)
  {
    typedef typename  BlockVector<K, A>::size_type size_type;
 
    for (size_type i=0; i<v.size(); i++)
      s << v[i] << std::endl;
 
    return s;
  }
 
namespace Imp {
 
#ifndef DOXYGEN
  template<class B, class A>
#else
  template<class B, class A=std::allocator<B> >
#endif
  class BlockVectorWindow : public Imp::block_vector_unmanaged<B,A>
  {
  public:
 
    //===== type definitions and constants
 
    using field_type = typename Imp::BlockTraits<B>::field_type;
 
    typedef B block_type;
 
    typedef A allocator_type;
 
    typedef typename A::size_type size_type;
 
    [[deprecated("Use free function blockLevel(). Will be removed after 2.8.")]]
    static constexpr unsigned int blocklevel = blockLevel<B>()+1;
 
    typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
 
    typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
 
 
    //===== constructors and such
    BlockVectorWindow () : Imp::block_vector_unmanaged<B,A>()
    {       }
 
    BlockVectorWindow (B* _p, size_type _n)
    {
      this->n = _n;
      this->p = _p;
    }
 
    BlockVectorWindow (const BlockVectorWindow& a)
    {
      this->n = a.n;
      this->p = a.p;
    }
 
    BlockVectorWindow& operator= (const BlockVectorWindow& a)
    {
      // check correct size
#ifdef DUNE_ISTL_WITH_CHECKING
      if (this->n!=a.N()) DUNE_THROW(ISTLError,"vector size mismatch");
#endif
 
      if (&a!=this)     // check if this and a are different objects
      {
        // copy data
        for (size_type i=0; i<this->n; i++) this->p[i]=a.p[i];
      }
      return *this;
    }
 
    BlockVectorWindow& operator= (const field_type& k)
    {
      (static_cast<Imp::block_vector_unmanaged<B,A>&>(*this)) = k;
      return *this;
    }
 
    operator BlockVector<B, A>() const {
      auto bv = BlockVector<B, A>(this->n);
 
      std::copy(this->begin(), this->end(), bv.begin());
 
      return bv;
    }
 
    //===== window manipulation methods
 
    void set (size_type _n, B* _p)
    {
      this->n = _n;
      this->p = _p;
    }
 
    void setsize (size_type _n)
    {
      this->n = _n;
    }
 
    void setptr (B* _p)
    {
      this->p = _p;
    }
 
    B* getptr ()
    {
      return this->p;
    }
 
    size_type getsize () const
    {
      return this->n;
    }
  };
 
 
 
  template<class B, class A=std::allocator<B> >
  class compressed_block_vector_unmanaged : public compressed_base_array_unmanaged<B,A>
  {
  public:
 
    //===== type definitions and constants
 
    using field_type = typename Imp::BlockTraits<B>::field_type;
 
    typedef B block_type;
 
    typedef A allocator_type;
 
    typedef typename compressed_base_array_unmanaged<B,A>::iterator Iterator;
 
    typedef typename compressed_base_array_unmanaged<B,A>::const_iterator ConstIterator;
 
    typedef typename A::size_type size_type;
 
    //===== assignment from scalar
 
    compressed_block_vector_unmanaged& operator= (const field_type& k)
    {
      for (size_type i=0; i<this->n; i++)
        (this->p)[i] = k;
      return *this;
    }
 
 
    //===== vector space arithmetic
 
    template<class V>
    compressed_block_vector_unmanaged& operator+= (const V& y)
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (!includesindexset(y)) DUNE_THROW(ISTLError,"index set mismatch");
#endif
      for (size_type i=0; i<y.n; ++i) this->operator[](y.j[i]) += y.p[i];
      return *this;
    }
 
    template<class V>
    compressed_block_vector_unmanaged& operator-= (const V& y)
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (!includesindexset(y)) DUNE_THROW(ISTLError,"index set mismatch");
#endif
      for (size_type i=0; i<y.n; ++i) this->operator[](y.j[i]) -= y.p[i];
      return *this;
    }
 
    template<class V>
    compressed_block_vector_unmanaged& axpy (const field_type& a, const V& y)
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (!includesindexset(y)) DUNE_THROW(ISTLError,"index set mismatch");
#endif
      for (size_type i=0; i<y.n; ++i)
        Impl::asVector((*this)[y.j[i]]).axpy(a,Impl::asVector(y.p[i]));
      return *this;
    }
 
    compressed_block_vector_unmanaged& operator*= (const field_type& k)
    {
      for (size_type i=0; i<this->n; ++i) (this->p)[i] *= k;
      return *this;
    }
 
    compressed_block_vector_unmanaged& operator/= (const field_type& k)
    {
      for (size_type i=0; i<this->n; ++i) (this->p)[i] /= k;
      return *this;
    }
 
 
    //===== Euclidean scalar product
 
    field_type operator* (const compressed_block_vector_unmanaged& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (!includesindexset(y) || !y.includesindexset(*this) )
        DUNE_THROW(ISTLError,"index set mismatch");
#endif
      field_type sum=0;
      for (size_type i=0; i<this->n; ++i)
        sum += (this->p)[i] * y[(this->j)[i]];
      return sum;
    }
 
 
    //===== norms
 
    typename FieldTraits<field_type>::real_type one_norm () const
    {
      typename FieldTraits<field_type>::real_type sum=0;
      for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].one_norm();
      return sum;
    }
 
    typename FieldTraits<field_type>::real_type one_norm_real () const
    {
      typename FieldTraits<field_type>::real_type sum=0;
      for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].one_norm_real();
      return sum;
    }
 
    typename FieldTraits<field_type>::real_type two_norm () const
    {
      using std::sqrt;
      typename FieldTraits<field_type>::real_type sum=0;
      for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].two_norm2();
      return sqrt(sum);
    }
 
    typename FieldTraits<field_type>::real_type two_norm2 () const
    {
      typename FieldTraits<field_type>::real_type sum=0;
      for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].two_norm2();
      return sum;
    }
 
    template <typename ft = field_type,
              typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
    typename FieldTraits<ft>::real_type infinity_norm() const {
      using real_type = typename FieldTraits<ft>::real_type;
      using std::max;
 
      real_type norm = 0;
      for (auto const &x : *this) {
        real_type const a = x.infinity_norm();
        norm = max(a, norm);
      }
      return norm;
    }
 
    template <typename ft = field_type,
              typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
    typename FieldTraits<ft>::real_type infinity_norm_real() const {
      using real_type = typename FieldTraits<ft>::real_type;
      using std::max;
 
      real_type norm = 0;
      for (auto const &x : *this) {
        real_type const a = x.infinity_norm_real();
        norm = max(a, norm);
      }
      return norm;
    }
 
    template <typename ft = field_type,
              typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
    typename FieldTraits<ft>::real_type infinity_norm() const {
      using real_type = typename FieldTraits<ft>::real_type;
      using std::max;
 
      real_type norm = 0;
      real_type isNaN = 1;
      for (auto const &x : *this) {
        real_type const a = x.infinity_norm();
        norm = max(a, norm);
        isNaN += a;
      }
      return norm * (isNaN / isNaN);
    }
 
    template <typename ft = field_type,
              typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
    typename FieldTraits<ft>::real_type infinity_norm_real() const {
      using real_type = typename FieldTraits<ft>::real_type;
      using std::max;
 
      real_type norm = 0;
      real_type isNaN = 1;
      for (auto const &x : *this) {
        real_type const a = x.infinity_norm_real();
        norm = max(a, norm);
        isNaN += a;
      }
      return norm * (isNaN / isNaN);
    }
 
    //===== sizes
 
    size_type N () const
    {
      return this->n;
    }
 
    size_type dim () const
    {
      size_type d=0;
      for (size_type i=0; i<this->n; i++)
        d += (this->p)[i].dim();
      return d;
    }
 
  protected:
    compressed_block_vector_unmanaged () : compressed_base_array_unmanaged<B,A>()
    {       }
 
    template<class V>
    bool includesindexset (const V& y)
    {
      typename V::ConstIterator e=this->end();
      for (size_type i=0; i<y.n; i++)
        if (this->find(y.j[i])==e)
          return false;
      return true;
    }
  };
 
 
  template<class B, class A=std::allocator<B> >
  class CompressedBlockVectorWindow : public compressed_block_vector_unmanaged<B,A>
  {
  public:
 
    //===== type definitions and constants
 
    using field_type = typename Imp::BlockTraits<B>::field_type;
 
    typedef B block_type;
 
    typedef A allocator_type;
 
    typedef typename A::size_type size_type;
 
    [[deprecated("Use free function blockLevel(). Will be removed after 2.8.")]]
    static constexpr unsigned int blocklevel = blockLevel<B>()+1;
 
    typedef typename compressed_block_vector_unmanaged<B,A>::Iterator Iterator;
 
    typedef typename compressed_block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
 
 
    //===== constructors and such
    CompressedBlockVectorWindow () : compressed_block_vector_unmanaged<B,A>()
    {       }
 
    CompressedBlockVectorWindow (B* _p, size_type* _j, size_type _n)
    {
      this->n = _n;
      this->p = _p;
      this->j = _j;
    }
 
    CompressedBlockVectorWindow (const CompressedBlockVectorWindow& a)
    {
      this->n = a.n;
      this->p = a.p;
      this->j = a.j;
    }
 
    CompressedBlockVectorWindow& operator= (const CompressedBlockVectorWindow& a)
    {
      // check correct size
#ifdef DUNE_ISTL_WITH_CHECKING
      if (this->n!=a.N()) DUNE_THROW(ISTLError,"vector size mismatch");
#endif
 
      if (&a!=this)     // check if this and a are different objects
      {
        // copy data
        for (size_type i=0; i<this->n; i++) this->p[i]=a.p[i];
        for (size_type i=0; i<this->n; i++) this->j[i]=a.j[i];
      }
      return *this;
    }
 
    CompressedBlockVectorWindow& operator= (const field_type& k)
    {
      (static_cast<compressed_block_vector_unmanaged<B,A>&>(*this)) = k;
      return *this;
    }
 
 
    //===== window manipulation methods
 
    void set (size_type _n, B* _p, size_type* _j)
    {
      this->n = _n;
      this->p = _p;
      this->j = _j;
    }
 
    void setsize (size_type _n)
    {
      this->n = _n;
    }
 
    void setptr (B* _p)
    {
      this->p = _p;
    }
 
    void setindexptr (size_type* _j)
    {
      this->j = _j;
    }
 
    B* getptr ()
    {
      return this->p;
    }
 
    size_type* getindexptr ()
    {
      return this->j;
    }
 
    const B* getptr () const
    {
      return this->p;
    }
 
    const size_type* getindexptr () const
    {
      return this->j;
    }
    size_type getsize () const
    {
      return this->n;
    }
  };
 
} // end namespace 'Imp'
 
 
  template<typename B, typename A>
  struct AutonomousValueType<Imp::BlockVectorWindow<B,A>>
  {
    using type = BlockVector<B, A>;
  };
 
 
} // end namespace 'Dune'
 
#endif