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