ssor_precond.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/base_header.hpp>
#include <kernel/solver/base.hpp>
 
namespace FEAT
{
  namespace Solver
  {
    namespace Intern
    {
      int cuda_ssor_forward_apply(int m, double * y, const double * x, const double * csrVal, const int * csrColInd, int ncolors, double omega,
          int * colored_row_ptr, int * rows_per_color, int * inverse_row_ptr);
      int cuda_ssor_backward_apply(int m, double * y, const double * x, const double * csrVal, const int * csrColInd, int ncolors, double omega,
          int * colored_row_ptr, int * rows_per_color, int * inverse_row_ptr);
      template <int BlockSize_>
      int cuda_ssor_forward_bcsr_apply(int m, double * y, const double * x, const double * csrVal, const int * csrColInd, int ncolors, double omega,
          int * colored_row_ptr, int * rows_per_color, int * inverse_row_ptr);
      template <int BlockSize_>
      int cuda_ssor_backward_bcsr_apply(int m, double * y, const double * x, const double * csrVal, const int * csrColInd, int ncolors, double omega,
          int * colored_row_ptr, int * rows_per_color, int * inverse_row_ptr);
      void cuda_sor_done_symbolic(int * colored_row_ptr, int * rows_per_color, int * inverse_row_ptr);
      void cuda_sor_init_symbolic(int m, int nnz, const double * csrVal, const int * csrRowPtr, const int * csrColInd, int & ncolors,
        int* & colored_row_ptr, int* & rows_per_color, int* & inverse_row_ptr);
    }
 
    template<typename Matrix_>
    class SSORPrecondBase
    {
    public:
      typedef typename Matrix_::DataType DataType;
      typedef typename Matrix_::VectorTypeL VectorType;
 
      virtual void set_omega(DataType omega) = 0;
 
      virtual void init_symbolic() = 0;
 
      virtual void done_symbolic() = 0;
 
      virtual String name() const = 0;
 
      virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) = 0;
 
      virtual ~SSORPrecondBase() {}
    }; // SSORPrecondBase class
 
    template<PreferredBackend backend_, typename Matrix_, typename Filter_>
    class SSORPrecondWithBackend;
 
