KroneckerTensorProduct.h revision c981c48f5bc9aefeffc0bcb0cc3934c2fae179dd 
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
16namespace Eigen {
17
18namespace 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 */
27template<typename Derived_A, typename Derived_B, typename Derived_AB>
28void 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 */
47template<typename Derived_A, typename Derived_B, typename Derived_AB>
48void 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 */
90template<typename A,typename B,typename CScalar,int CRows,int CCols, int COptions, int CMaxRows, int CMaxCols>
91void 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 */
109template<typename A,typename B,typename C>
110void 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 */
123template<typename A,typename B,typename C>
124void 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 */
136template<typename A,typename B,typename C>
137void 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 */
149template<typename A,typename B,typename C>
150void 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