bcrsmatrix.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_BCRSMATRIX_HH
#define DUNE_ISTL_BCRSMATRIX_HH

#include <cmath>
#include <complex>
#include <set>
#include <iostream>
#include <algorithm>
#include <numeric>
#include <vector>
#include <map>
#include <memory>

#include "istlexception.hh"
#include "bvector.hh"
#include "matrixutils.hh"
#include "matrixindexset.hh"
#include <dune/common/stdstreams.hh>
#include <dune/common/iteratorfacades.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/ftraits.hh>
#include <dune/common/scalarvectorview.hh>
#include <dune/common/scalarmatrixview.hh>

#include <dune/istl/blocklevel.hh>

namespace Dune {

  template<typename M>
  struct MatrixDimension;


  template<typename size_type>
  struct CompressionStatistics
  {
    double avg;
    size_type maximum;
    size_type overflow_total;

    double mem_ratio;
  };


  template<class M_>
  class ImplicitMatrixBuilder
  {

  public:

    typedef M_ Matrix;

    typedef typename Matrix::block_type block_type;

    typedef typename Matrix::size_type size_type;


    class row_object
    {

    public:

      block_type& operator[](size_type j) const
      {
        return _m.entry(_i,j);
      }

#ifndef DOXYGEN

      row_object(Matrix& m, size_type i)
        : _m(m)
        , _i(i)
      {}

#endif

    private:

      Matrix& _m;
      size_type _i;

    };


    ImplicitMatrixBuilder(Matrix& m)
      : _m(m)
    {
      if (m.buildMode() != Matrix::implicit)
        DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix in implicit build mode");
      if (m.buildStage() != Matrix::building)
        DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix with set size that has not been compressed() yet");
    }


    ImplicitMatrixBuilder(Matrix& m, size_type rows, size_type cols, size_type avg_cols_per_row, double overflow_fraction)
      : _m(m)
    {
      if (m.buildStage() != Matrix::notAllocated)
        DUNE_THROW(BCRSMatrixError,"You can only set up a matrix for this ImplicitBuilder if it has no memory allocated yet");
      m.setBuildMode(Matrix::implicit);
      m.setImplicitBuildModeParameters(avg_cols_per_row,overflow_fraction);
      m.setSize(rows,cols);
    }

    row_object operator[](size_type i) const
    {
      return row_object(_m,i);
    }

    size_type N() const
    {
      return _m.N();
    }

    size_type M() const
    {
      return _m.M();
    }

  private:

    Matrix& _m;

  };

  template<class B, class A=std::allocator<B> >
  class BCRSMatrix
  {
    friend struct MatrixDimension<BCRSMatrix>;
  public:
    enum BuildStage {
      notbuilt=0,
      notAllocated=0,
      building=1,
      rowSizesBuilt=2,
      built=3
    };

    //===== 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 Imp::CompressedBlockVectorWindow<B,size_type> row_type;

    typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;

    enum BuildMode {
      row_wise,
      random,
      implicit,
      unknown
    };

    //===== random access interface to rows of the matrix

    row_type& operator[] (size_type i)
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (build_mode == implicit && ready != built)
        DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
      if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
      if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
      return r[i];
    }

