power_diag_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
 
// includes, FEAT
#include <kernel/lafem/power_vector.hpp>
#include <kernel/lafem/sparse_layout.hpp>
#include <kernel/lafem/meta_element.hpp>
#include <kernel/lafem/container.hpp>
 
 
namespace FEAT
{
  namespace LAFEM
  {
    template<
      typename SubType_,
      int blocks_>
    class PowerDiagMatrix
    {
      // Note: the case = 1 is specialized below
      static_assert(blocks_ > 1, "invalid block size");
 
      // declare this class template as a friend for recursive inheritance
      template<typename, int>
      friend class PowerDiagMatrix;
 
      typedef PowerDiagMatrix<SubType_, blocks_-1> RestClass;
 
    public:
      typedef SubType_ SubMatrixType;
      typedef typename SubMatrixType::DataType DataType;
      typedef typename SubMatrixType::IndexType IndexType;
      static constexpr SparseLayoutId layout_id = SubMatrixType::layout_id;
      typedef PowerVector<typename SubMatrixType::VectorTypeL, blocks_> VectorTypeL;
      typedef PowerVector<typename SubMatrixType::VectorTypeR, blocks_> VectorTypeR;
      template <typename DT2_ = DataType, typename IT2_ = IndexType>
      using ContainerType = PowerDiagMatrix<typename SubType_::template ContainerType<DT2_, IT2_>, blocks_>;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
      static constexpr int num_row_blocks = blocks_;
      static constexpr int num_col_blocks = blocks_;
 
    protected:
      SubMatrixType _first;
      RestClass _rest;
 
      explicit PowerDiagMatrix(SubMatrixType&& the_first, RestClass&& the_rest) :
        _first(std::move(the_first)),
        _rest(std::move(the_rest))
      {
      }
 
    public:
      PowerDiagMatrix()
      {
      }
 
      explicit PowerDiagMatrix(const SparseLayout<IndexType, layout_id>& layout) :
        _first(layout),
        _rest(layout)
      {
      }
 
      PowerDiagMatrix(PowerDiagMatrix&& other) :
        _first(std::move(other._first)),
        _rest(std::move(other._rest))
      {
      }
 
      explicit PowerDiagMatrix(FileMode mode, String filename)
      {
        read_from(mode, filename);
      }
 
      explicit PowerDiagMatrix(FileMode mode, std::istream& file, String directory = "")
      {
        String line;
        do {
          if (file.eof())
            XABORTM("Wrong Input-file");
          std::getline(file, line);
          line.trim_me();
        } while (line.find("%%") == 0 || line == "");
 
        SubMatrixType tmp_first(mode, directory + line);
        _first = std::move(tmp_first);
 
        RestClass tmp_rest(mode, file, directory);
        _rest = std::move(tmp_rest);
      }
 
      void read_from(FileMode mode, String filename)
      {
        String directory;
        auto found = filename.rfind("/");
        if (found != std::string::npos)
        {
          directory = filename.substr(0, found + 1);
        }
 
        std::ifstream file(filename.c_str(), std::ifstream::in);
        if (! file.is_open())
          XABORTM("Unable to open Matrix file " + filename);
 
        String line;
        std::getline(file, line);
        if (line.find("%%MatrixMarket powerdiagmatrix coordinate real general") == String::npos)
          XABORTM("Input-file is not a complatible file");
 
        PowerDiagMatrix other(mode, file, directory);
 
        _first = std::move(other._first);
        _rest = std::move(other._rest);
 
        file.close();
      }
 
      PowerDiagMatrix& operator=(PowerDiagMatrix&& other)
      {
        if(this != &other)
        {
          _first = std::move(other._first);
          _rest = std::move(other._rest);
        }
        return *this;
      }
 
      PowerDiagMatrix(const PowerDiagMatrix&) = delete;
      PowerDiagMatrix& operator=(const PowerDiagMatrix&) = delete;
 
      virtual ~PowerDiagMatrix()
      {
      }
 
      void write_out(FileMode mode, String filename) const
      {
        std::ofstream file(filename.c_str(), std::ofstream::out);
        if (! file.is_open())
          XABORTM("Unable to open Matrix file " + filename);
 
        String suffix, directory;
        auto found = filename.rfind(".");
        if (found != std::string::npos)
        {
          suffix = filename.substr(found);
          filename.erase(found);
        }
        found = filename.rfind("/");
        if (found != std::string::npos)
        {
          directory = filename.substr(0, found + 1);
          filename.erase(0, found + 1);
        }
 
        file << "%%MatrixMarket powerdiagmatrix coordinate real general\n";
        for (Index i(1); i <= blocks_; ++i)
        {
          file << filename << "_pd" << i << suffix << "\n";
        }
 
        file.close();
 
        this->write_out_submatrices(mode, directory, filename, suffix);
      }
 
