sparse_matrix_bcsr.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/base_header.hpp>
#include <kernel/util/assertion.hpp>
#include <kernel/lafem/forward.hpp>
#include <kernel/lafem/container.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/sparse_layout.hpp>
#include <kernel/lafem/arch/scale.hpp>
#include <kernel/lafem/arch/axpy.hpp>
#include <kernel/lafem/arch/apply.hpp>
#include <kernel/lafem/arch/lumping.hpp>
#include <kernel/lafem/arch/norm.hpp>
#include <kernel/lafem/arch/scale_row_col.hpp>
#include <kernel/lafem/arch/row_norm.hpp>
#include <kernel/lafem/arch/diagonal.hpp>
#include <kernel/adjacency/graph.hpp>
#include <kernel/util/tiny_algebra.hpp>
#include <kernel/util/statistics.hpp>
#include <kernel/util/time_stamp.hpp>
#include <kernel/util/likwid_marker.hpp>
 
#include <fstream>
#include <set>
 
namespace FEAT
{
  namespace LAFEM
  {
    namespace Intern
    {
      template<typename DT_, typename IT_, int size_>
      struct BCSRVectorHelper
      {
        static_assert(size_ > 1, "invalid block size");
        typedef DenseVectorBlocked<DT_, IT_, size_> VectorType;
      };
 
      template<typename DT_, typename IT_>
      struct BCSRVectorHelper<DT_, IT_, 1>
      {
        typedef DenseVector< DT_, IT_> VectorType;
      };
 
      template<typename DT_, int BlockHeight_, int BlockWidth_, Perspective perspective_>
      struct BCSRPerspectiveHelper
      {
        typedef Tiny::Matrix<DT_, BlockHeight_, BlockWidth_> Type;
      };
 
      template<typename DT_, int BlockHeight_, int BlockWidth_>
      struct BCSRPerspectiveHelper<DT_, BlockHeight_, BlockWidth_, Perspective::pod>
      {
        typedef DT_ Type;
      };
    } // namespace Intern
 
    template <typename DT_, typename IT_, int BlockHeight_, int BlockWidth_>
    class SparseMatrixBCSR : public Container<DT_, IT_>
    {
      static_assert(BlockHeight_ > 0, "invalid block size");
      static_assert(BlockWidth_ > 0, "invalid block size");
 
    public:
      typedef DT_ DataType;
      typedef IT_ IndexType;
      static constexpr int BlockHeight = BlockHeight_;
      static constexpr int BlockWidth = BlockWidth_;
      typedef Tiny::Matrix<DataType, BlockHeight, BlockWidth> ValueType;
 
      static constexpr SparseLayoutId layout_id = SparseLayoutId::lt_csr;
 
      typedef const IT_* ImageIterator;
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      using ContainerType = SparseMatrixBCSR<DT2_, IT2_, BlockHeight_, BlockWidth_>;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
      typedef typename Intern::BCSRVectorHelper<DT_, IT_, BlockHeight_>::VectorType VectorTypeL;
      typedef typename Intern::BCSRVectorHelper<DT_, IT_, BlockWidth_>::VectorType VectorTypeR;
 
      static constexpr bool is_global = false;
      static constexpr bool is_local = true;
 
      class ScatterAxpy
      {
      public:
        typedef LAFEM::SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_> MatrixType;
        typedef DT_ DataType;
        typedef IT_ IndexType;
        typedef Tiny::Matrix<DT_, BlockHeight_, BlockWidth_> ValueType;
 
      private:
#ifdef DEBUG
        const IT_ _deadcode;
#endif
        Index _num_rows;
        Index _num_cols;
        IT_* _row_ptr;
        IT_* _col_idx;
        IT_* _col_ptr;
        ValueType* _data;
 
      public:
        explicit ScatterAxpy(MatrixType& matrix) :
#ifdef DEBUG
          _deadcode(~IT_(0)),
#endif
          _num_rows(matrix.rows()),
          _num_cols(matrix.columns()),
          _row_ptr(matrix.row_ptr()),
          _col_idx(matrix.col_ind()),
          _col_ptr(nullptr),
          _data(matrix.val())
        {
          // allocate column-pointer array
          _col_ptr = new IT_[_num_cols];
#ifdef DEBUG
          for(Index i(0); i < _num_cols; ++i)
          {
            _col_ptr[i] = _deadcode;
          }
#endif
        }
 
        virtual ~ScatterAxpy()
        {
          if(_col_ptr != nullptr)
          {
            delete [] _col_ptr;
          }
        }
 
        template<typename LocalMatrix_, typename RowMapping_, typename ColMapping_>
        void operator()(const LocalMatrix_& loc_mat, const RowMapping_& row_map,
                        const ColMapping_& col_map, DT_ alpha = DT_(1))
        {
          // loop over all local row entries
          for(int i(0); i < row_map.get_num_local_dofs(); ++i)
          {
            // fetch row index
            const Index ix = row_map.get_index(i);
 
            // build column pointer for this row entry contribution
            for(IT_ k(_row_ptr[ix]); k < _row_ptr[ix + 1]; ++k)
            {
              _col_ptr[_col_idx[k]] = k;
            }
 
            // loop over all local column entries
            for(int j(0); j < col_map.get_num_local_dofs(); ++j)
            {
              // fetch column index
              const Index jx = col_map.get_index(j);
 
              // ensure that the column pointer is valid for this index
              ASSERTM(_col_ptr[jx] != _deadcode, "invalid column index");
 
              // incorporate data into global matrix
              _data[_col_ptr[jx]] += alpha * loc_mat[i][j];
 
              // continue with next column entry
            }
 
#ifdef DEBUG
            // reformat column-pointer array
            for(IT_ k(_row_ptr[ix]); k < _row_ptr[ix + 1]; ++k)
            {
              _col_ptr[_col_idx[k]] = _deadcode;
            }
#endif
            // continue with next row entry
          }
        }
      }; // class ScatterAxpy
 
    private:
      Index & _size()
      {
        return this->_scalar_index.at(0);
      }
 
      Index & _rows()
      {
        return this->_scalar_index.at(1);
      }
 
      Index & _columns()
      {
        return this->_scalar_index.at(2);
      }
 
      Index & _used_elements()
      {
        return this->_scalar_index.at(3);
      }
 
    public:
 
      explicit SparseMatrixBCSR() :
        Container<DT_, IT_> (0)
      {
        this->_scalar_index.push_back(0);
        this->_scalar_index.push_back(0);
        this->_scalar_index.push_back(0);
      }
 
      explicit SparseMatrixBCSR(Index rows_in, Index columns_in) :
        Container<DT_, IT_> (rows_in * columns_in)
      {
        this->_scalar_index.push_back(rows_in);
        this->_scalar_index.push_back(columns_in);
        this->_scalar_index.push_back(0);
      }
 
      explicit SparseMatrixBCSR(Index rows_in, Index columns_in, Index used_elements_in) :
        Container<DT_, IT_> (rows_in * columns_in)
      {
        XASSERT(rows_in != Index(0) && columns_in != Index(0));
 
        this->_scalar_index.push_back(rows_in);
        this->_scalar_index.push_back(columns_in);
        this->_scalar_index.push_back(used_elements_in);
 
        this->_indices.push_back(MemoryPool::template allocate_memory<IT_>(_used_elements()));
        this->_indices_size.push_back(_used_elements());
 
        this->_indices.push_back(MemoryPool::template allocate_memory<IT_>(_rows() + 1));
        this->_indices_size.push_back(_rows() + 1);
 
        this->_elements.push_back(MemoryPool::template allocate_memory<DT_>(used_elements<Perspective::pod>()));
        this->_elements_size.push_back(used_elements<Perspective::pod>());
 
      }
 
      explicit SparseMatrixBCSR(const SparseLayout<IT_, layout_id> & layout_in) :
        Container<DT_, IT_> (layout_in._scalar_index.at(0))
      {
        this->_indices.assign(layout_in._indices.begin(), layout_in._indices.end());
        this->_indices_size.assign(layout_in._indices_size.begin(), layout_in._indices_size.end());
        this->_scalar_index.assign(layout_in._scalar_index.begin(), layout_in._scalar_index.end());
 
        for (auto i : this->_indices)
          MemoryPool::increase_memory(i);
 
        this->_elements.push_back(MemoryPool::template allocate_memory<DT_>(used_elements<Perspective::pod>()));
        this->_elements_size.push_back(used_elements<Perspective::pod>());
      }
 
      explicit SparseMatrixBCSR(const Adjacency::Graph & graph) :
        SparseMatrixBCSR(graph.get_num_nodes_domain(), graph.get_num_nodes_image(), graph.get_num_indices())
      {
        const Index num_nnze = graph.get_num_indices();
        const Index num_rows = graph.get_num_nodes_domain();
        const Index * dom_ptr(graph.get_domain_ptr());
        const Index * img_idx(graph.get_image_idx());
        IT_ * prow_ptr(this->row_ptr());
        IT_ * pcol_idx(this->col_ind());
 
        FEAT_PRAGMA_OMP(parallel for)
        for(Index i = 0; i <= num_rows; ++i)
          prow_ptr[i] = IT_(dom_ptr[i]);
 
        FEAT_PRAGMA_OMP(parallel for)
        for(Index i = 0; i < num_nnze; ++i)
          pcol_idx[i] = IT_(img_idx[i]);
      }
 
