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
20 /*!
21 * Kronecker tensor product helper function for dense matrices
22 *
23 * \param A Dense matrix A
24 * \param B Dense matrix B
25 * \param AB_ Kronecker tensor product of A and B
26 */
27 template<typename Derived_A, typename Derived_B, typename Derived_AB>
kroneckerProduct_full(const Derived_A & A,const Derived_B & B,Derived_AB & 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
40 /*!
41 * Kronecker tensor product helper function for matrices, where at least one is sparse
42 *
43 * \param A Matrix A
44 * \param B Matrix B
45 * \param AB_ Kronecker tensor product of A and B
46 */
47 template<typename Derived_A, typename Derived_B, typename Derived_AB>
kroneckerProduct_sparse(const Derived_A & A,const Derived_B & B,Derived_AB & 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
83 /*!
84 * Computes Kronecker tensor product of two dense matrices
85 *
86 * \param a Dense matrix a
87 * \param b Dense matrix b
88 * \param c Kronecker tensor product of a and b
89 */
90 template<typename A,typename B,typename CScalar,int CRows,int CCols, int COptions, int CMaxRows, int CMaxCols>
kroneckerProduct(const MatrixBase<A> & a,const MatrixBase<B> & b,Matrix<CScalar,CRows,CCols,COptions,CMaxRows,CMaxCols> & c)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
97 /*!
98 * Computes Kronecker tensor product of two dense matrices
99 *
100 * Remark: this function uses the const cast hack and has been
101 * implemented to make the function call possible, where the
102 * output matrix is a submatrix, e.g.
103 * kroneckerProduct(A,B,AB.block(2,5,6,6));
104 *
105 * \param a Dense matrix a
106 * \param b Dense matrix b
107 * \param c Kronecker tensor product of a and b
108 */
109 template<typename A,typename B,typename C>
kroneckerProduct(const MatrixBase<A> & a,const MatrixBase<B> & b,MatrixBase<C> const & 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
116 /*!
117 * Computes Kronecker tensor product of a dense and a sparse matrix
118 *
119 * \param a Dense matrix a
120 * \param b Sparse matrix b
121 * \param c Kronecker tensor product of a and b
122 */
123 template<typename A,typename B,typename C>
kroneckerProduct(const MatrixBase<A> & a,const SparseMatrixBase<B> & b,SparseMatrixBase<C> & 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
129 /*!
130 * Computes Kronecker tensor product of a sparse and a dense matrix
131 *
132 * \param a Sparse matrix a
133 * \param b Dense matrix b
134 * \param c Kronecker tensor product of a and b
135 */
136 template<typename A,typename B,typename C>
kroneckerProduct(const SparseMatrixBase<A> & a,const MatrixBase<B> & b,SparseMatrixBase<C> & 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
142 /*!
143 * Computes Kronecker tensor product of two sparse matrices
144 *
145 * \param a Sparse matrix a
146 * \param b Sparse matrix b
147 * \param c Kronecker tensor product of a and b
148 */
149 template<typename A,typename B,typename C>
kroneckerProduct(const SparseMatrixBase<A> & a,const SparseMatrixBase<B> & b,SparseMatrixBase<C> & 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
158