From fc85f91df06699cbc643409df4c3aacaabc6f484 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Tue, 28 Feb 2012 16:19:40 +0100
Subject: [PATCH] fix MKL interface with LLT::rankUpdate

---
 Eigen/src/Cholesky/LLT.h     | 142 ++++++++++++++++++-----------------
 Eigen/src/Cholesky/LLT_MKL.h |  35 ++-------
 2 files changed, 80 insertions(+), 97 deletions(-)

diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 0ce496c35..dd705c57e 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -184,7 +184,7 @@ template<typename _MatrixType, int _UpLo> class LLT
     inline Index cols() const { return m_matrix.cols(); }
 
     template<typename VectorType>
-    LLT& rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
+    LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
 
   protected:
     /** \internal
@@ -200,6 +200,76 @@ namespace internal {
 
 template<typename Scalar, int UpLo> struct llt_inplace;
 
+template<typename MatrixType, typename VectorType>
+static typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma)
+{
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef typename MatrixType::Index Index;
+  typedef typename MatrixType::ColXpr ColXpr;
+  typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
+  typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
+  typedef Matrix<Scalar,Dynamic,1> TempVectorType;
+  typedef typename TempVectorType::SegmentReturnType TempVecSegment;
+
+  int n = mat.cols();
+  eigen_assert(mat.rows()==n && vec.size()==n);
+
+  TempVectorType temp;
+
+  if(sigma>0)
+  {
+    // This version is based on Givens rotations.
+    // It is faster than the other one below, but only works for updates,
+    // i.e., for sigma > 0
+    temp = sqrt(sigma) * vec;
+
+    for(int i=0; i<n; ++i)
+    {
+      JacobiRotation<Scalar> g;
+      g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
+
+      int rs = n-i-1;
+      if(rs>0)
+      {
+        ColXprSegment x(mat.col(i).tail(rs));
+        TempVecSegment y(temp.tail(rs));
+        apply_rotation_in_the_plane(x, y, g);
+      }
+    }
+  }
+  else
+  {
+    temp = vec;
+    RealScalar beta = 1;
+    for(int j=0; j<n; ++j)
+    {
+      RealScalar Ljj = real(mat.coeff(j,j));
+      RealScalar dj = abs2(Ljj);
+      Scalar wj = temp.coeff(j);
+      RealScalar swj2 = sigma*abs2(wj);
+      RealScalar gamma = dj*beta + swj2;
+
+      RealScalar x = dj + swj2/beta;
+      if (x<=RealScalar(0))
+        return j;
+      RealScalar nLjj = sqrt(x);
+      mat.coeffRef(j,j) = nLjj;
+      beta += swj2/dj;
+
+      // Update the terms of L
+      Index rs = n-j-1;
+      if(rs)
+      {
+        temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
+        if(gamma != 0)
+          mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs);
+      }
+    }
+  }
+  return -1;
+}
+
 template<typename Scalar> struct llt_inplace<Scalar, Lower>
 {
   typedef typename NumTraits<Scalar>::Real RealScalar;
@@ -265,72 +335,10 @@ template<typename Scalar> struct llt_inplace<Scalar, Lower>
   template<typename MatrixType, typename VectorType>
   static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
   {
-    typedef typename MatrixType::Index Index;
-    typedef typename MatrixType::ColXpr ColXpr;
-    typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
-    typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
-    typedef Matrix<Scalar,Dynamic,1> TempVectorType;
-    typedef typename TempVectorType::SegmentReturnType TempVecSegment;
-
-    int n = mat.cols();
-    eigen_assert(mat.rows()==n && vec.size()==n);
-
-    TempVectorType temp;
-
-    if(sigma>0)
-    {
-      // This version is based on Givens rotations.
-      // It is faster than the other one below, but only works for updates,
-      // i.e., for sigma > 0
-      temp = sqrt(sigma) * vec;
-
-      for(int i=0; i<n; ++i)
-      {
-        JacobiRotation<Scalar> g;
-        g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
-
-        int rs = n-i-1;
-        if(rs>0)
-        {
-          ColXprSegment x(mat.col(i).tail(rs));
-          TempVecSegment y(temp.tail(rs));
-          apply_rotation_in_the_plane(x, y, g);
-        }
-      }
-    }
-    else
-    {
-      temp = vec;
-      RealScalar beta = 1;
-      for(int j=0; j<n; ++j)
-      {
-        RealScalar Ljj = real(mat.coeff(j,j));
-        RealScalar dj = abs2(Ljj);
-        Scalar wj = temp.coeff(j);
-        RealScalar swj2 = sigma*abs2(wj);
-        RealScalar gamma = dj*beta + swj2;
-
-        RealScalar x = dj + swj2/beta;
-        if (x<=RealScalar(0))
-          return j;
-        RealScalar nLjj = sqrt(x);
-        mat.coeffRef(j,j) = nLjj;
-        beta += swj2/dj;
-
-        // Update the terms of L
-        Index rs = n-j-1;
-        if(rs)
-        {
-          temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
-          if(gamma != 0)
-            mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs);
-        }
-      }
-    }
-    return -1;
+    return llt_rank_update_lower(mat, vec, sigma);
   }
 };
-
+  
 template<typename Scalar> struct llt_inplace<Scalar, Upper>
 {
   typedef typename NumTraits<Scalar>::Real RealScalar;
@@ -404,9 +412,9 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
   * then after it we have LL^* = A + sigma * v v^* where \a v must be a vector
   * of same dimension.
   */
