From ccdcebcf03a529b429feda831ff4b44f8433e045 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Thu, 15 Jan 2009 18:52:14 +0000
Subject: [PATCH] Sparse module: add support for sparse selfadjoint * dense

---
 Eigen/src/Sparse/SparseMatrixBase.h |  4 +--
 Eigen/src/Sparse/SparseProduct.h    | 49 +++++++++++++++++++++++++++--
 test/sparse.h                       |  7 ++++-
 test/sparse_basic.cpp               | 33 +++++++++++++++++++
 4 files changed, 87 insertions(+), 6 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();
 }
 
diff --git a/test/sparse.h b/test/sparse.h
index feb8bc82d..919afb760 100644
--- a/test/sparse.h
+++ b/test/sparse.h
@@ -40,7 +40,8 @@
 enum {
   ForceNonZeroDiag = 1,
   MakeLowerTriangular = 2,
-  MakeUpperTriangular = 4
+  MakeUpperTriangular = 4,
+  ForceRealDiag = 8
 };
 
 /* Initializes both a sparse and dense matrix with same random values,
@@ -73,6 +74,10 @@ initSparse(double density,
         v = Scalar(0);
       else if ((flags & MakeUpperTriangular) && j<i)
         v = Scalar(0);
+      
+      if ((flags&ForceRealDiag) && (i==j))
+        v = ei_real(v);
+        
       if (v!=Scalar(0))
       {
         sparseMat.fill(i,j) = v;
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp
index 40d452398..31244000d 100644
--- a/test/sparse_basic.cpp
+++ b/test/sparse_basic.cpp
@@ -269,6 +269,39 @@ template<typename Scalar> void sparse_basic(int rows, int cols)
     VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
     VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
   }
+  
+  // test self adjoint products
+  {
+    DenseMatrix b = DenseMatrix::Random(rows, rows);
+    DenseMatrix x = DenseMatrix::Random(rows, rows);
+    DenseMatrix refX = DenseMatrix::Random(rows, rows);
+    DenseMatrix refUp = DenseMatrix::Zero(rows, rows);
+    DenseMatrix refLo = DenseMatrix::Zero(rows, rows);
+    DenseMatrix refS = DenseMatrix::Zero(rows, rows);
+    SparseMatrix<Scalar> mUp(rows, rows);
+    SparseMatrix<Scalar> mLo(rows, rows);
+    SparseMatrix<Scalar> mS(rows, rows);
+    do {
+      initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular);
+    } while (refUp.isZero());
+    refLo = refUp.transpose().conjugate();
+    mLo = mUp.transpose().conjugate();
+    refS = refUp + refLo;
+    refS.diagonal() *= 0.5;
+    mS = mUp + mLo;
+    for (int k=0; k<mS.outerSize(); ++k)
+      for (typename SparseMatrix<Scalar>::InnerIterator it(mS,k); it; ++it)
+        if (it.index() == k)
+          it.valueRef() *= 0.5;
+    
+    VERIFY_IS_APPROX(refS.adjoint(), refS);
+    VERIFY_IS_APPROX(mS.transpose().conjugate(), mS);
+    VERIFY_IS_APPROX(mS, refS);
+    VERIFY_IS_APPROX(x=mS*b, refX=refS*b);
+    VERIFY_IS_APPROX(x=mUp.template marked<UpperTriangular|SelfAdjoint>()*b, refX=refS*b);
+    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_basic()