SparseView.h (8127B)
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 protected: 94 typedef typename XprType::Scalar Scalar; 95 public: 96 97 EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& sve, Index outer) 98 : EvalIterator(sve.m_argImpl,outer), m_view(sve.m_view) 99 { 100 incrementToNonZero(); 101 } 102 103 EIGEN_STRONG_INLINE InnerIterator& operator++() 104 { 105 EvalIterator::operator++(); 106 incrementToNonZero(); 107 return *this; 108 } 109 110 using EvalIterator::value; 111 112 protected: 113 const XprType &m_view; 114 115 private: 116 void incrementToNonZero() 117 { 118 while((bool(*this)) && internal::isMuchSmallerThan(value(), m_view.reference(), m_view.epsilon())) 119 { 120 EvalIterator::operator++(); 121 } 122 } 123 }; 124 125 enum { 126 CoeffReadCost = evaluator<ArgType>::CoeffReadCost, 127 Flags = XprType::Flags 128 }; 129 130 explicit unary_evaluator(const XprType& xpr) : m_argImpl(xpr.nestedExpression()), m_view(xpr) {} 131 132 protected: 133 evaluator<ArgType> m_argImpl; 134 const XprType &m_view; 135 }; 136 137 template<typename ArgType> 138 struct unary_evaluator<SparseView<ArgType>, IndexBased> 139 : public evaluator_base<SparseView<ArgType> > 140 { 141 public: 142 typedef SparseView<ArgType> XprType; 143 protected: 144 enum { IsRowMajor = (XprType::Flags&RowMajorBit)==RowMajorBit }; 145 typedef typename XprType::Scalar Scalar; 146 typedef typename XprType::StorageIndex StorageIndex; 147 public: 148 149 class InnerIterator 150 { 151 public: 152 153 EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& sve, Index outer) 154 : m_sve(sve), m_inner(0), m_outer(outer), m_end(sve.m_view.innerSize()) 155 { 156 incrementToNonZero(); 157 } 158 159 EIGEN_STRONG_INLINE InnerIterator& operator++() 160 { 161 m_inner++; 162 incrementToNonZero(); 163 return *this; 164 } 165 166 EIGEN_STRONG_INLINE Scalar value() const 167 { 168 return (IsRowMajor) ? m_sve.m_argImpl.coeff(m_outer, m_inner) 169 : m_sve.m_argImpl.coeff(m_inner, m_outer); 170 } 171 172 EIGEN_STRONG_INLINE StorageIndex index() const { return m_inner; } 173 inline Index row() const { return IsRowMajor ? m_outer : index(); } 174 inline Index col() const { return IsRowMajor ? index() : m_outer; } 175 176 EIGEN_STRONG_INLINE operator bool() const { return m_inner < m_end && m_inner>=0; } 177 178 protected: 179 const unary_evaluator &m_sve; 180 Index m_inner; 181 const Index m_outer; 182 const Index m_end; 183 184 private: 185 void incrementToNonZero() 186 { 187 while((bool(*this)) && internal::isMuchSmallerThan(value(), m_sve.m_view.reference(), m_sve.m_view.epsilon())) 188 { 189 m_inner++; 190 } 191 } 192 }; 193 194 enum { 195 CoeffReadCost = evaluator<ArgType>::CoeffReadCost, 196 Flags = XprType::Flags 197 }; 198 199 explicit unary_evaluator(const XprType& xpr) : m_argImpl(xpr.nestedExpression()), m_view(xpr) {} 200 201 protected: 202 evaluator<ArgType> m_argImpl; 203 const XprType &m_view; 204 }; 205 206 } // end namespace internal 207 208 /** \ingroup SparseCore_Module 209 * 210 * \returns a sparse expression of the dense expression \c *this with values smaller than 211 * \a reference * \a epsilon removed. 212 * 213 * This method is typically used when prototyping to convert a quickly assembled dense Matrix \c D to a SparseMatrix \c S: 214 * \code 215 * MatrixXd D(n,m); 216 * SparseMatrix<double> S; 217 * S = D.sparseView(); // suppress numerical zeros (exact) 218 * S = D.sparseView(reference); 219 * S = D.sparseView(reference,epsilon); 220 * \endcode 221 * where \a reference is a meaningful non zero reference value, 222 * and \a epsilon is a tolerance factor defaulting to NumTraits<Scalar>::dummy_precision(). 223 * 224 * \sa SparseMatrixBase::pruned(), class SparseView */ 225 template<typename Derived> 226 const SparseView<Derived> MatrixBase<Derived>::sparseView(const Scalar& reference, 227 const typename NumTraits<Scalar>::Real& epsilon) const 228 { 229 return SparseView<Derived>(derived(), reference, epsilon); 230 } 231 232 /** \returns an expression of \c *this with values smaller than 233 * \a reference * \a epsilon removed. 234 * 235 * This method is typically used in conjunction with the product of two sparse matrices 236 * to automatically prune the smallest values as follows: 237 * \code 238 * C = (A*B).pruned(); // suppress numerical zeros (exact) 239 * C = (A*B).pruned(ref); 240 * C = (A*B).pruned(ref,epsilon); 241 * \endcode 242 * where \c ref is a meaningful non zero reference value. 243 * */ 244 template<typename Derived> 245 const SparseView<Derived> 246 SparseMatrixBase<Derived>::pruned(const Scalar& reference, 247 const RealScalar& epsilon) const 248 { 249 return SparseView<Derived>(derived(), reference, epsilon); 250 } 251 252 } // end namespace Eigen 253 254 #endif