00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
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 }
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 }
00171
00172 #endif // KRONECKER_TENSOR_PRODUCT_H