From 08e419dcb151da41f169304f751e5467cf0c7b4a Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Sat, 4 Jul 2009 11:16:27 +0200
Subject: * update sparse module wrt new diagonal matrix impl * fix a bug is
 SparseMatrix

---
 Eigen/src/Core/DiagonalMatrix.h          |  6 +++-
 Eigen/src/Core/util/Constants.h          | 12 +-------
 Eigen/src/Core/util/Meta.h               | 12 +++++++-
 Eigen/src/Sparse/SparseDiagonalProduct.h | 52 +++++++++++++++++++++++---------
 Eigen/src/Sparse/SparseMatrix.h          |  2 +-
 Eigen/src/Sparse/SparseMatrixBase.h      | 17 ++++++++++-
 Eigen/src/Sparse/SparseProduct.h         | 17 ++---------
 test/sparse_product.cpp                  |  1 -
 8 files changed, 75 insertions(+), 44 deletions(-)

diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h
index 413ad5c9f..a092b7794 100644
--- a/Eigen/src/Core/DiagonalMatrix.h
+++ b/Eigen/src/Core/DiagonalMatrix.h
@@ -38,7 +38,8 @@ class DiagonalBase : public MultiplierBase<Derived>
       ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
       MaxRowsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
       MaxColsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
-      IsVectorAtCompileTime = 0
+      IsVectorAtCompileTime = 0,
+      Flags = 0
     };
     
     typedef Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime, 0, MaxRowsAtCompileTime, MaxColsAtCompileTime> DenseMatrixType;
@@ -52,6 +53,9 @@ class DiagonalBase : public MultiplierBase<Derived>
 
     inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); }
     inline DiagonalVectorType& diagonal() { return derived().diagonal(); }
+
+    inline int rows() const { return diagonal().size(); }
+    inline int cols() const { return diagonal().size(); }
     
     template<typename MatrixDerived>
     const DiagonalProduct<MatrixDerived, Derived, DiagonalOnTheLeft>
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h
index 35e4767b8..454b44e52 100644
--- a/Eigen/src/Core/util/Constants.h
+++ b/Eigen/src/Core/util/Constants.h
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
 //
 // Eigen is free software; you can redistribute it and/or
@@ -185,9 +185,6 @@ const unsigned int HereditaryBits = RowMajorBit
                                   | EvalBeforeNestingBit
                                   | EvalBeforeAssigningBit
                                   | SparseBit;
-
-// diagonal means both upper and lower triangular
-const unsigned DiagonalBits = UpperTriangularBit | LowerTriangularBit;
     
 // Possible values for the Mode parameter of part()
 const unsigned int UpperTriangular = UpperTriangularBit;
@@ -198,13 +195,6 @@ const unsigned int SelfAdjoint = SelfAdjointBit;
 const unsigned int UnitUpperTriangular = UpperTriangularBit | UnitDiagBit;
 const unsigned int UnitLowerTriangular = LowerTriangularBit | UnitDiagBit;
 
-template<typename T> struct ei_is_diagonal
-{
-  enum {
-    ret = ( (unsigned int)(T::Flags) & DiagonalBits ) == DiagonalBits
-  };
-};
-
 enum { DiagonalOnTheLeft, DiagonalOnTheRight };
 
 enum { Aligned, Unaligned };
diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h
index 7552dd70c..2f5363275 100644
--- a/Eigen/src/Core/util/Meta.h
+++ b/Eigen/src/Core/util/Meta.h
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
 //
 // Eigen is free software; you can redistribute it and/or
@@ -186,6 +186,16 @@ struct ei_result_of<ei_scalar_product_op<Scalar>(ArgType0,ArgType1)> {
 typedef typename ei_scalar_product_traits<typename ei_cleantype<ArgType0>::type, typename ei_cleantype<ArgType1>::type>::ReturnType type;
 };
 
+template<typename T> struct ei_is_diagonal
+{ enum { ret = false }; };
 
