Eigen/src/SparseCore/SparseDenseProduct.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// 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_SPARSEDENSEPRODUCT_H
#define EIGEN_SPARSEDENSEPRODUCT_H

namespace Eigen { 

namespace internal {

template <> struct product_promote_storage_type<Sparse,Dense, OuterProduct> { typedef Sparse ret; };
template <> struct product_promote_storage_type<Dense,Sparse, OuterProduct> { typedef Sparse ret; };

template<typename SparseLhsType, typename DenseRhsType, typename DenseResType,
         typename AlphaType,
         int LhsStorageOrder = ((SparseLhsType::Flags&RowMajorBit)==RowMajorBit) ? RowMajor : ColMajor,
         bool ColPerCol = ((DenseRhsType::Flags&RowMajorBit)==0) || DenseRhsType::ColsAtCompileTime==1>
struct sparse_time_dense_product_impl;

template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, typename DenseResType::Scalar, RowMajor, true>
{
  typedef typename internal::remove_all<SparseLhsType>::type Lhs;
  typedef typename internal::remove_all<DenseRhsType>::type Rhs;
  typedef typename internal::remove_all<DenseResType>::type Res;
  typedef typename Lhs::Index Index;
  typedef typename evaluator<Lhs>::InnerIterator LhsInnerIterator;
  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
  {
    typename evaluator<Lhs>::type lhsEval(lhs);
    for(Index c=0; c<rhs.cols(); ++c)
    {
      Index n = lhs.outerSize();
      for(Index j=0; j<n; ++j)
      {
        typename Res::Scalar tmp(0);
        for(LhsInnerIterator it(lhsEval,j); it ;++it)
          tmp += it.value() * rhs.coeff(it.index(),c);
        res.coeffRef(j,c) = alpha * tmp;
      }
    }
  }
};

template<typename T1, typename T2/*, int _Options, typename _StrideType*/>
struct scalar_product_traits<T1, Ref<T2/*, _Options, _StrideType*/> >
{
  enum {
    Defined = 1
  };
  typedef typename CwiseUnaryOp<scalar_multiple2_op<T1, typename T2::Scalar>, T2>::PlainObject ReturnType;
};
template<typename SparseLhsType, typename DenseRhsType, typename DenseResType, typename AlphaType>
struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, AlphaType, ColMajor, true>
{
  typedef typename internal::remove_all<SparseLhsType>::type Lhs;
  typedef typename internal::remove_all<DenseRhsType>::type Rhs;
  typedef typename internal::remove_all<DenseResType>::type Res;
  typedef typename Lhs::Index Index;
  typedef typename evaluator<Lhs>::InnerIterator LhsInnerIterator;
  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha)
  {
    typename evaluator<Lhs>::type lhsEval(lhs);
    for(Index c=0; c<rhs.cols(); ++c)
    {
      for(Index j=0; j<lhs.outerSize(); ++j)
      {
//        typename Res::Scalar rhs_j = alpha * rhs.coeff(j,c);
        typename internal::scalar_product_traits<AlphaType, typename Rhs::Scalar>::ReturnType rhs_j(alpha * rhs.coeff(j,c));
        for(LhsInnerIterator it(lhsEval,j); it ;++it)
          res.coeffRef(it.index(),c) += it.value() * rhs_j;
      }
    }
  }
};

template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, typename DenseResType::Scalar, RowMajor, false>
{
  typedef typename internal::remove_all<SparseLhsType>::type Lhs;
  typedef typename internal::remove_all<DenseRhsType>::type Rhs;
  typedef typename internal::remove_all<DenseResType>::type Res;
  typedef typename Lhs::Index Index;
  typedef typename evaluator<Lhs>::InnerIterator LhsInnerIterator;
  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
  {
    typename evaluator<Lhs>::type lhsEval(lhs);
    for(Index j=0; j<lhs.outerSize(); ++j)
    {
      typename Res::RowXpr res_j(res.row(j));
      for(LhsInnerIterator it(lhsEval,j); it ;++it)
        res_j += (alpha*it.value()) * rhs.row(it.index());
    }
  }
};

template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, typename DenseResType::Scalar, ColMajor, false>
{
  typedef typename internal::remove_all<SparseLhsType>::type Lhs;
  typedef typename internal::remove_all<DenseRhsType>::type Rhs;
  typedef typename internal::remove_all<DenseResType>::type Res;
  typedef typename Lhs::Index Index;
  typedef typename evaluator<Lhs>::InnerIterator LhsInnerIterator;
  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
  {
    typename evaluator<Lhs>::type lhsEval(lhs);
    for(Index j=0; j<lhs.outerSize(); ++j)
    {
      typename Rhs::ConstRowXpr rhs_j(rhs.row(j));
      for(LhsInnerIterator it(lhsEval,j); it ;++it)
        res.row(it.index()) += (alpha*it.value()) * rhs_j;
    }
  }
};

template<typename SparseLhsType, typename DenseRhsType, typename DenseResType,typename AlphaType>
inline void sparse_time_dense_product(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha)
{
  sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, AlphaType>::run(lhs, rhs, res, alpha);
}

} // end namespace internal

