sparse_matrix_banded.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/util/math.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/norm.hpp>
#include <kernel/adjacency/graph.hpp>
 
#include <fstream>
#include <set>
 
 
namespace FEAT
{
  namespace LAFEM
  {
    template <typename DT_, typename IT_ = Index>
    class SparseMatrixBanded : public Container<DT_, IT_>
    {
    public:
      class ScatterAxpy
      {
      public:
        typedef LAFEM::SparseMatrixBanded<DT_, IT_> MatrixType;
        typedef DT_ DataType;
        typedef IT_ IndexType;
 
      private:
#ifdef DEBUG
        const IT_ _deadcode;
#endif
        Index _num_rows;
        Index _num_cols;
        Index _num_of_offsets;
        IT_* _offsets;
        IT_* _col_ptr;
        DT_* _data;
 
      public:
        explicit ScatterAxpy(MatrixType& matrix) :
#ifdef DEBUG
          _deadcode(~IT_(0)),
#endif
          _num_rows(matrix.rows()),
          _num_cols(matrix.columns()),
          _num_of_offsets(matrix.num_of_offsets()),
          _offsets(matrix.offsets()),
          _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(0); k < _num_of_offsets; ++k)
            {
              if(_offsets[k] + ix + 1 >= _num_rows && _offsets[k] + ix + 1 < 2 * _num_rows)
              {
                _col_ptr[_offsets[k] + ix + 1 - _num_rows] = IT_(k * _num_rows + ix);
              }
            }
 
            // 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(0); k < _num_of_offsets; ++k)
            {
              if(_offsets[k] + ix + 1 >= _num_rows && _offsets[k] + ix + 1 < 2 * _num_rows)
              {
                _col_ptr[_offsets[k] + ix + 1 - _num_rows] = _deadcode;
              }
            }
#endif
            // continue with next row entry
          }
        }
      }; // class ScatterAxpy
 
      class GatherAxpy
      {
      public:
        typedef LAFEM::SparseMatrixBanded<DT_, IT_> MatrixType;
        typedef DT_ DataType;
        typedef IT_ IndexType;
 
      private:
#ifdef DEBUG
        const IT_ _deadcode;
#endif
        Index _num_rows;
        Index _num_cols;
        Index _num_of_offsets;
        const IT_* _offsets;
        IT_* _col_ptr;
        const DT_* _data;
 
      public:
        explicit GatherAxpy(const MatrixType& matrix) :
#ifdef DEBUG
          _deadcode(~IT_(0)),
#endif
          _num_rows(matrix.rows()),
          _num_cols(matrix.columns()),
          _num_of_offsets(matrix.num_of_offsets()),
          _offsets(matrix.offsets()),
          _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 ~GatherAxpy()
        {
          if(_col_ptr != nullptr)
          {
            delete [] _col_ptr;
          }
        }
 
        template<typename LocalMatrix_, typename RowMapping_, typename ColMapping_>
        void operator()(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(0); k < _num_of_offsets; ++k)
            {
              if(_offsets[k] + ix + 1 >= _num_rows && _offsets[k] + ix + 1 < 2 * _num_rows)
              {
                _col_ptr[_offsets[k] + ix + 1 - _num_rows] = IT_(k * _num_rows + ix);
              }
            }
 
            // 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");
 
              // update local matrix data
              loc_mat[i][j] += alpha * _data[_col_ptr[jx]];
 
              // continue with next column entry
            }
 
#ifdef DEBUG
            // reformat column-pointer array
            for(IT_ k(0); k < _num_of_offsets; ++k)
            {
              if(_offsets[k] + ix + 1 >= _num_rows && _offsets[k] + ix + 1 < 2 * _num_rows)
              {
                _col_ptr[_offsets[k] + ix + 1 - _num_rows] = _deadcode;
              }
            }