    template<typename Filter_, typename DT_, typename IT_>
    class SSORPrecondWithBackend<PreferredBackend::generic, LAFEM::SparseMatrixCSR<DT_, IT_>, Filter_> :
      public SSORPrecondBase<typename LAFEM::SparseMatrixCSR<DT_, IT_>>
    {
    public:
      typedef LAFEM::SparseMatrixCSR<DT_, IT_> MatrixType;
      typedef DT_ DataType;
      typedef IT_ IndexType;
      typedef Filter_ FilterType;
      typedef typename MatrixType::VectorTypeL VectorType;
 
    protected:
      const MatrixType& _matrix;
      const FilterType& _filter;
      DataType _omega;
 
    public:
      explicit SSORPrecondWithBackend(const MatrixType& matrix, const FilterType& filter, const DataType omega = DataType(1)) :
        _matrix(matrix),
        _filter(filter),
        _omega(omega)
      {
        if (_matrix.columns() != _matrix.rows())
        {
          XABORTM("Matrix is not square!");
        }
      }
 
      explicit SSORPrecondWithBackend(const String& section_name, const PropertyMap* section,
        const MatrixType& matrix, const FilterType& filter) :
        _matrix(matrix),
        _filter(filter),
        _omega(1)
      {
        if (_matrix.columns() != _matrix.rows())
        {
          XABORTM("Matrix is not square!");
        }
 
        auto omega_p = section->query("omega");
        if(omega_p.second && !omega_p.first.parse(this->_omega))
          throw ParseError(section_name + ".omega", omega_p.first, "a positive float");
      }
 
      virtual String name() const override
      {
        return "SSOR";
      }
 
      virtual void set_omega(DataType omega) override
      {
        XASSERT(omega > DataType(0));
        _omega = omega;
      }
 
       virtual void init_symbolic() override
      {
      }
 
      virtual void done_symbolic() override
      {
      }
 
      virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
      {
        XASSERTM(_matrix.rows() == vec_cor.size(), "matrix / vector size mismatch!");
        XASSERTM(_matrix.rows() == vec_def.size(), "matrix / vector size mismatch!");
 
        TimeStamp ts_start;
 
        // copy in-vector to out-vector
        vec_cor.copy(vec_def);
 
        _apply_intern(_matrix, vec_cor, vec_def);
 
        vec_cor.scale(vec_cor, _omega * (DataType(2.0) - _omega));
 
        this->_filter.filter_cor(vec_cor);
 
        TimeStamp ts_stop;
        Statistics::add_time_precon(ts_stop.elapsed(ts_start));
        Statistics::add_flops(_matrix.used_elements() *2 + 6 * vec_cor.size());
 
        return Status::success;
      }
 
    protected:
      void _apply_intern(const LAFEM::SparseMatrixCSR<DataType, IndexType>& matrix, VectorType& vec_cor, const VectorType& vec_def)
      {
        // create pointers
        DataType * pout(vec_cor.elements());
        const DataType * pin(vec_def.elements());
        const DataType * pval(matrix.val());
        const IndexType * pcol_ind(matrix.col_ind());
        const IndexType * prow_ptr(matrix.row_ptr());
        const IndexType n((IndexType(matrix.rows())));
 
        // __forward-insertion__
        // iteration over all rows
        for (Index i(0); i < n; ++i)
        {
          IndexType col;
          DataType d(0);
          // iteration over all elements on the left side of the main-diagonal
          for (col = prow_ptr[i]; pcol_ind[col] < i; ++col)
          {
            d += pval[col] * pout[pcol_ind[col]];
          }
          pout[i] = (pin[i] - _omega * d) / pval[col];
        }
 
        // __backward-insertion__
        // iteration over all rows
        for (Index i(n); i > 0;)
        {
          --i;
          IndexType col;
          DataType d(0);
          // iteration over all elements on the right side of the main-diagonal
          for (col = prow_ptr[i+1] - IndexType(1); pcol_ind[col] > i; --col)
          {
            d += pval[col] * pout[pcol_ind[col]];
          }
          pout[i] -= _omega * d / pval[col];
        }
      }
    }; // class SSORPrecondWithBackend<generic, SparseMatrixCSR>
 
