null_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/base.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
 
 
namespace FEAT
{
  namespace LAFEM
  {
    namespace Intern
    {
      template<typename DT_, typename IT_, int size_>
      struct BlockedVectorHelper
      {
        static_assert(size_ > 1, "invalid block size");
        typedef DenseVectorBlocked<DT_, IT_, size_> VectorType;
      };
 
      template<typename DT_, typename IT_>
      struct BlockedVectorHelper<DT_, IT_, 1>
      {
        typedef DenseVector<DT_, IT_> VectorType;
      };
    } // namespace Intern
 
    template<typename DT_, typename IT_, int BlockHeight_ = 1, int BlockWidth_ = 1>
    class NullMatrix
    {
      static_assert(BlockHeight_ > 0, "invalid block size");
      static_assert(BlockWidth_ > 0, "invalid block size");
 
    public:
      typedef DT_ DataType;
      typedef IT_ IndexType;
      static constexpr int BlockHeight = BlockHeight_;
      static constexpr int BlockWidth = BlockWidth_;
      typedef Tiny::Matrix<DataType, BlockHeight, BlockWidth> ValueType;
 
      typedef const IT_* ImageIterator;
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      using ContainerType = NullMatrix<DT2_, IT2_, BlockHeight_, BlockWidth_>;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByDI = ContainerType<DataType2_, IndexType2_>;
 
      typedef typename Intern::BlockedVectorHelper<DT_, IT_, BlockHeight_>::VectorType VectorTypeL;
      typedef typename Intern::BlockedVectorHelper<DT_, IT_, BlockWidth_>::VectorType VectorTypeR;
 
      static constexpr bool is_global = false;
      static constexpr bool is_local = true;
 
    private:
      Index _num_rows, _num_cols;
 
    public:
 
      explicit NullMatrix() :
        _num_rows(0u), _num_cols(0u)
      {
      }
 
      explicit NullMatrix(Index rows_in, Index columns_in) :
        _num_rows(rows_in), _num_cols(columns_in)
      {
      }
 
      NullMatrix(NullMatrix && other) :
        _num_rows(other._num_rows), _num_cols(other._num_cols)
      {
      }
 
      NullMatrix & operator= (NullMatrix && other)
      {
        this->_num_rows = other._num_rows;
        this->_num_cols = other._num_cols;
 
        return *this;
      }
 
      NullMatrix clone(CloneMode DOXY(clone_mode) = CloneMode::Weak) const
      {
        return NullMatrix(_num_rows, _num_cols);
      }
 
      template<typename DT2_, typename IT2_>
      void clone(
        const NullMatrix<DT2_, IT2_, BlockHeight_, BlockWidth_> & other,
        CloneMode DOXY(clone_mode) = CloneMode::Weak)
      {
        this->_num_rows = other.rows();
        this->_num_cols = other.columns();
      }
 
