• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
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. */
KroneckerProductBase(const Lhs & A,const Rhs & B)37     KroneckerProductBase(const Lhs& A, const Rhs& B)
38       : m_A(A), m_B(B)
39     {}
40 
rows()41     inline Index rows() const { return m_A.rows() * m_B.rows(); }
cols()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      */
coeff(Index row,Index col)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      */
coeff(Index i)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. */
KroneckerProduct(const Lhs & A,const Rhs & B)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. */
KroneckerProductSparse(const Lhs & A,const Rhs & B)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>
evalTo(Dest & dst)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>
evalTo(Dest & dst)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