      void write_out_submatrices(FileMode mode, String directory, String prefix, String suffix, Index length = blocks_) const
      {
        _first.write_out(mode, directory + prefix + "_pd" + stringify(length + 1 - blocks_) + suffix);
        _rest.write_out_submatrices(mode, directory, prefix, suffix, length);
      }
 
      PowerDiagMatrix clone(LAFEM::CloneMode mode = LAFEM::CloneMode::Weak) const
      {
        return PowerDiagMatrix(_first.clone(mode), _rest.clone(mode));
      }
 
      void clone(const PowerDiagMatrix & other, CloneMode clone_mode = CloneMode::Weak)
      {
        (*this)=other.clone(clone_mode);
      }
 
      std::size_t bytes() const
      {
        return _first.bytes() + _rest.bytes();
      }
 
      PowerDiagMatrix transpose() const
      {
        return PowerDiagMatrix(_first.transpose(), _rest.transpose());
      }
 
      template<int i_, int j_>
      SubMatrixType& at()
      {
        static_assert(i_ == j_, "invalid sub-matrix index");
        static_assert((0 <= i_) && (i_ < blocks_), "invalid sub-matrix index");
        return PowerElement<i_, SubMatrixType>::get(*this);
      }
 
      template<int i_, int j_>
      const SubMatrixType& at() const
      {
        static_assert(i_ == j_, "invalid sub-matrix index");
        static_assert((0 <= i_) && (i_ < blocks_), "invalid sub-matrix index");
        return PowerElement<i_, SubMatrixType>::get(*this);
      }
 
      SubMatrixType& get(int i, int j)
      {
        XASSERTM((i == j) && (0 <= i) && (i < blocks_), "invalid sub-matrix index");
        return (i == 0) ? _first : _rest.get(i-1, j-1);
      }
 
      const SubMatrixType& get(int i, int j) const
      {
        XASSERTM((i == j) && (0 <= i) && (i < blocks_), "invalid sub-matrix index");
        return (i == 0) ? _first : _rest.get(i-1, j-1);
      }
 
      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
      }
 
      SubMatrixType& first()
      {
        return _first;
      }
 
      const SubMatrixType& first() const
      {
        return _first;
      }
 
      RestClass& rest()
      {
        return _rest;
      }
 
      const RestClass& rest() const
      {
        return _rest;
      }
 
      int row_blocks() const
      {
        return num_row_blocks;
      }
 
      int col_blocks() const
      {
        return num_col_blocks;
      }
 
      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
      {
        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_>();
      }
 
      static String name()
      {
        return String("PowerDiagMatrix<") + SubMatrixType::name() + "," + stringify(blocks_) + ">";
      }
 
      Index size() const
      {
        return rows() * columns();
      }
 
      void format(DataType value = DataType(0))
      {
        first().format(value);
        rest().format(value);
      }
 
      void clear()
      {
        first().clear();
        rest().clear();
      }
 
      void extract_diag(VectorTypeL& diag) const
      {
        first().extract_diag(diag.first());
        rest().extract_diag(diag.rest());
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x) const
      {
        first().apply(r.first(), x.first());
        rest().apply(r.rest(), x.rest());
      }
 
      void apply(DenseVector<DataType , IndexType>& r, const DenseVector<DataType , IndexType>& 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!");
 
        DenseVector<DataType, IndexType> r_first(r, first().rows(), 0);
        DenseVector<DataType, IndexType> r_rest(r, rest().rows(), first().rows());
 
        DenseVector<DataType, IndexType> x_first(x, first().columns(), 0);
        DenseVector<DataType, IndexType> x_rest(x, rest().columns(), first().columns());
 
        first().apply(r_first, x_first);
        rest().apply(r_rest, x_rest);
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x) const
      {
        first().apply_transposed(r.first(), x.first());
        rest().apply_transposed(r.rest(), x.rest());
      }
 
      void apply_transposed(DenseVector<DataType , IndexType>& r, const DenseVector<DataType , IndexType>& 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!");
 
        DenseVector<DataType, IndexType> r_first(r, first().columns(), 0);
        DenseVector<DataType, IndexType> r_rest(r, rest().columns(), first().columns());
 
        DenseVector<DataType, IndexType> x_first(x, first().rows(), 0);
        DenseVector<DataType, IndexType> x_rest(x, rest().rows(), first().rows());
 
        first().apply_transposed(r_first, x_first);
        rest().apply_transposed(r_rest, x_rest);
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x, const VectorTypeL& y, DataType alpha = DataType(1)) const
      {
        first().apply(r.first(), x.first(), y.first(), alpha);
        rest().apply(r.rest(), x.rest(), y.rest(), alpha);
      }
 