      template <typename DT2_, typename IT2_>
      void convert(const NullMatrix<DT2_, IT2_, BlockHeight_, BlockWidth_> & other)
      {
        this->_num_rows = other.rows();
        this->_num_cols = other.columns();
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index rows() const
      {
        Index result(this->_num_rows);
        if (perspective_ == Perspective::pod)
        {
            result *= Index(BlockHeight_);
        }
        return result;
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index columns() const
      {
        Index result(this->_num_cols);
        if (perspective_ == Perspective::pod)
        {
            result *= Index(BlockWidth_);
        }
        return result;
      }
 
      template <Perspective perspective_ = Perspective::native>
      Index used_elements() const
      {
        return Index(0);
      }
 
      static String name()
      {
        return "NullMatrix";
      }
 
      void resize(Index rows_in, Index columns_in)
      {
        this->_num_rows = rows_in;
        this->_num_cols = columns_in;
      }
 
      std::size_t bytes() const
      {
        return std::size_t(0);
      }
 
      void copy(const NullMatrix & DOXY(x), bool DOXY(full) = false)
      {
        // nothing to do here
      }
 
      void format(const DT_ DOXY(value) = DT_(0))
      {
        // nothing to do here
      }
 
 
      void axpy(
                const NullMatrix & x,
                const NullMatrix & y,
                const DT_ DOXY(alpha) = DT_(1))
      {
        XASSERTM(x.rows() == y.rows(), "Matrix rows do not match!");
        XASSERTM(x.rows() == this->rows(), "Matrix rows do not match!");
        XASSERTM(x.columns() == y.columns(), "Matrix columns do not match!");
        XASSERTM(x.columns() == this->columns(), "Matrix columns do not match!");
        XASSERTM(x.used_elements() == y.used_elements(), "Matrix used_elements do not match!");
        XASSERTM(x.used_elements() == this->used_elements(), "Matrix used_elements do not match!");
 
        // nothing to do here
      }
 
      void scale(const NullMatrix & x, const DT_ DOXY(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!");
 
        // nothing to do here
      }
 
      DT_ norm_frobenius() const
      {
        return DT_(0);
      }
 
      void row_norm2(VectorTypeL& row_norms) const
      {
        XASSERTM(row_norms.size() == this->rows(), "Matrix/Vector dimension mismatch");
 
        row_norms.format();
      }
 
      void row_norm2sqr(VectorTypeL& row_norms) const
      {
        XASSERTM(row_norms.size() == this->rows(), "Matrix/Vector dimension mismatch");
 
        row_norms.format();
      }
 
      void row_norm2sqr(VectorTypeL& row_norms, const VectorTypeR& scal) const
      {
        XASSERTM(row_norms.size() == this->rows(), "Matrix/Vector dimension mismatch");
        XASSERTM(scal.size() == this->columns(), "Matrix/scalings dimension mismatch");
 
        row_norms.format();
      }
 
      NullMatrix<DT_, IT_, BlockWidth_, BlockHeight_> transpose() const
      {
        return NullMatrix<DT_, IT_, BlockWidth_, BlockHeight_>(_num_cols, _num_rows);
      }
 
      void transpose(const NullMatrix<DT_, IT_, BlockWidth_, BlockHeight_> & x)
      {
        x = this->transpose();
      }
 
      void apply(DenseVector<DT_, IT_> & r, const DenseVector<DT_, IT_> & x) const
      {
        XASSERTM(r.size() == this->rows<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns<Perspective::pod>(), "Vector size of x does not match!");
 
        r.format();
      }
 
      void apply_transposed(DenseVector<DT_, IT_> & r, const DenseVector<DT_, IT_> & x) const
      {
        XASSERTM(r.size() == this->columns<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows<Perspective::pod>(), "Vector size of x does not match!");
 
        r.format();
      }
 
      void apply(DenseVectorBlocked<DT_, IT_, BlockHeight_> & r, const DenseVector<DT_, IT_> & x) const
      {
        XASSERTM(r.size() == this->rows(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns<Perspective::pod>(), "Vector size of x does not match!");
 
        r.format();
      }
 
      void apply_transposed(DenseVectorBlocked<DT_, IT_, BlockWidth_> & r, const DenseVector<DT_, IT_> & x) const
      {
        XASSERTM(r.size() == this->columns(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows<Perspective::pod>(), "Vector size of x does not match!");
 
        r.format();
      }
 
      void apply(DenseVector<DT_, IT_> & r, const DenseVectorBlocked<DT_, IT_, BlockWidth_> & x) const
      {
        XASSERTM(r.size() == this->rows<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns(), "Vector size of x does not match!");
 
        r.format();
      }
 
      void apply_transposed(DenseVector<DT_, IT_> & r, const DenseVectorBlocked<DT_, IT_, BlockHeight_> & x) const
      {
        XASSERTM(r.size() == this->columns<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows(), "Vector size of x does not match!");
 
        r.format();
      }
 
      void apply(DenseVectorBlocked<DT_, IT_, BlockHeight_> & r, const DenseVectorBlocked<DT_, IT_, BlockWidth_> & 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!");
 
        r.format();
      }
 
      void apply_transposed(DenseVectorBlocked<DT_, IT_, BlockWidth_> & r, const DenseVectorBlocked<DT_, IT_, BlockHeight_> & 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!");
 
        r.format();
      }
 
      void apply(
                 DenseVector<DT_, IT_> & r,
                 const DenseVector<DT_, IT_> & x,
                 const DenseVector<DT_, IT_> & y,
                 const DT_ DOXY(alpha) = DT_(1)) const
      {
        XASSERTM(r.size() == this->rows<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns<Perspective::pod>(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->rows<Perspective::pod>(), "Vector size of y does not match!");
 
        r.copy(y);
      }
 
      void apply_transposed(
        DenseVector<DT_, IT_> & r,
        const DenseVector<DT_, IT_> & x,
        const DenseVector<DT_, IT_> & y,
        const DT_ DOXY(alpha) = DT_(1)) const
      {
        XASSERTM(r.size() == this->columns<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows<Perspective::pod>(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->columns<Perspective::pod>(), "Vector size of y does not match!");
 
        r.copy(y);
      }
 
      void apply(
                 DenseVectorBlocked<DT_, IT_, BlockHeight_> & r,
                 const DenseVector<DT_, IT_> & x,
                 const DenseVectorBlocked<DT_, IT_, BlockHeight_> & y,
                 const DT_ DOXY(alpha) = DT_(1)) const
      {
        XASSERTM(r.size() == this->rows(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns<Perspective::pod>(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->rows(), "Vector size of y does not match!");
 
        r.copy(y);
      }
      void apply_transposed(
        DenseVectorBlocked<DT_, IT_, BlockWidth_> & r,
        const DenseVector<DT_, IT_> & x,
        const DenseVectorBlocked<DT_, IT_, BlockWidth_> & y,
        const DT_ DOXY(alpha) = DT_(1)) const
      {
        XASSERTM(r.size() == this->columns(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows<Perspective::pod>(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->columns(), "Vector size of y does not match!");
 
        r.copy(y);
      }
 
 
      void apply(
                 DenseVector<DT_, IT_> & r,
                 const DenseVectorBlocked<DT_, IT_, BlockWidth_> & x,
                 const DenseVector<DT_, IT_> & y,
                 const DT_ DOXY(alpha) = DT_(1)) const
      {
        XASSERTM(r.size() == this->rows<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->columns(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->rows<Perspective::pod>(), "Vector size of y does not match!");
 
        r.copy(y);
      }
 
      void apply_transposed(
        DenseVector<DT_, IT_> & r,
        const DenseVectorBlocked<DT_, IT_, BlockHeight_> & x,
        const DenseVector<DT_, IT_> & y,
        const DT_ DOXY(alpha) = DT_(1)) const
      {
        XASSERTM(r.size() == this->columns<Perspective::pod>(), "Vector size of r does not match!");
        XASSERTM(x.size() == this->rows(), "Vector size of x does not match!");
        XASSERTM(y.size() == this->columns<Perspective::pod>(), "Vector size of y does not match!");
 
        r.copy(y);
      }
 
      void apply(
                 DenseVectorBlocked<DT_, IT_, BlockHeight_> & r,
                 const DenseVectorBlocked<DT_, IT_, BlockWidth_> & x,
                 const DenseVectorBlocked<DT_, IT_, BlockHeight_> & y,
                 const DT_ DOXY(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!");
 
        r.copy(y);
      }
 
      void apply_transposed(
        DenseVectorBlocked<DT_, IT_, BlockWidth_> & r,
        const DenseVectorBlocked<DT_, IT_, BlockHeight_> & x,
        const DenseVectorBlocked<DT_, IT_, BlockWidth_> & y,
        const DT_ DOXY(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!");
 
        r.copy(y);
      }
 
      void apply(
                 DenseVectorBlocked<DT_, IT_, BlockHeight_> & r,
                 const DenseVectorBlocked<DT_, IT_, BlockWidth_> & x,
                 const DenseVector<DT_, IT_> & y,
                 const DT_ DOXY(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<Perspective::pod>(), "Vector size of y does not match!");
 
        r.copy(y);
      }
 
      void apply_transposed(
        DenseVectorBlocked<DT_, IT_, BlockWidth_> & r,
        const DenseVectorBlocked<DT_, IT_, BlockHeight_> & x,
        const DenseVector<DT_, IT_> & y,
        const DT_ DOXY(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<Perspective::pod>(), "Vector size of y does not match!");
 
        r.copy(y);
      }
 
 
      void lump_rows(VectorTypeL& lump) const
      {
        XASSERTM(lump.size() == rows(), "lump vector size does not match matrix row count!");
        lump.format();
      }
 
      VectorTypeL lump_rows() const
      {
        VectorTypeL lump = create_vector_l();
        lump_rows(lump);
        return lump;
      }
 
 
      void extract_diag(VectorTypeL & diag) const
      {
        XASSERTM(diag.size() == rows(), "diag size does not match matrix row count!");
        XASSERTM(rows() == columns(), "matrix is not square!");
        diag.format();
      }
 
      VectorTypeL extract_diag() const
      {
        VectorTypeL diag = create_vector_l();
        extract_diag(diag);
        return diag;
      }
 
      // Returns a new compatible L-Vector.
      VectorTypeL create_vector_l() const
      {
        return VectorTypeL(this->rows());
      }
 
      // Returns a new compatible R-Vector.
      VectorTypeR create_vector_r() const
      {
        return VectorTypeR(this->columns());
      }
 
      Index get_length_of_line(const Index) const
      {
        return Index(0);
      }
 
      void set_line(const Index, DT_ * const, IT_ * const, const Index, const Index = 1) const
      {
        // nothing to do here
      }
 
      void set_line_reverse(const Index, DT_ * const, const Index = 1)
      {
        // nothing to do here
      }
 
      uint64_t get_checkpoint_size()
      {
        return size_t(0);
      }
 
      void restore_from_checkpoint_data(std::vector<char> & /*data*/)
      {
        // nothing to do here
      }
 
      void set_checkpoint_data(std::vector<char>& /*data*/)
      {
        // nothing to do here
      }
 
      /* ******************************************************************* */
      /*  A D J A C T O R   I N T E R F A C E   I M P L E M E N T A T I O N  */
      /* ******************************************************************* */
    public:
      inline Index get_num_nodes_domain() const
      {
        return rows();
      }
 
      inline Index get_num_nodes_image() const
      {
        return columns();
      }
 
      inline ImageIterator image_begin(Index domain_node) const
      {
        XASSERTM(domain_node < rows(), "Domain node index out of range");
        return nullptr;
      }
 
      inline ImageIterator image_end(Index domain_node) const
      {
        XASSERTM(domain_node < rows(), "Domain node index out of range");
        return nullptr;
      }
    }; // class NullMatrix
  } // namespace LAFEM
} // namespace FEAT