Make kroneckerProduct take two arguments and return an expression, which is more straight-forward.

Eigen/src/Core/MatrixBase.h

@@ -165,6 +165,9 @@ template<typename Derived> class MatrixBase
 
     template<typename ProductDerived, typename Lhs, typename Rhs>
     Derived& lazyAssign(const MatrixPowerProductBase<ProductDerived, Lhs,Rhs>& other);
+
+    template<typename Lhs, typename Rhs>
+    Derived& lazyAssign(const KroneckerProduct<Lhs,Rhs>& other);
 #endif // not EIGEN_PARSED_BY_DOXYGEN
 
     template<typename OtherDerived>

Eigen/src/Core/util/ForwardDeclarations.h

@@ -282,6 +282,9 @@ struct stem_function
 };
 }
 
+// KroneckerProduct module
+template<typename Lhs, typename Rhs> class KroneckerProduct;
+template<typename Lhs, typename Rhs> class KroneckerProductSparse;
 
 #ifdef EIGEN2_SUPPORT
 template<typename ExpressionType> class Cwise;

Eigen/src/SparseCore/SparseMatrixBase.h

@@ -259,6 +259,9 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived>
     template<typename Lhs, typename Rhs>
     inline Derived& operator=(const SparseSparseProduct<Lhs,Rhs>& product);
 
