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