Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / Sparse / SparseProduct.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra. Eigen itself is part of the KDE project.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24
25 #ifndef EIGEN_SPARSEPRODUCT_H
26 #define EIGEN_SPARSEPRODUCT_H
27
28 template<typename Lhs, typename Rhs> struct ei_sparse_product_mode
29 {
30   enum {
31
32     value = ((Lhs::Flags&Diagonal)==Diagonal || (Rhs::Flags&Diagonal)==Diagonal)
33           ? DiagonalProduct
34           :  (Rhs::Flags&Lhs::Flags&SparseBit)==SparseBit
35           ? SparseTimeSparseProduct
36           : (Lhs::Flags&SparseBit)==SparseBit
37           ? SparseTimeDenseProduct
38           : DenseTimeSparseProduct };
39 };
40
41 template<typename Lhs, typename Rhs, int ProductMode>
42 struct SparseProductReturnType
43 {
44   typedef const typename ei_nested<Lhs,Rhs::RowsAtCompileTime>::type LhsNested;
45   typedef const typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
46
47   typedef SparseProduct<LhsNested, RhsNested, ProductMode> Type;
48 };
49
50 template<typename Lhs, typename Rhs>
51 struct SparseProductReturnType<Lhs,Rhs,DiagonalProduct>
52 {
53   typedef const typename ei_nested<Lhs,Rhs::RowsAtCompileTime>::type LhsNested;
54   typedef const typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
55
56   typedef SparseDiagonalProduct<LhsNested, RhsNested> Type;
57 };
58
59 // sparse product return type specialization
60 template<typename Lhs, typename Rhs>
61 struct SparseProductReturnType<Lhs,Rhs,SparseTimeSparseProduct>
62 {
63   typedef typename ei_traits<Lhs>::Scalar Scalar;
64   enum {
65     LhsRowMajor = ei_traits<Lhs>::Flags & RowMajorBit,
66     RhsRowMajor = ei_traits<Rhs>::Flags & RowMajorBit,
67     TransposeRhs = (!LhsRowMajor) && RhsRowMajor,
68     TransposeLhs = LhsRowMajor && (!RhsRowMajor)
69   };
70
71   // FIXME if we transpose let's evaluate to a LinkedVectorMatrix since it is the
72   // type of the temporary to perform the transpose op
73   typedef typename ei_meta_if<TransposeLhs,
74     SparseMatrix<Scalar,0>,
75     const typename ei_nested<Lhs,Rhs::RowsAtCompileTime>::type>::ret LhsNested;
76
77   typedef typename ei_meta_if<TransposeRhs,
78     SparseMatrix<Scalar,0>,
79     const typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type>::ret RhsNested;
80
81   typedef SparseProduct<LhsNested, RhsNested, SparseTimeSparseProduct> Type;
82 };
83
84 template<typename LhsNested, typename RhsNested, int ProductMode>
85 struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >
86 {
87   // clean the nested types:
88   typedef typename ei_cleantype<LhsNested>::type _LhsNested;
89   typedef typename ei_cleantype<RhsNested>::type _RhsNested;
90   typedef typename _LhsNested::Scalar Scalar;
91
92   enum {
93     LhsCoeffReadCost = _LhsNested::CoeffReadCost,
94     RhsCoeffReadCost = _RhsNested::CoeffReadCost,
95     LhsFlags = _LhsNested::Flags,
96     RhsFlags = _RhsNested::Flags,
97
98     RowsAtCompileTime = _LhsNested::RowsAtCompileTime,
99     ColsAtCompileTime = _RhsNested::ColsAtCompileTime,
100     InnerSize = EIGEN_ENUM_MIN(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
101
102     MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
103     MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
104
105 //     LhsIsRowMajor = (LhsFlags & RowMajorBit)==RowMajorBit,
106 //     RhsIsRowMajor = (RhsFlags & RowMajorBit)==RowMajorBit,
107
108     EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit),
109     ResultIsSparse = ProductMode==SparseTimeSparseProduct || ProductMode==DiagonalProduct,
110
111     RemovedBits = ~( (EvalToRowMajor ? 0 : RowMajorBit) | (ResultIsSparse ? 0 : SparseBit) ),
112
113     Flags = (int(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
114           | EvalBeforeAssigningBit
115           | EvalBeforeNestingBit,
116
117     CoeffReadCost = Dynamic
118   };
119
120   typedef typename ei_meta_if<ResultIsSparse,
121     SparseMatrixBase<SparseProduct<LhsNested, RhsNested, ProductMode> >,
122     MatrixBase<SparseProduct<LhsNested, RhsNested, ProductMode> > >::ret Base;
123 };
124
125 template<typename LhsNested, typename RhsNested, int ProductMode>
126 class SparseProduct : ei_no_assignment_operator,
127   public ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >::Base
128 {
129   public:
130
131     EIGEN_GENERIC_PUBLIC_INTERFACE(SparseProduct)
132
133   private:
134
135     typedef typename ei_traits<SparseProduct>::_LhsNested _LhsNested;
136     typedef typename ei_traits<SparseProduct>::_RhsNested _RhsNested;
137
138   public:
139
140     template<typename Lhs, typename Rhs>
141     EIGEN_STRONG_INLINE SparseProduct(const Lhs& lhs, const Rhs& rhs)
142       : m_lhs(lhs), m_rhs(rhs)
143     {
144       ei_assert(lhs.cols() == rhs.rows());
145
146       enum {
147         ProductIsValid = _LhsNested::ColsAtCompileTime==Dynamic
148                       || _RhsNested::RowsAtCompileTime==Dynamic
149                       || int(_LhsNested::ColsAtCompileTime)==int(_RhsNested::RowsAtCompileTime),
150         AreVectors = _LhsNested::IsVectorAtCompileTime && _RhsNested::IsVectorAtCompileTime,
151         SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(_LhsNested,_RhsNested)
152       };
153       // note to the lost user:
154       //    * for a dot product use: v1.dot(v2)
155       //    * for a coeff-wise product use: v1.cwise()*v2
156       EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
157         INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
158       EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
159         INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
160       EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
161     }
162
163     EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); }
164     EIGEN_STRONG_INLINE int cols() const { return m_rhs.cols(); }
165
166     EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
167     EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
168
169   protected:
170     LhsNested m_lhs;
171     RhsNested m_rhs;
172 };
173
174 // perform a pseudo in-place sparse * sparse product assuming all matrices are col major
175 template<typename Lhs, typename Rhs, typename ResultType>
176 static void ei_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
177 {
178   typedef typename ei_traits<typename ei_cleantype<Lhs>::type>::Scalar Scalar;
179
180   // make sure to call innerSize/outerSize since we fake the storage order.
181   int rows = lhs.innerSize();
182   int cols = rhs.outerSize();
183   //int size = lhs.outerSize();
184   ei_assert(lhs.outerSize() == rhs.innerSize());
185
186   // allocate a temporary buffer
187   AmbiVector<Scalar> tempVector(rows);
188
189   // estimate the number of non zero entries
190   float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols()));
191   float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
192   float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f);
193
194   res.resize(rows, cols);
195   res.startFill(int(ratioRes*rows*cols));
196   for (int j=0; j<cols; ++j)
197   {
198     // let's do a more accurate determination of the nnz ratio for the current column j of res
199     //float ratioColRes = std::min(ratioLhs * rhs.innerNonZeros(j), 1.f);
200     // FIXME find a nice way to get the number of nonzeros of a sub matrix (here an inner vector)
201     float ratioColRes = ratioRes;
202     tempVector.init(ratioColRes);
203     tempVector.setZero();
204     for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
205     {
206       // FIXME should be written like this: tmp += rhsIt.value() * lhs.col(rhsIt.index())
207       tempVector.restart();
208       Scalar x = rhsIt.value();
209       for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
210       {
211         tempVector.coeffRef(lhsIt.index()) += lhsIt.value() * x;
212       }
213     }
214     for (typename AmbiVector<Scalar>::Iterator it(tempVector); it; ++it)
215       if (ResultType::Flags&RowMajorBit)
216         res.fill(j,it.index()) = it.value();
217       else
218         res.fill(it.index(), j) = it.value();
219   }
220   res.endFill();
221 }
222
223 template<typename Lhs, typename Rhs, typename ResultType,
224   int LhsStorageOrder = ei_traits<Lhs>::Flags&RowMajorBit,
225   int RhsStorageOrder = ei_traits<Rhs>::Flags&RowMajorBit,
226   int ResStorageOrder = ei_traits<ResultType>::Flags&RowMajorBit>
227 struct ei_sparse_product_selector;
228
229 template<typename Lhs, typename Rhs, typename ResultType>
230 struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
231 {
232   typedef typename ei_traits<typename ei_cleantype<Lhs>::type>::Scalar Scalar;
233
234   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
235   {
236     typename ei_cleantype<ResultType>::type _res(res.rows(), res.cols());
237     ei_sparse_product_impl<Lhs,Rhs,ResultType>(lhs, rhs, _res);
238     res.swap(_res);
239   }
240 };
241
242 template<typename Lhs, typename Rhs, typename ResultType>
243 struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
244 {
245   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
246   {
247     // we need a col-major matrix to hold the result
248     typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
249     SparseTemporaryType _res(res.rows(), res.cols());
250     ei_sparse_product_impl<Lhs,Rhs,SparseTemporaryType>(lhs, rhs, _res);
251     res = _res;
252   }
253 };
254
255 template<typename Lhs, typename Rhs, typename ResultType>
256 struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
257 {
258   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
259   {
260     // let's transpose the product to get a column x column product
261     typename ei_cleantype<ResultType>::type _res(res.rows(), res.cols());
262     ei_sparse_product_impl<Rhs,Lhs,ResultType>(rhs, lhs, _res);
263     res.swap(_res);
264   }
265 };
266
267 template<typename Lhs, typename Rhs, typename ResultType>
268 struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
269 {
270   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
271   {
272     // let's transpose the product to get a column x column product
273     typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
274     SparseTemporaryType _res(res.cols(), res.rows());
275     ei_sparse_product_impl<Rhs,Lhs,SparseTemporaryType>(rhs, lhs, _res);
276     res = _res.transpose();
277   }
278 };
279
280 // NOTE eventually let's transpose one argument even in this case since it might be expensive if
281 // the result is not dense.
282 // template<typename Lhs, typename Rhs, typename ResultType, int ResStorageOrder>
283 // struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,ResStorageOrder>
284 // {
285 //   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
286 //   {
287 //     // trivial product as lhs.row/rhs.col dot products
288 //     // loop over the preferred order of the result
289 //   }
290 // };
291
292 // NOTE the 2 others cases (col row *) must never occurs since they are caught
293 // by ProductReturnType which transform it to (col col *) by evaluating rhs.
294
295
296 // template<typename Derived>
297 // template<typename Lhs, typename Rhs>
298 // inline Derived& SparseMatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs>& product)
299 // {
300 // //   std::cout << "sparse product to dense\n";
301 //   ei_sparse_product_selector<
302 //     typename ei_cleantype<Lhs>::type,
303 //     typename ei_cleantype<Rhs>::type,
304 //     typename ei_cleantype<Derived>::type>::run(product.lhs(),product.rhs(),derived());
305 //   return derived();
306 // }
307
308 // sparse = sparse * sparse
309 template<typename Derived>
310 template<typename Lhs, typename Rhs>
311 inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs,SparseTimeSparseProduct>& product)
312 {
313   ei_sparse_product_selector<
314     typename ei_cleantype<Lhs>::type,
315     typename ei_cleantype<Rhs>::type,
316     Derived>::run(product.lhs(),product.rhs(),derived());
317   return derived();
318 }
319
320 // dense = sparse * dense
321 // template<typename Derived>
322 // template<typename Lhs, typename Rhs>
323 // Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeDenseProduct>& product)
324 // {
325 //   typedef typename ei_cleantype<Lhs>::type _Lhs;
326 //   typedef typename _Lhs::InnerIterator LhsInnerIterator;
327 //   enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit };
328 //   derived().setZero();
329 //   for (int j=0; j<product.lhs().outerSize(); ++j)
330 //     for (LhsInnerIterator i(product.lhs(),j); i; ++i)
331 //       derived().row(LhsIsRowMajor ? j : i.index()) += i.value() * product.rhs().row(LhsIsRowMajor ? i.index() : j);
332 //   return derived();
333 // }
334
335 template<typename Derived>
336 template<typename Lhs, typename Rhs>
337 Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeDenseProduct>& product)
338 {
339   typedef typename ei_cleantype<Lhs>::type _Lhs;
340   typedef typename ei_cleantype<Rhs>::type _Rhs;
341   typedef typename _Lhs::InnerIterator LhsInnerIterator;
342   enum {
343     LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit,
344     LhsIsSelfAdjoint = (_Lhs::Flags&SelfAdjointBit)==SelfAdjointBit,
345     ProcessFirstHalf = LhsIsSelfAdjoint
346       && (   ((_Lhs::Flags&(UpperTriangularBit|LowerTriangularBit))==0)
347           || ( (_Lhs::Flags&UpperTriangularBit) && !LhsIsRowMajor)
348           || ( (_Lhs::Flags&LowerTriangularBit) && LhsIsRowMajor) ),
349     ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf)
350   };
351   derived().setZero();
352   for (int j=0; j<product.lhs().outerSize(); ++j)
353   {
354     LhsInnerIterator i(product.lhs(),j);
355     if (ProcessSecondHalf && i && (i.index()==j))
356     {
357       derived().row(j) += i.value() * product.rhs().row(j);
358       ++i;
359     }
360     Block<Derived,1,Derived::ColsAtCompileTime> res(derived().row(LhsIsRowMajor ? j : 0));
361     for (; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
362     {
363       if (LhsIsSelfAdjoint)
364       {
365         int a = LhsIsRowMajor ? j : i.index();
366         int b = LhsIsRowMajor ? i.index() : j;
367         Scalar v = i.value();
368         derived().row(a) += (v) * product.rhs().row(b);
369         derived().row(b) += ei_conj(v) * product.rhs().row(a);
370       }
371       else if (LhsIsRowMajor)
372         res += i.value() * product.rhs().row(i.index());
373       else
374         derived().row(i.index()) += i.value() * product.rhs().row(j);
375     }
376     if (ProcessFirstHalf && i && (i.index()==j))
377       derived().row(j) += i.value() * product.rhs().row(j);
378   }
379   return derived();
380 }
381
382 // dense = dense * sparse
383 template<typename Derived>
384 template<typename Lhs, typename Rhs>
385 Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,DenseTimeSparseProduct>& product)
386 {
387   typedef typename ei_cleantype<Rhs>::type _Rhs;
388   typedef typename _Rhs::InnerIterator RhsInnerIterator;
389   enum { RhsIsRowMajor = (_Rhs::Flags&RowMajorBit)==RowMajorBit };
390   derived().setZero();
391   for (int j=0; j<product.rhs().outerSize(); ++j)
392     for (RhsInnerIterator i(product.rhs(),j); i; ++i)
393       derived().col(RhsIsRowMajor ? i.index() : j) += i.value() * product.lhs().col(RhsIsRowMajor ? j : i.index());
394   return derived();
395 }
396
397 // sparse * sparse
398 template<typename Derived>
399 template<typename OtherDerived>
400 EIGEN_STRONG_INLINE const typename SparseProductReturnType<Derived,OtherDerived>::Type
401 SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other) const
402 {
403   return typename SparseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
404 }
405
406 // sparse * dense
407 template<typename Derived>
408 template<typename OtherDerived>
409 EIGEN_STRONG_INLINE const typename SparseProductReturnType<Derived,OtherDerived>::Type
410 SparseMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
411 {
412   return typename SparseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
413 }
414
415 #endif // EIGEN_SPARSEPRODUCT_H