SparseSparseProduct.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-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_SPARSESPARSEPRODUCT_H
00026 #define EIGEN_SPARSESPARSEPRODUCT_H
00027 
00028 namespace internal {
00029 
00030 template<typename Lhs, typename Rhs, typename ResultType>
00031 static void sparse_product_impl2(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00032 {
00033   typedef typename remove_all<Lhs>::type::Scalar Scalar;
00034   typedef typename remove_all<Lhs>::type::Index Index;
00035 
00036   // make sure to call innerSize/outerSize since we fake the storage order.
00037   Index rows = lhs.innerSize();
00038   Index cols = rhs.outerSize();
00039   eigen_assert(lhs.outerSize() == rhs.innerSize());
00040 
00041   std::vector<bool> mask(rows,false);
00042   Matrix<Scalar,Dynamic,1> values(rows);
00043   Matrix<Index,Dynamic,1>    indices(rows);
00044 
00045   // estimate the number of non zero entries
00046   float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols()));
00047   float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
00048   float ratioRes = (std::min)(ratioLhs * avgNnzPerRhsColumn, 1.f);
00049 
00050 //  int t200 = rows/(log2(200)*1.39);
00051 //  int t = (rows*100)/139;
00052 
00053   res.resize(rows, cols);
00054   res.reserve(Index(ratioRes*rows*cols));
00055   // we compute each column of the result, one after the other
00056   for (Index j=0; j<cols; ++j)
00057   {
00058 
00059     res.startVec(j);
00060     Index nnz = 0;
00061     for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
00062     {
00063       Scalar y = rhsIt.value();
00064       Index k = rhsIt.index();
00065       for (typename Lhs::InnerIterator lhsIt(lhs, k); lhsIt; ++lhsIt)
00066       {
00067         Index i = lhsIt.index();
00068         Scalar x = lhsIt.value();
00069         if(!mask[i])
00070         {
00071           mask[i] = true;
00072 //           values[i] = x * y;
00073 //           indices[nnz] = i;
00074           ++nnz;
00075         }
00076         else
00077           values[i] += x * y;
00078       }
00079     }
00080     // FIXME reserve nnz non zeros
00081     // FIXME implement fast sort algorithms for very small nnz
00082     // if the result is sparse enough => use a quick sort
00083     // otherwise => loop through the entire vector
00084     // In order to avoid to perform an expensive log2 when the
00085     // result is clearly very sparse we use a linear bound up to 200.
00086 //     if((nnz<200 && nnz<t200) || nnz * log2(nnz) < t)
00087 //     {
00088 //       if(nnz>1) std::sort(indices.data(),indices.data()+nnz);
00089 //       for(int k=0; k<nnz; ++k)
00090 //       {
00091 //         int i = indices[k];
00092 //         res.insertBackNoCheck(j,i) = values[i];
00093 //         mask[i] = false;
00094 //       }
00095 //     }
00096 //     else
00097 //     {
00098 //       // dense path
00099 //       for(int i=0; i<rows; ++i)
00100 //       {
00101 //         if(mask[i])
00102 //         {
00103 //           mask[i] = false;
00104 //           res.insertBackNoCheck(j,i) = values[i];
00105 //         }
00106 //       }
00107 //     }
00108 
00109   }
00110   res.finalize();
00111 }
00112 
00113 // perform a pseudo in-place sparse * sparse product assuming all matrices are col major
00114 template<typename Lhs, typename Rhs, typename ResultType>
00115 static void sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00116 {
00117 //   return sparse_product_impl2(lhs,rhs,res);
00118 
00119   typedef typename remove_all<Lhs>::type::Scalar Scalar;
00120   typedef typename remove_all<Lhs>::type::Index Index;
00121 
00122   // make sure to call innerSize/outerSize since we fake the storage order.
00123   Index rows = lhs.innerSize();
00124   Index cols = rhs.outerSize();
00125   //int size = lhs.outerSize();
00126   eigen_assert(lhs.outerSize() == rhs.innerSize());
00127 
00128   // allocate a temporary buffer
00129   AmbiVector<Scalar,Index> tempVector(rows);
00130 
00131   // estimate the number of non zero entries
00132   float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols()));
00133   float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
00134   float ratioRes = (std::min)(ratioLhs * avgNnzPerRhsColumn, 1.f);
00135 
00136   // mimics a resizeByInnerOuter:
00137   if(ResultType::IsRowMajor)
00138     res.resize(cols, rows);
00139   else
00140     res.resize(rows, cols);
00141 
00142   res.reserve(Index(ratioRes*rows*cols));
00143   for (Index j=0; j<cols; ++j)
00144   {
00145     // let's do a more accurate determination of the nnz ratio for the current column j of res
00146     //float ratioColRes = (std::min)(ratioLhs * rhs.innerNonZeros(j), 1.f);
00147     // FIXME find a nice way to get the number of nonzeros of a sub matrix (here an inner vector)
00148     float ratioColRes = ratioRes;
00149     tempVector.init(ratioColRes);
00150     tempVector.setZero();
00151     for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
00152     {
00153       // FIXME should be written like this: tmp += rhsIt.value() * lhs.col(rhsIt.index())
00154       tempVector.restart();
00155       Scalar x = rhsIt.value();
00156       for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
00157       {
00158         tempVector.coeffRef(lhsIt.index()) += lhsIt.value() * x;
00159       }
00160     }
00161     res.startVec(j);
00162     for (typename AmbiVector<Scalar,Index>::Iterator it(tempVector); it; ++it)
00163       res.insertBackByOuterInner(j,it.index()) = it.value();
00164   }
00165   res.finalize();
00166 }
00167 
00168 template<typename Lhs, typename Rhs, typename ResultType,
00169   int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit,
00170   int RhsStorageOrder = traits<Rhs>::Flags&RowMajorBit,
00171   int ResStorageOrder = traits<ResultType>::Flags&RowMajorBit>
00172 struct sparse_product_selector;
00173 
00174 template<typename Lhs, typename Rhs, typename ResultType>
00175 struct sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
00176 {
00177   typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
00178 
00179   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00180   {
00181 //     std::cerr << __LINE__ << "\n";
00182     typename remove_all<ResultType>::type _res(res.rows(), res.cols());
00183     sparse_product_impl<Lhs,Rhs,ResultType>(lhs, rhs, _res);
00184     res.swap(_res);
00185   }
00186 };
00187 
00188 template<typename Lhs, typename Rhs, typename ResultType>
00189 struct sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
00190 {
00191   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00192   {
00193 //     std::cerr << __LINE__ << "\n";
00194     // we need a col-major matrix to hold the result
00195     typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
00196     SparseTemporaryType _res(res.rows(), res.cols());
00197     sparse_product_impl<Lhs,Rhs,SparseTemporaryType>(lhs, rhs, _res);
00198     res = _res;
00199   }
00200 };
00201 
00202 template<typename Lhs, typename Rhs, typename ResultType>
00203 struct sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
00204 {
00205   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00206   {
00207 //     std::cerr << __LINE__ << "\n";
00208     // let's transpose the product to get a column x column product
00209     typename remove_all<ResultType>::type _res(res.rows(), res.cols());
00210     sparse_product_impl<Rhs,Lhs,ResultType>(rhs, lhs, _res);
00211     res.swap(_res);
00212   }
00213 };
00214 
00215 template<typename Lhs, typename Rhs, typename ResultType>
00216 struct sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
00217 {
00218   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00219   {
00220 //     std::cerr << "here...\n";
00221     typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00222     ColMajorMatrix colLhs(lhs);
00223     ColMajorMatrix colRhs(rhs);
00224 //     std::cerr << "more...\n";
00225     sparse_product_impl<ColMajorMatrix,ColMajorMatrix,ResultType>(colLhs, colRhs, res);
00226 //     std::cerr << "OK.\n";
00227 
00228     // let's transpose the product to get a column x column product
00229 
00230 //     typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
00231 //     SparseTemporaryType _res(res.cols(), res.rows());
00232 //     sparse_product_impl<Rhs,Lhs,SparseTemporaryType>(rhs, lhs, _res);
00233 //     res = _res.transpose();
00234   }
00235 };
00236 
00237 // NOTE the 2 others cases (col row *) must never occur since they are caught
00238 // by ProductReturnType which transforms it to (col col *) by evaluating rhs.
00239 
00240 } // end namespace internal
00241 
00242 // sparse = sparse * sparse
00243 template<typename Derived>
00244 template<typename Lhs, typename Rhs>
00245 inline Derived& SparseMatrixBase<Derived>::operator=(const SparseSparseProduct<Lhs,Rhs>& product)
00246 {
00247 //   std::cerr << "there..." << typeid(Lhs).name() << "  " << typeid(Lhs).name() << " " << (Derived::Flags&&RowMajorBit) << "\n";
00248   internal::sparse_product_selector<
00249     typename internal::remove_all<Lhs>::type,
00250     typename internal::remove_all<Rhs>::type,
00251     Derived>::run(product.lhs(),product.rhs(),derived());
00252   return derived();
00253 }
00254 
00255 namespace internal {
00256 
00257 template<typename Lhs, typename Rhs, typename ResultType,
00258   int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit,
00259   int RhsStorageOrder = traits<Rhs>::Flags&RowMajorBit,
00260   int ResStorageOrder = traits<ResultType>::Flags&RowMajorBit>
00261 struct sparse_product_selector2;
00262 
00263 template<typename Lhs, typename Rhs, typename ResultType>
00264 struct sparse_product_selector2<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
00265 {
00266   typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
00267 
00268   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00269   {
00270     sparse_product_impl2<Lhs,Rhs,ResultType>(lhs, rhs, res);
00271   }
00272 };
00273 
00274 template<typename Lhs, typename Rhs, typename ResultType>
00275 struct sparse_product_selector2<Lhs,Rhs,ResultType,RowMajor,ColMajor,ColMajor>
00276 {
00277   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00278   {
00279       // prevent warnings until the code is fixed
00280       EIGEN_UNUSED_VARIABLE(lhs);
00281       EIGEN_UNUSED_VARIABLE(rhs);
00282       EIGEN_UNUSED_VARIABLE(res);
00283 
00284 //     typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
00285 //     RowMajorMatrix rhsRow = rhs;
00286 //     RowMajorMatrix resRow(res.rows(), res.cols());
00287 //     sparse_product_impl2<RowMajorMatrix,Lhs,RowMajorMatrix>(rhsRow, lhs, resRow);
00288 //     res = resRow;
00289   }
00290 };
00291 
00292 template<typename Lhs, typename Rhs, typename ResultType>
00293 struct sparse_product_selector2<Lhs,Rhs,ResultType,ColMajor,RowMajor,ColMajor>
00294 {
00295   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00296   {
00297     typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
00298     RowMajorMatrix lhsRow = lhs;
00299     RowMajorMatrix resRow(res.rows(), res.cols());
00300     sparse_product_impl2<Rhs,RowMajorMatrix,RowMajorMatrix>(rhs, lhsRow, resRow);
00301     res = resRow;
00302   }
00303 };
00304 
00305 template<typename Lhs, typename Rhs, typename ResultType>
00306 struct sparse_product_selector2<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
00307 {
00308   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00309   {
00310     typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
00311     RowMajorMatrix resRow(res.rows(), res.cols());
00312     sparse_product_impl2<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
00313     res = resRow;
00314   }
00315 };
00316 
00317 
00318 template<typename Lhs, typename Rhs, typename ResultType>
00319 struct sparse_product_selector2<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
00320 {
00321   typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
00322 
00323   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00324   {
00325     typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00326     ColMajorMatrix resCol(res.rows(), res.cols());
00327     sparse_product_impl2<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
00328     res = resCol;
00329   }
00330 };
00331 
00332 template<typename Lhs, typename Rhs, typename ResultType>
00333 struct sparse_product_selector2<Lhs,Rhs,ResultType,RowMajor,ColMajor,RowMajor>
00334 {
00335   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00336   {
00337     typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00338     ColMajorMatrix lhsCol = lhs;
00339     ColMajorMatrix resCol(res.rows(), res.cols());
00340     sparse_product_impl2<ColMajorMatrix,Rhs,ColMajorMatrix>(lhsCol, rhs, resCol);
00341     res = resCol;
00342   }
00343 };
00344 
00345 template<typename Lhs, typename Rhs, typename ResultType>
00346 struct sparse_product_selector2<Lhs,Rhs,ResultType,ColMajor,RowMajor,RowMajor>
00347 {
00348   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00349   {
00350     typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00351     ColMajorMatrix rhsCol = rhs;
00352     ColMajorMatrix resCol(res.rows(), res.cols());
00353     sparse_product_impl2<Lhs,ColMajorMatrix,ColMajorMatrix>(lhs, rhsCol, resCol);
00354     res = resCol;
00355   }
00356 };
00357 
00358 template<typename Lhs, typename Rhs, typename ResultType>
00359 struct sparse_product_selector2<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
00360 {
00361   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00362   {
00363     typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00364 //     ColMajorMatrix lhsTr(lhs);
00365 //     ColMajorMatrix rhsTr(rhs);
00366 //     ColMajorMatrix aux(res.rows(), res.cols());
00367 //     sparse_product_impl2<Rhs,Lhs,ColMajorMatrix>(rhs, lhs, aux);
00368 // //     ColMajorMatrix aux2 = aux.transpose();
00369 //     res = aux;
00370     typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00371     ColMajorMatrix lhsCol(lhs);
00372     ColMajorMatrix rhsCol(rhs);
00373     ColMajorMatrix resCol(res.rows(), res.cols());
00374     sparse_product_impl2<ColMajorMatrix,ColMajorMatrix,ColMajorMatrix>(lhsCol, rhsCol, resCol);
00375     res = resCol;
00376   }
00377 };
00378 
00379 } // end namespace internal
00380 
00381 template<typename Derived>
00382 template<typename Lhs, typename Rhs>
00383 inline void SparseMatrixBase<Derived>::_experimentalNewProduct(const Lhs& lhs, const Rhs& rhs)
00384 {
00385   //derived().resize(lhs.rows(), rhs.cols());
00386   internal::sparse_product_selector2<
00387     typename internal::remove_all<Lhs>::type,
00388     typename internal::remove_all<Rhs>::type,
00389     Derived>::run(lhs,rhs,derived());
00390 }
00391 
00392 // sparse * sparse
00393 template<typename Derived>
00394 template<typename OtherDerived>
00395 inline const typename SparseSparseProductReturnType<Derived,OtherDerived>::Type
00396 SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other) const
00397 {
00398   return typename SparseSparseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
00399 }
00400 
00401 #endif // EIGEN_SPARSESPARSEPRODUCT_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:39