    const row_type& operator[] (size_type i) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (build_mode == implicit && ready != built)
        DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
      if (built!=ready) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
      if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
      return r[i];
    }


    //===== iterator interface to rows of the matrix

    template<class T>
    class RealRowIterator
      : public RandomAccessIteratorFacade<RealRowIterator<T>, T>
    {

    public:
      typedef typename std::remove_const<T>::type ValueType;

      friend class RandomAccessIteratorFacade<RealRowIterator<const ValueType>, const ValueType>;
      friend class RandomAccessIteratorFacade<RealRowIterator<ValueType>, ValueType>;
      friend class RealRowIterator<const ValueType>;
      friend class RealRowIterator<ValueType>;

      RealRowIterator (row_type* _p, size_type _i)
        : p(_p), i(_i)
      {}

      RealRowIterator ()
        : p(0), i(0)
      {}

      RealRowIterator(const RealRowIterator<ValueType>& it)
        : p(it.p), i(it.i)
      {}


      size_type index () const
      {
        return i;
      }

      std::ptrdiff_t distanceTo(const RealRowIterator<ValueType>& other) const
      {
        assert(other.p==p);
        return (other.i-i);
      }

      std::ptrdiff_t distanceTo(const RealRowIterator<const ValueType>& other) const
      {
        assert(other.p==p);
        return (other.i-i);
      }

      bool equals (const RealRowIterator<ValueType>& other) const
      {
        assert(other.p==p);
        return i==other.i;
      }

      bool equals (const RealRowIterator<const ValueType>& other) const
      {
        assert(other.p==p);
        return i==other.i;
      }

    private:
      void increment()
      {
        ++i;
      }

      void decrement()
      {
        --i;
      }

      void advance(std::ptrdiff_t diff)
      {
        i+=diff;
      }

      T& elementAt(std::ptrdiff_t diff) const
      {
        return p[i+diff];
      }

      row_type& dereference () const
      {
        return p[i];
      }

      row_type* p;
      size_type i;
    };

    typedef RealRowIterator<row_type> iterator;
    typedef RealRowIterator<row_type> Iterator;

    Iterator begin ()
    {
      return Iterator(r,0);
    }

    Iterator end ()
    {
      return Iterator(r,n);
    }

    Iterator beforeEnd ()
    {
      return Iterator(r,n-1);
    }

    Iterator beforeBegin ()
    {
      return Iterator(r,-1);
    }

    typedef Iterator RowIterator;

    typedef typename row_type::Iterator ColIterator;

    typedef RealRowIterator<const row_type> const_iterator;
    typedef RealRowIterator<const row_type> ConstIterator;


    ConstIterator begin () const
    {
      return ConstIterator(r,0);
    }

    ConstIterator end () const
    {
      return ConstIterator(r,n);
    }

    ConstIterator beforeEnd() const
    {
      return ConstIterator(r,n-1);
    }

    ConstIterator beforeBegin () const
    {
      return ConstIterator(r,-1);
    }

    typedef ConstIterator ConstRowIterator;

    typedef typename row_type::ConstIterator ConstColIterator;

    //===== constructors & resizers

    // we use a negative compressionBufferSize to indicate that the implicit
    // mode parameters have not been set yet

    BCRSMatrix ()
      : build_mode(unknown), ready(notAllocated), n(0), m(0), nnz_(0),
        allocationSize_(0), r(0), a(0),
        avg(0), compressionBufferSize_(-1.0)
    {}

    BCRSMatrix (size_type _n, size_type _m, size_type _nnz, BuildMode bm)
      : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
        allocationSize_(0), r(0), a(0),
        avg(0), compressionBufferSize_(-1.0)
    {
      allocate(_n, _m, _nnz,true,false);
    }

    BCRSMatrix (size_type _n, size_type _m, BuildMode bm)
      : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
        allocationSize_(0), r(0), a(0),
        avg(0), compressionBufferSize_(-1.0)
    {
      allocate(_n, _m,0,true,false);
    }


    BCRSMatrix (size_type _n, size_type _m, size_type _avg, double compressionBufferSize, BuildMode bm)
      : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
        allocationSize_(0), r(0), a(0),
        avg(_avg), compressionBufferSize_(compressionBufferSize)
    {
      if (bm != implicit)
        DUNE_THROW(BCRSMatrixError,"Only call this constructor when using the implicit build mode");
      // Prevent user from setting a negative compression buffer size:
      // 1) It doesn't make sense
      // 2) We use a negative value to indicate that the parameters
      //    have not been set yet
      if (compressionBufferSize_ < 0.0)
        DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
      implicit_allocate(_n,_m);
    }

    BCRSMatrix (const BCRSMatrix& Mat)
      : build_mode(Mat.build_mode), ready(notAllocated), n(0), m(0), nnz_(0),
        allocationSize_(0), r(0), a(0),
        avg(Mat.avg), compressionBufferSize_(Mat.compressionBufferSize_)
    {
      if (!(Mat.ready == notAllocated || Mat.ready == built))
        DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copy-constructed when source matrix is completely empty (size not set) or fully built)");

      // deep copy in global array
      size_type _nnz = Mat.nonzeroes();

      j_ = Mat.j_; // enable column index sharing, release array in case of row-wise allocation
      allocate(Mat.n, Mat.m, _nnz, true, true);

      // build window structure
      copyWindowStructure(Mat);
    }

    ~BCRSMatrix ()
    {
      deallocate();
    }

    void setBuildMode(BuildMode bm)
    {
      if (ready == notAllocated)
        {
          build_mode = bm;
          return;
        }
      if (ready == building && (build_mode == unknown || build_mode == random || build_mode == row_wise) && (bm == row_wise || bm == random))
        build_mode = bm;
      else
        DUNE_THROW(InvalidStateException, "Matrix structure cannot be changed at this stage anymore (ready == "<<ready<<").");
    }

    void setSize(size_type rows, size_type columns, size_type nnz=0)
    {
      // deallocate already setup memory
      deallocate();

      if (build_mode == implicit)
      {
        if (nnz>0)
          DUNE_THROW(Dune::BCRSMatrixError,"number of non-zeroes may not be set in implicit mode, use setImplicitBuildModeParameters() instead");

        // implicit allocates differently
        implicit_allocate(rows,columns);
      }
      else
      {
        // allocate matrix memory
        allocate(rows, columns, nnz, true, false);
      }
    }

    void setImplicitBuildModeParameters(size_type _avg, double compressionBufferSize)
    {
      // Prevent user from setting a negative compression buffer size:
      // 1) It doesn't make sense
      // 2) We use a negative value to indicate that the parameters
      //    have not been set yet
      if (compressionBufferSize < 0.0)
        DUNE_THROW(BCRSMatrixError,"You cannot set a negative compressionBufferSize value");

      // make sure the parameters aren't changed after memory has been allocated
      if (ready != notAllocated)
        DUNE_THROW(InvalidStateException,"You cannot modify build mode parameters at this stage anymore");
      avg = _avg;
      compressionBufferSize_ = compressionBufferSize;
    }

    BCRSMatrix& operator= (const BCRSMatrix& Mat)
    {
      // return immediately when self-assignment
      if (&Mat==this) return *this;

      if (!((ready == notAllocated || ready == built) && (Mat.ready == notAllocated || Mat.ready == built)))
        DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copied when both target and source are empty or fully built)");

      // make it simple: ALWAYS throw away memory for a and j_
      // and deallocate rows only if n != Mat.n
      deallocate(n!=Mat.n);

      // reallocate the rows if required
      if (n>0 && n!=Mat.n) {
        // free rows
        for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
          std::allocator_traits<decltype(rowAllocator_)>::destroy(rowAllocator_, riter);
        rowAllocator_.deallocate(r,n);
      }

      nnz_ = Mat.nonzeroes();

      // allocate a, share j_
      j_ = Mat.j_;
      allocate(Mat.n, Mat.m, nnz_, n!=Mat.n, true);

      // build window structure
      copyWindowStructure(Mat);
      return *this;
    }

    BCRSMatrix& operator= (const field_type& k)
    {

      if (!(ready == notAllocated || ready == built))
        DUNE_THROW(InvalidStateException,"Scalar assignment only works on fully built BCRSMatrix)");

      for (size_type i=0; i<n; i++) r[i] = k;
      return *this;
    }

    //===== row-wise creation interface

    class CreateIterator
    {
    public:
      CreateIterator (BCRSMatrix& _Mat, size_type _i)
        : Mat(_Mat), i(_i), nnz(0), current_row(nullptr, Mat.j_.get(), 0)
      {
        if (Mat.build_mode == unknown && Mat.ready == building)
          {
            Mat.build_mode = row_wise;
          }
        if (i==0 && Mat.ready != building)
          DUNE_THROW(BCRSMatrixError,"creation only allowed for uninitialized matrix");
        if(Mat.build_mode!=row_wise)
          DUNE_THROW(BCRSMatrixError,"creation only allowed if row wise allocation was requested in the constructor");
        if(i==0 && _Mat.N()==0)
          // empty Matrix is always built.
           Mat.ready = built;
      }

      CreateIterator& operator++()
      {
        // this should only be called if matrix is in creation
        if (Mat.ready != building)
          DUNE_THROW(BCRSMatrixError,"matrix already built up");

        // row i is defined through the pattern
        // get memory for the row and initialize the j_ array
        // this depends on the allocation mode

        // compute size of the row
        size_type s = pattern.size();

        if(s>0) {
          // update number of nonzeroes including this row
          nnz += s;

          // alloc memory / set window
          if (Mat.nnz_ > 0)
          {
            // memory is allocated in one long array

            // check if that memory is sufficient
            if (nnz > Mat.nnz_)
              DUNE_THROW(BCRSMatrixError,"allocated nnz too small");

            // set row i
            Mat.r[i].set(s,nullptr,current_row.getindexptr());
            current_row.setindexptr(current_row.getindexptr()+s);
          }else{
            // memory is allocated individually per row
            // allocate and set row i
            B* b = Mat.allocator_.allocate(s);
            // use placement new to call constructor that allocates
            // additional memory.
            new (b) B[s];
            size_type* j = Mat.sizeAllocator_.allocate(s);
            Mat.r[i].set(s,b,j);
          }
        }else
          // setup empty row
          Mat.r[i].set(0,nullptr,nullptr);

        // initialize the j array for row i from pattern
        std::visit([&](const auto& row_indices){
          std::copy(row_indices.cbegin(), row_indices.cend(), Mat.r[i].getindexptr());
        }, pattern.storage());

        // now go to next row
        i++;
        pattern.clear();

        // check if this was last row
        if (i==Mat.n)
        {
          Mat.ready = built;
          if(Mat.nnz_ > 0)
          {
            // Set nnz to the exact number of nonzero blocks inserted
            // as some methods rely on it
            Mat.nnz_ = nnz;
            // allocate data array
            Mat.allocateData();
            Mat.setDataPointers();
          }
        }
        // done
        return *this;
      }

      bool operator!= (const CreateIterator& it) const
      {
        return (i!=it.i) || (&Mat!=&it.Mat);
      }

      bool operator== (const CreateIterator& it) const
      {
        return (i==it.i) && (&Mat==&it.Mat);
      }

      size_type index () const
      {
        return i;
      }

      void insert (size_type j)
      {
        pattern.insert(j);
      }

      bool contains (size_type j)
      {
        return pattern.contains(j);
      }
      size_type size() const
      {
        return pattern.size();
      }

    private:
      BCRSMatrix& Mat;     // the matrix we are defining
      size_type i;               // current row to be defined
      size_type nnz;             // count total number of nonzeros
      using PatternType = typename Impl::RowIndexSet;
      PatternType pattern;      // used to compile entries in a row
      row_type current_row;     // row pointing to the current row to setup
    };

    friend class CreateIterator;

    CreateIterator createbegin ()
    {
      return CreateIterator(*this,0);
    }

    CreateIterator createend ()
    {
      return CreateIterator(*this,n);
    }


    //===== random creation interface

    void setrowsize (size_type i, size_type s)
    {
      if (build_mode!=random)
        DUNE_THROW(BCRSMatrixError,"requires random build mode");
      if (ready != building)
        DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");

      r[i].setsize(s);
    }

    size_type getrowsize (size_type i) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
      if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
      return r[i].getsize();
    }

    void incrementrowsize (size_type i, size_type s = 1)
    {
      if (build_mode!=random)
        DUNE_THROW(BCRSMatrixError,"requires random build mode");
      if (ready != building)
        DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");

      r[i].setsize(r[i].getsize()+s);
    }

    void endrowsizes ()
    {
      if (build_mode!=random)
        DUNE_THROW(BCRSMatrixError,"requires random build mode");
      if (ready != building)
        DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");

      // compute total size, check positivity
      size_type total=0;
      for (size_type i=0; i<n; i++)
      {
        total += r[i].getsize();
      }

      if(nnz_ == 0)
        // allocate/check memory
        allocate(n,m,total,false,false);
      else if(nnz_ < total)
        DUNE_THROW(BCRSMatrixError,"Specified number of nonzeros ("<<nnz_<<") not "
                                                             <<"sufficient for calculated nonzeros ("<<total<<"! ");

      // set the window pointers correctly
      setColumnPointers(begin());

      // initialize j_ array with m (an invalid column index)
      // this indicates an unused entry
      for (size_type k=0; k<nnz_; k++)
        j_.get()[k] = m;
      ready = rowSizesBuilt;
    }


    void addindex (size_type row, size_type col)
    {
      if (build_mode!=random)
        DUNE_THROW(BCRSMatrixError,"requires random build mode");
      if (ready==built)
        DUNE_THROW(BCRSMatrixError,"matrix already built up");
      if (ready==building)
        DUNE_THROW(BCRSMatrixError,"matrix row sizes not built up yet");
      if (ready==notAllocated)
        DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");

      if (col >= m)
        DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");

      // get row range
      size_type* const first = r[row].getindexptr();
      size_type* const last = first + r[row].getsize();

      // find correct insertion position for new column index
      size_type* pos = std::lower_bound(first,last,col);

      // check if index is already in row
      if (pos!=last && *pos == col) return;

      // find end of already inserted column indices
      size_type* end = std::lower_bound(pos,last,m);
      if (end==last)
        DUNE_THROW(BCRSMatrixError,"row is too small");

      // insert new column index at correct position
      std::copy_backward(pos,end,end+1);
      *pos = col;
    }


    template<typename It>
    void setIndicesNoSort(size_type row, It begin, It end)
    {
      size_type row_size = r[row].size();
      size_type* col_begin = r[row].getindexptr();
      size_type* col_end;
      // consistency check between allocated row size and number of passed column indices
      if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
        DUNE_THROW(BCRSMatrixError,"Given size of row " << row
                   << " (" << row_size
                   << ") does not match number of passed entries (" << (col_end - col_begin) << ")");
    }



    template<typename It>
    void setIndices(size_type row, It begin, It end)
    {
      size_type row_size = r[row].size();
      size_type* col_begin = r[row].getindexptr();
      size_type* col_end;
      // consistency check between allocated row size and number of passed column indices
      if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
        DUNE_THROW(BCRSMatrixError,"Given size of row " << row
                   << " (" << row_size
                   << ") does not match number of passed entries (" << (col_end - col_begin) << ")");
      std::sort(col_begin,col_end);
    }

    void endindices ()
    {
      if (build_mode!=random)
        DUNE_THROW(BCRSMatrixError,"requires random build mode");
      if (ready==built)
        DUNE_THROW(BCRSMatrixError,"matrix already built up");
      if (ready==building)
        DUNE_THROW(BCRSMatrixError,"row sizes are not built up yet");
      if (ready==notAllocated)
        DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");

      // check if there are undefined indices
      RowIterator endi=end();
      for (RowIterator i=begin(); i!=endi; ++i)
      {
        ColIterator endj = (*i).end();
        for (ColIterator j=(*i).begin(); j!=endj; ++j) {
          if (j.index() >= m) {
            dwarn << "WARNING: size of row "<< i.index()<<" is "<<j.offset()<<". But was specified as being "<< (*i).end().offset()
                  <<". This means you are wasting valuable space and creating additional cache misses!"<<std::endl;
            nnz_ -= ((*i).end().offset() - j.offset());
            r[i.index()].setsize(j.offset());
            break;
          }
        }
      }

      allocateData();
      setDataPointers();

      // if not, set matrix to built
      ready = built;
    }

    //===== implicit creation interface


    B& entry(size_type row, size_type col)
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (build_mode!=implicit)
        DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
      if (ready==built)
        DUNE_THROW(BCRSMatrixError,"matrix already built up, use operator[] for entry access now");
      if (ready==notAllocated)
        DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
      if (ready!=building)
        DUNE_THROW(InvalidStateException,"You may only use entry() during the 'building' stage");

      if (row >= n)
        DUNE_THROW(BCRSMatrixError,"row index exceeds matrix size");
      if (col >= m)
        DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
