GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-24 02:29:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

improve the new experimental sparse product

Browse Source
This commit is contained in:
Gael Guennebaud 2010-01-05 19:56:59 +01:00
parent eda4e98c61
commit 023e0dfb4e
2 changed files with 66 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Eigen/src/Sparse/SparseMatrix.h
View File

@ -450,6 +450,12 @@ class SparseMatrix
return *this; return *this;
} }
template<typename Lhs, typename Rhs>
inline SparseMatrix& operator=(const SparseProduct<Lhs,Rhs>& product)
{
return Base::operator=(product);
}
template<typename OtherDerived> template<typename OtherDerived>
EIGEN_DONT_INLINE SparseMatrix& operator=(const SparseMatrixBase<OtherDerived>& other) EIGEN_DONT_INLINE SparseMatrix& operator=(const SparseMatrixBase<OtherDerived>& other)
{ {

101
Eigen/src/Sparse/SparseProduct.h
View File

@ -147,13 +147,16 @@ static void ei_sparse_product_impl2(const Lhs& lhs, const Rhs& rhs, ResultType&
ei_assert(lhs.outerSize() == rhs.innerSize()); ei_assert(lhs.outerSize() == rhs.innerSize());
std::vector<bool> mask(rows,false); std::vector<bool> mask(rows,false);
Matrix<Scalar,Dynamic,1> values(rows);
Matrix<int,Dynamic,1> indices(rows);
// estimate the number of non zero entries // estimate the number of non zero entries
float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols())); float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols()));
float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols); float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f); float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f);
float ratio; int t200 = rows/(log2(200)*1.39);
int t = (rows*100)/139;
res.resize(rows, cols); res.resize(rows, cols);
res.reserve(int(ratioRes*rows*cols)); res.reserve(int(ratioRes*rows*cols));
@ -162,6 +165,7 @@ static void ei_sparse_product_impl2(const Lhs& lhs, const Rhs& rhs, ResultType&
{ {
res.startVec(j); res.startVec(j);
int nnz = 0;
for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt) for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
{ {
Scalar y = rhsIt.value(); Scalar y = rhsIt.value();
@ -173,42 +177,42 @@ static void ei_sparse_product_impl2(const Lhs& lhs, const Rhs& rhs, ResultType&
if(!mask[i]) if(!mask[i])
{ {
mask[i] = true; mask[i] = true;
values[i] = x * y; // values[i] = x * y;
res.insertBackNoCheck(j,i); // indices[nnz] = i;
++nnz;
} }
else else
res._valuePtr()[mask[i]] += x* y; values[i] += x * y;
} }
} }
// FIXME reserve nnz non zeros
// FIXME implement fast sort algorithms for very small nnz
// if the result is sparse enough => use a quick sort // if the result is sparse enough => use a quick sort
// otherwise => loop through the entire vector // otherwise => loop through the entire vector
SparseInnerVectorSet<ResultType,1> vec(res,j); // In order to avoid to perform an expensive log2 when the
// result is clearly very sparse we use a linear bound up to 200.
int nnz = vec.nonZeros(); // if((nnz<200 && nnz<t200) || nnz * log2(nnz) < t)
if(rows/1.39 > nnz * log2(nnz)) // {
{ // if(nnz>1) std::sort(indices.data(),indices.data()+nnz);
std::sort(vec._innerIndexPtr(), vec._innerIndexPtr()+vec.nonZeros()); // for(int k=0; k<nnz; ++k)
for (typename ResultType::InnerIterator it(res, j); it; ++it) // {
{ // int i = indices[k];
it.valueRef() = values[it.index()]; // res.insertBackNoCheck(j,i) = values[i];
mask[it.index()] = false; // mask[i] = false;
} // }
} // }
else // else
{ // {
// dense path // // dense path
int count = 0; // for(int i=0; i<rows; ++i)
for(int i=0; i<rows; ++i) // {
{ // if(mask[i])
if(mask[i]) // {
{ // mask[i] = false;
mask[i] = false; // res.insertBackNoCheck(j,i) = values[i];
vec._innerIndexPtr()[count] = i; // }
vec._valuePtr()[count] = i; // }
++count; // }
}
}
}
} }
res.finalize(); res.finalize();
@ -218,6 +222,8 @@ static void ei_sparse_product_impl2(const Lhs& lhs, const Rhs& rhs, ResultType&
template<typename Lhs, typename Rhs, typename ResultType> template<typename Lhs, typename Rhs, typename ResultType>
static void ei_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res) static void ei_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{ {
// return ei_sparse_product_impl2(lhs,rhs,res);
typedef typename ei_traits<typename ei_cleantype<Lhs>::type>::Scalar Scalar; typedef typename ei_traits<typename ei_cleantype<Lhs>::type>::Scalar Scalar;
// make sure to call innerSize/outerSize since we fake the storage order. // make sure to call innerSize/outerSize since we fake the storage order.
@ -274,6 +280,7 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{ {
// std::cerr << __LINE__ << "\n";
typename ei_cleantype<ResultType>::type _res(res.rows(), res.cols()); typename ei_cleantype<ResultType>::type _res(res.rows(), res.cols());
ei_sparse_product_impl<Lhs,Rhs,ResultType>(lhs, rhs, _res); ei_sparse_product_impl<Lhs,Rhs,ResultType>(lhs, rhs, _res);
res.swap(_res); res.swap(_res);
@ -285,6 +292,7 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
{ {
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{ {
// std::cerr << __LINE__ << "\n";
// we need a col-major matrix to hold the result // we need a col-major matrix to hold the result
typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType; typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
SparseTemporaryType _res(res.rows(), res.cols()); SparseTemporaryType _res(res.rows(), res.cols());
@ -298,6 +306,7 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
{ {
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{ {
// std::cerr << __LINE__ << "\n";
// let's transpose the product to get a column x column product // let's transpose the product to get a column x column product
typename ei_cleantype<ResultType>::type _res(res.rows(), res.cols()); typename ei_cleantype<ResultType>::type _res(res.rows(), res.cols());
ei_sparse_product_impl<Rhs,Lhs,ResultType>(rhs, lhs, _res); ei_sparse_product_impl<Rhs,Lhs,ResultType>(rhs, lhs, _res);
@ -310,11 +319,20 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
{ {
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{ {
// std::cerr << "here...\n";
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
ColMajorMatrix colLhs(lhs);
ColMajorMatrix colRhs(rhs);
// std::cerr << "more...\n";
ei_sparse_product_impl<ColMajorMatrix,ColMajorMatrix,ResultType>(colLhs, colRhs, res);
// std::cerr << "OK.\n";
// let's transpose the product to get a column x column product // let's transpose the product to get a column x column product
typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
SparseTemporaryType _res(res.cols(), res.rows()); // typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
ei_sparse_product_impl<Rhs,Lhs,SparseTemporaryType>(rhs, lhs, _res); // SparseTemporaryType _res(res.cols(), res.rows());
res = _res.transpose(); // ei_sparse_product_impl<Rhs,Lhs,SparseTemporaryType>(rhs, lhs, _res);
// res = _res.transpose();
} }
}; };
@ -327,6 +345,7 @@ template<typename Derived>
template<typename Lhs, typename Rhs> template<typename Lhs, typename Rhs>
inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs>& product) inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs>& product)
{ {
// std::cerr << "there..." << typeid(Lhs).name() << " " << typeid(Lhs).name() << " " << (Derived::Flags&&RowMajorBit) << "\n";
ei_sparse_product_selector< ei_sparse_product_selector<
typename ei_cleantype<Lhs>::type, typename ei_cleantype<Lhs>::type,
typename ei_cleantype<Rhs>::type, typename ei_cleantype<Rhs>::type,
@ -348,7 +367,7 @@ struct ei_sparse_product_selector2<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{ {
ei_sparse_product_impl2<Lhs,Rhs,ResultType>(lhs, rhs, res, 0); ei_sparse_product_impl2<Lhs,Rhs,ResultType>(lhs, rhs, res);
} }
}; };
@ -357,11 +376,11 @@ struct ei_sparse_product_selector2<Lhs,Rhs,ResultType,RowMajor,ColMajor,ColMajor
{ {
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{ {
typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix; // typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
RowMajorMatrix rhsRow = rhs; // RowMajorMatrix rhsRow = rhs;
RowMajorMatrix resRow(res.rows(), res.cols()); // RowMajorMatrix resRow(res.rows(), res.cols());
ei_sparse_product_impl2<RowMajorMatrix,Lhs,RowMajorMatrix>(rhsRow, lhs, resRow); // ei_sparse_product_impl2<RowMajorMatrix,Lhs,RowMajorMatrix>(rhsRow, lhs, resRow);
res = resRow; // res = resRow;
} }
}; };