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 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 
12 #ifndef KRONECKER_TENSOR_PRODUCT_H
13 #define KRONECKER_TENSOR_PRODUCT_H
14 
15 
16 namespace Eigen {
17 
18 namespace internal {
19 
27 template<typename Derived_A, typename Derived_B, typename Derived_AB>
28 void kroneckerProduct_full(const Derived_A& A, const Derived_B& B, Derived_AB & AB)
29 {
30  const unsigned int Ar = A.rows(),
31  Ac = A.cols(),
32  Br = B.rows(),
33  Bc = B.cols();
34  for (unsigned int i=0; i<Ar; ++i)
35  for (unsigned int j=0; j<Ac; ++j)
36  AB.block(i*Br,j*Bc,Br,Bc) = A(i,j)*B;
37 }
38 
39 
47 template<typename Derived_A, typename Derived_B, typename Derived_AB>
48 void kroneckerProduct_sparse(const Derived_A &A, const Derived_B &B, Derived_AB &AB)
49 {
50  const unsigned int Ar = A.rows(),
51  Ac = A.cols(),
52  Br = B.rows(),
53  Bc = B.cols();
54  AB.resize(Ar*Br,Ac*Bc);
55  AB.resizeNonZeros(0);
56  AB.reserve(A.nonZeros()*B.nonZeros());
57 
58  for (int kA=0; kA<A.outerSize(); ++kA)
59  {
60  for (int kB=0; kB<B.outerSize(); ++kB)
61  {
62  for (typename Derived_A::InnerIterator itA(A,kA); itA; ++itA)
63  {
64  for (typename Derived_B::InnerIterator itB(B,kB); itB; ++itB)
65  {
66  const unsigned int iA = itA.row(),
67  jA = itA.col(),
68  iB = itB.row(),
69  jB = itB.col(),
70  i = iA*Br + iB,
71  j = jA*Bc + jB;
72  AB.insert(i,j) = itA.value() * itB.value();
73  }
74  }
75  }
76  }
77 }
78 
79 } // end namespace internal
80 
81 
82 
90 template<typename A,typename B,typename CScalar,int CRows,int CCols, int COptions, int CMaxRows, int CMaxCols>
91 void kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b, Matrix<CScalar,CRows,CCols,COptions,CMaxRows,CMaxCols>& c)
92 {
93  c.resize(a.rows()*b.rows(),a.cols()*b.cols());
94  internal::kroneckerProduct_full(a.derived(), b.derived(), c);
95 }
96 
109 template<typename A,typename B,typename C>
110 void kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b, MatrixBase<C> const & c_)
111 {
112  MatrixBase<C>& c = const_cast<MatrixBase<C>& >(c_);
113  internal::kroneckerProduct_full(a.derived(), b.derived(), c.derived());
114 }
115 
123 template<typename A,typename B,typename C>
124 void kroneckerProduct(const MatrixBase<A>& a, const SparseMatrixBase<B>& b, SparseMatrixBase<C>& c)
125 {
126  internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived());
127 }
128 
136 template<typename A,typename B,typename C>
137 void kroneckerProduct(const SparseMatrixBase<A>& a, const MatrixBase<B>& b, SparseMatrixBase<C>& c)
138 {
139  internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived());
140 }
141 
149 template<typename A,typename B,typename C>
150 void kroneckerProduct(const SparseMatrixBase<A>& a, const SparseMatrixBase<B>& b, SparseMatrixBase<C>& c)
151 {
152  internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived());
153 }
154 
155 } // end namespace Eigen
156 
157 #endif // KRONECKER_TENSOR_PRODUCT_H