src/KroneckerProduct/KroneckerTensorProduct.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2011 Kolja Brix <brix@igpm.rwth-aachen.de>
0005 // Copyright (C) 2011 Andreas Platen <andiplaten@gmx.de>
0006 // Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
0007 //
0008 // This Source Code Form is subject to the terms of the Mozilla
0009 // Public License v. 2.0. If a copy of the MPL was not distributed
0010 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0011
0012 #ifndef KRONECKER_TENSOR_PRODUCT_H
0013 #define KRONECKER_TENSOR_PRODUCT_H
0014
0015 namespace Eigen {
0016
0017 /*!
0018  * \ingroup KroneckerProduct_Module
0019  *
0020  * \brief The base class of dense and sparse Kronecker product.
0021  *
0022  * \tparam Derived is the derived type.
0023  */
0024 template<typename Derived>
0025 class KroneckerProductBase : public ReturnByValue<Derived>
0026 {
0027   private:
0028     typedef typename internal::traits<Derived> Traits;
0029     typedef typename Traits::Scalar Scalar;
0030
0031   protected:
0032     typedef typename Traits::Lhs Lhs;
0033     typedef typename Traits::Rhs Rhs;
0034
0035   public:
0036     /*! \brief Constructor. */
0037     KroneckerProductBase(const Lhs& A, const Rhs& B)
0038       : m_A(A), m_B(B)
0039     {}
0040
0041     inline Index rows() const { return m_A.rows() * m_B.rows(); }
0042     inline Index cols() const { return m_A.cols() * m_B.cols(); }
0043
0044     /*!
0045      * This overrides ReturnByValue::coeff because this function is
0046      * efficient enough.
0047      */
0048     Scalar coeff(Index row, Index col) const
0049     {
0050       return m_A.coeff(row / m_B.rows(), col / m_B.cols()) *
0051              m_B.coeff(row % m_B.rows(), col % m_B.cols());
0052     }
0053
0054     /*!
0055      * This overrides ReturnByValue::coeff because this function is
0056      * efficient enough.
0057      */
0058     Scalar coeff(Index i) const
0059     {
0060       EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived);
0061       return m_A.coeff(i / m_A.size()) * m_B.coeff(i % m_A.size());
0062     }
0063
0064   protected:
0065     typename Lhs::Nested m_A;
0066     typename Rhs::Nested m_B;
0067 };
0068
0069 /*!
0070  * \ingroup KroneckerProduct_Module
0071  *
0072  * \brief Kronecker tensor product helper class for dense matrices
0073  *
0074  * This class is the return value of kroneckerProduct(MatrixBase,
0075  * MatrixBase). Use the function rather than construct this class
0076  * directly to avoid specifying template prarameters.
0077  *
0078  * \tparam Lhs  Type of the left-hand side, a matrix expression.
0079  * \tparam Rhs  Type of the rignt-hand side, a matrix expression.
0080  */
0081 template<typename Lhs, typename Rhs>
0082 class KroneckerProduct : public KroneckerProductBase<KroneckerProduct<Lhs,Rhs> >
0083 {
0084   private:
0085     typedef KroneckerProductBase<KroneckerProduct> Base;
0086     using Base::m_A;
0087     using Base::m_B;
0088
0089   public:
0090     /*! \brief Constructor. */
0091     KroneckerProduct(const Lhs& A, const Rhs& B)
0092       : Base(A, B)
0093     {}
0094
0095     /*! \brief Evaluate the Kronecker tensor product. */
0096     template<typename Dest> void evalTo(Dest& dst) const;
0097 };
0098
0099 /*!
0100  * \ingroup KroneckerProduct_Module
0101  *
0102  * \brief Kronecker tensor product helper class for sparse matrices
0103  *
0104  * If at least one of the operands is a sparse matrix expression,
0105  * then this class is returned and evaluates into a sparse matrix.
0106  *
0107  * This class is the return value of kroneckerProduct(EigenBase,
0108  * EigenBase). Use the function rather than construct this class
0109  * directly to avoid specifying template prarameters.
0110  *
0111  * \tparam Lhs  Type of the left-hand side, a matrix expression.
0112  * \tparam Rhs  Type of the rignt-hand side, a matrix expression.
0113  */
0114 template<typename Lhs, typename Rhs>
0115 class KroneckerProductSparse : public KroneckerProductBase<KroneckerProductSparse<Lhs,Rhs> >
0116 {
0117   private:
0118     typedef KroneckerProductBase<KroneckerProductSparse> Base;
0119     using Base::m_A;
0120     using Base::m_B;
0121
0122   public:
0123     /*! \brief Constructor. */
0124     KroneckerProductSparse(const Lhs& A, const Rhs& B)
0125       : Base(A, B)
0126     {}
0127
0128     /*! \brief Evaluate the Kronecker tensor product. */
0129     template<typename Dest> void evalTo(Dest& dst) const;
0130 };
0131
0132 template<typename Lhs, typename Rhs>
0133 template<typename Dest>
0134 void KroneckerProduct<Lhs,Rhs>::evalTo(Dest& dst) const
0135 {
0136   const int BlockRows = Rhs::RowsAtCompileTime,
0137             BlockCols = Rhs::ColsAtCompileTime;
0138   const Index Br = m_B.rows(),
0139               Bc = m_B.cols();
0140   for (Index i=0; i < m_A.rows(); ++i)
0141     for (Index j=0; j < m_A.cols(); ++j)
0142       Block<Dest,BlockRows,BlockCols>(dst,i*Br,j*Bc,Br,Bc) = m_A.coeff(i,j) * m_B;
0143 }
0144
0145 template<typename Lhs, typename Rhs>
0146 template<typename Dest>
0147 void KroneckerProductSparse<Lhs,Rhs>::evalTo(Dest& dst) const
0148 {
0149   Index Br = m_B.rows(), Bc = m_B.cols();
0150   dst.resize(this->rows(), this->cols());
0151   dst.resizeNonZeros(0);
0152
0153   // 1 - evaluate the operands if needed:
0154   typedef typename internal::nested_eval<Lhs,Dynamic>::type Lhs1;
0155   typedef typename internal::remove_all<Lhs1>::type Lhs1Cleaned;
0156   const Lhs1 lhs1(m_A);
0157   typedef typename internal::nested_eval<Rhs,Dynamic>::type Rhs1;
0158   typedef typename internal::remove_all<Rhs1>::type Rhs1Cleaned;
0159   const Rhs1 rhs1(m_B);
0160
0161   // 2 - construct respective iterators
0162   typedef Eigen::InnerIterator<Lhs1Cleaned> LhsInnerIterator;
0163   typedef Eigen::InnerIterator<Rhs1Cleaned> RhsInnerIterator;
0164
0165   // compute number of non-zeros per innervectors of dst
0166   {
0167     // TODO VectorXi is not necessarily big enough!
0168     VectorXi nnzA = VectorXi::Zero(Dest::IsRowMajor ? m_A.rows() : m_A.cols());
0169     for (Index kA=0; kA < m_A.outerSize(); ++kA)
0170       for (LhsInnerIterator itA(lhs1,kA); itA; ++itA)
0171         nnzA(Dest::IsRowMajor ? itA.row() : itA.col())++;
0172
0173     VectorXi nnzB = VectorXi::Zero(Dest::IsRowMajor ? m_B.rows() : m_B.cols());
0174     for (Index kB=0; kB < m_B.outerSize(); ++kB)
0175       for (RhsInnerIterator itB(rhs1,kB); itB; ++itB)
0176         nnzB(Dest::IsRowMajor ? itB.row() : itB.col())++;
0177
0178     Matrix<int,Dynamic,Dynamic,ColMajor> nnzAB = nnzB * nnzA.transpose();
0179     dst.reserve(VectorXi::Map(nnzAB.data(), nnzAB.size()));
0180   }
0181
0182   for (Index kA=0; kA < m_A.outerSize(); ++kA)
0183   {
0184     for (Index kB=0; kB < m_B.outerSize(); ++kB)
0185     {
0186       for (LhsInnerIterator itA(lhs1,kA); itA; ++itA)
0187       {
0188         for (RhsInnerIterator itB(rhs1,kB); itB; ++itB)
0189         {
0190           Index i = itA.row() * Br + itB.row(),
0191                 j = itA.col() * Bc + itB.col();
0192           dst.insert(i,j) = itA.value() * itB.value();
0193         }
0194       }
0195     }
0196   }
0197 }
0198
0199 namespace internal {
0200
0201 template<typename _Lhs, typename _Rhs>
0202 struct traits<KroneckerProduct<_Lhs,_Rhs> >
0203 {
0204   typedef typename remove_all<_Lhs>::type Lhs;
0205   typedef typename remove_all<_Rhs>::type Rhs;
0206   typedef typename ScalarBinaryOpTraits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar;
0207   typedef typename promote_index_type<typename Lhs::StorageIndex, typename Rhs::StorageIndex>::type StorageIndex;
0208
0209   enum {
0210     Rows = size_at_compile_time<traits<Lhs>::RowsAtCompileTime, traits<Rhs>::RowsAtCompileTime>::ret,
0211     Cols = size_at_compile_time<traits<Lhs>::ColsAtCompileTime, traits<Rhs>::ColsAtCompileTime>::ret,
0212     MaxRows = size_at_compile_time<traits<Lhs>::MaxRowsAtCompileTime, traits<Rhs>::MaxRowsAtCompileTime>::ret,
0213     MaxCols = size_at_compile_time<traits<Lhs>::MaxColsAtCompileTime, traits<Rhs>::MaxColsAtCompileTime>::ret
0214   };
0215
0216   typedef Matrix<Scalar,Rows,Cols> ReturnType;
0217 };
0218
0219 template<typename _Lhs, typename _Rhs>
0220 struct traits<KroneckerProductSparse<_Lhs,_Rhs> >
0221 {
0222   typedef MatrixXpr XprKind;
0223   typedef typename remove_all<_Lhs>::type Lhs;
0224   typedef typename remove_all<_Rhs>::type Rhs;
0225   typedef typename ScalarBinaryOpTraits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar;
0226   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;
0227   typedef typename promote_index_type<typename Lhs::StorageIndex, typename Rhs::StorageIndex>::type StorageIndex;
0228
0229   enum {
0230     LhsFlags = Lhs::Flags,
0231     RhsFlags = Rhs::Flags,
0232
0233     RowsAtCompileTime = size_at_compile_time<traits<Lhs>::RowsAtCompileTime, traits<Rhs>::RowsAtCompileTime>::ret,
0234     ColsAtCompileTime = size_at_compile_time<traits<Lhs>::ColsAtCompileTime, traits<Rhs>::ColsAtCompileTime>::ret,
0235     MaxRowsAtCompileTime = size_at_compile_time<traits<Lhs>::MaxRowsAtCompileTime, traits<Rhs>::MaxRowsAtCompileTime>::ret,
0236     MaxColsAtCompileTime = size_at_compile_time<traits<Lhs>::MaxColsAtCompileTime, traits<Rhs>::MaxColsAtCompileTime>::ret,
0237
0238     EvalToRowMajor = (int(LhsFlags) & int(RhsFlags) & RowMajorBit),
0239     RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit),
0240
0241     Flags = ((int(LhsFlags) | int(RhsFlags)) & HereditaryBits & RemovedBits)
0242           | EvalBeforeNestingBit,
0243     CoeffReadCost = HugeCost
0244   };
0245
0246   typedef SparseMatrix<Scalar, 0, StorageIndex> ReturnType;
0247 };
0248
0249 } // end namespace internal
0250
0251 /*!
0252  * \ingroup KroneckerProduct_Module
0253  *
0254  * Computes Kronecker tensor product of two dense matrices
0255  *
0256  * \warning If you want to replace a matrix by its Kronecker product
0257  *          with some matrix, do \b NOT do this:
0258  * \code
0259  * A = kroneckerProduct(A,B); // bug!!! caused by aliasing effect
0260  * \endcode
0261  * instead, use eval() to work around this:
0262  * \code
0263  * A = kroneckerProduct(A,B).eval();
0264  * \endcode
0265  *
0266  * \param a  Dense matrix a
0267  * \param b  Dense matrix b
0268  * \return   Kronecker tensor product of a and b
0269  */
0270 template<typename A, typename B>
0271 KroneckerProduct<A,B> kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b)
0272 {
0273   return KroneckerProduct<A, B>(a.derived(), b.derived());
0274 }
0275
0276 /*!
0277  * \ingroup KroneckerProduct_Module
0278  *
0279  * Computes Kronecker tensor product of two matrices, at least one of
0280  * which is sparse
0281  *
0282  * \warning If you want to replace a matrix by its Kronecker product
0283  *          with some matrix, do \b NOT do this:
0284  * \code
0285  * A = kroneckerProduct(A,B); // bug!!! caused by aliasing effect
0286  * \endcode
0287  * instead, use eval() to work around this:
0288  * \code
0289  * A = kroneckerProduct(A,B).eval();
0290  * \endcode
0291  *
0292  * \param a  Dense/sparse matrix a
0293  * \param b  Dense/sparse matrix b
0294  * \return   Kronecker tensor product of a and b, stored in a sparse
0295  *           matrix
0296  */
0297 template<typename A, typename B>
0298 KroneckerProductSparse<A,B> kroneckerProduct(const EigenBase<A>& a, const EigenBase<B>& b)
0299 {
0300   return KroneckerProductSparse<A,B>(a.derived(), b.derived());
0301 }
0302
0303 } // end namespace Eigen
0304
0305 #endif // KRONECKER_TENSOR_PRODUCT_H