1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2011 Kolja Brix <brix (at) igpm.rwth-aachen.de> 5 // Copyright (C) 2011 Andreas Platen <andiplaten (at) 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> 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> 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> 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> 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> 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> 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> 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
