1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008-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_SPARSEPRODUCT_H
11#define EIGEN_SPARSEPRODUCT_H
12
13namespace Eigen {
14
15/** \returns an expression of the product of two sparse matrices.
16 * By default a conservative product preserving the symbolic non zeros is performed.
17 * The automatic pruning of the small values can be achieved by calling the pruned() function
18 * in which case a totally different product algorithm is employed:
19 * \code
20 * C = (A*B).pruned(); // supress numerical zeros (exact)
21 * C = (A*B).pruned(ref);
22 * C = (A*B).pruned(ref,epsilon);
23 * \endcode
24 * where \c ref is a meaningful non zero reference value.
25 * */
26template<typename Derived>
27template<typename OtherDerived>
28inline const Product<Derived,OtherDerived,AliasFreeProduct>
29SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other) const
30{
31 return Product<Derived,OtherDerived,AliasFreeProduct>(derived(), other.derived());
32}
33
34namespace internal {
35
36// sparse * sparse
37template<typename Lhs, typename Rhs, int ProductType>
38struct generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType>
39{
40 template<typename Dest>
41 static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
42 {
43 evalTo(dst, lhs, rhs, typename evaluator_traits<Dest>::Shape());
44 }
45
46 // dense += sparse * sparse
47 template<typename Dest,typename ActualLhs>
48 static void addTo(Dest& dst, const ActualLhs& lhs, const Rhs& rhs, typename enable_if<is_same<typename evaluator_traits<Dest>::Shape,DenseShape>::value,int*>::type* = 0)
49 {
50 typedef typename nested_eval<ActualLhs,Dynamic>::type LhsNested;
51 typedef typename nested_eval<Rhs,Dynamic>::type RhsNested;
52 LhsNested lhsNested(lhs);
53 RhsNested rhsNested(rhs);
54 internal::sparse_sparse_to_dense_product_selector<typename remove_all<LhsNested>::type,
55 typename remove_all<RhsNested>::type, Dest>::run(lhsNested,rhsNested,dst);
56 }
57
58 // dense -= sparse * sparse
59 template<typename Dest>
60 static void subTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, typename enable_if<is_same<typename evaluator_traits<Dest>::Shape,DenseShape>::value,int*>::type* = 0)
61 {
62 addTo(dst, -lhs, rhs);
63 }
64
65protected:
66
67 // sparse = sparse * sparse
68 template<typename Dest>
69 static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, SparseShape)
70 {
71 typedef typename nested_eval<Lhs,Dynamic>::type LhsNested;
72 typedef typename nested_eval<Rhs,Dynamic>::type RhsNested;
73 LhsNested lhsNested(lhs);
74 RhsNested rhsNested(rhs);
75 internal::conservative_sparse_sparse_product_selector<typename remove_all<LhsNested>::type,
76 typename remove_all<RhsNested>::type, Dest>::run(lhsNested,rhsNested,dst);
77 }
78
79 // dense = sparse * sparse
80 template<typename Dest>
81 static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, DenseShape)
82 {
83 dst.setZero();
84 addTo(dst, lhs, rhs);
85 }
86};
87
88// sparse * sparse-triangular
89template<typename Lhs, typename Rhs, int ProductType>
90struct generic_product_impl<Lhs, Rhs, SparseShape, SparseTriangularShape, ProductType>
91 : public generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType>
92{};
93
94// sparse-triangular * sparse
95template<typename Lhs, typename Rhs, int ProductType>
96struct generic_product_impl<Lhs, Rhs, SparseTriangularShape, SparseShape, ProductType>
97 : public generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType>
98{};
99
100// dense = sparse-product (can be sparse*sparse, sparse*perm, etc.)
101template< typename DstXprType, typename Lhs, typename Rhs>
102struct Assignment<DstXprType, Product<Lhs,Rhs,AliasFreeProduct>, internal::assign_op<typename DstXprType::Scalar,typename Product<Lhs,Rhs,AliasFreeProduct>::Scalar>, Sparse2Dense>
103{
104 typedef Product<Lhs,Rhs,AliasFreeProduct> SrcXprType;
105 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &)
106 {
107 Index dstRows = src.rows();
108 Index dstCols = src.cols();
109 if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
110 dst.resize(dstRows, dstCols);
111
112 generic_product_impl<Lhs, Rhs>::evalTo(dst,src.lhs(),src.rhs());
113 }
114};
115
116// dense += sparse-product (can be sparse*sparse, sparse*perm, etc.)
117template< typename DstXprType, typename Lhs, typename Rhs>
118struct Assignment<DstXprType, Product<Lhs,Rhs,AliasFreeProduct>, internal::add_assign_op<typename DstXprType::Scalar,typename Product<Lhs,Rhs,AliasFreeProduct>::Scalar>, Sparse2Dense>
119{
120 typedef Product<Lhs,Rhs,AliasFreeProduct> SrcXprType;
121 static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &)
122 {
123 generic_product_impl<Lhs, Rhs>::addTo(dst,src.lhs(),src.rhs());
124 }
125};
126
127// dense -= sparse-product (can be sparse*sparse, sparse*perm, etc.)
128template< typename DstXprType, typename Lhs, typename Rhs>
129struct Assignment<DstXprType, Product<Lhs,Rhs,AliasFreeProduct>, internal::sub_assign_op<typename DstXprType::Scalar,typename Product<Lhs,Rhs,AliasFreeProduct>::Scalar>, Sparse2Dense>
130{
131 typedef Product<Lhs,Rhs,AliasFreeProduct> SrcXprType;
132 static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &)
133 {
134 generic_product_impl<Lhs, Rhs>::subTo(dst,src.lhs(),src.rhs());
135 }
136};
137
138template<typename Lhs, typename Rhs, int Options>
139struct unary_evaluator<SparseView<Product<Lhs, Rhs, Options> >, IteratorBased>
140 : public evaluator<typename Product<Lhs, Rhs, DefaultProduct>::PlainObject>
141{
142 typedef SparseView<Product<Lhs, Rhs, Options> > XprType;
143 typedef typename XprType::PlainObject PlainObject;
144 typedef evaluator<PlainObject> Base;
145
146 explicit unary_evaluator(const XprType& xpr)
147 : m_result(xpr.rows(), xpr.cols())
148 {
149 using std::abs;
150 ::new (static_cast<Base*>(this)) Base(m_result);
151 typedef typename nested_eval<Lhs,Dynamic>::type LhsNested;
152 typedef typename nested_eval<Rhs,Dynamic>::type RhsNested;
153 LhsNested lhsNested(xpr.nestedExpression().lhs());
154 RhsNested rhsNested(xpr.nestedExpression().rhs());
155
156 internal::sparse_sparse_product_with_pruning_selector<typename remove_all<LhsNested>::type,
157 typename remove_all<RhsNested>::type, PlainObject>::run(lhsNested,rhsNested,m_result,
158 abs(xpr.reference())*xpr.epsilon());
159 }
160
161protected:
162 PlainObject m_result;
163};
164
165} // end namespace internal
166
167} // end namespace Eigen
168
169#endif // EIGEN_SPARSEPRODUCT_H
170