+    template<typename Lhs, typename Rhs>
+    inline Derived& operator=(const KroneckerProductSparse<Lhs,Rhs>& product);
+
     friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m)
     {
       typedef typename Derived::Nested Nested;

unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h

@@ -3,138 +3,230 @@
 //
 // Copyright (C) 2011 Kolja Brix <brix@igpm.rwth-aachen.de>
 // Copyright (C) 2011 Andreas Platen <andiplaten@gmx.de>
+// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-
 #ifndef KRONECKER_TENSOR_PRODUCT_H
 #define KRONECKER_TENSOR_PRODUCT_H
 
+#define EIGEN_SIZE_PRODUCT(a,b) (!((int)a && (int)b) ? 0 \
+                : ((int)a == Dynamic || (int)b == Dynamic) ? Dynamic \
+                : (int)a * (int)b)
 
 namespace Eigen { 
 
 namespace internal {
 
-/*!
+template<typename _Lhs, typename _Rhs>
- * Kronecker tensor product helper function for dense matrices
+struct traits<KroneckerProduct<_Lhs,_Rhs> >
- *
- * \param A   Dense matrix A
- * \param B   Dense matrix B
- * \param AB_ Kronecker tensor product of A and B
- */
-template<typename Derived_A, typename Derived_B, typename Derived_AB>
-void kroneckerProduct_full(const Derived_A& A, const Derived_B& B, Derived_AB & AB)
 {
-  const unsigned int Ar = A.rows(),
+  typedef MatrixXpr XprKind;
-                     Ac = A.cols(),
+  typedef typename remove_all<_Lhs>::type Lhs;
-                     Br = B.rows(),
+  typedef typename remove_all<_Rhs>::type Rhs;
-                     Bc = B.cols();
+  typedef typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar;
-  AB.resize(Ar*Br,Ac*Bc);
+  typedef Dense StorageKind;
+  typedef typename promote_index_type<typename Lhs::Index, typename Rhs::Index>::type Index;
 
-  for (unsigned int i=0; i<Ar; ++i)
+  enum {
-    for (unsigned int j=0; j<Ac; ++j)
+    RowsAtCompileTime = EIGEN_SIZE_PRODUCT(traits<Lhs>::RowsAtCompileTime, traits<Rhs>::RowsAtCompileTime),
-      AB.block(i*Br,j*Bc,Br,Bc) = A(i,j)*B;
+    ColsAtCompileTime = EIGEN_SIZE_PRODUCT(traits<Lhs>::ColsAtCompileTime, traits<Rhs>::ColsAtCompileTime),
-}
+    MaxRowsAtCompileTime = EIGEN_SIZE_PRODUCT(traits<Lhs>::MaxRowsAtCompileTime, traits<Rhs>::MaxRowsAtCompileTime),
+    MaxColsAtCompileTime = EIGEN_SIZE_PRODUCT(traits<Lhs>::MaxColsAtCompileTime, traits<Rhs>::MaxColsAtCompileTime),
+    Flags = (MaxRowsAtCompileTime==1 ? RowMajorBit : 0)
+          | EvalBeforeNestingBit | EvalBeforeAssigningBit | NestByRefBit,
+    CoeffReadCost = Lhs::CoeffReadCost + Rhs::CoeffReadCost + NumTraits<Scalar>::MulCost
+  };
+};
 
-
+template<typename _Lhs, typename _Rhs>
-/*!
+struct traits<KroneckerProductSparse<_Lhs,_Rhs> >
- * Kronecker tensor product helper function for matrices, where at least one is sparse
- *
- * \param A   Matrix A
- * \param B   Matrix B
- * \param AB_ Kronecker tensor product of A and B
- */
-template<typename Derived_A, typename Derived_B, typename Derived_AB>
-void kroneckerProduct_sparse(const Derived_A &A, const Derived_B &B, Derived_AB &AB)
 {
-  const unsigned int Ar = A.rows(),
+  typedef MatrixXpr XprKind;
-                     Ac = A.cols(),
+  typedef typename remove_all<_Lhs>::type Lhs;
-                     Br = B.rows(),
+  typedef typename remove_all<_Rhs>::type Rhs;
-                     Bc = B.cols();
+  typedef typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar;
-  AB.resize(Ar*Br,Ac*Bc);
+  typedef Sparse StorageKind;
-  AB.resizeNonZeros(0);
+  typedef typename promote_index_type<typename Lhs::Index, typename Rhs::Index>::type Index;
-  AB.reserve(A.nonZeros()*B.nonZeros());
 
-  for (int kA=0; kA<A.outerSize(); ++kA)
+  enum {
-  {
+    LhsFlags = Lhs::Flags,
-    for (int kB=0; kB<B.outerSize(); ++kB)
+    RhsFlags = Rhs::Flags,
-    {
+
-      for (typename Derived_A::InnerIterator itA(A,kA); itA; ++itA)
+    RowsAtCompileTime = EIGEN_SIZE_PRODUCT(traits<Lhs>::RowsAtCompileTime, traits<Rhs>::RowsAtCompileTime),
-      {
+    ColsAtCompileTime = EIGEN_SIZE_PRODUCT(traits<Lhs>::ColsAtCompileTime, traits<Rhs>::ColsAtCompileTime),
-        for (typename Derived_B::InnerIterator itB(B,kB); itB; ++itB)
+    MaxRowsAtCompileTime = EIGEN_SIZE_PRODUCT(traits<Lhs>::MaxRowsAtCompileTime, traits<Rhs>::MaxRowsAtCompileTime),
-        {
+    MaxColsAtCompileTime = EIGEN_SIZE_PRODUCT(traits<Lhs>::MaxColsAtCompileTime, traits<Rhs>::MaxColsAtCompileTime),
-          const unsigned int iA = itA.row(),
+
-                             jA = itA.col(),
+    EvalToRowMajor = (LhsFlags & RhsFlags & RowMajorBit),
-                             iB = itB.row(),
+    RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit),
-                             jB = itB.col(),
+
-                             i  = iA*Br + iB,
+    Flags = ((LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
-                             j  = jA*Bc + jB;
+          | EvalBeforeNestingBit | EvalBeforeAssigningBit,
-          AB.insert(i,j) = itA.value() * itB.value();
+    CoeffReadCost = Dynamic
-        }
+  };
-      }
+};
-    }
-  }
-}
 
 } // end namespace internal
 
