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