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 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25 
26 
27 #ifndef KRONECKER_TENSOR_PRODUCT_H
28 #define KRONECKER_TENSOR_PRODUCT_H
29 
30 
31 namespace Eigen {
32 
33 namespace internal {
34 
42 template<typename Derived_A, typename Derived_B, typename Derived_AB>
43 void kroneckerProduct_full(const Derived_A& A, const Derived_B& B, Derived_AB & AB)
44 {
45  const unsigned int Ar = A.rows(),
46  Ac = A.cols(),
47  Br = B.rows(),
48  Bc = B.cols();
49  for (unsigned int i=0; i<Ar; ++i)
50  for (unsigned int j=0; j<Ac; ++j)
51  AB.block(i*Br,j*Bc,Br,Bc) = A(i,j)*B;
52 }
53 
54 
62 template<typename Derived_A, typename Derived_B, typename Derived_AB>
63 void kroneckerProduct_sparse(const Derived_A &A, const Derived_B &B, Derived_AB &AB)
64 {
65  const unsigned int Ar = A.rows(),
66  Ac = A.cols(),
67  Br = B.rows(),
68  Bc = B.cols();
69  AB.resize(Ar*Br,Ac*Bc);
70  AB.resizeNonZeros(0);
71  AB.reserve(A.nonZeros()*B.nonZeros());
72 
73  for (int kA=0; kA<A.outerSize(); ++kA)
74  {
75  for (int kB=0; kB<B.outerSize(); ++kB)
76  {
77  for (typename Derived_A::InnerIterator itA(A,kA); itA; ++itA)
78  {
79  for (typename Derived_B::InnerIterator itB(B,kB); itB; ++itB)
80  {
81  const unsigned int iA = itA.row(),
82  jA = itA.col(),
83  iB = itB.row(),
84  jB = itB.col(),
85  i = iA*Br + iB,
86  j = jA*Bc + jB;
87  AB.insert(i,j) = itA.value() * itB.value();
88  }
89  }
90  }
91  }
92 }
93 
94 } // end namespace internal
95 
96 
97 
105 template<typename A,typename B,typename CScalar,int CRows,int CCols, int COptions, int CMaxRows, int CMaxCols>
106 void kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b, Matrix<CScalar,CRows,CCols,COptions,CMaxRows,CMaxCols>& c)
107 {
108  c.resize(a.rows()*b.rows(),a.cols()*b.cols());
109  internal::kroneckerProduct_full(a.derived(), b.derived(), c);
110 }
111 
124 template<typename A,typename B,typename C>
125 void kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b, MatrixBase<C> const & c_)
126 {
127  MatrixBase<C>& c = const_cast<MatrixBase<C>& >(c_);
128  internal::kroneckerProduct_full(a.derived(), b.derived(), c.derived());
129 }
130 
138 template<typename A,typename B,typename C>
139 void kroneckerProduct(const MatrixBase<A>& a, const SparseMatrixBase<B>& b, SparseMatrixBase<C>& c)
140 {
141  internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived());
142 }
143 
151 template<typename A,typename B,typename C>
152 void kroneckerProduct(const SparseMatrixBase<A>& a, const MatrixBase<B>& b, SparseMatrixBase<C>& c)
153 {
154  internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived());
155 }
156 
164 template<typename A,typename B,typename C>
165 void kroneckerProduct(const SparseMatrixBase<A>& a, const SparseMatrixBase<B>& b, SparseMatrixBase<C>& c)
166 {
167  internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived());
168 }
169 
170 } // end namespace Eigen
171 
172 #endif // KRONECKER_TENSOR_PRODUCT_H