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/container.hpp> // required for LAFEM::CloneMode
#include <kernel/global/gate.hpp>
#include <kernel/global/vector.hpp>
#include <kernel/global/synch_mat.hpp>
 
namespace FEAT
{
  namespace Global
  {
    template<typename LocalMatrix_, typename RowMirror_, typename ColMirror_>
    class Matrix
    {
    public:
      typedef LocalMatrix_ LocalMatrixType;
      typedef RowMirror_ RowMirrorType;
      typedef ColMirror_ ColMirrorType;
 
      typedef typename LocalMatrix_::DataType DataType;
      typedef typename LocalMatrix_::IndexType IndexType;
 
      typedef typename LocalMatrix_::VectorTypeL LocalVectorTypeL;
      typedef typename LocalMatrix_::VectorTypeR LocalVectorTypeR;
 
      typedef Vector<LocalVectorTypeL, RowMirror_> VectorTypeL;
      typedef Vector<LocalVectorTypeR, ColMirror_> VectorTypeR;
 
      typedef Gate<LocalVectorTypeL, RowMirror_> GateRowType;
      typedef Gate<LocalVectorTypeR, ColMirror_> GateColType;
 
      template <typename LocalMatrix2_, typename RowMirror2_ = RowMirror_, typename ColMirror2_ = ColMirror_>
      using ContainerType = Matrix<LocalMatrix2_, RowMirror2_, ColMirror2_>;
 
      template <typename DataType2_, typename IndexType2_>
      using ContainerTypeByMDI = Matrix<
        typename LocalMatrix_::template ContainerType<DataType2_, IndexType2_>,
        typename RowMirror_::template MirrorType<DataType2_, IndexType2_>,
        typename ColMirror_::template MirrorType<DataType2_, IndexType2_> >;
 
      static constexpr bool is_global = true;
      static constexpr bool is_local = false;
 
    protected:
      const GateRowType* _row_gate;
      const GateColType* _col_gate;
      LocalMatrix_ _matrix;
 
    public:
      Matrix() :
        _row_gate(nullptr),
        _col_gate(nullptr),
        _matrix()
      {
      }
 
      template<typename... Args_>
      explicit Matrix(const GateRowType* row_gate, const GateColType* col_gate, Args_&&... args) :
        _row_gate(row_gate),
        _col_gate(col_gate),
        _matrix(std::forward<Args_>(args)...)
      {
        if((_row_gate != nullptr) && (_col_gate != nullptr))
        {
          // the gates must be defined on the same communicator
          XASSERT(_row_gate->get_comm() == _col_gate->get_comm());
        }
      }
 
      LocalMatrix_& local()
      {
        return _matrix;
      }
 
      const LocalMatrix_& local() const
      {
        return _matrix;
      }
 
      template<typename OtherGlobalMatrix_>
      void convert(GateRowType* row_gate, GateColType* col_gate, const OtherGlobalMatrix_ & other)
      {
        this->_row_gate = row_gate;
        this->_col_gate = col_gate;
        this->_matrix.convert(other.local());
      }
 
      const GateRowType* get_row_gate() const
      {
        return _row_gate;
      }
 
      const GateColType* get_col_gate() const
      {
        return _col_gate;
      }
 
      const Dist::Comm* get_comm() const
      {
        return _row_gate != nullptr ? _row_gate->get_comm() : nullptr;
      }
 
      Matrix clone(LAFEM::CloneMode mode = LAFEM::CloneMode::Weak) const
      {
        return Matrix(_row_gate, _col_gate, _matrix.clone(mode));
      }
 
      VectorTypeL create_vector_l() const
      {
        return VectorTypeL(_row_gate, _matrix.create_vector_l());
      }
 
      VectorTypeR create_vector_r() const
      {
        return VectorTypeR(_col_gate, _matrix.create_vector_r());
      }
 
      template<LAFEM::Perspective perspective_ = LAFEM::Perspective::pod>
      Index columns() const
      {
        return _col_gate->template get_num_global_dofs<perspective_>();
      }
 
      template<LAFEM::Perspective perspective_ = LAFEM::Perspective::pod>
      Index rows() const
      {
        return _row_gate->template get_num_global_dofs<perspective_>();
      }
 
      template<LAFEM::Perspective perspective_ = LAFEM::Perspective::pod>
      Index used_elements() const
      {
        Index my_used_elements(_matrix.template used_elements<perspective_>());
        const Dist::Comm* comm = (_row_gate != nullptr ? _row_gate->get_comm() : (_col_gate != nullptr ? _col_gate->get_comm() : nullptr));
        if((comm != nullptr) && (comm->size() > 1))
          comm->allreduce(&my_used_elements, &my_used_elements, std::size_t(1), Dist::op_sum);
        return my_used_elements;
      }
 
      std::size_t bytes() const
      {
        size_t my_bytes(_matrix.bytes());
        const Dist::Comm* comm = (_row_gate != nullptr ? _row_gate->get_comm() : (_col_gate != nullptr ? _col_gate->get_comm() : nullptr));
        if((comm != nullptr) && (comm->size() > 1))
          comm->allreduce(&my_bytes, &my_bytes, std::size_t(1), Dist::op_sum);
        return my_bytes;
      }
 
