Eigen/src/SparseCore/SparseView.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2010 Daniel Lowengrub <lowdanie@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_SPARSEVIEW_H
#define EIGEN_SPARSEVIEW_H

// IWYU pragma: private
#include "./InternalHeaderCheck.h"

namespace Eigen {

namespace internal {

template <typename MatrixType>
struct traits<SparseView<MatrixType> > : traits<MatrixType> {
  typedef typename MatrixType::StorageIndex StorageIndex;
  typedef Sparse StorageKind;
  enum { Flags = int(traits<MatrixType>::Flags) & (RowMajorBit) };
};

}  // end namespace internal

/** \ingroup SparseCore_Module
 * \class SparseView
 *
 * \brief Expression of a dense or sparse matrix with zero or too small values removed
 *
 * \tparam MatrixType the type of the object of which we are removing the small entries
 *
 * This class represents an expression of a given dense or sparse matrix with
 * entries smaller than \c reference * \c epsilon are removed.
 * It is the return type of MatrixBase::sparseView() and SparseMatrixBase::pruned()
 * and most of the time this is the only way it is used.
 *
 * \sa MatrixBase::sparseView(), SparseMatrixBase::pruned()
 */
template <typename MatrixType>
class SparseView : public SparseMatrixBase<SparseView<MatrixType> > {
  typedef typename MatrixType::Nested MatrixTypeNested;
  typedef internal::remove_all_t<MatrixTypeNested> MatrixTypeNested_;
  typedef SparseMatrixBase<SparseView> Base;

 public:
  EIGEN_SPARSE_PUBLIC_INTERFACE(SparseView)
  typedef internal::remove_all_t<MatrixType> NestedExpression;

  explicit SparseView(const MatrixType& mat, const Scalar& reference = Scalar(0),
                      const RealScalar& epsilon = NumTraits<Scalar>::dummy_precision())
      : m_matrix(mat), m_reference(reference), m_epsilon(epsilon) {}

  inline Index rows() const { return m_matrix.rows(); }
  inline Index cols() const { return m_matrix.cols(); }

  inline Index innerSize() const { return m_matrix.innerSize(); }
  inline Index outerSize() const { return m_matrix.outerSize(); }

  /** \returns the nested expression */
  const internal::remove_all_t<MatrixTypeNested>& nestedExpression() const { return m_matrix; }

  Scalar reference() const { return m_reference; }
  RealScalar epsilon() const { return m_epsilon; }

 protected:
  MatrixTypeNested m_matrix;
  Scalar m_reference;
  RealScalar m_epsilon;
};

namespace internal {

// TODO find a way to unify the two following variants
// This is tricky because implementing an inner iterator on top of an IndexBased evaluator is
// not easy because the evaluators do not expose the sizes of the underlying expression.

template <typename ArgType>
struct unary_evaluator<SparseView<ArgType>, IteratorBased> : public evaluator_base<SparseView<ArgType> > {
  typedef typename evaluator<ArgType>::InnerIterator EvalIterator;

 public:
  typedef SparseView<ArgType> XprType;

  class InnerIterator : public EvalIterator {
   protected:
    typedef typename XprType::Scalar Scalar;

   public:
    EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& sve, Index outer)
        : EvalIterator(sve.m_argImpl, outer), m_view(sve.m_view) {
      incrementToNonZero();
    }

    EIGEN_STRONG_INLINE InnerIterator& operator++() {
      EvalIterator::operator++();
      incrementToNonZero();
      return *this;
    }

    using EvalIterator::value;

   protected:
    const XprType& m_view;

   private:
    void incrementToNonZero() {
      while ((bool(*this)) && internal::isMuchSmallerThan(value(), m_view.reference(), m_view.epsilon())) {
        EvalIterator::operator++();
      }
    }
  };

  enum { CoeffReadCost = evaluator<ArgType>::CoeffReadCost, Flags = XprType::Flags };

  explicit unary_evaluator(const XprType& xpr) : m_argImpl(xpr.nestedExpression()), m_view(xpr) {}