#endif

      size_type* begin = r[row].getindexptr();
      size_type* end = begin + r[row].getsize();

      size_type* pos = std::find(begin, end, col);

      //treat the case that there was a match in the array
      if (pos != end)
        if (*pos == col)
        {
          std::ptrdiff_t offset = pos - r[row].getindexptr();
          B* aptr = r[row].getptr() + offset;

          return *aptr;
        }

      //determine whether overflow has to be taken into account or not
      if (r[row].getsize() == avg)
        return overflow[std::make_pair(row,col)];
      else
      {
        //modify index array
        *end = col;

        //do simultaneous operations on data array a
        std::ptrdiff_t offset = end - r[row].getindexptr();
        B* apos = r[row].getptr() + offset;

        //increase rowsize
        r[row].setsize(r[row].getsize()+1);

        //return reference to the newly created entry
        return *apos;
      }
    }


    CompressionStatistics compress()
    {
      if (build_mode!=implicit)
        DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
      if (ready==built)
        DUNE_THROW(BCRSMatrixError,"matrix already built up, no more need for compression");
      if (ready==notAllocated)
        DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
      if (ready!=building)
        DUNE_THROW(InvalidStateException,"You may only call compress() at the end of the 'building' stage");

      //calculate statistics
      CompressionStatistics stats;
      stats.overflow_total = overflow.size();
      stats.maximum = 0;

      //get insertion iterators pointing to one before start (for later use of ++it)
      size_type* jiit = j_.get();
      B* aiit = a;

      //get iterator to the smallest overflow element
      typename OverflowType::iterator oit = overflow.begin();

      //store a copy of index pointers on which to perform sorting
      std::vector<size_type*> perm;

      //iterate over all rows and copy elements into their position in the compressed array
      for (size_type i=0; i<n; i++)
      {
        //get old pointers into a and j and size without overflow changes
        size_type* begin = r[i].getindexptr();
        //B* apos = r[i].getptr();
        size_type size = r[i].getsize();

        perm.resize(size);

        typename std::vector<size_type*>::iterator it = perm.begin();
        for (size_type* iit = begin; iit < begin + size; ++iit, ++it)
          *it = iit;

        //sort permutation array
        std::sort(perm.begin(),perm.end(), [](const size_type* l, const size_type* r){ return *l < *r; });

        //change row window pointer to their new positions
        r[i].setindexptr(jiit);
        r[i].setptr(aiit);

        for (it = perm.begin(); it != perm.end(); ++it)
        {
          //check whether there are elements in the overflow area which take precedence
          while ((oit!=overflow.end()) && (oit->first < std::make_pair(i,**it)))
          {
            //check whether there is enough memory to write to
            if (jiit > begin)
              DUNE_THROW(Dune::ImplicitModeCompressionBufferExhausted,
                         "Allocated memory for BCRSMatrix exhausted during compress()!"
                         "Please increase either the average number of entries per row or the compressionBufferSize value."
                         );
            //copy an element from the overflow area to the insertion position in a and j
            *jiit = oit->first.second;
            ++jiit;
            *aiit = oit->second;
            ++aiit;
            ++oit;
            r[i].setsize(r[i].getsize()+1);
          }

          //check whether there is enough memory to write to
          if (jiit > begin)
              DUNE_THROW(Dune::ImplicitModeCompressionBufferExhausted,
                         "Allocated memory for BCRSMatrix exhausted during compress()!"
                         "Please increase either the average number of entries per row or the compressionBufferSize value."
                         );

          //copy element from array
          *jiit = **it;
          ++jiit;
          B* apos = *it - j_.get() + a;
          *aiit = *apos;
          ++aiit;
        }

        //copy remaining elements from the overflow area
        while ((oit!=overflow.end()) && (oit->first.first == i))
        {
          //check whether there is enough memory to write to
          if (jiit > begin)
              DUNE_THROW(Dune::ImplicitModeCompressionBufferExhausted,
                         "Allocated memory for BCRSMatrix exhausted during compress()!"
                         "Please increase either the average number of entries per row or the compressionBufferSize value."
                         );

          //copy and element from the overflow area to the insertion position in a and j
          *jiit = oit->first.second;
          ++jiit;
          *aiit = oit->second;
          ++aiit;
          ++oit;
          r[i].setsize(r[i].getsize()+1);
        }

        // update maximum row size
        if (r[i].getsize()>stats.maximum)
          stats.maximum = r[i].getsize();
      }

      // overflow area may be cleared
      overflow.clear();

      //determine average number of entries and memory usage
      if ( n == 0)
      {
        stats.avg = 0;
        stats.mem_ratio = 1;
      }
      else
      {
        std::ptrdiff_t diff = (r[n-1].getindexptr() + r[n-1].getsize() - j_.get());
        nnz_ = diff;
        stats.avg = (double) (nnz_) / (double) n;
        stats.mem_ratio = (double) (nnz_) / (double) allocationSize_;
      }

      //matrix is now built
      ready = built;

      return stats;
    }

    //===== vector space arithmetic

    BCRSMatrix& operator*= (const field_type& k)
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
#endif

      if (nnz_ > 0)
      {
        // process 1D array
        for (size_type i=0; i<nnz_; i++)
          a[i] *= k;
      }
      else
      {
        RowIterator endi=end();
        for (RowIterator i=begin(); i!=endi; ++i)
        {
          ColIterator endj = (*i).end();
          for (ColIterator j=(*i).begin(); j!=endj; ++j)
            (*j) *= k;
        }
      }

      return *this;
    }

    BCRSMatrix& operator/= (const field_type& k)
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
#endif

      if (nnz_ > 0)
      {
        // process 1D array
        for (size_type i=0; i<nnz_; i++)
          a[i] /= k;
      }
      else
      {
        RowIterator endi=end();
        for (RowIterator i=begin(); i!=endi; ++i)
        {
          ColIterator endj = (*i).end();
          for (ColIterator j=(*i).begin(); j!=endj; ++j)
            (*j) /= k;
        }
      }

      return *this;
    }


    BCRSMatrix& operator+= (const BCRSMatrix& b)
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built || b.ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
      if(N()!=b.N() || M() != b.M())
        DUNE_THROW(RangeError, "Matrix sizes do not match!");
