1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2011 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_CONSERVATIVESPARSESPARSEPRODUCT_H
11 #define EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
12 
13 namespace Eigen {
14 
15 namespace internal {
16 
17 template<typename Lhs, typename Rhs, typename ResultType>
conservative_sparse_sparse_product_impl(const Lhs & lhs,const Rhs & rhs,ResultType & res)18 static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
19 {
20   typedef typename remove_all<Lhs>::type::Scalar Scalar;
21   typedef typename remove_all<Lhs>::type::Index Index;
22 
23   // make sure to call innerSize/outerSize since we fake the storage order.
24   Index rows = lhs.innerSize();
25   Index cols = rhs.outerSize();
26   eigen_assert(lhs.outerSize() == rhs.innerSize());
27 
28   std::vector<bool> mask(rows,false);
29   Matrix<Scalar,Dynamic,1> values(rows);
30   Matrix<Index,Dynamic,1>  indices(rows);
31 
32   // estimate the number of non zero entries
33   // given a rhs column containing Y non zeros, we assume that the respective Y columns
34   // of the lhs differs in average of one non zeros, thus the number of non zeros for
35   // the product of a rhs column with the lhs is X+Y where X is the average number of non zero
36   // per column of the lhs.
37   // Therefore, we have nnz(lhs*rhs) = nnz(lhs) + nnz(rhs)
38   Index estimated_nnz_prod = lhs.nonZeros() + rhs.nonZeros();
39 
40   res.setZero();
41   res.reserve(Index(estimated_nnz_prod));
42   // we compute each column of the result, one after the other
43   for (Index j=0; j<cols; ++j)
44   {
45 
46     res.startVec(j);
47     Index nnz = 0;
48     for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
49     {
50       Scalar y = rhsIt.value();
51       Index k = rhsIt.index();
52       for (typename Lhs::InnerIterator lhsIt(lhs, k); lhsIt; ++lhsIt)
53       {
54         Index i = lhsIt.index();
55         Scalar x = lhsIt.value();
56         if(!mask[i])
57         {
58           mask[i] = true;
59           values[i] = x * y;
60           indices[nnz] = i;
61           ++nnz;
62         }
63         else
64           values[i] += x * y;
65       }
66     }
67 
68     // unordered insertion
69     for(Index k=0; k<nnz; ++k)
70     {
71       Index i = indices[k];
72       res.insertBackByOuterInnerUnordered(j,i) = values[i];
73       mask[i] = false;
74     }
75 
76 #if 0
77     // alternative ordered insertion code:
78 
79     Index t200 = rows/(log2(200)*1.39);
80     Index t = (rows*100)/139;
81 
82     // FIXME reserve nnz non zeros
83     // FIXME implement fast sort algorithms for very small nnz
84     // if the result is sparse enough => use a quick sort
85     // otherwise => loop through the entire vector
86     // In order to avoid to perform an expensive log2 when the
87     // result is clearly very sparse we use a linear bound up to 200.
88     //if((nnz<200 && nnz<t200) || nnz * log2(nnz) < t)
89     //res.startVec(j);
90     if(true)
91     {
92       if(nnz>1) std::sort(indices.data(),indices.data()+nnz);
93       for(Index k=0; k<nnz; ++k)
94       {
95         Index i = indices[k];
96         res.insertBackByOuterInner(j,i) = values[i];
97         mask[i] = false;
98       }
99     }
100     else
101     {
102       // dense path
103       for(Index i=0; i<rows; ++i)
104       {
105         if(mask[i])
106         {
107           mask[i] = false;
108           res.insertBackByOuterInner(j,i) = values[i];
109         }
110       }
111     }
112 #endif
113 
114   }
115   res.finalize();
116 }
117 
118 
119 } // end namespace internal
120 
121 namespace internal {
122 
123 template<typename Lhs, typename Rhs, typename ResultType,
124   int LhsStorageOrder = (traits<Lhs>::Flags&RowMajorBit) ? RowMajor : ColMajor,
125   int RhsStorageOrder = (traits<Rhs>::Flags&RowMajorBit) ? RowMajor : ColMajor,
126   int ResStorageOrder = (traits<ResultType>::Flags&RowMajorBit) ? RowMajor : ColMajor>
127 struct conservative_sparse_sparse_product_selector;
128 
129 template<typename Lhs, typename Rhs, typename ResultType>
130 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
131 {
132   typedef typename remove_all<Lhs>::type LhsCleaned;
133   typedef typename LhsCleaned::Scalar Scalar;
134 
135   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
136   {
137     typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::Index> RowMajorMatrix;
138     typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> ColMajorMatrix;
139     ColMajorMatrix resCol(lhs.rows(),rhs.cols());
140     internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
141     // sort the non zeros:
142     RowMajorMatrix resRow(resCol);
143     res = resRow;
144   }
145 };
146 
147 template<typename Lhs, typename Rhs, typename ResultType>
148 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,ColMajor>
149 {
150   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
151   {
152      typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::Index> RowMajorMatrix;
153      RowMajorMatrix rhsRow = rhs;
154      RowMajorMatrix resRow(lhs.rows(), rhs.cols());
155      internal::conservative_sparse_sparse_product_impl<RowMajorMatrix,Lhs,RowMajorMatrix>(rhsRow, lhs, resRow);
156      res = resRow;
157   }
158 };
159 
160 template<typename Lhs, typename Rhs, typename ResultType>
161 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,RowMajor,ColMajor>
162 {
163   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
164   {
165     typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::Index> RowMajorMatrix;
166     RowMajorMatrix lhsRow = lhs;
167     RowMajorMatrix resRow(lhs.rows(), rhs.cols());
168     internal::conservative_sparse_sparse_product_impl<Rhs,RowMajorMatrix,RowMajorMatrix>(rhs, lhsRow, resRow);
169     res = resRow;
170   }
171 };
172 
173 template<typename Lhs, typename Rhs, typename ResultType>
174 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
175 {
176   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
177   {
178     typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::Index> RowMajorMatrix;
179     RowMajorMatrix resRow(lhs.rows(), rhs.cols());
180     internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
181     res = resRow;
182   }
183 };
184 
185 
186 template<typename Lhs, typename Rhs, typename ResultType>
187 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
188 {
189   typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
190 
191   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
192   {
193     typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> ColMajorMatrix;
194     ColMajorMatrix resCol(lhs.rows(), rhs.cols());
195     internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
196     res = resCol;
197   }
198 };
199 
200 template<typename Lhs, typename Rhs, typename ResultType>
201 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,RowMajor>
202 {
203   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
204   {
205     typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> ColMajorMatrix;
206     ColMajorMatrix lhsCol = lhs;
207     ColMajorMatrix resCol(lhs.rows(), rhs.cols());
208     internal::conservative_sparse_sparse_product_impl<ColMajorMatrix,Rhs,ColMajorMatrix>(lhsCol, rhs, resCol);
209     res = resCol;
210   }
211 };
212 
213 template<typename Lhs, typename Rhs, typename ResultType>
214 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,RowMajor,RowMajor>
215 {
216   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
217   {
218     typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> ColMajorMatrix;
219     ColMajorMatrix rhsCol = rhs;
220     ColMajorMatrix resCol(lhs.rows(), rhs.cols());
221     internal::conservative_sparse_sparse_product_impl<Lhs,ColMajorMatrix,ColMajorMatrix>(lhs, rhsCol, resCol);
222     res = resCol;
223   }
224 };
225 
226 template<typename Lhs, typename Rhs, typename ResultType>
227 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
228 {
229   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
230   {
231     typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::Index> RowMajorMatrix;
232     typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> ColMajorMatrix;
233     RowMajorMatrix resRow(lhs.rows(),rhs.cols());
234     internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
235     // sort the non zeros:
236     ColMajorMatrix resCol(resRow);
237     res = resCol;
238   }
239 };
240 
241 } // end namespace internal
242 
243 } // end namespace Eigen
244 
245 #endif // EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
246