bvector.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- 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 ST=std::size_t >
  class block_vector_unmanaged : public base_array_unmanaged<B,ST>
  {
  public:

    //===== type definitions and constants
    using field_type = typename Imp::BlockTraits<B>::field_type;

    typedef B block_type;

    typedef ST size_type;

    typedef typename base_array_unmanaged<B,ST>::iterator Iterator;

    typedef typename base_array_unmanaged<B,ST>::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 OtherST>
    auto operator* (const block_vector_unmanaged<OtherB,OtherST>& 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 OtherST>
    auto dot(const block_vector_unmanaged<OtherB,OtherST>& 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,ST>()
    {       }
  };


  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,typename A::size_type>
  {
  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 Imp::block_vector_unmanaged<B,size_type>::Iterator Iterator;

    typedef typename Imp::block_vector_unmanaged<B,size_type>::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,size_type>&>(*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,typename A::size_type>
  {
  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 Imp::block_vector_unmanaged<B,size_type>::Iterator Iterator;

    typedef typename Imp::block_vector_unmanaged<B,size_type>::ConstIterator ConstIterator;


    //===== constructors and such
    BlockVectorWindow () : Imp::block_vector_unmanaged<B,size_type>()
    {       }

    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,size_type>&>(*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 ST=std::size_t >
  class compressed_block_vector_unmanaged : public compressed_base_array_unmanaged<B,ST>
  {
  public:

    //===== type definitions and constants

    using field_type = typename Imp::BlockTraits<B>::field_type;

    typedef B block_type;

    typedef typename compressed_base_array_unmanaged<B,ST>::iterator Iterator;

    typedef typename compressed_base_array_unmanaged<B,ST>::const_iterator ConstIterator;

    typedef ST 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,ST>()
    {       }

    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 ST=std::size_t >
  class CompressedBlockVectorWindow : public compressed_block_vector_unmanaged<B,ST>
  {
  public:

    //===== type definitions and constants

    using field_type = typename Imp::BlockTraits<B>::field_type;

    typedef B block_type;

    typedef ST size_type;

    typedef typename compressed_block_vector_unmanaged<B,ST>::Iterator Iterator;

    typedef typename compressed_block_vector_unmanaged<B,ST>::ConstIterator ConstIterator;


    //===== constructors and such
    CompressedBlockVectorWindow () : compressed_block_vector_unmanaged<B,ST>()
    {       }

    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,ST>&>(*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