      void apply(DenseVector<DataType, IndexType>& r, const DenseVector<DataType, IndexType>& x,
                 const DenseVector<DataType, IndexType>& y, DataType alpha = DataType(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!");
 
        DenseVector<DataType, IndexType> r_first(r, first().rows(), 0);
        DenseVector<DataType, IndexType> r_rest(r, rest().rows(), first().rows());
 
        DenseVector<DataType, IndexType> x_first(x, first().columns(), 0);
        DenseVector<DataType, IndexType> x_rest(x, rest().columns(), first().columns());
 
        DenseVector<DataType, IndexType> y_first(y, first().rows(), 0);
        DenseVector<DataType, IndexType> y_rest(y, rest().rows(), first().rows());
 
        first().apply(r_first, x_first, y_first, alpha);
        rest().apply(r_rest, x_rest, y_rest, alpha);
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x, const VectorTypeR& y, DataType alpha = DataType(1)) const
      {
        first().apply_transposed(r.first(), x.first(), y.first(), alpha);
        rest().apply_transposed(r.rest(), x.rest(), y.rest(), alpha);
      }
 
      void apply_transposed(DenseVector<DataType, IndexType>& r, const DenseVector<DataType, IndexType>& x,
        const DenseVector<DataType, IndexType>& y, DataType alpha = DataType(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!");
 
        DenseVector<DataType, IndexType> r_first(r, first().columns(), 0);
        DenseVector<DataType, IndexType> r_rest(r, rest().columns(), first().columns());
 
        DenseVector<DataType, IndexType> x_first(x, first().rows(), 0);
        DenseVector<DataType, IndexType> x_rest(x, rest().rows(), first().rows());
 
        DenseVector<DataType, IndexType> y_first(y, first().columns(), 0);
        DenseVector<DataType, IndexType> y_rest(y, rest().columns(), first().columns());
 
        first().apply_transposed(r_first, x_first, y_first, alpha);
        rest().apply_transposed(r_rest, x_rest, y_rest, alpha);
      }
 
      DataType max_rel_diff(const PowerDiagMatrix& x) const
      {
        DataType max_rel_diff = Math::max(this->first().max_rel_diff(x.first()), this->rest().max_rel_diff(x.rest()));
        return max_rel_diff;
      }
 
      bool same_layout(const PowerDiagMatrix& x) const
      {
        return (this->name() == x.name()) && (this->first().same_layout(x.first())) && (this->rest().same_layout(x.rest()));
      }
 
      VectorTypeL create_vector_l() const
      {
        return VectorTypeL(first().create_vector_l(), rest().create_vector_l());
      }
 
      VectorTypeR create_vector_r() const
      {
        return VectorTypeR(first().create_vector_r(), rest().create_vector_r());
      }
 
      void scale_rows(const PowerDiagMatrix& a, const VectorTypeL& w)
      {
        first().scale_rows(a.first(), w.first());
        rest().scale_rows(a.rest(), w.rest());
      }
 
      void scale_cols(const PowerDiagMatrix& a, const VectorTypeR& w)
      {
        first().scale_cols(a.first(), w.first());
        rest().scale_cols(a.rest(), w.rest());
      }
 
      Index get_length_of_line(const Index row) const
      {
        const Index brows(this->first().template rows<Perspective::pod>());
 
        if (row < brows)
        {
          return this->first().get_length_of_line(row);
        }
        else
        {
          return this->rest().get_length_of_line(row - brows);
        }
      }
 
      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 brows(this->first().template rows<Perspective::pod>());
        const Index bcolumns(this->first().template columns<Perspective::pod>());
 
        if (row < brows)
        {
          this->first().set_line(row, pval_set, pcol_set, col_start, stride);
        }
        else
        {
          this->rest().set_line(row - brows, pval_set, pcol_set, col_start + bcolumns, stride);
        }
      }
 
      void set_line_reverse(const Index row, const DataType * const pval_set, const Index stride = 1)
      {
        const Index brows(this->first().template rows<Perspective::pod>());
 
        if (row < brows)
        {
          this->first().set_line_reverse(row, pval_set, stride);
        }
        else
        {
          this->rest().set_line_reverse(row - brows, 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>();
        const Index first_cols = first().template columns<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 + IT2_(first_cols));
      }
 
      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);
      }
 
 
      template <typename SubType2_>
      void convert(const PowerDiagMatrix<SubType2_, blocks_> & other)
      {
        this->first().convert(other.first());
        this->rest().convert(other.rest());
      }
 
