adp_solver_base.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/solver/base.hpp>
#include <kernel/global/alg_dof_parti_system.hpp>
#include <kernel/global/matrix.hpp>
#include <kernel/global/filter.hpp>
#include <kernel/global/vector.hpp>
#include <kernel/global/pmdcdsc_matrix.hpp>
 
namespace FEAT
{
  namespace Solver
  {
    template<typename Matrix_, typename Filter_, typename SolverBase_ = Solver::SolverBase<typename Matrix_::VectorTypeL>>
#ifndef DOXYGEN
    class ADPSolverBase;
#else // DOXYGEN
    class ADPSolverBase :
      public SolverBase_
    {
    protected:
      explicit ADPSolverBase(const GlobalMatrixType& matrix, const GlobalFilterType& filter);
 
      const Dist::Comm* _get_comm() const;
 
      Index _get_adp_vector_size() const;
 
      Index _get_adp_matrix_num_rows() const;
 
      Index _get_adp_matrix_num_cols() const;
 
      Index _get_adp_matrix_num_nzes() const;
 
      Index _get_num_global_nonzeros() const;
 
      Index _get_num_global_dofs() const;
 
      Index _get_num_owned_dofs() const;
 
      Index _get_global_dof_offset() const;
 
      String _get_adp_block_information() const;
 
      template<typename RPT_, typename CIT_>
      void _upload_symbolic(RPT_* row_ptr, CIT_* col_idx);
 
      template<typename DTV_, typename RPT_, typename CIT_>
      void _upload_numeric(DTV_* val, const RPT_* row_ptr, const CIT_* col_idx);
 
      template<typename DTV_>
      void _upload_vector(DTV_* val, const LocalVectorType& vector);
 
      template<typename DTV_>
      void _download_vector(LocalVectorType& vector, const DTV_* val);
    }; // class ADPSolverBase
#endif // DOXYGEN
 
 
    template<typename LocalMatrix_, typename Mirror_, typename LocalFilter_, typename SolverBase_>
    class ADPSolverBase<Global::Matrix<LocalMatrix_, Mirror_, Mirror_>, Global::Filter<LocalFilter_, Mirror_>, SolverBase_> :
      public SolverBase_
    {
    public:
      typedef LocalMatrix_ LocalMatrixType;
      typedef typename LocalMatrixType::VectorTypeL LocalVectorType;
      typedef LocalFilter_ LocalFilterType;
      typedef Mirror_ MirrorType;
 
      typedef Global::Matrix<LocalMatrixType, MirrorType, MirrorType> GlobalMatrixType;
      typedef Global::Vector<LocalVectorType, MirrorType> GlobalVectorType;
      typedef Global::Filter<LocalFilterType, MirrorType> GlobalFilterType;
 
      typedef SolverBase_ BaseClass;
 
      typedef GlobalMatrixType MatrixType;
      typedef GlobalVectorType VectorType;
      typedef GlobalFilterType FilterType;
 
      typedef Global::AlgDofParti<LocalVectorType, MirrorType> AlgDofPartiType;
      typedef Global::AlgDofPartiSystem<LocalMatrixType, LocalFilterType, MirrorType> AlgDofPartiSystemType;
 
      typedef typename AlgDofPartiType::DataType DataType;
      typedef typename AlgDofPartiType::IndexType IndexType;
 
      typedef LAFEM::SparseMatrixCSR<DataType, IndexType> ADPMatrixType;
      typedef LAFEM::DenseVector<DataType, IndexType> ADPVectorType;
 
    protected:
      const GlobalMatrixType& _system_matrix;
      const GlobalFilterType& _system_filter;
      std::shared_ptr<AlgDofPartiSystemType> _adp_system;
 
      explicit ADPSolverBase(const GlobalMatrixType& matrix, const GlobalFilterType& filter) :
        BaseClass(),
        _system_matrix(matrix),
        _system_filter(filter),
        _adp_system()
      {
      }
 
    public:
      const Dist::Comm* _get_comm() const
      {
        return this->_system_matrix.get_comm();
      }
 
      virtual void init_symbolic() override
      {
        BaseClass::init_symbolic();
 
        // assemble the algebraic dof partitioning
        this->_adp_system = std::make_shared<AlgDofPartiSystemType>(this->_system_matrix, this->_system_filter);
        this->_adp_system->init_symbolic();
      }
 
      virtual void done_symbolic() override
      {
        _adp_system.reset();
 
        BaseClass::done_symbolic();
      }
 
    protected:
      Index _get_adp_vector_size() const
      {
        return this->_adp_system->get_adp_vector_size();
      }
 
      Index _get_adp_matrix_num_rows() const
      {
        return this->_adp_system->get_adp_matrix_rows();
      }
 
      Index _get_adp_matrix_num_cols() const
      {
        return this->_adp_system->get_adp_matrix_cols();
      }
 
      Index _get_adp_matrix_num_nzes() const
      {
        return this->_adp_system->get_adp_matrix_nzes();
      }
 
      Index _get_num_global_nonzeros() const
      {
        return this->_adp_system->get_num_global_nonzeros();
      }
 
      Index _get_num_global_dofs() const
      {
        return this->_adp_system->get_num_global_dofs();
      }
 
      Index _get_num_owned_dofs() const
      {
        return this->_adp_system->get_num_owned_dofs();
      }
 
      Index _get_global_dof_offset() const
      {
        return this->_adp_system->get_global_dof_offset();
      }
 
      String _get_adp_block_information() const
      {
        return this->_adp_system->get_alg_dof_parti()->get_block_information();
      }
 
      template<typename RPT_, typename CIT_>
      void _upload_symbolic(RPT_* row_ptr, CIT_* col_idx)
      {
        this->_adp_system->upload_matrix_symbolic(row_ptr, col_idx);
      }
 
      template<typename DTV_, typename RPT_, typename CIT_>
      void _upload_numeric(DTV_* val, const RPT_* row_ptr, const CIT_* col_idx)
      {
        this->_adp_system->upload_matrix_numeric(val, row_ptr, col_idx);
        this->_adp_system->upload_filter();
        this->_adp_system->filter_matrix(val, row_ptr, col_idx);
      }
 
      template<typename DTV_>
      void _upload_vector(DTV_* val, const GlobalVectorType& vector)
      {
        this->_adp_system->upload_vector(val, vector.local());
      }
 
      template<typename DTV_>
      void _upload_vector(DTV_* val, const LocalVectorType& vector)
      {
        this->_adp_system->upload_vector(val, vector);
      }
 
      template<typename DTV_>
      void _download_vector(GlobalVectorType& vector, const DTV_* val)
      {
        this->_adp_system->download_vector(vector.local(), val);
      }
 
      template<typename DTV_>
      void _download_vector(LocalVectorType& vector, const DTV_* val)
      {
        this->_adp_system->download_vector(vector, val);
      }
    }; // class ADPSolverBase<Global::Matrix<...>, ...>
 
 
    template<typename MatrixB_, typename MatrixD_, typename GlobalFilter_, typename SolverBase_>
    class ADPSolverBase<Global::PMDCDSCMatrix<MatrixB_, MatrixD_>, GlobalFilter_, SolverBase_> :
      public SolverBase_
    {
    public:
      typedef Global::PMDCDSCMatrix<MatrixB_, MatrixD_> GlobalMatrixType;
      typedef typename GlobalMatrixType::VectorTypeL GlobalVectorType;
      typedef GlobalFilter_ GlobalFilterType;
 
      typedef typename GlobalMatrixType::LocalMatrixTypeS LocalMatrixTypeS;
      typedef typename LocalMatrixTypeS::VectorTypeL LocalVectorType;
 
      typedef SolverBase_ BaseClass;
 
      typedef GlobalVectorType VectorType;
 
      typedef typename LocalMatrixTypeS::DataType DataType;
      typedef typename LocalMatrixTypeS::IndexType IndexType;
 
    protected:
      const GlobalMatrixType& _system_matrix;
      const GlobalFilterType& _system_filter;
 
    private:
      Index _global_dof_count;
      Index _global_dof_offset;
      Index _owned_dof_count;
      Index _owned_num_nzes;
      Index _global_num_nzes;
 
    protected:
      explicit ADPSolverBase(const GlobalMatrixType& matrix, const GlobalFilterType& filter) :
        BaseClass(),
        _system_matrix(matrix),
        _system_filter(filter),
        _global_dof_count(0),
        _global_dof_offset(0),
        _owned_dof_count(0),
        _owned_num_nzes(0),
        _global_num_nzes(0)
      {
      }
 
    public:
      const Dist::Comm* _get_comm() const
      {
        return this->_system_matrix.get_comm();
      }
 
      virtual void init_symbolic() override
      {
        BaseClass::init_symbolic();
 
        // compute the ADP dimensions
        _system_matrix.adp_compute_counts(_global_dof_offset, _global_dof_count, _owned_dof_count, _owned_num_nzes, _global_num_nzes);
      }
 
    protected:
      Index _get_adp_vector_size() const
      {
        return _get_num_owned_dofs();
      }
 
      Index _get_adp_matrix_num_rows() const
      {
        return _get_num_owned_dofs();
      }
 
      Index _get_adp_matrix_num_cols() const
      {
        return _get_num_global_dofs();
      }
 
      Index _get_adp_matrix_num_nzes() const
      {
        return _owned_num_nzes;
      }
 
      Index _get_num_global_nonzeros() const
      {
        return _global_num_nzes;
      }
 
      Index _get_num_global_dofs() const
      {
        return _global_dof_count;
      }
 
      Index _get_num_owned_dofs() const
      {
        return _owned_dof_count;
      }
 
      Index _get_global_dof_offset() const
      {
        return _global_dof_offset;
      }
 
      String _get_adp_block_information() const
      {
        XABORTM("Block information for ADPSolverBase<PMDCDSCMatrix> not available yet!");
        return "-N/A-";
      }
 
      template<typename RPT_, typename CIT_>
      void _upload_symbolic(RPT_* row_ptr, CIT_* col_idx)
      {
        _system_matrix.adp_upload_symbolic(row_ptr, col_idx, _global_dof_offset);
      }
 
      template<typename DTV_, typename RPT_, typename CIT_>
      void _upload_numeric(DTV_* val, const RPT_* row_ptr, const CIT_* col_idx)
      {
        _system_matrix.adp_upload_numeric(val, row_ptr, col_idx);
        this->_filter_mat(val, row_ptr, col_idx, this->_system_filter.local());
      }
 
      template<typename DTV_>
      void _upload_vector(DTV_* val, const LAFEM::DenseVector<DataType, IndexType>& vector)
      {
        const Index n = _get_adp_vector_size();
        const DataType* vx = vector.elements();
 
        XASSERT(vector.size() == n);
 
        FEAT_PRAGMA_OMP(parallel for)
        for(Index i = 0; i < n; ++i)
          val[i] = DTV_(vx[i]);
      }
 
      template<typename DTV_>
      void _download_vector(LAFEM::DenseVector<DataType, IndexType>& vector, const DTV_* val)
      {
        const Index n = _get_adp_vector_size();
        DataType* vx = vector.elements();
 
        XASSERT(vector.size() == n);
 
        FEAT_PRAGMA_OMP(parallel for)
          for(Index i = 0; i < n; ++i)
            vx[i] = DataType(val[i]);
      }
 
      template<typename DTV_, typename RPT_, typename CIT_>
      void _filter_mat(DTV_*, const RPT_*, const CIT_*,
        const LAFEM::NoneFilter<DataType, IndexType>&)
      {
        // nothing to do here
      }
 
      template<typename DTV_, typename RPT_, typename CIT_>
      void _filter_mat(DTV_* val, const RPT_* row_ptr, const CIT_* col_idx,
        const LAFEM::UnitFilter<DataType, IndexType>& filter)
      {
        const IndexType n = filter.used_elements();
        const IndexType* fil_idx = filter.get_indices();
        for(IndexType i = 0; i < n; ++i)
        {
          const Index row = fil_idx[i];
          for(RPT_ j = row_ptr[row]; j < row_ptr[row+1]; ++j)
            val[j] = DTV_(Index(col_idx[j]) == row ? 1 : 0);
        }
      }
 
      template<typename DTV_, typename RPT_, typename CIT_>
      void _filter_mat(DTV_*, const RPT_*, const CIT_*,
        const LAFEM::MeanFilter<DataType, IndexType>& filter)
      {
        XASSERTM(filter.empty(), "LAFEM::MeanFilter is not supported by ADPSolverBase yet!");
      }
 
      template<typename DTV_, typename RPT_, typename CIT_>
      void _filter_mat(DTV_*, const RPT_*, const CIT_*,
        const Global::MeanFilter<DataType, IndexType>& filter)
      {
        XASSERTM(filter.empty(), "LAFEM::MeanFilter is not supported by ADPSolverBase yet!");
      }
 
      template<typename DTV_, typename RPT_, typename CIT_, typename SubFilter_>
      void _filter_mat(DTV_* val, const RPT_* row_ptr, const CIT_* col_idx,
        const LAFEM::FilterSequence<SubFilter_>& filter)
      {
        for(const auto& sf : filter)
          this->_filter_mat(val, row_ptr, col_idx, sf.second);
      }
 
      template<typename DTV_, typename RPT_, typename CIT_, typename First_, typename... Rest_>
      void _filter_mat(DTV_* val, const RPT_* row_ptr, const CIT_* col_idx,
        const LAFEM::FilterChain<First_, Rest_...>& filter)
      {
        this->_filter_mat(val, row_ptr, col_idx, filter.first());
        this->_filter_mat(val, row_ptr, col_idx, filter.rest());
      }
 
      template<typename DTV_, typename RPT_, typename CIT_, typename First_>
      void _filter_mat(DTV_* val, const RPT_* row_ptr, const CIT_* col_idx,
        const LAFEM::FilterChain<First_>& filter)
      {
        this->_filter_mat(val, row_ptr, col_idx, filter.first());
      }
    }; // class ADPSolverBase<Global::PMDCDSCMatrix<...>, ...>
 
 
 
