ConservativeSparseSparseProduct.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
00011 #define EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal {
00016 
00017 template<typename Lhs, typename Rhs, typename ResultType>
00018 static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00019 {
00020   typedef typename remove_all<Lhs>::type::Scalar Scalar;
00021   typedef typename remove_all<Lhs>::type::Index Index;
00022 
00023   // make sure to call innerSize/outerSize since we fake the storage order.
00024   Index rows = lhs.innerSize();
00025   Index cols = rhs.outerSize();
00026   eigen_assert(lhs.outerSize() == rhs.innerSize());
00027 
00028   std::vector<bool> mask(rows,false);
00029   Matrix<Scalar,Dynamic,1> values(rows);
00030   Matrix<Index,Dynamic,1>  indices(rows);
00031 
00032   // estimate the number of non zero entries
00033   // given a rhs column containing Y non zeros, we assume that the respective Y columns
00034   // of the lhs differs in average of one non zeros, thus the number of non zeros for
00035   // the product of a rhs column with the lhs is X+Y where X is the average number of non zero
00036   // per column of the lhs.
00037   // Therefore, we have nnz(lhs*rhs) = nnz(lhs) + nnz(rhs)
00038   Index estimated_nnz_prod = lhs.nonZeros() + rhs.nonZeros();
00039 
00040   res.setZero();
00041   res.reserve(Index(estimated_nnz_prod));
00042   // we compute each column of the result, one after the other
00043   for (Index j=0; j<cols; ++j)
00044   {
00045 
00046     res.startVec(j);
00047     Index nnz = 0;
00048     for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
00049     {
00050       Scalar y = rhsIt.value();
00051       Index k = rhsIt.index();
00052       for (typename Lhs::InnerIterator lhsIt(lhs, k); lhsIt; ++lhsIt)
00053       {
00054         Index i = lhsIt.index();
00055         Scalar x = lhsIt.value();
00056         if(!mask[i])
00057         {
00058           mask[i] = true;
00059           values[i] = x * y;
00060           indices[nnz] = i;
00061           ++nnz;
00062         }
00063         else
00064           values[i] += x * y;
00065       }
00066     }
00067 
00068     // unordered insertion
00069     for(int k=0; k<nnz; ++k)
00070     {
00071       int i = indices[k];
00072       res.insertBackByOuterInnerUnordered(j,i) = values[i];
00073       mask[i] = false;
00074     }
00075 
00076 #if 0
00077     // alternative ordered insertion code:
00078 
00079     int t200 = rows/(log2(200)*1.39);
00080     int t = (rows*100)/139;
00081 
00082     // FIXME reserve nnz non zeros
00083     // FIXME implement fast sort algorithms for very small nnz
00084     // if the result is sparse enough => use a quick sort
00085     // otherwise => loop through the entire vector
00086     // In order to avoid to perform an expensive log2 when the
00087     // result is clearly very sparse we use a linear bound up to 200.
00088     //if((nnz<200 && nnz<t200) || nnz * log2(nnz) < t)
00089     //res.startVec(j);
00090     if(true)
00091     {
00092       if(nnz>1) std::sort(indices.data(),indices.data()+nnz);
00093       for(int k=0; k<nnz; ++k)
00094       {
00095         int i = indices[k];
00096         res.insertBackByOuterInner(j,i) = values[i];
00097         mask[i] = false;
00098       }
00099     }
00100     else
00101     {
00102       // dense path
00103       for(int i=0; i<rows; ++i)
00104       {
00105         if(mask[i])
00106         {
00107           mask[i] = false;
00108           res.insertBackByOuterInner(j,i) = values[i];
00109         }
00110       }
00111     }
00112 #endif
00113 
00114   }
00115   res.finalize();
00116 }
00117 
00118 
00119 } // end namespace internal
00120 
00121 namespace internal {
00122 
00123 template<typename Lhs, typename Rhs, typename ResultType,
00124   int LhsStorageOrder = (traits<Lhs>::Flags&RowMajorBit) ? RowMajor : ColMajor,
00125   int RhsStorageOrder = (traits<Rhs>::Flags&RowMajorBit) ? RowMajor : ColMajor,
00126   int ResStorageOrder = (traits<ResultType>::Flags&RowMajorBit) ? RowMajor : ColMajor>
00127 struct conservative_sparse_sparse_product_selector;
00128 
00129 template<typename Lhs, typename Rhs, typename ResultType>
00130 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
00131 {
00132   typedef typename remove_all<Lhs>::type LhsCleaned;
00133   typedef typename LhsCleaned::Scalar Scalar;
00134 
00135   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00136   {
00137     typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
00138     typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00139     ColMajorMatrix resCol(lhs.rows(),rhs.cols());
00140     internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
00141     // sort the non zeros:
00142     RowMajorMatrix resRow(resCol);
00143     res = resRow;
00144   }
00145 };
00146 
00147 template<typename Lhs, typename Rhs, typename ResultType>
00148 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,ColMajor>
00149 {
00150   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00151   {
00152      typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
00153      RowMajorMatrix rhsRow = rhs;
00154      RowMajorMatrix resRow(lhs.rows(), rhs.cols());
00155      internal::conservative_sparse_sparse_product_impl<RowMajorMatrix,Lhs,RowMajorMatrix>(rhsRow, lhs, resRow);
00156      res = resRow;
00157   }
00158 };
00159 
00160 template<typename Lhs, typename Rhs, typename ResultType>
00161 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,RowMajor,ColMajor>
00162 {
00163   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00164   {
00165     typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
00166     RowMajorMatrix lhsRow = lhs;
00167     RowMajorMatrix resRow(lhs.rows(), rhs.cols());
00168     internal::conservative_sparse_sparse_product_impl<Rhs,RowMajorMatrix,RowMajorMatrix>(rhs, lhsRow, resRow);
00169     res = resRow;
00170   }
00171 };
00172 
00173 template<typename Lhs, typename Rhs, typename ResultType>
00174 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
00175 {
00176   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00177   {
00178     typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
00179     RowMajorMatrix resRow(lhs.rows(), rhs.cols());
00180     internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
00181     res = resRow;
00182   }
00183 };
00184 
00185 
00186 template<typename Lhs, typename Rhs, typename ResultType>
00187 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
00188 {
00189   typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
00190 
00191   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00192   {
00193     typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00194     ColMajorMatrix resCol(lhs.rows(), rhs.cols());
00195     internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
00196     res = resCol;
00197   }
00198 };
00199 
00200 template<typename Lhs, typename Rhs, typename ResultType>
00201 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,RowMajor>
00202 {
00203   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00204   {
00205     typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00206     ColMajorMatrix lhsCol = lhs;
00207     ColMajorMatrix resCol(lhs.rows(), rhs.cols());
00208     internal::conservative_sparse_sparse_product_impl<ColMajorMatrix,Rhs,ColMajorMatrix>(lhsCol, rhs, resCol);
00209     res = resCol;
00210   }
00211 };
00212 
00213 template<typename Lhs, typename Rhs, typename ResultType>
00214 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,RowMajor,RowMajor>
00215 {
00216   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00217   {
00218     typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00219     ColMajorMatrix rhsCol = rhs;
00220     ColMajorMatrix resCol(lhs.rows(), rhs.cols());
00221     internal::conservative_sparse_sparse_product_impl<Lhs,ColMajorMatrix,ColMajorMatrix>(lhs, rhsCol, resCol);
00222     res = resCol;
00223   }
00224 };
00225 
00226 template<typename Lhs, typename Rhs, typename ResultType>
00227 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
00228 {
00229   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00230   {
00231     typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
00232     typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00233     RowMajorMatrix resRow(lhs.rows(),rhs.cols());
00234     internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
00235     // sort the non zeros:
00236     ColMajorMatrix resCol(resRow);
00237     res = resCol;
00238   }
00239 };
00240 
00241 } // end namespace internal
00242 
00243 } // end namespace Eigen
00244 
00245 #endif // EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:36:56