diagonallocalmatrix.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_PDELAB_GRIDOPERATOR_COMMON_DIAGONALLOCALMATRIX_HH
#define DUNE_PDELAB_GRIDOPERATOR_COMMON_DIAGONALLOCALMATRIX_HH
 
#include <dune/pdelab/gridoperator/common/localmatrix.hh>
 
namespace Dune {
  namespace PDELab {
 
        template<typename T, typename W = T>
     class DiagonalLocalMatrix
        {
        public:
 
 
          typedef std::vector<T> BaseContainer;
 
          typedef typename BaseContainer::value_type  value_type;
 
          typedef typename BaseContainer::size_type    size_type;
 
          typedef typename BaseContainer::reference    reference;
 
          typedef typename BaseContainer::const_reference const_reference;
 
 
          typedef W weight_type;
 
          typedef WeightedMatrixAccumulationView<DiagonalLocalMatrix> WeightedAccumulationView;
 
          struct iterator
            : public Dune::BidirectionalIteratorFacade<iterator,value_type>
          {
 
            iterator()
              : _m(nullptr)
              , _i(0)
            {}
 
            iterator(DiagonalLocalMatrix& m, size_type i)
              : _m(&m)
              , _i(i)
            {}
 
            bool equals(const iterator& other) const
            {
              return _m == other._m && _i == other._i;
            }
 
            value_type& dereference() const
            {
              return _m->getEntry(_i,_i);
            }
 
            void increment()
            {
              ++_i;
            }
 
            void decrement()
            {
              --_i;
            }
 
            size_type row() const
            {
              return _i;
            }
 
            size_type col() const
            {
              return _i;
            }
 
            DiagonalLocalMatrix* _m;
            size_type _i;
 
          };
 
          iterator begin()
          {
            return iterator(*this,0);
          }
 
          iterator end()
          {
            return iterator(*this,nrows());
          }
 
          DiagonalLocalMatrix () {}
 
          DiagonalLocalMatrix (size_type r, size_type c)
            : _container(r)
            , _rows(r)
            , _cols(c)
          {
            assert(r == c);
          }
 
          DiagonalLocalMatrix (size_type r, size_type c, const T& t)
            : _container(r,t)
            , _rows(r)
            , _cols(c)
          {
            assert(r == c);
          }
 
          void resize (size_type r, size_type c)
          {
            assert(r == c);
            _container.resize(r);
            _rows = r;
            _cols = c;
          }
 
          DiagonalLocalMatrix& operator=(const T& t)
          {
            std::fill(_container.begin(),_container.end(),t);
            return *this;
          }
 
          void assign (size_type r, size_type c, const T& t)
          {
            assert(r == c);
            _container.assign(r,t);
            _rows = r;
            _cols = c;
          }
 
 
          template<typename LFSU, typename LFSV>
          T& operator() (const LFSV& lfsv, size_type i, const LFSU& lfsu, size_type j)
          {
            return getEntry(lfsv.localIndex(i),lfsu.localIndex(j));
          }
 
 
          template<typename LFSU, typename LFSV>
          const T& operator() (const LFSV& lfsv, size_type i, const LFSU& lfsu, size_type j) const
          {
            return getEntry(lfsv.localIndex(i),lfsu.localIndex(j));
          }
 
          DiagonalLocalMatrix& operator *= (const T& x)
          {
            using namespace std::placeholders;
            std::transform(
              _container.begin(),
              _container.end(),
              _container.begin(),
              std::bind(std::multiplies<T>(),x,_1)
              );
            return *this;
          }
 
          size_type nrows () const
          {
            return _rows;
          }
 
          size_type ncols () const
          {
            return _cols;
          }
 
          template<class X, class R>
          void umv (const X& x, R& y) const
          {
            for (size_type i=0; i<_rows; ++i)
              {
                for (size_type j=0; j<_cols; j++)
                  accessBaseContainer(y)[i] += getEntry(i,j) * accessBaseContainer(x)[j];
              }
          }
 
          template<class X, class R>
          void usmv (const value_type& alpha, const X& x, R& y) const
          {
            for (size_type i=0; i<_rows; ++i)
              {
                for (size_type j=0; j<_cols; j++)
                  accessBaseContainer(y)[i] += alpha * getEntry(i,j) * accessBaseContainer(x)[j];
              }
          }
 
          WeightedAccumulationView weightedAccumulationView(weight_type weight)
          {
            return WeightedAccumulationView(*this,weight);
          }
 
 
          BaseContainer& base()
          {
            return _container;
          }
 
 
          const BaseContainer& base() const
          {
            return _container;
          }
 
 
          value_type& getEntry(size_type i, size_type j)
          {
            assert(i == j);
            return _container[i];
          }
 
 
          const value_type& getEntry(size_type i, size_type j) const
          {
            assert(i == j);
            return _container[i];
          }
 
        private:
 
          std::vector<T> _container;
          size_type _rows, _cols;
     };
 
  } // namespace PDELab
} // namespace Dune
 
#endif // DUNE_PDELAB_GRIDOPERATOR_COMMON_DIAGONALLOCALMATRIX_HH