diff --git a/unsupported/Eigen/src/NonLinear/chkder.h b/unsupported/Eigen/src/NonLinear/chkder.h
index 913727227..692b95f31 100644
--- a/unsupported/Eigen/src/NonLinear/chkder.h
+++ b/unsupported/Eigen/src/NonLinear/chkder.h
@@ -20,10 +20,9 @@ void ei_chkder(
     int i,j;
 
     const int m = fvec.size(), n = x.size();
-    int ldfjac = m;
 
     if (mode != 2) {
-        xp.resize(ldfjac);
+        xp.resize(m);
         /*        mode = 1. */
         for (j = 0; j < n; ++j) {
             temp = eps * ei_abs(x[j]);
@@ -34,7 +33,7 @@ void ei_chkder(
     }
     else {
         /*        mode = 2. */
-        err.setZero(ldfjac); 
+        err.setZero(m); 
         for (j = 0; j < n; ++j) {
             temp = ei_abs(x[j]);
             if (temp == 0.)
@@ -43,20 +42,13 @@ void ei_chkder(
         }
         for (i = 0; i < m; ++i) {
             temp = 1.;
-            if (fvec[i] != 0. && fvecp[i] != 0. && ei_abs(fvecp[i] - 
-                        fvec[i]) >= epsf * ei_abs(fvec[i]))
-            {
-                temp = eps * ei_abs((fvecp[i] - fvec[i]) / eps - err[i]) 
-                    / (ei_abs(fvec[i]) +
-                            ei_abs(fvecp[i]));
-            }
+            if (fvec[i] != 0. && fvecp[i] != 0. && ei_abs(fvecp[i] - fvec[i]) >= epsf * ei_abs(fvec[i]))
+                temp = eps * ei_abs((fvecp[i] - fvec[i]) / eps - err[i]) / (ei_abs(fvec[i]) + ei_abs(fvecp[i]));
             err[i] = 1.;
-            if (temp > epsilon<Scalar>() && temp < eps) {
+            if (temp > epsilon<Scalar>() && temp < eps)
                 err[i] = (chkder_log10e * ei_log(temp) - epslog) / epslog;
-            }
-            if (temp >= eps) {
+            if (temp >= eps)
                 err[i] = 0.;
-            }
         }
     }
 } /* chkder_ */
diff --git a/unsupported/Eigen/src/NonLinear/hybrd.h b/unsupported/Eigen/src/NonLinear/hybrd.h
index e935b3c35..3b17e444b 100644
--- a/unsupported/Eigen/src/NonLinear/hybrd.h
+++ b/unsupported/Eigen/src/NonLinear/hybrd.h
@@ -28,8 +28,7 @@ int ei_hybrd(
     fvec.resize(n);
     qtf.resize(n);
     R.resize(lr);
-    int ldfjac = n;
-    fjac.resize(ldfjac, n);
+    fjac.resize(n, n);
 
     /* Local variables */
     int i, j, l, iwa[1];
@@ -58,7 +57,7 @@ int ei_hybrd(
     /*     check the input parameters for errors. */
 
     if (n <= 0 || xtol < 0. || maxfev <= 0 || nb_of_subdiagonals < 0 || nb_of_superdiagonals < 0 ||
-            factor <= 0. || ldfjac < n || lr < n * (n + 1) / 2) {
+            factor <= 0. || lr < n * (n + 1) / 2) {
         goto L300;
     }
     if (mode == 2)
@@ -105,7 +104,7 @@ L30:
 
     /*        compute the qr factorization of the jacobian. */
 
-    ei_qrfac<Scalar>(n, n, fjac.data(), ldfjac, false, iwa, 1, wa1.data(), wa2.data(), wa3.data());
+    ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, 1, wa1.data(), wa2.data(), wa3.data());
 
     /*        on the first iteration and if mode is 1, scale according */
     /*        to the norms of the columns of the initial jacobian. */
@@ -179,7 +178,7 @@ L110:
 
     /*        accumulate the orthogonal factor in fjac. */
 
-    ei_qform<Scalar>(n, n, fjac.data(), ldfjac, wa1.data());
+    ei_qform<Scalar>(n, n, fjac.data(), fjac.rows(), wa1.data());
 
     /*        rescale if necessary. */
 
@@ -358,7 +357,7 @@ L260:
     /*           compute the qr factorization of the updated jacobian. */
 
     ei_r1updt<Scalar>(n, n, R.data(), lr, wa1.data(), wa2.data(), wa3.data(), &sing);
-    ei_r1mpyq<Scalar>(n, n, fjac.data(), ldfjac, wa2.data(), wa3.data());
+    ei_r1mpyq<Scalar>(n, n, fjac.data(), fjac.rows(), wa2.data(), wa3.data());
     ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data());
 
     /*           end of the inner loop. */
diff --git a/unsupported/Eigen/src/NonLinear/hybrj.h b/unsupported/Eigen/src/NonLinear/hybrj.h
index 2e46d09e5..78c8a00a8 100644
--- a/unsupported/Eigen/src/NonLinear/hybrj.h
+++ b/unsupported/Eigen/src/NonLinear/hybrj.h
@@ -23,8 +23,7 @@ int ei_hybrj(
     fvec.resize(n);
     qtf.resize(n);
     R.resize(lr);
-    int ldfjac = n;
-    fjac.resize(ldfjac, n);
+    fjac.resize(n, n);
 
     /* Local variables */
     int i, j, l, iwa[1];
@@ -52,7 +51,7 @@ int ei_hybrj(
 
     /*     check the input parameters for errors. */
 
-    if (n <= 0 || ldfjac < n || xtol < 0. || maxfev <= 0 || factor <= 
+    if (n <= 0 || xtol < 0. || maxfev <= 0 || factor <= 
             0. || lr < n * (n + 1) / 2) {
         goto L300;
     }
@@ -94,7 +93,7 @@ L30:
 
     /*        compute the qr factorization of the jacobian. */
 
-    ei_qrfac<Scalar>(n, n,fjac.data(), ldfjac, false, iwa, 1, wa1.data(), wa2.data(), wa3.data());
+    ei_qrfac<Scalar>(n, n,fjac.data(), fjac.rows(), false, iwa, 1, wa1.data(), wa2.data(), wa3.data());
 
     /*        on the first iteration and if mode is 1, scale according */
     /*        to the norms of the columns of the initial jacobian. */
@@ -168,7 +167,7 @@ L110:
 
     /*        accumulate the orthogonal factor in fjac. */
 
-    ei_qform<Scalar>(n, n, fjac.data(), ldfjac, wa1.data());
+    ei_qform<Scalar>(n, n, fjac.data(), fjac.rows(), wa1.data());
 
     /*        rescale if necessary. */
 
@@ -355,7 +354,7 @@ L260:
     /*           compute the qr factorization of the updated jacobian. */
 
     ei_r1updt<Scalar>(n, n, R.data(), lr, wa1.data(), wa2.data(), wa3.data(), &sing);
-    ei_r1mpyq<Scalar>(n, n, fjac.data(), ldfjac, wa2.data(), wa3.data());
+    ei_r1mpyq<Scalar>(n, n, fjac.data(), fjac.rows(), wa2.data(), wa3.data());
     ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data());
 
     /*           end of the inner loop. */
diff --git a/unsupported/Eigen/src/NonLinear/lmder.h b/unsupported/Eigen/src/NonLinear/lmder.h
index 8f0db3e0a..57407911b 100644
--- a/unsupported/Eigen/src/NonLinear/lmder.h
+++ b/unsupported/Eigen/src/NonLinear/lmder.h
@@ -19,10 +19,9 @@ int ei_lmder(
 {
     const int m = fvec.size(), n = x.size();
     Matrix< Scalar, Dynamic, 1 > qtf(n), wa1(n), wa2(n), wa3(n), wa4(m);
-    int ldfjac = m;
 
     ipvt.resize(n);
-    fjac.resize(ldfjac, n);
+    fjac.resize(m, n);
     diag.resize(n);
 
     /* Local variables */
@@ -46,7 +45,7 @@ int ei_lmder(
 
     /*     check the input parameters for errors. */
 
-    if (n <= 0 || m < n || ldfjac < m || ftol < 0. || xtol < 0. || 
+    if (n <= 0 || m < n || ftol < 0. || xtol < 0. || 
             gtol < 0. || maxfev <= 0 || factor <= 0.) {
         goto L300;
     }
@@ -98,7 +97,7 @@ L40:
 
     /*        compute the qr factorization of the jacobian. */
 
-    ei_qrfac<Scalar>(m, n, fjac.data(), ldfjac, true, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data());
+    ei_qrfac<Scalar>(m, n, fjac.data(), fjac.rows(), 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 */
@@ -200,7 +199,7 @@ L200:
     /*           determine the levenberg-marquardt parameter. */
 
     ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
-    ei_lmpar<Scalar>(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta,
+    ei_lmpar<Scalar>(n, fjac.data(), fjac.rows(), ipvt.data(), diag.data(), qtf.data(), delta,
             par, wa1.data(), wa2.data(), wa3.data(), wa4.data());
     ipvt.cwise()-=1;
 
diff --git a/unsupported/Eigen/src/NonLinear/lmdif.h b/unsupported/Eigen/src/NonLinear/lmdif.h
index a9883228e..ab9f9ad04 100644
--- a/unsupported/Eigen/src/NonLinear/lmdif.h
+++ b/unsupported/Eigen/src/NonLinear/lmdif.h
@@ -22,10 +22,9 @@ int ei_lmdif(
     Matrix< Scalar, Dynamic, 1 >
         wa1(n), wa2(n), wa3(n),
         wa4(m);
-    int ldfjac = m;
 
     ipvt.resize(n);
-    fjac.resize(ldfjac, n);
+    fjac.resize(m, n);
     diag.resize(n);
     qtf.resize(n);
 
@@ -48,7 +47,7 @@ int ei_lmdif(
 
     /*     check the input parameters for errors. */
 
-    if (n <= 0 || m < n || ldfjac < m || ftol < 0. || xtol < 0. || 
+    if (n <= 0 || m < n || ftol < 0. || xtol < 0. || 
             gtol < 0. || maxfev <= 0 || factor <= 0.) {
         goto L300;
     }
@@ -100,7 +99,7 @@ L40:
 
     /*        compute the qr factorization of the jacobian. */
 
-    ei_qrfac<Scalar>(m, n, fjac.data(), ldfjac, true, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data());
+    ei_qrfac<Scalar>(m, n, fjac.data(), fjac.rows(), 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 */
@@ -137,21 +136,21 @@ L80:
 
     wa4 = fvec;
     for (j = 0; j < n; ++j) {
-        if (fjac[j + j * ldfjac] == 0.) {
+        if (fjac(j,j) == 0.) {
             goto L120;
         }
         sum = 0.;
         for (i = j; i < m; ++i) {
-            sum += fjac[i + j * ldfjac] * wa4[i];
+            sum += fjac(i,j) * wa4[i];
             /* L100: */
         }
-        temp = -sum / fjac[j + j * ldfjac];
+        temp = -sum / fjac(j,j);
         for (i = j; i < m; ++i) {
-            wa4[i] += fjac[i + j * ldfjac] * temp;
+            wa4[i] += fjac(i,j) * temp;
             /* L110: */
         }
 L120:
-        fjac[j + j * ldfjac] = wa1[j];
+        fjac(j,j) = wa1[j];
         qtf[j] = wa4[j];
         /* L130: */
     }
@@ -164,19 +163,13 @@ L120:
     }
     for (j = 0; j < n; ++j) {
         l = ipvt[j];
-        if (wa2[l] == 0.) {
-            goto L150;
+        if (wa2[l] != 0.) {
+            sum = 0.;
+            for (i = 0; i <= j; ++i)
+                sum += fjac(i,j) * (qtf[i] / fnorm);
+            /* Computing MAX */
+            gnorm = std::max(gnorm, ei_abs(sum / wa2[l]));
         }
-        sum = 0.;
-        for (i = 0; i <= j; ++i) {
-            sum += fjac[i + j * ldfjac] * (qtf[i] / fnorm);
-            /* L140: */
-        }
-        /* Computing MAX */
-        gnorm = std::max(gnorm, ei_abs(sum / wa2[l]));
-L150:
-        /* L160: */
-        ;
     }
 L170:
 
@@ -205,7 +198,7 @@ L200:
     /*           determine the levenberg-marquardt parameter. */
 
     ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
-    ei_lmpar<Scalar>(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta,
+    ei_lmpar<Scalar>(n, fjac.data(), fjac.rows(), ipvt.data(), diag.data(), qtf.data(), delta,
             par, wa1.data(), wa2.data(), wa3.data(), wa4.data());
     ipvt.cwise()-=1;
 
@@ -245,7 +238,7 @@ L200:
         l = ipvt[j];
         temp = wa1[l];
         for (i = 0; i <= j; ++i) {
-            wa3[i] += fjac[i + j * ldfjac] * temp;
+            wa3[i] += fjac(i,j) * temp;
             /* L220: */
         }
         /* L230: */
diff --git a/unsupported/Eigen/src/NonLinear/lmstr.h b/unsupported/Eigen/src/NonLinear/lmstr.h
index cd50b72a8..f8257feb3 100644
--- a/unsupported/Eigen/src/NonLinear/lmstr.h
+++ b/unsupported/Eigen/src/NonLinear/lmstr.h
@@ -22,10 +22,9 @@ int ei_lmstr(
         qtf(n),
         wa1(n), wa2(n), wa3(n),
         wa4(m);
-    int ldfjac = m;
 
     ipvt.resize(n);
-    fjac.resize(ldfjac, n);
+    fjac.resize(m, n);
     diag.resize(n);
 
     /* Local variables */
@@ -48,7 +47,7 @@ int ei_lmstr(
 
     /*     check the input parameters for errors. */
 
-    if (n <= 0 || m < n || ldfjac < n || ftol < 0. || xtol < 0. || 
+    if (n <= 0 || m < n || ftol < 0. || xtol < 0. || 
             gtol < 0. || maxfev <= 0 || factor <= 0.) {
         goto L340;
     }
@@ -110,7 +109,7 @@ L40:
             goto L340;
         }
         temp = fvec[i];
-        ei_rwupdt<Scalar>(n, fjac.data(), ldfjac, wa3.data(), qtf.data(), &temp, wa1.data(), wa2.data());
+        ei_rwupdt<Scalar>(n, fjac.data(), fjac.rows(), wa3.data(), qtf.data(), &temp, wa1.data(), wa2.data());
         ++iflag;
         /* L70: */
     }
@@ -126,15 +125,12 @@ L40:
         }
         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);
-        /* L80: */
     }
     if (! sing) {
         goto L130;
     }
     ipvt.cwise()+=1;
-    ei_qrfac<Scalar>(n, n, fjac.data(), ldfjac, true, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data());
+    ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), 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.) {
@@ -237,7 +233,7 @@ L240:
     /*           determine the levenberg-marquardt parameter. */
 
     ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
-    ei_lmpar<Scalar>(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta, par,
+    ei_lmpar<Scalar>(n, fjac.data(), fjac.rows(), ipvt.data(), diag.data(), qtf.data(), delta, par,
             wa1.data(), wa2.data(), wa3.data(), wa4.data());
     ipvt.cwise()-=1;