    template<typename Filter_, typename DT_, typename IT_, int BlockHeight_, int BlockWidth_>
    class SSORPrecondWithBackend<PreferredBackend::generic, LAFEM::SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>, Filter_> :
      public SSORPrecondBase<typename LAFEM::SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_>>
    {
      static_assert(BlockHeight_ == BlockWidth_, "only square blocks are supported!");
    public:
      typedef LAFEM::SparseMatrixBCSR<DT_, IT_, BlockHeight_, BlockWidth_> MatrixType;
      typedef Filter_ FilterType;
      typedef typename MatrixType::VectorTypeL VectorType;
      typedef typename MatrixType::DataType DataType;
      typedef typename MatrixType::IndexType IndexType;
 
    protected:
      const MatrixType& _matrix;
      const FilterType& _filter;
      DataType _omega;
 
      void _apply_intern(const MatrixType & matrix, VectorType& vec_cor, const VectorType& vec_def)
      {
        // create pointers
        auto* pout(vec_cor.elements());
        const auto* pin(vec_def.elements());
        const auto* pval(matrix.val());
        const IndexType * pcol_ind(matrix.col_ind());
        const IndexType * prow_ptr(matrix.row_ptr());
        const IndexType n((IndexType(matrix.rows())));
        typename MatrixType::ValueType inverse;
 
        // __forward-insertion__
        // iteration over all rows
        for (Index i(0); i < n; ++i)
        {
          IndexType col;
          typename VectorType::ValueType d(0);
          // iteration over all elements on the left side of the main-diagonal
          for (col = prow_ptr[i]; pcol_ind[col] < i; ++col)
          {
            d += pval[col] * pout[pcol_ind[col]];
          }
          //pout[i] = (pin[i] - _omega * d) / pval[col];
          inverse.set_inverse(pval[col]);
          pout[i] = inverse * (pin[i] - _omega * d);
        }
 
        // __backward-insertion__
        // iteration over all rows
        for (Index i(n); i > 0;)
        {
          --i;
          IndexType col;
          typename VectorType::ValueType d(0);
          // iteration over all elements on the right side of the main-diagonal
          for (col = prow_ptr[i+1] - IndexType(1); pcol_ind[col] > i; --col)
          {
            d += pval[col] * pout[pcol_ind[col]];
          }
          //pout[i] -= _omega * d / pval[col];
          inverse.set_inverse(pval[col]);
          pout[i] = pout[i] - (_omega * inverse * d);
        }
      }
 
    public:
      explicit SSORPrecondWithBackend(const MatrixType& matrix, const FilterType& filter, const DataType omega = DataType(1)) :
        _matrix(matrix),
        _filter(filter),
        _omega(omega)
      {
        if (_matrix.columns() != _matrix.rows())
        {
          XABORTM("Matrix is not square!");
        }
      }
 
      explicit SSORPrecondWithBackend(const String& /*section_name*/, const PropertyMap* section,
        const MatrixType& matrix, const FilterType& filter) :
        _matrix(matrix),
        _filter(filter),
        _omega(1)
      {
        if (_matrix.columns() != _matrix.rows())
        {
          XABORTM("Matrix is not square!");
        }
 
        // Check if we have set _omega
        auto omega_p = section->query("omega");
        if(omega_p.second)
        {
          set_omega(DataType(std::stod(omega_p.first)));
        }
      }
 
      virtual String name() const override
      {
        return "SSOR";
      }
 
      virtual void set_omega(DataType omega) override
      {
        XASSERT(omega > DataType(0));
        _omega = omega;
      }
 
       virtual void init_symbolic() override
      {
      }
 
      virtual void done_symbolic() override
      {
      }
 
      virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
      {
        XASSERTM(_matrix.rows() == vec_cor.size(), "matrix / vector size mismatch!");
        XASSERTM(_matrix.rows() == vec_def.size(), "matrix / vector size mismatch!");
 
        TimeStamp ts_start;
 
        // copy in-vector to out-vector
        vec_cor.copy(vec_def);
 
        _apply_intern(_matrix, vec_cor, vec_def);
 
        this->_filter.filter_cor(vec_cor);
 
        TimeStamp ts_stop;
        Statistics::add_time_precon(ts_stop.elapsed(ts_start));
        Statistics::add_flops(_matrix.template used_elements<LAFEM::Perspective::pod>() * 2 + 6 * vec_cor.template size<LAFEM::Perspective::pod>());
 
        return Status::success;
      }
    }; // class SSORPrecondWithBackend<generic, SparseMatrixBCSR>
 
#ifdef FEAT_HAVE_CUDA
    template<typename Filter_>
    class SSORPrecondWithBackend<PreferredBackend::cuda, LAFEM::SparseMatrixCSR<double, unsigned int>, Filter_> :
      public SSORPrecondBase<LAFEM::SparseMatrixCSR<double, unsigned int>>
    {
    public:
      typedef LAFEM::SparseMatrixCSR<double, unsigned int> MatrixType;
      typedef Filter_ FilterType;
      typedef typename MatrixType::VectorTypeL VectorType;
      typedef typename MatrixType::DataType DataType;
 
    protected:
      const MatrixType& _matrix;
      const FilterType& _filter;
      double _omega;
      // row ptr permutation, sorted by color(each color sorted by amount of rows), start/end index per row
      int * _colored_row_ptr;
      // amount of rows per color (sorted by amount of rows)
      int * _rows_per_color;
      // mapping of idx to native row number
      int * _inverse_row_ptr;
      // number of colors
      int _ncolors;
 
    public:
      explicit SSORPrecondWithBackend(const MatrixType& matrix, const FilterType& filter, const DataType omega = DataType(1)) :
        _matrix(matrix),
        _filter(filter),
        _omega(omega),
        _colored_row_ptr(nullptr),
        _rows_per_color(nullptr),
        _inverse_row_ptr(nullptr),
        _ncolors(0)
      {
        if (_matrix.columns() != _matrix.rows())
        {
          XABORTM("Matrix is not square!");
        }
      }
 
      explicit SSORPrecondWithBackend(const String& /*section_name*/, const PropertyMap* section,
      const MatrixType& matrix, const FilterType& filter) :
        _matrix(matrix),
        _filter(filter),
        _omega(1),
        _colored_row_ptr(nullptr),
        _rows_per_color(nullptr),
        _inverse_row_ptr(nullptr),
        _ncolors(0)
      {
        if (_matrix.columns() != _matrix.rows())
        {
          XABORTM("Matrix is not square!");
        }
 
        // Check if we have set _omega
        auto omega_p = section->query("omega");
        if(omega_p.second)
        {
          set_omega(DataType(std::stod(omega_p.first)));
        }
      }
 
      virtual String name() const override
      {
        return "SSOR";
      }
 
      virtual void set_omega(DataType omega) override
      {
        XASSERT(omega > DataType(0));
        _omega = omega;
      }
 
      virtual void init_symbolic() override
      {
        Intern::cuda_sor_init_symbolic((int)_matrix.rows(), (int)_matrix.used_elements(), _matrix.val(), (const int*)_matrix.row_ptr(), (const int*)_matrix.col_ind(), _ncolors,
          _colored_row_ptr, _rows_per_color, _inverse_row_ptr);
      }
 
      virtual void done_symbolic() override
      {
        Intern::cuda_sor_done_symbolic(_colored_row_ptr, _rows_per_color, _inverse_row_ptr);
      }
 
      virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
      {
        XASSERTM(_matrix.rows() == vec_cor.size(), "matrix / vector size mismatch!");
        XASSERTM(_matrix.rows() == vec_def.size(), "matrix / vector size mismatch!");
 
        TimeStamp ts_start;
 
        int status = Intern::cuda_ssor_forward_apply((int)vec_cor.size(), vec_cor.elements(), vec_def.elements(), (const double*)_matrix.val(), (const int*)_matrix.col_ind(), _ncolors, _omega, _colored_row_ptr, _rows_per_color, _inverse_row_ptr);
        status |= Intern::cuda_ssor_backward_apply((int)vec_cor.size(), vec_cor.elements(), vec_def.elements(), (const double*)_matrix.val(), (const int*)_matrix.col_ind(), _ncolors, _omega, _colored_row_ptr, _rows_per_color, _inverse_row_ptr);
 
        vec_cor.scale(vec_cor, _omega * (2.0 - _omega));
 
        this->_filter.filter_cor(vec_cor);
 
        TimeStamp ts_stop;
        Statistics::add_time_precon(ts_stop.elapsed(ts_start));
        Statistics::add_flops(_matrix.used_elements() *2 + 6 * vec_cor.size());
 
        return (status == 0) ? Status::success :  Status::aborted;
      }
    }; // class SSORPrecondWithBackend<cuda, SparseMatrixCSR>
 