 protected:
  evaluator<ArgType> m_argImpl;
  const XprType& m_view;
};

template <typename ArgType>
struct unary_evaluator<SparseView<ArgType>, IndexBased> : public evaluator_base<SparseView<ArgType> > {
 public:
  typedef SparseView<ArgType> XprType;

 protected:
  enum { IsRowMajor = (XprType::Flags & RowMajorBit) == RowMajorBit };
  typedef typename XprType::Scalar Scalar;
  typedef typename XprType::StorageIndex StorageIndex;

 public:
  class InnerIterator {
   public:
    EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& sve, Index outer)
        : m_sve(sve), m_inner(0), m_outer(outer), m_end(sve.m_view.innerSize()) {
      incrementToNonZero();
    }

    EIGEN_STRONG_INLINE InnerIterator& operator++() {
      m_inner++;
      incrementToNonZero();
      return *this;
    }

    EIGEN_STRONG_INLINE Scalar value() const {
      return (IsRowMajor) ? m_sve.m_argImpl.coeff(m_outer, m_inner) : m_sve.m_argImpl.coeff(m_inner, m_outer);
    }

    EIGEN_STRONG_INLINE StorageIndex index() const { return m_inner; }
    inline Index row() const { return IsRowMajor ? m_outer : index(); }
    inline Index col() const { return IsRowMajor ? index() : m_outer; }

    EIGEN_STRONG_INLINE operator bool() const { return m_inner < m_end && m_inner >= 0; }

   protected:
    const unary_evaluator& m_sve;
    Index m_inner;
    const Index m_outer;
    const Index m_end;

   private:
    void incrementToNonZero() {
      while ((bool(*this)) && internal::isMuchSmallerThan(value(), m_sve.m_view.reference(), m_sve.m_view.epsilon())) {
        m_inner++;
      }
    }
  };

  enum { CoeffReadCost = evaluator<ArgType>::CoeffReadCost, Flags = XprType::Flags };

  explicit unary_evaluator(const XprType& xpr) : m_argImpl(xpr.nestedExpression()), m_view(xpr) {}

 protected:
  evaluator<ArgType> m_argImpl;
  const XprType& m_view;
};

}  // end namespace internal

/** \ingroup SparseCore_Module
 *
 * \returns a sparse expression of the dense expression \c *this with values smaller than
 * \a reference * \a epsilon removed.
 *
 * This method is typically used when prototyping to convert a quickly assembled dense Matrix \c D to a SparseMatrix \c
 * S: \code MatrixXd D(n,m); SparseMatrix<double> S; S = D.sparseView();             // suppress numerical zeros (exact)
 * S = D.sparseView(reference);
 * S = D.sparseView(reference,epsilon);
 * \endcode
 * where \a reference is a meaningful non zero reference value,
 * and \a epsilon is a tolerance factor defaulting to NumTraits<Scalar>::dummy_precision().
 *
 * \sa SparseMatrixBase::pruned(), class SparseView */
template <typename Derived>
const SparseView<Derived> MatrixBase<Derived>::sparseView(const Scalar& reference,
                                                          const typename NumTraits<Scalar>::Real& epsilon) const {
  return SparseView<Derived>(derived(), reference, epsilon);
}

/** \returns an expression of \c *this with values smaller than
 * \a reference * \a epsilon removed.
 *
 * This method is typically used in conjunction with the product of two sparse matrices
 * to automatically prune the smallest values as follows:
 * \code
 * C = (A*B).pruned();             // suppress numerical zeros (exact)
 * C = (A*B).pruned(ref);
 * C = (A*B).pruned(ref,epsilon);
 * \endcode
 * where \c ref is a meaningful non zero reference value.
 * */
template <typename Derived>
const SparseView<Derived> SparseMatrixBase<Derived>::pruned(const Scalar& reference, const RealScalar& epsilon) const {
  return SparseView<Derived>(derived(), reference, epsilon);
}

}  // end namespace Eigen

#endif