+/*!
+ * \brief Kronecker tensor product helper class for dense matrices
+ *
+ * This class is the return value of kroneckerProduct(MatrixBase,
+ * MatrixBase). Use the function rather than construct this class
+ * directly to avoid specifying template prarameters.
+ *
+ * \tparam Lhs  Type of the left-hand side, a matrix expression.
+ * \tparam Rhs  Type of the rignt-hand side, a matrix expression.
+ */
+template<typename Lhs, typename Rhs>
+class KroneckerProduct : public MatrixBase<KroneckerProduct<Lhs,Rhs> >
+{
+  public:
+    typedef MatrixBase<KroneckerProduct> Base;
+    EIGEN_DENSE_PUBLIC_INTERFACE(KroneckerProduct)
+
+    /*! \brief Constructor. */
+    KroneckerProduct(const Lhs& A, const Rhs& B)
+      : m_A(A), m_B(B)
+    {}
+
+    /*! \brief Evaluate the Kronecker tensor product. */
+    template<typename Dest> void evalTo(Dest& dst) const;
+    
+    inline Index rows() const { return m_A.rows() * m_B.rows(); }
+    inline Index cols() const { return m_A.cols() * m_B.cols(); }
+
+    typename Base::CoeffReturnType coeff(Index row, Index col) const
+    {
+      return m_A.coeff(row / m_A.cols(), col / m_A.rows()) *
+             m_B.coeff(row % m_A.cols(), col % m_A.rows());
+    }
+
+    typename Base::CoeffReturnType coeff(Index i) const
+    {
+      EIGEN_STATIC_ASSERT_VECTOR_ONLY(KroneckerProduct);
+      return m_A.coeff(i / m_A.size()) * m_B.coeff(i % m_A.size());
+    }
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+    struct Unusable {};
+    Unusable& coeffRef(Index) { return *reinterpret_cast<Unusable*>(this); }
+    Unusable& coeffRef(Index,Index) { return *reinterpret_cast<Unusable*>(this); }
+#endif
+
+  private:
+    typename Lhs::Nested m_A;
+    typename Rhs::Nested m_B;
+};
+
+template<typename Lhs, typename Rhs>
+class KroneckerProductSparse : public SparseMatrixBase<KroneckerProductSparse<Lhs,Rhs> >
+{
+  public:
+    typedef SparseMatrixBase<KroneckerProductSparse> Base;
+    EIGEN_DENSE_PUBLIC_INTERFACE(KroneckerProductSparse)
+
+    /*! \brief Constructor. */
+    KroneckerProductSparse(const Lhs& A, const Rhs& B)
+      : m_A(A), m_B(B)
+    {}
+
+    /*! \brief Evaluate the Kronecker tensor product. */
+    template<typename Dest> void evalTo(Dest& dst) const;
+    
+    inline Index rows() const { return m_A.rows() * m_B.rows(); }
+    inline Index cols() const { return m_A.cols() * m_B.cols(); }
+
+  private:
+    typename Lhs::Nested m_A;
+    typename Rhs::Nested m_B;
+};
+
+template<typename Lhs, typename Rhs>
+template<typename Dest>
+void KroneckerProduct<Lhs,Rhs>::evalTo(Dest& dst) const
+{
+  const int BlockRows = Rhs::RowsAtCompileTime,
+            BlockCols = Rhs::ColsAtCompileTime;
+  const Index Br = m_B.rows(),
+              Bc = m_B.cols();
+  for (Index i=0; i < m_A.rows(); ++i)
+    for (Index j=0; j < m_A.cols(); ++j)
+      Block<Dest,BlockRows,BlockCols>(dst,i*Br,j*Bc,Br,Bc) = m_A.coeff(i,j) * m_B;
+}
+
+template<typename Lhs, typename Rhs>
+template<typename Dest>
+void KroneckerProductSparse<Lhs,Rhs>::evalTo(Dest& dst) const
+{
+  const Index Br = m_B.rows(),
+              Bc = m_B.cols();
+  dst.resize(rows(),cols());
+  dst.resizeNonZeros(0);
+  dst.reserve(m_A.nonZeros() * m_B.nonZeros());
+
+  for (Index kA=0; kA < m_A.outerSize(); ++kA)
+  {
+    for (Index kB=0; kB < m_B.outerSize(); ++kB)
+    {
+      for (typename Lhs::InnerIterator itA(m_A,kA); itA; ++itA)
+      {
+        for (typename Rhs::InnerIterator itB(m_B,kB); itB; ++itB)
+        {
+          const Index i = itA.row() * Br + itB.row(),
+                      j = itA.col() * Bc + itB.col();
+          dst.insert(i,j) = itA.value() * itB.value();
+        }
+      }
+    }
+  }
+}
+
 /*!
  * Computes Kronecker tensor product of two dense matrices
  *
- * Remark: this function uses the const cast hack and has been
- *         implemented to make the function call possible, where the
- *         output matrix is a submatrix, e.g.
- *           kroneckerProduct(A,B,AB.block(2,5,6,6));
- *
  * \param a  Dense matrix a
  * \param b  Dense matrix b
- * \param c  Kronecker tensor product of a and b
+ * \return   Kronecker tensor product of a and b
  */
