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
14namespace Eigen {
15
16namespace internal {
17
18template<typename MatrixType>
19struct 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 */
44template<typename MatrixType>
45class 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;
50public:
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
71protected:
72 MatrixTypeNested m_matrix;
73 Scalar m_reference;
74 RealScalar m_epsilon;
75};
76
77namespace 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
83template<typename ArgType>
84struct 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
136template<typename ArgType>
137struct 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 */
224template<typename Derived>
225const 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 * */
243template<typename Derived>
244const SparseView<Derived>
245SparseMatrixBase<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