      explicit SparseMatrixBCSR(FileMode mode, const String& filename) :
        Container< DT_, IT_>(0)
      {
        read_from(mode, filename);
      }
 
      explicit SparseMatrixBCSR(FileMode mode, std::istream& file) :
        Container<DT_, IT_>(0)
      {
        read_from(mode, file);
      }
 
      explicit SparseMatrixBCSR(const Index rows_in, const Index columns_in,
                                      DenseVector< IT_, IT_> & col_ind_in, DenseVector< DT_, IT_> & val_in, DenseVector< IT_, IT_> & row_ptr_in) :
        Container<DT_, IT_>(rows_in * columns_in)
      {
        XASSERT(col_ind_in.size() > 0);
        XASSERT(val_in.size() > 0);
        XASSERT(row_ptr_in.size() > 0);
 
        XASSERT(rows_in != Index(0) && columns_in != Index(0));
        XASSERTM(val_in.size() % (BlockHeight_ * BlockWidth_) == 0, "input values size is not a multiple of container blocksize!");
 
        this->_scalar_index.push_back(rows_in);
        this->_scalar_index.push_back(columns_in);
        this->_scalar_index.push_back(val_in.size() / Index(BlockHeight_ * BlockWidth_));
 
        this->_elements.push_back(val_in.elements());
        this->_elements_size.push_back(val_in.size());
        this->_indices.push_back(col_ind_in.elements());
        this->_indices_size.push_back(col_ind_in.size());
        this->_indices.push_back(row_ptr_in.elements());
        this->_indices_size.push_back(row_ptr_in.size());
 
        for (Index i(0) ; i < this->_elements.size() ; ++i)
          MemoryPool::increase_memory(this->_elements.at(i));
        for (Index i(0) ; i < this->_indices.size() ; ++i)
          MemoryPool::increase_memory(this->_indices.at(i));
      }
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      explicit SparseMatrixBCSR(std::vector<char> input) :
        Container<DT_, IT_>(0)
      {
        deserialize<DT2_, IT2_>(input);
      }
 
      SparseMatrixBCSR(SparseMatrixBCSR && other) :
        Container<DT_, IT_>(std::forward<SparseMatrixBCSR>(other))
      {
      }
 
      SparseMatrixBCSR & operator= (SparseMatrixBCSR && other)
      {
        this->move(std::forward<SparseMatrixBCSR>(other));
 
        return *this;
      }
 
      SparseMatrixBCSR clone(CloneMode clone_mode = CloneMode::Weak) const
      {
        SparseMatrixBCSR t;
        t.clone(*this, clone_mode);
        return t;
      }
 
      template<typename DT2_, typename IT2_>
      void clone(const SparseMatrixBCSR<DT2_, IT2_, BlockHeight_, BlockWidth_> & other, CloneMode clone_mode = CloneMode::Weak)
      {
        Container<DT_, IT_>::clone(other, clone_mode);
      }
 
      template <typename DT2_, typename IT2_>
      void convert(const SparseMatrixBCSR<DT2_, IT2_, BlockHeight_, BlockWidth_> & other)
      {
        this->assign(other);
      }
 
      template <typename DT2_, typename IT2_>
      void convert(const SparseMatrixCSR<DT2_, IT2_>& other)
      {
        this->clear();
        // setup internal sizes
        Index n_rows = Index(other.rows()) / Index(BlockHeight);
        Index n_cols = Index(other.columns()) / Index(BlockWidth);
 
        // temprorary arrays for row and column ptr
        std::vector<IndexType> tmp_row_ptr(n_rows+1u);
        std::vector<IndexType> tmp_col_ptr;
 
        // ptrs of csr matrix
        const IT2_* _crs_row = other.row_ptr();
        const IT2_* _crs_col = other.col_ind();
        const DT2_ * _crs_val = other.val();
 
        // reserve to sensical size, we will have at most used_element blocks
        tmp_col_ptr.reserve(other.used_elements());
 
        // for each block row, gather which blocks should be actually allocated
        std::set<IndexType> tmp_uni_col;
 
        tmp_row_ptr[0] = IndexType(0);
        for(Index b_row = 0; b_row < n_rows; ++b_row)
        {
          for(Index l_row = b_row*Index(BlockHeight); l_row < (b_row+1)*Index(BlockHeight); ++l_row)
          {
            for(IT2_ idx = _crs_row[l_row]; idx < _crs_row[l_row+1]; ++idx)
            {
              tmp_uni_col.insert(IndexType(_crs_col[idx])/IndexType(BlockWidth));
            }
          }
          tmp_row_ptr[b_row+1] = IndexType(tmp_uni_col.size()) + tmp_row_ptr[b_row];
          std::for_each(tmp_uni_col.begin(), tmp_uni_col.end(), [&](const auto& ele){tmp_col_ptr.push_back(ele);});
          // std::printf("Brow %i:  ", int(b_row));
          // std::for_each(tmp_uni_col.begin(), tmp_uni_col.end(), [&](const auto& ele){std::cout << ele << ", ";});
          // std::cout << "\n";
          tmp_uni_col.clear();
        }
 
        Index elements_size = tmp_col_ptr.size();
        // std::printf("Elements size %i\n", int(elements_size));
        // std::printf("Tmp col ptr:\n");
        // std::for_each(tmp_col_ptr.begin(), tmp_col_ptr.end(), [&](const auto& ele){std::cout << ele << ", ";});
        // std::cout << "\n";
 
        // allocate actual matrix data
        this->_scalar_index.push_back(n_rows*n_cols);
        this->_scalar_index.push_back(n_rows);
        this->_scalar_index.push_back(n_cols);
        this->_scalar_index.push_back(elements_size);
 
        this->_indices.push_back(MemoryPool::template allocate_memory<IT_>(_used_elements()));
        this->_indices_size.push_back(_used_elements());
 
        this->_indices.push_back(MemoryPool::template allocate_memory<IT_>(_rows() + 1));
        this->_indices_size.push_back(_rows() + 1);
 
        this->_elements.push_back(MemoryPool::template allocate_memory<DT_>(used_elements<Perspective::pod>()));
        this->_elements_size.push_back(used_elements<Perspective::pod>());
 
        auto* row_ptr = this->_indices.at(1);
        auto* col = this->_indices.at(0);
        auto* val = this->_elements.at(0);
        // now, fill up array data
        std::copy(tmp_row_ptr.begin(), tmp_row_ptr.end(), row_ptr);
        std::copy(tmp_col_ptr.begin(), tmp_col_ptr.end(), col);
        // std::printf("Row ptr:\n");
        // std::for_each(row_ptr, row_ptr+n_rows+1, [&](const auto& ele){std::cout << ele << ", ";});
        // std::cout << "\n";
        // std::printf("Col ptr:\n");
        // std::for_each(col, col+elements_size, [&](const auto& ele){std::cout << ele << ", ";});
        // std::cout << "\n";
 
        // clear temprorary arrays
        tmp_row_ptr.clear();
        tmp_col_ptr.clear();
 
        // format data
        this->format();
 
        // run through crs data and simply write to the correct location, this can be done
        // in parallel
        FEAT_PRAGMA_OMP(parallel for)
        for(Index b_row = Index(0); b_row < n_rows; ++b_row)
        {
          for(Index l_row = Index(0); l_row < Index(BlockHeight); ++l_row)
          {
            IndexType b_idx = row_ptr[b_row];
            for(IT2_ idx = _crs_row[b_row*Index(BlockHeight) + l_row]; idx < _crs_row[b_row*Index(BlockHeight) + l_row + Index(1)]; ++idx)
            {
              //advance block row index until we find our current row
              while(col[b_idx] != IndexType(_crs_col[idx]/BlockWidth))
              {
                ++b_idx;
              }
              ASSERTM(col[b_idx] == IndexType(_crs_col[idx]/BlockWidth), "Blocked Column Array does not fit csr columns");
              // write data
              val[b_idx*BlockHeight*BlockWidth + BlockWidth*l_row + (_crs_col[idx]%BlockWidth)] = DT_(_crs_val[idx]);
            }
          }
        }
        //done
      }
 
      void convert(const Adjacency::Graph & graph)
      {
        this->move(SparseMatrixBCSR(graph));
      }
 
      SparseMatrixBCSR & operator= (const SparseLayout<IT_, layout_id> & layout_in)
      {
        for (Index i(0) ; i < this->_elements.size() ; ++i)
          MemoryPool::release_memory(this->_elements.at(i));
        for (Index i(0) ; i < this->_indices.size() ; ++i)
          MemoryPool::release_memory(this->_indices.at(i));
 
        this->_elements.clear();
        this->_indices.clear();
        this->_elements_size.clear();
        this->_indices_size.clear();
        this->_scalar_index.clear();
 
        this->_indices.assign(layout_in._indices.begin(), layout_in._indices.end());
        this->_indices_size.assign(layout_in._indices_size.begin(), layout_in._indices_size.end());
        this->_scalar_index.assign(layout_in._scalar_index.begin(), layout_in._scalar_index.end());
 
        for (auto i : this->_indices)
          MemoryPool::increase_memory(i);
 
        this->_elements.push_back(MemoryPool::template allocate_memory<DT_>(used_elements<Perspective::pod>()));
        this->_elements_size.push_back(used_elements<Perspective::pod>());
 
        return *this;
      }
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      void deserialize(std::vector<char> input)
      {
        this->template _deserialize<DT2_, IT2_>(FileMode::fm_bcsr, input);
      }
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      std::vector<char> serialize(const LAFEM::SerialConfig& config = SerialConfig()) const
      {
        return this->template _serialize<DT2_, IT2_>(FileMode::fm_bcsr, config);
      }
 
