diff --git a/unsupported/Eigen/src/NonLinear/lmder.h b/unsupported/Eigen/src/NonLinear/lmder.h
index d279168aa..7ae2d82bf 100644
--- a/unsupported/Eigen/src/NonLinear/lmder.h
+++ b/unsupported/Eigen/src/NonLinear/lmder.h
@@ -198,9 +198,7 @@ L200:
 
     /*           determine the levenberg-marquardt parameter. */
 
-    ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
     ei_lmpar<Scalar>(fjac, ipvt, diag, qtf, delta, par, wa1, wa2);
-    ipvt.cwise()-=1;
 
     /*           store the direction p and x + p. calculate the norm of p. */
 
diff --git a/unsupported/Eigen/src/NonLinear/lmdif.h b/unsupported/Eigen/src/NonLinear/lmdif.h
index 7a003887f..67244f306 100644
--- a/unsupported/Eigen/src/NonLinear/lmdif.h
+++ b/unsupported/Eigen/src/NonLinear/lmdif.h
@@ -194,9 +194,7 @@ L200:
 
     /*           determine the levenberg-marquardt parameter. */
 
-    ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
     ei_lmpar<Scalar>(fjac, ipvt, diag, qtf, delta, par, wa1, wa2);
-    ipvt.cwise()-=1;
 
     /*           store the direction p and x + p. calculate the norm of p. */
 
diff --git a/unsupported/Eigen/src/NonLinear/lmpar.h b/unsupported/Eigen/src/NonLinear/lmpar.h
index 28fcf390b..664c78bff 100644
--- a/unsupported/Eigen/src/NonLinear/lmpar.h
+++ b/unsupported/Eigen/src/NonLinear/lmpar.h
@@ -2,7 +2,7 @@
 template <typename Scalar>
 void ei_lmpar(
         Matrix< Scalar, Dynamic, Dynamic > &r,
-        const VectorXi &ipvt,
+        VectorXi &ipvt, // TODO : const once ipvt mess fixed
         const Matrix< Scalar, Dynamic, 1 >  &diag,
         const Matrix< Scalar, Dynamic, 1 >  &qtb,
         Scalar delta,
@@ -50,7 +50,7 @@ void ei_lmpar(
     }
 
     for (j = 0; j < n; ++j) {
-        l = ipvt[j]-1;
+        l = ipvt[j];
         x[l] = wa1[j];
     }
 
@@ -73,7 +73,7 @@ void ei_lmpar(
     if (nsing < n-1)
         goto L120;
     for (j = 0; j < n; ++j) {
-        l = ipvt[j]-1;
+        l = ipvt[j];
         wa1[j] = diag[l] * (wa2[l] / dxnorm);
     }
     for (j = 0; j < n; ++j) {
@@ -92,7 +92,7 @@ L120:
         sum = 0.;
         for (i = 0; i <= j; ++i)
             sum += r(i,j) * qtb[i];
-        l = ipvt[j]-1;
+        l = ipvt[j];
         wa1[j] = sum / diag[l];
     }
     gnorm = wa1.stableNorm();
@@ -122,7 +122,9 @@ L150:
     temp = ei_sqrt(par);
     wa1 = temp * diag;
 
+    ipvt.cwise()+=1; // qrsolv() expects the fortran convention (as qrfac provides)
     ei_qrsolv<Scalar>(n, r.data(), r.rows(), ipvt.data(), wa1.data(), qtb.data(), x.data(), sdiag.data(), wa2.data());
+    ipvt.cwise()-=1;
 
     wa2 = diag.cwise() * x;
     dxnorm = wa2.blueNorm();
@@ -141,7 +143,7 @@ L150:
     /*        compute the newton correction. */
 
     for (j = 0; j < n; ++j) {
-        l = ipvt[j]-1;
+        l = ipvt[j];
         wa1[j] = diag[l] * (wa2[l] / dxnorm);
         /* L180: */
     }
diff --git a/unsupported/Eigen/src/NonLinear/lmstr.h b/unsupported/Eigen/src/NonLinear/lmstr.h
index 45593a505..783217a54 100644
--- a/unsupported/Eigen/src/NonLinear/lmstr.h
+++ b/unsupported/Eigen/src/NonLinear/lmstr.h
@@ -216,9 +216,7 @@ L240:
 
     /*           determine the levenberg-marquardt parameter. */
 
-    ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
     ei_lmpar<Scalar>(fjac, ipvt, diag, qtf, delta, par, wa1, wa2);
-    ipvt.cwise()-=1;
 
     /*           store the direction p and x + p. calculate the norm of p. */