#endif
      RowIterator endi=end();
      ConstRowIterator j=b.begin();
      for (RowIterator i=begin(); i!=endi; ++i, ++j) {
        i->operator+=(*j);
      }

      return *this;
    }

    BCRSMatrix& operator-= (const BCRSMatrix& b)
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built || b.ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
      if(N()!=b.N() || M() != b.M())
        DUNE_THROW(RangeError, "Matrix sizes do not match!");
#endif
      RowIterator endi=end();
      ConstRowIterator j=b.begin();
      for (RowIterator i=begin(); i!=endi; ++i, ++j) {
        i->operator-=(*j);
      }

      return *this;
    }

    BCRSMatrix& axpy(field_type alpha, const BCRSMatrix& b)
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built || b.ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
      if(N()!=b.N() || M() != b.M())
        DUNE_THROW(RangeError, "Matrix sizes do not match!");
#endif
      RowIterator endi=end();
      ConstRowIterator j=b.begin();
      for(RowIterator i=begin(); i!=endi; ++i, ++j)
        i->axpy(alpha, *j);

      return *this;
    }

    //===== linear maps

    template<class X, class Y>
    void mv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
      if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,
                                 "Size mismatch: M: " << N() << "x" << M() << " x: " << x.N());
      if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,
                                 "Size mismatch: M: " << N() << "x" << M() << " y: " << y.N());
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
      {
        y[i.index()]=0;
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
        {
          auto&& xj = Impl::asVector(x[j.index()]);
          auto&& yi = Impl::asVector(y[i.index()]);
          Impl::asMatrix(*j).umv(xj, yi);
        }
      }
    }

    template<class X, class Y>
    void umv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
      if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
      if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
        {
          auto&& xj = Impl::asVector(x[j.index()]);
          auto&& yi = Impl::asVector(y[i.index()]);
          Impl::asMatrix(*j).umv(xj,yi);
        }
      }
    }

    template<class X, class Y>
    void mmv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
      if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
      if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
        {
          auto&& xj = Impl::asVector(x[j.index()]);
          auto&& yi = Impl::asVector(y[i.index()]);
          Impl::asMatrix(*j).mmv(xj,yi);
        }
      }
    }

    template<class X, class Y, class F>
    void usmv (F&& alpha, const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
      if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
      if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
        {
          auto&& xj = Impl::asVector(x[j.index()]);
          auto&& yi = Impl::asVector(y[i.index()]);
          Impl::asMatrix(*j).usmv(alpha,xj,yi);
        }
      }
    }

    template<class X, class Y>
    void mtv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
      if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
      if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
      for(size_type i=0; i<y.N(); ++i)
        y[i]=0;
      umtv(x,y);
    }

    template<class X, class Y>
    void umtv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
      if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
      if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
        {
          auto&& xi = Impl::asVector(x[i.index()]);
          auto&& yj = Impl::asVector(y[j.index()]);
          Impl::asMatrix(*j).umtv(xi,yj);
        }
      }
    }

    template<class X, class Y>
    void mmtv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
      if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
        {
          auto&& xi = Impl::asVector(x[i.index()]);
          auto&& yj = Impl::asVector(y[j.index()]);
          Impl::asMatrix(*j).mmtv(xi,yj);
        }
      }
    }

    template<class X, class Y>
    void usmtv (const field_type& alpha, const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
      if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
      if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
        {
          auto&& xi = Impl::asVector(x[i.index()]);
          auto&& yj = Impl::asVector(y[j.index()]);
          Impl::asMatrix(*j).usmtv(alpha,xi,yj);
        }
      }
    }

    template<class X, class Y>
    void umhv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
      if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
      if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
        {
          auto&& xi = Impl::asVector(x[i.index()]);
          auto&& yj = Impl::asVector(y[j.index()]);
          Impl::asMatrix(*j).umhv(xi,yj);
        }
      }
    }

    template<class X, class Y>
    void mmhv (const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
      if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
      if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
        {
          auto&& xi = Impl::asVector(x[i.index()]);
          auto&& yj = Impl::asVector(y[j.index()]);
          Impl::asMatrix(*j).mmhv(xi,yj);
        }
      }
    }

    template<class X, class Y>
    void usmhv (const field_type& alpha, const X& x, Y& y) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
      if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
      if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