    template<typename Filter_, int BlockHeight_, int BlockWidth_>
    class SSORPrecondWithBackend<PreferredBackend::cuda, LAFEM::SparseMatrixBCSR<double, unsigned int, BlockHeight_, BlockWidth_>, Filter_> :
      public SSORPrecondBase<typename LAFEM::SparseMatrixBCSR<double, unsigned int, BlockHeight_, BlockWidth_>>
    {
      static_assert(BlockHeight_ == BlockWidth_, "only square blocks are supported!");
    public:
      typedef LAFEM::SparseMatrixBCSR<double, unsigned int, BlockHeight_, BlockWidth_> MatrixType;
      typedef Filter_ FilterType;
      typedef typename MatrixType::VectorTypeL VectorType;
      typedef typename MatrixType::DataType DataType;
 
    protected:
      const MatrixType& _matrix;
      const FilterType& _filter;
      double _omega;
      // row ptr permutation, sorted by color(each color sorted by amount of rows), start/end index per row
      int * _colored_row_ptr;
      // amount of rows per color (sorted by amount of rows)
      int * _rows_per_color;
      // mapping of idx to native row number
      int * _inverse_row_ptr;
      // number of colors
      int _ncolors;
 
    public:
      explicit SSORPrecondWithBackend(const MatrixType& matrix, const FilterType& filter, const DataType omega = DataType(1)) :
        _matrix(matrix),
        _filter(filter),
        _omega(omega),
        _colored_row_ptr(nullptr),
        _rows_per_color(nullptr),
        _inverse_row_ptr(nullptr),
        _ncolors(0)
      {
        if (_matrix.columns() != _matrix.rows())
        {
          XABORTM("Matrix is not square!");
        }
      }
 
      explicit SSORPrecondWithBackend(const String& /*section_name*/, const PropertyMap* section,
      const MatrixType& matrix, const FilterType& filter) :
        _matrix(matrix),
        _filter(filter),
        _omega(1),
        _colored_row_ptr(nullptr),
        _rows_per_color(nullptr),
        _inverse_row_ptr(nullptr),
        _ncolors(0)
      {
        if (_matrix.columns() != _matrix.rows())
        {
          XABORTM("Matrix is not square!");
        }
 
        // Check if we have set _omega
        auto omega_p = section->query("omega");
        if(omega_p.second)
        {
          set_omega(DataType(std::stod(omega_p.first)));
        }
      }
 
      virtual String name() const override
      {
        return "SSOR";
      }
 
      virtual void set_omega(DataType omega) override
      {
        XASSERT(omega > DataType(0));
        _omega = omega;
      }
 
      virtual void init_symbolic() override
      {
        Intern::cuda_sor_init_symbolic((int)_matrix.rows(), (int)_matrix.used_elements(), _matrix.template val<LAFEM::Perspective::pod>(), (const int*)_matrix.row_ptr(), (const int*)_matrix.col_ind(), _ncolors,
          _colored_row_ptr, _rows_per_color, _inverse_row_ptr);
      }
 
      virtual void done_symbolic() override
      {
        Intern::cuda_sor_done_symbolic(_colored_row_ptr, _rows_per_color, _inverse_row_ptr);
      }
 
      virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
      {
        XASSERTM(_matrix.rows() == vec_cor.size(), "matrix / vector size mismatch!");
        XASSERTM(_matrix.rows() == vec_def.size(), "matrix / vector size mismatch!");
 
        TimeStamp ts_start;
 
        int status = Intern::cuda_ssor_forward_bcsr_apply<BlockHeight_>((int)vec_cor.size(), vec_cor. template elements<LAFEM::Perspective::pod>(), vec_def.template elements<LAFEM::Perspective::pod>(), _matrix.template val<LAFEM::Perspective::pod>(), (const int*)_matrix.col_ind(), _ncolors, _omega, _colored_row_ptr, _rows_per_color, _inverse_row_ptr);
        status |= Intern::cuda_ssor_backward_bcsr_apply<BlockHeight_>((int)vec_cor.size(), vec_cor.template elements<LAFEM::Perspective::pod>(), vec_def.template elements<LAFEM::Perspective::pod>(), _matrix.template val<LAFEM::Perspective::pod>(), (const int*)_matrix.col_ind(), _ncolors, _omega, _colored_row_ptr, _rows_per_color, _inverse_row_ptr);
 
        vec_cor.scale(vec_cor, _omega * (2.0 - _omega));
 
        this->_filter.filter_cor(vec_cor);
 
        TimeStamp ts_stop;
        Statistics::add_time_precon(ts_stop.elapsed(ts_start));
        Statistics::add_flops(_matrix.template used_elements<LAFEM::Perspective::pod>() *2 + 6 * vec_cor.template size<LAFEM::Perspective::pod>());
 
        return (status == 0) ? Status::success :  Status::aborted;
      }
    }; // class SSORPrecondWithBackend<cuda, SparseMatrixBCSR>
#endif //FEAT_HAVE_CUDA
 