-template<typename A,typename B,typename C>
+template<typename A, typename B>
-void kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b, const MatrixBase<C>& c)
+KroneckerProduct<A,B> kroneckerProduct(const MatrixBase<A>& a, const MatrixBase<B>& b)
 {
-  internal::kroneckerProduct_full(a.derived(), b.derived(), c.const_cast_derived());
+  return KroneckerProduct<A, B>(a.derived(), b.derived());
 }
 
 /*!
- * Computes Kronecker tensor product of a dense and a sparse matrix
+ * Computes Kronecker tensor product of two matrices, at least one of
+ * which is sparse.
  *
- * \param a  Dense matrix a
+ * \param a  Dense/sparse matrix a
- * \param b  Sparse matrix b
+ * \param b  Dense/sparse matrix b
- * \param c  Kronecker tensor product of a and b
+ * \return   Kronecker tensor product of a and b, stored in a sparse
+ *           matrix
  */
-template<typename A,typename B,typename C>
+template<typename A, typename B>
-void kroneckerProduct(const MatrixBase<A>& a, const SparseMatrixBase<B>& b, SparseMatrixBase<C>& c)
+KroneckerProductSparse<A,B> kroneckerProduct(const EigenBase<A>& a, const EigenBase<B>& b)
 {
-  internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived());
+  return KroneckerProductSparse<A,B>(a.derived(), b.derived());
 }
 
-/*!
+template<typename Derived>
- * Computes Kronecker tensor product of a sparse and a dense matrix
+template<typename Lhs, typename Rhs>
- *
+Derived& MatrixBase<Derived>::lazyAssign(const KroneckerProduct<Lhs,Rhs>& other)
- * \param a  Sparse matrix a
- * \param b  Dense matrix b
- * \param c  Kronecker tensor product of a and b
- */
-template<typename A,typename B,typename C>
-void kroneckerProduct(const SparseMatrixBase<A>& a, const MatrixBase<B>& b, SparseMatrixBase<C>& c)
 {
-  internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived());
+  other.evalTo(derived());
+  return derived();
 }
 
-/*!
+template<typename Derived>
- * Computes Kronecker tensor product of two sparse matrices
+template<typename Lhs, typename Rhs>
- *
+Derived& SparseMatrixBase<Derived>::operator=(const KroneckerProductSparse<Lhs,Rhs>& product)
- * \param a  Sparse matrix a
- * \param b  Sparse matrix b
- * \param c  Kronecker tensor product of a and b
- */
-template<typename A,typename B,typename C>
-void kroneckerProduct(const SparseMatrixBase<A>& a, const SparseMatrixBase<B>& b, SparseMatrixBase<C>& c)
 {
-  internal::kroneckerProduct_sparse(a.derived(), b.derived(), c.derived());
+  product.evalTo(derived());
+  return derived();
 }
 
 } // end namespace Eigen

