tuple_matrix.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
 
#include <kernel/lafem/tuple_vector.hpp>
 
 
namespace FEAT
{
  namespace LAFEM
  {
    template<typename First_, typename... Rest_>
    class TupleMatrixRow
    {
    private:
      template<typename,typename...>
      friend class TupleMatrixRow;
 
      typedef TupleMatrixRow<Rest_...> RestClass;
 
    public:
      static constexpr int num_blocks = TupleMatrixRow<Rest_...>::num_blocks + 1;
 
      typedef typename First_::DataType DataType;
      typedef typename First_::IndexType IndexType;
 
      // ensure that all sub-matrix have the same data-type
      static_assert(std::is_same<DataType, typename RestClass::DataType>::value,
                    "sub-matrices have different data-types");
      static_assert(std::is_same<IndexType, typename RestClass::IndexType>::value,
                    "sub-matrices have different index-types");
 
      template <typename DT2_ = DataType, typename IT2_ = IndexType>
      using ContainerType = TupleMatrixRow<
        typename First_::template ContainerType<DT2_, IT2_>,
        typename Rest_::template ContainerType<DT2_, IT2_>...>;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
      typedef typename First_::VectorTypeL VectorTypeL;
      typedef TupleVector<typename First_::VectorTypeR, typename Rest_::VectorTypeR...> VectorTypeR;
 
    protected:
      First_ _first;
      RestClass _rest;
 
      static String sub_name_list()
      {
        return First_::name() + "," + RestClass::sub_name_list();
      }
 
    public:
      TupleMatrixRow()
      {
      }
 
      explicit TupleMatrixRow(First_&& the_first, RestClass&& the_rest) :
        _first(std::forward<First_>(the_first)),
        _rest(std::forward<RestClass>(the_rest))
      {
      }
 
      explicit TupleMatrixRow(First_&& the_first, Rest_&&... the_rest) :
        _first(std::forward<First_>(the_first)),
        _rest(std::forward<Rest_>(the_rest)...)
      {
      }
 
      TupleMatrixRow(TupleMatrixRow&& other) :
        _first(std::forward<First_>(other._first)),
        _rest(std::forward<RestClass>(other._rest))
      {
      }
 
      TupleMatrixRow& operator=(TupleMatrixRow&& other)
      {
        if(this != &other)
        {
          _first = std::forward<First_>(other._first);
          _rest = std::forward<RestClass>(other._rest);
        }
        return *this;
      }
 
      TupleMatrixRow(const TupleMatrixRow&) = delete;
      TupleMatrixRow& operator=(const TupleMatrixRow&) = delete;
 
      TupleMatrixRow clone(LAFEM::CloneMode mode = LAFEM::CloneMode::Weak) const
      {
        return TupleMatrixRow(_first.clone(mode), _rest.clone(mode));
      }
 
      void clone(const TupleMatrixRow& other, CloneMode clone_mode)
      {
        _first.clone(other._first, clone_mode);
        _rest.clone(other._rest, clone_mode);
      }
 
      void clone(const TupleMatrixRow& other)
      {
        _first.clone(other._first);
        _rest.clone(other._rest);
      }
 
      First_& first()
      {
        return _first;
      }
 
      const First_& first() const
      {
        return _first;
      }
 
      TupleMatrixRow<Rest_...>& rest()
      {
        return _rest;
      }
 
      const TupleMatrixRow<Rest_...>& rest() const
      {
        return _rest;
      }
 
      template<int i_>
      typename TupleElement<i_, First_, Rest_...>::Type& at()
      {
        static_assert((0 <= i_) && (i_ < num_blocks), "invalid sub-matrix index");
        return TupleElement<i_, First_, Rest_...>::get(*this);
      }
 