      void read_from(FileMode mode, const String& filename)
      {
        std::ios_base::openmode bin = std::ifstream::in | std::ifstream::binary;
        if(mode == FileMode::fm_mtx)
          bin = std::ifstream::in;
        std::ifstream file(filename.c_str(), bin);
        if (! file.is_open())
          XABORTM("Unable to open Matrix file " + filename);
        read_from(mode, file);
        file.close();
      }
      void read_from(FileMode mode, std::istream& file)
      {
        switch(mode)
        {
          /*case FileMode::fm_mtx:
           *   read_from_mtx(filename);
           *   break;*/
          case FileMode::fm_bcsr:
          case FileMode::fm_binary:
            this->template _deserialize<double, std::uint64_t>(FileMode::fm_bcsr, file);
            break;
          default:
            XABORTM("Filemode not supported!");
        }
      }
 
      /*
       * \brief Read in matrix from MatrixMarket mtx file.
       *
       * \param[in] filename The file that shall be read in.
       */
      /*void read_from_mtx(String filename)
      {
        std::ifstream file(filename.c_str(), std::ifstream::in);
        if (! file.is_open())
          XABORTM("Unable to open Matrix file " + filename);
        read_from_mtx(file);
        file.close();
      }*/
 
      /*
       * \brief Read in matrix from MatrixMarket mtx stream.
       *
       * \param[in] file The stream that shall be read in.
       */
      /*void read_from_mtx(std::istream& file)
      {
        this->clear();
        this->_scalar_index.push_back(0);
        this->_scalar_index.push_back(0);
        this->_scalar_index.push_back(0);
        this->_scalar_index.push_back(0);
 
        itd::map<IT_, std::map<IT_, DT_> > entries; // map<row, map<column, value> >
 
        Index ue(0);
        String line;
        std::getline(file, line);
        const bool general((line.find("%%MatrixMarket matrix coordinate real general") != String::npos) ? true : false);
        const bool symmetric((line.find("%%MatrixMarket matrix coordinate real symmetric") != String::npos) ? true : false);
 
        if (symmetric == false && general == false)
        {
          XABORTM("Input-file is not a compatible mtx-file");
        }
 
        while(!file.eof())
        {
          std::getline(file,line);
          if (file.eof())
            XABORTM("Input-file is empty");
 
          String::size_type begin(line.find_first_not_of(" "));
          if (line.at(begin) != '%')
            break;
        }
        {
          String::size_type begin(line.find_first_not_of(" "));
          line.erase(0, begin);
          String::size_type end(line.find_first_of(" "));
          String srow(line, 0, end);
          Index row((Index)atol(srow.c_str()));
          line.erase(0, end);
 
          begin = line.find_first_not_of(" ");
          line.erase(0, begin);
          end = line.find_first_of(" ");
          String scol(line, 0, end);
          Index col((Index)atol(scol.c_str()));
          line.erase(0, end);
          _rows() = row;
          _columns() = col;
          _size() = this->rows() * this->columns();
        }
 
        while(!file.eof())
        {
          std::getline(file, line);
          if (file.eof())
            break;
 
          String::size_type begin(line.find_first_not_of(" "));
          line.erase(0, begin);
          String::size_type end(line.find_first_of(" "));
          String srow(line, 0, end);
          IT_ row((IT_)atol(srow.c_str()));
          --row;
          line.erase(0, end);
 
          begin = line.find_first_not_of(" ");
          line.erase(0, begin);
          end = line.find_first_of(" ");
          String scol(line, 0, end);
          IT_ col((IT_)atol(scol.c_str()));
          --col;
          line.erase(0, end);
 
          begin = line.find_first_not_of(" ");
          line.erase(0, begin);
          end = line.find_first_of(" ");
          String sval(line, 0, end);
          DT_ tval((DT_)atof(sval.c_str()));
 
          entries[IT_(row)].insert(std::pair<IT_, DT_>(col, tval));
          ++ue;
          if (symmetric == true && row != col)
          {
            entries[IT_(col)].insert(std::pair<IT_, DT_>(row, tval));
            ++ue;
          }
        }
        _size() = this->rows() * this->columns();
        _used_elements() = ue;
 
        DT_ * tval = new DT_[ue];
        IT_ * tcol_ind = new IT_[ue];
        IT_ * trow_ptr = new IT_[rows() + 1];
 
        IT_ idx(0);
        Index row_idx(0);
        for (auto row : entries)
        {
          trow_ptr[row_idx] = idx;
          for (auto col : row.second )
          {
            tcol_ind[idx] = col.first;
            tval[idx] = col.second;
            ++idx;
          }
          row.second.clear();
          ++row_idx;
        }
        trow_ptr[rows()] = IT_(ue);
        entries.clear();
 
        this->_elements.push_back(MemoryPool<Mem_>::template allocate_memory<DT_>(_used_elements()));
        this->_elements_size.push_back(_used_elements());
        this->_indices.push_back(MemoryPool<Mem_>::template allocate_memory<IT_>(_used_elements()));
        this->_indices_size.push_back(_used_elements());
        this->_indices.push_back(MemoryPool<Mem_>::template allocate_memory<IT_>(rows() + 1));
        this->_indices_size.push_back(rows() + 1);
 
        MemoryPool<Mem_>::template upload<DT_>(this->_elements.at(0), tval, _used_elements());
        MemoryPool<Mem_>::template upload<IT_>(this->_indices.at(0), tcol_ind, _used_elements());
        MemoryPool<Mem_>::template upload<IT_>(this->_indices.at(1), trow_ptr, rows() + 1);
 
        delete[] tval;
        delete[] tcol_ind;
        delete[] trow_ptr;
      }*/
      void write_out(FileMode mode, const String& filename) const
      {
        std::ios_base::openmode bin = std::ofstream::out | std::ofstream::binary;
        if(mode == FileMode::fm_mtx)
          bin = std::ofstream::out;
        std::ofstream file;
        char* buff = nullptr;
        if(mode == FileMode::fm_mtx)
        {
          buff = new char[LAFEM::FileOutStreamBufferSize];
          file.rdbuf()->pubsetbuf(buff, LAFEM::FileOutStreamBufferSize);
        }
        file.open(filename.c_str(), bin);
        if(! file.is_open())
          XABORTM("Unable to open Matrix file " + filename);
        write_out(mode, file);
        file.close();
        delete[] buff;
      }
 
      void write_out(FileMode mode, std::ostream& file) const
      {
        switch(mode)
        {
          case FileMode::fm_bcsr:
          case FileMode::fm_binary:
            this->template _serialize<double, std::uint64_t>(FileMode::fm_bcsr, file);
            break;
          case FileMode::fm_mtx:
          {
            file << "%%MatrixMarket matrix coordinate real general\n";
            file << this->template rows<Perspective::pod>() << " " << this->template columns<Perspective::pod>() << " " << this->template used_elements<Perspective::pod>() << "\n";
 
            for (Index row(0) ; row < rows() ; ++row)
            {
              const IT_ end(this->row_ptr()[row + 1]);
              for (IT_ i(this->row_ptr()[row]) ; i < end ; ++i)
              {
                auto block = this->val()[i];
                for (int y(0) ; y < BlockHeight_ ; ++y)
                {
                  for (int x(0) ; x < BlockWidth_ ; ++x)
                  {
                    file << ((int)row * BlockHeight_) + y + 1 << " " << ((int)this->col_ind()[i] * BlockWidth_) + x + 1 << " " << stringify_fp_sci(block[y][x]) << "\n";
                  }
                }
              }
            }
            break;
          }
          default:
            XABORTM("Filemode not supported!");
        }
      }
 
      ValueType operator()(Index row, Index col) const
      {
        ASSERT(row < rows());
        ASSERT(col < columns());
 
        MemoryPool::synchronize();
 
        for (Index i(this->row_ptr()[row]) ; i < this->row_ptr()[row+1] ; ++i)
        {
          if (this->col_ind()[i] == col)
          {
            return this->val<Perspective::native>()[i];
          }
          if (this->col_ind()[i] > col)
            break; //return zero element
        }
 
        return ValueType(0.);
      }
 
