src/SparseCore/SparseView.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 // Copyright (C) 2010 Daniel Lowengrub <lowdanie@gmail.com>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_SPARSEVIEW_H
0012 #define EIGEN_SPARSEVIEW_H
0013
0014 namespace Eigen {
0015
0016 namespace internal {
0017
0018 template<typename MatrixType>
0019 struct traits<SparseView<MatrixType> > : traits<MatrixType>
0020 {
0021   typedef typename MatrixType::StorageIndex StorageIndex;
0022   typedef Sparse StorageKind;
0023   enum {
0024     Flags = int(traits<MatrixType>::Flags) & (RowMajorBit)
0025   };
0026 };
0027
0028 } // end namespace internal
0029
0030 /** \ingroup SparseCore_Module
0031   * \class SparseView
0032   *
0033   * \brief Expression of a dense or sparse matrix with zero or too small values removed
0034   *
0035   * \tparam MatrixType the type of the object of which we are removing the small entries
0036   *
0037   * This class represents an expression of a given dense or sparse matrix with
0038   * entries smaller than \c reference * \c epsilon are removed.
0039   * It is the return type of MatrixBase::sparseView() and SparseMatrixBase::pruned()
0040   * and most of the time this is the only way it is used.
0041   *
0042   * \sa MatrixBase::sparseView(), SparseMatrixBase::pruned()
0043   */
0044 template<typename MatrixType>
0045 class SparseView : public SparseMatrixBase<SparseView<MatrixType> >
0046 {
0047   typedef typename MatrixType::Nested MatrixTypeNested;
0048   typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
0049   typedef SparseMatrixBase<SparseView > Base;
0050 public:
0051   EIGEN_SPARSE_PUBLIC_INTERFACE(SparseView)
0052   typedef typename internal::remove_all<MatrixType>::type NestedExpression;
0053
0054   explicit SparseView(const MatrixType& mat, const Scalar& reference = Scalar(0),
0055                       const RealScalar &epsilon = NumTraits<Scalar>::dummy_precision())
0056     : m_matrix(mat), m_reference(reference), m_epsilon(epsilon) {}
0057
0058   inline Index rows() const { return m_matrix.rows(); }
0059   inline Index cols() const { return m_matrix.cols(); }
0060
0061   inline Index innerSize() const { return m_matrix.innerSize(); }
0062   inline Index outerSize() const { return m_matrix.outerSize(); }
0063
0064   /** \returns the nested expression */
0065   const typename internal::remove_all<MatrixTypeNested>::type&
0066   nestedExpression() const { return m_matrix; }
0067
0068   Scalar reference() const { return m_reference; }
0069   RealScalar epsilon() const { return m_epsilon; }
0070
0071 protected:
0072   MatrixTypeNested m_matrix;
0073   Scalar m_reference;
0074   RealScalar m_epsilon;
0075 };
0076
0077 namespace internal {
0078
0079 // TODO find a way to unify the two following variants
0080 // This is tricky because implementing an inner iterator on top of an IndexBased evaluator is
0081 // not easy because the evaluators do not expose the sizes of the underlying expression.
0082
0083 template<typename ArgType>
0084 struct unary_evaluator<SparseView<ArgType>, IteratorBased>
0085   : public evaluator_base<SparseView<ArgType> >
0086 {
0087     typedef typename evaluator<ArgType>::InnerIterator EvalIterator;
0088   public:
0089     typedef SparseView<ArgType> XprType;
0090
0091     class InnerIterator : public EvalIterator
0092     {
0093       protected:
0094         typedef typename XprType::Scalar Scalar;
0095       public:
0096
0097         EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& sve, Index outer)
0098           : EvalIterator(sve.m_argImpl,outer), m_view(sve.m_view)
0099         {
0100           incrementToNonZero();
0101         }
0102
0103         EIGEN_STRONG_INLINE InnerIterator& operator++()
0104         {
0105           EvalIterator::operator++();
0106           incrementToNonZero();
0107           return *this;
0108         }
0109
0110         using EvalIterator::value;
0111
0112       protected:
0113         const XprType &m_view;
0114
0115       private:
0116         void incrementToNonZero()
0117         {
0118           while((bool(*this)) && internal::isMuchSmallerThan(value(), m_view.reference(), m_view.epsilon()))
0119           {
0120             EvalIterator::operator++();
0121           }
0122         }
0123     };
0124
0125     enum {
0126       CoeffReadCost = evaluator<ArgType>::CoeffReadCost,
0127       Flags = XprType::Flags
0128     };
0129
0130     explicit unary_evaluator(const XprType& xpr) : m_argImpl(xpr.nestedExpression()), m_view(xpr) {}
0131
0132   protected:
0133     evaluator<ArgType> m_argImpl;
0134     const XprType &m_view;
0135 };
0136
0137 template<typename ArgType>
0138 struct unary_evaluator<SparseView<ArgType>, IndexBased>
0139   : public evaluator_base<SparseView<ArgType> >
0140 {
0141   public:
0142     typedef SparseView<ArgType> XprType;
0143   protected:
0144     enum { IsRowMajor = (XprType::Flags&RowMajorBit)==RowMajorBit };
0145     typedef typename XprType::Scalar Scalar;
0146     typedef typename XprType::StorageIndex StorageIndex;
0147   public:
0148
0149     class InnerIterator
0150     {
0151       public:
0152
0153         EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& sve, Index outer)
0154           : m_sve(sve), m_inner(0), m_outer(outer), m_end(sve.m_view.innerSize())
0155         {
0156           incrementToNonZero();
0157         }
0158
0159         EIGEN_STRONG_INLINE InnerIterator& operator++()
0160         {
0161           m_inner++;
0162           incrementToNonZero();
0163           return *this;
0164         }
0165
0166         EIGEN_STRONG_INLINE Scalar value() const
0167         {
0168           return (IsRowMajor) ? m_sve.m_argImpl.coeff(m_outer, m_inner)
0169                               : m_sve.m_argImpl.coeff(m_inner, m_outer);
0170         }
0171
0172         EIGEN_STRONG_INLINE StorageIndex index() const { return m_inner; }
0173         inline Index row() const { return IsRowMajor ? m_outer : index(); }
0174         inline Index col() const { return IsRowMajor ? index() : m_outer; }
0175
0176         EIGEN_STRONG_INLINE operator bool() const { return m_inner < m_end && m_inner>=0; }
0177
0178       protected:
0179         const unary_evaluator &m_sve;
0180         Index m_inner;
0181         const Index m_outer;
0182         const Index m_end;
0183
0184       private:
0185         void incrementToNonZero()
0186         {
0187           while((bool(*this)) && internal::isMuchSmallerThan(value(), m_sve.m_view.reference(), m_sve.m_view.epsilon()))
0188           {
0189             m_inner++;
0190           }
0191         }
0192     };
0193
0194     enum {
0195       CoeffReadCost = evaluator<ArgType>::CoeffReadCost,
0196       Flags = XprType::Flags
0197     };
0198
0199     explicit unary_evaluator(const XprType& xpr) : m_argImpl(xpr.nestedExpression()), m_view(xpr) {}
0200
0201   protected:
0202     evaluator<ArgType> m_argImpl;
0203     const XprType &m_view;
0204 };
0205
0206 } // end namespace internal
0207
0208 /** \ingroup SparseCore_Module
0209   *
0210   * \returns a sparse expression of the dense expression \c *this with values smaller than
0211   * \a reference * \a epsilon removed.
0212   *
0213   * This method is typically used when prototyping to convert a quickly assembled dense Matrix \c D to a SparseMatrix \c S:
0214   * \code
0215   * MatrixXd D(n,m);
0216   * SparseMatrix<double> S;
0217   * S = D.sparseView();             // suppress numerical zeros (exact)
0218   * S = D.sparseView(reference);
0219   * S = D.sparseView(reference,epsilon);
0220   * \endcode
0221   * where \a reference is a meaningful non zero reference value,
0222   * and \a epsilon is a tolerance factor defaulting to NumTraits<Scalar>::dummy_precision().
0223   *
0224   * \sa SparseMatrixBase::pruned(), class SparseView */
0225 template<typename Derived>
0226 const SparseView<Derived> MatrixBase<Derived>::sparseView(const Scalar& reference,
0227                                                           const typename NumTraits<Scalar>::Real& epsilon) const
0228 {
0229   return SparseView<Derived>(derived(), reference, epsilon);
0230 }
0231
0232 /** \returns an expression of \c *this with values smaller than
0233   * \a reference * \a epsilon removed.
0234   *
0235   * This method is typically used in conjunction with the product of two sparse matrices
0236   * to automatically prune the smallest values as follows:
0237   * \code
0238   * C = (A*B).pruned();             // suppress numerical zeros (exact)
0239   * C = (A*B).pruned(ref);
0240   * C = (A*B).pruned(ref,epsilon);
0241   * \endcode
0242   * where \c ref is a meaningful non zero reference value.
0243   * */
0244 template<typename Derived>
0245 const SparseView<Derived>
0246 SparseMatrixBase<Derived>::pruned(const Scalar& reference,
0247                                   const RealScalar& epsilon) const
0248 {
0249   return SparseView<Derived>(derived(), reference, epsilon);
0250 }
0251
0252 } // end namespace Eigen
0253
0254 #endif