sparse.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_PDELAB_BACKEND_SIMPLE_SPARSE_HH
#define DUNE_PDELAB_BACKEND_SIMPLE_SPARSE_HH

#include <vector>
#include <algorithm>
#include <functional>
#include <numeric>
#include <unordered_set>

#include <dune/common/typetraits.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/pdelab/backend/tags.hh>
#include <dune/pdelab/backend/backendselector.hh>
#include <dune/pdelab/backend/common/uncachedmatrixview.hh>
#include <dune/pdelab/backend/simple/descriptors.hh>

namespace Dune {
  namespace PDELab {
    namespace simple {

      template<typename _RowOrdering, typename _ColOrdering>
      class SparseMatrixPattern
        : public std::vector< std::unordered_set<std::size_t> >
      {

      public:

        typedef _RowOrdering RowOrdering;
        typedef _ColOrdering ColOrdering;

        typedef std::unordered_set<std::size_t> col_type;

        template<typename RI, typename CI>
        void add_link(const RI& ri, const CI& ci)
        {
          this->resize(_row_ordering.blockCount());
          (*this)[ri.back()].insert(ci.back());
        }

        SparseMatrixPattern(const RowOrdering& row_ordering, const ColOrdering& col_ordering)
          : _row_ordering(row_ordering)
          , _col_ordering(col_ordering)
        {}

      private:

        const RowOrdering& _row_ordering;
        const ColOrdering& _col_ordering;

      };

      template<template<typename> class C, typename ET, typename I>
      struct SparseMatrixData
      {
        typedef ET ElementType;
        typedef I  index_type;
        typedef std::size_t size_type;
        std::size_t _rows;
        std::size_t _cols;
        std::size_t _non_zeros;
        C<ElementType> _data;
        C<index_type>  _colindex;
        C<index_type>  _rowoffset;
      };

      template<typename GFSV, typename GFSU, template<typename> class C, typename ET, typename I>
      class SparseMatrixContainer
      {

      public:

        typedef SparseMatrixData<C,ET,I> Container;
        typedef ET ElementType;

        typedef ElementType field_type;
        typedef typename Container::size_type size_type;
        typedef I index_type;

        typedef GFSU TrialGridFunctionSpace;
        typedef GFSV TestGridFunctionSpace;

        typedef typename GFSV::Ordering::Traits::ContainerIndex RowIndex;
        typedef typename GFSU::Ordering::Traits::ContainerIndex ColIndex;

        template<typename RowCache, typename ColCache>
        using LocalView = UncachedMatrixView<SparseMatrixContainer,RowCache,ColCache>;

        template<typename RowCache, typename ColCache>
        using ConstLocalView = ConstUncachedMatrixView<const SparseMatrixContainer,RowCache,ColCache>;

        typedef OrderingBase<
          typename GFSV::Ordering::Traits::DOFIndex,
          typename GFSV::Ordering::Traits::ContainerIndex
          > RowOrdering;

        typedef OrderingBase<
          typename GFSU::Ordering::Traits::DOFIndex,
          typename GFSU::Ordering::Traits::ContainerIndex
          > ColOrdering;

        typedef SparseMatrixPattern<RowOrdering,ColOrdering> Pattern;

        template<typename GO>
        SparseMatrixContainer(const GO& go)
          : _container(make_shared<Container>())
        {
          allocate_matrix(_container, go, ElementType(0));
        }

        template<typename GO>
        SparseMatrixContainer(const GO& go, const ElementType& e)
          : _container(make_shared<Container>())
        {
          allocate_matrix(_container, go, e);
        }

        explicit SparseMatrixContainer(tags::unattached_container = tags::unattached_container())
        {}

        explicit SparseMatrixContainer(tags::attached_container)
        : _container(make_shared<Container>())
        {}

        SparseMatrixContainer(const SparseMatrixContainer& rhs)
          : _container(make_shared<Container>(*(rhs._container)))
        {}

        SparseMatrixContainer& operator=(const SparseMatrixContainer& rhs)
        {
          if (this == &rhs)
            return *this;
          if (attached())
          {
            (*_container) = (*(rhs._container));
          }
          else
          {
            _container = make_shared<Container>(*(rhs._container));
          }
          return *this;
        }

        void detach()
        {
          _container.reset();
        }

        void attach(shared_ptr<Container> container)
        {
          _container = container;
        }

        bool attached() const
        {
          return bool(_container);
        }

        const shared_ptr<Container>& storage() const
        {
          return _container;
        }

        size_type N() const
        {
          return _container->_rows;
        }

