KroneckerTensorProduct.h revision c981c48f5bc9aefeffc0bcb0cc3934c2fae179dd
1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2011 Kolja Brix <brix@igpm.rwth-aachen.de> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2011 Andreas Platen <andiplaten@gmx.de> 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef KRONECKER_TENSOR_PRODUCT_H 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define KRONECKER_TENSOR_PRODUCT_H 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*! 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Kronecker tensor product helper function for dense matrices 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param A Dense matrix A 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param B Dense matrix B 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param AB_ Kronecker tensor product of A and B 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived_A, typename Derived_B, typename Derived_AB> 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid kroneckerProduct_full(const Derived_A& A, const Derived_B& B, Derived_AB & AB) 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const unsigned int Ar = A.rows(), 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Ac = A.cols(), 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Br = B.rows(), 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Bc = B.cols(); 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (unsigned int i=0; i<Ar; ++i) 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (unsigned int j=0; j<Ac; ++j) 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath AB.block(i*Br,j*Bc,Br,Bc) = A(i,j)*B; 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*! 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Kronecker tensor product helper function for matrices, where at least one is sparse 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param A Matrix A 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param B Matrix B 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param AB_ Kronecker tensor product of A and B 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived_A, typename Derived_B, typename Derived_AB> 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid kroneckerProduct_sparse(const Derived_A &A, const Derived_B &B, Derived_AB &AB) 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const unsigned int Ar = A.rows(), 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Ac = A.cols(), 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Br = B.rows(), 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Bc = B.cols(); 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath AB.resize(Ar*Br,Ac*Bc); 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath AB.resizeNonZeros(0); 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath AB.reserve(A.nonZeros()*B.nonZeros()); 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int kA=0; kA<A.outerSize(); ++kA) 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int kB=0; kB<B.outerSize(); ++kB) 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (typename Derived_A::InnerIterator itA(A,kA); itA; ++itA) 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (typename Derived_B::InnerIterator itB(B,kB); itB; ++itB) 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const unsigned int iA = itA.row(), 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jA = itA.col(), 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath iB = itB.row(), 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jB = itB.col(), 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath i = iA*Br + iB, 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath j = jA*Bc + jB; 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath AB.insert(i,j) = itA.value() * itB.value(); 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*! 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Computes Kronecker tensor product of two dense matrices 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param a Dense matrix a 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param b Dense matrix b 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param c Kronecker tensor product of a and b 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename A,typename B,typename CScalar,int CRows,int CCols, int COptions, int CMaxRows, int CMaxCols> 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b, Matrix<CScalar,CRows,CCols,COptions,CMaxRows,CMaxCols>& c) 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath c.resize(a.rows()*b.rows(),a.cols()*b.cols()); 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::kroneckerProduct_full(a.derived(), b.derived(), c); 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*! 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Computes Kronecker tensor product of two dense matrices 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Remark: this function uses the const cast hack and has been 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * implemented to make the function call possible, where the 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * output matrix is a submatrix, e.g. 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * kroneckerProduct(A,B,AB.block(2,5,6,6)); 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param a Dense matrix a 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param b Dense matrix b 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param c Kronecker tensor product of a and b 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename A,typename B,typename C> 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b, MatrixBase<C> const & c_) 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixBase<C>& c = const_cast<MatrixBase<C>& >(c_); 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::kroneckerProduct_full(a.derived(), b.derived(), c.derived()); 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*! 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Computes Kronecker tensor product of a dense and a sparse matrix 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param a Dense matrix a 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param b Sparse matrix b 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param c Kronecker tensor product of a and b 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename A,typename B,typename C> 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid kroneckerProduct(const MatrixBase<A>& a, const SparseMatrixBase<B>& b, SparseMatrixBase<C>& c) 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived()); 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*! 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Computes Kronecker tensor product of a sparse and a dense matrix 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param a Sparse matrix a 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param b Dense matrix b 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param c Kronecker tensor product of a and b 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename A,typename B,typename C> 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid kroneckerProduct(const SparseMatrixBase<A>& a, const MatrixBase<B>& b, SparseMatrixBase<C>& c) 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived()); 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*! 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Computes Kronecker tensor product of two sparse matrices 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param a Sparse matrix a 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param b Sparse matrix b 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param c Kronecker tensor product of a and b 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename A,typename B,typename C> 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid kroneckerProduct(const SparseMatrixBase<A>& a, const SparseMatrixBase<B>& b, SparseMatrixBase<C>& c) 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived()); 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // KRONECKER_TENSOR_PRODUCT_H 158