diff --git a/unsupported/Eigen/src/NonLinear/lmder.h b/unsupported/Eigen/src/NonLinear/lmder.h
index bb92f446a..4499b381e 100644
--- a/unsupported/Eigen/src/NonLinear/lmder.h
+++ b/unsupported/Eigen/src/NonLinear/lmder.h
@@ -105,6 +105,7 @@ L40:
     /*        compute the qr factorization of the jacobian. */
 
     qrfac(m, n, fjac.data(), ldfjac, TRUE_, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data());
+    ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1)
 
     /*        on the first iteration and if mode is 1, scale according */
     /*        to the norms of the columns of the initial jacobian. */
@@ -172,7 +173,7 @@ L120:
         goto L170;
     }
     for (j = 0; j < n; ++j) {
-        l = ipvt[j]-1;
+        l = ipvt[j];
         if (wa2[l] == 0.) {
             goto L150;
         }
@@ -213,8 +214,10 @@ L200:
 
     /*           determine the levenberg-marquardt parameter. */
 
+    ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
     lmpar(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta,
             &par, wa1.data(), wa2.data(), wa3.data(), wa4.data());
+    ipvt.cwise()-=1;
 
     /*           store the direction p and x + p. calculate the norm of p. */
 
@@ -252,7 +255,7 @@ L200:
 
     for (j = 0; j < n; ++j) {
         wa3[j] = 0.;
-        l = ipvt[j]-1;
+        l = ipvt[j];
         temp = wa1[l];
         for (i = 0; i <= j; ++i) {
             wa3[i] += fjac(i,j) * temp;
diff --git a/unsupported/Eigen/src/NonLinear/lmstr.h b/unsupported/Eigen/src/NonLinear/lmstr.h
index 1bf9caff3..a221b897c 100644
--- a/unsupported/Eigen/src/NonLinear/lmstr.h
+++ b/unsupported/Eigen/src/NonLinear/lmstr.h
@@ -129,7 +129,7 @@ L40:
         if (fjac(j,j) == 0.) {
             sing = TRUE_;
         }
-        ipvt[j] = j+1;
+        ipvt[j] = j;
         wa2[j] = fjac.col(j).start(j).stableNorm();
 //        wa2[j] = ei_enorm<Scalar>(j, &fjac[j * ldfjac + 1]);
 //            sum += fjac[i + j * ldfjac] * (qtf[i] / fnorm);
@@ -138,7 +138,9 @@ L40:
     if (! sing) {
         goto L130;
     }
+    ipvt.cwise()+=1;
     qrfac(n, n, fjac.data(), ldfjac, TRUE_, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data());
+    ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1)
     for (j = 0; j < n; ++j) {
         if (fjac(j,j) == 0.) {
             goto L110;
@@ -198,7 +200,7 @@ L170:
         goto L210;
     }
     for (j = 0; j < n; ++j) {
-        l = ipvt[j]-1;
+        l = ipvt[j];
         if (wa2[l] == 0.) {
             goto L190;
         }
@@ -239,8 +241,10 @@ L240:
 
     /*           determine the levenberg-marquardt parameter. */
 
+    ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
     lmpar(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta, &par,
             wa1.data(), wa2.data(), wa3.data(), wa4.data());
+    ipvt.cwise()-=1;
 
     /*           store the direction p and x + p. calculate the norm of p. */
 
@@ -278,7 +282,7 @@ L240:
 
     for (j = 0; j < n; ++j) {
         wa3[j] = 0.;
-        l = ipvt[j]-1;
+        l = ipvt[j];
         temp = wa1[l];
         for (i = 0; i <= j; ++i) {
             wa3[i] += fjac(i,j) * temp;
diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp
index bb16feb50..e7eeff58c 100644
--- a/unsupported/test/NonLinear.cpp
+++ b/unsupported/test/NonLinear.cpp
@@ -246,6 +246,7 @@ void testLmder()
   covar_ftol = epsilon<double>();
   covfac = fnorm*fnorm/(m-n);
   VectorXd wa(n);
+  ipvt.cwise()+=1; // covar() expects the fortran convention (as qrfac provides)
   covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data());
 
   MatrixXd cov_ref(n,n);
@@ -761,6 +762,7 @@ void testLmdif()
   covar_ftol = epsilon<double>();
   covfac = fnorm*fnorm/(m-n);
   VectorXd wa(n);
+//  ipvt.cwise()+=1; // covar() expects the fortran convention (as qrfac provides)
   covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data());
 
   MatrixXd cov_ref(n,n);