    template<typename Matrix_, typename Filter_, typename SolverBase_>
    class ADPSolverBase :
      public SolverBase_
    {
      static_assert(Matrix_::is_local, "invalid instantiation of ADPSolverBase for non-local matrix type!");
 
    public:
      typedef Matrix_ MatrixType;
      typedef typename MatrixType::VectorTypeL VectorType;
      typedef Filter_ FilterType;
 
      typedef SolverBase_ BaseClass;
 
      typedef typename MatrixType::DataType DataType;
      typedef typename MatrixType::IndexType IndexType;
 
    protected:
      Dist::Comm _comm_self;
      const MatrixType& _system_matrix;
      const FilterType& _system_filter;
 
    protected:
      explicit ADPSolverBase(const MatrixType& matrix, const FilterType& filter) :
        _comm_self(Dist::Comm::self()),
        _system_matrix(matrix),
        _system_filter(filter)
      {
      }
 
      const Dist::Comm* _get_comm() const
      {
        return &this->_comm_self;
      }
 
      Index _get_adp_vector_size() const
      {
        return _system_matrix.rows();
      }
 
      Index _get_adp_matrix_num_rows() const
      {
        return _system_matrix.rows();
      }
 
      Index _get_adp_matrix_num_cols() const
      {
        return _system_matrix.columns();
      }
 
      Index _get_adp_matrix_num_nzes() const
      {
        return _system_matrix.template used_elements<LAFEM::Perspective::pod>();
      }
 
      Index _get_num_global_nonzeros() const
      {
        return _system_matrix.template used_elements<LAFEM::Perspective::pod>();
      }
 
      Index _get_num_global_dofs() const
      {
        return _system_matrix.rows();
      }
 
      Index _get_num_owned_dofs() const
      {
        return _system_matrix.rows();
      }
 
      Index _get_global_dof_offset() const
      {
        return Index(0);
      }
 
      String _get_adp_block_information() const
      {
        XABORTM("Block information for ADPSolverBase<...> not available yet!");
        return "-N/A-";
      }
 
      template<typename RPT_, typename CIT_>
      void _upload_symbolic(RPT_* row_ptr, CIT_* col_idx)
      {
        const Index num_rows = _system_matrix.rows();
        row_ptr[0] = RPT_(0);
        for(Index i(0); i < num_rows; ++i)
        {
          row_ptr[i+1u] = row_ptr[i] + RPT_(_system_matrix.get_row_col_indices(i, &col_idx[row_ptr[i]], CIT_(0)));
        }
      }
 
      template<typename DTV_, typename RPT_, typename CIT_>
      void _upload_numeric(DTV_* val, const RPT_* row_ptr, const CIT_* col_idx)
      {
        const Index num_rows = _system_matrix.rows();
        for(Index i(0); i < num_rows; ++i)
        {
          _system_matrix.get_row_values(i, &val[row_ptr[i]]);
        }
 
        this->_filter_mat(val, row_ptr, col_idx, this->_system_filter);
      }
 
      template<typename DTV_>
      void _upload_vector(DTV_* val, const VectorType& vector)
      {
        vector.set_vec(val);
      }
 
      template<typename DTV_>
      void _download_vector(VectorType& vector, const DTV_* val)
      {
        vector.set_vec_inv(val);
      }
 
      template<typename DTV_, typename RPT_, typename CIT_>
      void _filter_mat(DTV_*, const RPT_*, const CIT_*,
        const LAFEM::NoneFilter<DataType, IndexType>&)
      {
        // nothing to do here
      }
 
      template<typename DTV_, typename RPT_, typename CIT_>
      void _filter_mat(DTV_* val, const RPT_* row_ptr, const CIT_* col_idx,
        const LAFEM::UnitFilter<DataType, IndexType>& filter)
      {
        const IndexType n = filter.used_elements();
        const IndexType* fil_idx = filter.get_indices();
        for(IndexType i = 0; i < n; ++i)
        {
          const Index row = fil_idx[i];
          for(RPT_ j = row_ptr[row]; j < row_ptr[row+1]; ++j)
            val[j] = DTV_(Index(col_idx[j]) == row ? 1 : 0);
        }
      }
 
      template<typename DTV_, typename RPT_, typename CIT_>
      void _filter_mat(DTV_*, const RPT_*, const CIT_*,
        const LAFEM::MeanFilter<DataType, IndexType>& filter)
      {
        XASSERTM(filter.empty(), "LAFEM::MeanFilter is not supported by ADPSolverBase yet!");
      }
 
      template<typename DTV_, typename RPT_, typename CIT_>
      void _filter_mat(DTV_*, const RPT_*, const CIT_*,
        const Global::MeanFilter<DataType, IndexType>& filter)
      {
        XASSERTM(filter.empty(), "LAFEM::MeanFilter is not supported by ADPSolverBase yet!");
      }
 
      template<typename DTV_, typename RPT_, typename CIT_, typename SubFilter_>
      void _filter_mat(DTV_* val, const RPT_* row_ptr, const CIT_* col_idx,
        const LAFEM::FilterSequence<SubFilter_>& filter)
      {
        for(const auto& sf : filter)
          this->_filter_mat(val, row_ptr, col_idx, sf.second);
      }
 
      template<typename DTV_, typename RPT_, typename CIT_, typename First_, typename... Rest_>
      void _filter_mat(DTV_* val, const RPT_* row_ptr, const CIT_* col_idx,
        const LAFEM::FilterChain<First_, Rest_...>& filter)
      {
        this->_filter_mat(val, row_ptr, col_idx, filter.first());
        this->_filter_mat(val, row_ptr, col_idx, filter.rest());
      }
 
      template<typename DTV_, typename RPT_, typename CIT_, typename First_>
      void _filter_mat(DTV_* val, const RPT_* row_ptr, const CIT_* col_idx,
        const LAFEM::FilterChain<First_>& filter)
      {
        this->_filter_mat(val, row_ptr, col_idx, filter.first());
      }
    }; // class ADPSolverBase<...>
  } // namespace Solver
} // namespace FEAT