unsupported/test/kronecker_product.cpp

@@ -3,6 +3,7 @@
 //
 // Copyright (C) 2011 Kolja Brix <brix@igpm.rwth-aachen.de>
 // Copyright (C) 2011 Andreas Platen <andiplaten@gmx.de>
+// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -89,64 +90,55 @@ void test_kronecker_product()
   MatrixXd             DM_b(3,2);
   SparseMatrix<double> SM_a(2,3);
   SparseMatrix<double> SM_b(3,2);
-  SM_a.insert(0,0) = DM_a(0,0) = -0.4461540300782201;
+  SM_a.insert(0,0) = DM_a.coeffRef(0,0) = -0.4461540300782201;
-  SM_a.insert(0,1) = DM_a(0,1) = -0.8057364375283049;
+  SM_a.insert(0,1) = DM_a.coeffRef(0,1) = -0.8057364375283049;
-  SM_a.insert(0,2) = DM_a(0,2) =  0.3896572459516341;
+  SM_a.insert(0,2) = DM_a.coeffRef(0,2) =  0.3896572459516341;
-  SM_a.insert(1,0) = DM_a(1,0) = -0.9076572187376921;
+  SM_a.insert(1,0) = DM_a.coeffRef(1,0) = -0.9076572187376921;
-  SM_a.insert(1,1) = DM_a(1,1) =  0.6469156566545853;
+  SM_a.insert(1,1) = DM_a.coeffRef(1,1) =  0.6469156566545853;
-  SM_a.insert(1,2) = DM_a(1,2) = -0.3658010398782789;
+  SM_a.insert(1,2) = DM_a.coeffRef(1,2) = -0.3658010398782789;
-  SM_b.insert(0,0) = DM_b(0,0) =  0.9004440976767099;
+  SM_b.insert(0,0) = DM_b.coeffRef(0,0) =  0.9004440976767099;
-  SM_b.insert(0,1) = DM_b(0,1) = -0.2368830858139832;
+  SM_b.insert(0,1) = DM_b.coeffRef(0,1) = -0.2368830858139832;
-  SM_b.insert(1,0) = DM_b(1,0) = -0.9311078389941825;
+  SM_b.insert(1,0) = DM_b.coeffRef(1,0) = -0.9311078389941825;
-  SM_b.insert(1,1) = DM_b(1,1) =  0.5310335762980047;
+  SM_b.insert(1,1) = DM_b.coeffRef(1,1) =  0.5310335762980047;
-  SM_b.insert(2,0) = DM_b(2,0) = -0.1225112806872035;
+  SM_b.insert(2,0) = DM_b.coeffRef(2,0) = -0.1225112806872035;
-  SM_b.insert(2,1) = DM_b(2,1) =  0.5903998022741264;
+  SM_b.insert(2,1) = DM_b.coeffRef(2,1) =  0.5903998022741264;
   SparseMatrix<double,RowMajor> SM_row_a(SM_a), SM_row_b(SM_b);
 
   // test kroneckerProduct(DM_block,DM,DM_fixedSize)
-  Matrix<double, 6, 6> DM_fix_ab;
+  Matrix<double, 6, 6> DM_fix_ab = kroneckerProduct(DM_a.topLeftCorner<2,3>(),DM_b);
-  DM_fix_ab(0,0)=37.0;
-  kroneckerProduct(DM_a.block(0,0,2,3),DM_b,DM_fix_ab);
   CALL_SUBTEST(check_kronecker_product(DM_fix_ab));
 
   // test kroneckerProduct(DM,DM,DM_block)
   MatrixXd DM_block_ab(10,15);
