diff options
Diffstat (limited to 'unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h')
| -rw-r--r-- | unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h | 169 |
1 files changed, 115 insertions, 54 deletions
diff --git a/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h b/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h index 532896c3b..582fa8512 100644 --- a/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h +++ b/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h @@ -12,58 +12,93 @@ #ifndef KRONECKER_TENSOR_PRODUCT_H #define KRONECKER_TENSOR_PRODUCT_H -namespace Eigen { - -template<typename Scalar, int Options, typename Index> class SparseMatrix; +namespace Eigen { /*! - * \brief Kronecker tensor product helper class for dense matrices + * \ingroup KroneckerProduct_Module * - * This class is the return value of kroneckerProduct(MatrixBase, - * MatrixBase). Use the function rather than construct this class - * directly to avoid specifying template prarameters. + * \brief The base class of dense and sparse Kronecker product. * - * \tparam Lhs Type of the left-hand side, a matrix expression. - * \tparam Rhs Type of the rignt-hand side, a matrix expression. + * \tparam Derived is the derived type. */ -template<typename Lhs, typename Rhs> -class KroneckerProduct : public ReturnByValue<KroneckerProduct<Lhs,Rhs> > +template<typename Derived> +class KroneckerProductBase : public ReturnByValue<Derived> { private: - typedef ReturnByValue<KroneckerProduct> Base; - typedef typename Base::Scalar Scalar; - typedef typename Base::Index Index; + typedef typename internal::traits<Derived> Traits; + typedef typename Traits::Scalar Scalar; + + protected: + typedef typename Traits::Lhs Lhs; + typedef typename Traits::Rhs Rhs; public: /*! \brief Constructor. */ - KroneckerProduct(const Lhs& A, const Rhs& B) + KroneckerProductBase(const Lhs& A, const Rhs& B) : m_A(A), m_B(B) {} - /*! \brief Evaluate the Kronecker tensor product. */ - template<typename Dest> void evalTo(Dest& dst) const; - inline Index rows() const { return m_A.rows() * m_B.rows(); } inline Index cols() const { return m_A.cols() * m_B.cols(); } + /*! + * This overrides ReturnByValue::coeff because this function is + * efficient enough. + */ Scalar coeff(Index row, Index col) const { return m_A.coeff(row / m_B.rows(), col / m_B.cols()) * m_B.coeff(row % m_B.rows(), col % m_B.cols()); } + /*! + * This overrides ReturnByValue::coeff because this function is + * efficient enough. + */ Scalar coeff(Index i) const { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(KroneckerProduct); + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived); return m_A.coeff(i / m_A.size()) * m_B.coeff(i % m_A.size()); } - private: + protected: typename Lhs::Nested m_A; typename Rhs::Nested m_B; }; /*! + * \ingroup KroneckerProduct_Module + * + * \brief Kronecker tensor product helper class for dense matrices + * + * This class is the return value of kroneckerProduct(MatrixBase, + * MatrixBase). Use the function rather than construct this class + * directly to avoid specifying template prarameters. + * + * \tparam Lhs Type of the left-hand side, a matrix expression. + * \tparam Rhs Type of the rignt-hand side, a matrix expression. + */ +template<typename Lhs, typename Rhs> +class KroneckerProduct : public KroneckerProductBase<KroneckerProduct<Lhs,Rhs> > +{ + private: + typedef KroneckerProductBase<KroneckerProduct> Base; + using Base::m_A; + using Base::m_B; + + public: + /*! \brief Constructor. */ + KroneckerProduct(const Lhs& A, const Rhs& B) + : Base(A, B) + {} + + /*! \brief Evaluate the Kronecker tensor product. */ + template<typename Dest> void evalTo(Dest& dst) const; +}; + +/*! + * \ingroup KroneckerProduct_Module + * * \brief Kronecker tensor product helper class for sparse matrices * * If at least one of the operands is a sparse matrix expression, @@ -77,34 +112,21 @@ class KroneckerProduct : public ReturnByValue<KroneckerProduct<Lhs,Rhs> > * \tparam Rhs Type of the rignt-hand side, a matrix expression. */ template<typename Lhs, typename Rhs> -class KroneckerProductSparse : public EigenBase<KroneckerProductSparse<Lhs,Rhs> > +class KroneckerProductSparse : public KroneckerProductBase<KroneckerProductSparse<Lhs,Rhs> > { private: - typedef typename internal::traits<KroneckerProductSparse>::Index Index; + typedef KroneckerProductBase<KroneckerProductSparse> Base; + using Base::m_A; + using Base::m_B; public: /*! \brief Constructor. */ KroneckerProductSparse(const Lhs& A, const Rhs& B) - : m_A(A), m_B(B) + : Base(A, B) {} /*! \brief Evaluate the Kronecker tensor product. */ template<typename Dest> void evalTo(Dest& dst) const; - - inline Index rows() const { return m_A.rows() * m_B.rows(); } - inline Index cols() const { return m_A.cols() * m_B.cols(); } - - template<typename Scalar, int Options, typename Index> - operator SparseMatrix<Scalar, Options, Index>() - { - SparseMatrix<Scalar, Options, Index> result; - evalTo(result.derived()); - return result; - } - - private: - typename Lhs::Nested m_A; - typename Rhs::Nested m_B; }; template<typename Lhs, typename Rhs> @@ -124,22 +146,49 @@ template<typename Lhs, typename Rhs> template<typename Dest> void KroneckerProductSparse<Lhs,Rhs>::evalTo(Dest& dst) const { - const Index Br = m_B.rows(), - Bc = m_B.cols(); - dst.resize(rows(),cols()); + Index Br = m_B.rows(), Bc = m_B.cols(); + dst.resize(this->rows(), this->cols()); dst.resizeNonZeros(0); - dst.reserve(m_A.nonZeros() * m_B.nonZeros()); + + // 1 - evaluate the operands if needed: + typedef typename internal::nested_eval<Lhs,Dynamic>::type Lhs1; + typedef typename internal::remove_all<Lhs1>::type Lhs1Cleaned; + const Lhs1 lhs1(m_A); + typedef typename internal::nested_eval<Rhs,Dynamic>::type Rhs1; + typedef typename internal::remove_all<Rhs1>::type Rhs1Cleaned; + const Rhs1 rhs1(m_B); + + // 2 - construct respective iterators + typedef Eigen::InnerIterator<Lhs1Cleaned> LhsInnerIterator; + typedef Eigen::InnerIterator<Rhs1Cleaned> RhsInnerIterator; + + // compute number of non-zeros per innervectors of dst + { + // TODO VectorXi is not necessarily big enough! + VectorXi nnzA = VectorXi::Zero(Dest::IsRowMajor ? m_A.rows() : m_A.cols()); + for (Index kA=0; kA < m_A.outerSize(); ++kA) + for (LhsInnerIterator itA(lhs1,kA); itA; ++itA) + nnzA(Dest::IsRowMajor ? itA.row() : itA.col())++; + + VectorXi nnzB = VectorXi::Zero(Dest::IsRowMajor ? m_B.rows() : m_B.cols()); + for (Index kB=0; kB < m_B.outerSize(); ++kB) + for (RhsInnerIterator itB(rhs1,kB); itB; ++itB) + nnzB(Dest::IsRowMajor ? itB.row() : itB.col())++; + + Matrix<int,Dynamic,Dynamic,ColMajor> nnzAB = nnzB * nnzA.transpose(); + dst.reserve(VectorXi::Map(nnzAB.data(), nnzAB.size())); + } for (Index kA=0; kA < m_A.outerSize(); ++kA) { for (Index kB=0; kB < m_B.outerSize(); ++kB) { - for (typename Lhs::InnerIterator itA(m_A,kA); itA; ++itA) + for (LhsInnerIterator itA(lhs1,kA); itA; ++itA) { - for (typename Rhs::InnerIterator itB(m_B,kB); itB; ++itB) + for (RhsInnerIterator itB(rhs1,kB); itB; ++itB) { - const Index i = itA.row() * Br + itB.row(), - j = itA.col() * Bc + itB.col(); + Index i = itA.row() * Br + itB.row(), + j = itA.col() * Bc + itB.col(); dst.insert(i,j) = itA.value() * itB.value(); } } @@ -154,14 +203,14 @@ struct traits<KroneckerProduct<_Lhs,_Rhs> > { typedef typename remove_all<_Lhs>::type Lhs; typedef typename remove_all<_Rhs>::type Rhs; - typedef typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar; + typedef typename ScalarBinaryOpTraits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar; + typedef typename promote_index_type<typename Lhs::StorageIndex, typename Rhs::StorageIndex>::type StorageIndex; enum { Rows = size_at_compile_time<traits<Lhs>::RowsAtCompileTime, traits<Rhs>::RowsAtCompileTime>::ret, Cols = size_at_compile_time<traits<Lhs>::ColsAtCompileTime, traits<Rhs>::ColsAtCompileTime>::ret, MaxRows = size_at_compile_time<traits<Lhs>::MaxRowsAtCompileTime, traits<Rhs>::MaxRowsAtCompileTime>::ret, - MaxCols = size_at_compile_time<traits<Lhs>::MaxColsAtCompileTime, traits<Rhs>::MaxColsAtCompileTime>::ret, - CoeffReadCost = Lhs::CoeffReadCost + Rhs::CoeffReadCost + NumTraits<Scalar>::MulCost + MaxCols = size_at_compile_time<traits<Lhs>::MaxColsAtCompileTime, traits<Rhs>::MaxColsAtCompileTime>::ret }; typedef Matrix<Scalar,Rows,Cols> ReturnType; @@ -173,9 +222,9 @@ struct traits<KroneckerProductSparse<_Lhs,_Rhs> > typedef MatrixXpr XprKind; typedef typename remove_all<_Lhs>::type Lhs; typedef typename remove_all<_Rhs>::type Rhs; - typedef typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar; - typedef typename promote_storage_type<typename traits<Lhs>::StorageKind, typename traits<Rhs>::StorageKind>::ret StorageKind; - typedef typename promote_index_type<typename Lhs::Index, typename Rhs::Index>::type Index; + typedef typename ScalarBinaryOpTraits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar; + typedef typename cwise_promote_storage_type<typename traits<Lhs>::StorageKind, typename traits<Rhs>::StorageKind, scalar_product_op<typename Lhs::Scalar, typename Rhs::Scalar> >::ret StorageKind; + typedef typename promote_index_type<typename Lhs::StorageIndex, typename Rhs::StorageIndex>::type StorageIndex; enum { LhsFlags = Lhs::Flags, @@ -190,9 +239,11 @@ struct traits<KroneckerProductSparse<_Lhs,_Rhs> > RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit), Flags = ((LhsFlags | RhsFlags) & HereditaryBits & RemovedBits) - | EvalBeforeNestingBit | EvalBeforeAssigningBit, - CoeffReadCost = Dynamic + | EvalBeforeNestingBit, + CoeffReadCost = HugeCost }; + + typedef SparseMatrix<Scalar, 0, StorageIndex> ReturnType; }; } // end namespace internal @@ -228,6 +279,16 @@ KroneckerProduct<A,B> kroneckerProduct(const MatrixBase<A>& a, const MatrixBase< * Computes Kronecker tensor product of two matrices, at least one of * which is sparse * + * \warning If you want to replace a matrix by its Kronecker product + * with some matrix, do \b NOT do this: + * \code + * A = kroneckerProduct(A,B); // bug!!! caused by aliasing effect + * \endcode + * instead, use eval() to work around this: + * \code + * A = kroneckerProduct(A,B).eval(); + * \endcode + * * \param a Dense/sparse matrix a * \param b Dense/sparse matrix b * \return Kronecker tensor product of a and b, stored in a sparse |