1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr> 5 // Copyright (C) 2010 Daniel Lowengrub <lowdanie@gmail.com> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_SPARSEVIEW_H 12 #define EIGEN_SPARSEVIEW_H 13 14 namespace Eigen { 15 16 namespace internal { 17 18 template<typename MatrixType> 19 struct traits<SparseView<MatrixType> > : traits<MatrixType> 20 { 21 typedef typename MatrixType::StorageIndex StorageIndex; 22 typedef Sparse StorageKind; 23 enum { 24 Flags = int(traits<MatrixType>::Flags) & (RowMajorBit) 25 }; 26 }; 27 28 } // end namespace internal 29 30 /** \ingroup SparseCore_Module 31 * \class SparseView 32 * 33 * \brief Expression of a dense or sparse matrix with zero or too small values removed 34 * 35 * \tparam MatrixType the type of the object of which we are removing the small entries 36 * 37 * This class represents an expression of a given dense or sparse matrix with 38 * entries smaller than \c reference * \c epsilon are removed. 39 * It is the return type of MatrixBase::sparseView() and SparseMatrixBase::pruned() 40 * and most of the time this is the only way it is used. 41 * 42 * \sa MatrixBase::sparseView(), SparseMatrixBase::pruned() 43 */ 44 template<typename MatrixType> 45 class SparseView : public SparseMatrixBase<SparseView<MatrixType> > 46 { 47 typedef typename MatrixType::Nested MatrixTypeNested; 48 typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested; 49 typedef SparseMatrixBase<SparseView > Base; 50 public: 51 EIGEN_SPARSE_PUBLIC_INTERFACE(SparseView) 52 typedef typename internal::remove_all<MatrixType>::type NestedExpression; 53 54 explicit SparseView(const MatrixType& mat, const Scalar& reference = Scalar(0), 55 const RealScalar &epsilon = NumTraits<Scalar>::dummy_precision()) 56 : m_matrix(mat), m_reference(reference), m_epsilon(epsilon) {} 57 58 inline Index rows() const { return m_matrix.rows(); } 59 inline Index cols() const { return m_matrix.cols(); } 60 61 inline Index innerSize() const { return m_matrix.innerSize(); } 62 inline Index outerSize() const { return m_matrix.outerSize(); } 63 64 /** \returns the nested expression */ 65 const typename internal::remove_all<MatrixTypeNested>::type& 66 nestedExpression() const { return m_matrix; } 67 68 Scalar reference() const { return m_reference; } 69 RealScalar epsilon() const { return m_epsilon; } 70 71 protected: 72 MatrixTypeNested m_matrix; 73 Scalar m_reference; 74 RealScalar m_epsilon; 75 }; 76 77 namespace internal { 78 79 // TODO find a way to unify the two following variants 80 // This is tricky because implementing an inner iterator on top of an IndexBased evaluator is 81 // not easy because the evaluators do not expose the sizes of the underlying expression. 82 83 template<typename ArgType> 84 struct unary_evaluator<SparseView<ArgType>, IteratorBased> 85 : public evaluator_base<SparseView<ArgType> > 86 { 87 typedef typename evaluator<ArgType>::InnerIterator EvalIterator; 88 public: 89 typedef SparseView<ArgType> XprType; 90 91 class InnerIterator : public EvalIterator 92 { 93 typedef typename XprType::Scalar Scalar; 94 public: 95 96 EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& sve, Index outer) 97 : EvalIterator(sve.m_argImpl,outer), m_view(sve.m_view) 98 { 99 incrementToNonZero(); 100 } 101 102 EIGEN_STRONG_INLINE InnerIterator& operator++() 103 { 104 EvalIterator::operator++(); 105 incrementToNonZero(); 106 return *this; 107 } 108 109 using EvalIterator::value; 110 111 protected: 112 const XprType &m_view; 113 114 private: 115 void incrementToNonZero() 116 { 117 while((bool(*this)) && internal::isMuchSmallerThan(value(), m_view.reference(), m_view.epsilon())) 118 { 119 EvalIterator::operator++(); 120 } 121 } 122 }; 123 124 enum { 125 CoeffReadCost = evaluator<ArgType>::CoeffReadCost, 126 Flags = XprType::Flags 127 }; 128 129 explicit unary_evaluator(const XprType& xpr) : m_argImpl(xpr.nestedExpression()), m_view(xpr) {} 130 131 protected: 132 evaluator<ArgType> m_argImpl; 133 const XprType &m_view; 134 }; 135 136 template<typename ArgType> 137 struct unary_evaluator<SparseView<ArgType>, IndexBased> 138 : public evaluator_base<SparseView<ArgType> > 139 { 140 public: 141 typedef SparseView<ArgType> XprType; 142 protected: 143 enum { IsRowMajor = (XprType::Flags&RowMajorBit)==RowMajorBit }; 144 typedef typename XprType::Scalar Scalar; 145 typedef typename XprType::StorageIndex StorageIndex; 146 public: 147 148 class InnerIterator 149 { 150 public: 151 152 EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& sve, Index outer) 153 : m_sve(sve), m_inner(0), m_outer(outer), m_end(sve.m_view.innerSize()) 154 { 155 incrementToNonZero(); 156 } 157 158 EIGEN_STRONG_INLINE InnerIterator& operator++() 159 { 160 m_inner++; 161 incrementToNonZero(); 162 return *this; 163 } 164 165 EIGEN_STRONG_INLINE Scalar value() const 166 { 167 return (IsRowMajor) ? m_sve.m_argImpl.coeff(m_outer, m_inner) 168 : m_sve.m_argImpl.coeff(m_inner, m_outer); 169 } 170 171 EIGEN_STRONG_INLINE StorageIndex index() const { return m_inner; } 172 inline Index row() const { return IsRowMajor ? m_outer : index(); } 173 inline Index col() const { return IsRowMajor ? index() : m_outer; } 174 175 EIGEN_STRONG_INLINE operator bool() const { return m_inner < m_end && m_inner>=0; } 176 177 protected: 178 const unary_evaluator &m_sve; 179 Index m_inner; 180 const Index m_outer; 181 const Index m_end; 182 183 private: 184 void incrementToNonZero() 185 { 186 while((bool(*this)) && internal::isMuchSmallerThan(value(), m_sve.m_view.reference(), m_sve.m_view.epsilon())) 187 { 188 m_inner++; 189 } 190 } 191 }; 192 193 enum { 194 CoeffReadCost = evaluator<ArgType>::CoeffReadCost, 195 Flags = XprType::Flags 196 }; 197 198 explicit unary_evaluator(const XprType& xpr) : m_argImpl(xpr.nestedExpression()), m_view(xpr) {} 199 200 protected: 201 evaluator<ArgType> m_argImpl; 202 const XprType &m_view; 203 }; 204 205 } // end namespace internal 206 207 /** \ingroup SparseCore_Module 208 * 209 * \returns a sparse expression of the dense expression \c *this with values smaller than 210 * \a reference * \a epsilon removed. 211 * 212 * This method is typically used when prototyping to convert a quickly assembled dense Matrix \c D to a SparseMatrix \c S: 213 * \code 214 * MatrixXd D(n,m); 215 * SparseMatrix<double> S; 216 * S = D.sparseView(); // suppress numerical zeros (exact) 217 * S = D.sparseView(reference); 218 * S = D.sparseView(reference,epsilon); 219 * \endcode 220 * where \a reference is a meaningful non zero reference value, 221 * and \a epsilon is a tolerance factor defaulting to NumTraits<Scalar>::dummy_precision(). 222 * 223 * \sa SparseMatrixBase::pruned(), class SparseView */ 224 template<typename Derived> 225 const SparseView<Derived> MatrixBase<Derived>::sparseView(const Scalar& reference, 226 const typename NumTraits<Scalar>::Real& epsilon) const 227 { 228 return SparseView<Derived>(derived(), reference, epsilon); 229 } 230 231 /** \returns an expression of \c *this with values smaller than 232 * \a reference * \a epsilon removed. 233 * 234 * This method is typically used in conjunction with the product of two sparse matrices 235 * to automatically prune the smallest values as follows: 236 * \code 237 * C = (A*B).pruned(); // suppress numerical zeros (exact) 238 * C = (A*B).pruned(ref); 239 * C = (A*B).pruned(ref,epsilon); 240 * \endcode 241 * where \c ref is a meaningful non zero reference value. 242 * */ 243 template<typename Derived> 244 const SparseView<Derived> 245 SparseMatrixBase<Derived>::pruned(const Scalar& reference, 246 const RealScalar& epsilon) const 247 { 248 return SparseView<Derived>(derived(), reference, epsilon); 249 } 250 251 } // end namespace Eigen 252 253 #endif 254