-  DM_block_ab(0,0)=37.0;
+  DM_block_ab.block<6,6>(2,5) = kroneckerProduct(DM_a,DM_b);
-  kroneckerProduct(DM_a,DM_b,DM_block_ab.block(2,5,6,6));
+  CALL_SUBTEST(check_kronecker_product(DM_block_ab.block<6,6>(2,5)));
-  CALL_SUBTEST(check_kronecker_product(DM_block_ab.block(2,5,6,6)));
 
   // test kroneckerProduct(DM,DM,DM)
-  MatrixXd DM_ab(1,5);
+  MatrixXd DM_ab = kroneckerProduct(DM_a,DM_b);
-  DM_ab(0,0)=37.0;
-  kroneckerProduct(DM_a,DM_b,DM_ab);
   CALL_SUBTEST(check_kronecker_product(DM_ab));
 
   // test kroneckerProduct(SM,DM,SM)
-  SparseMatrix<double> SM_ab(1,20);
+  SparseMatrix<double> SM_ab = kroneckerProduct(SM_a,DM_b);
-  SM_ab.insert(0,0)=37.0;
-  kroneckerProduct(SM_a,DM_b,SM_ab);
   CALL_SUBTEST(check_kronecker_product(SM_ab));
-  SparseMatrix<double,RowMajor> SM_ab2(10,3);
+  SparseMatrix<double,RowMajor> SM_ab2 = kroneckerProduct(SM_a,DM_b);
-  SM_ab2.insert(0,0)=37.0;
-  kroneckerProduct(SM_a,DM_b,SM_ab2);
   CALL_SUBTEST(check_kronecker_product(SM_ab2));
 
   // test kroneckerProduct(DM,SM,SM)
   SM_ab.insert(0,0)=37.0;
-  kroneckerProduct(DM_a,SM_b,SM_ab);
+  SM_ab = kroneckerProduct(DM_a,SM_b);
   CALL_SUBTEST(check_kronecker_product(SM_ab));
   SM_ab2.insert(0,0)=37.0;
-  kroneckerProduct(DM_a,SM_b,SM_ab2);
+  SM_ab2 = kroneckerProduct(DM_a,SM_b);
   CALL_SUBTEST(check_kronecker_product(SM_ab2));
 
   // test kroneckerProduct(SM,SM,SM)
   SM_ab.resize(2,33);
   SM_ab.insert(0,0)=37.0;
-  kroneckerProduct(SM_a,SM_b,SM_ab);
+  SM_ab = kroneckerProduct(SM_a,SM_b);
   CALL_SUBTEST(check_kronecker_product(SM_ab));
   SM_ab2.resize(5,11);
   SM_ab2.insert(0,0)=37.0;
-  kroneckerProduct(SM_a,SM_b,SM_ab2);
+  SM_ab2 = kroneckerProduct(SM_a,SM_b);
   CALL_SUBTEST(check_kronecker_product(SM_ab2));
 
   // test kroneckerProduct(SM,SM,SM) with sparse pattern
@@ -163,17 +155,16 @@ void test_kronecker_product()
   SM_b.finalize();
   SM_ab.resize(1,1);
   SM_ab.insert(0,0)=37.0;
-  kroneckerProduct(SM_a,SM_b,SM_ab);
+  SM_ab = kroneckerProduct(SM_a,SM_b);
   CALL_SUBTEST(check_sparse_kronecker_product(SM_ab));
 
   // test dimension of result of kroneckerProduct(DM,DM,DM)
   MatrixXd DM_a2(2,1);
   MatrixXd DM_b2(5,4);
-  MatrixXd DM_ab2;
+  MatrixXd DM_ab2 = kroneckerProduct(DM_a2,DM_b2);
-  kroneckerProduct(DM_a2,DM_b2,DM_ab2);
   CALL_SUBTEST(check_dimension(DM_ab2,2*5,1*4));
   DM_a2.resize(10,9);
   DM_b2.resize(4,8);
-  kroneckerProduct(DM_a2,DM_b2,DM_ab2);
+  DM_ab2 = kroneckerProduct(DM_a2,DM_b2);
   CALL_SUBTEST(check_dimension(DM_ab2,10*4,9*8));
 }