      template <typename SubType2_>
      void convert_reverse(PowerDiagMatrix<SubType2_, blocks_> & other) const
      {
        this->first().convert_reverse(other.first());
        this->rest().convert_reverse(other.rest());
      }
    };
 
    template<typename SubType_>
    class PowerDiagMatrix<SubType_, 1>
    {
      template<typename, int>
      friend class PowerDiagMatrix;
 
    public:
      typedef SubType_ SubMatrixType;
      typedef typename SubMatrixType::DataType DataType;
      typedef typename SubMatrixType::IndexType IndexType;
      static constexpr SparseLayoutId layout_id = SubMatrixType::layout_id;
      typedef PowerVector<typename SubMatrixType::VectorTypeL, 1> VectorTypeL;
      typedef PowerVector<typename SubMatrixType::VectorTypeR, 1> VectorTypeR;
      template <typename DT2_ = DataType, typename IT2_ = IndexType>
      using ContainerType = PowerDiagMatrix<typename SubType_::template ContainerType<DT2_, IT2_>, 1>;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
      static constexpr int num_row_blocks = 1;
      static constexpr int num_col_blocks = 1;
 
    protected:
      SubMatrixType _first;
 
      explicit PowerDiagMatrix(SubMatrixType&& the_first) :
        _first(std::move(the_first))
      {
      }
 
    public:
      PowerDiagMatrix()
      {
      }
 
      explicit PowerDiagMatrix(const SparseLayout<IndexType, layout_id>& layout) :
        _first(layout)
      {
      }
 
      PowerDiagMatrix(PowerDiagMatrix&& other) :
        _first(std::move(other._first))
      {
      }
 
      PowerDiagMatrix& operator=(PowerDiagMatrix&& other)
      {
        if(this != &other)
        {
          _first = std::move(other._first);
        }
        return *this;
      }
 
      explicit PowerDiagMatrix(FileMode mode, const String& filename)
      {
        read_from(mode, filename);
      }
 
      explicit PowerDiagMatrix(FileMode mode, std::istream& file, const String& directory = "")
      {
        String line;
        do {
          if (file.eof())
            XABORTM("Wrong Input-file");
          std::getline(file, line);
          line.trim_me();
        } while (line.find("%%") == 0 || line == "");
 
        SubMatrixType tmp_first(mode, directory + line);
        _first = std::move(tmp_first);
      }
 
      void read_from(FileMode mode, const String& filename)
      {
        String directory;
        auto found = filename.rfind("/");
        if (found != std::string::npos)
        {
          directory = filename.substr(0, found + 1);
        }
 
        std::ifstream file(filename.c_str(), std::ifstream::in);
        if (! file.is_open())
          XABORTM("Unable to open Matrix file " + filename);
 
        PowerDiagMatrix other(mode, file, directory);
 
        _first = std::move(other._first);
 
        file.close();
      }
 
      PowerDiagMatrix(const PowerDiagMatrix&) = delete;
      PowerDiagMatrix& operator=(const PowerDiagMatrix&) = delete;
 
      virtual ~PowerDiagMatrix()
      {
      }
 
      void write_out(FileMode mode, String filename) const
      {
        std::ofstream file(filename.c_str(), std::ofstream::out);
        if (! file.is_open())
          XABORTM("Unable to open Matrix file " + filename);
 
        String suffix, directory;
        auto found = filename.rfind(".");
        if (found != std::string::npos)
        {
          suffix = filename.substr(found);
          filename.erase(found);
        }
        found = filename.rfind("/");
        if (found != std::string::npos)
        {
          directory = filename.substr(0, found + 1);
          filename.erase(0, found + 1);
        }
 
        file << "%%MatrixMarket powerdiagmatrix coordinate real general\n";
        file << filename << "_pd" << 1 << suffix << "\n";
 
        file.close();
 
        this->write_out_submatrices(mode, directory, filename, suffix);
      }
 
      void write_out_submatrices(FileMode mode, const String& directory, const String& prefix, const String& suffix, Index length = 1) const
      {
        _first.write_out(mode, directory + prefix + "_pd" + stringify(length) + suffix);
      }
 
      PowerDiagMatrix clone(LAFEM::CloneMode mode = LAFEM::CloneMode::Weak) const
      {
        return PowerDiagMatrix(_first.clone(mode));
      }
 
      void clone(const PowerDiagMatrix & other, CloneMode clone_mode = CloneMode::Weak)
      {
        (*this)=other.clone(clone_mode);
      }
 
