sparse_matrix_cscr.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_row_col.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/lafem/arch/diagonal.hpp>
#include <kernel/lafem/arch/lumping.hpp>
#include <kernel/lafem/arch/row_norm.hpp>
#include <kernel/adjacency/graph.hpp>
#include <kernel/adjacency/permutation.hpp>
#include <kernel/util/statistics.hpp>
#include <kernel/util/time_stamp.hpp>
#include <kernel/lafem/vector_mirror.hpp>
 
#include <fstream>
 
 
namespace FEAT
{
  namespace LAFEM
  {
    template <typename DT_, typename IT_ = Index>
    class SparseMatrixCSCR : public Container<DT_, IT_>
    {
    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 & _used_rows()
      {
        return this->_scalar_index.at(4);
      }
 
    public:
      typedef DT_ DataType;
      typedef IT_ IndexType;
      typedef DT_ ValueType;
      typedef DenseVector<DT_, IT_> VectorTypeL;
      typedef DenseVector<DT_, IT_> VectorTypeR;
      static constexpr SparseLayoutId layout_id = SparseLayoutId::lt_cscr;
      typedef const IT_* ImageIterator;
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      using ContainerType = SparseMatrixCSCR<DT2_, IT2_>;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
      static constexpr bool is_global = false;
      static constexpr bool is_local = true;
 
      explicit SparseMatrixCSCR() :
        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 SparseMatrixCSCR(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);
        this->_scalar_index.push_back(0);
      }
 
      explicit SparseMatrixCSCR(Index rows_in, Index columns_in, Index used_elements_in, Index used_rows_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->_scalar_index.push_back(used_rows_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_>(_used_rows() + 1));
        this->_indices_size.push_back(_used_rows() + 1);
 
        this->_indices.push_back(MemoryPool::template allocate_memory<IT_>(_used_rows()));
        this->_indices_size.push_back(_used_rows());
 
        this->_elements.push_back(MemoryPool::template allocate_memory<DT_>(_used_elements()));
        this->_elements_size.push_back(_used_elements());
      }
 
