diff --git a/dune/matrix-vector/singlenonzerorowmatrix.hh b/dune/matrix-vector/singlenonzerorowmatrix.hh index 1dd1bc7567ac4ef10eee9b5ce433576e8cae1a70..5de2968a495103c8bc4040c297220f1af0fb0b31 100644 --- a/dune/matrix-vector/singlenonzerorowmatrix.hh +++ b/dune/matrix-vector/singlenonzerorowmatrix.hh @@ -8,7 +8,6 @@ #include <dune/fufem/arithmetic.hh> #include <dune/fufem/indexedsliceiterator.hh> - /** * \brief A static matrix that has only a single nonzero row * @@ -17,39 +16,81 @@ * all entries accordingly. So it will not reduce memory * requirements but allows to implement methods more efficiently. */ -template<class K, int ROWS, int COLS> -class SingleNonZeroRowMatrix : - public Dune::FieldMatrix<K, ROWS, COLS> +template<class K, int ROWS, int COLS> //TODO allow to set storage type as in singlenonzerocolumnmatrix. Update docu accordingly. +class SingleNonZeroRowMatrix { - typedef typename Dune::FieldMatrix<K, ROWS, COLS> Base; +public: + typedef typename Dune::FieldMatrix<K, 1, COLS> SingleRowMatrix; + +protected: + class RowProxy + { + public: + typedef std::size_t size_type; + typedef typename SingleRowMatrix::ConstColIterator ConstIterator; + typedef ConstIterator const_iterator; + + RowProxy(const SingleRowMatrix* nzRow, size_type rowIndex, size_type nzRowIndex) : + nzRow_(nzRow), + isNonZero_(rowIndex==nzRowIndex) + {} + + ConstIterator begin() const + { + if (isNonZero_) + return (*nzRow_)[0].begin(); + else + return (*nzRow_)[0].end(); + } + + ConstIterator end() const + { + return (*nzRow_)[0].end(); + } + + protected: + const SingleRowMatrix* nzRow_; + bool isNonZero_; + }; public: + enum { rows = ROWS, cols = COLS}; - typedef typename Dune::FieldMatrix<K, 1, Base::cols> SingleRowMatrix; + typedef RowProxy row_type; + typedef row_type const_row_reference; + typedef typename std::size_t size_type; + typedef typename RowProxy::ConstIterator ConstColIterator; /** * \brief Create from single row matrix and row index */ - SingleNonZeroRowMatrix(const SingleRowMatrix& c, int row) : - Base(0.0), - row_(row) + SingleNonZeroRowMatrix(const SingleRowMatrix& r, size_type nzRowIndex) : + nzRow_(r), + nzRowIndex_(nzRowIndex) + {} + + size_type N() const + { + return ROWS; + } + + size_type M() const { - for(int i=0; i<Base::cols; ++i) - (*this)[row_][i] = c[0][i]; + return COLS; } template<class X , class Y > void umv(const X& x, Y& y) const { - for(int i=0; i<Base::cols; ++i) - y[row_] += (*this)[row_][i] * x[i]; + for(size_type i=0; i<M(); ++i) + y[nzRowIndex_] += nzRow_[0][i] * x[i]; } template<class X , class Y > void umtv(const X& x, Y& y) const { - for(int i=0; i<Base::cols; ++i) - y[i] = (*this)[row_][i] * x[row_]; + for(size_type i=0; i<N(); ++i) + y[i] = nzRow_[0][i] * x[nzRowIndex_]; } template<class X , class Y > @@ -59,19 +100,19 @@ public: umtv(x, y); } - int nonZeroRowIndex() const + const_row_reference operator[] (size_type rowIndex) const { - return row_; + return const_row_reference(&nzRow_, rowIndex, nzRowIndex_); } - const Base::row_type& nonZeroRow() const + size_type nonZeroRowIndex() const { - return (*this)[row_]; + return nzRowIndex_; } protected: - - const int row_; + SingleRowMatrix nzRow_; + const size_type nzRowIndex_; }; @@ -81,27 +122,14 @@ namespace Arithmetic template<class K, int ROWS, int COLS> struct MatrixTraits<SingleNonZeroRowMatrix<K, ROWS, COLS> > { - enum { isMatrix=true }; - enum { rows=ROWS}; - enum { cols=COLS}; - typedef typename SingleNonZeroRowMatrix<K, ROWS, COLS>::row_type::ConstColIterator ConstColIterator; - - static ConstColIterator rowBegin(const typename SingleNonZeroRowMatrix<K, ROWS, COLS>::row_type& m, int row) - { - if (row == m.nonZeroRowIndex()) - return m.nonZeroRow().begin(); - else - return m.nonZeroRow().end(); - } - - static ConstColIterator rowEnd(const typename SingleNonZeroRowMatrix<K, ROWS, COLS>::row_type& m, int row) - { - return m.nonZeroRow().end(); - } + enum { isMatrix = true }; + enum { isDense = false }; + enum { sizeIsStatic = true }; //TODO only one of these should be used + enum { hasStaticSize = true }; //TODO only one of these should be used + enum { rows = ROWS}; + enum { cols = COLS}; }; - -}; - +} #endif