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 
13 namespace 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();             // suppress 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   * */
26 template<typename Derived>
27 template<typename OtherDerived>
28 inline const Product<Derived,OtherDerived,AliasFreeProduct>
29 SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other) const
30 {
31   return Product<Derived,OtherDerived,AliasFreeProduct>(derived(), other.derived());
32 }
33 
34 namespace internal {
35 
36 // sparse * sparse
37 template<typename Lhs, typename Rhs, int ProductType>
38 struct 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 
65 protected:
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
89 template<typename Lhs, typename Rhs, int ProductType>
90 struct generic_product_impl<Lhs, Rhs, SparseShape, SparseTriangularShape, ProductType>
91  : public generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType>
92 {};
93 
94 // sparse-triangular * sparse
95 template<typename Lhs, typename Rhs, int ProductType>
96 struct 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.)
101 template< typename DstXprType, typename Lhs, typename Rhs>
102 struct 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.)
117 template< typename DstXprType, typename Lhs, typename Rhs>
118 struct 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.)
128 template< typename DstXprType, typename Lhs, typename Rhs>
129 struct 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 
138 template<typename Lhs, typename Rhs, int Options>
139 struct 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 
161 protected:
162   PlainObject m_result;
163 };
164 
165 } // end namespace internal
166 
167 // sparse matrix = sparse-product (can be sparse*sparse, sparse*perm, etc.)
168 template<typename Scalar, int _Options, typename _StorageIndex>
169 template<typename Lhs, typename Rhs>
170 SparseMatrix<Scalar,_Options,_StorageIndex>& SparseMatrix<Scalar,_Options,_StorageIndex>::operator=(const Product<Lhs,Rhs,AliasFreeProduct>& src)
171 {
172   // std::cout << "in Assignment : " << DstOptions << "\n";
173   SparseMatrix dst(src.rows(),src.cols());
174   internal::generic_product_impl<Lhs, Rhs>::evalTo(dst,src.lhs(),src.rhs());
175   this->swap(dst);
176   return *this;
177 }
178 
179 } // end namespace Eigen
180 
181 #endif // EIGEN_SPARSEPRODUCT_H
182