      std::size_t bytes() const
      {
        return _first.bytes();
      }
 
      PowerDiagMatrix transpose() const
      {
        return PowerDiagMatrix(_first.transpose());
      }
 
      template<int i, int j>
      SubMatrixType& at()
      {
        static_assert(i == 0, "invalid sub-matrix index");
        static_assert(j == 0, "invalid sub-matrix index");
        return _first;
      }
 
      template<int i, int j>
      const SubMatrixType& at() const
      {
        static_assert(i == 0, "invalid sub-matrix index");
        static_assert(j == 0, "invalid sub-matrix index");
        return _first;
      }
 
      SubMatrixType& get(int i, int j)
      {
        XASSERTM((i == 0) && (j == 0), "invalid sub-matrix index");
        return _first;
      }
 
      const SubMatrixType& get(int i, int j) const
      {
        XASSERTM((i == 0) && (j == 0), "invalid sub-matrix index");
        return _first;
      }
 
      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);
      }
 
      SubMatrixType& first()
      {
        return _first;
      }
 
      const SubMatrixType& first() const
      {
        return _first;
      }
 
      int row_blocks() const
      {
        return 1;
      }
 
      int col_blocks() const
      {
        return 1;
      }
 
      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_>();
      }
 
      static String name()
      {
        return String("PowerDiagMatrix<") + SubMatrixType::name() + "," + stringify(1) + ">";
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index size() const
      {
        return rows<perspective_>() * columns<perspective_>();
      }
 
      void format(DataType value = DataType(0))
      {
        first().format(value);
      }
 
      void clear()
      {
        first().clear();
      }
 
      void extract_diag(VectorTypeL& diag) const
      {
        first().extract_diag(diag.first());
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x) const
      {
        first().apply(r.first(), x.first());
      }
 
      void apply(DenseVector<DataType, IndexType>& r, const DenseVector<DataType, IndexType>& x) const
      {
        first().apply(r, x);
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x, const VectorTypeL& y, DataType alpha = DataType(1)) const
      {
        first().apply(r.first(), x.first(), y.first(), alpha);
      }
 
      void apply(DenseVector<DataType, IndexType>& r, const DenseVector<DataType, IndexType>& x,
                 const DenseVector<DataType, IndexType>& y, DataType alpha = DataType(1)) const
      {
        first().apply(r, x, y, alpha);
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x) const
      {
        first().apply_transposed(r.first(), x.first());
      }
 
      void apply_transposed(DenseVector<DataType, IndexType>& r, const DenseVector<DataType, IndexType>& x) const
      {
        first().apply_transposed(r, x);
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x, const VectorTypeR& y, DataType alpha = DataType(1)) const
      {
        first().apply_transposed(r.first(), x.first(), y.first(), alpha);
      }
 
      void apply_transposed(DenseVector<DataType, IndexType>& r, const DenseVector<DataType, IndexType>& x,
        const DenseVector<DataType, IndexType>& y, DataType alpha = DataType(1)) const
      {
        first().apply_transposed(r, x, y, alpha);
      }
 
      DataType max_rel_diff(const PowerDiagMatrix& x) const
      {
        return first().max_rel_diff(x.first());
      }
 
      VectorTypeL create_vector_l() const
      {
        return VectorTypeL(first().create_vector_l());
      }
 
      VectorTypeR create_vector_r() const
      {
        return VectorTypeR(first().create_vector_r());
      }
 
      void scale_rows(const PowerDiagMatrix& a, const VectorTypeL& w)
      {
        first().scale_rows(a.first(), w.first());
      }
 
      void scale_cols(const PowerDiagMatrix& a, const VectorTypeR& w)
      {
        first().scale_cols(a.first(), w.first());
      }
 
      Index get_length_of_line(const Index row) const
      {
        return this->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
      {
        this->first().set_line(row, pval_set, pcol_set, col_start, stride);
      }
 
      void set_line_reverse(const Index row, const DataType * const pval_set, const Index stride = 1)
      {
        this->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);
      }
 
      template <typename SubType2_>
      void convert(const PowerDiagMatrix<SubType2_, 1> & other)
      {
        this->first().convert(other.first());
      }
 
      template <typename SubType2_>
      void convert_reverse(PowerDiagMatrix<SubType2_, 1> & other) const
      {
        this->first().convert_reverse(other.first());
      }
 
      bool same_layout(const PowerDiagMatrix& x) const
      {
        return (this->name() == x.name()) && (this->first().same_layout(x.first()));
      }
    };
  } // namespace LAFEM
} // namespace FEAT