    template<PreferredBackend backend_, typename Matrix_, typename Filter_>
    class SSORPrecondWithBackend :
      public SSORPrecondBase<Matrix_>
    {
    public:
      template<typename DT_>
      explicit SSORPrecondWithBackend(const Matrix_&, const Filter_&, const DT_)
      {
      }
 
      explicit SSORPrecondWithBackend(const Matrix_&, const Filter_&)
      {
      }
 
      explicit SSORPrecondWithBackend(const String&, const PropertyMap*, const Matrix_&, const Filter_&)
      {
      }
 
      virtual Status apply(typename Matrix_::VectorTypeL &, const typename Matrix_::VectorTypeL &) override
      {
          XABORTM("not implemented yet!");
      }
 
      virtual String name() const override
      {
          XABORTM("not implemented yet!");
      }
 
      virtual void init_symbolic() override
      {
        XABORTM("not implemented yet!");
      }
 
      virtual void done_symbolic() override
      {
        XABORTM("not implemented yet!");
      }
 
      virtual void set_omega(typename Matrix_::DataType /*omega*/) override
      {
        XABORTM("not implemented yet!");
      }
    };
 
    template<typename Matrix_, typename Filter_>
    class SSORPrecond : public SolverBase<typename Matrix_::VectorTypeL>
    {
    private:
      std::shared_ptr<SSORPrecondBase<Matrix_>> _impl;
 