        size_type M() const
        {
          return _container->_cols;
        }

        SparseMatrixContainer& operator=(const ElementType& e)
        {
          std::fill(_container->_data.begin(),_container->_data.end(),e);
          return *this;
        }

        SparseMatrixContainer& operator*=(const ElementType& e)
        {
          using namespace std::placeholders;
          std::transform(_container->_data.begin(),_container->_data.end(),_container->_data.begin(),std::bind(std::multiplies<ET>(),e,_1));
          return *this;
        }

        template<typename V>
        void mv(const V& x, V& y) const
        {
          assert(y.N() == N());
          assert(x.N() == M());
          for (std::size_t r = 0; r < N(); ++r)
          {
            y.base()[r] = sparse_inner_product(r,x);
          }
        }

        template<typename V>
        void usmv(const ElementType alpha, const V& x, V& y) const
        {
          assert(y.N() == N());
          assert(x.N() == M());
          for (std::size_t r = 0; r < N(); ++r)
          {
            y.base()[r] += alpha * sparse_inner_product(r,x);
          }
        }

        ElementType& operator()(const RowIndex& ri, const ColIndex& ci)
        {
          // entries are in ascending order
          auto begin = _container->_colindex.begin() + _container->_rowoffset[ri[0]];
          auto end = _container->_colindex.begin() + _container->_rowoffset[ri[0]+1];
          auto it = std::lower_bound(begin,end,ci[0]);
          assert (it != end);
          return _container->_data[it - _container->_colindex.begin()];
        }

        const ElementType& operator()(const RowIndex& ri, const ColIndex& ci) const
        {
          // entries are in ascending order
          auto begin = _container->_colindex.begin() + _container->_rowoffset[ri[0]];
          auto end = _container->_colindex.begin() + _container->_rowoffset[ri[0]+1];
          auto it = std::lower_bound(begin,end,ci[0]);
          assert(it != end);
          return _container->_data[it - _container->_colindex.begin()];
        }

        const Container& base() const
        {
          return *_container;
        }

        Container& base()
        {
          return *_container;
        }

        void flush()
        {}

        void finalize()
        {}

        void clear_row(const RowIndex& ri, const ElementType& diagonal_entry)
        {
          std::fill(
            _container->_data.begin() + _container->_rowoffset[ri[0]],
            _container->_data.begin() + _container->_rowoffset[ri[0]+1], ElementType(0));
          (*this)(ri,ri) = diagonal_entry;
        }

      protected:
        template<typename GO>
        static void allocate_matrix(shared_ptr<Container> & c, const GO & go, const ElementType& e)
        {
          typedef typename Pattern::col_type col_type;
          Pattern pattern(go.testGridFunctionSpace().ordering(),go.trialGridFunctionSpace().ordering());
          go.fill_pattern(pattern);

          c->_rows = go.testGridFunctionSpace().size();
          c->_cols = go.trialGridFunctionSpace().size();
          // compute row offsets
          c->_rowoffset.resize(c->_rows+1);
          size_type offset = 0;
          auto calc_offset = [=](const col_type & entry) mutable -> size_t { offset += entry.size(); return offset; };
          std::transform(pattern.begin(), pattern.end(),
            c->_rowoffset.begin()+1,
            calc_offset);
          // compute non-zeros
          c->_non_zeros = c->_rowoffset.back();
          // allocate col/data vectors
          c->_data.resize(c->_non_zeros, e);
          c->_colindex.resize(c->_non_zeros);
          // copy pattern
          auto colit = c->_colindex.begin();
          c->_rowoffset[0] = 0;
          for (auto & row : pattern)
          {
            auto last = std::copy(row.begin(),row.end(),colit);
            std::sort(colit, last);
            colit = last;
          }
        }

        template<typename V>
        ElementType sparse_inner_product (std::size_t row, const V & x) const {
            std::size_t begin = _container->_rowoffset[row];
            std::size_t end = _container->_rowoffset[row+1];
            auto accu = [](const ElementType & a, const ElementType & b) -> ElementType { return a+b; };
            auto mult = [=](const ElementType & a, const index_type & i) -> ElementType { return a * x.base()[i]; };
            return std::inner_product(
              &_container->_data[begin],
              &_container->_data[end],
              &_container->_colindex[begin],
              ElementType(0),
              accu, mult
              );
        }

        shared_ptr< Container > _container;
      };

    } // namespace simple
  } // namespace PDELab
} // namespace Dune

#endif // DUNE_PDELAB_BACKEND_SIMPLE_SPARSE_HH