      template<int i_>
      typename TupleElement<i_, First_, Rest_...>::Type const& at() const
      {
        static_assert((0 <= i_) && (i_ < num_blocks), "invalid sub-matrix index");
        return TupleElement<i_, First_, Rest_...>::get(*this);
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index rows() const
      {
        const Index rows_f = _first.template rows<perspective_>();
        const Index rows_r = _rest.template rows<perspective_>();
        XASSERTM(rows_f == rows_r, "row count mismatch");
        return rows_f;
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index columns() const
      {
        return _first.template columns<perspective_>() + _rest.template columns<perspective_>();
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index used_elements() const
      {
        return _first.template used_elements<perspective_>() + _rest.template used_elements<perspective_>();
      }
 
      std::size_t bytes() const
      {
        return _first.bytes() + _rest.bytes();
      }
 
      static String name()
      {
        return String("TupleMatrixRow<") + sub_name_list() + ">";
      }
 
      void format(DataType value = DataType(0))
      {
        first().format(value);
        rest().format(value);
      }
 
      void format(Random & rng, DataType min, DataType max)
      {
        first().format(rng, min, max);
        rest().format(rng, min, max);
      }
 
      void clear()
      {
        first().clear();
        rest().clear();
      }
 
      VectorTypeL create_vector_l() const
      {
        return first().create_vector_l();
      }
 
      VectorTypeR create_vector_r() const
      {
        return VectorTypeR(first().create_vector_r(), rest().create_vector_r());
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x) const
      {
        first().apply(r, x.first());
        rest().apply(r, x.rest(), r, DataType(1));
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x) const
      {
        first().apply_transposed(r.first(), x);
        rest().apply_transposed(r.rest(), x);
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x, const VectorTypeL& y, DataType alpha = DataType(1)) const
      {
        first().apply(r, x.first(), y, alpha);
        rest().apply(r, x.rest(), r, alpha);
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x, const VectorTypeR& y, DataType alpha = DataType(1)) const
      {
        first().apply_transposed(r.first(), x, y.first(), alpha);
        rest().apply_transposed(r.rest(), x, y.rest(), alpha);
      }
 
      Index get_length_of_line(const Index row) const
      {
        return first().get_length_of_line(row) + rest().get_length_of_line(row);
      }
 
      void set_line(const Index row, DataType * const pval_set, IndexType * const pcol_set, const Index col_start, const Index stride = 1) const
      {
        const Index first_length(first().get_length_of_line(row));
 
        first().set_line(row, pval_set, pcol_set, col_start, stride);
        rest().set_line(row, pval_set + stride * first_length, pcol_set + stride * first_length, col_start + first().template columns<Perspective::pod>(), stride);
      }
 
      void set_line_reverse(const Index row, DataType * const pval_set, const Index stride = 1)
      {
        const Index first_length(first().get_length_of_line(row));
 
        first().set_line_reverse(row, pval_set, stride);
        rest().set_line_reverse(row, pval_set + stride * first_length, stride);
      }
 
      Index row_degree(const Index row) const
      {
        return first().row_degree(row) + rest().row_degree(row);
      }
 
      template<typename IT2_>
      Index get_row_col_indices(const Index row, IT2_* const pcol_idx, const IT2_ col_offset) const
      {
        const Index first_cols = first().template columns<Perspective::pod>();
        const Index first_degree = first().get_row_col_indices(row, pcol_idx, col_offset);
        return first_degree + rest().get_row_col_indices(row,  pcol_idx + first_degree, col_offset + IT2_(first_cols));
      }
 
      template<typename DT2_>
      Index get_row_values(const Index row, DT2_ * const pvals) const
      {
        const Index first_degree = first().get_row_values(row, pvals);
        return first_degree + rest().get_row_values(row,  pvals + first_degree);
      }
 
      template<typename DT2_>
      Index set_row_values(const Index row, const DT2_ * const pvals)
      {
        const Index first_degree = first().set_row_values(row, pvals);
        return first_degree + rest().set_row_values(row,  pvals + first_degree);
      }
 
      std::uint64_t get_checkpoint_size(SerialConfig& config)
      {
        return sizeof(std::uint64_t) + _first.get_checkpoint_size(config) + _rest.get_checkpoint_size(config); //sizeof(std::uint64_t) bits needed to store lenght of checkpointed _first
      }
 
      void restore_from_checkpoint_data(std::vector<char> & data)
      {
        std::uint64_t isize = *(std::uint64_t*) data.data(); //get size of checkpointed _first
        std::vector<char>::iterator start = std::begin(data) + sizeof(std::uint64_t);
        std::vector<char>::iterator last_of_first = std::begin(data) + sizeof(std::uint64_t) + (int) isize;
        std::vector<char> buffer_first(start, last_of_first);
        _first.restore_from_checkpoint_data(buffer_first);
 
        data.erase(std::begin(data), last_of_first);
        _rest.restore_from_checkpoint_data(data);
      }
 
      std::uint64_t set_checkpoint_data(std::vector<char>& data, SerialConfig& config)
      {
        std::size_t old_size = data.size();
        data.insert(std::end(data), sizeof(std::uint64_t), 0); //add placeholder
        std::uint64_t ireal_size = _first.set_checkpoint_data(data, config); //add data of _first to the overall checkpoint and save its size
        char* csize = reinterpret_cast<char*>(&ireal_size);
        for(std::size_t i(0) ; i < sizeof(std::uint64_t) ; ++i)  //overwrite the guessed datalength
        {
          data[old_size + i] = csize[i];
        }
 
        return sizeof(std::uint64_t) + ireal_size + _rest.set_checkpoint_data(data, config); //generate and add checkpoint data for the _rest
      }
 
      template<typename OtherFirst_, typename... OtherRest_>
      void convert(const TupleMatrixRow<OtherFirst_, OtherRest_...>& other)
      {
        first().convert(other.first());
        rest().convert(other.rest());
      }
 
      template<typename OtherFirst_, typename... OtherRest_>
      void convert_reverse(TupleMatrixRow<OtherFirst_, OtherRest_...>& other) const
      {
        first().convert_reverse(other.first());
        rest().convert_reverse(other.rest());
      }
    }; // template class TupleMatrixRow
 
    // specialization for a single sub-matrix
    template<typename First_>
    class TupleMatrixRow<First_>
    {
    private:
      template<typename,typename...>
      friend class TupleMatrixRow;
 
    public:
      static constexpr int num_blocks = 1;
 
      typedef typename First_::DataType DataType;
      typedef typename First_::IndexType IndexType;
 
      template <typename DT2_ = DataType, typename IT2_ = IndexType>
      using ContainerType = TupleMatrixRow<typename First_::template ContainerType<DT2_, IT2_> >;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
      typedef typename First_::VectorTypeL VectorTypeL;
      typedef TupleVector<typename First_::VectorTypeR> VectorTypeR;
 
    protected:
      First_ _first;
 
      static String sub_name_list()
      {
        return First_::name();
      }
 
    public:
      TupleMatrixRow()
      {
      }
 
      explicit TupleMatrixRow(First_&& the_first) :
        _first(std::forward<First_>(the_first))
      {
      }
 
      TupleMatrixRow(TupleMatrixRow&& other) :
        _first(std::forward<First_>(other._first))
      {
      }
 
      TupleMatrixRow& operator=(TupleMatrixRow&& other)
      {
        if(this != &other)
        {
          _first = std::forward<First_>(other._first);
        }
        return *this;
      }
 
      TupleMatrixRow(const TupleMatrixRow&) = delete;
      TupleMatrixRow& operator=(const TupleMatrixRow&) = delete;
 
      TupleMatrixRow clone(LAFEM::CloneMode mode = LAFEM::CloneMode::Weak) const
      {
        return TupleMatrixRow(_first.clone(mode));
      }
 
      void clone(const TupleMatrixRow& other, CloneMode clone_mode)
      {
        _first.clone(other._first, clone_mode);
      }
 
      void clone(const TupleMatrixRow& other)
      {
        _first.clone(other._first);
      }
 
      First_& first()
      {
        return _first;
      }
 
      const First_& first() const
      {
        return _first;
      }
 
      template<int i_>
      First_& at()
      {
        static_assert(i_ == 0, "invalid sub-matrix index");
        return first();
      }
 
      template<int i_>
      const First_& at() const
      {
        static_assert(i_ == 0, "invalid sub-matrix index");
        return first();
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index rows() const
      {
        return _first.template rows<perspective_>();
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index columns() const
      {
        return _first.template columns<perspective_>();
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index used_elements() const
      {
        return _first.template used_elements<perspective_>();
      }
 
      std::size_t bytes() const
      {
        return _first.bytes();
      }
 
      static String name()
      {
        return String("TupleMatrixRow<") + sub_name_list() + ">";
      }
 
      void format(DataType value = DataType(0))
      {
        first().format(value);
      }
 
      void format(Random & rng, DataType min, DataType max)
      {
        first().format(rng, min, max);
      }
 
      void clear()
      {
        first().clear();
      }
 
      VectorTypeL create_vector_l() const
      {
        return first().create_vector_l();
      }
 
      VectorTypeR create_vector_r() const
      {
        return VectorTypeR(first().create_vector_r());
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x) const
      {
        first().apply(r, x.first());
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x, const VectorTypeL& y, DataType alpha = DataType(1)) const
      {
        first().apply(r, x.first(), y, alpha);
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x) const
      {
        first().apply_transposed(r.first(), x);
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x, const VectorTypeR& y, DataType alpha = DataType(1)) const
      {
        first().apply_transposed(r.first(), x, y.first(), alpha);
      }
 
      Index get_length_of_line(const Index row) const
      {
        return first().get_length_of_line(row);
      }
 
      void set_line(const Index row, DataType * const pval_set, IndexType * const pcol_set, const Index col_start, const Index stride = 1) const
      {
        first().set_line(row, pval_set, pcol_set, col_start, stride);
      }
 
      void set_line_reverse(const Index row, DataType * const pval_set, const Index stride = 1)
      {
        first().set_line_reverse(row, pval_set, stride);
      }
 
      Index row_degree(const Index row) const
      {
        return first().row_degree(row);
      }
 
      template<typename IT2_>
      Index get_row_col_indices(const Index row, IT2_* const pcol_idx, const IT2_ col_offset) const
      {
        return first().get_row_col_indices(row, pcol_idx, col_offset);
      }
 
      template<typename DT2_>
      Index get_row_values(const Index row, DT2_ * const pvals) const
      {
        return first().get_row_values(row, pvals);
      }
 
      template<typename DT2_>
      Index set_row_values(const Index row, const DT2_ * const pvals)
      {
        return first().set_row_values(row, pvals);
      }
 
      std::uint64_t get_checkpoint_size(SerialConfig& config)
      {
        return _first.get_checkpoint_size(config);
      }
 
      void restore_from_checkpoint_data(std::vector<char> & data)
      {
        _first.restore_from_checkpoint_data(data);
      }
 
      std::uint64_t set_checkpoint_data(std::vector<char>& data, SerialConfig& config)
      {
        return _first.set_checkpoint_data(data, config);
      }
 
      template<typename OtherFirst_>
      void convert(const TupleMatrixRow<OtherFirst_>& other)
      {
        first().convert(other.first());
      }
 
      template<typename OtherFirst_>
      void convert_reverse(TupleMatrixRow<OtherFirst_>& other) const
      {
        first().convert_reverse(other.first());
      }
    }; // template class TupleMatrixRow
 
    namespace Intern
    {
      // helper class: checks whether 'T_' is an instance of 'TupleMatrixRow<...>'
      template<typename T_>
      struct IsTMR_
      {
        static constexpr bool value = false;
      };
 
      template<typename... T_>
      struct IsTMR_<TupleMatrixRow<T_...>>
      {
        static constexpr bool value = true;
      };
    } // namespace Intern
 
    template<int i_, int j_, typename FirstRow_, typename... RestRows_>
    class TupleMatrixElement
    {
    public:
      static_assert(i_ > 0, "invalid row index"); // i_=0 is specialized below
      static_assert(j_ >= 0, "invalid column index");
 
      typedef typename TupleMatrixElement<i_-1, j_, RestRows_...>::Type Type;
 
      template<typename Meta_>
      static Type& get(Meta_& meta)
      {
        return TupleMatrixElement<i_-1, j_, RestRows_...>::get(meta.rest());
      }
 
      template<typename Meta_>
      static const Type& get(const Meta_& meta)
      {
        return TupleMatrixElement<i_-1, j_, RestRows_...>::get(meta.rest());
      }
    };
 
    // specialization for i_ = 0;
    // Note: This specialization explicitly assumes that 'FirstRow_' is an instance of
    // 'TupleMatrixRow', which is ensured by a static_assert in the TupleMatrix class.
    template<int j_, typename First_, typename... Rest_, typename... RestRows_>
    class TupleMatrixElement<0, j_, TupleMatrixRow<First_, Rest_...>, RestRows_... >
    {
    public:
      static_assert(j_ >= 0, "invalid column index");
 
      typedef typename TupleElement<j_, First_, Rest_...>::Type Type;
 
      template<typename Meta_>
      static Type& get(Meta_& meta)
      {
        return TupleElement<j_, First_, Rest_...>::get(meta.first());
      }
 
      template<typename Meta_>
      static const Type& get(const Meta_& meta)
      {
        return TupleElement<j_, First_, Rest_...>::get(meta.first());
      }
    };
 
    template<typename FirstRow_, typename... RestRows_>
    class TupleMatrix
    {
    private:
      template<typename,typename...>
      friend class TupleMatrix;
 
      typedef TupleMatrix<RestRows_...> RestClass;
 
      // sanity check: ensure that 'FirstRow_' is an instance of TupleMatrixRow
      // note: this is essential for the implementation of TupleMatrixElement!
      static_assert(Intern::IsTMR_<FirstRow_>::value,
        "template parameter of TupleMatrix must be TupleMatrixRow<...>");
 
    public:
      static constexpr int num_row_blocks = TupleMatrix<RestRows_...>::num_row_blocks + 1;
      static constexpr int num_col_blocks = FirstRow_::num_blocks;
 
      // sanity check: number of column blocks must match for all rows
      static_assert(num_col_blocks == RestClass::num_col_blocks, "column blocks mismatch");
 
      typedef typename FirstRow_::DataType DataType;
      typedef typename FirstRow_::IndexType IndexType;
 
      // ensure that all sub-matrix have the same data-type
      static_assert(std::is_same<DataType, typename RestClass::DataType>::value,
                    "sub-matrices have different data-types");
      static_assert(std::is_same<IndexType, typename RestClass::IndexType>::value,
                    "sub-matrices have different index-types");
 
      template <typename DT2_ = DataType, typename IT2_ = IndexType>
      using ContainerType = TupleMatrix<
        typename FirstRow_::template ContainerType<DT2_, IT2_>,
        typename RestRows_::template ContainerType<DT2_, IT2_>...>;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
      typedef TupleVector<typename FirstRow_::VectorTypeL, typename RestRows_::VectorTypeL...> VectorTypeL;
      typedef typename FirstRow_::VectorTypeR VectorTypeR;
 
      static_assert(std::is_same<VectorTypeR, typename RestClass::VectorTypeR>::value, "vector type mismatch");
 
    protected:
      FirstRow_ _first;
      RestClass _rest;
 
      static String sub_name_list()
      {
        return FirstRow_::name() + "," + RestClass::sub_name_list();
      }
 
    public:
      TupleMatrix()
      {
      }
 
      explicit TupleMatrix(FirstRow_&& the_first, RestClass&& the_rest) :
        _first(std::forward<FirstRow_>(the_first)),
        _rest(std::forward<RestClass>(the_rest))
      {
      }
 
      explicit TupleMatrix(FirstRow_&& the_first, RestRows_&&... the_rest) :
        _first(std::forward<FirstRow_>(the_first)),
        _rest(std::forward<RestRows_>(the_rest)...)
      {
      }
 
      TupleMatrix(TupleMatrix&& other) :
        _first(std::forward<FirstRow_>(other._first)),
        _rest(std::forward<RestClass>(other._rest))
      {
      }
 
      TupleMatrix& operator=(TupleMatrix&& other)
      {
        if(this != &other)
        {
          _first = std::forward<FirstRow_>(other._first);
          _rest = std::forward<RestClass>(other._rest);
        }
        return *this;
      }
 
      TupleMatrix(const TupleMatrix&) = delete;
      TupleMatrix& operator=(const TupleMatrix&) = delete;
 
      TupleMatrix clone(LAFEM::CloneMode mode = LAFEM::CloneMode::Weak) const
      {
        return TupleMatrix(_first.clone(mode), _rest.clone(mode));
      }
 
      void clone(const TupleMatrix& other, CloneMode clone_mode)
      {
        _first.clone(other._first, clone_mode);
        _rest.clone(other._rest, clone_mode);
      }
 
      void clone(const TupleMatrix& other)
      {
        _first.clone(other._first);
        _rest.clone(other._rest);
      }
 
      FirstRow_& first()
      {
        return _first;
      }
 
      const FirstRow_& first() const
      {
        return _first;
      }
 
      TupleMatrix<RestRows_...>& rest()
      {
        return _rest;
      }
 
      const TupleMatrix<RestRows_...>& rest() const
      {
        return _rest;
      }
 
      template<int i_, int j_>
      typename TupleMatrixElement<i_, j_, FirstRow_, RestRows_...>::Type& at()
      {
        static_assert((0 <= i_) && (i_ < num_row_blocks), "invalid sub-matrix row index");
        static_assert((0 <= j_) && (j_ < num_col_blocks), "invalid sub-matrix column index");
        return TupleMatrixElement<i_, j_, FirstRow_, RestRows_...>::get(*this);
      }
 
      template<int i_, int j_>
      const typename TupleMatrixElement<i_, j_, FirstRow_, RestRows_...>::Type& at() const
      {
        static_assert((0 <= i_) && (i_ < num_row_blocks), "invalid sub-matrix row index");
        static_assert((0 <= j_) && (j_ < num_col_blocks), "invalid sub-matrix column index");
        return TupleMatrixElement<i_, j_, FirstRow_, RestRows_...>::get(*this);
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index rows() const
      {
        return _first.template rows<perspective_>() + _rest.template rows<perspective_>();
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index columns() const
      {
        const Index cols_f = _first.template columns<perspective_>();
        const Index cols_r = _rest.template columns<perspective_>();
        XASSERTM(cols_f == cols_r, "column count mismatch");
        return cols_f;
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index used_elements() const
      {
        return _first.template used_elements<perspective_>() + _rest.template used_elements<perspective_>();
      }
 
      std::size_t bytes() const
      {
        return _first.bytes() + _rest.bytes();
      }
 
      static String name()
      {
        return String("TupleMatrix<") + sub_name_list() + ">";
      }
 
      void format(DataType value = DataType(0))
      {
        first().format(value);
        rest().format(value);
      }
 
      void format(Random & rng, DataType min, DataType max)
      {
        first().format(rng, min, max);
        rest().format(rng, min, max);
      }
 
      void clear()
      {
        first().clear();
        rest().clear();
      }
 
      VectorTypeL create_vector_l() const
      {
        return VectorTypeL(first().create_vector_l(), rest().create_vector_l());
      }
 
      VectorTypeR create_vector_r() const
      {
        return first().create_vector_r();
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x) const
      {
        first().apply(r.first(), x);
        rest().apply(r.rest(), x);
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x) const
      {
        first().apply_transposed(r, x.first());
        rest().apply_transposed(r, x.rest(), r, DataType(1));
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x, const VectorTypeL& y, DataType alpha = DataType(1)) const
      {
        first().apply(r.first(), x, y.first(), alpha);
        rest().apply(r.rest(), x, y.rest(), alpha);
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x, const VectorTypeR& y, DataType alpha = DataType(1)) const
      {
        first().apply_transposed(r, x.first(), y, alpha);
        rest().apply_transposed(r, x.rest(), r, alpha);
      }
 
      Index get_length_of_line(const Index row) const
      {
        const Index first_rows = first().template rows<Perspective::pod>();
        if(row < first_rows)
          return first().get_length_of_line(row);
        else
          return rest().get_length_of_line(row - first_rows);
      }
 
      void set_line(const Index row, DataType * const pval_set, IndexType * const pcol_set, const Index col_start, const Index stride = 1) const
      {
        const Index first_rows = first().template rows<Perspective::pod>();
        if(row < first_rows)
          first().set_line(row, pval_set, pcol_set, col_start, stride);
        else
          rest().set_line(row - first_rows, pval_set, pcol_set, col_start, stride);
      }
 
      void set_line_reverse(const Index row, DataType * const pval_set, const Index stride = 1)
      {
        const Index first_rows = first().template rows<Perspective::pod>();
        if(row < first_rows)
          first().set_line_reverse(row, pval_set, stride);
        else
          rest().set_line_reverse(row - first_rows, pval_set, stride);
      }
 
      Index row_degree(const Index row) const
      {
        const Index first_rows = first().template rows<Perspective::pod>();
        if(row < first_rows)
          return first().row_degree(row);
        else
          return rest().row_degree(row - first_rows);
      }
 
      template<typename IT2_>
      Index get_row_col_indices(const Index row, IT2_* const pcol_idx, const IT2_ col_offset) const
      {
        const Index first_rows = first().template rows<Perspective::pod>();
        if(row < first_rows)
          return first().get_row_col_indices(row, pcol_idx, col_offset);
        else
          return rest().get_row_col_indices(row - first_rows, pcol_idx, col_offset);
      }
 
      template<typename DT2_>
      Index get_row_values(const Index row, DT2_ * const pvals) const
      {
        const Index first_rows = first().template rows<Perspective::pod>();
        if(row < first_rows)
          return first().get_row_values(row, pvals);
        else
          return rest().get_row_values(row - first_rows, pvals);
      }
 
      template<typename DT2_>
      Index set_row_values(const Index row, const DT2_ * const pvals)
      {
        const Index first_rows = first().template rows<Perspective::pod>();
        if(row < first_rows)
          return first().set_row_values(row, pvals);
        else
          return rest().set_row_values(row - first_rows, pvals);
      }
 
      std::uint64_t get_checkpoint_size(SerialConfig& config)
      {
        return sizeof(std::uint64_t) + _first.get_checkpoint_size(config) + _rest.get_checkpoint_size(config); //sizeof(std::uint64_t) bits needed to store lenght of checkpointed _first
      }
 
      void restore_from_checkpoint_data(std::vector<char> & data)
      {
        std::uint64_t isize = *(std::uint64_t*) data.data(); //get size of checkpointed _first
        std::vector<char>::iterator start = std::begin(data) + sizeof(std::uint64_t);
        std::vector<char>::iterator last_of_first = std::begin(data) + sizeof(std::uint64_t) + (int) isize;
        std::vector<char> buffer_first(start, last_of_first);
        _first.restore_from_checkpoint_data(buffer_first);
 
        data.erase(std::begin(data), last_of_first);
        _rest.restore_from_checkpoint_data(data);
      }
 
      std::uint64_t set_checkpoint_data(std::vector<char>& data, SerialConfig& config)
      {
        std::size_t old_size = data.size();
        data.insert(std::end(data), sizeof(std::uint64_t), 0); //add placeholder
        std::uint64_t ireal_size = _first.set_checkpoint_data(data, config); //add data of _first to the overall checkpoint and save its size
        char* csize = reinterpret_cast<char*>(&ireal_size);
        for(std::size_t i(0) ; i < sizeof(std::uint64_t) ; ++i)  //overwrite the guessed datalength
        {
          data[old_size +i] = csize[i];
        }
        return sizeof(std::uint64_t) + ireal_size + _rest.set_checkpoint_data(data, config); //generate and add checkpoint data for the _rest
      }
 
      template<typename OtherFirstRow_, typename... OtherRestRows_>
      void convert(const TupleMatrix<OtherFirstRow_, OtherRestRows_...>& other)
      {
        first().convert(other.first());
        rest().convert(other.rest());
      }
 
      template<typename OtherFirstRow_, typename... OtherRestRows_>
      void convert_reverse(TupleMatrix<OtherFirstRow_, OtherRestRows_...>& other) const
      {
        first().convert_reverse(other.first());
        rest().convert_reverse(other.rest());
      }
    }; // template class TupleMatrix
 
    // specialization for 1 row
    template<typename FirstRow_>
    class TupleMatrix<FirstRow_>
    {
    private:
      template<typename,typename...>
      friend class TupleMatrix;
 
      // sanity check: ensure that 'FirstRow_' is an instance of TupleMatrixRow
      static_assert(Intern::IsTMR_<FirstRow_>::value,
        "template parameter of TupleMatrix must be TupleMatrixRow<...>");
 
    public:
      static constexpr int num_blocks = 1;
 
      static constexpr int num_row_blocks = num_blocks;
      static constexpr int num_col_blocks = FirstRow_::num_blocks;
 
 
      typedef typename FirstRow_::DataType DataType;
      typedef typename FirstRow_::IndexType IndexType;
 
      template < typename DT2_ = DataType, typename IT2_ = IndexType>
      using ContainerType = TupleMatrix<typename FirstRow_::template ContainerType<DT2_, IT2_> >;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
      typedef TupleVector<typename FirstRow_::VectorTypeL> VectorTypeL;
      typedef typename FirstRow_::VectorTypeR VectorTypeR;
 
    protected:
      FirstRow_ _first;
 
      static String sub_name_list()
      {
        return FirstRow_::name();
      }
 
    public:
      TupleMatrix()
      {
      }
 
      explicit TupleMatrix(FirstRow_&& the_first) :
        _first(std::forward<FirstRow_>(the_first))
      {
      }
 
      TupleMatrix(TupleMatrix&& other) :
        _first(std::forward<FirstRow_>(other._first))
      {
      }
 
      TupleMatrix& operator=(TupleMatrix&& other)
      {
        if(this != &other)
        {
          _first = std::forward<FirstRow_>(other._first);
        }
        return *this;
      }
 
      TupleMatrix(const TupleMatrix&) = delete;
      TupleMatrix& operator=(const TupleMatrix&) = delete;
 
      TupleMatrix clone(LAFEM::CloneMode mode = LAFEM::CloneMode::Weak) const
      {
        return TupleMatrix(_first.clone(mode));
      }
 
      void clone(const TupleMatrix& other, CloneMode clone_mode)
      {
        _first.clone(other._first, clone_mode);
      }
 
      void clone(const TupleMatrix& other)
      {
        _first.clone(other._first);
      }
 
      FirstRow_& first()
      {
        return _first;
      }
 
      const FirstRow_& first() const
      {
        return _first;
      }
 
      template<int i_, int j_>
      typename TupleMatrixElement<i_, j_, FirstRow_>::Type& at()
      {
        static_assert(0 == i_, "invalid sub-matrix row index");
        static_assert((0 <= j_) && (j_ < num_col_blocks), "invalid sub-matrix column index");
        return TupleMatrixElement<i_, j_, FirstRow_>::get(*this);
      }
 
      template<int i_, int j_>
      const typename TupleMatrixElement<i_, j_, FirstRow_>::Type& at() const
      {
        static_assert(0 == i_, "invalid sub-matrix row index");
        static_assert((0 <= j_) && (j_ < num_col_blocks), "invalid sub-matrix column index");
        return TupleMatrixElement<i_, j_, FirstRow_>::get(*this);
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index rows() const
      {
        return _first.template rows<perspective_>();
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index columns() const
      {
        return _first.template columns<perspective_>();
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index used_elements() const
      {
        return _first.template used_elements<perspective_>();
      }
 
      std::size_t bytes() const
      {
        return _first.bytes();
      }
 
      static String name()
      {
        return String("TupleMatrix<") + sub_name_list() + ">";
      }
 
      void format(DataType value = DataType(0))
      {
        first().format(value);
      }
 
      void format(Random & rng, DataType min, DataType max)
      {
        first().format(rng, min, max);
      }
 
      void clear()
      {
        first().clear();
      }
 
      VectorTypeL create_vector_l() const
      {
        return VectorTypeL(first().create_vector_l());
      }
 
      VectorTypeR create_vector_r() const
      {
        return first().create_vector_r();
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x) const
      {
        first().apply(r.first(), x);
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x, const VectorTypeL& y, DataType alpha = DataType(1)) const
      {
        first().apply(r.first(), x, y.first(), alpha);
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x) const
      {
        first().apply_transposed(r, x.first());
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x, const VectorTypeR& y, DataType alpha = DataType(1)) const
      {
        first().apply(r, x.first(), y, alpha);
      }
 
      Index get_length_of_line(const Index row) const
      {
        return first().get_length_of_line(row);
      }
 
      void set_line(const Index row, DataType * const pval_set, IndexType * const pcol_set, const Index col_start, const Index stride = 1) const
      {
        first().set_line(row, pval_set, pcol_set, col_start, stride);
      }
 
      void set_line_reverse(const Index row, DataType * const pval_set, const Index stride = 1)
      {
        first().set_line_reverse(row, pval_set, stride);
      }
 
      Index row_degree(const Index row) const
      {
        return first().row_degree(row);
      }
 
      template<typename IT2_>
      Index get_row_col_indices(const Index row, IT2_* const pcol_idx, const IT2_ col_offset) const
      {
        return first().get_row_col_indices(row, pcol_idx, col_offset);
      }
 
      template<typename DT2_>
      Index get_row_values(const Index row, DT2_ * const pvals) const
      {
        return first().get_row_values(row, pvals);
      }
 
      template<typename DT2_>
      Index set_row_values(const Index row, const DT2_ * const pvals)
      {
        return first().set_row_values(row, pvals);
      }
 
      std::uint64_t get_checkpoint_size(SerialConfig& config)
      {
        return _first.get_checkpoint_size(config);
      }
 
      void restore_from_checkpoint_data(std::vector<char> & data)
      {
        _first.restore_from_checkpoint_data(data);
      }
 
      std::uint64_t set_checkpoint_data(std::vector<char>& data, SerialConfig& config)
      {
        return _first.set_checkpoint_data(data, config);
      }
 
      template<typename OtherFirstRow_>
      void convert(const TupleMatrix<OtherFirstRow_>& other)
      {
        first().convert(other.first());
      }
 
      template<typename OtherFirstRow_>
      void convert_reverse(TupleMatrix<OtherFirstRow_>& other) const
      {
        first().convert_reverse(other.first());
      }
    }; // template class TupleMatrix
  } // namespace LAFEM
} // namespace FEAT