KroneckerTensorProduct.h
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 
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:
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 
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 
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 
81 template<typename Lhs, typename Rhs>
82 class KroneckerProduct : public KroneckerProductBase<KroneckerProduct<Lhs,Rhs> >
83 {
84  private:
86  using Base::m_A;
87  using Base::m_B;
88 
89  public:
91  KroneckerProduct(const Lhs& A, const Rhs& B)
92  : Base(A, B)
93  {}
94 
96  template<typename Dest> void evalTo(Dest& dst) const;
97 };
98 
114 template<typename Lhs, typename Rhs>
115 class KroneckerProductSparse : public KroneckerProductBase<KroneckerProductSparse<Lhs,Rhs> >
116 {
117  private:
119  using Base::m_A;
120  using Base::m_B;
121 
122  public:
124  KroneckerProductSparse(const Lhs& A, const Rhs& B)
125  : Base(A, B)
126  {}
127 
129  template<typename Dest> void evalTo(Dest& dst) const;
130 };
131 
132 template<typename Lhs, typename Rhs>
133 template<typename Dest>
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>
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 
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 
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
Scalar coeff(Index row, Index col) const
Definition: KroneckerTensorProduct.h:48
Namespace containing all symbols from the Eigen library.
Definition: AdolcForward:45
KroneckerProduct< A, B > kroneckerProduct(const MatrixBase< A > &a, const MatrixBase< B > &b)
Definition: KroneckerTensorProduct.h:271
void evalTo(Dest &dst) const
Evaluate the Kronecker tensor product.
Definition: KroneckerTensorProduct.h:134
KroneckerProductSparse(const Lhs &A, const Rhs &B)
Constructor.
Definition: KroneckerTensorProduct.h:124
KroneckerProductBase(const Lhs &A, const Rhs &B)
Constructor.
Definition: KroneckerTensorProduct.h:37
KroneckerProduct(const Lhs &A, const Rhs &B)
Constructor.
Definition: KroneckerTensorProduct.h:91
void evalTo(Dest &dst) const
Evaluate the Kronecker tensor product.
Definition: KroneckerTensorProduct.h:147
Kronecker tensor product helper class for dense matrices.
Definition: KroneckerTensorProduct.h:82
The base class of dense and sparse Kronecker product.
Definition: KroneckerTensorProduct.h:25
Kronecker tensor product helper class for sparse matrices.
Definition: KroneckerTensorProduct.h:115
Scalar coeff(Index i) const
Definition: KroneckerTensorProduct.h:58