namespace internal {

template<typename Lhs, typename Rhs, int ProductType>
struct generic_product_impl<Lhs, Rhs, SparseShape, DenseShape, ProductType>
{
  template<typename Dest>
  static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
  {
    typedef typename nested_eval<Lhs,Dynamic>::type LhsNested;
    typedef typename nested_eval<Rhs,Dynamic>::type RhsNested;
    LhsNested lhsNested(lhs);
    RhsNested rhsNested(rhs);
    
    dst.setZero();
    internal::sparse_time_dense_product(lhsNested, rhsNested, dst, typename Dest::Scalar(1));
  }
};

template<typename Lhs, typename Rhs, int ProductType>
struct generic_product_impl<Lhs, Rhs, SparseTriangularShape, DenseShape, ProductType>
  : generic_product_impl<Lhs, Rhs, SparseShape, DenseShape, ProductType>
{};

template<typename Lhs, typename Rhs, int ProductType>
struct generic_product_impl<Lhs, Rhs, DenseShape, SparseShape, ProductType>
{
  template<typename Dest>
  static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
  {
    typedef typename nested_eval<Lhs,Dynamic>::type LhsNested;
    typedef typename nested_eval<Rhs,Dynamic>::type RhsNested;
    LhsNested lhsNested(lhs);
    RhsNested rhsNested(rhs);
    
    dst.setZero();
    // transpose everything
    Transpose<Dest> dstT(dst);
    internal::sparse_time_dense_product(rhsNested.transpose(), lhsNested.transpose(), dstT, typename Dest::Scalar(1));
  }
};

template<typename Lhs, typename Rhs, int ProductType>
struct generic_product_impl<Lhs, Rhs, DenseShape, SparseTriangularShape, ProductType>
  : generic_product_impl<Lhs, Rhs, DenseShape, SparseShape, ProductType>
{};

template<typename LhsT, typename RhsT, bool NeedToTranspose>
struct sparse_dense_outer_product_evaluator
{
protected:
  typedef typename conditional<NeedToTranspose,RhsT,LhsT>::type Lhs1;
  typedef typename conditional<NeedToTranspose,LhsT,RhsT>::type ActualRhs;
  typedef Product<LhsT,RhsT,DefaultProduct> ProdXprType;
  
  // if the actual left-hand side is a dense vector,
  // then build a sparse-view so that we can seamlessly iterate over it.
  typedef typename conditional<is_same<typename internal::traits<Lhs1>::StorageKind,Sparse>::value,
            Lhs1, SparseView<Lhs1> >::type ActualLhs;
  typedef typename conditional<is_same<typename internal::traits<Lhs1>::StorageKind,Sparse>::value,
            Lhs1 const&, SparseView<Lhs1> >::type LhsArg;
            
