src/SparseCore/ConservativeSparseSparseProduct.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
0011 #define EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
0012
0013 namespace Eigen {
0014
0015 namespace internal {
0016
0017 template<typename Lhs, typename Rhs, typename ResultType>
0018 static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res, bool sortedInsertion = false)
0019 {
0020   typedef typename remove_all<Lhs>::type::Scalar LhsScalar;
0021   typedef typename remove_all<Rhs>::type::Scalar RhsScalar;
0022   typedef typename remove_all<ResultType>::type::Scalar ResScalar;
0023
0024   // make sure to call innerSize/outerSize since we fake the storage order.
0025   Index rows = lhs.innerSize();
0026   Index cols = rhs.outerSize();
0027   eigen_assert(lhs.outerSize() == rhs.innerSize());
0028
0029   ei_declare_aligned_stack_constructed_variable(bool,   mask,     rows, 0);
0030   ei_declare_aligned_stack_constructed_variable(ResScalar, values,   rows, 0);
0031   ei_declare_aligned_stack_constructed_variable(Index,  indices,  rows, 0);
0032
0033   std::memset(mask,0,sizeof(bool)*rows);
0034
0035   evaluator<Lhs> lhsEval(lhs);
0036   evaluator<Rhs> rhsEval(rhs);
0037
0038   // estimate the number of non zero entries
0039   // given a rhs column containing Y non zeros, we assume that the respective Y columns
0040   // of the lhs differs in average of one non zeros, thus the number of non zeros for
0041   // the product of a rhs column with the lhs is X+Y where X is the average number of non zero
0042   // per column of the lhs.
0043   // Therefore, we have nnz(lhs*rhs) = nnz(lhs) + nnz(rhs)
0044   Index estimated_nnz_prod = lhsEval.nonZerosEstimate() + rhsEval.nonZerosEstimate();
0045
0046   res.setZero();
0047   res.reserve(Index(estimated_nnz_prod));
0048   // we compute each column of the result, one after the other
0049   for (Index j=0; j<cols; ++j)
0050   {
0051
0052     res.startVec(j);
0053     Index nnz = 0;
0054     for (typename evaluator<Rhs>::InnerIterator rhsIt(rhsEval, j); rhsIt; ++rhsIt)
0055     {
0056       RhsScalar y = rhsIt.value();
0057       Index k = rhsIt.index();
0058       for (typename evaluator<Lhs>::InnerIterator lhsIt(lhsEval, k); lhsIt; ++lhsIt)
0059       {
0060         Index i = lhsIt.index();
0061         LhsScalar x = lhsIt.value();
0062         if(!mask[i])
0063         {
0064           mask[i] = true;
0065           values[i] = x * y;
0066           indices[nnz] = i;
0067           ++nnz;
0068         }
0069         else
0070           values[i] += x * y;
0071       }
0072     }
0073     if(!sortedInsertion)
0074     {
0075       // unordered insertion
0076       for(Index k=0; k<nnz; ++k)
0077       {
0078         Index i = indices[k];
0079         res.insertBackByOuterInnerUnordered(j,i) = values[i];
0080         mask[i] = false;
0081       }
0082     }
0083     else
0084     {
0085       // alternative ordered insertion code:
0086       const Index t200 = rows/11; // 11 == (log2(200)*1.39)
0087       const Index t = (rows*100)/139;
0088
0089       // FIXME reserve nnz non zeros
0090       // FIXME implement faster sorting algorithms for very small nnz
0091       // if the result is sparse enough => use a quick sort
0092       // otherwise => loop through the entire vector
0093       // In order to avoid to perform an expensive log2 when the
0094       // result is clearly very sparse we use a linear bound up to 200.
0095       if((nnz<200 && nnz<t200) || nnz * numext::log2(int(nnz)) < t)
0096       {
0097         if(nnz>1) std::sort(indices,indices+nnz);
0098         for(Index k=0; k<nnz; ++k)
0099         {
0100           Index i = indices[k];
0101           res.insertBackByOuterInner(j,i) = values[i];
0102           mask[i] = false;
0103         }
0104       }
0105       else
0106       {
0107         // dense path
0108         for(Index i=0; i<rows; ++i)
0109         {
0110           if(mask[i])
0111           {
0112             mask[i] = false;
0113             res.insertBackByOuterInner(j,i) = values[i];
0114           }
0115         }
0116       }
0117     }
0118   }
0119   res.finalize();
0120 }
0121
0122
0123 } // end namespace internal
0124
0125 namespace internal {
0126
0127 template<typename Lhs, typename Rhs, typename ResultType,
0128   int LhsStorageOrder = (traits<Lhs>::Flags&RowMajorBit) ? RowMajor : ColMajor,
0129   int RhsStorageOrder = (traits<Rhs>::Flags&RowMajorBit) ? RowMajor : ColMajor,
0130   int ResStorageOrder = (traits<ResultType>::Flags&RowMajorBit) ? RowMajor : ColMajor>
0131 struct conservative_sparse_sparse_product_selector;
0132
0133 template<typename Lhs, typename Rhs, typename ResultType>
0134 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
0135 {
0136   typedef typename remove_all<Lhs>::type LhsCleaned;
0137   typedef typename LhsCleaned::Scalar Scalar;
0138
0139   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
0140   {
0141     typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::StorageIndex> RowMajorMatrix;
0142     typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::StorageIndex> ColMajorMatrixAux;
0143     typedef typename sparse_eval<ColMajorMatrixAux,ResultType::RowsAtCompileTime,ResultType::ColsAtCompileTime,ColMajorMatrixAux::Flags>::type ColMajorMatrix;
0144
0145     // If the result is tall and thin (in the extreme case a column vector)
0146     // then it is faster to sort the coefficients inplace instead of transposing twice.
0147     // FIXME, the following heuristic is probably not very good.
0148     if(lhs.rows()>rhs.cols())
0149     {
0150       ColMajorMatrix resCol(lhs.rows(),rhs.cols());
0151       // perform sorted insertion
0152       internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol, true);
0153       res = resCol.markAsRValue();
0154     }
0155     else
0156     {
0157       ColMajorMatrixAux resCol(lhs.rows(),rhs.cols());
0158       // resort to transpose to sort the entries
0159       internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrixAux>(lhs, rhs, resCol, false);
0160       RowMajorMatrix resRow(resCol);
0161       res = resRow.markAsRValue();
0162     }
0163   }
0164 };
0165
0166 template<typename Lhs, typename Rhs, typename ResultType>
0167 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,ColMajor>
0168 {
0169   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
0170   {
0171     typedef SparseMatrix<typename Rhs::Scalar,RowMajor,typename ResultType::StorageIndex> RowMajorRhs;
0172     typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::StorageIndex> RowMajorRes;
0173     RowMajorRhs rhsRow = rhs;
0174     RowMajorRes resRow(lhs.rows(), rhs.cols());
0175     internal::conservative_sparse_sparse_product_impl<RowMajorRhs,Lhs,RowMajorRes>(rhsRow, lhs, resRow);
0176     res = resRow;
0177   }
0178 };
0179
0180 template<typename Lhs, typename Rhs, typename ResultType>
0181 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,RowMajor,ColMajor>
0182 {
0183   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
0184   {
0185     typedef SparseMatrix<typename Lhs::Scalar,RowMajor,typename ResultType::StorageIndex> RowMajorLhs;
0186     typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::StorageIndex> RowMajorRes;
0187     RowMajorLhs lhsRow = lhs;
0188     RowMajorRes resRow(lhs.rows(), rhs.cols());
0189     internal::conservative_sparse_sparse_product_impl<Rhs,RowMajorLhs,RowMajorRes>(rhs, lhsRow, resRow);
0190     res = resRow;
0191   }
0192 };
0193
0194 template<typename Lhs, typename Rhs, typename ResultType>
0195 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
0196 {
0197   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
0198   {
0199     typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::StorageIndex> RowMajorMatrix;
0200     RowMajorMatrix resRow(lhs.rows(), rhs.cols());
0201     internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
0202     res = resRow;
0203   }
0204 };
0205
0206
0207 template<typename Lhs, typename Rhs, typename ResultType>
0208 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
0209 {
0210   typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
0211
0212   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
0213   {
0214     typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::StorageIndex> ColMajorMatrix;
0215     ColMajorMatrix resCol(lhs.rows(), rhs.cols());
0216     internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
0217     res = resCol;
0218   }
0219 };
0220
0221 template<typename Lhs, typename Rhs, typename ResultType>
0222 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,RowMajor>
0223 {
0224   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
0225   {
0226     typedef SparseMatrix<typename Lhs::Scalar,ColMajor,typename ResultType::StorageIndex> ColMajorLhs;
0227     typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::StorageIndex> ColMajorRes;
0228     ColMajorLhs lhsCol = lhs;
0229     ColMajorRes resCol(lhs.rows(), rhs.cols());
0230     internal::conservative_sparse_sparse_product_impl<ColMajorLhs,Rhs,ColMajorRes>(lhsCol, rhs, resCol);
0231     res = resCol;
0232   }
0233 };
0234
0235 template<typename Lhs, typename Rhs, typename ResultType>
0236 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,RowMajor,RowMajor>
0237 {
0238   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
0239   {
0240     typedef SparseMatrix<typename Rhs::Scalar,ColMajor,typename ResultType::StorageIndex> ColMajorRhs;
0241     typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::StorageIndex> ColMajorRes;
0242     ColMajorRhs rhsCol = rhs;
0243     ColMajorRes resCol(lhs.rows(), rhs.cols());
0244     internal::conservative_sparse_sparse_product_impl<Lhs,ColMajorRhs,ColMajorRes>(lhs, rhsCol, resCol);
0245     res = resCol;
0246   }
0247 };
0248
0249 template<typename Lhs, typename Rhs, typename ResultType>
0250 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
0251 {
0252   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
0253   {
0254     typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::StorageIndex> RowMajorMatrix;
0255     typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::StorageIndex> ColMajorMatrix;
0256     RowMajorMatrix resRow(lhs.rows(),rhs.cols());
0257     internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
0258     // sort the non zeros:
0259     ColMajorMatrix resCol(resRow);
0260     res = resCol;
0261   }
0262 };
0263
0264 } // end namespace internal
0265
0266
0267 namespace internal {
0268
0269 template<typename Lhs, typename Rhs, typename ResultType>
0270 static void sparse_sparse_to_dense_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
0271 {
0272   typedef typename remove_all<Lhs>::type::Scalar LhsScalar;
0273   typedef typename remove_all<Rhs>::type::Scalar RhsScalar;
0274   Index cols = rhs.outerSize();
0275   eigen_assert(lhs.outerSize() == rhs.innerSize());
0276
0277   evaluator<Lhs> lhsEval(lhs);
0278   evaluator<Rhs> rhsEval(rhs);
0279
0280   for (Index j=0; j<cols; ++j)
0281   {
0282     for (typename evaluator<Rhs>::InnerIterator rhsIt(rhsEval, j); rhsIt; ++rhsIt)
0283     {
0284       RhsScalar y = rhsIt.value();
0285       Index k = rhsIt.index();
0286       for (typename evaluator<Lhs>::InnerIterator lhsIt(lhsEval, k); lhsIt; ++lhsIt)
0287       {
0288         Index i = lhsIt.index();
0289         LhsScalar x = lhsIt.value();
0290         res.coeffRef(i,j) += x * y;
0291       }
0292     }
0293   }
0294 }
0295
0296
0297 } // end namespace internal
0298
0299 namespace internal {
0300
0301 template<typename Lhs, typename Rhs, typename ResultType,
0302   int LhsStorageOrder = (traits<Lhs>::Flags&RowMajorBit) ? RowMajor : ColMajor,
0303   int RhsStorageOrder = (traits<Rhs>::Flags&RowMajorBit) ? RowMajor : ColMajor>
0304 struct sparse_sparse_to_dense_product_selector;
0305
0306 template<typename Lhs, typename Rhs, typename ResultType>
0307 struct sparse_sparse_to_dense_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor>
0308 {
0309   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
0310   {
0311     internal::sparse_sparse_to_dense_product_impl<Lhs,Rhs,ResultType>(lhs, rhs, res);
0312   }
0313 };
0314
0315 template<typename Lhs, typename Rhs, typename ResultType>
0316 struct sparse_sparse_to_dense_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor>
0317 {
0318   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
0319   {
0320     typedef SparseMatrix<typename Lhs::Scalar,ColMajor,typename ResultType::StorageIndex> ColMajorLhs;
0321     ColMajorLhs lhsCol(lhs);
0322     internal::sparse_sparse_to_dense_product_impl<ColMajorLhs,Rhs,ResultType>(lhsCol, rhs, res);
0323   }
0324 };
0325
0326 template<typename Lhs, typename Rhs, typename ResultType>
0327 struct sparse_sparse_to_dense_product_selector<Lhs,Rhs,ResultType,ColMajor,RowMajor>
0328 {
0329   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
0330   {
0331     typedef SparseMatrix<typename Rhs::Scalar,ColMajor,typename ResultType::StorageIndex> ColMajorRhs;
0332     ColMajorRhs rhsCol(rhs);
0333     internal::sparse_sparse_to_dense_product_impl<Lhs,ColMajorRhs,ResultType>(lhs, rhsCol, res);
0334   }
0335 };
0336
0337 template<typename Lhs, typename Rhs, typename ResultType>
0338 struct sparse_sparse_to_dense_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor>
0339 {
0340   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
0341   {
0342     Transpose<ResultType> trRes(res);
0343     internal::sparse_sparse_to_dense_product_impl<Rhs,Lhs,Transpose<ResultType> >(rhs, lhs, trRes);
0344   }
0345 };
0346
0347
0348 } // end namespace internal
0349
0350 } // end namespace Eigen
0351
0352 #endif // EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H