Solve the issue found by Timothy in solveTriangular:

  => row-major rhs are now evaluated to a column-major
     temporary before the computations.
Add solveInPlace in Cholesky*

2 files changed, 67 insertions, 4 deletions

diff --git a/Eigen/src/Cholesky/Cholesky.h b/Eigen/src/Cholesky/Cholesky.h
index a64ab7c70..b94fea8dc 100644
--- a/Eigen/src/Cholesky/Cholesky.h
+++ b/Eigen/src/Cholesky/Cholesky.h
@@ -74,6 +74,9 @@ template<typename MatrixType> class Cholesky
     template<typename Derived>
     typename Derived::Eval solve(const MatrixBase<Derived> &b) const;
 
+    template<typename Derived>
+    bool solveInPlace(MatrixBase<Derived> &bAndX) const;
+
     void compute(const MatrixType& matrix);
 
   protected:
@@ -141,8 +144,37 @@ typename Derived::Eval Cholesky<MatrixType>::solve(const MatrixBase<Derived> &b)
 {
   const int size = m_matrix.rows();
   ei_assert(size==b.rows());
+  typename ei_eval_to_column_major<Derived>::type x(b);
+  solveInPlace(x);
+  return x;
+  //return m_matrix.adjoint().template part<Upper>().solveTriangular(matrixL().solveTriangular(b));
+}
 
-  return m_matrix.adjoint().template part<Upper>().solveTriangular(matrixL().solveTriangular(b));
+/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
+  * The result is stored in \a bAndx
+  *
+  * \returns true in case of success, false otherwise.
+  *
+  * In other words, it computes \f$ b = A^{-1} b \f$ with
+  * \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left.
+  * \param bAndX stores both the matrix \f$ b \f$ and the result \f$ x \f$
+  *
+  * Example: \include Cholesky_solve.cpp
+  * Output: \verbinclude Cholesky_solve.out
+  *
+  * \sa MatrixBase::cholesky(), Cholesky::solve()
+  */
+template<typename MatrixType>
+template<typename Derived>
+bool Cholesky<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
+{
+  const int size = m_matrix.rows();
+  ei_assert(size==bAndX.rows());
+  if (!m_isPositiveDefinite)
+    return false;
+  matrixL().solveTriangularInPlace(bAndX);
+  m_matrix.adjoint().template part<Upper>().solveTriangularInPlace(bAndX);
+  return true;
 }
 
 /** \cholesky_module
diff --git a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
index 4040869b0..fd111fb1f 100644
--- a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
+++ b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
@@ -71,6 +71,9 @@ template<typename MatrixType> class CholeskyWithoutSquareRoot
     template<typename Derived>
     typename Derived::Eval solve(const MatrixBase<Derived> &b) const;
 
+		template<typename Derived>
+    bool solveInPlace(MatrixBase<Derived> &bAndX) const;
+
     void compute(const MatrixType& matrix);
 
   protected:
@@ -101,7 +104,7 @@ void CholeskyWithoutSquareRoot<MatrixType>::compute(const MatrixType& a)
     m_matrix = a;
     return;
   }
-  
+
   // Let's preallocate a temporay vector to evaluate the matrix-vector product into it.
   // Unlike the standard Cholesky decomposition, here we cannot evaluate it to the destination
   // matrix because it a sub-row which is not compatible suitable for efficient packet evaluation.
@@ -144,8 +147,8 @@ void CholeskyWithoutSquareRoot<MatrixType>::compute(const MatrixType& a)
   * \param b the column vector \f$ b \f$, which can also be a matrix.
   *
   * See Cholesky::solve() for a example.
-  * 
-  * \sa MatrixBase::choleskyNoSqrt()
+  *
+  * \sa CholeskyWithoutSquareRoot::solveInPlace(), MatrixBase::choleskyNoSqrt()
   */
 template<typename MatrixType>
 template<typename Derived>
@@ -161,6 +164,34 @@ typename Derived::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(const Matrix
      );
 }
 
+/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
+  * The result is stored in \a bAndx
+  *
+  * \returns true in case of success, false otherwise.
+  *
+  * In other words, it computes \f$ b = A^{-1} b \f$ with
+  * \f$ {L^{*}}^{-1} D^{-1} L^{-1} b \f$ from right to left.
+  * \param bAndX stores both the matrix \f$ b \f$ and the result \f$ x \f$
+  *
+  * Example: \include Cholesky_solve.cpp
+  * Output: \verbinclude Cholesky_solve.out
+  *
+  * \sa MatrixBase::cholesky(), CholeskyWithoutSquareRoot::solve()
+  */
+template<typename MatrixType>
+template<typename Derived>
+bool CholeskyWithoutSquareRoot<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
+{
+  const int size = m_matrix.rows();
+  ei_assert(size==bAndX.rows());
+  if (!m_isPositiveDefinite)
+    return false;
+  matrixL().solveTriangularInPlace(bAndX);
+	bAndX *= m_matrix.cwise().inverse().template part<Diagonal>();
+  m_matrix.adjoint().template part<UnitUpper>().solveTriangularInPlace(bAndX);
+  return true;
+}
+
 /** \cholesky_module
   * \returns the Cholesky decomposition without square root of \c *this
   */