1// This file is part of Eigen, a lightweight C++ template library 2// for linear algebra. 3// 4// Copyright (C) 2011 Kolja Brix <brix@igpm.rwth-aachen.de> 5// Copyright (C) 2011 Andreas Platen <andiplaten@gmx.de> 6// Copyright (C) 2012 Chen-Pang He <jdh8@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 15namespace Eigen { 16 17template<typename Scalar, int Options, typename Index> class SparseMatrix; 18 19/*! 20 * \brief Kronecker tensor product helper class for dense matrices 21 * 22 * This class is the return value of kroneckerProduct(MatrixBase, 23 * MatrixBase). Use the function rather than construct this class 24 * directly to avoid specifying template prarameters. 25 * 26 * \tparam Lhs Type of the left-hand side, a matrix expression. 27 * \tparam Rhs Type of the rignt-hand side, a matrix expression. 28 */ 29template<typename Lhs, typename Rhs> 30class KroneckerProduct : public ReturnByValue<KroneckerProduct<Lhs,Rhs> > 31{ 32 private: 33 typedef ReturnByValue<KroneckerProduct> Base; 34 typedef typename Base::Scalar Scalar; 35 typedef typename Base::Index Index; 36 37 public: 38 /*! \brief Constructor. */ 39 KroneckerProduct(const Lhs& A, const Rhs& B) 40 : m_A(A), m_B(B) 41 {} 42 43 /*! \brief Evaluate the Kronecker tensor product. */ 44 template<typename Dest> void evalTo(Dest& dst) const; 45 46 inline Index rows() const { return m_A.rows() * m_B.rows(); } 47 inline Index cols() const { return m_A.cols() * m_B.cols(); } 48 49 Scalar coeff(Index row, Index col) const 50 { 51 return m_A.coeff(row / m_B.rows(), col / m_B.cols()) * 52 m_B.coeff(row % m_B.rows(), col % m_B.cols()); 53 } 54 55 Scalar coeff(Index i) const 56 { 57 EIGEN_STATIC_ASSERT_VECTOR_ONLY(KroneckerProduct); 58 return m_A.coeff(i / m_A.size()) * m_B.coeff(i % m_A.size()); 59 } 60 61 private: 62 typename Lhs::Nested m_A; 63 typename Rhs::Nested m_B; 64}; 65 66/*! 67 * \brief Kronecker tensor product helper class for sparse matrices 68 * 69 * If at least one of the operands is a sparse matrix expression, 70 * then this class is returned and evaluates into a sparse matrix. 71 * 72 * This class is the return value of kroneckerProduct(EigenBase, 73 * EigenBase). Use the function rather than construct this class 74 * directly to avoid specifying template prarameters. 75 * 76 * \tparam Lhs Type of the left-hand side, a matrix expression. 77 * \tparam Rhs Type of the rignt-hand side, a matrix expression. 78 */ 79template<typename Lhs, typename Rhs> 80class KroneckerProductSparse : public EigenBase<KroneckerProductSparse<Lhs,Rhs> > 81{ 82 private: 83 typedef typename internal::traits<KroneckerProductSparse>::Index Index; 84 85 public: 86 /*! \brief Constructor. */ 87 KroneckerProductSparse(const Lhs& A, const Rhs& B) 88 : m_A(A), m_B(B) 89 {} 90 91 /*! \brief Evaluate the Kronecker tensor product. */ 92 template<typename Dest> void evalTo(Dest& dst) const; 93 94 inline Index rows() const { return m_A.rows() * m_B.rows(); } 95 inline Index cols() const { return m_A.cols() * m_B.cols(); } 96 97 template<typename Scalar, int Options, typename Index> 98 operator SparseMatrix<Scalar, Options, Index>() 99 { 100 SparseMatrix<Scalar, Options, Index> result; 101 evalTo(result.derived()); 102 return result; 103 } 104 105 private: 106 typename Lhs::Nested m_A; 107 typename Rhs::Nested m_B; 108}; 109 110template<typename Lhs, typename Rhs> 111template<typename Dest> 112void KroneckerProduct<Lhs,Rhs>::evalTo(Dest& dst) const 113{ 114 const int BlockRows = Rhs::RowsAtCompileTime, 115 BlockCols = Rhs::ColsAtCompileTime; 116 const Index Br = m_B.rows(), 117 Bc = m_B.cols(); 118 for (Index i=0; i < m_A.rows(); ++i) 119 for (Index j=0; j < m_A.cols(); ++j) 120 Block<Dest,BlockRows,BlockCols>(dst,i*Br,j*Bc,Br,Bc) = m_A.coeff(i,j) * m_B; 121} 122 123template<typename Lhs, typename Rhs> 124template<typename Dest> 125void KroneckerProductSparse<Lhs,Rhs>::evalTo(Dest& dst) const 126{ 127 const Index Br = m_B.rows(), 128 Bc = m_B.cols(); 129 dst.resize(rows(),cols()); 130 dst.resizeNonZeros(0); 131 dst.reserve(m_A.nonZeros() * m_B.nonZeros()); 132 133 for (Index kA=0; kA < m_A.outerSize(); ++kA) 134 { 135 for (Index kB=0; kB < m_B.outerSize(); ++kB) 136 { 137 for (typename Lhs::InnerIterator itA(m_A,kA); itA; ++itA) 138 { 139 for (typename Rhs::InnerIterator itB(m_B,kB); itB; ++itB) 140 { 141 const Index i = itA.row() * Br + itB.row(), 142 j = itA.col() * Bc + itB.col(); 143 dst.insert(i,j) = itA.value() * itB.value(); 144 } 145 } 146 } 147 } 148} 149 150namespace internal { 151 152template<typename _Lhs, typename _Rhs> 153struct traits<KroneckerProduct<_Lhs,_Rhs> > 154{ 155 typedef typename remove_all<_Lhs>::type Lhs; 156 typedef typename remove_all<_Rhs>::type Rhs; 157 typedef typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar; 158 159 enum { 160 Rows = size_at_compile_time<traits<Lhs>::RowsAtCompileTime, traits<Rhs>::RowsAtCompileTime>::ret, 161 Cols = size_at_compile_time<traits<Lhs>::ColsAtCompileTime, traits<Rhs>::ColsAtCompileTime>::ret, 162 MaxRows = size_at_compile_time<traits<Lhs>::MaxRowsAtCompileTime, traits<Rhs>::MaxRowsAtCompileTime>::ret, 163 MaxCols = size_at_compile_time<traits<Lhs>::MaxColsAtCompileTime, traits<Rhs>::MaxColsAtCompileTime>::ret, 164 CoeffReadCost = Lhs::CoeffReadCost + Rhs::CoeffReadCost + NumTraits<Scalar>::MulCost 165 }; 166 167 typedef Matrix<Scalar,Rows,Cols> ReturnType; 168}; 169 170template<typename _Lhs, typename _Rhs> 171struct traits<KroneckerProductSparse<_Lhs,_Rhs> > 172{ 173 typedef MatrixXpr XprKind; 174 typedef typename remove_all<_Lhs>::type Lhs; 175 typedef typename remove_all<_Rhs>::type Rhs; 176 typedef typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar; 177 typedef typename promote_storage_type<typename traits<Lhs>::StorageKind, typename traits<Rhs>::StorageKind>::ret StorageKind; 178 typedef typename promote_index_type<typename Lhs::Index, typename Rhs::Index>::type Index; 179 180 enum { 181 LhsFlags = Lhs::Flags, 182 RhsFlags = Rhs::Flags, 183 184 RowsAtCompileTime = size_at_compile_time<traits<Lhs>::RowsAtCompileTime, traits<Rhs>::RowsAtCompileTime>::ret, 185 ColsAtCompileTime = size_at_compile_time<traits<Lhs>::ColsAtCompileTime, traits<Rhs>::ColsAtCompileTime>::ret, 186 MaxRowsAtCompileTime = size_at_compile_time<traits<Lhs>::MaxRowsAtCompileTime, traits<Rhs>::MaxRowsAtCompileTime>::ret, 187 MaxColsAtCompileTime = size_at_compile_time<traits<Lhs>::MaxColsAtCompileTime, traits<Rhs>::MaxColsAtCompileTime>::ret, 188 189 EvalToRowMajor = (LhsFlags & RhsFlags & RowMajorBit), 190 RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit), 191 192 Flags = ((LhsFlags | RhsFlags) & HereditaryBits & RemovedBits) 193 | EvalBeforeNestingBit | EvalBeforeAssigningBit, 194 CoeffReadCost = Dynamic 195 }; 196}; 197 198} // end namespace internal 199 200/*! 201 * \ingroup KroneckerProduct_Module 202 * 203 * Computes Kronecker tensor product of two dense matrices 204 * 205 * \warning If you want to replace a matrix by its Kronecker product 206 * with some matrix, do \b NOT do this: 207 * \code 208 * A = kroneckerProduct(A,B); // bug!!! caused by aliasing effect 209 * \endcode 210 * instead, use eval() to work around this: 211 * \code 212 * A = kroneckerProduct(A,B).eval(); 213 * \endcode 214 * 215 * \param a Dense matrix a 216 * \param b Dense matrix b 217 * \return Kronecker tensor product of a and b 218 */ 219template<typename A, typename B> 220KroneckerProduct<A,B> kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b) 221{ 222 return KroneckerProduct<A, B>(a.derived(), b.derived()); 223} 224 225/*! 226 * \ingroup KroneckerProduct_Module 227 * 228 * Computes Kronecker tensor product of two matrices, at least one of 229 * which is sparse 230 * 231 * \param a Dense/sparse matrix a 232 * \param b Dense/sparse matrix b 233 * \return Kronecker tensor product of a and b, stored in a sparse 234 * matrix 235 */ 236template<typename A, typename B> 237KroneckerProductSparse<A,B> kroneckerProduct(const EigenBase<A>& a, const EigenBase<B>& b) 238{ 239 return KroneckerProductSparse<A,B>(a.derived(), b.derived()); 240} 241 242} // end namespace Eigen 243 244#endif // KRONECKER_TENSOR_PRODUCT_H 245