      SparseLayout<IT_, layout_id> layout() const
      {
        return SparseLayout<IT_, layout_id>(this->_indices, this->_indices_size, this->_scalar_index);
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index rows() const
      {
        if (perspective_ == Perspective::pod)
          return this->_scalar_index.at(1) * Index(BlockHeight_);
        else
          return this->_scalar_index.at(1);
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index columns() const
      {
        if (perspective_ == Perspective::pod)
          return this->_scalar_index.at(2) * Index(BlockWidth_);
        else
          return this->_scalar_index.at(2);
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index used_elements() const
      {
        if (perspective_ == Perspective::pod)
          return this->_scalar_index.at(3) * Index(BlockHeight_ * BlockWidth_);
        else
          return this->_scalar_index.at(3);
      }
 
      IT_ * col_ind()
      {
        if (this->size() == 0)
          return nullptr;
 
        return this->_indices.at(0);
      }
 
      IT_ const * col_ind() const
      {
        if (this->size() == 0)
          return nullptr;
 
        return this->_indices.at(0);
      }
 
      template <Perspective perspective_ = Perspective::native>
      auto val() const -> const typename Intern::BCSRPerspectiveHelper<DT_, BlockHeight_, BlockWidth_, perspective_>::Type *
      {
        if (this->size() == 0)
          return nullptr;
 
        return (const typename Intern::BCSRPerspectiveHelper<DT_, BlockHeight_, BlockWidth_, perspective_>::Type *)(this->_elements.at(0));
      }
 
      template <Perspective perspective_ = Perspective::native>
      auto val() -> typename Intern::BCSRPerspectiveHelper<DT_, BlockHeight_, BlockWidth_, perspective_>::Type *
      {
        if (this->size() == 0)
          return nullptr;
 
        return (typename Intern::BCSRPerspectiveHelper<DT_, BlockHeight_, BlockWidth_, perspective_>::Type *)(this->_elements.at(0));
      }
 
      IT_ * row_ptr()
      {
        if (this->_indices.size() == 0)
          return nullptr;
 
        return this->_indices.at(1);
      }
 
      IT_ const * row_ptr() const
      {
        if (this->_indices.size() == 0)
          return nullptr;
 
        return this->_indices.at(1);
      }
 
      static String name()
      {
        return "SparseMatrixBCSR";
      }
 
      void copy(const SparseMatrixBCSR & x, bool full = false)
      {
        this->_copy_content(x, full);
      }
 
 
      void axpy(
                const SparseMatrixBCSR & x,
                const DT_ alpha = DT_(1))
      {
        XASSERTM(x.rows() == this->rows(), "Matrix rows do not match!");
        XASSERTM(x.columns() == this->columns(), "Matrix columns do not match!");
        XASSERTM(x.used_elements() == this->used_elements(), "Matrix used_elements do not match!");
 
        TimeStamp ts_start;
 
        Statistics::add_flops(this->used_elements<Perspective::pod>() * 2);
        Arch::Axpy::value(this->template val<Perspective::pod>(),
            alpha,
            x.template val<Perspective::pod>(),
            this->used_elements<Perspective::pod>());
 
        TimeStamp ts_stop;
        Statistics::add_time_axpy(ts_stop.elapsed(ts_start));
      }
 
      void scale(const SparseMatrixBCSR & x, const DT_ alpha)
      {
        XASSERTM(x.rows() == this->rows(), "Row count does not match!");
        XASSERTM(x.columns() == this->columns(), "Column count does not match!");
        XASSERTM(x.used_elements() == this->used_elements(), "Nonzero count does not match!");
 
        TimeStamp ts_start;
        Statistics::add_flops(this->used_elements<Perspective::pod>());
        Arch::Scale::value(this->template val<Perspective::pod>(), x.template val<Perspective::pod>(), alpha, this->used_elements<Perspective::pod>());
        TimeStamp ts_stop;
        Statistics::add_time_axpy(ts_stop.elapsed(ts_start));
      }
 
      DT_ norm_frobenius() const
      {
        TimeStamp ts_start;
        Statistics::add_flops(this->used_elements<Perspective::pod>() * 2);
        DT_ result = Arch::Norm2::value(this->template val<Perspective::pod>(),
                                              this->used_elements<Perspective::pod>());
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
        return result;
      }
 
      void row_norm2(VectorTypeL& row_norms) const
      {
        XASSERTM(row_norms.size() == this->rows(), "Matrix/Vector dimension mismatch");
 
        TimeStamp ts_start;
        Statistics::add_flops(this->used_elements<Perspective::pod>() * 2);
 
        Arch::RowNorm::bcsr_norm2(row_norms.template elements<Perspective::pod>(),
          this->template val<Perspective::pod>(),
          col_ind(), row_ptr(), rows(), BlockHeight, BlockWidth);
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
      }
 
      void row_norm2sqr(VectorTypeL& row_norms) const
      {
        XASSERTM(row_norms.size() == this->rows(), "Matrix/Vector dimension mismatch");
 
        TimeStamp ts_start;
        Statistics::add_flops(this->used_elements<Perspective::pod>() * 2);
 
        Arch::RowNorm::bcsr_norm2sqr(row_norms.template elements<Perspective::pod>(),
          this->template val<Perspective::pod>(),
          col_ind(), row_ptr(), rows(), BlockHeight, BlockWidth);
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
      }
 
      void row_norm2sqr(VectorTypeL& row_norms, const VectorTypeR& scal) const
      {
        XASSERTM(row_norms.size() == this->rows(), "Matrix/Vector dimension mismatch");
        XASSERTM(scal.size() == this->columns(), "Matrix/scalings dimension mismatch");
 
        TimeStamp ts_start;
        Statistics::add_flops(this->used_elements<Perspective::pod>() * 2);
 
        Arch::RowNorm::bcsr_scaled_norm2sqr(row_norms.template elements<Perspective::pod>(),
           scal.template elements<Perspective::pod>(), this->template val<Perspective::pod>(),
          col_ind(), row_ptr(), rows(), BlockHeight, BlockWidth);
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
      }
 
      DT_ max_abs_element() const
      {
        TimeStamp ts_start;
 
        Index max_abs_index = Arch::MaxAbsIndex::value(this->template val<Perspective::pod>(), this->template used_elements<Perspective::pod>());
        ASSERT(max_abs_index < this->template used_elements<Perspective::pod>());
        DT_ result;
        MemoryPool::template copy<DT_>(&result, this->template val<Perspective::pod>() + max_abs_index, 1);
        result = Math::abs(result);
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      DT_ min_abs_element() const
      {
        TimeStamp ts_start;
 
        Index min_abs_index = Arch::MinAbsIndex::value(this->template val<Perspective::pod>(), this->template used_elements<Perspective::pod>());
        ASSERT(min_abs_index < this->template used_elements<Perspective::pod>());
        DT_ result;
        MemoryPool::template copy<DT_>(&result, this->template val<Perspective::pod>() + min_abs_index, 1);
        result = Math::abs(result);
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      DT_ max_element() const
      {
        TimeStamp ts_start;
 
        Index max_index = Arch::MaxIndex::value(this->template val<Perspective::pod>(), this->template used_elements<Perspective::pod>());
        ASSERT(max_index < this->template used_elements<Perspective::pod>());
        DT_ result;
        MemoryPool::template copy<DT_>(&result, this->template val<Perspective::pod>() + max_index, 1);
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      DT_ min_element() const
      {
        TimeStamp ts_start;
 
        Index min_index = Arch::MinIndex::value(this->template val<Perspective::pod>(), this->template used_elements<Perspective::pod>());
        ASSERT(min_index < this->template used_elements<Perspective::pod>());
        DT_ result;
        MemoryPool::template copy<DT_>(&result, this->template val<Perspective::pod>() + min_index, 1);
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      DT_ max_rel_diff(const SparseMatrixBCSR& x) const
      {
        XASSERTM(x.used_elements() == this->used_elements(), "Nonzero count does not match!");
        TimeStamp ts_start;
 
        DataType max_rel_diff = Arch::MaxRelDiff::value(this->template val<Perspective::pod>(), x.template val<Perspective::pod>(), this->template used_elements<Perspective::pod>());
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return max_rel_diff;
      }
 
      bool same_layout(const SparseMatrixBCSR& x) const
      {
        if (this->row_ptr ()== x.row_ptr() && this->col_ind() == x.col_ind())
          return true;
        if(this->size() == 0 && x.size() == 0 && this->get_elements().size() == 0 && this->get_indices().size() == 0 && x.get_elements().size() == 0 && x.get_indices().size() == 0)
          return true;
        if (this->rows() != x.rows())
          return false;
        if (this->columns() != x.columns())
          return false;
        if (this->used_elements() != x.used_elements())
          return false;
 
        IT_ * col_ind_a;
        IT_ * col_ind_b;
        IT_ * row_ptr_a;
        IT_ * row_ptr_b;
 
        col_ind_a = const_cast<IT_*>(this->col_ind());
        row_ptr_a = const_cast<IT_*>(this->row_ptr());
        col_ind_b = const_cast<IT_*>(x.col_ind());
        row_ptr_b = const_cast<IT_*>(x.row_ptr());
 
        bool ret(true);
        for (Index i(0) ; i < this->used_elements() ; ++i)
        {
          if (col_ind_a[i] != col_ind_b[i])
          {
            ret = false;
            break;
          }
        }
        if (ret)
        {
          for (Index i(0) ; i < this->rows() + 1; ++i)
          {
            if (row_ptr_a[i] != row_ptr_b[i])
            {
              ret = false;
              break;
            }
          }
        }
 
        return ret;
      }
 
      SparseMatrixBCSR<DT_, IT_, BlockWidth_, BlockHeight_> transpose() const
      {
        SparseMatrixBCSR<DT_, IT_, BlockWidth_, BlockHeight_> x_t;
        x_t.transpose(*this);
        return x_t;
      }
 
      void transpose(const SparseMatrixBCSR<DT_, IT_, BlockWidth_, BlockHeight_> & x)
      {
        using XType = SparseMatrixBCSR<DT_, IT_, BlockWidth_, BlockHeight_>;
        if (x.used_elements() == 0)
        {
          this->move(SparseMatrixBCSR(x.rows(), x.columns()));
          return;
        }
 
        const Index txrows(x.rows());
        const Index txcolumns(x.columns());
        const Index txused_elements(x.used_elements());
 
        const IT_ * ptxcol_ind(x.col_ind());
        const IT_ * ptxrow_ptr(x.row_ptr());
        const typename XType::ValueType * ptxval(x.val());
 
        DenseVector<IT_, IT_> tcol_ind(txused_elements);
        DenseVector<DT_, IT_> tval(txused_elements * BlockHeight_ * BlockWidth_);
        DenseVector<IT_, IT_> trow_ptr(txcolumns + 1, IT_(0));
 
        IT_ * ptcol_ind(tcol_ind.elements());
        ValueType * ptval((ValueType*)tval.elements());
        IT_ * ptrow_ptr(trow_ptr.elements());
 
        ptrow_ptr[0] = 0;
 
        for (Index i(0); i < txused_elements; ++i)
        {
          ++ptrow_ptr[ptxcol_ind[i] + 1];
        }
 
        for (Index i(1); i < txcolumns - 1; ++i)
        {
          ptrow_ptr[i + 1] += ptrow_ptr[i];
        }
 
        for (Index i(0); i < txrows; ++i)
        {
          for (IT_ k(ptxrow_ptr[i]); k < ptxrow_ptr[i+1]; ++k)
          {
            const IT_ l(ptxcol_ind[k]);
            const IT_ j(ptrow_ptr[l]);
            ptval[j].set_transpose(ptxval[k]);
            ptcol_ind[j] = IT_(i);
            ++ptrow_ptr[l];
          }
        }
 
        for (Index i(txcolumns); i > 0; --i)
        {
          ptrow_ptr[i] = ptrow_ptr[i - 1];
        }
        ptrow_ptr[0] = 0;
 
        this->move(SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_> (txcolumns, txrows, tcol_ind, tval, trow_ptr));
      }
 
      void scale_rows(const SparseMatrixBCSR & x, const DenseVectorBlocked<DT_,IT_, BlockHeight_> & s)
      {
        XASSERTM(x.rows() == this->rows(), "Row count does not match!");
        XASSERTM(x.columns() == this->columns(), "Column count does not match!");
        XASSERTM(x.used_elements() == this->used_elements(), "Nonzero count does not match!");
        XASSERTM(s.size() == this->rows(), "Vector size does not match!");
 
        TimeStamp ts_start;
 
        Statistics::add_flops(this->used_elements<Perspective::pod>());
        Arch::ScaleRows::template bcsr<BlockHeight_, BlockWidth_>(this->val<Perspective::pod>(), x.template val<Perspective::pod>(),
          this->col_ind(), this->row_ptr(), s.template elements<Perspective::pod>(), this->rows(), this->columns(), this->used_elements());
 
        TimeStamp ts_stop;
        Statistics::add_time_axpy(ts_stop.elapsed(ts_start));
      }
 
      void scale_cols(const SparseMatrixBCSR & x, const DenseVectorBlocked<DT_,IT_, BlockWidth_> & s)
      {
        XASSERTM(x.rows() == this->rows(), "Row count does not match!");
        XASSERTM(x.columns() == this->columns(), "Column count does not match!");
        XASSERTM(x.used_elements() == this->used_elements(), "Nonzero count does not match!");
        XASSERTM(s.size() == this->columns(), "Vector size does not match!");
 
        TimeStamp ts_start;
 
        Statistics::add_flops(this->used_elements<Perspective::pod>());
        Arch::ScaleCols::template bcsr<BlockHeight_, BlockWidth_>(this->val<Perspective::pod>(), x.template val<Perspective::pod>(),
          this->col_ind(), this->row_ptr(), s.template elements<Perspective::pod>(), this->rows(), this->columns(), this->used_elements());
 
        TimeStamp ts_stop;
        Statistics::add_time_axpy(ts_stop.elapsed(ts_start));
      }
 
      void apply(DenseVector<DT_, IT_> & r, const DenseVector< DT_, IT_> & x) const
      {
        XASSERTM(r.size() == this->rows<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns<Perspective::pod>(), "Vector size of x does not match!");
 
        TimeStamp ts_start;
        FEAT_KERNEL_MARKER_START("BCSR-apply");
 
        if (this->used_elements() == 0)
        {
          r.format();
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops(this->used_elements<Perspective::pod>() * 2);
        Arch::Apply::template bcsr<BlockHeight_, BlockWidth_>(
            r.elements(), DT_(1), x.elements(), DT_(0), r.elements(), this->template val<Perspective::pod>(), this->col_ind(),
            this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
      }
 
      void apply_transposed(DenseVector<DT_, IT_> & r, const DenseVector< DT_, IT_> & x) const
      {
        XASSERTM(r.size() == this->columns<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows<Perspective::pod>(), "Vector size of x does not match!");
 
        TimeStamp ts_start;
        FEAT_KERNEL_MARKER_START("BCSR-apply");
 
        if (this->used_elements() == 0)
        {
          r.format();
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops(this->used_elements<Perspective::pod>() * 2);
        Arch::Apply::template bcsr_transposed<BlockHeight_, BlockWidth_>(
          r.elements(), DT_(1), x.elements(), DT_(0), r.elements(), this->template val<Perspective::pod>(), this->col_ind(),
          this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
      }
 
      void apply(DenseVectorBlocked<DT_, IT_, BlockHeight_> & r, const DenseVector< DT_, IT_> & x) const
      {
        XASSERTM(r.size() == this->rows(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns<Perspective::pod>(), "Vector size of x does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0)
        {
          r.format();
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops(this->used_elements<Perspective::pod>() * 2);
 
        Arch::Apply::template bcsr<BlockHeight_, BlockWidth_>(
            r.template elements<Perspective::pod>(), DT_(1), x.elements(), DT_(0), r.template elements<Perspective::pod>(), this->template val<Perspective::pod>(), this->col_ind(),
            this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
 
      void apply_transposed(DenseVectorBlocked<DT_, IT_, BlockWidth_> & r, const DenseVector< DT_, IT_> & x) const
      {
        XASSERTM(r.size() == this->columns(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows<Perspective::pod>(), "Vector size of x does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0)
        {
          r.format();
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops(this->used_elements<Perspective::pod>() * 2);
 
        Arch::Apply::template bcsr_transposed<BlockHeight_, BlockWidth_>(
          r.template elements<Perspective::pod>(), DT_(1), x.elements(), DT_(0), r.template elements<Perspective::pod>(), this->template val<Perspective::pod>(), this->col_ind(),
          this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply(DenseVector<DT_, IT_> & r, const DenseVectorBlocked<DT_, IT_, BlockWidth_> & x) const
      {
        XASSERTM(r.size() == this->rows<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns(), "Vector size of x does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0)
        {
          r.format();
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops(this->used_elements<Perspective::pod>() * 2);
 
        Arch::Apply::template bcsr<BlockHeight_, BlockWidth_>(
            r.elements(), DT_(1), x.template elements<Perspective::pod>(), DT_(0), r.elements(), this->template val<Perspective::pod>(), this->col_ind(),
            this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply_transposed(DenseVector<DT_, IT_> & r, const DenseVectorBlocked<DT_, IT_, BlockHeight_> & x) const
      {
        XASSERTM(r.size() == this->columns<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows(), "Vector size of x does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0)
        {
          r.format();
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops(this->used_elements<Perspective::pod>() * 2);
 
        Arch::Apply::template bcsr_transposed<BlockHeight_, BlockWidth_>(
          r.elements(), DT_(1), x.template elements<Perspective::pod>(), DT_(0), r.elements(), this->template val<Perspective::pod>(), this->col_ind(),
          this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply(DenseVectorBlocked<DT_, IT_, BlockHeight_> & r, const DenseVectorBlocked<DT_, IT_, BlockWidth_> & x) const
      {
        XASSERTM(r.size() == this->rows(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns(), "Vector size of x does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0)
        {
          r.format();
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops(this->used_elements<Perspective::pod>() * 2);
 
        Arch::Apply::template bcsr<BlockHeight_, BlockWidth_>(
            r.template elements<Perspective::pod>(), DT_(1), x.template elements<Perspective::pod>(), DT_(0), r.template elements<Perspective::pod>(), this->template val<Perspective::pod>(),
            this->col_ind(), this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply_transposed(DenseVectorBlocked<DT_, IT_, BlockWidth_> & r, const DenseVectorBlocked<DT_, IT_, BlockHeight_> & x) const
      {
        XASSERTM(r.size() == this->columns(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows(), "Vector size of x does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0)
        {
          r.format();
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops(this->used_elements<Perspective::pod>() * 2);
 
        Arch::Apply::template bcsr_transposed<BlockHeight_, BlockWidth_>(
          r.template elements<Perspective::pod>(), DT_(1), x.template elements<Perspective::pod>(), DT_(0), r.template elements<Perspective::pod>(), this->template val<Perspective::pod>(),
          this->col_ind(), this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply(
                 DenseVector<DT_, IT_> & r,
                 const DenseVector< DT_, IT_> & x,
                 const DenseVector< DT_, IT_> & y,
                 const DT_ alpha = DT_(1)) const
      {
        XASSERTM(r.size() == this->rows<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns<Perspective::pod>(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->rows<Perspective::pod>(), "Vector size of y does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0 || Math::abs(alpha) < Math::eps<DT_>())
        {
          r.copy(y);
          //r.scale(beta);
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops( (this->used_elements<Perspective::pod>() + this->rows<Perspective::pod>()) * 2 );
 
        Arch::Apply::template bcsr<BlockHeight_, BlockWidth_>(
            r.elements(), alpha, x.elements(), DT_(1), y.elements(), this->template val<Perspective::pod>(), this->col_ind(),
            this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply_transposed(
        DenseVector<DT_, IT_> & r,
        const DenseVector< DT_, IT_> & x,
        const DenseVector< DT_, IT_> & y,
        const DT_ alpha = DT_(1)) const
      {
        XASSERTM(r.size() == this->columns<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows<Perspective::pod>(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->columns<Perspective::pod>(), "Vector size of y does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0 || Math::abs(alpha) < Math::eps<DT_>())
        {
          r.copy(y);
          //r.scale(beta);
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops( (this->used_elements<Perspective::pod>() + this->rows<Perspective::pod>()) * 2 );
 
        Arch::Apply::template bcsr_transposed<BlockHeight_, BlockWidth_>(
          r.elements(), alpha, x.elements(), DT_(1), y.elements(), this->template val<Perspective::pod>(), this->col_ind(),
          this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply(
                 DenseVectorBlocked<DT_, IT_, BlockHeight_> & r,
                 const DenseVector< DT_, IT_> & x,
                 const DenseVectorBlocked<DT_, IT_, BlockHeight_> & y,
                 const DT_ alpha = DT_(1)) const
      {
        XASSERTM(r.size() == this->rows(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns<Perspective::pod>(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->rows(), "Vector size of y does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0 || Math::abs(alpha) < Math::eps<DT_>())
        {
          r.copy(y);
          //r.scale(beta);
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops( (this->used_elements<Perspective::pod>() + this->rows<Perspective::pod>()) * 2 );
        Arch::Apply::template bcsr<BlockHeight_, BlockWidth_>(
            r.template elements<Perspective::pod>(), alpha, x.elements(), DT_(1), y.template elements<Perspective::pod>(), this->template val<Perspective::pod>(), this->col_ind(),
            this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply_transposed(
        DenseVectorBlocked<DT_, IT_, BlockWidth_> & r,
        const DenseVector< DT_, IT_> & x,
        const DenseVectorBlocked<DT_, IT_, BlockWidth_> & y,
        const DT_ alpha = DT_(1)) const
      {
        XASSERTM(r.size() == this->columns(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows<Perspective::pod>(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->columns(), "Vector size of y does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if(this->used_elements() == 0 || Math::abs(alpha) < Math::eps<DT_>())
        {
          r.copy(y);
          //r.scale(beta);
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops((this->used_elements<Perspective::pod>() + this->rows<Perspective::pod>()) * 2);
        Arch::Apply::template bcsr_transposed<BlockHeight_, BlockWidth_>(
          r.template elements<Perspective::pod>(), alpha, x.elements(), DT_(1), y.template elements<Perspective::pod>(), this->template val<Perspective::pod>(), this->col_ind(),
          this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply(
                 DenseVector<DT_, IT_> & r,
                 const DenseVectorBlocked<DT_, IT_, BlockWidth_> & x,
                 const DenseVector< DT_, IT_> & y,
                 const DT_ alpha = DT_(1)) const
      {
        XASSERTM(r.size() == this->rows<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->rows<Perspective::pod>(), "Vector size of y does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0 || Math::abs(alpha) < Math::eps<DT_>())
        {
          r.copy(y);
          //r.scale(beta);
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops( (this->used_elements<Perspective::pod>() + this->rows<Perspective::pod>()) * 2 );
        Arch::Apply::template bcsr<BlockHeight_, BlockWidth_>(
            r.elements(), alpha, x.template elements<Perspective::pod>(), DT_(1), y.elements(), this->template val<Perspective::pod>(), this->col_ind(),
            this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply_transposed(
        DenseVector<DT_, IT_> & r,
        const DenseVectorBlocked<DT_, IT_, BlockHeight_> & x,
        const DenseVector< DT_, IT_> & y,
        const DT_ alpha = DT_(1)) const
      {
        XASSERTM(r.size() == this->columns<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->columns<Perspective::pod>(), "Vector size of y does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0 || Math::abs(alpha) < Math::eps<DT_>())
        {
          r.copy(y);
          //r.scale(beta);
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops( (this->used_elements<Perspective::pod>() + this->rows<Perspective::pod>()) * 2 );
        Arch::Apply::template bcsr_transposed<BlockHeight_, BlockWidth_>(
          r.elements(), alpha, x.template elements<Perspective::pod>(), DT_(1), y.elements(), this->template val<Perspective::pod>(), this->col_ind(),
          this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply(
                 DenseVectorBlocked<DT_, IT_, BlockHeight_> & r,
                 const DenseVectorBlocked<DT_, IT_, BlockWidth_> & x,
                 const DenseVectorBlocked<DT_, IT_, BlockHeight_> & y,
                 const DT_ alpha = DT_(1)) const
      {
        XASSERTM(r.size() == this->rows(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->rows(), "Vector size of y does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0 || Math::abs(alpha) < Math::eps<DT_>())
        {
          r.copy(y);
          //r.scale(beta);
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops( (this->used_elements<Perspective::pod>() + this->rows<Perspective::pod>()) * 2 );
        Arch::Apply::template bcsr<BlockHeight_, BlockWidth_>(
            r.template elements<Perspective::pod>(), alpha, x.template elements<Perspective::pod>(), DT_(1), y.template elements<Perspective::pod>(), this->template val<Perspective::pod>(),
            this->col_ind(), this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply_transposed(
        DenseVectorBlocked<DT_, IT_, BlockWidth_> & r,
        const DenseVectorBlocked<DT_, IT_, BlockHeight_> & x,
        const DenseVectorBlocked<DT_, IT_, BlockWidth_> & y,
        const DT_ alpha = DT_(1)) const
      {
        XASSERTM(r.size() == this->columns(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->columns(), "Vector size of y does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0 || Math::abs(alpha) < Math::eps<DT_>())
        {
          r.copy(y);
          //r.scale(beta);
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops( (this->used_elements<Perspective::pod>() + this->rows<Perspective::pod>()) * 2 );
        Arch::Apply::template bcsr_transposed<BlockHeight_, BlockWidth_>(
          r.template elements<Perspective::pod>(), alpha, x.template elements<Perspective::pod>(), DT_(1), y.template elements<Perspective::pod>(), this->template val<Perspective::pod>(),
          this->col_ind(), this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply(
                 DenseVectorBlocked<DT_, IT_, BlockHeight_> & r,
                 const DenseVectorBlocked<DT_, IT_, BlockWidth_> & x,
                 const DenseVector< DT_, IT_> & y,
                 const DT_ alpha = DT_(1)) const
      {
        XASSERTM(r.size() == this->rows(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->rows<Perspective::pod>(), "Vector size of y does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0 || Math::abs(alpha) < Math::eps<DT_>())
        {
          r.convert(y);
          //r.scale(beta);
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops( (this->used_elements<Perspective::pod>() + this->rows<Perspective::pod>()) * 2 );
        Arch::Apply::template bcsr<BlockHeight_, BlockWidth_>(
            r.template elements<Perspective::pod>(), alpha, x.template elements<Perspective::pod>(), DT_(1), y.template elements<Perspective::pod>(), this->template val<Perspective::pod>(),
            this->col_ind(), this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void apply_transposed(
        DenseVectorBlocked<DT_, IT_, BlockWidth_> & r,
        const DenseVectorBlocked<DT_, IT_, BlockHeight_> & x,
        const DenseVector< DT_, IT_> & y,
        const DT_ alpha = DT_(1)) const
      {
        XASSERTM(r.size() == this->columns(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->columns<Perspective::pod>(), "Vector size of y does not match!");
 
        TimeStamp ts_start;
 
        FEAT_KERNEL_MARKER_START("BCSR-apply");
        if (this->used_elements() == 0 || Math::abs(alpha) < Math::eps<DT_>())
        {
          r.convert(y);
          //r.scale(beta);
          return;
        }
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        Statistics::add_flops( (this->used_elements<Perspective::pod>() + this->rows<Perspective::pod>()) * 2 );
        Arch::Apply::template bcsr_transposed<BlockHeight_, BlockWidth_>(
          r.template elements<Perspective::pod>(), alpha, x.template elements<Perspective::pod>(), DT_(1), y.template elements<Perspective::pod>(), this->template val<Perspective::pod>(),
          this->col_ind(), this->row_ptr(), this->rows(), this->columns(), this->used_elements());
 
        FEAT_KERNEL_MARKER_STOP("BCSR-apply");
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void add_double_mat_product(
        const LAFEM::SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>& d,
        const LAFEM::SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>& a,
        const LAFEM::SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>& b,
        const DT_ alpha = DT_(1),
        const bool allow_incomplete = false)
      {
        // validate matrix dimensions
        XASSERT(BlockHeight_ == BlockWidth_);
        XASSERT(this->rows() == d.rows());
        XASSERT(d.columns() == a.rows());
        XASSERT(a.columns() == b.rows());
        XASSERT(b.columns() == this->columns());
 
        // fetch matrix arrays:
        ValueType* data_x = this->val();
        const ValueType* data_d = d.val();
        const ValueType* data_a = a.val();
        const ValueType* data_b = b.val();
        const IT_* row_ptr_x = this->row_ptr();
        const IT_* col_idx_x = this->col_ind();
        const IT_* row_ptr_d = d.row_ptr();
        const IT_* col_idx_d = d.col_ind();
        const IT_* row_ptr_a = a.row_ptr();
        const IT_* col_idx_a = a.col_ind();
        const IT_* row_ptr_b = b.row_ptr();
        const IT_* col_idx_b = b.col_ind();
 
        // loop over all rows of D and X, resp.
        for(IT_ i(0); i < IT_(this->rows()); ++i)
        {
          // loop over all non-zeros D_ik in row i of D
          for(IT_ ik(row_ptr_d[i]); ik  < row_ptr_d[i+1]; ++ik)
          {
            // get column index k
            const IT_ k = col_idx_d[ik];
 
            // loop over all non-zeros A_kl in row k of A
            for(IT_ kl(row_ptr_a[k]); kl < row_ptr_a[k+1]; ++kl)
            {
              // get column index l
              const IT_ l = col_idx_a[kl];
 
              // pre-compute factor (alpha * D_ik * A_kl)
              ValueType omega;
              omega.set_mat_mat_mult(data_d[ik], data_a[kl]);
              omega *= alpha;
 
              // loop over all non-zeros B_lj in row j of B and
              // loop over all non-zeros X_ij in row i of X and
              // perform a "sparse axpy" of B_l onto X_i, i.e.:
              //   X_i. += (alpha * D_ik * A_kl) * B_l.
              IT_ ij = row_ptr_x[i];
              IT_ lj = row_ptr_b[l];
              while(lj < row_ptr_b[l+1])
              {
                if(ij >= row_ptr_x[i+1])
                {
                  // we have reached the end of row X_i, but there is at least
                  // one entry in row B_l left, so the pattern of X is incomplete
                  // We let the caller decide whether this is a valid case or not:
                  if(allow_incomplete)
                    break;  // continue with next row
                  else
                    XABORTM("Incomplete output matrix structure");
                }
                else if(col_idx_x[ij] == col_idx_b[lj])
                {
                  // okay: B_lj contributes to X_ij
                  ValueType temp;
                  temp.set_mat_mat_mult(omega, data_b[lj]);
                  data_x[ij] += temp;
                  ++ij;
                  ++lj;
                }
                else if(col_idx_x[ij] < col_idx_b[lj])
                {
                  // entry X_ij exists, but B_lj is missing:
                  // this is a perfectly valid case, so continue with the next non-zero of X_i
                  ++ij;
                }
                else //if(col_idx_x[ij] > col_idx_b[lj])
                {
                  // If we come out here, then the sparsity pattern of X is incomplete:
                  // B_lj is meant to be added onto X_ij, but the entry X_ij is missing
                  // We let the caller decide whether this is a valid case or not:
                  if(allow_incomplete)
                    ++lj;
                  else
                    XABORTM("Incomplete output matrix structure");
                }
              }
            }
          }
        }
      }
 
      void add_double_mat_product(
        const LAFEM::SparseMatrixCSR<DT_, IT_>& d,
        const LAFEM::SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>& a,
        const LAFEM::SparseMatrixCSR<DT_, IT_>& b,
        const DT_ alpha = DT_(1),
        const bool allow_incomplete = false)
      {
        // validate matrix dimensions
        XASSERT(BlockHeight_ == BlockWidth_);
        XASSERT(this->rows() == d.rows());
        XASSERT(d.columns() == a.rows());
        XASSERT(a.columns() == b.rows());
        XASSERT(b.columns() == this->columns());
 
        // fetch matrix arrays:
        ValueType* data_x = this->val();
        const DT_* data_d = d.val();
        const ValueType* data_a = a.val();
        const DT_* data_b = b.val();
        const IT_* row_ptr_x = this->row_ptr();
        const IT_* col_idx_x = this->col_ind();
        const IT_* row_ptr_d = d.row_ptr();
        const IT_* col_idx_d = d.col_ind();
        const IT_* row_ptr_a = a.row_ptr();
        const IT_* col_idx_a = a.col_ind();
        const IT_* row_ptr_b = b.row_ptr();
        const IT_* col_idx_b = b.col_ind();
 
        // loop over all rows of D and X, resp.
        for(IT_ i(0); i < IT_(this->rows()); ++i)
        {
          // loop over all non-zeros D_ik in row i of D
          for(IT_ ik(row_ptr_d[i]); ik  < row_ptr_d[i+1]; ++ik)
          {
            // get column index k
            const IT_ k = col_idx_d[ik];
 
            // loop over all non-zeros A_kl in row k of A
            for(IT_ kl(row_ptr_a[k]); kl < row_ptr_a[k+1]; ++kl)
            {
              // get column index l
              const IT_ l = col_idx_a[kl];
 
              // pre-compute factor (alpha * D_ik * A_kl)
              ValueType omega = (alpha * data_d[ik]) * data_a[kl];
 
              // loop over all non-zeros B_lj in row j of B and
              // loop over all non-zeros X_ij in row i of X and
              // perform a "sparse axpy" of B_l onto X_i, i.e.:
              //   X_i. += (alpha * D_ik * A_kl) * B_l.
              IT_ ij = row_ptr_x[i];
              IT_ lj = row_ptr_b[l];
              while(lj < row_ptr_b[l+1])
              {
                if(ij >= row_ptr_x[i+1])
                {
                  // we have reached the end of row X_i, but there is at least
                  // one entry in row B_l left, so the pattern of X is incomplete
                  // We let the caller decide whether this is a valid case or not:
                  if(allow_incomplete)
                    break;  // continue with next row
                  else
                    XABORTM("Incomplete output matrix structure");
                }
                else if(col_idx_x[ij] == col_idx_b[lj])
                {
                  // okay: B_lj contributes to X_ij
                  Tiny::axpy(data_x[ij], omega, data_b[lj]);
                  ++ij;
                  ++lj;
                }
                else if(col_idx_x[ij] < col_idx_b[lj])
                {
                  // entry X_ij exists, but B_lj is missing:
                  // this is a perfectly valid case, so continue with the next non-zero of X_i
                  ++ij;
                }
                else //if(col_idx_x[ij] > col_idx_b[lj])
                {
                  // If we come out here, then the sparsity pattern of X is incomplete:
                  // B_lj is meant to be added onto X_ij, but the entry X_ij is missing
                  // We let the caller decide whether this is a valid case or not:
                  if(allow_incomplete)
                    ++lj;
                  else
                    XABORTM("Incomplete output matrix structure");
                }
              }
            }
          }
        }
      }
 
      void lump_rows(VectorTypeL& lump) const
      {
        XASSERTM(lump.size() == rows(), "lump vector size does not match matrix row count!");
 
        Arch::Lumping::bcsr(
          lump.template elements<Perspective::pod>(),
          this->template val<Perspective::pod>(),
          col_ind(), row_ptr(), rows(), BlockHeight, BlockWidth);
      }
 
      VectorTypeL lump_rows() const
      {
        VectorTypeL lump = create_vector_l();
        lump_rows(lump);
        return lump;
      }
 
 
      void extract_diag(VectorTypeL & diag, DenseVector<IT_, IT_> & diag_indices) const
      {
        XASSERTM(diag.size() == rows(), "diag size does not match matrix row count!");
        XASSERTM(rows() == columns(), "matrix is not square!");
 
        for (Index row(0); row < rows(); row++)
        {
          const Index index = diag_indices.elements()[row];
          typename VectorTypeL::ValueType t(0);
          if (index != used_elements())
          {
            ValueType m = this->template val<LAFEM::Perspective::native>()[index];
            for (int i(0) ; i < BlockHeight_ ; ++i)
            {
              t[i] = m[i][i];
            }
          }
          diag(row, t);
        }
      }
 
      void extract_diag(VectorTypeL & diag) const
      {
        auto diag_indices = extract_diag_indices();
        extract_diag(diag, diag_indices);
      }
 
 
      VectorTypeL extract_diag() const
      {
        VectorTypeL diag = create_vector_l();
        extract_diag(diag);
        return diag;
      }
 
      void extract_diag_indices(DenseVector<IT_, IT_> & diag_indices) const
      {
        XASSERTM(diag_indices.size() == rows(), "diag size does not match matrix row count!");
        XASSERTM(rows() == columns(), "matrix is not square!");
 
        Arch::Diagonal::csr(diag_indices.elements(), col_ind(), row_ptr(), rows());
      }
 
      DenseVector<IT_, IT_> extract_diag_indices() const
      {
        DenseVector<IT_, IT_> diag_indices(this->template rows<LAFEM::Perspective::native>());
        extract_diag_indices(diag_indices);
        return diag_indices;
      }
 
      void permute(Adjacency::Permutation & perm_row, Adjacency::Permutation & perm_col)
      {
        if (perm_row.size() == 0 && perm_col.size() == 0)
          return;
 
        XASSERTM(perm_row.size() == this->rows(), "Container rows does not match permutation size");
        XASSERTM(perm_col.size() == this->columns(), "Container columns does not match permutation size");
 
        // http://de.mathworks.com/help/matlab/math/sparse-matrix-operations.html#f6-13070
        IT_ * temp_row_ptr = new IT_[rows() + 1];
        IT_ * temp_col_ind = new IT_[used_elements()];
        ValueType * temp_val = new ValueType[used_elements()];
 
        Index * perm_pos;
        perm_pos = perm_row.get_perm_pos();
 
        //permute rows from source to temp_*
        Index new_start(0);
        temp_row_ptr[0] = 0;
        for (Index row(0) ; row < this->rows() ; ++row)
        {
          Index row_size(this->row_ptr()[perm_pos[row] + 1] - this->row_ptr()[perm_pos[row]]);
 
          //iterate over all elements in single one new and old row
          for (Index i(new_start), j(this->row_ptr()[perm_pos[row]]) ; i < new_start + row_size ; ++i, ++j)
          {
            temp_col_ind[i] = this->col_ind()[j];
            temp_val[i] = this->val()[j];
          }
 
          new_start += row_size;
          temp_row_ptr[row+1] = (IT_)new_start;
        }
 
        //use inverse col permutation as lookup table: i -> new location of i
        Adjacency::Permutation perm_col_inv = perm_col.inverse();
        perm_pos = perm_col_inv.get_perm_pos();
 
        //permute columns from temp_* to source
        ::memcpy(this->row_ptr(), temp_row_ptr, (rows() + 1) * sizeof(IT_));
        ::memcpy(this->val(), temp_val, used_elements() * sizeof(ValueType));
        for (Index i(0) ; i < used_elements() ; ++i)
        {
          this->col_ind()[i] = (IT_)perm_pos[temp_col_ind[i]];
        }
 
        delete[] temp_row_ptr;
        delete[] temp_col_ind;
        delete[] temp_val;
 
        //sort columns in every row by column index
        IT_ swap_key;
        ValueType swap_val;
        for (Index row(0) ; row < rows() ; ++row)
        {
          Index offset(this->row_ptr()[row]);
          Index row_size(this->row_ptr()[row+1] - this->row_ptr()[row]);
          for (Index i(1), j ; i < row_size ; ++i)
          {
            swap_key = this->col_ind()[i + offset];
            swap_val = this->val()[i + offset];
            j = i;
            while (j > 0 && this->col_ind()[j - 1 + offset] > swap_key)
            {
              this->col_ind()[j + offset] = this->col_ind()[j - 1 + offset];
              this->val()[j + offset] = this->val()[j - 1 + offset];
              --j;
            }
            this->col_ind()[j + offset] = swap_key;
            this->val()[j + offset] = swap_val;
          }
        }
      }
 
      template<int BHD_>
      void add_trace_double_mat_mult(
        typename LAFEM::SparseMatrixBCSR< DT_, IT_, BHD_, BlockHeight>::VectorTypeL& v,
        const LAFEM::SparseMatrixBCSR<DT_, IT_, BHD_, BlockHeight>& d,
        const LAFEM::DenseVectorBlocked<DT_, IT_, BlockHeight>& a,
        const DT_ alpha = DT_(1)) const
      {
        typedef LAFEM::SparseMatrixBCSR<DT_, IT_, BHD_, BlockHeight> MatrixD;
 
        static constexpr IT_ bhd = IT_(BHD_);
        static constexpr IT_ bwd = IT_(BlockHeight);
 
        static constexpr IT_ bh = IT_(BlockHeight);
        static constexpr IT_ bw = IT_(BlockWidth);
 
        const auto& b(*this);
 
        // Check matrix dimensions
        XASSERT(v.size() == d.rows());
        XASSERT(d.columns() == a.size());
        XASSERT(a.size() == b.rows());
        XASSERT(b.columns() == this->columns());
 
        // Fetch matrix arrays
        // Use POD here to avoid hassle with 1x1 matrices and stuff
        DT_* data_v = v.template elements<LAFEM::Perspective::pod>();
        const DT_* data_d = d.template val<LAFEM::Perspective::pod>();
        const DT_* data_a = a.template elements<LAFEM::Perspective::pod>();
        const DT_* data_b = b.template val<LAFEM::Perspective::pod>();
        const IT_* row_ptr_d = d.row_ptr();
        const IT_* col_idx_d = d.col_ind();
        const IT_* row_ptr_b = b.row_ptr();
        const IT_* col_idx_b = b.col_ind();
 
        // Loop over all rows i of v
        for(IT_ row(0); row < IT_(v.size()); ++row)
        {
          // For all non-zeros D_ik in row i of D
          for(IT_ pos_d(row_ptr_d[row]); pos_d  < row_ptr_d[row+1]; ++pos_d)
          {
            // Get column index k
            const IT_ col_d(col_idx_d[pos_d]);
 
            // Pre-compute block factor (alpha * D_ik * A_kk)
            typename MatrixD::ValueType omega(0);
            for(IT_ i(0); i < bhd; ++i)
            {
              for(IT_ j(0); j < bwd; ++j)
              {
                omega[int(i)][int(j)] = alpha * data_d[bhd*bwd*pos_d + bwd*i + j] * data_a[bwd*col_d + j];
              }
            }
 
            // Loop over all non-zero blocks B_kj in row j of B and perform a "sparse axpy" of B_l onto v_i, i.e.:
            // v_i += alpha * (D_ik * A_kk) * B_k.
            IT_ pos_b(row_ptr_b[col_d]);
            while(pos_b < row_ptr_b[col_d+1])
            {
              if(col_idx_b[pos_b] == row)
              {
                // B_kj contributes to v_i here
                for(IT_ i(0); i < bhd; ++i)
                {
                  for(IT_ j(0); j < bh; ++j)
                  {
                    data_v[bhd*row + i] += omega[int(i)][int(j)]*data_b[bh*bw*pos_b + bw*j + i];
                  }
                }
                break;
              }
              ++pos_b;
            }
          }
        }
      }
 
      friend std::ostream & operator<< (std::ostream & lhs, const SparseMatrixBCSR & b)
      {
        lhs << "[\n";
        for (Index i(0) ; i < b.rows() ; ++i)
        {
          for (int k(0) ; k < BlockHeight_ ; ++k)
          {
            lhs << "[";
            for (Index j(0) ; j < b.columns() ; ++j)
            {
              for (int l(0) ; l < BlockWidth_ ; ++l)
                lhs << "  " << stringify(b(i, j).v[k][l]);
            }
            lhs << "]\n";
          }
        }
        lhs << "]\n";
 
        return lhs;
      }
 
      // Returns a new compatible L-Vector.
      VectorTypeL create_vector_l() const
      {
        return VectorTypeL(this->rows());
      }
 
      // Returns a new compatible R-Vector.
      VectorTypeR create_vector_r() const
      {
        return VectorTypeR(this->columns());
      }
 
      Index get_length_of_line(const Index row) const
      {
        const auto * prow_ptr(this->row_ptr());
        const auto trow = row / Index(BlockHeight_);
        return Index(prow_ptr[trow + 1] - prow_ptr[trow]) * Index(BlockWidth_);
      }
 
      void set_line(const Index row, DT_ * const pval_set, IT_ * const pcol_set,
                    const Index col_start, const Index stride = 1) const
      {
        const auto * prow_ptr(this->row_ptr());
        const auto * pcol_ind(this->col_ind());
        const auto * pval(this->val());
 
        const auto trow = int(row) / BlockHeight_;
        const auto lrow = int(row) - trow * BlockHeight_;
 
        const Index start((Index(prow_ptr[trow])));
        const Index end((Index(prow_ptr[trow + 1] - prow_ptr[trow])));
        for (Index i(0); i < end; ++i)
        {
          for (Index ti(0); ti < Index(BlockWidth_); ++ti)
          {
            pval_set[(i * Index(BlockWidth_) + ti) * stride] = pval[start + i](lrow, int(ti));
            pcol_set[(i * Index(BlockWidth_) + ti) * stride] = pcol_ind[start + i] * IT_(BlockWidth_) + IT_(ti) + IT_(col_start);
          }
        }
      }
 
      void set_line_reverse(const Index row, const DT_ * const pval_set, const Index stride = 1)
      {
        const auto * prow_ptr(this->row_ptr());
        auto * pval(this->val());
 
        const auto trow = int(row) / BlockHeight_;
        const auto lrow = int(row) - trow * BlockHeight_;
 
        const Index start((Index(prow_ptr[trow])));
        const Index end((Index(prow_ptr[trow + 1] - prow_ptr[trow])));
        for (Index i(0); i < end; ++i)
        {
          for (Index ti(0); ti < Index(BlockWidth_); ++ti)
          {
            pval[start + i](lrow, int(ti)) = pval_set[(i * Index(BlockWidth_) + ti) * stride];
          }
        }
      }
 
      /* ******************************************************************* */
      /*  A D J A C T O R   I N T E R F A C E   I M P L E M E N T A T I O N  */
      /* ******************************************************************* */
    public:
      inline Index get_num_nodes_domain() const
      {
        return rows();
      }
 
      inline Index get_num_nodes_image() const
      {
        return columns();
      }
 
      inline ImageIterator image_begin(Index domain_node) const
      {
        XASSERTM(domain_node < rows(), "Domain node index out of range");
        return &this->_indices.at(0)[this->_indices.at(1)[domain_node]];
      }
 
      inline ImageIterator image_end(Index domain_node) const
      {
        XASSERTM(domain_node < rows(), "Domain node index out of range");
        return &this->_indices.at(0)[this->_indices.at(1)[domain_node + 1]];
      }
    }; // class SparseMatrixBCSR
  } // namespace LAFEM
} // namespace FEAT