1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2011 Kolja Brix <brix (at) igpm.rwth-aachen.de> 5 // Copyright (C) 2011 Andreas Platen <andiplaten (at) gmx.de> 6 // Copyright (C) 2012 Chen-Pang He <jdh8 (at) ms63.hinet.net> 7 // 8 // This Source Code Form is subject to the terms of the Mozilla 9 // Public License v. 2.0. If a copy of the MPL was not distributed 10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 11 12 #ifndef KRONECKER_TENSOR_PRODUCT_H 13 #define KRONECKER_TENSOR_PRODUCT_H 14 15 namespace Eigen { 16 17 /*! 18 * \ingroup KroneckerProduct_Module 19 * 20 * \brief The base class of dense and sparse Kronecker product. 21 * 22 * \tparam Derived is the derived type. 23 */ 24 template<typename Derived> 25 class KroneckerProductBase : public ReturnByValue<Derived> 26 { 27 private: 28 typedef typename internal::traits<Derived> Traits; 29 typedef typename Traits::Scalar Scalar; 30 31 protected: 32 typedef typename Traits::Lhs Lhs; 33 typedef typename Traits::Rhs Rhs; 34 35 public: 36 /*! \brief Constructor. */ 37 KroneckerProductBase(const Lhs& A, const Rhs& B) 38 : m_A(A), m_B(B) 39 {} 40 41 inline Index rows() const { return m_A.rows() * m_B.rows(); } 42 inline Index cols() const { return m_A.cols() * m_B.cols(); } 43 44 /*! 45 * This overrides ReturnByValue::coeff because this function is 46 * efficient enough. 47 */ 48 Scalar coeff(Index row, Index col) const 49 { 50 return m_A.coeff(row / m_B.rows(), col / m_B.cols()) * 51 m_B.coeff(row % m_B.rows(), col % m_B.cols()); 52 } 53 54 /*! 55 * This overrides ReturnByValue::coeff because this function is 56 * efficient enough. 57 */ 58 Scalar coeff(Index i) const 59 { 60 EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived); 61 return m_A.coeff(i / m_A.size()) * m_B.coeff(i % m_A.size()); 62 } 63 64 protected: 65 typename Lhs::Nested m_A; 66 typename Rhs::Nested m_B; 67 }; 68 69 /*! 70 * \ingroup KroneckerProduct_Module 71 * 72 * \brief Kronecker tensor product helper class for dense matrices 73 * 74 * This class is the return value of kroneckerProduct(MatrixBase, 75 * MatrixBase). Use the function rather than construct this class 76 * directly to avoid specifying template prarameters. 77 * 78 * \tparam Lhs Type of the left-hand side, a matrix expression. 79 * \tparam Rhs Type of the rignt-hand side, a matrix expression. 80 */ 81 template<typename Lhs, typename Rhs> 82 class KroneckerProduct : public KroneckerProductBase<KroneckerProduct<Lhs,Rhs> > 83 { 84 private: 85 typedef KroneckerProductBase<KroneckerProduct> Base; 86 using Base::m_A; 87 using Base::m_B; 88 89 public: 90 /*! \brief Constructor. */ 91 KroneckerProduct(const Lhs& A, const Rhs& B) 92 : Base(A, B) 93 {} 94 95 /*! \brief Evaluate the Kronecker tensor product. */ 96 template<typename Dest> void evalTo(Dest& dst) const; 97 }; 98 99 /*! 100 * \ingroup KroneckerProduct_Module 101 * 102 * \brief Kronecker tensor product helper class for sparse matrices 103 * 104 * If at least one of the operands is a sparse matrix expression, 105 * then this class is returned and evaluates into a sparse matrix. 106 * 107 * This class is the return value of kroneckerProduct(EigenBase, 108 * EigenBase). Use the function rather than construct this class 109 * directly to avoid specifying template prarameters. 110 * 111 * \tparam Lhs Type of the left-hand side, a matrix expression. 112 * \tparam Rhs Type of the rignt-hand side, a matrix expression. 113 */ 114 template<typename Lhs, typename Rhs> 115 class KroneckerProductSparse : public KroneckerProductBase<KroneckerProductSparse<Lhs,Rhs> > 116 { 117 private: 118 typedef KroneckerProductBase<KroneckerProductSparse> Base; 119 using Base::m_A; 120 using Base::m_B; 121 122 public: 123 /*! \brief Constructor. */ 124 KroneckerProductSparse(const Lhs& A, const Rhs& B) 125 : Base(A, B) 126 {} 127 128 /*! \brief Evaluate the Kronecker tensor product. */ 129 template<typename Dest> void evalTo(Dest& dst) const; 130 }; 131 132 template<typename Lhs, typename Rhs> 133 template<typename Dest> 134 void KroneckerProduct<Lhs,Rhs>::evalTo(Dest& dst) const 135 { 136 const int BlockRows = Rhs::RowsAtCompileTime, 137 BlockCols = Rhs::ColsAtCompileTime; 138 const Index Br = m_B.rows(), 139 Bc = m_B.cols(); 140 for (Index i=0; i < m_A.rows(); ++i) 141 for (Index j=0; j < m_A.cols(); ++j) 142 Block<Dest,BlockRows,BlockCols>(dst,i*Br,j*Bc,Br,Bc) = m_A.coeff(i,j) * m_B; 143 } 144 145 template<typename Lhs, typename Rhs> 146 template<typename Dest> 147 void KroneckerProductSparse<Lhs,Rhs>::evalTo(Dest& dst) const 148 { 149 Index Br = m_B.rows(), Bc = m_B.cols(); 150 dst.resize(this->rows(), this->cols()); 151 dst.resizeNonZeros(0); 152 153 // 1 - evaluate the operands if needed: 154 typedef typename internal::nested_eval<Lhs,Dynamic>::type Lhs1; 155 typedef typename internal::remove_all<Lhs1>::type Lhs1Cleaned; 156 const Lhs1 lhs1(m_A); 157 typedef typename internal::nested_eval<Rhs,Dynamic>::type Rhs1; 158 typedef typename internal::remove_all<Rhs1>::type Rhs1Cleaned; 159 const Rhs1 rhs1(m_B); 160 161 // 2 - construct respective iterators 162 typedef Eigen::InnerIterator<Lhs1Cleaned> LhsInnerIterator; 163 typedef Eigen::InnerIterator<Rhs1Cleaned> RhsInnerIterator; 164 165 // compute number of non-zeros per innervectors of dst 166 { 167 // TODO VectorXi is not necessarily big enough! 168 VectorXi nnzA = VectorXi::Zero(Dest::IsRowMajor ? m_A.rows() : m_A.cols()); 169 for (Index kA=0; kA < m_A.outerSize(); ++kA) 170 for (LhsInnerIterator itA(lhs1,kA); itA; ++itA) 171 nnzA(Dest::IsRowMajor ? itA.row() : itA.col())++; 172 173 VectorXi nnzB = VectorXi::Zero(Dest::IsRowMajor ? m_B.rows() : m_B.cols()); 174 for (Index kB=0; kB < m_B.outerSize(); ++kB) 175 for (RhsInnerIterator itB(rhs1,kB); itB; ++itB) 176 nnzB(Dest::IsRowMajor ? itB.row() : itB.col())++; 177 178 Matrix<int,Dynamic,Dynamic,ColMajor> nnzAB = nnzB * nnzA.transpose(); 179 dst.reserve(VectorXi::Map(nnzAB.data(), nnzAB.size())); 180 } 181 182 for (Index kA=0; kA < m_A.outerSize(); ++kA) 183 { 184 for (Index kB=0; kB < m_B.outerSize(); ++kB) 185 { 186 for (LhsInnerIterator itA(lhs1,kA); itA; ++itA) 187 { 188 for (RhsInnerIterator itB(rhs1,kB); itB; ++itB) 189 { 190 Index i = itA.row() * Br + itB.row(), 191 j = itA.col() * Bc + itB.col(); 192 dst.insert(i,j) = itA.value() * itB.value(); 193 } 194 } 195 } 196 } 197 } 198 199 namespace internal { 200 201 template<typename _Lhs, typename _Rhs> 202 struct traits<KroneckerProduct<_Lhs,_Rhs> > 203 { 204 typedef typename remove_all<_Lhs>::type Lhs; 205 typedef typename remove_all<_Rhs>::type Rhs; 206 typedef typename ScalarBinaryOpTraits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar; 207 typedef typename promote_index_type<typename Lhs::StorageIndex, typename Rhs::StorageIndex>::type StorageIndex; 208 209 enum { 210 Rows = size_at_compile_time<traits<Lhs>::RowsAtCompileTime, traits<Rhs>::RowsAtCompileTime>::ret, 211 Cols = size_at_compile_time<traits<Lhs>::ColsAtCompileTime, traits<Rhs>::ColsAtCompileTime>::ret, 212 MaxRows = size_at_compile_time<traits<Lhs>::MaxRowsAtCompileTime, traits<Rhs>::MaxRowsAtCompileTime>::ret, 213 MaxCols = size_at_compile_time<traits<Lhs>::MaxColsAtCompileTime, traits<Rhs>::MaxColsAtCompileTime>::ret 214 }; 215 216 typedef Matrix<Scalar,Rows,Cols> ReturnType; 217 }; 218 219 template<typename _Lhs, typename _Rhs> 220 struct traits<KroneckerProductSparse<_Lhs,_Rhs> > 221 { 222 typedef MatrixXpr XprKind; 223 typedef typename remove_all<_Lhs>::type Lhs; 224 typedef typename remove_all<_Rhs>::type Rhs; 225 typedef typename ScalarBinaryOpTraits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar; 226 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; 227 typedef typename promote_index_type<typename Lhs::StorageIndex, typename Rhs::StorageIndex>::type StorageIndex; 228 229 enum { 230 LhsFlags = Lhs::Flags, 231 RhsFlags = Rhs::Flags, 232 233 RowsAtCompileTime = size_at_compile_time<traits<Lhs>::RowsAtCompileTime, traits<Rhs>::RowsAtCompileTime>::ret, 234 ColsAtCompileTime = size_at_compile_time<traits<Lhs>::ColsAtCompileTime, traits<Rhs>::ColsAtCompileTime>::ret, 235 MaxRowsAtCompileTime = size_at_compile_time<traits<Lhs>::MaxRowsAtCompileTime, traits<Rhs>::MaxRowsAtCompileTime>::ret, 236 MaxColsAtCompileTime = size_at_compile_time<traits<Lhs>::MaxColsAtCompileTime, traits<Rhs>::MaxColsAtCompileTime>::ret, 237 238 EvalToRowMajor = (LhsFlags & RhsFlags & RowMajorBit), 239 RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit), 240 241 Flags = ((LhsFlags | RhsFlags) & HereditaryBits & RemovedBits) 242 | EvalBeforeNestingBit, 243 CoeffReadCost = HugeCost 244 }; 245 246 typedef SparseMatrix<Scalar, 0, StorageIndex> ReturnType; 247 }; 248 249 } // end namespace internal 250 251 /*! 252 * \ingroup KroneckerProduct_Module 253 * 254 * Computes Kronecker tensor product of two dense matrices 255 * 256 * \warning If you want to replace a matrix by its Kronecker product 257 * with some matrix, do \b NOT do this: 258 * \code 259 * A = kroneckerProduct(A,B); // bug!!! caused by aliasing effect 260 * \endcode 261 * instead, use eval() to work around this: 262 * \code 263 * A = kroneckerProduct(A,B).eval(); 264 * \endcode 265 * 266 * \param a Dense matrix a 267 * \param b Dense matrix b 268 * \return Kronecker tensor product of a and b 269 */ 270 template<typename A, typename B> 271 KroneckerProduct<A,B> kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b) 272 { 273 return KroneckerProduct<A, B>(a.derived(), b.derived()); 274 } 275 276 /*! 277 * \ingroup KroneckerProduct_Module 278 * 279 * Computes Kronecker tensor product of two matrices, at least one of 280 * which is sparse 281 * 282 * \warning If you want to replace a matrix by its Kronecker product 283 * with some matrix, do \b NOT do this: 284 * \code 285 * A = kroneckerProduct(A,B); // bug!!! caused by aliasing effect 286 * \endcode 287 * instead, use eval() to work around this: 288 * \code 289 * A = kroneckerProduct(A,B).eval(); 290 * \endcode 291 * 292 * \param a Dense/sparse matrix a 293 * \param b Dense/sparse matrix b 294 * \return Kronecker tensor product of a and b, stored in a sparse 295 * matrix 296 */ 297 template<typename A, typename B> 298 KroneckerProductSparse<A,B> kroneckerProduct(const EigenBase<A>& a, const EigenBase<B>& b) 299 { 300 return KroneckerProductSparse<A,B>(a.derived(), b.derived()); 301 } 302 303 } // end namespace Eigen 304 305 #endif // KRONECKER_TENSOR_PRODUCT_H 306