From 5546e963c922aeafd5a9f1abc47a3f0ff9081ff8 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Sat, 16 Aug 2008 07:09:39 +0000
Subject: [PATCH] bugfix in CholeskyWithoutSquareRoot::solve found by Timothy
 Hunter

---
 Eigen/Cholesky                                 | 2 ++
 Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h | 7 +++----
 test/cholesky.cpp                              | 9 ++++++---
 3 files changed, 11 insertions(+), 7 deletions(-)

diff --git a/Eigen/Cholesky b/Eigen/Cholesky
index a0d4311f9..f5bcec52d 100644
--- a/Eigen/Cholesky
+++ b/Eigen/Cholesky
@@ -25,6 +25,8 @@ namespace Eigen {
   * \endcode
   */
 
+#include "src/Array/CwiseOperators.h"
+#include "src/Array/Functors.h"
 #include "src/Cholesky/Cholesky.h"
 #include "src/Cholesky/CholeskyWithoutSquareRoot.h"
 
diff --git a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
index 69dcdb8a2..6cf5d84c2 100644
--- a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
+++ b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
@@ -149,10 +149,9 @@ typename Derived::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(const Matrix
 
   return m_matrix.adjoint().template part<UnitUpper>()
     .solveTriangular(
-      (matrixL()
-        .solveTriangular(b))
-        .cwise()/m_matrix.diagonal()
-      );
+      (  m_matrix.cwise().inverse().diagonal().asDiagonal()
+       * matrixL().solveTriangular(b))
+     );
 }
 
 /** \cholesky_module
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index eab1febf2..55a74c477 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -39,16 +39,19 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
 
   MatrixType a = MatrixType::Random(rows,cols);
-  VectorType b = VectorType::Random(rows);
+  VectorType vecB = VectorType::Random(rows);
+  MatrixType matB = MatrixType::Random(rows,cols);
   SquareMatrixType covMat =  a * a.adjoint();
 
   CholeskyWithoutSquareRoot<SquareMatrixType> cholnosqrt(covMat);
   VERIFY_IS_APPROX(covMat, cholnosqrt.matrixL() * cholnosqrt.vectorD().asDiagonal() * cholnosqrt.matrixL().adjoint());
-  VERIFY_IS_APPROX(covMat * cholnosqrt.solve(b), b);
+  VERIFY_IS_APPROX(covMat * cholnosqrt.solve(vecB), vecB);
+  VERIFY_IS_APPROX(covMat * cholnosqrt.solve(matB), matB);
 
   Cholesky<SquareMatrixType> chol(covMat);
   VERIFY_IS_APPROX(covMat, chol.matrixL() * chol.matrixL().adjoint());
-  VERIFY_IS_APPROX(covMat * chol.solve(b), b);
+  VERIFY_IS_APPROX(covMat * chol.solve(vecB), vecB);
+  VERIFY_IS_APPROX(covMat * chol.solve(matB), matB);
 }
 
 void test_cholesky()