      void extract_diag(VectorTypeL& diag, bool sync = true) const
      {
        _matrix.extract_diag(diag.local());
        if(sync)
        {
          diag.sync_0();
        }
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x) const
      {
        _matrix.apply(r.local(), x.local());
        r.sync_0();
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x) const
      {
        _matrix.apply_transposed(r.local(), x.local());
        r.sync_0();
      }
 
      auto apply_async(VectorTypeL& r, const VectorTypeR& x) const -> decltype(r.sync_0_async())
      {
        _matrix.apply(r.local(), x.local());
        return r.sync_0_async();
      }
 
      auto apply_transposed_async(VectorTypeR& r, const VectorTypeL& x) const -> decltype(r.sync_0_async())
      {
        _matrix.apply_transposed(r.local(), x.local());
        return r.sync_0_async();
      }
 
      void apply(VectorTypeL& r, const VectorTypeR& x, const VectorTypeL& y, const DataType alpha = DataType(1)) const
      {
        // copy y to r
        r.copy(y);
 
        // convert from type-1 to type-0
        r.from_1_to_0();
 
        // r <- r + alpha*A*x
        _matrix.apply(r.local(), x.local(), r.local(), alpha);
 
        // synchronize r
        r.sync_0();
      }
 
      void apply_transposed(VectorTypeR& r, const VectorTypeL& x, const VectorTypeR& y, const DataType alpha = DataType(1)) const
      {
        // copy y to r
        r.copy(y);
 
        // convert from type-1 to type-0
        r.from_1_to_0();
 
        // r <- r + alpha*A^T*x
        _matrix.apply_transposed(r.local(), x.local(), r.local(), alpha);
 
        // synchronize r
        r.sync_0();
      }
 
      auto apply_async(VectorTypeL& r, const VectorTypeR& x, const VectorTypeL& y, const DataType alpha = DataType(1)) const -> decltype(r.sync_0_async())
      {
        // copy y to r
        r.copy(y);
 
        // convert from type-1 to type-0
        r.from_1_to_0();
 
        // r <- r + alpha*A*x
        _matrix.apply(r.local(), x.local(), r.local(), alpha);
 
        // synchronize r
        return r.sync_0_async();
      }
 
      auto apply_transposed_async(VectorTypeR& r, const VectorTypeL& x, const VectorTypeR& y, const DataType alpha = DataType(1)) const -> decltype(r.sync_0_async())
      {
        // copy y to r
        r.copy(y);
 
        // convert from type-1 to type-0
        r.from_1_to_0();
 
        // r <- r + alpha*A^T*x
        _matrix.apply_transposed(r.local(), x.local(), r.local(), alpha);
 
        // synchronize r
        return r.sync_0_async();
      }
 
      void lump_rows(VectorTypeL& lump, bool sync = true) const
      {
        XASSERTM(lump.local().size() == _matrix.rows(), "lump vector size does not match matrix row count!");
 
        _matrix.lump_rows(lump.local());
 
        if(sync)
        {
          lump.sync_0();
        }
      }
 
      VectorTypeL lump_rows(bool sync = true) const
      {
        VectorTypeL lump = create_vector_l();
 
        lump_rows(lump, sync);
 
        return lump;
      }
 
      LocalMatrix_ convert_to_1() const
      {
        ASSERTM(_col_gate, "Column gate is not set!");
        ASSERTM(_row_gate, "Row gate is not set!");
        LocalMatrix_ locmat = _matrix.clone(LAFEM::CloneMode::Weak);
        SynchMatrix<LocalMatrix_, RowMirror_> synch(*_row_gate->_comm, _row_gate->_ranks, _row_gate->_mirrors, _col_gate->_mirrors);
        synch.init(locmat);
        synch.exec(locmat);
        return locmat;
      }
 
      void convert_to_1(LocalMatrix_& matrix_sync_1, bool full_copy = true) const
      {
        matrix_sync_1.copy(this->local(), full_copy);
        if((_row_gate != nullptr) && (_col_gate != nullptr))
        {
          SynchMatrix<LocalMatrix_, RowMirror_> synch(*_row_gate->_comm, _row_gate->_ranks, _row_gate->_mirrors, _col_gate->_mirrors);
          synch.init(matrix_sync_1);
          synch.exec(matrix_sync_1);
        }
      }
 
      std::uint64_t get_checkpoint_size(LAFEM::SerialConfig& config)
      {
        return _matrix.get_checkpoint_size(config);
      }
 
      void restore_from_checkpoint_data(std::vector<char>& data)
      {
        _matrix.restore_from_checkpoint_data(data);
      }
 
      std::uint64_t set_checkpoint_data(std::vector<char>& data, LAFEM::SerialConfig& config)
      {
        return _matrix.set_checkpoint_data(data, config);
      }
    }; // class Matrix<...>
  } // namespace Global
} // namespace FEAT