1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2009-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_SPARSE_DIAGONAL_PRODUCT_H
11#define EIGEN_SPARSE_DIAGONAL_PRODUCT_H
12
13namespace Eigen {
14
15// The product of a diagonal matrix with a sparse matrix can be easily
16// implemented using expression template.
17// We have two consider very different cases:
18// 1 - diag * row-major sparse
19// => each inner vector <=> scalar * sparse vector product
20// => so we can reuse CwiseUnaryOp::InnerIterator
21// 2 - diag * col-major sparse
22// => each inner vector <=> densevector * sparse vector cwise product
23// => again, we can reuse specialization of CwiseBinaryOp::InnerIterator
24// for that particular case
25// The two other cases are symmetric.
26
27namespace internal {
28
29enum {
30 SDP_AsScalarProduct,
31 SDP_AsCwiseProduct
32};
33
34template<typename SparseXprType, typename DiagonalCoeffType, int SDP_Tag>
35struct sparse_diagonal_product_evaluator;
36
37template<typename Lhs, typename Rhs, int ProductTag>
38struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, DiagonalShape, SparseShape>
39 : public sparse_diagonal_product_evaluator<Rhs, typename Lhs::DiagonalVectorType, Rhs::Flags&RowMajorBit?SDP_AsScalarProduct:SDP_AsCwiseProduct>
40{
41 typedef Product<Lhs, Rhs, DefaultProduct> XprType;
42 enum { CoeffReadCost = HugeCost, Flags = Rhs::Flags&RowMajorBit, Alignment = 0 }; // FIXME CoeffReadCost & Flags
43
44 typedef sparse_diagonal_product_evaluator<Rhs, typename Lhs::DiagonalVectorType, Rhs::Flags&RowMajorBit?SDP_AsScalarProduct:SDP_AsCwiseProduct> Base;
45 explicit product_evaluator(const XprType& xpr) : Base(xpr.rhs(), xpr.lhs().diagonal()) {}
46};
47
48template<typename Lhs, typename Rhs, int ProductTag>
49struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, SparseShape, DiagonalShape>
50 : public sparse_diagonal_product_evaluator<Lhs, Transpose<const typename Rhs::DiagonalVectorType>, Lhs::Flags&RowMajorBit?SDP_AsCwiseProduct:SDP_AsScalarProduct>
51{
52 typedef Product<Lhs, Rhs, DefaultProduct> XprType;
53 enum { CoeffReadCost = HugeCost, Flags = Lhs::Flags&RowMajorBit, Alignment = 0 }; // FIXME CoeffReadCost & Flags
54
55 typedef sparse_diagonal_product_evaluator<Lhs, Transpose<const typename Rhs::DiagonalVectorType>, Lhs::Flags&RowMajorBit?SDP_AsCwiseProduct:SDP_AsScalarProduct> Base;
56 explicit product_evaluator(const XprType& xpr) : Base(xpr.lhs(), xpr.rhs().diagonal().transpose()) {}
57};
58
59template<typename SparseXprType, typename DiagonalCoeffType>
60struct sparse_diagonal_product_evaluator<SparseXprType, DiagonalCoeffType, SDP_AsScalarProduct>
61{
62protected:
63 typedef typename evaluator<SparseXprType>::InnerIterator SparseXprInnerIterator;
64 typedef typename SparseXprType::Scalar Scalar;
65
66public:
67 class InnerIterator : public SparseXprInnerIterator
68 {
69 public:
70 InnerIterator(const sparse_diagonal_product_evaluator &xprEval, Index outer)
71 : SparseXprInnerIterator(xprEval.m_sparseXprImpl, outer),
72 m_coeff(xprEval.m_diagCoeffImpl.coeff(outer))
73 {}
74
75 EIGEN_STRONG_INLINE Scalar value() const { return m_coeff * SparseXprInnerIterator::value(); }
76 protected:
77 typename DiagonalCoeffType::Scalar m_coeff;
78 };
79
80 sparse_diagonal_product_evaluator(const SparseXprType &sparseXpr, const DiagonalCoeffType &diagCoeff)
81 : m_sparseXprImpl(sparseXpr), m_diagCoeffImpl(diagCoeff)
82 {}
83
84 Index nonZerosEstimate() const { return m_sparseXprImpl.nonZerosEstimate(); }
85
86protected:
87 evaluator<SparseXprType> m_sparseXprImpl;
88 evaluator<DiagonalCoeffType> m_diagCoeffImpl;
89};
90
91
92template<typename SparseXprType, typename DiagCoeffType>
93struct sparse_diagonal_product_evaluator<SparseXprType, DiagCoeffType, SDP_AsCwiseProduct>
94{
95 typedef typename SparseXprType::Scalar Scalar;
96 typedef typename SparseXprType::StorageIndex StorageIndex;
97
98 typedef typename nested_eval<DiagCoeffType,SparseXprType::IsRowMajor ? SparseXprType::RowsAtCompileTime
99 : SparseXprType::ColsAtCompileTime>::type DiagCoeffNested;
100
101 class InnerIterator
102 {
103 typedef typename evaluator<SparseXprType>::InnerIterator SparseXprIter;
104 public:
105 InnerIterator(const sparse_diagonal_product_evaluator &xprEval, Index outer)
106 : m_sparseIter(xprEval.m_sparseXprEval, outer), m_diagCoeffNested(xprEval.m_diagCoeffNested)
107 {}
108
109 inline Scalar value() const { return m_sparseIter.value() * m_diagCoeffNested.coeff(index()); }
110 inline StorageIndex index() const { return m_sparseIter.index(); }
111 inline Index outer() const { return m_sparseIter.outer(); }
112 inline Index col() const { return SparseXprType::IsRowMajor ? m_sparseIter.index() : m_sparseIter.outer(); }
113 inline Index row() const { return SparseXprType::IsRowMajor ? m_sparseIter.outer() : m_sparseIter.index(); }
114
115 EIGEN_STRONG_INLINE InnerIterator& operator++() { ++m_sparseIter; return *this; }
116 inline operator bool() const { return m_sparseIter; }
117
118 protected:
119 SparseXprIter m_sparseIter;
120 DiagCoeffNested m_diagCoeffNested;
121 };
122
123 sparse_diagonal_product_evaluator(const SparseXprType &sparseXpr, const DiagCoeffType &diagCoeff)
124 : m_sparseXprEval(sparseXpr), m_diagCoeffNested(diagCoeff)
125 {}
126
127 Index nonZerosEstimate() const { return m_sparseXprEval.nonZerosEstimate(); }
128
129protected:
130 evaluator<SparseXprType> m_sparseXprEval;
131 DiagCoeffNested m_diagCoeffNested;
132};
133
134} // end namespace internal
135
136} // end namespace Eigen
137
138#endif // EIGEN_SPARSE_DIAGONAL_PRODUCT_H
139