$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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(Index k=0; k<nnz; ++k) 00070 { 00071 Index 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 Index t200 = rows/(log2(200)*1.39); 00080 Index 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(Index k=0; k<nnz; ++k) 00094 { 00095 Index i = indices[k]; 00096 res.insertBackByOuterInner(j,i) = values[i]; 00097 mask[i] = false; 00098 } 00099 } 00100 else 00101 { 00102 // dense path 00103 for(Index 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,typename ResultType::Index> RowMajorMatrix; 00138 typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> 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,typename ResultType::Index> 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,typename ResultType::Index> 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,typename ResultType::Index> 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,typename ResultType::Index> 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,typename ResultType::Index> 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,typename ResultType::Index> 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,typename ResultType::Index> RowMajorMatrix; 00232 typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> 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