  typedef typename evaluator<ActualLhs>::type LhsEval;
  typedef typename evaluator<ActualRhs>::type RhsEval;
  typedef typename evaluator<ActualLhs>::InnerIterator LhsIterator;
  typedef typename ProdXprType::Scalar Scalar;
  typedef typename ProdXprType::Index Index;
  
public:
  enum {
    Flags = NeedToTranspose ? RowMajorBit : 0,
    CoeffReadCost = Dynamic
  };
  
  class InnerIterator : public LhsIterator
  {
  public:
    InnerIterator(const sparse_dense_outer_product_evaluator &xprEval, Index outer)
      : LhsIterator(xprEval.m_lhsXprImpl, 0),
        m_outer(outer),
        m_empty(false),
        m_factor(get(xprEval.m_rhsXprImpl, outer, typename internal::traits<ActualRhs>::StorageKind() ))
    {}
    
    EIGEN_STRONG_INLINE Index outer() const { return m_outer; }
    EIGEN_STRONG_INLINE Index row()   const { return NeedToTranspose ? m_outer : LhsIterator::index(); }
    EIGEN_STRONG_INLINE Index col()   const { return NeedToTranspose ? LhsIterator::index() : m_outer; }

    EIGEN_STRONG_INLINE Scalar value() const { return LhsIterator::value() * m_factor; }
    EIGEN_STRONG_INLINE operator bool() const { return LhsIterator::operator bool() && (!m_empty); }
    
  protected:
    Scalar get(const RhsEval &rhs, Index outer, Dense = Dense()) const
    {
      return rhs.coeff(outer);
    }
    
    Scalar get(const RhsEval &rhs, Index outer, Sparse = Sparse())
    {
      typename RhsEval::InnerIterator it(rhs, outer);
      if (it && it.index()==0 && it.value()!=Scalar(0))
        return it.value();
      m_empty = true;
      return Scalar(0);
    }
    
    Index m_outer;
    bool m_empty;
    Scalar m_factor;
  };
  
  sparse_dense_outer_product_evaluator(const Lhs1 &lhs, const ActualRhs &rhs)
     : m_lhs(lhs), m_lhsXprImpl(m_lhs), m_rhsXprImpl(rhs)
  {}
  
  // transpose case
  sparse_dense_outer_product_evaluator(const ActualRhs &rhs, const Lhs1 &lhs)
     : m_lhs(lhs), m_lhsXprImpl(m_lhs), m_rhsXprImpl(rhs)
  {}
    
protected:
  const LhsArg m_lhs;
  typename evaluator<ActualLhs>::nestedType m_lhsXprImpl;
  typename evaluator<ActualRhs>::nestedType m_rhsXprImpl;
};

// sparse * dense outer product
template<typename Lhs, typename Rhs>
struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, OuterProduct, SparseShape, DenseShape, typename traits<Lhs>::Scalar, typename traits<Rhs>::Scalar> 
  : sparse_dense_outer_product_evaluator<Lhs,Rhs, Lhs::IsRowMajor>
{
  typedef sparse_dense_outer_product_evaluator<Lhs,Rhs, Lhs::IsRowMajor> Base;
  
  typedef Product<Lhs, Rhs> XprType;
  typedef typename XprType::PlainObject PlainObject;

  explicit product_evaluator(const XprType& xpr)
    : Base(xpr.lhs(), xpr.rhs())
  {}
  
};

template<typename Lhs, typename Rhs>
struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, OuterProduct, DenseShape, SparseShape, typename traits<Lhs>::Scalar, typename traits<Rhs>::Scalar> 
  : sparse_dense_outer_product_evaluator<Lhs,Rhs, Rhs::IsRowMajor>
{
  typedef sparse_dense_outer_product_evaluator<Lhs,Rhs, Rhs::IsRowMajor> Base;
  
  typedef Product<Lhs, Rhs> XprType;
  typedef typename XprType::PlainObject PlainObject;

  explicit product_evaluator(const XprType& xpr)
    : Base(xpr.lhs(), xpr.rhs())
  {}
  
};

} // end namespace internal

} // end namespace Eigen

#endif // EIGEN_SPARSEDENSEPRODUCT_H