+template<typename T> struct ei_is_diagonal<DiagonalBase<T> >
+{ enum { ret = true }; };
+
+template<typename T> struct ei_is_diagonal<DiagonalWrapper<T> >
+{ enum { ret = true }; };
+
+template<typename T, int S> struct ei_is_diagonal<DiagonalMatrix<T,S> >
+{ enum { ret = true }; };
 
 #endif // EIGEN_META_H
diff --git a/Eigen/src/Sparse/SparseDiagonalProduct.h b/Eigen/src/Sparse/SparseDiagonalProduct.h
index 3a148b09e..5fb149a2c 100644
--- a/Eigen/src/Sparse/SparseDiagonalProduct.h
+++ b/Eigen/src/Sparse/SparseDiagonalProduct.h
@@ -25,8 +25,9 @@
 #ifndef EIGEN_SPARSE_DIAGONAL_PRODUCT_H
 #define EIGEN_SPARSE_DIAGONAL_PRODUCT_H
 
-// the product a diagonal matrix with a sparse matrix can be easily
-// implemented using expression template. We have two very different cases:
+// The product of a diagonal matrix with a sparse matrix can be easily
+// implemented using expression template.
+// We have two consider very different cases:
 // 1 - diag * row-major sparse
 //     => each inner vector <=> scalar * sparse vector product
 //     => so we can reuse CwiseUnaryOp::InnerIterator
@@ -37,13 +38,21 @@
 // The two other cases are symmetric.
 
 template<typename Lhs, typename Rhs>
-struct ei_traits<SparseDiagonalProduct<Lhs, Rhs> > : ei_traits<SparseProduct<Lhs, Rhs, DiagonalProduct> >
+struct ei_traits<SparseDiagonalProduct<Lhs, Rhs> >
 {
   typedef typename ei_cleantype<Lhs>::type _Lhs;
   typedef typename ei_cleantype<Rhs>::type _Rhs;
+  typedef typename _Lhs::Scalar Scalar;
   enum {
+    RowsAtCompileTime = _Lhs::RowsAtCompileTime,
+    ColsAtCompileTime = _Rhs::ColsAtCompileTime,
+
+    MaxRowsAtCompileTime = _Lhs::MaxRowsAtCompileTime,
+    MaxColsAtCompileTime = _Rhs::MaxColsAtCompileTime,
+    
     SparseFlags = ei_is_diagonal<_Lhs>::ret ? int(_Rhs::Flags) : int(_Lhs::Flags),
-    Flags = SparseBit | (SparseFlags&RowMajorBit)
+    Flags = SparseBit | (SparseFlags&RowMajorBit),
+    CoeffReadCost = Dynamic
   };
 };
 
@@ -51,11 +60,16 @@ enum {SDP_IsDiagonal, SDP_IsSparseRowMajor, SDP_IsSparseColMajor};
 template<typename Lhs, typename Rhs, typename SparseDiagonalProductType, int RhsMode, int LhsMode> 
 class ei_sparse_diagonal_product_inner_iterator_selector;
 
