From 1927b4dff513f66866de205a46c66a1f2c877d01 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 25 Jun 2010 10:26:24 +0200 Subject: [PATCH] Fix use of nesting types in SparseTranspose and split the big SparseProduct.h file --- Eigen/Sparse | 2 + Eigen/src/Sparse/SparseDenseProduct.h | 111 ++++++ Eigen/src/Sparse/SparseMatrix.h | 2 +- Eigen/src/Sparse/SparseMatrixBase.h | 4 +- Eigen/src/Sparse/SparseProduct.h | 483 +------------------------ Eigen/src/Sparse/SparseSparseProduct.h | 390 ++++++++++++++++++++ Eigen/src/Sparse/SparseTranspose.h | 24 +- Eigen/src/Sparse/SparseUtil.h | 8 +- Eigen/src/Sparse/SparseVector.h | 2 +- test/sparse_product.cpp | 10 +- 10 files changed, 541 insertions(+), 495 deletions(-) create mode 100644 Eigen/src/Sparse/SparseDenseProduct.h create mode 100644 Eigen/src/Sparse/SparseSparseProduct.h diff --git a/Eigen/Sparse b/Eigen/Sparse index 864f194f4..1167ce78f 100644 --- a/Eigen/Sparse +++ b/Eigen/Sparse @@ -45,6 +45,8 @@ struct Sparse {}; #include "src/Sparse/SparseRedux.h" #include "src/Sparse/SparseFuzzy.h" #include "src/Sparse/SparseProduct.h" +#include "src/Sparse/SparseSparseProduct.h" +#include "src/Sparse/SparseDenseProduct.h" #include "src/Sparse/SparseDiagonalProduct.h" #include "src/Sparse/SparseTriangularView.h" #include "src/Sparse/SparseSelfAdjointView.h" diff --git a/Eigen/src/Sparse/SparseDenseProduct.h b/Eigen/src/Sparse/SparseDenseProduct.h new file mode 100644 index 000000000..8eaa60a26 --- /dev/null +++ b/Eigen/src/Sparse/SparseDenseProduct.h @@ -0,0 +1,111 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2010 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_SPARSEDENSEPRODUCT_H +#define EIGEN_SPARSEDENSEPRODUCT_H + +template +struct ei_traits > + : ei_traits, Lhs, Rhs> > +{ + typedef Dense StorageKind; + typedef MatrixXpr XprKind; +}; + +template +class SparseTimeDenseProduct + : public ProductBase, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseTimeDenseProduct) + + SparseTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template void scaleAndAddTo(Dest& dest, Scalar alpha) const + { + typedef typename ei_cleantype::type _Lhs; + typedef typename ei_cleantype::type _Rhs; + typedef typename _Lhs::InnerIterator LhsInnerIterator; + enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit }; + for(Index j=0; j dest_j(dest.row(LhsIsRowMajor ? j : 0)); + for(LhsInnerIterator it(m_lhs,j); it ;++it) + { + if(LhsIsRowMajor) dest_j += (alpha*it.value()) * m_rhs.row(it.index()); + else if(Rhs::ColsAtCompileTime==1) dest.coeffRef(it.index()) += it.value() * rhs_j; + else dest.row(it.index()) += (alpha*it.value()) * m_rhs.row(j); + } + } + } + + private: + SparseTimeDenseProduct& operator=(const SparseTimeDenseProduct&); +}; + + +// dense = dense * sparse +template +struct ei_traits > + : ei_traits, Lhs, Rhs> > +{ + typedef Dense StorageKind; +}; + +template +class DenseTimeSparseProduct + : public ProductBase, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseProduct) + + DenseTimeSparseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template void scaleAndAddTo(Dest& dest, Scalar alpha) const + { + typedef typename ei_cleantype::type _Rhs; + typedef typename _Rhs::InnerIterator RhsInnerIterator; + enum { RhsIsRowMajor = (_Rhs::Flags&RowMajorBit)==RowMajorBit }; + for(Index j=0; j +template +inline const SparseTimeDenseProduct +SparseMatrixBase::operator*(const MatrixBase &other) const +{ + return SparseTimeDenseProduct(derived(), other.derived()); +} + +#endif // EIGEN_SPARSEDENSEPRODUCT_H diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index 53335e00f..ecf8efbe7 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -429,7 +429,7 @@ class SparseMatrix #ifndef EIGEN_PARSED_BY_DOXYGEN template - inline SparseMatrix& operator=(const SparseProduct& product) + inline SparseMatrix& operator=(const SparseSparseProduct& product) { return Base::operator=(product); } diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index 8507fd334..12a1cb538 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -252,7 +252,7 @@ template class SparseMatrixBase : public EigenBase } template - inline Derived& operator=(const SparseProduct& product); + inline Derived& operator=(const SparseSparseProduct& product); template inline void _experimentalNewProduct(const Lhs& lhs, const Rhs& rhs); @@ -348,7 +348,7 @@ template class SparseMatrixBase : public EigenBase // sparse * sparse template - const typename SparseProductReturnType::Type + const typename SparseSparseProductReturnType::Type operator*(const SparseMatrixBase &other) const; // sparse * diagonal diff --git a/Eigen/src/Sparse/SparseProduct.h b/Eigen/src/Sparse/SparseProduct.h index 0927fc4c5..ef6ef3c02 100644 --- a/Eigen/src/Sparse/SparseProduct.h +++ b/Eigen/src/Sparse/SparseProduct.h @@ -26,18 +26,16 @@ #define EIGEN_SPARSEPRODUCT_H template -struct SparseProductReturnType +struct SparseSparseProductReturnType { typedef typename ei_traits::Scalar Scalar; enum { LhsRowMajor = ei_traits::Flags & RowMajorBit, RhsRowMajor = ei_traits::Flags & RowMajorBit, - TransposeRhs = /*false,*/ (!LhsRowMajor) && RhsRowMajor, - TransposeLhs = /*false*/ LhsRowMajor && (!RhsRowMajor) + TransposeRhs = (!LhsRowMajor) && RhsRowMajor, + TransposeLhs = LhsRowMajor && (!RhsRowMajor) }; - // FIXME if we transpose let's evaluate to a LinkedVectorMatrix since it is the - // type of the temporary to perform the transpose op typedef typename ei_meta_if, const typename ei_nested::type>::ret LhsNested; @@ -46,11 +44,11 @@ struct SparseProductReturnType SparseMatrix, const typename ei_nested::type>::ret RhsNested; - typedef SparseProduct Type; + typedef SparseSparseProduct Type; }; template -struct ei_traits > +struct ei_traits > { typedef MatrixXpr XprKind; // clean the nested types: @@ -66,13 +64,13 @@ struct ei_traits > LhsFlags = _LhsNested::Flags, RhsFlags = _RhsNested::Flags, - RowsAtCompileTime = _LhsNested::RowsAtCompileTime, - ColsAtCompileTime = _RhsNested::ColsAtCompileTime, - InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime), - + RowsAtCompileTime = _LhsNested::RowsAtCompileTime, + ColsAtCompileTime = _RhsNested::ColsAtCompileTime, MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime, MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime, + InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime), + EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit), RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit), @@ -85,28 +83,26 @@ struct ei_traits > }; typedef Sparse StorageKind; - - typedef SparseMatrixBase > Base; }; template -class SparseProduct : ei_no_assignment_operator, - public ei_traits >::Base +class SparseSparseProduct : ei_no_assignment_operator, + public SparseMatrixBase > { public: - typedef typename ei_traits >::Base Base; - EIGEN_DENSE_PUBLIC_INTERFACE(SparseProduct) + typedef SparseMatrixBase Base; + EIGEN_DENSE_PUBLIC_INTERFACE(SparseSparseProduct) private: - typedef typename ei_traits::_LhsNested _LhsNested; - typedef typename ei_traits::_RhsNested _RhsNested; + typedef typename ei_traits::_LhsNested _LhsNested; + typedef typename ei_traits::_RhsNested _RhsNested; public: template - EIGEN_STRONG_INLINE SparseProduct(const Lhs& lhs, const Rhs& rhs) + EIGEN_STRONG_INLINE SparseSparseProduct(const Lhs& lhs, const Rhs& rhs) : m_lhs(lhs), m_rhs(rhs) { ei_assert(lhs.cols() == rhs.rows()); @@ -139,451 +135,4 @@ class SparseProduct : ei_no_assignment_operator, RhsNested m_rhs; }; -template -static void ei_sparse_product_impl2(const Lhs& lhs, const Rhs& rhs, ResultType& res) -{ - typedef typename ei_cleantype::type::Scalar Scalar; - typedef typename ei_cleantype::type::Index Index; - - // make sure to call innerSize/outerSize since we fake the storage order. - Index rows = lhs.innerSize(); - Index cols = rhs.outerSize(); - ei_assert(lhs.outerSize() == rhs.innerSize()); - - std::vector mask(rows,false); - Matrix values(rows); - Matrix indices(rows); - - // estimate the number of non zero entries - float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols())); - float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols); - float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f); - -// int t200 = rows/(log2(200)*1.39); -// int t = (rows*100)/139; - - res.resize(rows, cols); - res.reserve(Index(ratioRes*rows*cols)); - // we compute each column of the result, one after the other - for (Index j=0; j use a quick sort - // otherwise => loop through the entire vector - // In order to avoid to perform an expensive log2 when the - // result is clearly very sparse we use a linear bound up to 200. -// if((nnz<200 && nnz1) std::sort(indices.data(),indices.data()+nnz); -// for(int k=0; k -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_cleantype::type::Scalar Scalar; - typedef typename ei_cleantype::type::Index Index; - - // make sure to call innerSize/outerSize since we fake the storage order. - Index rows = lhs.innerSize(); - Index cols = rhs.outerSize(); - //int size = lhs.outerSize(); - ei_assert(lhs.outerSize() == rhs.innerSize()); - - // allocate a temporary buffer - AmbiVector tempVector(rows); - - // estimate the number of non zero entries - float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols())); - float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols); - float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f); - - res.resize(rows, cols); - res.reserve(Index(ratioRes*rows*cols)); - for (Index j=0; j::Iterator it(tempVector); it; ++it) - res.insertBackByOuterInner(j,it.index()) = it.value(); - } - res.finalize(); -} - -template::Flags&RowMajorBit, - int RhsStorageOrder = ei_traits::Flags&RowMajorBit, - int ResStorageOrder = ei_traits::Flags&RowMajorBit> -struct ei_sparse_product_selector; - -template -struct ei_sparse_product_selector -{ - typedef typename ei_traits::type>::Scalar Scalar; - - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { -// std::cerr << __LINE__ << "\n"; - typename ei_cleantype::type _res(res.rows(), res.cols()); - ei_sparse_product_impl(lhs, rhs, _res); - res.swap(_res); - } -}; - -template -struct ei_sparse_product_selector -{ - 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 - typedef SparseMatrix SparseTemporaryType; - SparseTemporaryType _res(res.rows(), res.cols()); - ei_sparse_product_impl(lhs, rhs, _res); - res = _res; - } -}; - -template -struct ei_sparse_product_selector -{ - 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 - typename ei_cleantype::type _res(res.rows(), res.cols()); - ei_sparse_product_impl(rhs, lhs, _res); - res.swap(_res); - } -}; - -template -struct ei_sparse_product_selector -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { -// std::cerr << "here...\n"; - typedef SparseMatrix ColMajorMatrix; - ColMajorMatrix colLhs(lhs); - ColMajorMatrix colRhs(rhs); -// std::cerr << "more...\n"; - ei_sparse_product_impl(colLhs, colRhs, res); -// std::cerr << "OK.\n"; - - // let's transpose the product to get a column x column product - -// typedef SparseMatrix SparseTemporaryType; -// SparseTemporaryType _res(res.cols(), res.rows()); -// ei_sparse_product_impl(rhs, lhs, _res); -// res = _res.transpose(); - } -}; - -// NOTE the 2 others cases (col row *) must never occurs since they are caught -// by ProductReturnType which transform it to (col col *) by evaluating rhs. - - -// sparse = sparse * sparse -template -template -inline Derived& SparseMatrixBase::operator=(const SparseProduct& product) -{ -// std::cerr << "there..." << typeid(Lhs).name() << " " << typeid(Lhs).name() << " " << (Derived::Flags&&RowMajorBit) << "\n"; - ei_sparse_product_selector< - typename ei_cleantype::type, - typename ei_cleantype::type, - Derived>::run(product.lhs(),product.rhs(),derived()); - return derived(); -} - - -template::Flags&RowMajorBit, - int RhsStorageOrder = ei_traits::Flags&RowMajorBit, - int ResStorageOrder = ei_traits::Flags&RowMajorBit> -struct ei_sparse_product_selector2; - -template -struct ei_sparse_product_selector2 -{ - typedef typename ei_traits::type>::Scalar Scalar; - - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - ei_sparse_product_impl2(lhs, rhs, res); - } -}; - -template -struct ei_sparse_product_selector2 -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - // prevent warnings until the code is fixed - EIGEN_UNUSED_VARIABLE(lhs); - EIGEN_UNUSED_VARIABLE(rhs); - EIGEN_UNUSED_VARIABLE(res); - -// typedef SparseMatrix RowMajorMatrix; -// RowMajorMatrix rhsRow = rhs; -// RowMajorMatrix resRow(res.rows(), res.cols()); -// ei_sparse_product_impl2(rhsRow, lhs, resRow); -// res = resRow; - } -}; - -template -struct ei_sparse_product_selector2 -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix RowMajorMatrix; - RowMajorMatrix lhsRow = lhs; - RowMajorMatrix resRow(res.rows(), res.cols()); - ei_sparse_product_impl2(rhs, lhsRow, resRow); - res = resRow; - } -}; - -template -struct ei_sparse_product_selector2 -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix RowMajorMatrix; - RowMajorMatrix resRow(res.rows(), res.cols()); - ei_sparse_product_impl2(rhs, lhs, resRow); - res = resRow; - } -}; - - -template -struct ei_sparse_product_selector2 -{ - typedef typename ei_traits::type>::Scalar Scalar; - - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix ColMajorMatrix; - ColMajorMatrix resCol(res.rows(), res.cols()); - ei_sparse_product_impl2(lhs, rhs, resCol); - res = resCol; - } -}; - -template -struct ei_sparse_product_selector2 -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix ColMajorMatrix; - ColMajorMatrix lhsCol = lhs; - ColMajorMatrix resCol(res.rows(), res.cols()); - ei_sparse_product_impl2(lhsCol, rhs, resCol); - res = resCol; - } -}; - -template -struct ei_sparse_product_selector2 -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix ColMajorMatrix; - ColMajorMatrix rhsCol = rhs; - ColMajorMatrix resCol(res.rows(), res.cols()); - ei_sparse_product_impl2(lhs, rhsCol, resCol); - res = resCol; - } -}; - -template -struct ei_sparse_product_selector2 -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix ColMajorMatrix; -// ColMajorMatrix lhsTr(lhs); -// ColMajorMatrix rhsTr(rhs); -// ColMajorMatrix aux(res.rows(), res.cols()); -// ei_sparse_product_impl2(rhs, lhs, aux); -// // ColMajorMatrix aux2 = aux.transpose(); -// res = aux; - typedef SparseMatrix ColMajorMatrix; - ColMajorMatrix lhsCol(lhs); - ColMajorMatrix rhsCol(rhs); - ColMajorMatrix resCol(res.rows(), res.cols()); - ei_sparse_product_impl2(lhsCol, rhsCol, resCol); - res = resCol; - } -}; - -template -template -inline void SparseMatrixBase::_experimentalNewProduct(const Lhs& lhs, const Rhs& rhs) -{ - //derived().resize(lhs.rows(), rhs.cols()); - ei_sparse_product_selector2< - typename ei_cleantype::type, - typename ei_cleantype::type, - Derived>::run(lhs,rhs,derived()); -} - - - -template -struct ei_traits > - : ei_traits, Lhs, Rhs> > -{ - typedef Dense StorageKind; - typedef MatrixXpr XprKind; -}; - -template -class SparseTimeDenseProduct - : public ProductBase, Lhs, Rhs> -{ - public: - EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseTimeDenseProduct) - - SparseTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) - {} - - template void scaleAndAddTo(Dest& dest, Scalar alpha) const - { - typedef typename ei_cleantype::type _Lhs; - typedef typename ei_cleantype::type _Rhs; - typedef typename _Lhs::InnerIterator LhsInnerIterator; - enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit }; - for(Index j=0; j dest_j(dest.row(LhsIsRowMajor ? j : 0)); - for(LhsInnerIterator it(m_lhs,j); it ;++it) - { - if(LhsIsRowMajor) dest_j += (alpha*it.value()) * m_rhs.row(it.index()); - else if(Rhs::ColsAtCompileTime==1) dest.coeffRef(it.index()) += it.value() * rhs_j; - else dest.row(it.index()) += (alpha*it.value()) * m_rhs.row(j); - } - } - } - - private: - SparseTimeDenseProduct& operator=(const SparseTimeDenseProduct&); -}; - - -// dense = dense * sparse -template -struct ei_traits > - : ei_traits, Lhs, Rhs> > -{ - typedef Dense StorageKind; -}; - -template -class DenseTimeSparseProduct - : public ProductBase, Lhs, Rhs> -{ - public: - EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseProduct) - - DenseTimeSparseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) - {} - - template void scaleAndAddTo(Dest& dest, Scalar alpha) const - { - typedef typename ei_cleantype::type _Rhs; - typedef typename _Rhs::InnerIterator RhsInnerIterator; - enum { RhsIsRowMajor = (_Rhs::Flags&RowMajorBit)==RowMajorBit }; - for(Index j=0; j -template -inline const typename SparseProductReturnType::Type -SparseMatrixBase::operator*(const SparseMatrixBase &other) const -{ - return typename SparseProductReturnType::Type(derived(), other.derived()); -} - -// sparse * dense -template -template -inline const SparseTimeDenseProduct -SparseMatrixBase::operator*(const MatrixBase &other) const -{ - return SparseTimeDenseProduct(derived(), other.derived()); -} - #endif // EIGEN_SPARSEPRODUCT_H diff --git a/Eigen/src/Sparse/SparseSparseProduct.h b/Eigen/src/Sparse/SparseSparseProduct.h new file mode 100644 index 000000000..c8724c118 --- /dev/null +++ b/Eigen/src/Sparse/SparseSparseProduct.h @@ -0,0 +1,390 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2010 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_SPARSESPARSEPRODUCT_H +#define EIGEN_SPARSESPARSEPRODUCT_H + +template +static void ei_sparse_product_impl2(const Lhs& lhs, const Rhs& rhs, ResultType& res) +{ + typedef typename ei_cleantype::type::Scalar Scalar; + typedef typename ei_cleantype::type::Index Index; + + // make sure to call innerSize/outerSize since we fake the storage order. + Index rows = lhs.innerSize(); + Index cols = rhs.outerSize(); + ei_assert(lhs.outerSize() == rhs.innerSize()); + + std::vector mask(rows,false); + Matrix values(rows); + Matrix indices(rows); + + // estimate the number of non zero entries + float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols())); + float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols); + float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f); + +// int t200 = rows/(log2(200)*1.39); +// int t = (rows*100)/139; + + res.resize(rows, cols); + res.reserve(Index(ratioRes*rows*cols)); + // we compute each column of the result, one after the other + for (Index j=0; j use a quick sort + // otherwise => loop through the entire vector + // In order to avoid to perform an expensive log2 when the + // result is clearly very sparse we use a linear bound up to 200. +// if((nnz<200 && nnz1) std::sort(indices.data(),indices.data()+nnz); +// for(int k=0; k +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_cleantype::type::Scalar Scalar; + typedef typename ei_cleantype::type::Index Index; + + // make sure to call innerSize/outerSize since we fake the storage order. + Index rows = lhs.innerSize(); + Index cols = rhs.outerSize(); + //int size = lhs.outerSize(); + ei_assert(lhs.outerSize() == rhs.innerSize()); + + // allocate a temporary buffer + AmbiVector tempVector(rows); + + // estimate the number of non zero entries + float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols())); + float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols); + float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f); + + res.resize(rows, cols); + res.reserve(Index(ratioRes*rows*cols)); + for (Index j=0; j::Iterator it(tempVector); it; ++it) + res.insertBackByOuterInner(j,it.index()) = it.value(); + } + res.finalize(); +} + +template::Flags&RowMajorBit, + int RhsStorageOrder = ei_traits::Flags&RowMajorBit, + int ResStorageOrder = ei_traits::Flags&RowMajorBit> +struct ei_sparse_product_selector; + +template +struct ei_sparse_product_selector +{ + typedef typename ei_traits::type>::Scalar Scalar; + + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { +// std::cerr << __LINE__ << "\n"; + typename ei_cleantype::type _res(res.rows(), res.cols()); + ei_sparse_product_impl(lhs, rhs, _res); + res.swap(_res); + } +}; + +template +struct ei_sparse_product_selector +{ + 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 + typedef SparseMatrix SparseTemporaryType; + SparseTemporaryType _res(res.rows(), res.cols()); + ei_sparse_product_impl(lhs, rhs, _res); + res = _res; + } +}; + +template +struct ei_sparse_product_selector +{ + 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 + typename ei_cleantype::type _res(res.rows(), res.cols()); + ei_sparse_product_impl(rhs, lhs, _res); + res.swap(_res); + } +}; + +template +struct ei_sparse_product_selector +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { +// std::cerr << "here...\n"; + typedef SparseMatrix ColMajorMatrix; + ColMajorMatrix colLhs(lhs); + ColMajorMatrix colRhs(rhs); +// std::cerr << "more...\n"; + ei_sparse_product_impl(colLhs, colRhs, res); +// std::cerr << "OK.\n"; + + // let's transpose the product to get a column x column product + +// typedef SparseMatrix SparseTemporaryType; +// SparseTemporaryType _res(res.cols(), res.rows()); +// ei_sparse_product_impl(rhs, lhs, _res); +// res = _res.transpose(); + } +}; + +// NOTE the 2 others cases (col row *) must never occurs since they are caught +// by ProductReturnType which transform it to (col col *) by evaluating rhs. + + +// sparse = sparse * sparse +template +template +inline Derived& SparseMatrixBase::operator=(const SparseSparseProduct& product) +{ +// std::cerr << "there..." << typeid(Lhs).name() << " " << typeid(Lhs).name() << " " << (Derived::Flags&&RowMajorBit) << "\n"; + ei_sparse_product_selector< + typename ei_cleantype::type, + typename ei_cleantype::type, + Derived>::run(product.lhs(),product.rhs(),derived()); + return derived(); +} + + +template::Flags&RowMajorBit, + int RhsStorageOrder = ei_traits::Flags&RowMajorBit, + int ResStorageOrder = ei_traits::Flags&RowMajorBit> +struct ei_sparse_product_selector2; + +template +struct ei_sparse_product_selector2 +{ + typedef typename ei_traits::type>::Scalar Scalar; + + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + ei_sparse_product_impl2(lhs, rhs, res); + } +}; + +template +struct ei_sparse_product_selector2 +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + // prevent warnings until the code is fixed + EIGEN_UNUSED_VARIABLE(lhs); + EIGEN_UNUSED_VARIABLE(rhs); + EIGEN_UNUSED_VARIABLE(res); + +// typedef SparseMatrix RowMajorMatrix; +// RowMajorMatrix rhsRow = rhs; +// RowMajorMatrix resRow(res.rows(), res.cols()); +// ei_sparse_product_impl2(rhsRow, lhs, resRow); +// res = resRow; + } +}; + +template +struct ei_sparse_product_selector2 +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix RowMajorMatrix; + RowMajorMatrix lhsRow = lhs; + RowMajorMatrix resRow(res.rows(), res.cols()); + ei_sparse_product_impl2(rhs, lhsRow, resRow); + res = resRow; + } +}; + +template +struct ei_sparse_product_selector2 +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix RowMajorMatrix; + RowMajorMatrix resRow(res.rows(), res.cols()); + ei_sparse_product_impl2(rhs, lhs, resRow); + res = resRow; + } +}; + + +template +struct ei_sparse_product_selector2 +{ + typedef typename ei_traits::type>::Scalar Scalar; + + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix ColMajorMatrix; + ColMajorMatrix resCol(res.rows(), res.cols()); + ei_sparse_product_impl2(lhs, rhs, resCol); + res = resCol; + } +}; + +template +struct ei_sparse_product_selector2 +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix ColMajorMatrix; + ColMajorMatrix lhsCol = lhs; + ColMajorMatrix resCol(res.rows(), res.cols()); + ei_sparse_product_impl2(lhsCol, rhs, resCol); + res = resCol; + } +}; + +template +struct ei_sparse_product_selector2 +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix ColMajorMatrix; + ColMajorMatrix rhsCol = rhs; + ColMajorMatrix resCol(res.rows(), res.cols()); + ei_sparse_product_impl2(lhs, rhsCol, resCol); + res = resCol; + } +}; + +template +struct ei_sparse_product_selector2 +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix ColMajorMatrix; +// ColMajorMatrix lhsTr(lhs); +// ColMajorMatrix rhsTr(rhs); +// ColMajorMatrix aux(res.rows(), res.cols()); +// ei_sparse_product_impl2(rhs, lhs, aux); +// // ColMajorMatrix aux2 = aux.transpose(); +// res = aux; + typedef SparseMatrix ColMajorMatrix; + ColMajorMatrix lhsCol(lhs); + ColMajorMatrix rhsCol(rhs); + ColMajorMatrix resCol(res.rows(), res.cols()); + ei_sparse_product_impl2(lhsCol, rhsCol, resCol); + res = resCol; + } +}; + +template +template +inline void SparseMatrixBase::_experimentalNewProduct(const Lhs& lhs, const Rhs& rhs) +{ + //derived().resize(lhs.rows(), rhs.cols()); + ei_sparse_product_selector2< + typename ei_cleantype::type, + typename ei_cleantype::type, + Derived>::run(lhs,rhs,derived()); +} + +// sparse * sparse +template +template +inline const typename SparseSparseProductReturnType::Type +SparseMatrixBase::operator*(const SparseMatrixBase &other) const +{ + return typename SparseSparseProductReturnType::Type(derived(), other.derived()); +} + +#endif // EIGEN_SPARSESPARSEPRODUCT_H diff --git a/Eigen/src/Sparse/SparseTranspose.h b/Eigen/src/Sparse/SparseTranspose.h index d6b631d4d..79e0a04db 100644 --- a/Eigen/src/Sparse/SparseTranspose.h +++ b/Eigen/src/Sparse/SparseTranspose.h @@ -28,6 +28,7 @@ template class TransposeImpl : public SparseMatrixBase > { + typedef typename ei_cleantype::type _MatrixTypeNested; public: EIGEN_SPARSE_PUBLIC_INTERFACE(Transpose) @@ -36,24 +37,12 @@ template class TransposeImpl class ReverseInnerIterator; inline Index nonZeros() const { return derived().nestedExpression().nonZeros(); } - - // FIXME should be keep them ? - inline Scalar& coeffRef(Index row, Index col) - { return const_cast_derived().nestedExpression().coeffRef(col, row); } - - inline const Scalar coeff(Index row, Index col) const - { return derived().nestedExpression().coeff(col, row); } - - inline const Scalar coeff(Index index) const - { return derived().nestedExpression().coeff(index); } - - inline Scalar& coeffRef(Index index) - { return const_cast_derived().nestedExpression().coeffRef(index); } }; -template class TransposeImpl::InnerIterator : public MatrixType::InnerIterator +template class TransposeImpl::InnerIterator + : public _MatrixTypeNested::InnerIterator { - typedef typename MatrixType::InnerIterator Base; + typedef typename _MatrixTypeNested::InnerIterator Base; public: EIGEN_STRONG_INLINE InnerIterator(const TransposeImpl& trans, Index outer) @@ -63,9 +52,10 @@ template class TransposeImpl::InnerItera inline Index col() const { return Base::row(); } }; -template class TransposeImpl::ReverseInnerIterator : public MatrixType::ReverseInnerIterator +template class TransposeImpl::ReverseInnerIterator + : public _MatrixTypeNested::ReverseInnerIterator { - typedef typename MatrixType::ReverseInnerIterator Base; + typedef typename _MatrixTypeNested::ReverseInnerIterator Base; public: EIGEN_STRONG_INLINE ReverseInnerIterator(const TransposeImpl& xpr, Index outer) diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h index 76a29e3aa..81941994f 100644 --- a/Eigen/src/Sparse/SparseUtil.h +++ b/Eigen/src/Sparse/SparseUtil.h @@ -92,15 +92,11 @@ template class SparseSelfAdjointView; template class SparseDiagonalProduct; template class SparseView; -template class SparseProduct; +template class SparseSparseProduct; template class SparseTimeDenseProduct; template class DenseTimeSparseProduct; -template::StorageKind, - typename RhsStorage = typename ei_traits::StorageKind> struct ei_sparse_product_mode; - -template struct SparseProductReturnType; +template struct SparseSparseProductReturnType; template struct ei_eval { diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h index 8bda8f395..dfe79a3b2 100644 --- a/Eigen/src/Sparse/SparseVector.h +++ b/Eigen/src/Sparse/SparseVector.h @@ -247,7 +247,7 @@ class SparseVector #ifndef EIGEN_PARSED_BY_DOXYGEN template - inline SparseVector& operator=(const SparseProduct& product) + inline SparseVector& operator=(const SparseSparseProduct& product) { return Base::operator=(product); } diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp index c07735e43..44239545c 100644 --- a/test/sparse_product.cpp +++ b/test/sparse_product.cpp @@ -36,6 +36,9 @@ template void sparse_product(const SparseMatrixType& typedef Matrix DenseMatrix; typedef Matrix DenseVector; + Scalar s1 = ei_random(); + Scalar s2 = ei_random(); + // test matrix-matrix product { DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); @@ -49,11 +52,16 @@ template void sparse_product(const SparseMatrixType& initSparse(density, refMat2, m2); initSparse(density, refMat3, m3); initSparse(density, refMat4, m4); + VERIFY_IS_APPROX(m4=m2*m3, refMat4=refMat2*refMat3); VERIFY_IS_APPROX(m4=m2.transpose()*m3, refMat4=refMat2.transpose()*refMat3); VERIFY_IS_APPROX(m4=m2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); VERIFY_IS_APPROX(m4=m2*m3.transpose(), refMat4=refMat2*refMat3.transpose()); + VERIFY_IS_APPROX(m4 = m2*m3/s1, refMat4 = refMat2*refMat3/s1); + VERIFY_IS_APPROX(m4 = m2*m3*s1, refMat4 = refMat2*refMat3*s1); + VERIFY_IS_APPROX(m4 = s2*m2*m3*s1, refMat4 = s2*refMat2*refMat3*s1); + // sparse * dense VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3); VERIFY_IS_APPROX(dm4=m2*refMat3.transpose(), refMat4=refMat2*refMat3.transpose()); @@ -62,7 +70,7 @@ template void sparse_product(const SparseMatrixType& VERIFY_IS_APPROX(dm4=m2*(refMat3+refMat3), refMat4=refMat2*(refMat3+refMat3)); VERIFY_IS_APPROX(dm4=m2.transpose()*(refMat3+refMat5)*0.5, refMat4=refMat2.transpose()*(refMat3+refMat5)*0.5); - + // dense * sparse VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3); VERIFY_IS_APPROX(dm4=refMat2*m3.transpose(), refMat4=refMat2*refMat3.transpose());