some cleaning in Cholesky and removed evil ei_sqrt of complex

Eigen/src/Cholesky/Cholesky.h

@@ -32,12 +32,15 @@
   * \param MatrixType the type of the matrix of which we are computing the Cholesky decomposition
   *
   * This class performs a standard Cholesky decomposition of a symmetric, positive definite
-  * matrix A such that A = U'U = LL', where U is upper triangular.
+  * matrix A such that A = LL^* = U^*U, where L is lower triangular.
   *
-  * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like  A'A x = b,
+  * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like  D^*D x = b,
   * for that purpose, we recommend the Cholesky decomposition without square root which is more stable
   * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
-  * situation like generalised eigen problem with hermitian matrices.
+  * situations like generalised eigen problems with hermitian matrices.
+  *
+  * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
+  * the strict lower part does not have to store correct values.
   *
   * \sa class CholeskyWithoutSquareRoot
   */
@@ -46,6 +49,7 @@ template<typename MatrixType> class Cholesky
   public:
 
     typedef typename MatrixType::Scalar Scalar;
+    typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
 
     Cholesky(const MatrixType& matrix)
@@ -61,75 +65,60 @@ template<typename MatrixType> class Cholesky
 
     bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
 
-    template<typename DerivedVec>
+    template<typename Derived>
-    typename DerivedVec::Eval solve(MatrixBase<DerivedVec> &vecB);
+    typename Derived::Eval solve(MatrixBase<Derived> &b);
 
-    /** Compute / recompute the Cholesky decomposition A = U'U = LL' of \a matrix
-      */
     void compute(const MatrixType& matrix);
 
   protected:
     /** \internal
-      * Used to compute and store the cholesky decomposition.
+      * Used to compute and store L
-      * The strict upper part correspond to the coefficients of the input
+      * The strict upper part is not used and even not initialized.
-      * symmetric matrix, while the lower part store U'=L.
       */
     MatrixType m_matrix;
     bool m_isPositiveDefinite;
 };
 
-template<typename MatrixType>
+/** Compute / recompute the Cholesky decomposition A = LL^* = U^*U of \a matrix
-void Cholesky<MatrixType>::compute(const MatrixType& matrix)
-{
-  assert(matrix.rows()==matrix.cols());
-  const int size = matrix.rows();
-  m_matrix = matrix.conjugate();
-
-  #if 0
-  m_isPositiveDefinite = true;
-  for (int i = 0; i < size; ++i)
-  {
-    m_isPositiveDefinite = m_isPositiveDefinite && m_matrix(i,i) > Scalar(0);
-    m_matrix(i,i) = ei_sqrt(m_matrix(i,i));
-    if (i+1<size)
-      m_matrix.col(i).end(size-i-1) /= m_matrix(i,i);
-    for (int j = i+1; j < size; ++j)
-    {
-      m_matrix.col(j).end(size-j) -= m_matrix(j,i) * m_matrix.col(i).end(size-j);
-    }
-  }
-  #else
-  // this version looks faster for large matrices
-//   m_isPositiveDefinite = m_matrix(0,0) > Scalar(0);
-  m_matrix(0,0) = ei_sqrt(m_matrix(0,0));
-  m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix(0,0);
-  for (int j = 1; j < size; ++j)
-  {
-//     Scalar tmp = m_matrix(j,j) - m_matrix.row(j).start(j).norm2();
-    Scalar tmp = m_matrix(j,j) - (m_matrix.row(j).start(j) * m_matrix.row(j).start(j).adjoint())(0,0);
-//     m_isPositiveDefinite = m_isPositiveDefinite && tmp > Scalar(0);
-//     m_matrix(j,j) = ei_sqrt(tmp<Scalar(0) ? Scalar(0) : tmp);
-    m_matrix(j,j) = ei_sqrt(tmp);
-    tmp = 1. / m_matrix(j,j);
-    for (int i = j+1; i < size; ++i)
-      m_matrix(i,j) = tmp * (m_matrix(j,i) -
-          (m_matrix.row(i).start(j) * m_matrix.row(j).start(j).adjoint())(0,0) );
-  }
-  #endif
-}
-
-/** Solve A*x = b with A symmeric positive definite using the available Cholesky decomposition.
   */
 template<typename MatrixType>
