Sparse module: add support for sparse selfadjoint * dense

2 files changed, 48 insertions, 5 deletions

diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h
index 953dd30a7..dd4eeff16 100644
--- a/Eigen/src/Sparse/SparseMatrixBase.h
+++ b/Eigen/src/Sparse/SparseMatrixBase.h
@@ -88,7 +88,7 @@ template<typename Derived> class SparseMatrixBase
     /** \internal the return type of MatrixBase::imag() */
     typedef CwiseUnaryOp<ei_scalar_imag_op<Scalar>, Derived> ImagReturnType;
     /** \internal the return type of MatrixBase::adjoint() */
-    typedef Eigen::Transpose<NestByValue<typename ei_cleantype<ConjugateReturnType>::type> >
+    typedef SparseTranspose</*NestByValue<*/typename ei_cleantype<ConjugateReturnType>::type> /*>*/
             AdjointReturnType;
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
@@ -322,7 +322,7 @@ template<typename Derived> class SparseMatrixBase
     SparseTranspose<Derived> transpose() { return derived(); }
     const SparseTranspose<Derived> transpose() const { return derived(); }
     // void transposeInPlace();
-    // const AdjointReturnType adjoint() const;
+    const AdjointReturnType adjoint() const { return conjugate()/*.nestByValue()*/; }
 
     SparseInnerVector<Derived> innerVector(int outer);
     const SparseInnerVector<Derived> innerVector(int outer) const;
diff --git a/Eigen/src/Sparse/SparseProduct.h b/Eigen/src/Sparse/SparseProduct.h
index 29f5208fa..5a2c294a2 100644
--- a/Eigen/src/Sparse/SparseProduct.h
+++ b/Eigen/src/Sparse/SparseProduct.h
@@ -294,17 +294,60 @@ inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs
 }
 
 // dense = sparse * dense
+// template<typename Derived>
+// template<typename Lhs, typename Rhs>
+// Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeDenseProduct>& product)
+// {
+//   typedef typename ei_cleantype<Lhs>::type _Lhs;
+//   typedef typename _Lhs::InnerIterator LhsInnerIterator;
+//   enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit };
+//   derived().setZero();
+//   for (int j=0; j<product.lhs().outerSize(); ++j)
+//     for (LhsInnerIterator i(product.lhs(),j); i; ++i)
+//       derived().row(LhsIsRowMajor ? j : i.index()) += i.value() * product.rhs().row(LhsIsRowMajor ? i.index() : j);
+//   return derived();
+// }
+
 template<typename Derived>
 template<typename Lhs, typename Rhs>
 Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeDenseProduct>& product)
 {
   typedef typename ei_cleantype<Lhs>::type _Lhs;
   typedef typename _Lhs::InnerIterator LhsInnerIterator;
-  enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit };
+  enum {
+    LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit,
+    LhsIsSelfAdjoint = (_Lhs::Flags&SelfAdjointBit)==SelfAdjointBit,
+    ProcessFirstHalf = LhsIsSelfAdjoint
+      && (   ((_Lhs::Flags&(UpperTriangularBit|LowerTriangularBit))==0)
+          || ( (_Lhs::Flags&UpperTriangularBit) && !LhsIsRowMajor)
+          || ( (_Lhs::Flags&LowerTriangularBit) && LhsIsRowMajor) ),
+    ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf)
+  };
   derived().setZero();
   for (int j=0; j<product.lhs().outerSize(); ++j)
-    for (LhsInnerIterator i(product.lhs(),j); i; ++i)
-      derived().row(LhsIsRowMajor ? j : i.index()) += i.value() * product.rhs().row(LhsIsRowMajor ? i.index() : j);
+  {
+    LhsInnerIterator i(product.lhs(),j);
+    if (ProcessSecondHalf && i && (i.index()==j))
+    {
+      derived().row(j) += i.value() * product.rhs().row(j);
+      ++i;
+    }
+    for (; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
+    {
+      if (LhsIsSelfAdjoint)
+      {
+        int a = LhsIsRowMajor ? j : i.index();
+        int b = LhsIsRowMajor ? i.index() : j;
+        Scalar v = i.value();
+        derived().row(a) += (v) * product.rhs().row(b);
+        derived().row(b) += ei_conj(v) * product.rhs().row(a);
+      }
+      else
+        derived().row(LhsIsRowMajor ? j : i.index()) += i.value() * product.rhs().row(LhsIsRowMajor ? i.index() : j);
+    }
+    if (ProcessFirstHalf && i && (i.index()==j))
+      derived().row(j) += i.value() * product.rhs().row(j);
+  }
   return derived();
 }