Simplify Rank2Update.h

1 files changed, 14 insertions, 54 deletions

diff --git a/blas/Rank2Update.h b/blas/Rank2Update.h
index 2767916e2..8e1a676ee 100644
--- a/blas/Rank2Update.h
+++ b/blas/Rank2Update.h
@@ -16,34 +16,16 @@ namespace internal {
  * This is the low-level version of SelfadjointRank2Update.h
  */
 template<typename Scalar, typename Index, int UpLo>
-struct rank2_update_selector;
-
-template<typename Scalar, typename Index>
-struct rank2_update_selector<Scalar,Index,Upper>
-{
-  static void run(Index size, Scalar* mat, Index stride, const Scalar* _u, const Scalar* _v, Scalar alpha)
-  {
-    Map<const Matrix<Scalar,Dynamic,1> > u(_u, size), v(_v, size);
-    for (Index i=0; i<size; ++i)
-    {
-      Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i, i+1) +=
-			conj(alpha) * conj(_u[i]) * v.head(i+1)
-		      + alpha * conj(_v[i]) * u.head(i+1);
-    }
-  }
-};
-
-template<typename Scalar, typename Index>
-struct rank2_update_selector<Scalar,Index,Lower>
+struct rank2_update_selector
 {
-  static void run(Index size, Scalar* mat, Index stride, const Scalar* _u, const Scalar* _v, Scalar alpha)
+  static void run(Index size, Scalar* mat, Index stride, const Scalar* u, const Scalar* v, Scalar alpha)
   {
-    Map<const Matrix<Scalar,Dynamic,1> > u(_u, size), v(_v, size);
+    typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
     for (Index i=0; i<size; ++i)
     {
-      Map<Matrix<Scalar,Dynamic,1> >(mat+(stride+1)*i, size-i) +=
-			conj(alpha) * conj(_u[i]) * v.tail(size-i)
-		      + alpha * conj(_v[i]) * u.tail(size-i);
+      Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1)) +=
+			conj(alpha) * conj(u[i]) * OtherMap(v+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1))
+		      + alpha * conj(v[i]) * OtherMap(u+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1));
     }
   }
 };
@@ -52,42 +34,20 @@ struct rank2_update_selector<Scalar,Index,Lower>
  * The matrix is in packed form.
  */
 template<typename Scalar, typename Index, int UpLo>
-struct packed_rank2_update_selector;
-
-template<typename Scalar, typename Index>
-struct packed_rank2_update_selector<Scalar,Index,Upper>
-{
-  static void run(Index size, Scalar* mat, const Scalar* _u, const Scalar* _v, Scalar alpha)
-  {
-    Map<const Matrix<Scalar,Dynamic,1> > u(_u, size), v(_v, size);
-    Index offset = 0;
-    for (Index i=0; i<size; ++i)
-    {
-      Map<Matrix<Scalar,Dynamic,1> >(mat+offset, i+1) +=
-			conj(alpha) * conj(_u[i]) * v.head(i+1)
-		      + alpha * conj(_v[i]) * u.head(i+1);
-      //FIXME This should be handled outside.
-      mat[offset+i] = real(mat[offset+i]);
-      offset += i+1;
-    }
-  }
-};
-
-template<typename Scalar, typename Index>
-struct packed_rank2_update_selector<Scalar,Index,Lower>
+struct packed_rank2_update_selector
 {
-  static void run(Index size, Scalar* mat, const Scalar* _u, const Scalar* _v, Scalar alpha)
+  static void run(Index size, Scalar* mat, const Scalar* u, const Scalar* v, Scalar alpha)
   {
-    Map<const Matrix<Scalar,Dynamic,1> > u(_u, size), v(_v, size);
+    typedef Map<const Matrix<Scalar,Dynamic,1> > OtherMap;
     Index offset = 0;
     for (Index i=0; i<size; ++i)
     {
-      Map<Matrix<Scalar,Dynamic,1> >(mat+offset, size-i) +=
-			conj(alpha) * conj(_u[i]) * v.tail(size-i)
-		      + alpha * conj(_v[i]) * u.tail(size-i);
+      Map<Matrix<Scalar,Dynamic,1> >(mat+offset, UpLo==Lower ? size-i : (i+1)) +=
+			conj(alpha) * conj(u[i]) * OtherMap(v+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1))
+		      + alpha * conj(v[i]) * OtherMap(u+(UpLo==Lower ? i : 0), UpLo==Lower ? size-i : (i+1));
       //FIXME This should be handled outside.
-      mat[offset] = real(mat[offset]);
-      offset += size-i;
+      mat[offset+(UpLo==Lower ? 0 : i)] = real(mat[offset+(UpLo==Lower ? 0 : i)]);
+      offset += UpLo==Lower ? size-i : (i+1);
     }
   }
 };