src/SparseCore/SparseDiagonalProduct.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_SPARSE_DIAGONAL_PRODUCT_H
0011 #define EIGEN_SPARSE_DIAGONAL_PRODUCT_H
0012
0013 namespace Eigen {
0014
0015 // The product of a diagonal matrix with a sparse matrix can be easily
0016 // implemented using expression template.
0017 // We have two consider very different cases:
0018 // 1 - diag * row-major sparse
0019 //     => each inner vector <=> scalar * sparse vector product
0020 //     => so we can reuse CwiseUnaryOp::InnerIterator
0021 // 2 - diag * col-major sparse
0022 //     => each inner vector <=> densevector * sparse vector cwise product
0023 //     => again, we can reuse specialization of CwiseBinaryOp::InnerIterator
0024 //        for that particular case
0025 // The two other cases are symmetric.
0026
0027 namespace internal {
0028
0029 enum {
0030   SDP_AsScalarProduct,
0031   SDP_AsCwiseProduct
0032 };
0033
0034 template<typename SparseXprType, typename DiagonalCoeffType, int SDP_Tag>
0035 struct sparse_diagonal_product_evaluator;
0036
0037 template<typename Lhs, typename Rhs, int ProductTag>
0038 struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, DiagonalShape, SparseShape>
0039   : public sparse_diagonal_product_evaluator<Rhs, typename Lhs::DiagonalVectorType, Rhs::Flags&RowMajorBit?SDP_AsScalarProduct:SDP_AsCwiseProduct>
0040 {
0041   typedef Product<Lhs, Rhs, DefaultProduct> XprType;
0042   enum { CoeffReadCost = HugeCost, Flags = Rhs::Flags&RowMajorBit, Alignment = 0 }; // FIXME CoeffReadCost & Flags
0043
0044   typedef sparse_diagonal_product_evaluator<Rhs, typename Lhs::DiagonalVectorType, Rhs::Flags&RowMajorBit?SDP_AsScalarProduct:SDP_AsCwiseProduct> Base;
0045   explicit product_evaluator(const XprType& xpr) : Base(xpr.rhs(), xpr.lhs().diagonal()) {}
0046 };
0047
0048 template<typename Lhs, typename Rhs, int ProductTag>
0049 struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, SparseShape, DiagonalShape>
0050   : public sparse_diagonal_product_evaluator<Lhs, Transpose<const typename Rhs::DiagonalVectorType>, Lhs::Flags&RowMajorBit?SDP_AsCwiseProduct:SDP_AsScalarProduct>
0051 {
0052   typedef Product<Lhs, Rhs, DefaultProduct> XprType;
0053   enum { CoeffReadCost = HugeCost, Flags = Lhs::Flags&RowMajorBit, Alignment = 0 }; // FIXME CoeffReadCost & Flags
0054
0055   typedef sparse_diagonal_product_evaluator<Lhs, Transpose<const typename Rhs::DiagonalVectorType>, Lhs::Flags&RowMajorBit?SDP_AsCwiseProduct:SDP_AsScalarProduct> Base;
0056   explicit product_evaluator(const XprType& xpr) : Base(xpr.lhs(), xpr.rhs().diagonal().transpose()) {}
0057 };
0058
0059 template<typename SparseXprType, typename DiagonalCoeffType>
0060 struct sparse_diagonal_product_evaluator<SparseXprType, DiagonalCoeffType, SDP_AsScalarProduct>
0061 {
0062 protected:
0063   typedef typename evaluator<SparseXprType>::InnerIterator SparseXprInnerIterator;
0064   typedef typename SparseXprType::Scalar Scalar;
0065
0066 public:
0067   class InnerIterator : public SparseXprInnerIterator
0068   {
0069   public:
0070     InnerIterator(const sparse_diagonal_product_evaluator &xprEval, Index outer)
0071       : SparseXprInnerIterator(xprEval.m_sparseXprImpl, outer),
0072         m_coeff(xprEval.m_diagCoeffImpl.coeff(outer))
0073     {}
0074
0075     EIGEN_STRONG_INLINE Scalar value() const { return m_coeff * SparseXprInnerIterator::value(); }
0076   protected:
0077     typename DiagonalCoeffType::Scalar m_coeff;
0078   };
0079
0080   sparse_diagonal_product_evaluator(const SparseXprType &sparseXpr, const DiagonalCoeffType &diagCoeff)
0081     : m_sparseXprImpl(sparseXpr), m_diagCoeffImpl(diagCoeff)
0082   {}
0083
0084   Index nonZerosEstimate() const { return m_sparseXprImpl.nonZerosEstimate(); }
0085
0086 protected:
0087   evaluator<SparseXprType> m_sparseXprImpl;
0088   evaluator<DiagonalCoeffType> m_diagCoeffImpl;
0089 };
0090
0091
0092 template<typename SparseXprType, typename DiagCoeffType>
0093 struct sparse_diagonal_product_evaluator<SparseXprType, DiagCoeffType, SDP_AsCwiseProduct>
0094 {
0095   typedef typename SparseXprType::Scalar Scalar;
0096   typedef typename SparseXprType::StorageIndex StorageIndex;
0097
0098   typedef typename nested_eval<DiagCoeffType,SparseXprType::IsRowMajor ? SparseXprType::RowsAtCompileTime
0099                                                                        : SparseXprType::ColsAtCompileTime>::type DiagCoeffNested;
0100
0101   class InnerIterator
0102   {
0103     typedef typename evaluator<SparseXprType>::InnerIterator SparseXprIter;
0104   public:
0105     InnerIterator(const sparse_diagonal_product_evaluator &xprEval, Index outer)
0106       : m_sparseIter(xprEval.m_sparseXprEval, outer), m_diagCoeffNested(xprEval.m_diagCoeffNested)
0107     {}
0108
0109     inline Scalar value() const { return m_sparseIter.value() * m_diagCoeffNested.coeff(index()); }
0110     inline StorageIndex index() const  { return m_sparseIter.index(); }
0111     inline Index outer() const  { return m_sparseIter.outer(); }
0112     inline Index col() const    { return SparseXprType::IsRowMajor ? m_sparseIter.index() : m_sparseIter.outer(); }
0113     inline Index row() const    { return SparseXprType::IsRowMajor ? m_sparseIter.outer() : m_sparseIter.index(); }
0114
0115     EIGEN_STRONG_INLINE InnerIterator& operator++() { ++m_sparseIter; return *this; }
0116     inline operator bool() const  { return m_sparseIter; }
0117
0118   protected:
0119     SparseXprIter m_sparseIter;
0120     DiagCoeffNested m_diagCoeffNested;
0121   };
0122
0123   sparse_diagonal_product_evaluator(const SparseXprType &sparseXpr, const DiagCoeffType &diagCoeff)
0124     : m_sparseXprEval(sparseXpr), m_diagCoeffNested(diagCoeff)
0125   {}
0126
0127   Index nonZerosEstimate() const { return m_sparseXprEval.nonZerosEstimate(); }
0128
0129 protected:
0130   evaluator<SparseXprType> m_sparseXprEval;
0131   DiagCoeffNested m_diagCoeffNested;
0132 };
0133
0134 } // end namespace internal
0135
0136 } // end namespace Eigen
0137
0138 #endif // EIGEN_SPARSE_DIAGONAL_PRODUCT_H