KroneckerTensorProduct.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2011 Kolja Brix <brix@igpm.rwth-aachen.de>
00005 // Copyright (C) 2011 Andreas Platen <andiplaten@gmx.de>
00006 //
00007 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 
00027 #ifndef KRONECKER_TENSOR_PRODUCT_H
00028 #define KRONECKER_TENSOR_PRODUCT_H
00029 
00030 
00031 namespace Eigen { 
00032 
00033 namespace internal {
00034 
00042 template<typename Derived_A, typename Derived_B, typename Derived_AB>
00043 void kroneckerProduct_full(const Derived_A& A, const Derived_B& B, Derived_AB & AB)
00044 {
00045   const unsigned int Ar = A.rows(),
00046                      Ac = A.cols(),
00047                      Br = B.rows(),
00048                      Bc = B.cols();
00049   for (unsigned int i=0; i<Ar; ++i)
00050     for (unsigned int j=0; j<Ac; ++j)
00051       AB.block(i*Br,j*Bc,Br,Bc) = A(i,j)*B;
00052 }
00053 
00054 
00062 template<typename Derived_A, typename Derived_B, typename Derived_AB>
00063 void kroneckerProduct_sparse(const Derived_A &A, const Derived_B &B, Derived_AB &AB)
00064 {
00065   const unsigned int Ar = A.rows(),
00066                      Ac = A.cols(),
00067                      Br = B.rows(),
00068                      Bc = B.cols();
00069   AB.resize(Ar*Br,Ac*Bc);
00070   AB.resizeNonZeros(0);
00071   AB.reserve(A.nonZeros()*B.nonZeros());
00072 
00073   for (int kA=0; kA<A.outerSize(); ++kA)
00074   {
00075     for (int kB=0; kB<B.outerSize(); ++kB)
00076     {
00077       for (typename Derived_A::InnerIterator itA(A,kA); itA; ++itA)
00078       {
00079         for (typename Derived_B::InnerIterator itB(B,kB); itB; ++itB)
00080         {
00081           const unsigned int iA = itA.row(),
00082                              jA = itA.col(),
00083                              iB = itB.row(),
00084                              jB = itB.col(),
00085                              i  = iA*Br + iB,
00086                              j  = jA*Bc + jB;
00087           AB.insert(i,j) = itA.value() * itB.value();
00088         }
00089       }
00090     }
00091   }
00092 }
00093 
00094 } // end namespace internal
00095 
00096 
00097 
00105 template<typename A,typename B,typename CScalar,int CRows,int CCols, int COptions, int CMaxRows, int CMaxCols>
00106 void kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b, Matrix<CScalar,CRows,CCols,COptions,CMaxRows,CMaxCols>& c)
00107 {
00108   c.resize(a.rows()*b.rows(),a.cols()*b.cols());
00109   internal::kroneckerProduct_full(a.derived(), b.derived(), c);
00110 }
00111 
00124 template<typename A,typename B,typename C>
00125 void kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b, MatrixBase<C> const & c_)
00126 {
00127   MatrixBase<C>& c = const_cast<MatrixBase<C>& >(c_);
00128   internal::kroneckerProduct_full(a.derived(), b.derived(), c.derived());
00129 }
00130 
00138 template<typename A,typename B,typename C>
00139 void kroneckerProduct(const MatrixBase<A>& a, const SparseMatrixBase<B>& b, SparseMatrixBase<C>& c)
00140 {
00141   internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived());
00142 }
00143 
00151 template<typename A,typename B,typename C>
00152 void kroneckerProduct(const SparseMatrixBase<A>& a, const MatrixBase<B>& b, SparseMatrixBase<C>& c)
00153 {
00154   internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived());
00155 }
00156 
00164 template<typename A,typename B,typename C>
00165 void kroneckerProduct(const SparseMatrixBase<A>& a, const SparseMatrixBase<B>& b, SparseMatrixBase<C>& c)
00166 {
00167   internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived());
00168 }
00169 
00170 } // end namespace Eigen
00171 
00172 #endif // KRONECKER_TENSOR_PRODUCT_H