power_full_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_row_matrix.hpp>
#include <kernel/lafem/power_col_matrix.hpp>
#include <kernel/lafem/power_vector.hpp>
#include <kernel/lafem/sparse_layout.hpp>
#include <kernel/lafem/container.hpp>
 
 
#include <fstream>
 
namespace FEAT
{
  namespace LAFEM
  {
    namespace Intern
    {
      // helper class for extract_diag function of PowerFullMatrix
      template<int n_>
      struct ExtDiagPFM;
    }
 
    template<
      typename SubType_,
      int width_,
      int height_>
    class PowerFullMatrix
    {
      static_assert((width_ > 0) && (height_ > 0), "invalid matrix dimensions");
 
    public:
      typedef PowerColMatrix<PowerRowMatrix<SubType_, width_>, height_> ContClass;
      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, height_> VectorTypeL;
      typedef PowerVector<typename SubMatrixType::VectorTypeR, width_> VectorTypeR;
      template <typename DT2_ = DataType, typename IT2_ = IndexType>
      using ContainerType = PowerFullMatrix<
        typename SubType_::template ContainerType<DT2_, IT2_>, width_, height_>;
 
      static constexpr int num_row_blocks = height_;
      static constexpr int num_col_blocks = width_;
 
    protected:
      // the container
      ContClass _container;
 
    protected:
      explicit PowerFullMatrix(ContClass&& cont) :
        _container(std::move(cont))
      {
      }
 
    public:
      PowerFullMatrix()
      {
      }
 
      explicit PowerFullMatrix(const SparseLayout<IndexType, layout_id>& layout) :
        _container(layout)
      {
      }
 
      PowerFullMatrix(PowerFullMatrix&& other) :
        _container(std::move(other._container))
      {
      }
 
      explicit PowerFullMatrix(FileMode mode, const String& filename)
      {
        ContClass other(mode, filename);
        _container = std::move(other);
      }
 
      explicit PowerFullMatrix(FileMode mode, std::istream& file, const String& directory = "")
      {
        ContClass other(mode, file, directory);
        _container = std::move(other);
      }
 
      void read_from(FileMode mode, const String& filename)
      {
        ContClass other(mode, filename);
        _container = std::move(other);
      }
 
      PowerFullMatrix& operator=(PowerFullMatrix&& other)
      {
        if(this != &other)
        {
          _container = std::move(other._container);
        }
        return *this;
      }
 
      PowerFullMatrix(const PowerFullMatrix&) = delete;
      PowerFullMatrix& operator=(const PowerFullMatrix&) = delete;
 
      virtual ~PowerFullMatrix()
      {
      }
 
      void write_out(FileMode mode, const String& filename) const
      {
        _container.write_out(mode, filename);
      }
 
      PowerFullMatrix clone(LAFEM::CloneMode mode = LAFEM::CloneMode::Weak) const
      {
        return PowerFullMatrix(_container.clone(mode));
      }
 
      std::size_t bytes() const
      {
        return _container.bytes();
      }
 
      template<int i_, int j_>
      SubMatrixType& at()
      {
        static_assert((0 <= i_) && (i_ < height_), "invalid sub-matrix index");
        static_assert((0 <= j_) && (j_ < width_), "invalid sub-matrix index");
        return _container.template at<i_, 0>().template at<0, j_>();
      }
 
      template<int i_, int j_>
      const SubMatrixType& at() const
      {
        static_assert((0 <= i_) && (i_ < height_), "invalid sub-matrix index");
        static_assert((0 <= j_) && (j_ < width_), "invalid sub-matrix index");
        return _container.template at<i_, 0>().template at<0, j_>();
      }
 
      SubMatrixType& get(int i, int j)
      {
        XASSERTM((0 <= i) && (i < height_), "invalid sub-matrix row index");
        XASSERTM((0 <= j) && (j < width_), "invalid sub-matrix column index");
        return _container.get(i, 0).get(0, j);
      }
 
      const SubMatrixType& get(int i, int j) const
      {
        XASSERTM((0 <= i) && (i < height_), "invalid sub-matrix row index");
        XASSERTM((0 <= j) && (j < width_), "invalid sub-matrix column index");
        return _container.get(i, 0).get(0, j);
      }
 
      int row_blocks() const
      {
        return num_row_blocks;
      }
 
      int col_blocks() const
      {
        return num_col_blocks;
      }
 
