bcrsmatrix.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_BCRSMATRIX_HH
#define DUNE_BCRSMATRIX_HH
 
#include <cmath>
#include <complex>
#include <set>
#include <iostream>
#include <algorithm>
#include <numeric>
#include <vector>
#include <map>
 
#include "istlexception.hh"
#include "bvector.hh"
#include "matrixutils.hh"
#include <dune/common/shared_ptr.hh>
#include <dune/common/stdstreams.hh>
#include <dune/common/iteratorfacades.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/ftraits.hh>
#include <dune/common/static_assert.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
 
    typedef typename B::field_type field_type;
 
    typedef B block_type;
 
    typedef A allocator_type;
 
    typedef CompressedBlockVectorWindow<B,A> row_type;
 
    typedef typename A::size_type size_type;
 
    typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;
 
    enum {
      blocklevel = B::blocklevel+1
    };
 
    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");
      if (r[i].getptr()==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
#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 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 overflowsize 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), overflowsize(-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), overflowsize(-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), overflowsize(-1.0)
    {
      allocate(_n, _m,0,true,false);
    }
 
 
    BCRSMatrix (size_type _n, size_type _m, size_type _avg, double _overflowsize, BuildMode bm)
      : build_mode(bm), ready(notAllocated), n(0), m(0), nnz(0),
        allocationSize(0), r(0), a(0),
        avg(_avg), overflowsize(_overflowsize)
    {
      if (bm != implicit)
        DUNE_THROW(BCRSMatrixError,"Only call this constructor when using the implicit build mode");
      // Prevent user from setting a negative overflowsize:
      // 1) It doesn't make sense
      // 2) We use a negative overflow value to indicate that the parameters
      //    have not been set yet
      if (_overflowsize < 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), overflowsize(Mat.overflowsize)
    {
      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.nnz;
 
      // in case of row-wise allocation
      if (_nnz<=0)
      {
        _nnz = 0;
        for (size_type i=0; i<Mat.n; i++)
          _nnz += Mat.r[i].getsize();
      }
 
      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 _overflow)
    {
      // Prevent user from setting a negative overflowsize:
      // 1) It doesn't make sense
      // 2) We use a negative overflow value to indicate that the parameters
      //    have not been set yet
      if (_overflow < 0.0)
        DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
 
      // 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;
      overflowsize = _overflow;
    }
 
    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
      deallocate(false);
 
      // 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)
          rowAllocator_.destroy(riter);
        rowAllocator_.deallocate(r,n);
      }
 
      nnz=Mat.nnz;
      if (nnz<=0)
      {
        for (size_type i=0; i<Mat.n; i++)
          nnz += Mat.r[i].getsize();
      }
 
      // 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");
      }
 
      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*   a = Mat.allocator_.allocate(s);
            // use placement new to call constructor that allocates
            // additional memory.
            new (a) B[s];
            size_type* j = Mat.sizeAllocator_.allocate(s);
            Mat.r[i].set(s,a,j);
          }
        }else
          // setup empty row
          Mat.r[i].set(0,0,0);
 
        // initialize the j array for row i from pattern
        size_type k=0;
        size_type *j =  Mat.r[i].getindexptr();
        for (typename PatternType::const_iterator it=pattern.begin(); it!=pattern.end(); ++it)
          j[k++] = *it;
 
        // 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)
      {
        if (pattern.find(j)!=pattern.end())
          return true;
        else
          return false;
      }
      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
      typedef std::set<size_type,std::less<size_type> > PatternType;
      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 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;
            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 simulatenous 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 sortation
      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(),PointerCompare<size_type>());
 
        //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::ImplicitModeOverflowExhausted,
                         "Allocated memory for BCRSMatrix exhausted during compress()!"
                         "Please increase either the average number of entries per row or the overflow fraction."
                         );
            //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);
          }
 
          //check whether there is enough memory to write to
          if (jiit > begin)
              DUNE_THROW(Dune::ImplicitModeOverflowExhausted,
                         "Allocated memory for BCRSMatrix exhausted during compress()!"
                         "Please increase either the average number of entries per row or the overflow fraction."
                         );
 
          //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::ImplicitModeOverflowExhausted,
                         "Allocated memory for BCRSMatrix exhausted during compress()!"
                         "Please increase either the average number of entries per row or the overflow fraction."
                         );
 
          //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
      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=begin(); i!=endi; ++i)
      {
        y[i.index()]=0;
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).umv(x[j.index()],y[i.index()]);
      }
    }
 
    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=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).umv(x[j.index()],y[i.index()]);
      }
    }
 
    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=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).mmv(x[j.index()],y[i.index()]);
      }
    }
 
    template<class X, class Y>
    void usmv (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()!=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=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).usmv(alpha,x[j.index()],y[i.index()]);
      }
    }
 
    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=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).umtv(x[i.index()],y[j.index()]);
      }
    }
 
    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=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).mmtv(x[i.index()],y[j.index()]);
      }
    }
 
    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=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).usmtv(alpha,x[i.index()],y[j.index()]);
      }
    }
 
    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=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).umhv(x[i.index()],y[j.index()]);
      }
    }
 
    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=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).mmhv(x[i.index()],y[j.index()]);
      }
    }
 
    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=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          (*j).usmhv(alpha,x[i.index()],y[j.index()]);
      }
    }
 
 
    //===== 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
 
      double sum=0;
 
      ConstRowIterator endi=end();
      for (ConstRowIterator i=begin(); i!=endi; ++i)
      {
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          sum += (*j).frobenius_norm2();
      }
 
      return sum;
    }
 
    typename FieldTraits<field_type>::real_type frobenius_norm () const
    {
      return sqrt(frobenius_norm2());
    }
 
    typename FieldTraits<field_type>::real_type infinity_norm () const
    {
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
 
      double max=0;
      ConstRowIterator endi=end();
      for (ConstRowIterator i=begin(); i!=endi; ++i)
      {
        double sum=0;
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          sum += (*j).infinity_norm();
        max = std::max(max,sum);
      }
      return max;
    }
 
    typename FieldTraits<field_type>::real_type infinity_norm_real () const
    {
#ifdef DUNE_ISTL_WITH_CHECKING
      if (ready != built)
        DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
#endif
 
      double max=0;
      ConstRowIterator endi=end();
      for (ConstRowIterator i=begin(); i!=endi; ++i)
      {
        double sum=0;
        ConstColIterator endj = (*i).end();
        for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
          sum += (*j).infinity_norm_real();
        max = std::max(max,sum);
      }
      return max;
    }
 
 
    //===== sizes
 
    size_type N () const
    {
      return n;
    }
 
    size_type M () const
    {
      return m;
    }
 
    size_type nonzeroes () const
    {
      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
      if (r[i].size() && r[i].find(j)!=r[i].end())
        return true;
      else
        return false;
    }
 
 
  private:
    // 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 A::template rebind<B>::other allocator_;
 
    typename A::template rebind<row_type>::other rowAllocator_;
 
    typename A::template rebind<size_type>::other sizeAllocator_;
 
    // size of the matrix
    size_type n;       // number of rows
    size_type m;       // number of columns
    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
    Dune::shared_ptr<size_type> j;  // [allocationSize] column indices of entries
 
    // additional data is needed in implicit buildmode
    size_type avg;
    double overflowsize;
 
    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,0,0);
        }
      }
    }
 
 
    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,0,0);
        }
 
        // 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,0,0);
        }
 
        // 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)
              allocator_.destroy(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) {
              allocator_.destroy(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)
          rowAllocator_.destroy(riter);
        rowAllocator_.deallocate(r,n);
        r = nullptr;
      }
 
      // Mark matrix as not built at all.
      ready=notAllocated;
 
    }
 
    class Deallocator
    {
      typename A::template rebind<size_type>::other& sizeAllocator_;
 
    public:
      Deallocator(typename A::template rebind<size_type>::other& sizeAllocator)
        : sizeAllocator_(sizeAllocator)
      {}
 
      void operator()(size_type* p) { sizeAllocator_.deallocate(p,1); }
    };
 
 
    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);
        }else{
          r = 0;
        }
      }
 
      // allocate a and j array
      if (allocate_data)
        allocateData();
      if (allocationSize>0) {
        // allocate column indices only if not yet present (enable sharing)
        if (!j.get())
          j.reset(sizeAllocator_.allocate(allocationSize),Deallocator(sizeAllocator_));
      }else{
        j.reset();
        for(row_type* ri=r; ri!=r+rows; ++ri)
          rowAllocator_.construct(ri, row_type());
      }
 
      // 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 (overflowsize < 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)*overflowsize + 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;
    }
  };
 
 
} // end namespace
 
#endif