#endif
            // continue with next row entry
          }
        }
      }; // class GatherAxpy
 
    public: //shall be 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);
      }
 
      Index & _num_of_offsets()
      {
        return this->_scalar_index.at(4);
      }
 
    public:
      typedef DT_ DataType;
      typedef IT_ IndexType;
      typedef DenseVector<DataType, IT_> VectorTypeL;
      typedef DenseVector<DataType, IT_> VectorTypeR;
      static constexpr SparseLayoutId layout_id = SparseLayoutId::lt_banded;
      typedef DT_ ValueType;
      class ImageIterator
      {
      private:
        IT_ _k;
        const IT_ _s;
        const IT_ * const _offsets;
 
      public:
        ImageIterator() : _k(IT_(0)), _s(IT_(0)), _offsets(nullptr) {}
 
        ImageIterator(IT_ k, IT_ s, const IT_ * offsets) : _k(k), _s(s), _offsets(offsets) {}
 
        ImageIterator& operator=(const ImageIterator& other)
        {
          _k       = other._k;
          _s       = other._s;
          _offsets = other._offsets;
          return *this;
        }
 
        bool operator!=(const ImageIterator& other) const
        {
          return this->_k != other._k;
        }
 
        ImageIterator& operator++()
        {
          ++_k;
          return *this;
        }
 
        Index operator*() const
        {
          return Index(_s + _offsets[_k]);
        }
      };
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      using ContainerType = SparseMatrixBanded<DT2_, IT2_>;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
      explicit SparseMatrixBanded() :
        Container<DT_, IT_> (0)
      {
        this->_scalar_index.push_back(0);
        this->_scalar_index.push_back(0);
        this->_scalar_index.push_back(0);
        this->_scalar_index.push_back(0);
      }
 
      explicit SparseMatrixBanded(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_>(rows() * num_of_offsets()));
        this->_elements_size.push_back(rows() * num_of_offsets());
      }
 
      explicit SparseMatrixBanded(const Index rows_in, const Index columns_in,
                                  DenseVector<DT_, IT_> & val_in,
                                  DenseVector<IT_, IT_> & offsets_in) :
        Container<DT_, IT_>(rows_in * columns_in)
      {
        if (val_in.size() != rows_in * offsets_in.size())
        {
          XABORTM("Size of values does not match to number of offsets and row count!");
        }
 
        this->_scalar_index.push_back(rows_in);
        this->_scalar_index.push_back(columns_in);
 
        Index tused_elements(0);
 
        for (Index i(0); i < offsets_in.size(); ++i)
        {
          const Index toffset(Index(offsets_in(i)));
 
          if (toffset + Index(2) > rows_in + columns_in)
          {
            XABORTM("Offset out of matrix!");
          }
 
          tused_elements += columns_in + Math::min(rows_in, columns_in + rows_in - toffset - 1) - Math::max(columns_in + rows_in - toffset - 1, columns_in);
        }
 
        this->_scalar_index.push_back(tused_elements);
        this->_scalar_index.push_back(offsets_in.size());
 
        this->_elements.push_back(val_in.elements());
        this->_elements_size.push_back(val_in.size());
        this->_indices.push_back(offsets_in.elements());
        this->_indices_size.push_back(offsets_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 MT_>
      explicit SparseMatrixBanded(const MT_ & other) :
        Container<DT_, IT_>(other.size())
      {
        convert(other);
      }
 
      explicit SparseMatrixBanded(const Adjacency::Graph & graph) :
        Container<DT_, IT_>(0)
      {
        Index num_rows = graph.get_num_nodes_domain();
        Index num_cols = graph.get_num_nodes_image();
 
        const Index * dom_ptr(graph.get_domain_ptr());
        const Index * img_idx(graph.get_image_idx());
 
        std::set<IT_> moffsets;
 
        for(Index row(0); row < num_rows; ++row)
        {
          for(Index j(dom_ptr[row]); j < dom_ptr[row+1]; ++j)
          {
            moffsets.insert(IT_(num_rows - 1 + img_idx[j] - row));
          }
        }
 
        DenseVector<IT_, IT_> toffsets(Index(moffsets.size()));
        auto * ptoffsets = toffsets.elements();
 
        IT_ idx(0);
        for (auto off : moffsets)
        {
          ptoffsets[idx] = off;
          ++idx;
        }
 
        DenseVector<IT_, IT_> toffsets_mem;
        toffsets_mem.convert(toffsets);
        DenseVector<DT_, IT_> tval_mem(Index(moffsets.size()) * num_rows, DT_(0));
 
        this->move(SparseMatrixBanded(num_rows, num_cols, tval_mem, toffsets_mem));
      }
 
      explicit SparseMatrixBanded(FileMode mode, String filename) :
        Container<DT_, IT_>(0)
      {
        read_from(mode, filename);
      }
 
      explicit SparseMatrixBanded(FileMode mode, std::istream& file) :
        Container<DT_, IT_>(0)
      {
        read_from(mode, file);
      }
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      explicit SparseMatrixBanded(std::vector<char> input) :
        Container<DT_, IT_>(0)
      {
        deserialize<DT2_, IT2_>(input);
      }
 
      SparseMatrixBanded(SparseMatrixBanded && other) :
        Container<DT_, IT_>(std::forward<SparseMatrixBanded>(other))
      {
      }
 
      SparseMatrixBanded & operator= (SparseMatrixBanded && other)
      {
        this->move(std::forward<SparseMatrixBanded>(other));
 
        return *this;
      }
 
      SparseMatrixBanded clone(CloneMode clone_mode = CloneMode::Weak) const
      {
        SparseMatrixBanded t;
        t.clone(*this, clone_mode);
        return t;
      }
 
      template<typename DT2_, typename IT2_>
      void clone(const SparseMatrixBanded<DT2_, IT2_> & other, CloneMode clone_mode = CloneMode::Weak)
      {
        Container<DT_, IT_>::clone(other, clone_mode);
      }
 
      SparseMatrixBanded & 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_>(rows() * num_of_offsets()));
        this->_elements_size.push_back(rows() * num_of_offsets());
 
        return *this;
      }
 
      template <typename DT2_, typename IT2_>
      void convert(const SparseMatrixBanded<DT2_, IT2_> & other)
      {
        this->assign(other);
      }
 
      template <typename DT2_, typename IT2_>
      void convert(const SparseMatrixCSR<DT2_, IT2_> & csr)
      {
        XASSERT(csr.template used_elements<Perspective::pod>() > 0);
 
        std::set<IT_> offset_set;
        Index nrows(csr.rows());
        Index ncolumns(csr.columns());
        for (Index row(0) ; row < nrows ; ++row)
        {
          for (Index i(csr.row_ptr()[row]) ; i < csr.row_ptr()[row+1] ; ++i)
          {
            //compute band offset for each non zero entry
            IT_ off = IT_(csr.col_ind()[i] - row + nrows - 1);
            offset_set.insert(off);
          }
        }
 
        DenseVector<DT_, IT_> val_new(Index(offset_set.size()) * nrows, DT_(0));
        for (Index row(0) ; row < nrows ; ++row)
        {
          for (Index i(csr.row_ptr()[row]) ; i < csr.row_ptr()[row+1] ; ++i)
          {
            Index col = csr.col_ind()[i];
            Index offset = (col - row + nrows - 1);
            Index band = Index(std::distance(offset_set.begin(), offset_set.find(IT_(offset))));
            val_new.elements()[band * nrows + row] = csr.val()[i];
          }
        }
 
        DenseVector<IT_, IT_> offsets_new(Index(offset_set.size()));
        auto it_offset = offset_set.begin();
        for (Index i(0) ; i < offset_set.size() ; ++i, ++it_offset)
        {
          offsets_new(i, *it_offset);
        }
        offset_set.clear();
 
        SparseMatrixBanded<DT_, IT_> temp(nrows, ncolumns, val_new, offsets_new);
        this->move(std::move(temp));
      }
 
 
      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_bm:
        case FileMode::fm_binary:
          if (! std::is_same<DT_, double>::value)
            std::cout<<"Warning: You are writing out a banded matrix that is not double precision!\n";
 
        this->template _serialize<double, std::uint64_t>(FileMode::fm_bm, file);
          break;
        default:
          XABORTM("Filemode not supported!");
        }
      }
 
      DT_ operator()(Index row, Index col) const
      {
        ASSERT(row < rows());
        ASSERT(col < columns());
 
        MemoryPool::synchronize();
 
        const Index trows(this->_scalar_index.at(1));
 
        for (Index i(0); i < this->_scalar_index.at(4); ++i)
        {
          const Index toffset(this->offsets()[i]);
          if (row + toffset + Index(1) == col + trows)
          {
            return this->val()[i * trows + row];
          }
        }
        return DT_(0.);
      }
 
      SparseLayout<IT_, layout_id> layout() const
      {
        return SparseLayout<IT_, layout_id>(this->_indices, this->_indices_size, this->_scalar_index);
      }
 
      template <Perspective = Perspective::native>
      Index rows() const
      {
        return this->_scalar_index.at(1);
      }
 
      template <Perspective = Perspective::native>
      Index columns() const
      {
        return this->_scalar_index.at(2);
      }
 
      template <Perspective = Perspective::native>
      Index used_elements() const
      {
        return this->_scalar_index.at(3);
      }
 
      Index num_of_offsets() const
      {
        return this->_scalar_index.at(4);
      }
 
      DT_ * val()
      {
        return this->_elements.at(0);
      }
 
      DT_ const * val() const
      {
        return this->_elements.at(0);
      }
 
      IT_ * offsets()
      {
        return this->_indices.at(0);
      }
 
      IT_ const * offsets() const
      {
        return this->_indices.at(0);
      }
 
      static String name()
      {
        return "SparseMatrixBanded";
      }
 
      void copy(const SparseMatrixBanded & x, bool full = false)
      {
        this->_copy_content(x, full);
      }
 
 
      void axpy(
                const SparseMatrixBanded & 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.num_of_offsets() == this->num_of_offsets(), "Matrix num_of_offsets 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() * 2);
        Arch::Axpy::value(this->val(), alpha, x.val(), this->rows() * this->num_of_offsets());
 
        TimeStamp ts_stop;
        Statistics::add_time_axpy(ts_stop.elapsed(ts_start));
      }
 
      void scale(const SparseMatrixBanded & x, const DT_ alpha)
      {
        XASSERTM(x.rows() == this->rows(), "Matrix rows do not match!");
        XASSERTM(x.columns() == this->columns(), "Matrix columns do not match!");
        XASSERTM(x.num_of_offsets() == this->num_of_offsets(), "Matrix num_of_offsets 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());
        Arch::Scale::value(this->val(), x.val(), alpha, this->rows() * this->num_of_offsets());
 
        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() * 2);
        DT_ result = Arch::Norm2::value(this->val(), this->rows() * this->num_of_offsets());
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      void apply(DenseVector<DT_, IT_>& r, const DenseVector<DT_, IT_>& 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!");
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        TimeStamp ts_start;
        Statistics::add_flops( 2 * this->used_elements() );
 
        Arch::Apply::banded(r.elements(),
            DT_(1),
            x.elements(),
            DT_(0),
            r.elements(),
            this->val(),
            this->offsets(),
            this->num_of_offsets(),
            this->rows(),
            this->columns());
 
        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
      {
        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(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        TimeStamp ts_start;
        Statistics::add_flops( 2 * this->used_elements() );
 
        Arch::Apply::banded_transposed(r.elements(),
          DT_(1),
          x.elements(),
          DT_(0),
          r.elements(),
          this->val(),
          this->offsets(),
          this->num_of_offsets(),
          this->rows(),
          this->columns());
 
        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(), "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!");
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        if (this->used_elements() == 0 || Math::abs(alpha) < Math::eps<DT_>())
        {
          r.copy(y);
          //r.scale(beta);
          return;
        }
 
        TimeStamp ts_start;
        Statistics::add_flops( 2 * (this->used_elements() + this->rows()) );
 
        Arch::Apply::banded(r.elements(),
            alpha,
            x.elements(),
            DT_(1),
            y.elements(),
            this->val(),
            this->offsets(),
            this->num_of_offsets(),
            this->rows(),
            this->columns());
 
        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(), "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!");
 
        XASSERTM(r.template elements<Perspective::pod>() != x.template elements<Perspective::pod>(), "Vector x and r must not share the same memory!");
 
        if (this->used_elements() == 0 || Math::abs(alpha) < Math::eps<DT_>())
        {
          r.copy(y);
          //r.scale(beta);
          return;
        }
 
        TimeStamp ts_start;
        Statistics::add_flops( 2 * (this->used_elements() + this->rows()) );
 
        Arch::Apply::banded_transposed(r.elements(),
          alpha,
          x.elements(),
          DT_(1),
          y.elements(),
          this->val(),
          this->offsets(),
          this->num_of_offsets(),
          this->rows(),
          this->columns());
 
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      void extract_diag(VectorTypeL & diag) const
      {
        XASSERTM(diag.size() == rows(), "diag size does not match matrix row count!");
        XASSERTM(rows() == columns(), "matrix is not square!");
 
        for (Index i(0) ; i < num_of_offsets() ; ++i)
        {
          if (this->offsets()[i] == this->rows() - 1)
          {
            MemoryPool::copy(diag.elements(), val() + i * rows(), rows());
            break;
          }
        }
      }
 
      VectorTypeL extract_diag() const
      {
        VectorTypeL diag = create_vector_l();
        extract_diag(diag);
        return diag;
      }
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      void deserialize(std::vector<char> input)
      {
        this->template _deserialize<DT2_, IT2_>(FileMode::fm_bm, input);
      }
 
      DT_ max_rel_diff(const SparseMatrixBanded& 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->val(), x.val(), this->used_elements());
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return max_rel_diff;
      }
 
      bool same_layout(const SparseMatrixBanded& x) const
      {
        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->num_of_offsets() != x.num_of_offsets())
          return false;
        if (this->used_elements() != x.used_elements())
          return false;
 
        IT_ * offsets_a;
        IT_ * offsets_b;
 
        offsets_a = const_cast<IT_*>(this->offsets());
        offsets_b = const_cast<IT_*>(x.offsets());
 
        bool ret(true);
 
        for (Index i(0); i < this->num_of_offsets(); ++i)
        {
          if (offsets_a[i] != offsets_b[i])
          {
            ret = false;
            break;
          }
        }
 
        return ret;
      }
 
      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_bm, 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_bm:
        case FileMode::fm_binary:
          this->template _deserialize<double, std::uint64_t>(FileMode::fm_bm, file);
          break;
        default:
          XABORTM("Filemode not supported!");
        }
      }
 
      VectorTypeL create_vector_l() const
      {
        return VectorTypeL(this->rows());
      }
 
      VectorTypeR create_vector_r() const
      {
        return VectorTypeR(this->columns());
      }
 
      template <typename ITX_>
      inline Index _start_offset(const Index i, const ITX_ * const offsets,
                                const Index rows_in, const Index columns_in, const Index noo) const
      {
        if (i == Index(-1))
        {
          return rows_in;
        }
        else if (i == noo)
        {
          return Index(0);
        }
        else
        {
          return Math::max(columns_in + Index(1), rows_in + columns_in - Index(offsets[i])) - columns_in - Index(1);
        }
      }
 
      template <typename ITX_>
      inline Index _end_offset(const Index i, const ITX_ * const offsets,
                              const Index rows_in, const Index columns_in, const Index noo) const
      {
        if (i == Index (-1))
        {
          return rows_in - 1;
        }
        else if (i == noo)
        {
          return Index(-1);
        }
        else
        {
          return Math::min(rows_in, columns_in + rows_in - Index(offsets[i]) - Index(1)) - Index(1);
        }
      }
 
      Index start_offset(const Index i) const
      {
        return this->_start_offset(i, offsets(), rows(), columns(), num_of_offsets());
      }
 
      Index end_offset(const Index i) const
      {
        return this->_end_offset(i, offsets(), rows(), columns(), num_of_offsets());
      }
 
      /* ******************************************************************* */
      /*  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
      {
        XASSERT(domain_node < rows());
 
        const IT_ * toffsets(offsets());
        const IT_ tnum_of_offsets(num_of_offsets());
        const Index trows(rows());
 
        for(IT_ k(0); k < tnum_of_offsets; ++k)
        {
          if(toffsets[k] + domain_node + 1 >= trows)
            return ImageIterator(k, domain_node + 1 - trows, toffsets);
        }
      }
 
      inline ImageIterator image_end(Index domain_node) const
      {
        XASSERT(domain_node < rows());
 
        const IT_ * toffsets(offsets());
        const IT_ tnum_of_offsets(num_of_offsets());
        const Index trows(rows());
 
        for(IT_ k(0); k < tnum_of_offsets; ++k)
        {
          if(toffsets[k] + domain_node + 1 >= 2 * trows)
          {
            return ImageIterator(k, domain_node + 1 - trows, toffsets);
          }
        }
        return ImageIterator(tnum_of_offsets, domain_node + 1 - trows, toffsets);
      }
 
      friend std::ostream & operator<< (std::ostream & lhs, const SparseMatrixBanded & b)
      {
        lhs << "[\n";
        for (Index i(0) ; i < b.rows() ; ++i)
        {
          lhs << "[";
          for (Index j(0) ; j < b.columns() ; ++j)
          {
            lhs << "  " << stringify(b(i, j));
          }
          lhs << "]\n";
        }
        lhs << "]\n";
 
        return lhs;
      }
    };
 
  } // namespace LAFEM
} // namespace FEAT