-template<typename LhsNested, typename RhsNested>
-class SparseDiagonalProduct : public SparseMatrixBase<SparseDiagonalProduct<LhsNested,RhsNested> >, ei_no_assignment_operator
+template<typename Lhs, typename Rhs>
+class SparseDiagonalProduct
+  : public SparseMatrixBase<SparseDiagonalProduct<Lhs,Rhs> >,
+    ei_no_assignment_operator
 {
-    typedef typename ei_traits<SparseDiagonalProduct>::_LhsNested _LhsNested;
-    typedef typename ei_traits<SparseDiagonalProduct>::_RhsNested _RhsNested;
+    typedef typename Lhs::Nested LhsNested;
+    typedef typename Rhs::Nested RhsNested;
+
+    typedef typename ei_cleantype<LhsNested>::type _LhsNested;
+    typedef typename ei_cleantype<RhsNested>::type _RhsNested;
     
     enum {
       LhsMode = ei_is_diagonal<_LhsNested>::ret ? SDP_IsDiagonal
@@ -71,8 +85,8 @@ class SparseDiagonalProduct : public SparseMatrixBase<SparseDiagonalProduct<LhsN
     typedef ei_sparse_diagonal_product_inner_iterator_selector
                 <_LhsNested,_RhsNested,SparseDiagonalProduct,LhsMode,RhsMode> InnerIterator;
 
-    template<typename Lhs, typename Rhs>
-    EIGEN_STRONG_INLINE SparseDiagonalProduct(const Lhs& lhs, const Rhs& rhs)
+    template<typename _Lhs, typename _Rhs>
+    EIGEN_STRONG_INLINE SparseDiagonalProduct(const _Lhs& lhs, const _Rhs& rhs)
       : m_lhs(lhs), m_rhs(rhs)
     {
       ei_assert(lhs.cols() == rhs.rows() && "invalid sparse matrix * diagonal matrix product");
@@ -109,12 +123,12 @@ class ei_sparse_diagonal_product_inner_iterator_selector
   : public SparseCwiseBinaryOp<
       ei_scalar_product_op<typename Lhs::Scalar>,
       SparseInnerVectorSet<Rhs,1>,
-      typename Lhs::_CoeffsVectorType>::InnerIterator
+      typename Lhs::DiagonalVectorType>::InnerIterator
 {
     typedef typename SparseCwiseBinaryOp<
       ei_scalar_product_op<typename Lhs::Scalar>,
       SparseInnerVectorSet<Rhs,1>,
-      typename Lhs::_CoeffsVectorType>::InnerIterator Base;
+      typename Lhs::DiagonalVectorType>::InnerIterator Base;
   public:
     inline ei_sparse_diagonal_product_inner_iterator_selector(
               const SparseDiagonalProductType& expr, int outer)
@@ -141,12 +155,12 @@ class ei_sparse_diagonal_product_inner_iterator_selector
   : public SparseCwiseBinaryOp<
       ei_scalar_product_op<typename Rhs::Scalar>,
       SparseInnerVectorSet<Lhs,1>,
-      NestByValue<Transpose<typename Rhs::_CoeffsVectorType> > >::InnerIterator
+      NestByValue<Transpose<typename Rhs::DiagonalVectorType> > >::InnerIterator
 {
     typedef typename SparseCwiseBinaryOp<
       ei_scalar_product_op<typename Rhs::Scalar>,
       SparseInnerVectorSet<Lhs,1>,
-      NestByValue<Transpose<typename Rhs::_CoeffsVectorType> > >::InnerIterator Base;
+      NestByValue<Transpose<typename Rhs::DiagonalVectorType> > >::InnerIterator Base;
   public:
     inline ei_sparse_diagonal_product_inner_iterator_selector(
               const SparseDiagonalProductType& expr, int outer)
@@ -154,4 +168,14 @@ class ei_sparse_diagonal_product_inner_iterator_selector
     {}
 };
 
+// SparseMatrixBase functions
+
+template<typename Derived>
+template<typename OtherDerived>
+const SparseDiagonalProduct<Derived,OtherDerived>
+SparseMatrixBase<Derived>::operator*(const DiagonalBase<OtherDerived> &other) const
+{
+  return SparseDiagonalProduct<Derived,OtherDerived>(this->derived(), other.derived());
+}
+
 #endif // EIGEN_SPARSE_DIAGONAL_PRODUCT_H
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index 03bc60bb5..b751b5cb9 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -371,7 +371,7 @@ class SparseMatrix
       const int outerSize = IsRowMajor ? rows : cols;
       m_innerSize = IsRowMajor ? cols : rows;
       m_data.clear();
-      if (m_outerSize != outerSize)
+      if (m_outerSize != outerSize || m_outerSize==0)
       {
         delete[] m_outerIndex;
         m_outerIndex = new int [outerSize+1];
diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h
index 8d9406cfb..8c53757b0 100644
--- a/Eigen/src/Sparse/SparseMatrixBase.h
+++ b/Eigen/src/Sparse/SparseMatrixBase.h
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -37,6 +37,9 @@
   *
   */
 template<typename Derived> class SparseMatrixBase
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+  : public MultiplierBase<Derived>
+#endif // not EIGEN_PARSED_BY_DOXYGEN
 {
   public:
 
@@ -301,10 +304,22 @@ template<typename Derived> class SparseMatrixBase
     { return matrix*scalar; }
 
 
+    // sparse * sparse
     template<typename OtherDerived>
     const typename SparseProductReturnType<Derived,OtherDerived>::Type
     operator*(const SparseMatrixBase<OtherDerived> &other) const;
 
+    // sparse * diagonal
+    template<typename OtherDerived>
+    const SparseDiagonalProduct<Derived,OtherDerived>
+    operator*(const DiagonalBase<OtherDerived> &other) const;
+
+    // diagonal * sparse
+    template<typename OtherDerived> friend 
+    const SparseDiagonalProduct<OtherDerived,Derived>
+    operator*(const DiagonalBase<OtherDerived> &lhs, const SparseMatrixBase& rhs)
+    { return SparseDiagonalProduct<OtherDerived,Derived>(lhs.derived(), rhs.derived()); }
+
     // dense * sparse (return a dense object)
     template<typename OtherDerived> friend
     const typename SparseProductReturnType<OtherDerived,Derived>::Type
diff --git a/Eigen/src/Sparse/SparseProduct.h b/Eigen/src/Sparse/SparseProduct.h
index 5222439d8..a9ad8aad4 100644
--- a/Eigen/src/Sparse/SparseProduct.h
+++ b/Eigen/src/Sparse/SparseProduct.h
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -29,9 +29,7 @@ template<typename Lhs, typename Rhs> struct ei_sparse_product_mode
 {
   enum {
 
-    value = ei_is_diagonal<Lhs>::ret || ei_is_diagonal<Rhs>::ret
-          ? DiagonalProduct
-          :  (Rhs::Flags&Lhs::Flags&SparseBit)==SparseBit
+    value = (Rhs::Flags&Lhs::Flags&SparseBit)==SparseBit
           ? SparseTimeSparseProduct
           : (Lhs::Flags&SparseBit)==SparseBit
           ? SparseTimeDenseProduct
@@ -47,15 +45,6 @@ struct SparseProductReturnType
   typedef SparseProduct<LhsNested, RhsNested, ProductMode> Type;
 };
 
-template<typename Lhs, typename Rhs>
-struct SparseProductReturnType<Lhs,Rhs,DiagonalProduct>
-{
-  typedef const typename ei_nested<Lhs,Rhs::RowsAtCompileTime>::type LhsNested;
-  typedef const typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
-
-  typedef SparseDiagonalProduct<LhsNested, RhsNested> Type;
-};
-
 // sparse product return type specialization
 template<typename Lhs, typename Rhs>
 struct SparseProductReturnType<Lhs,Rhs,SparseTimeSparseProduct>
@@ -106,7 +95,7 @@ struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >
 //     RhsIsRowMajor = (RhsFlags & RowMajorBit)==RowMajorBit,
 
     EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit),
-    ResultIsSparse = ProductMode==SparseTimeSparseProduct || ProductMode==DiagonalProduct,
+    ResultIsSparse = ProductMode==SparseTimeSparseProduct,
 
     RemovedBits = ~( (EvalToRowMajor ? 0 : RowMajorBit) | (ResultIsSparse ? 0 : SparseBit) ),
 
diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp
index 122fd4325..743273f65 100644
--- a/test/sparse_product.cpp
+++ b/test/sparse_product.cpp
@@ -118,7 +118,6 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType&
     VERIFY_IS_APPROX(x=mLo.template marked<LowerTriangular|SelfAdjoint>()*b, refX=refS*b);
     VERIFY_IS_APPROX(x=mS.template marked<SelfAdjoint>()*b, refX=refS*b);
   }
-
 }
 
 void test_sparse_product()
-- 
cgit v1.2.3