      explicit SparseMatrixCSCR(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()));
        this->_elements_size.push_back(_used_elements());
      }
 
      template <typename MT_>
      explicit SparseMatrixCSCR(const MT_ & other) :
        Container<DT_, IT_>(other.size())
      {
        convert(other);
      }
 
      explicit SparseMatrixCSCR(FileMode mode, const String& filename) :
        Container<DT_, IT_>(0)
      {
        read_from(mode, filename);
      }
 
      explicit SparseMatrixCSCR(FileMode mode, std::istream& file) :
        Container<DT_, IT_>(0)
      {
        read_from(mode, file);
      }
 
      explicit SparseMatrixCSCR(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,
                               DenseVector<IT_, IT_> & row_numbers_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(row_numbers_in.size() > 0);
        XASSERT(row_numbers_in.size() + 1 == row_ptr_in.size());
 
        this->_scalar_index.push_back(rows_in);
        this->_scalar_index.push_back(columns_in);
        this->_scalar_index.push_back(val_in.size());
        this->_scalar_index.push_back(row_numbers_in.size());
 
        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());
        this->_indices.push_back(row_numbers_in.elements());
        this->_indices_size.push_back(row_numbers_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));
      }
 
      explicit SparseMatrixCSCR(const SparseMatrixCSR<DT_, IT_> & csr, const VectorMirror<DT_, IT_> & non_zero_rows) :
        Container<DT_, IT_>(csr.size())
      {
        XASSERT(non_zero_rows.num_indices() > 0);
 
        const IT_ * row_main(nullptr);
        const IT_ * indices(nullptr);
        IT_ * trow_main(nullptr);
        IT_ * tindices(nullptr);
        row_main = csr.row_ptr();
        indices = non_zero_rows.indices();
 
        Index tused_elements(0);
        for (Index i(0) ; i < non_zero_rows.num_indices() ; ++i)
        {
          const Index row(indices[i]);
          tused_elements += row_main[row+1] - row_main[row];
        }
 
        this->_scalar_index.push_back(csr.rows());
        this->_scalar_index.push_back(csr.columns());
        this->_scalar_index.push_back(tused_elements);
        this->_scalar_index.push_back(non_zero_rows.num_indices());
 
        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_>(_used_rows() + 1));
        this->_indices_size.push_back(_used_rows() + 1);
 
        this->_indices.push_back(MemoryPool::template allocate_memory<IT_>(_used_rows()));
        this->_indices_size.push_back(_used_rows());
 
        this->_elements.push_back(MemoryPool::template allocate_memory<DT_>(_used_elements()));
        this->_elements_size.push_back(_used_elements());
 
        Index row_count(0);
        IT_ offset(0);
        this->row_ptr()[0]=offset;
        for (Index i(0) ; i < non_zero_rows.num_indices() ; ++i)
        {
          const IT_ row(indices[i]);
          const IT_ start(row_main[row]);
          const IT_ end(row_main[row+1]);
          const IT_ row_size(end - start);
          MemoryPool::set_memory(this->row_ptr() + row_count + 1, offset + row_size);
          MemoryPool::set_memory(this->row_numbers() + row_count, row);
          MemoryPool::copy(this->val() + offset, csr.val() + start, row_size);
          MemoryPool::copy(this->col_ind() + offset, csr.col_ind() + start, row_size);
          ++row_count;
          offset += row_size;
        }
 
        delete[] trow_main;
        delete[] tindices;
      }
 
      explicit SparseMatrixCSCR(const Adjacency::Graph & graph) :
        Container<DT_, IT_>(0)
      {
        // get number of rows, columns and indices
        Index num_rows = graph.get_num_nodes_domain();
        Index num_cols = graph.get_num_nodes_image();
        Index num_nzes = graph.get_num_indices();
 
        if (num_nzes == 0)
        {
          this->move(SparseMatrixCSCR(num_rows, num_cols));
          return;
        }
 
        // get graph arrays
        const Index* dom_ptr = graph.get_domain_ptr();
        const Index* img_idx = graph.get_image_idx();
 
        // count number of non-empty rows
        Index num_nzrs = Index(0);
        for(Index i(0); i < num_rows; ++i)
        {
          num_nzrs += Index(dom_ptr[i] < dom_ptr[i+1] ? 1 : 0);
        }
 
        // allocate output matrix
        SparseMatrixCSCR<DT_, IT_> matrix(num_rows, num_cols, num_nzes, num_nzrs);
 
        // get matrix arrays
        IndexType* trow_ptr = matrix.row_ptr();
        IndexType* trow_idx = matrix.row_numbers();
        IndexType* tcol_idx = matrix.col_ind();
 
        // fill arrays
        trow_ptr[0] = IndexType(dom_ptr[0]);
        for(Index i(0), j(0); i < num_rows; ++i)
        {
          if(dom_ptr[i] < dom_ptr[i + 1])
          {
            ASSERT(j < num_nzrs);
            trow_idx[  j] = IndexType(i);
            trow_ptr[++j] = IndexType(dom_ptr[i+1]);
          }
        }
 
        for(Index k(0); k < num_nzes; ++k)
          tcol_idx[k] = IndexType(img_idx[k]);
 
        this->convert(matrix);
      }
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      explicit SparseMatrixCSCR(std::vector<char> input) :
        Container<DT_, IT_>(0)
      {
        deserialize<DT2_, IT2_>(input);
      }
 
      SparseMatrixCSCR(SparseMatrixCSCR && other) :
        Container<DT_, IT_>(std::forward<SparseMatrixCSCR>(other))
      {
      }
 
      SparseMatrixCSCR & operator= (SparseMatrixCSCR && other)
      {
        this->move(std::forward<SparseMatrixCSCR>(other));
 
        return *this;
      }
 
      SparseMatrixCSCR clone(CloneMode clone_mode = CloneMode::Weak) const
      {
        SparseMatrixCSCR t;
        t.clone(*this, clone_mode);
        return t;
      }
 
      template<typename DT2_, typename IT2_>
      void clone(const SparseMatrixCSCR<DT2_, IT2_> & other, CloneMode clone_mode = CloneMode::Weak)
      {
        Container<DT_, IT_>::clone(other, clone_mode);
      }
 
      template <typename DT2_, typename IT2_>
      void convert(const SparseMatrixCSCR<DT2_, IT2_> & other)
      {
        this->assign(other);
      }
 
      template <typename MT_>
      void convert(const MT_ & a)
      {
        typename MT_::template ContainerType<DT_, IT_> ta;
        ta.convert(a);
 
        const Index arows(ta.template rows<Perspective::pod>());
        const Index acolumns(ta.template columns<Perspective::pod>());
        const Index aused_elements(ta.template used_elements<Perspective::pod>());
        Index aused_rows(0);
 
        for (Index i(0); i < arows; ++i)
        {
          aused_rows += ta.get_length_of_line(i) > 0;
        }
 
        DenseVector<DT_, IT_> tval(aused_elements);
        DenseVector<IT_, IT_> tcol_ind(aused_elements);
        DenseVector<IT_, IT_> trow_ptr(aused_rows + 1);
        DenseVector<IT_, IT_> trow_numbers(aused_rows);
 
        DT_ * pval(tval.elements());
        IT_ * pcol_ind(tcol_ind.elements());
        IT_ * prow_ptr(trow_ptr.elements());
        IT_ * prow_numbers(trow_numbers.elements());
 
        Index row_index(0);
        Index offset(0);
        prow_ptr[0] = IT_(0);
        for (Index i(0); i < arows; ++i)
        {
          Index length = ta.get_length_of_line(i);
          if (length > 0)
          {
            prow_ptr[row_index + 1] = IT_(offset + length);
            prow_numbers[row_index] = IT_(i);
            ++row_index;
            offset += length;
          }
        }
 
 
        for (Index i(0); i < arows; ++i)
        {
          ta.set_line(i, pval + prow_ptr[i], pcol_ind + prow_ptr[i], 0);
        }
 
        this->move(SparseMatrixCSCR<DT_, IT_>(arows, acolumns, tcol_ind, tval, trow_ptr, trow_numbers));
      }
 
      SparseMatrixCSCR & 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()));
        this->_elements_size.push_back(_used_elements());
 
        return *this;
      }
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      void deserialize(std::vector<char> input)
      {
        this->template _deserialize<DT2_, IT2_>(FileMode::fm_cscr, 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_cscr, 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_cscr:
        case FileMode::fm_binary:
          this->template _deserialize<double, std::uint64_t>(FileMode::fm_cscr, file);
          break;
        default:
          XABORTM("Filemode not supported!");
        }
      }
 
      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_cscr:
        case FileMode::fm_binary:
          this->template _serialize<double, std::uint64_t>(FileMode::fm_cscr, file);
          break;
        default:
          XABORTM("Filemode not supported!");
        }
      }
 
      DT_ operator()(Index row, Index col) const
      {
        ASSERT(row < rows());
        ASSERT(col < columns());
 
        Index row_index(0);
        while (row_index < this->used_rows() && Index(this->row_numbers()[row_index]) < row)
        {
          ++row_index;
        }
 
        if (row_index == this->used_rows() || Index(this->row_numbers()[row_index]) > row)
          return DT_(0.);
 
        //row_numbers[row_index] == row
        for (Index i(Index(this->row_ptr()[row_index])) ; i < Index(this->row_ptr()[row_index + 1]) ; ++i)
        {
          if (Index(this->_indices.at(0)[i]) == col)
            return this->_elements.at(0)[i];
          if (Index(this->_indices.at(0)[i]) > col)
            return DT_(0.);
        }
        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_rows() const
      {
        return this->_scalar_index.at(4);
      }
 
      template <Perspective = Perspective::native>
      Index used_elements() const
      {
        return this->_scalar_index.at(3);
      }
 
      IT_ * col_ind()
      {
        if (this->_indices.size() == 0)
          return nullptr;
 
        return this->_indices.at(0);
      }
 
      IT_ const * col_ind() const
      {
        if (this->_indices.size() == 0)
          return nullptr;
 
        return this->_indices.at(0);
      }
 
      DT_ * val()
      {
        if (this->_elements.size() == 0)
          return nullptr;
 
        return this->_elements.at(0);
      }
 
      DT_ const * val() const
      {
        if (this->_elements.size() == 0)
          return nullptr;
 
        return 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);
      }
 
      IT_ * row_numbers()
      {
        if (this->_indices.size() == 0)
          return nullptr;
 
        return this->_indices.at(2);
      }
 
      IT_ const * row_numbers() const
      {
        if (this->_indices.size() == 0)
          return nullptr;
 
        return this->_indices.at(2);
      }
 
      static String name()
      {
        return "SparseMatrixCSCR";
      }
 
      void copy(const SparseMatrixCSCR & x, bool full = false)
      {
        this->_copy_content(x, full);
      }
 
 
      void axpy(
                const SparseMatrixCSCR & 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() * 2);
        Arch::Axpy::value(this->val(), alpha, x.val(), this->used_elements());
 
        TimeStamp ts_stop;
        Statistics::add_time_axpy(ts_stop.elapsed(ts_start));
      }
 
      void scale(const SparseMatrixCSCR & 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());
        Arch::Scale::value(this->val(), x.val(), alpha, this->used_elements());
 
        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->used_elements());
 
        TimeStamp ts_stop;
        Statistics::add_time_reduction(ts_stop.elapsed(ts_start));
 
        return result;
      }
 
      DT_ max_rel_diff(const SparseMatrixCSCR& 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 SparseMatrixCSCR& 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->used_elements() != x.used_elements())
          return false;
        if (this->used_rows() != x.used_rows())
          return false;
 
        IT_ * col_ind_a;
        IT_ * col_ind_b;
        IT_ * row_ptr_a;
        IT_ * row_ptr_b;
 
        bool ret(true);
 
        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());
        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->used_rows() + 1; ++i)
          {
            if (row_ptr_a[i] != row_ptr_b[i])
            {
              ret = false;
              break;
            }
          }
        }
 
        return ret;
      }
 
      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!");
 
        TimeStamp ts_start;
 
        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() * 2);
        Arch::Apply::cscr(r.elements(), DT_(1), x.elements(), DT_(0), r.elements(),
            this->val(), this->col_ind(), this->row_ptr(), this->row_numbers(), this->used_rows(), this->rows(), this->columns(), this->used_elements(), false);
 
        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!");
 
        TimeStamp ts_start;
 
        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() * 2);
        Arch::Apply::cscr(r.elements(), DT_(1), x.elements(), DT_(0), r.elements(),
          this->val(), this->col_ind(), this->row_ptr(), this->row_numbers(), this->used_rows(), this->rows(), this->columns(), this->used_elements(), true);
 
        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!");
 
        TimeStamp ts_start;
 
        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() + this->rows()) * 2 );
        Arch::Apply::cscr(r.elements(), alpha, x.elements(), DT_(1.), y.elements(),
            this->val(), this->col_ind(), this->row_ptr(), this->row_numbers(), this->used_rows(), this->rows(), this->columns(), this->used_elements(), false);
 
        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!");
 
        TimeStamp ts_start;
 
        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() + this->rows()) * 2 );
        Arch::Apply::cscr(r.elements(), alpha, x.elements(), DT_(1.), y.elements(),
          this->val(), this->col_ind(), this->row_ptr(), this->row_numbers(), this->used_rows(), this->rows(), this->columns(), this->used_elements(), true);
 
        TimeStamp ts_stop;
        Statistics::add_time_blas2(ts_stop.elapsed(ts_start));
      }
 
      VectorTypeL create_vector_l() const
      {
        return VectorTypeL(this->rows());
      }
 
      VectorTypeR create_vector_r() const
      {
        return VectorTypeR(this->columns());
      }
 
      Index get_length_of_line(const Index row) const
      {
        const IT_ * prow_ptr(this->row_ptr());
        return Index(prow_ptr[row + 1] - prow_ptr[row]);
      }
 
 
      void set_line(const Index row, DT_ * const pval_set, IT_ * const pcol_set,
                    const Index col_start, const Index stride = 1) const
      {
        const IT_ * prow_ptr(this->row_ptr());
        const IT_ * pcol_ind(this->col_ind());
        const DT_ * pval(this->val());
 
        Index row_idx(0);
        while (row_idx < this->used_rows() && row_idx < row)
        {
          ++row_idx;
        }
        XASSERTM(row_idx < this->used_rows(), "invalid row number provided!");
 
        const Index start((Index(prow_ptr[row_idx])));
        const Index end((Index(prow_ptr[row_idx + 1] - prow_ptr[row_idx])));
        for (Index i(0); i < end; ++i)
        {
          pval_set[i * stride] = pval[start + i];
          pcol_set[i * stride] = pcol_ind[start + i] + IT_(col_start);
        }
      }
 
 
      friend std::ostream & operator<< (std::ostream & lhs, const SparseMatrixCSCR & 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;
      }
    }; //SparseMatrixCSCR
 
#ifdef FEAT_EICKT
    extern template class SparseMatrixCSCR<float, std::uint32_t>;
    extern template class SparseMatrixCSCR<double, std::uint32_t>;
    extern template class SparseMatrixCSCR<float, std::uint64_t>;
    extern template class SparseMatrixCSCR<double, std::uint64_t>;
#endif
 
  } // namespace LAFEM
} // namespace FEAT