-template<typename MatrixType, int _UpLo>
+template<typename _MatrixType, int _UpLo>
 template<typename VectorType>
-LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
+LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
 {
   EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
   eigen_assert(v.size()==m_matrix.cols());
diff --git a/Eigen/src/Cholesky/LLT_MKL.h b/Eigen/src/Cholesky/LLT_MKL.h
index c5ed77d33..10aa746dc 100644
--- a/Eigen/src/Cholesky/LLT_MKL.h
+++ b/Eigen/src/Cholesky/LLT_MKL.h
@@ -50,7 +50,7 @@ template<> struct mkl_llt<EIGTYPE> \
     lapack_int size, lda, info, StorageOrder; \
     EIGTYPE* a; \
     eigen_assert(m.rows()==m.cols()); \
-/* Set up parameters for ?potrf */ \
+    /* Set up parameters for ?potrf */ \
     size = m.rows(); \
     StorageOrder = MatrixType::Flags&RowMajorBit?RowMajor:ColMajor; \
     matrix_order = StorageOrder==RowMajor ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \
@@ -70,33 +70,8 @@ template<> struct llt_inplace<EIGTYPE, Lower> \
     return mkl_llt<EIGTYPE>::potrf(m, 'L'); \
   } \
   template<typename MatrixType, typename VectorType> \
-  static void rankUpdate(MatrixType& mat, const VectorType& vec) \
-  { \
-    typedef typename MatrixType::ColXpr ColXpr; \
-    typedef typename internal::remove_all<ColXpr>::type ColXprCleaned; \
-    typedef typename ColXprCleaned::SegmentReturnType ColXprSegment; \
-    typedef typename MatrixType::Scalar Scalar; \
-    typedef Matrix<Scalar,Dynamic,1> TempVectorType; \
-    typedef typename TempVectorType::SegmentReturnType TempVecSegment; \
-\
-    int n = mat.cols(); \
-    eigen_assert(mat.rows()==n && vec.size()==n); \
-    TempVectorType temp(vec); \
-\
-    for(int i=0; i<n; ++i) \
-    { \
-      JacobiRotation<Scalar> g; \
-      g.makeGivens(mat(i,i), -temp(i), &mat(i,i)); \
-\
-      int rs = n-i-1; \
-      if(rs>0) \
-      { \
-        ColXprSegment x(mat.col(i).tail(rs)); \
-        TempVecSegment y(temp.tail(rs)); \
-        apply_rotation_in_the_plane(x, y, g); \
-      } \
-    } \
-  } \
+  static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma) \
+  { return llt_rank_update_lower(mat, vec, sigma); } \
 }; \
 template<> struct llt_inplace<EIGTYPE, Upper> \
 { \
@@ -106,10 +81,10 @@ template<> struct llt_inplace<EIGTYPE, Upper> \
     return mkl_llt<EIGTYPE>::potrf(m, 'U'); \
   } \
   template<typename MatrixType, typename VectorType> \
-  static void rankUpdate(MatrixType& mat, const VectorType& vec) \
+  static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma) \
   { \
     Transpose<MatrixType> matt(mat); \
-    return llt_inplace<EIGTYPE, Lower>::rankUpdate(matt, vec.conjugate()); \
+    return llt_inplace<EIGTYPE, Lower>::rankUpdate(matt, vec.conjugate(), sigma); \
   } \
 };