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 <cmath>
#include <complex>
#include <memory>
#include <limits>
 
#include <dune/common/promotiontraits.hh>
#include <dune/common/dotproduct.hh>
#include <dune/common/ftraits.hh>
 
#include "istlexception.hh"
#include "basearray.hh"
 
namespace Dune {
 
  template<class B, class A=std::allocator<B> >
  class BlockVectorWindow;
 
namespace Imp {
 
  template<class B, class A=std::allocator<B> >
  class block_vector_unmanaged : public base_array_unmanaged<B,A>
  {
  public:
 
    //===== type definitions and constants
 
    typedef typename B::field_type 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) (*this)[i].axpy(a,y[i]);
      return *this;
    }
 
 
    template<class OtherB, class OtherA>
    typename PromotionTraits<field_type,typename OtherB::field_type>::PromotedType operator* (const block_vector_unmanaged<OtherB,OtherA>& y) const
    {
      typedef typename PromotionTraits<field_type,typename 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>
    typename PromotionTraits<field_type,typename OtherB::field_type>::PromotedType dot(const block_vector_unmanaged<OtherB,OtherA>& y) const
    {
      typedef typename PromotionTraits<field_type,typename 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 += ((*this)[i]).dot(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 += (*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 += (*this)[i].one_norm_real();
      return sum;
    }
 
    typename FieldTraits<field_type>::real_type two_norm () const
    {
      typename FieldTraits<field_type>::real_type sum=0;
      for (size_type i=0; i<this->n; ++i) sum += (*this)[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)[i].two_norm2();
      return sum;
    }
 
    template <typename ft = field_type,
              typename std::enable_if<!has_nan<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<!has_nan<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<has_nan<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;
      }
      isNaN /= isNaN;
      return norm * isNaN;
    }
 
    template <typename ft = field_type,
              typename std::enable_if<has_nan<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;
      }
      isNaN /= isNaN;
      return norm * 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)[i].dim();
      return d;
    }
 
  protected:
    block_vector_unmanaged () : base_array_unmanaged<B,A>()
    {       }
  };
 
} // 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
 
    typedef typename B::field_type field_type;
 
    typedef B block_type;
 
    typedef A allocator_type;
 
    typedef typename A::size_type size_type;
 
    enum {
      blocklevel = B::blocklevel+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 () : Imp::block_vector_unmanaged<B,A>(),
                     capacity_(0)
    {}
 
    explicit BlockVector (size_type _n)
    {
      this->n = _n;
      capacity_ = _n;
      if (capacity_>0) {
        this->p = this->allocator_.allocate(capacity_);
        // actually construct the objects
        new(this->p)B[capacity_];
      } else
      {
        this->p = 0;
        this->n = 0;
        capacity_ = 0;
      }
    }
 
    BlockVector (std::initializer_list<B> const &l)
    {
      this->n = l.size();
      capacity_ = l.size();
      if (capacity_>0) {
        this->p = this->allocator_.allocate(capacity_);
        // actually construct the objects
        new(this->p)B[capacity_];
 
        std::copy_n(l.begin(), l.size(), this->p);
      } else
      {
        this->p = 0;
        this->n = 0;
        capacity_ = 0;
      }
    }
 
 
    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!)" );
      size_type capacity = _capacity;
      this->n = _n;
      if(this->n > capacity)
        capacity_ = _n;
      else
        capacity_ = capacity;
 
      if (capacity_>0) {
        this->p = this->allocator_.allocate(capacity_);
        new (this->p)B[capacity_];
      } else
      {
        this->p = 0;
        this->n = 0;
        capacity_ = 0;
      }
    }
 
 
    void reserve(size_type capacity, bool copyOldValues=true)
    {
      if(capacity >= Imp::block_vector_unmanaged<B,A>::N() && capacity != capacity_) {
        // save the old data
        B* pold = this->p;
 
        if(capacity>0) {
          // create new array with capacity
          this->p = this->allocator_.allocate(capacity);
          new (this->p)B[capacity];
 
          if(copyOldValues) {
            // copy the old values
            B* to = this->p;
            B* from = pold;
 
            for(size_type i=0; i < Imp::block_vector_unmanaged<B,A>::N(); ++i, ++from, ++to)
              *to = *from;
          }
          if(capacity_ > 0) {
            // Destruct old objects and free memory
            int i=capacity_;
            while (i)
              pold[--i].~B();
            this->allocator_.deallocate(pold,capacity_);
          }
        }else{
          if(capacity_ > 0)
            // free old data
            this->p = 0;
          capacity_ = 0;
        }
 
        capacity_ = capacity;
      }
    }
 
    size_type capacity() const
    {
      return capacity_;
    }
 
    void resize(size_type size, bool copyOldValues=true)
    {
      if (size > Imp::block_vector_unmanaged<B,A>::N())
        if(capacity_ < size)
          this->reserve(size, copyOldValues);
      this->n = size;
    }
 
 
 
 
    BlockVector (const BlockVector& a) :
      Imp::block_vector_unmanaged<B,A>(a)
    {
      // allocate memory with same size as a
      this->n = a.n;
      capacity_ = a.capacity_;
 
      if (capacity_>0) {
        this->p = this->allocator_.allocate(capacity_);
        new (this->p)B[capacity_];
      } else
      {
        this->n = 0;
        this->p = 0;
      }
 
      // and copy elements
      for (size_type i=0; i<this->n; i++) this->p[i]=a.p[i];
    }
 
    ~BlockVector ()
    {
      if (capacity_>0) {
        int i=capacity_;
        while (i)
          this->p[--i].~B();
        this->allocator_.deallocate(this->p,capacity_);
      }
    }
 
    BlockVector& operator= (const BlockVector& a)
    {
      if (&a!=this)     // check if this and a are different objects
      {
        // adjust size of vector
        if (capacity_!=a.capacity_)           // check if size is different
        {
          if (capacity_>0) {
            int i=capacity_;
            while (i)
              this->p[--i].~B();
            this->allocator_.deallocate(this->p,capacity_);                     // free old memory
          }
          capacity_ = a.capacity_;
          if (capacity_>0) {
            this->p = this->allocator_.allocate(capacity_);
            new (this->p)B[capacity_];
          } else
          {
            this->p = 0;
            capacity_ = 0;
          }
        }
        this->n = a.n;
        // copy data
        for (size_type i=0; i<this->n; i++)
          this->p[i]=a.p[i];
      }
      return *this;
    }
 
    BlockVector& operator= (const field_type& k)
    {
      // forward to operator= in base class
      (static_cast<Imp::block_vector_unmanaged<B,A>&>(*this)) = k;
      return *this;
    }
 
    template<class OtherAlloc>
    BlockVector& operator= (const BlockVectorWindow<B,OtherAlloc>& other)
    {
      resize(other.size());
      for(std::size_t i=0; i<other.size(); ++i)
        (*this)[i] = other[i];
      return *this;
    }
 
  protected:
    size_type capacity_;
 
    A allocator_;
 
  };
 
  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
 
    typedef typename B::field_type field_type;
 
    typedef B block_type;
 
    typedef A allocator_type;
 
    typedef typename A::size_type size_type;
 
    enum {
      blocklevel = B::blocklevel+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;
    }
 
 
    //===== 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 ()
    {
      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
 
    typedef typename B::field_type 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) (this->operator[](y.j[i])).axpy(a,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
    {
      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<!has_nan<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<!has_nan<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<has_nan<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;
      }
      isNaN /= isNaN;
      return norm * isNaN;
    }
 
    template <typename ft = field_type,
              typename std::enable_if<has_nan<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;
      }
      isNaN /= isNaN;
      return norm * 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
 
    typedef typename B::field_type field_type;
 
    typedef B block_type;
 
    typedef A allocator_type;
 
    typedef typename A::size_type size_type;
 
    enum {
      blocklevel = B::blocklevel+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'
 
} // end namespace 'Dune'
 
#endif