#endif
      ConstRowIterator endi=end();
      for (ConstRowIterator i=this->begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
        {
          auto&& xi = Impl::asVector(x[i.index()]);
          auto&& yj = Impl::asVector(y[j.index()]);
          Impl::asMatrix(*j).usmhv(alpha,xi,yj);
        }
      }
    }


    //===== norms

    typename FieldTraits<field_type>::real_type frobenius_norm2 () const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
#endif

      typename FieldTraits<field_type>::real_type sum=0;

      for (auto&& row : (*this))
        for (auto&& entry : row)
          sum += Impl::asMatrix(entry).frobenius_norm2();

      return sum;
    }

    typename FieldTraits<field_type>::real_type frobenius_norm () const
    {
      return sqrt(frobenius_norm2());
    }

    template <typename ft = field_type,
              typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
    typename FieldTraits<ft>::real_type infinity_norm() const {
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");

      using real_type = typename FieldTraits<ft>::real_type;
      using std::max;

      real_type norm = 0;
      for (auto const &x : *this) {
        real_type sum = 0;
        for (auto const &y : x)
          sum += Impl::asMatrix(y).infinity_norm();
        norm = max(sum, 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 {
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");

      using real_type = typename FieldTraits<ft>::real_type;
      using std::max;

      real_type norm = 0;
      for (auto const &x : *this) {
        real_type sum = 0;
        for (auto const &y : x)
          sum += Impl::asMatrix(y).infinity_norm_real();
        norm = max(sum, 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 {
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");

      using real_type = typename FieldTraits<ft>::real_type;
      using std::max;

      real_type norm = 0;
      real_type is_nan = 1;
      for (auto const &x : *this) {
        real_type sum = 0;
        for (auto const &y : x)
          sum += Impl::asMatrix(y).infinity_norm();
        norm = max(sum, norm);
        is_nan += sum;
      }

      return norm * (is_nan / is_nan);
    }

    template <typename ft = field_type,
              typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
    typename FieldTraits<ft>::real_type infinity_norm_real() const {
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");

      using real_type = typename FieldTraits<ft>::real_type;
      using std::max;

      real_type norm = 0;
      real_type is_nan = 1;

      for (auto const &x : *this) {
        real_type sum = 0;
        for (auto const &y : x)
          sum += Impl::asMatrix(y).infinity_norm_real();
        norm = max(sum, norm);
        is_nan += sum;
      }

      return norm * (is_nan / is_nan);
    }

    //===== sizes

    size_type N () const
    {
      return n;
    }

    size_type M () const
    {
      return m;
    }

    size_type nonzeroes () const
    {
      // in case of row-wise allocation
      if( nnz_ <= 0 )
        nnz_ = std::accumulate( begin(), end(), size_type( 0 ), [] ( size_type s, const row_type &row ) { return s+row.getsize(); } );
      return nnz_;
    }

    BuildStage buildStage() const
    {
      return ready;
    }

    BuildMode buildMode() const
    {
      return build_mode;
    }

    //===== query

    bool exists (size_type i, size_type j) const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (i<0 || i>=n) DUNE_THROW(BCRSMatrixError,"row index out of range");
      if (j<0 || j>=m) DUNE_THROW(BCRSMatrixError,"column index out of range");
#endif
      return (r[i].size() && r[i].find(j) != r[i].end());
    }


  protected:
    // state information
    BuildMode build_mode;     // row wise or whole matrix
    BuildStage ready;               // indicate the stage the matrix building is in

    // The allocator used for memory management
    typename std::allocator_traits<A>::template rebind_alloc<B> allocator_;

    typename std::allocator_traits<A>::template rebind_alloc<row_type> rowAllocator_;

    typename std::allocator_traits<A>::template rebind_alloc<size_type> sizeAllocator_;

    // size of the matrix
    size_type n;       // number of rows
    size_type m;       // number of columns
    mutable size_type nnz_;     // number of nonzeroes contained in the matrix
    size_type allocationSize_; //allocated size of a and j arrays, except in implicit mode: nnz_==allocationsSize_
    // zero means that memory is allocated separately for each row.

    // the rows are dynamically allocated
    row_type* r;     // [n] the individual rows having pointers into a,j arrays

    // dynamically allocated memory
    B*   a;      // [allocationSize] non-zero entries of the matrix in row-wise ordering
    // If a single array of column indices is used, it can be shared
    // between different matrices with the same sparsity pattern
    std::shared_ptr<size_type> j_;  // [allocationSize] column indices of entries

    // additional data is needed in implicit buildmode
    size_type avg;
    double compressionBufferSize_;

    typedef std::map<std::pair<size_type,size_type>, B> OverflowType;
    OverflowType overflow;

    void setWindowPointers(ConstRowIterator row)
    {
      row_type current_row(a,j_.get(),0); // Pointers to current row data
      for (size_type i=0; i<n; i++, ++row) {
        // set row i
        size_type s = row->getsize();

        if (s>0) {
          // setup pointers and size
          r[i].set(s,current_row.getptr(), current_row.getindexptr());
          // update pointer for next row
          current_row.setptr(current_row.getptr()+s);
          current_row.setindexptr(current_row.getindexptr()+s);
        } else{
          // empty row
          r[i].set(0,nullptr,nullptr);
        }
      }
    }


    void setColumnPointers(ConstRowIterator row)
    {
      size_type* jptr = j_.get();
      for (size_type i=0; i<n; ++i, ++row) {
        // set row i
        size_type s = row->getsize();

        if (s>0) {
          // setup pointers and size
          r[i].setsize(s);
          r[i].setindexptr(jptr);
        } else{
          // empty row
          r[i].set(0,nullptr,nullptr);
        }

        // advance position in global array
        jptr += s;
      }
    }


    void setDataPointers()
    {
      B* aptr = a;
      for (size_type i=0; i<n; ++i) {
        // set row i
        if (r[i].getsize() > 0) {
          // setup pointers and size
          r[i].setptr(aptr);
        } else{
          // empty row
          r[i].set(0,nullptr,nullptr);
        }

        // advance position in global array
        aptr += r[i].getsize();
      }
    }

    void copyWindowStructure(const BCRSMatrix& Mat)
    {
      setWindowPointers(Mat.begin());

      // copy data
      for (size_type i=0; i<n; i++) r[i] = Mat.r[i];

      // finish off
      build_mode = row_wise; // dummy
      ready = built;
    }

    void deallocate(bool deallocateRows=true)
    {

      if (notAllocated)
        return;

      if (allocationSize_>0)
      {
        // a,j_ have been allocated as one long vector
        j_.reset();
        if (a)
          {
            for(B *aiter=a+(allocationSize_-1), *aend=a-1; aiter!=aend; --aiter)
              std::allocator_traits<decltype(allocator_)>::destroy(allocator_, aiter);
            allocator_.deallocate(a,allocationSize_);
            a = nullptr;
          }
      }
      else if (r)
      {
        // check if memory for rows have been allocated individually
        for (size_type i=0; i<n; i++)
          if (r[i].getsize()>0)
          {
            for (B *col=r[i].getptr()+(r[i].getsize()-1),
                 *colend = r[i].getptr()-1; col!=colend; --col) {
              std::allocator_traits<decltype(allocator_)>::destroy(allocator_, col);
            }
            sizeAllocator_.deallocate(r[i].getindexptr(),1);
            allocator_.deallocate(r[i].getptr(),1);
            // clear out row data in case we don't want to deallocate the rows
            // otherwise we might run into a double free problem here later
            r[i].set(0,nullptr,nullptr);
          }
      }

      // deallocate the rows
      if (n>0 && deallocateRows && r) {
        for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
          std::allocator_traits<decltype(rowAllocator_)>::destroy(rowAllocator_, riter);
        rowAllocator_.deallocate(r,n);
        r = nullptr;
      }

      // Mark matrix as not built at all.
      ready=notAllocated;

    }

    void allocate(size_type rows, size_type columns, size_type allocationSize, bool allocateRows, bool allocate_data)
    {
      // Store size
      n = rows;
      m = columns;
      nnz_ = allocationSize;
      allocationSize_ = allocationSize;

      // allocate rows
      if(allocateRows) {
        if (n>0) {
          if (r)
            DUNE_THROW(InvalidStateException,"Rows have already been allocated, cannot allocate a second time");
          r = rowAllocator_.allocate(rows);
          // initialize row entries
          for(row_type* ri=r; ri!=r+rows; ++ri)
            std::allocator_traits<decltype(rowAllocator_)>::construct(rowAllocator_, ri, row_type());
        }else{
          r = 0;
        }
      }

      // allocate a and j_ array
      if (allocate_data)
        allocateData();
      // allocate column indices only if not yet present (enable sharing)
      if (allocationSize_>0) {
        // we copy allocator and size to the deleter since _j may outlive this class
        if (!j_.get())
          j_.reset(sizeAllocator_.allocate(allocationSize_),
            [alloc = sizeAllocator_, size = allocationSize_](auto ptr) mutable {
              alloc.deallocate(ptr, size);
            });
      }else{
        j_.reset();
      }

      // Mark the matrix as not built.
      ready = building;
    }

    void allocateData()
    {
      if (a)
        DUNE_THROW(InvalidStateException,"Cannot allocate data array (already allocated)");
      if (allocationSize_>0) {
        a = allocator_.allocate(allocationSize_);
        // use placement new to call constructor that allocates
        // additional memory.
        new (a) B[allocationSize_];
      } else {
        a = nullptr;
      }
    }

    void implicit_allocate(size_type _n, size_type _m)
    {
      if (build_mode != implicit)
        DUNE_THROW(InvalidStateException,"implicit_allocate() may only be called in implicit build mode");
      if (ready != notAllocated)
        DUNE_THROW(InvalidStateException,"memory has already been allocated");

      // check to make sure the user has actually set the parameters
      if (compressionBufferSize_ < 0)
        DUNE_THROW(InvalidStateException,"You have to set the implicit build mode parameters before starting to build the matrix");
      //calculate size of overflow region, add buffer for row sort!
      size_type osize = (size_type) (_n*avg)*compressionBufferSize_ + 4*avg;
      allocationSize_ = _n*avg + osize;

      allocate(_n, _m, allocationSize_,true,true);

      //set row pointers correctly
      size_type* jptr = j_.get() + osize;
      B* aptr = a + osize;
      for (size_type i=0; i<n; i++)
      {
        r[i].set(0,aptr,jptr);
        jptr = jptr + avg;
        aptr = aptr + avg;
      }

      ready = building;
    }
  };


  template<class B, class A>
  struct FieldTraits< BCRSMatrix<B, A> >
  {
    using field_type = typename BCRSMatrix<B, A>::field_type;
    using real_type = typename FieldTraits<field_type>::real_type;
  };

} // end namespace

#endif