    public:
      typedef Matrix_ MatrixType;
      typedef Filter_ FilterType;
      typedef typename MatrixType::VectorTypeL VectorType;
      typedef typename MatrixType::DataType DataType;
 
      typedef SolverBase<VectorType> BaseClass;
 
    public:
      SSORPrecond(const PreferredBackend backend, const MatrixType& matrix, const FilterType& filter, DataType omega = DataType(1))
      {
        switch (backend)
        {
          case PreferredBackend::cuda:
            _impl = std::make_shared<SSORPrecondWithBackend<PreferredBackend::cuda, Matrix_, Filter_>>(matrix, filter, omega);
            break;
          case PreferredBackend::mkl:
          case PreferredBackend::generic:
          default:
            _impl = std::make_shared<SSORPrecondWithBackend<PreferredBackend::generic, Matrix_, Filter_>>(matrix, filter, omega);
        }
      }
 
      SSORPrecond(const PreferredBackend backend, const String& section_name, const PropertyMap* section, const MatrixType& matrix, const FilterType& filter, DataType omega = DataType(1)) :
        BaseClass(section_name, section)
      {
        switch (backend)
        {
          case PreferredBackend::cuda:
            _impl = std::make_shared<SSORPrecondWithBackend<PreferredBackend::cuda, Matrix_, Filter_>>(section_name, section, matrix, filter, omega);
            break;
          case PreferredBackend::mkl:
          case PreferredBackend::generic:
          default:
            _impl = std::make_shared<SSORPrecondWithBackend<PreferredBackend::generic, Matrix_, Filter_>>(section_name, section, matrix, filter, omega);
        }
      }
 
      virtual ~SSORPrecond()
      {
      }
 
      virtual String name() const override
      {
        return _impl->name();
      }
 
      virtual void set_omega(DataType omega)
      {
        _impl->set_omega(omega);
      }
 
      virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
      {
        return _impl->apply(vec_cor, vec_def);
      }
 
      virtual void init_symbolic() override
      {
        _impl->init_symbolic();
      }
 
      virtual void done_symbolic() override
      {
        _impl->done_symbolic();
      }
    }; // class SSORPrecond
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<SSORPrecond<Matrix_, Filter_>> new_ssor_precond(
      const PreferredBackend backend, const Matrix_& matrix, const Filter_& filter,
      const typename Matrix_::DataType omega = typename Matrix_::DataType(1))
    {
      return std::make_shared<SSORPrecond<Matrix_, Filter_>>
        (backend, matrix, filter, omega);
    }
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<SSORPrecond<Matrix_, Filter_>> new_ssor_precond(
      const String& section_name, const PropertyMap* section, PreferredBackend backend,
      const Matrix_& matrix, const Filter_& filter)
    {
      return std::make_shared<SSORPrecond<Matrix_, Filter_>>
        (backend, section_name, section, matrix, filter);
    }
  } // namespace Solver
} // namespace FEAT