-template<typename DerivedVec>
+void Cholesky<MatrixType>::compute(const MatrixType& a)
-typename DerivedVec::Eval Cholesky<MatrixType>::solve(MatrixBase<DerivedVec> &vecB)
+{
+  assert(a.rows()==a.cols());
+  const int size = a.rows();
+  m_matrix.resize(size, size);
+
+  RealScalar x;
+  x = ei_real(a.coeff(0,0));
+  m_isPositiveDefinite = x > RealScalar(0) && ei_isMuchSmallerThan(ei_imag(m_matrix.coeff(0,0)), RealScalar(1));
+  m_matrix.coeffRef(0,0) = ei_sqrt(x);
+  m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / m_matrix.coeff(0,0);
+  for (int j = 1; j < size; ++j)
+  {
+    Scalar tmp = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).norm2();
+    x = ei_real(tmp);
+    m_isPositiveDefinite = m_isPositiveDefinite && x > RealScalar(0) && ei_isMuchSmallerThan(ei_imag(tmp), RealScalar(1));
+    m_matrix.coeffRef(j,j) = x = ei_sqrt(x);
+
+    int endSize = size-j-1;
+    if (endSize>0)
+      m_matrix.col(j).end(endSize) = (a.row(j).end(endSize).adjoint()
+        - m_matrix.block(j+1, 0, endSize, j) * m_matrix.row(j).start(j).adjoint()) / x;
+  }
+}
+
+/** \returns the solution of A x = \a b using the current decomposition of A.
+  * In other words, it returns \code A^-1 b \endcode computing
+  * \code L^-*  L^1 b \code from right to left.
+  */
+template<typename MatrixType>
+template<typename Derived>
+typename Derived::Eval Cholesky<MatrixType>::solve(MatrixBase<Derived> &b)
 {
   const int size = m_matrix.rows();
-  ei_assert(size==vecB.size());
+  ei_assert(size==b.size());
 
-  // FIXME .inverseProduct creates a temporary that is not nice since it is called twice
+  return m_matrix.adjoint().upper().inverseProduct(m_matrix.lower().inverseProduct(b));
-  // add a .inverseProductInPlace ??
-  return m_matrix.adjoint().upper().inverseProduct(m_matrix.lower().inverseProduct(vecB));
 }
 
 

Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h

@@ -32,13 +32,14 @@
   * \param MatrixType the type of the matrix of which we are computing the Cholesky decomposition
   *
   * This class performs a Cholesky decomposition without square root of a symmetric, positive definite
-  * matrix A such that A = U' D U = L D L', where U is upper triangular with a unit diagonal and D is a diagonal
+  * matrix A such that A = L D L^* = U^* D U, where L is lower triangular with a unit diagonal and D is a diagonal
   * matrix.
   *
   * Compared to a standard Cholesky decomposition, avoiding the square roots allows for faster and more
   * stable computation.
   *
-  * \todo what about complex matrices ?
+  * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
+  * the strict lower part does not have to store correct values.
   *
   * \sa class Cholesky
   */
@@ -47,6 +48,7 @@ template<typename MatrixType> class CholeskyWithoutSquareRoot
   public:
 
     typedef typename MatrixType::Scalar Scalar;
+    typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
 
     CholeskyWithoutSquareRoot(const MatrixType& matrix)
@@ -55,92 +57,85 @@ template<typename MatrixType> class CholeskyWithoutSquareRoot
       compute(matrix);
     }
 
+    /** \returns the lower triangular matrix L */
     Triangular<Lower|UnitDiagBit, MatrixType > matrixL(void) const
     {
       return m_matrix.lowerWithUnitDiag();
     }
 
+    /** \returns the coefficients of the diagonal matrix D */
     DiagonalCoeffs<MatrixType> vectorD(void) const
     {
       return m_matrix.diagonal();
     }
 
+    /** \returns whether the matrix is positive definite */
     bool isPositiveDefinite(void) const
     {
       return m_matrix.diagonal().minCoeff() > Scalar(0);
     }
 
-    template<typename DerivedVec>
+    template<typename Derived>
-    typename DerivedVec::Eval solve(MatrixBase<DerivedVec> &vecB);
+    typename Derived::Eval solve(MatrixBase<Derived> &b);
 
-    /** Compute / recompute the Cholesky decomposition A = U'DU = LDL' of \a matrix
-      */
     void compute(const MatrixType& matrix);
 
   protected:
     /** \internal
-      * Used to compute and store the cholesky decomposition A = U'DU = LDL'.
+      * Used to compute and store the cholesky decomposition A = L D L^* = U^* D U.
       * The strict upper part is used during the decomposition, the strict lower
-      * part correspond to the coefficients of U'=L (its diagonal is equal to 1 and
+      * part correspond to the coefficients of L (its diagonal is equal to 1 and
       * is not stored), and the diagonal entries correspond to D.
       */
     MatrixType m_matrix;
 };
 
