diff --git a/unsupported/Eigen/src/NonLinear/hybrd.h b/unsupported/Eigen/src/NonLinear/hybrd.h
index 774088516..e935b3c35 100644
--- a/unsupported/Eigen/src/NonLinear/hybrd.h
+++ b/unsupported/Eigen/src/NonLinear/hybrd.h
@@ -61,16 +61,9 @@ int ei_hybrd(
             factor <= 0. || ldfjac < n || lr < n * (n + 1) / 2) {
         goto L300;
     }
-    if (mode != 2) {
-        goto L20;
-    }
-    for (j = 0; j < n; ++j) {
-        if (diag[j] <= 0.) {
-            goto L300;
-        }
-        /* L10: */
-    }
-L20:
+    if (mode == 2)
+        for (j = 0; j < n; ++j)
+            if (diag[j] <= 0.) goto L300;
 
     /*     evaluate the function at the starting point */
     /*     and calculate its norm. */
@@ -135,10 +128,7 @@ L50:
     /*        on the first iteration, calculate the norm of the scaled x */
     /*        and initialize the step bound delta. */
 
-    for (j = 0; j < n; ++j) {
-        wa3[j] = diag[j] * x[j];
-        /* L60: */
-    }
+    wa3 = diag.cwise() * x;
     xnorm = wa3.stableNorm();
     delta = factor * xnorm;
     if (delta == 0.) {
@@ -148,10 +138,7 @@ L70:
 
     /*        form (q transpose)*fvec and store in qtf. */
 
-    for (i = 0; i < n; ++i) {
-        qtf[i] = fvec[i];
-        /* L80: */
-    }
+    qtf = fvec;
     for (j = 0; j < n; ++j) {
         if (fjac(j,j) == 0.) {
             goto L110;
@@ -200,8 +187,7 @@ L110:
         goto L170;
     }
     /* Computing MAX */
-    for (j = 0; j < n; ++j)
-        diag[j] = std::max(diag[j], wa2[j]);
+    diag = diag.cwise().max(wa2);
 L170:
 
     /*        beginning of the inner loop. */
@@ -228,12 +214,9 @@ L190:
 
     /*           store the direction p and x + p. calculate the norm of p. */
 
-    for (j = 0; j < n; ++j) {
-        wa1[j] = -wa1[j];
-        wa2[j] = x[j] + wa1[j];
-        wa3[j] = diag[j] * wa1[j];
-        /* L200: */
-    }
+    wa1 = -wa1;
+    wa2 = x + wa1;
+    wa3 = diag.cwise() * wa1;
     pnorm = wa3.stableNorm();
 
     /*           on the first iteration, adjust the initial step bound. */
@@ -310,12 +293,9 @@ L240:
 
     /*           successful iteration. update x, fvec, and their norms. */
 
-    for (j = 0; j < n; ++j) {
-        x[j] = wa2[j];
-        wa2[j] = diag[j] * x[j];
-        fvec[j] = wa4[j];
-        /* L250: */
-    }
+    x = wa2;
+    wa2 = diag.cwise() * x;
+    fvec = wa4;
     temp = wa2.stableNorm();
     fnorm = fnorm1;
     ++iter;
@@ -368,17 +348,11 @@ L260:
     /*           and update qtf if necessary. */
 
     for (j = 0; j < n; ++j) {
-        sum = 0.;
-        for (i = 0; i < n; ++i) {
-            sum += fjac(i,j) * wa4[i];
-            /* L270: */
-        }
+        sum = wa4.dot(fjac.col(j));
         wa2[j] = (sum - wa3[j]) / pnorm;
         wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
-        if (ratio >= Scalar(1e-4)) {
+        if (ratio >= Scalar(1e-4))
             qtf[j] = sum;
-        }
-        /* L280: */
     }
 
     /*           compute the qr factorization of the updated jacobian. */
diff --git a/unsupported/Eigen/src/NonLinear/hybrj.h b/unsupported/Eigen/src/NonLinear/hybrj.h
index 1b5fec348..2e46d09e5 100644
--- a/unsupported/Eigen/src/NonLinear/hybrj.h
+++ b/unsupported/Eigen/src/NonLinear/hybrj.h
@@ -56,16 +56,10 @@ int ei_hybrj(
             0. || lr < n * (n + 1) / 2) {
         goto L300;
     }
-    if (mode != 2) {
-        goto L20;
-    }
-    for (j = 0; j < n; ++j) {
-        if (diag[j] <= 0.) {
-            goto L300;
-        }
-        /* L10: */
-    }
-L20:
+    if (mode == 2)
+        for (j = 0; j < n; ++j)
+            if (diag[j] <= 0.)
+                goto L300;
 
     /*     evaluate the function at the starting point */
     /*     and calculate its norm. */
@@ -123,10 +117,7 @@ L50:
     /*        on the first iteration, calculate the norm of the scaled x */
     /*        and initialize the step bound delta. */
 
-    for (j = 0; j < n; ++j) {
-        wa3[j] = diag[j] * x[j];
-        /* L60: */
-    }
+    wa3 = diag.cwise() * x;
     xnorm = wa3.stableNorm();;
     delta = factor * xnorm;
     if (delta == 0.) {
@@ -136,10 +127,7 @@ L70:
 
     /*        form (q transpose)*fvec and store in qtf. */
 
-    for (i = 0; i < n; ++i) {
-        qtf[i] = fvec[i];
-        /* L80: */
-    }
+    qtf = fvec;
     for (j = 0; j < n; ++j) {
         if (fjac(j,j) == 0.) {
             goto L110;
@@ -219,7 +207,6 @@ L190:
         wa1[j] = -wa1[j];
         wa2[j] = x[j] + wa1[j];
         wa3[j] = diag[j] * wa1[j];
-        /* L200: */
     }
     pnorm = wa3.stableNorm();
 
@@ -297,12 +284,9 @@ L240:
 
     /*           successful iteration. update x, fvec, and their norms. */
 
-    for (j = 0; j < n; ++j) {
-        x[j] = wa2[j];
-        wa2[j] = diag[j] * x[j];
-        fvec[j] = wa4[j];
-        /* L250: */
-    }
+    x =wa2;
+    wa2 = diag.cwise() * x;
+    fvec = wa4;
     xnorm = wa2.stableNorm();
     fnorm = fnorm1;
     ++iter;
@@ -359,11 +343,7 @@ L260:
     /*           and update qtf if necessary. */
 
     for (j = 0; j < n; ++j) {
-        sum = 0.;
-        for (i = 0; i < n; ++i) {
-            sum += fjac(i,j) * wa4[i];
-            /* L270: */
-        }
+        sum = wa4.dot(fjac.col(j));
         wa2[j] = (sum - wa3[j]) / pnorm;
         wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
         if (ratio >= Scalar(1e-4)) {
diff --git a/unsupported/Eigen/src/NonLinear/lmder.h b/unsupported/Eigen/src/NonLinear/lmder.h
index a4422258e..8f0db3e0a 100644
--- a/unsupported/Eigen/src/NonLinear/lmder.h
+++ b/unsupported/Eigen/src/NonLinear/lmder.h
@@ -50,16 +50,10 @@ int ei_lmder(
             gtol < 0. || maxfev <= 0 || factor <= 0.) {
         goto L300;
     }
-    if (mode != 2) {
-        goto L20;
-    }
-    for (j = 0; j < n; ++j) {
-        if (diag[j] <= 0.) {
-            goto L300;
-        }
-        /* L10: */
-    }
-L20:
+    if (mode == 2)
+        for (j = 0; j < n; ++j)
+            if (diag[j] <= 0.)
+                goto L300;
 
     /*     evaluate the function at the starting point */
     /*     and calculate its norm. */
@@ -128,10 +122,7 @@ L60:
     /*        on the first iteration, calculate the norm of the scaled x */
     /*        and initialize the step bound delta. */
 
-    for (j = 0; j < n; ++j) {
-        wa3[j] = diag[j] * x[j];
-        /* L70: */
-    }
+    wa3 = diag.cwise() * x ;
     xnorm = wa3.stableNorm();
     delta = factor * xnorm;
     if (delta == 0.) {
@@ -142,10 +133,7 @@ L80:
     /*        form (q transpose)*fvec and store the first n components in */
     /*        qtf. */
 
-    for (i = 0; i < m; ++i) {
-        wa4[i] = fvec[i];
-        /* L90: */
-    }
+    wa4 = fvec;
     for (j = 0; j < n; ++j) {
         if (fjac(j,j) == 0.) {
             goto L120;
@@ -178,14 +166,11 @@ L120:
             goto L150;
         }
         sum = 0.;
-        for (i = 0; i <= j; ++i) {
+        for (i = 0; i <= j; ++i)
             sum += fjac(i,j) * (qtf[i] / fnorm);
-            /* L140: */
-        }
         /* Computing MAX */
         gnorm = std::max(gnorm, ei_abs(sum / wa2[l]));
 L150:
-        /* L160: */
         ;
     }
 L170:
@@ -221,12 +206,9 @@ L200:
 
     /*           store the direction p and x + p. calculate the norm of p. */
 
-    for (j = 0; j < n; ++j) {
-        wa1[j] = -wa1[j];
-        wa2[j] = x[j] + wa1[j];
-        wa3[j] = diag[j] * wa1[j];
-        /* L210: */
-    }
+    wa1 = -wa1;
+    wa2 = x + wa1;
+    wa3 = diag.cwise() * wa1;
     pnorm = wa3.stableNorm();
 
     /*           on the first iteration, adjust the initial step bound. */
@@ -253,6 +235,7 @@ L200:
     /*           compute the scaled predicted reduction and */
     /*           the scaled directional derivative. */
 
+    wa3.fill(0.);
     for (j = 0; j < n; ++j) {
         wa3[j] = 0.;
         l = ipvt[j];
@@ -312,15 +295,9 @@ L260:
 
     /*           successful iteration. update x, fvec, and their norms. */
 
-    for (j = 0; j < n; ++j) {
-        x[j] = wa2[j];
-        wa2[j] = diag[j] * x[j];
-        /* L270: */
-    }
-    for (i = 0; i < m; ++i) {
-        fvec[i] = wa4[i];
-        /* L280: */
-    }
+    x = wa2;
+    wa2 = diag.cwise() * x;
+    fvec = wa4;
     xnorm = wa2.stableNorm();
     fnorm = fnorm1;
     ++iter;
diff --git a/unsupported/Eigen/src/NonLinear/lmdif.h b/unsupported/Eigen/src/NonLinear/lmdif.h
index 267b1d6e3..a9883228e 100644
--- a/unsupported/Eigen/src/NonLinear/lmdif.h
+++ b/unsupported/Eigen/src/NonLinear/lmdif.h
@@ -52,16 +52,10 @@ int ei_lmdif(
             gtol < 0. || maxfev <= 0 || factor <= 0.) {
         goto L300;
     }
-    if (mode != 2) {
-        goto L20;
-    }
-    for (j = 0; j < n; ++j) {
-        if (diag[j] <= 0.) {
-            goto L300;
-        }
-        /* L10: */
-    }
-L20:
+    if (mode == 2)
+        for (j = 0; j < n; ++j)
+            if (diag[j] <= 0.)
+                goto L300;
 
     /*     evaluate the function at the starting point */
     /*     and calculate its norm. */
@@ -130,10 +124,7 @@ L60:
     /*        on the first iteration, calculate the norm of the scaled x */
     /*        and initialize the step bound delta. */
 
-    for (j = 0; j < n; ++j) {
-        wa3[j] = diag[j] * x[j];
-        /* L70: */
-    }
+    wa3 = diag.cwise() * x;
     xnorm = wa3.stableNorm();;
     delta = factor * xnorm;
     if (delta == 0.) {
@@ -144,10 +135,7 @@ L80:
     /*        form (q transpose)*fvec and store the first n components in */
     /*        qtf. */
 
-    for (i = 0; i < m; ++i) {
-        wa4[i] = fvec[i];
-        /* L90: */
-    }
+    wa4 = fvec;
     for (j = 0; j < n; ++j) {
         if (fjac[j + j * ldfjac] == 0.) {
             goto L120;
@@ -223,12 +211,9 @@ L200:
 
     /*           store the direction p and x + p. calculate the norm of p. */
 
-    for (j = 0; j < n; ++j) {
-        wa1[j] = -wa1[j];
-        wa2[j] = x[j] + wa1[j];
-        wa3[j] = diag[j] * wa1[j];
-        /* L210: */
-    }
+    wa1 = -wa1;
+    wa2 = x + wa1;
+    wa3 = diag.cwise() * wa1;
     pnorm = wa3.stableNorm();
 
     /*           on the first iteration, adjust the initial step bound. */
@@ -314,15 +299,9 @@ L260:
 
     /*           successful iteration. update x, fvec, and their norms. */
 
-    for (j = 0; j < n; ++j) {
-        x[j] = wa2[j];
-        wa2[j] = diag[j] * x[j];
-        /* L270: */
-    }
-    for (i = 0; i < m; ++i) {
-        fvec[i] = wa4[i];
-        /* L280: */
-    }
+    x = wa2;
+    wa2 = diag.cwise() * x;
+    fvec = wa4;
     xnorm = wa2.stableNorm();
     fnorm = fnorm1;
     ++iter;
diff --git a/unsupported/Eigen/src/NonLinear/lmstr.h b/unsupported/Eigen/src/NonLinear/lmstr.h
index 44fcce623..cd50b72a8 100644
--- a/unsupported/Eigen/src/NonLinear/lmstr.h
+++ b/unsupported/Eigen/src/NonLinear/lmstr.h
@@ -52,16 +52,11 @@ int ei_lmstr(
             gtol < 0. || maxfev <= 0 || factor <= 0.) {
         goto L340;
     }
-    if (mode != 2) {
-        goto L20;
-    }
-    for (j = 0; j < n; ++j) {
-        if (diag[j] <= 0.) {
-            goto L340;
-        }
-        /* L10: */
-    }
-L20:
+
+    if (mode == 2)
+        for (j = 0; j < n; ++j)
+            if (diag[j] <= 0.)
+                goto L300;
 
     /*     evaluate the function at the starting point */
     /*     and calculate its norm. */
@@ -248,12 +243,9 @@ L240:
 
     /*           store the direction p and x + p. calculate the norm of p. */
 
-    for (j = 0; j < n; ++j) {
-        wa1[j] = -wa1[j];
-        wa2[j] = x[j] + wa1[j];
-        wa3[j] = diag[j] * wa1[j];
-        /* L250: */
-    }
+    wa1 = -wa1;
+    wa2 = x + wa1;
+    wa3 = diag.cwise() * wa1;
     pnorm = wa3.stableNorm();
 
     /*           on the first iteration, adjust the initial step bound. */
@@ -340,15 +332,9 @@ L300:
 
     /*           successful iteration. update x, fvec, and their norms. */
 
-    for (j = 0; j < n; ++j) {
-        x[j] = wa2[j];
-        wa2[j] = diag[j] * x[j];
-        /* L310: */
-    }
-    for (i = 0; i < m; ++i) {
-        fvec[i] = wa4[i];
-        /* L320: */
-    }
+    x = wa2;
+    wa2 = diag.cwise() * x;
+    fvec = wa4;
     xnorm = wa2.stableNorm();
     fnorm = fnorm1;
     ++iter;