1 files changed, 9 insertions, 9 deletions

diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LMpar.h b/unsupported/Eigen/src/LevenbergMarquardt/LMpar.h
index af76a9799..9532042d9 100644
--- a/unsupported/Eigen/src/LevenbergMarquardt/LMpar.h
+++ b/unsupported/Eigen/src/LevenbergMarquardt/LMpar.h
@@ -40,11 +40,15 @@ namespace internal {
     Scalar temp, paru;
     Scalar gnorm;
     Scalar dxnorm;
-
+    
+    // Make a copy of the triangular factor. 
+    // This copy is modified during call the qrsolv
+    MatrixType s;
+    s = qr.matrixR();
 
     /* Function Body */
     const Scalar dwarf = (std::numeric_limits<Scalar>::min)();
-    const Index n = qr.matrixQR().cols();
+    const Index n = qr.matrixR().cols();
     eigen_assert(n==diag.size());
     eigen_assert(n==qtb.size());
 
@@ -58,8 +62,7 @@ namespace internal {
     wa1 = qtb;
     wa1.tail(n-rank).setZero();
     //FIXME There is no solve in place for sparse triangularView
-    //qr.matrixQR().topLeftCorner(rank, rank).template triangularView<Upper>().solveInPlace(wa1.head(rank));
-    wa1.head(rank) = qr.matrixQR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(qtb.head(rank));
+    wa1.head(rank) = s.topLeftCorner(rank,rank).template triangularView<Upper>().solve(qtb.head(rank));
 
     x = qr.colsPermutation()*wa1;
 
@@ -81,14 +84,14 @@ namespace internal {
     parl = 0.;
     if (rank==n) {
       wa1 = qr.colsPermutation().inverse() *  diag.cwiseProduct(wa2)/dxnorm;
-      qr.matrixQR().topLeftCorner(n, n).transpose().template triangularView<Lower>().solveInPlace(wa1);
+      s.topLeftCorner(n,n).transpose().template triangularView<Lower>().solveInPlace(wa1);
       temp = wa1.blueNorm();
       parl = fp / m_delta / temp / temp;
     }
 
     /* calculate an upper bound, paru, for the zero of the function. */
     for (j = 0; j < n; ++j)
-      wa1[j] = qr.matrixQR().col(j).head(j+1).dot(qtb.head(j+1)) / diag[qr.colsPermutation().indices()(j)];
+      wa1[j] = s.col(j).head(j+1).dot(qtb.head(j+1)) / diag[qr.colsPermutation().indices()(j)];
 
     gnorm = wa1.stableNorm();
     paru = gnorm / m_delta;
@@ -103,8 +106,6 @@ namespace internal {
       par = gnorm / dxnorm;
 
     /* beginning of an iteration. */
-    MatrixType s;
-    s = qr.matrixQR();
     while (true) {
       ++iter;
 
@@ -130,7 +131,6 @@ namespace internal {
       /* compute the newton correction. */
       wa1 = qr.colsPermutation().inverse() * diag.cwiseProduct(wa2/dxnorm);
       // we could almost use this here, but the diagonal is outside qr, in sdiag[]
-      // qr.matrixQR().topLeftCorner(n, n).transpose().template triangularView<Lower>().solveInPlace(wa1);
       for (j = 0; j < n; ++j) {
         wa1[j] /= sdiag[j];
         temp = wa1[j];