      Index get_length_of_line(const Index row) const
      {
        return _container.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
      {
        _container.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)
      {
        _container.set_line_reverse(row, pval_set, stride);
      }
 
      Index row_degree(const Index row) const
      {
        return _container.row_degree(row);
      }
 
      template<typename IT2_>
      Index get_row_col_indices(const Index row, IT2_* const pcol_idx, const IT2_ col_offset) const
      {
        return _container.get_row_col_indices(row, pcol_idx, col_offset);
      }
 
      template<typename DT2_>
      Index get_row_values(const Index row, DT2_ * const pvals) const
      {
        return _container.get_row_values(row, pvals);
      }
 
      template<typename DT2_>
      Index set_row_values(const Index row, const DT2_ * const pvals)
      {
        return _container.set_row_values(row, pvals);
      }
 
      ContClass& get_container()
      {
        return _container;
      }
 
      const ContClass& get_container() const
      {
        return _container;
      }
 
      VectorTypeL create_vector_l() const
      {
        return _container.create_vector_l();
      }
 
      VectorTypeR create_vector_r() const
      {
        return _container.create_vector_r();
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index rows() const
      {
        return _container.template rows<perspective_>();
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index columns() const
      {
        return _container.template columns<perspective_>();
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index used_elements() const
      {
        return _container.template used_elements<perspective_>();
      }
 
      static String name()
      {
        return String("PowerFullMatrix<") + SubMatrixType::name() + "," + stringify(width_) + "," + stringify(height_) + ">";
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index size() const
      {
        return rows<perspective_>() * columns<perspective_>();
      }
 
      void format(DataType value = DataType(0))
      {
        _container.format(value);
      }
 
      void extract_diag(VectorTypeL& diag) const
      {
        static_assert(width_ == height_, "cannot extract diagonal from rectangular matrix");
        Intern::ExtDiagPFM<width_>::extract_diag(*this, diag);
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x) const
      {
        _container.apply(r, x);
      }
 
      void apply(DenseVector<DataType, IndexType>& r, const DenseVector<DataType, IndexType>& x) const
      {
        _container.apply(r, x);
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x, const VectorTypeL& y, DataType alpha = DataType(1)) const
      {
        _container.apply(r, x, y, alpha);
      }
 
      void apply(DenseVector<DataType, IndexType>& r, const DenseVector<DataType, IndexType>& x,
                 const DenseVector<DataType, IndexType>& y, DataType alpha = DataType(1)) const
      {
        _container.apply(r, x, y, alpha);
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x) const
      {
        _container.apply_transposed(r, x);
      }
 
      void apply_transposed(DenseVector<DataType, IndexType>& r, const DenseVector<DataType, IndexType>& x) const
      {
        _container.apply_transposed(r, x);
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x, const VectorTypeR& y, DataType alpha = DataType(1)) const
      {
        _container.apply_transposed(r, x, y, alpha);
      }
 
      void apply_transposed(DenseVector<DataType, IndexType>& r, const DenseVector<DataType, IndexType>& x,
        const DenseVector<DataType, IndexType>& y, DataType alpha = DataType(1)) const
      {
        _container.apply_transposed(r, x, y, alpha);
      }
 
      DataType max_rel_diff(const PowerFullMatrix& x) const
      {
        return this->_container.max_rel_diff(x._container);
      }
 
      bool same_layout(const PowerFullMatrix& x) const
      {
        return this->_container.same_layout(x._container);
      }
 
      template <typename SubType2_>
      void convert(const PowerFullMatrix<SubType2_, width_, height_> & other)
      {
        _container.convert(other._container);
      }
 
      template <typename SubType2_>
      void convert_reverse(PowerFullMatrix<SubType2_, width_, height_> & other) const
      {
        _container.convert_reverse(other._container);
      }
 
      std::uint64_t get_checkpoint_size(SerialConfig& config)
      {
        return _container.get_checkpoint_size(config);
      }
 
      void restore_from_checkpoint_data(std::vector<char> & data)
      {
        _container.restore_from_checkpoint_data(data);
      }
 
      std::uint64_t set_checkpoint_data(std::vector<char>& data, SerialConfig& config)
      {
        return _container.set_checkpoint_data(data, config);
      }
    }; // class PowerFullMatrix
 
    namespace Intern
    {
      // JacobiHelper specialization for PowerFullMatrix
      template<int n_>
      struct ExtDiagPFM
      {
        static_assert(n_ > 1, "invalid block size");
        template<typename SMT_, int m_>
        static void extract_diag(
          const PowerFullMatrix<SMT_, m_, m_>& matrix,
          PowerVector<typename SMT_::VectorTypeL, m_>& diag)
        {
          ExtDiagPFM<n_-1>::extract_diag(matrix, diag);
          matrix.template at<n_-1,n_-1>().extract_diag(diag.template at<n_-1>());
        }
      };
 
      template<>
      struct ExtDiagPFM<1>
      {
        template<typename SMT_, int m_>
        static void extract_diag(
          const PowerFullMatrix<SMT_, m_, m_>& matrix,
          PowerVector<typename SMT_::VectorTypeL, m_>& diag)
        {
          matrix.template at<0,0>().extract_diag(diag.template at<0>());
        }
      };
    } // namespace Intern
  } // namespace LAFEM
} // namespace FEAT