-template<typename MatrixType>
+/** Compute / recompute the Cholesky decomposition A = L D L^* = U^* D U of \a matrix
-void CholeskyWithoutSquareRoot<MatrixType>::compute(const MatrixType& matrix)
-{
-  assert(matrix.rows()==matrix.cols());
-  const int size = matrix.rows();
-  m_matrix = matrix.conjugate();
-  #if 0
-  for (int i = 0; i < size; ++i)
-  {
-    Scalar tmp = Scalar(1) / m_matrix(i,i);
-    for (int j = i+1; j < size; ++j)
-    {
-      m_matrix(j,i) = m_matrix(i,j) * tmp;
-      m_matrix.row(j).end(size-j) -= m_matrix(j,i) * m_matrix.row(i).end(size-j);
-    }
-  }
-  #else
-  // this version looks faster for large matrices
-  m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix(0,0);
-  for (int j = 1; j < size; ++j)
-  {
-    Scalar tmp = m_matrix(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate())(0,0);
-    m_matrix(j,j) = tmp;
-    tmp = Scalar(1) / tmp;
-    for (int i = j+1; i < size; ++i)
-    {
-      m_matrix(j,i) = (m_matrix(j,i) - (m_matrix.row(i).start(j) * m_matrix.col(j).start(j).conjugate())(0,0) );
-      m_matrix(i,j) = tmp * m_matrix(j,i);
-    }
-  }
-  #endif
-}
-
-/** Solve A*x = b with A symmeric positive definite using the available Cholesky decomposition.
   */
 template<typename MatrixType>
-template<typename DerivedVec>
+void CholeskyWithoutSquareRoot<MatrixType>::compute(const MatrixType& a)
-typename DerivedVec::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(MatrixBase<DerivedVec> &vecB)
+{
+  assert(a.rows()==a.cols());
+  const int size = a.rows();
+  m_matrix.resize(size, size);
+
+  // Note that, in this algorithm the rows of the strict upper part of m_matrix is used to store
+  // column vector, thus the strange .conjugate() and .transpose()...
+
+  m_matrix.row(0) = a.row(0).conjugate();
+  m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix.coeff(0,0);
+  for (int j = 1; j < size; ++j)
+  {
+    RealScalar tmp = ei_real(a.coeff(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate()).coeff(0,0));
+    m_matrix.coeffRef(j,j) = tmp;
+
+    int endSize = size-j-1;
+    if (endSize>0)
+    {
+      m_matrix.row(j).end(endSize) = a.row(j).end(endSize).conjugate()
+        - (m_matrix.block(j+1,0, endSize, j) * m_matrix.col(j).start(j).conjugate()).transpose();
+      m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp;
+    }
+  }
+}
+
+/** \returns the solution of A x = \a b using the current decomposition of A.
+  * In other words, it returns \code A^-1 b \endcode computing
+  * \code (L^-*) (D^-1) (L^-1) b \code from right to left.
+  */
+template<typename MatrixType>
+template<typename Derived>
+typename Derived::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(MatrixBase<Derived> &vecB)
 {
   const int size = m_matrix.rows();
   ei_assert(size==vecB.size());
 
-  // FIXME .inverseProduct creates a temporary that is not nice since it is called twice
-  // maybe add a .inverseProductInPlace() ??
   return m_matrix.adjoint().upperWithUnitDiag()
     .inverseProduct(
       (m_matrix.lowerWithUnitDiag()
         .inverseProduct(vecB))
         .cwiseQuotient(m_matrix.diagonal())
       );
-
-//   return m_matrix.adjoint().upperWithUnitDiag()
-//     .inverseProduct(
-//       (m_matrix.lowerWithUnitDiag() * (m_matrix.diagonal().asDiagonal())).lower().inverseProduct(vecB));
 }
 
 

Eigen/src/Core/MathFunctions.h

@@ -181,43 +181,6 @@ inline std::complex<double> ei_exp(std::complex<double> x)  { return std::exp(x)
 inline std::complex<double> ei_sin(std::complex<double> x)  { return std::sin(x); }
 inline std::complex<double> ei_cos(std::complex<double> x)  { return std::cos(x); }
 
-template<typename T>
-inline std::complex<T> ei_sqrt(const std::complex<T>& x)
-{
-  if (std::real(x) == 0.0 && std::imag(x) == 0.0)
-    return std::complex<T>(0);
-  else
-  {
-    T a = ei_abs(std::real(x));
-    T b = ei_abs(std::imag(x));
-    T c;
-
-    if (a >= b)
-    {
-      T t = b / a;
-      c = ei_sqrt(a) * ei_sqrt(0.5 * (1.0 + ei_sqrt(1.0 + t * t)));
-    }
-    else
-    {
-      T t = a / b;
-      c = ei_sqrt(b) * ei_sqrt(0.5 * (t + ei_sqrt (1.0 + t * t)));
-    }
-
-    T d = std::imag(x) / (2.0 * c);
-    if (std::real(x) >= 0.0)
-    {
-      return std::complex<T>(c, d);
-    }
-    else
-    {
-      std::complex<T> res(d, c);
-      if (std::imag(x)<0.0)
-        res = -res;
-      return res;
-    }
-  }
-}
-
 template<> inline std::complex<double> ei_random()
 {
   return std::complex<double>(ei_random<double>(), ei_random<double>());

Eigen/src/LU/Inverse.h

@@ -98,7 +98,7 @@ template<typename MatrixType, bool CheckExistence> class Inverse : ei_no_assignm
 
   protected:
     bool m_exists;
-    MatrixType m_inverse;
+    typename MatrixType::Eval m_inverse;
 };
 
 template<typename MatrixType, bool CheckExistence>