diff --git a/unsupported/Eigen/src/NonLinear/lmder.h b/unsupported/Eigen/src/NonLinear/lmder.h
index d08530947..6090bd48a 100644
--- a/unsupported/Eigen/src/NonLinear/lmder.h
+++ b/unsupported/Eigen/src/NonLinear/lmder.h
@@ -121,7 +121,7 @@ L60:
     /*        on the first iteration, calculate the norm of the scaled x */
     /*        and initialize the step bound delta. */
 
-    wa3 = diag.cwise() * x ;
+    wa3 = diag.cwise() * x;
     xnorm = wa3.stableNorm();
     delta = factor * xnorm;
     if (delta == 0.) {
@@ -161,16 +161,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) * (qtf[i] / fnorm);
-        /* Computing MAX */
-        gnorm = std::max(gnorm, ei_abs(sum / wa2[l]));
-L150:
-        ;
     }
 L170:
 
diff --git a/unsupported/Eigen/src/NonLinear/lmdif.h b/unsupported/Eigen/src/NonLinear/lmdif.h
index 59bc51cbd..42aefb32d 100644
--- a/unsupported/Eigen/src/NonLinear/lmdif.h
+++ b/unsupported/Eigen/src/NonLinear/lmdif.h
@@ -121,7 +121,7 @@ L60:
     /*        and initialize the step bound delta. */
 
     wa3 = diag.cwise() * x;
-    xnorm = wa3.stableNorm();;
+    xnorm = wa3.stableNorm();
     delta = factor * xnorm;
     if (delta == 0.) {
         delta = factor;
diff --git a/unsupported/Eigen/src/NonLinear/lmstr.h b/unsupported/Eigen/src/NonLinear/lmstr.h
index b228a2e9e..dab5f8c4d 100644
--- a/unsupported/Eigen/src/NonLinear/lmstr.h
+++ b/unsupported/Eigen/src/NonLinear/lmstr.h
@@ -7,6 +7,7 @@ int ei_lmstr(
         int &njev,
         Matrix< Scalar, Dynamic, Dynamic > &fjac,
         VectorXi &ipvt,
+        Matrix< Scalar, Dynamic, 1 >  &qtf,
         Matrix< Scalar, Dynamic, 1 >  &diag,
         int mode=1,
         Scalar factor = 100.,
@@ -18,25 +19,23 @@ int ei_lmstr(
         )
 {
     const int m = fvec.size(), n = x.size();
-    Matrix< Scalar, Dynamic, 1 >
-        qtf(n),
-        wa1(n), wa2(n), wa3(n),
-        wa4(m);
+    Matrix< Scalar, Dynamic, 1 > wa1(n), wa2(n), wa3(n), wa4(m);
 
     ipvt.resize(n);
     fjac.resize(m, n);
     diag.resize(n);
+    qtf.resize(n);
 
     /* Local variables */
     int i, j, l;
     Scalar par, sum;
-    int sing;
-    int iter;
+    int sing, iter;
     Scalar temp, temp1, temp2;
     int iflag;
     Scalar delta;
     Scalar ratio;
-    Scalar fnorm, gnorm, pnorm, xnorm, fnorm1, actred, dirder, prered;
+    Scalar fnorm, gnorm;
+    Scalar pnorm, xnorm, fnorm1, actred, dirder, prered;
     int info;
 
     /* Function Body */
@@ -51,7 +50,6 @@ int ei_lmstr(
             gtol < 0. || maxfev <= 0 || factor <= 0.) {
         goto L340;
     }
-
     if (mode == 2)
         for (j = 0; j < n; ++j)
             if (diag[j] <= 0.) goto L300;
@@ -64,7 +62,7 @@ int ei_lmstr(
     if (iflag < 0) {
         goto L340;
     }
-    fnorm = fvec.stableNorm();;
+    fnorm = fvec.stableNorm();
 
     /*     initialize levenberg-marquardt parameter and iteration counter. */
 
@@ -176,19 +174,13 @@ L170:
     }
     for (j = 0; j < n; ++j) {
         l = ipvt[j];
-        if (wa2[l] == 0.) {
-            goto L190;
+        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) * (qtf[i] / fnorm);
-            /* L180: */
-        }
-        /* Computing MAX */
-        gnorm = std::max(gnorm, ei_abs(sum/wa2[l]));
-L190:
-        /* L200: */
-        ;
     }
 L210:
 
@@ -249,8 +241,8 @@ L240:
     /*           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];
         temp = wa1[l];
         for (i = 0; i <= j; ++i) {
@@ -260,7 +252,7 @@ L240:
         /* L270: */
     }
     temp1 = ei_abs2(wa3.stableNorm() / fnorm);
-    temp2 = ei_abs2( ei_sqrt(par) * pnorm / fnorm);
+    temp2 = ei_abs2(ei_sqrt(par) * pnorm / fnorm);
     /* Computing 2nd power */
     prered = temp1 + temp2 / Scalar(.5);
     dirder = -(temp1 + temp2);
@@ -284,9 +276,8 @@ L240:
     if (actred < 0.) {
         temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
     }
-    if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) {
+    if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1))
         temp = Scalar(.1);
-    }
     /* Computing MIN */
     delta = temp * std::min(delta, pnorm / Scalar(.1));
 
diff --git a/unsupported/Eigen/src/NonLinear/lmstr1.h b/unsupported/Eigen/src/NonLinear/lmstr1.h
index 896bba394..737b569ad 100644
--- a/unsupported/Eigen/src/NonLinear/lmstr1.h
+++ b/unsupported/Eigen/src/NonLinear/lmstr1.h
@@ -10,7 +10,7 @@ int ei_lmstr1(
     const int n = x.size(), m=fvec.size();
     int info, nfev=0, njev=0;
     Matrix< Scalar, Dynamic, Dynamic > fjac(m, n);
-    Matrix< Scalar, Dynamic, 1> diag;
+    Matrix< Scalar, Dynamic, 1> diag, qtf;
 
     /* check the input parameters for errors. */
     if (n <= 0 || m < n || tol < 0.) {
@@ -22,7 +22,7 @@ int ei_lmstr1(
     info = ei_lmstr<Functor,Scalar>(
         x, fvec,
         nfev, njev,
-        fjac, ipvt, diag,
+        fjac, ipvt, qtf, diag,
         1,
         100.,
         (n+1)*100,
diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp
index 241714e3f..40d8bd89b 100644
--- a/unsupported/test/NonLinear.cpp
+++ b/unsupported/test/NonLinear.cpp
@@ -461,7 +461,7 @@ void testLmstr()
   const int m=15, n=3;
   int info, nfev=0, njev=0;
   double fnorm;
-  VectorXd x(n), fvec(m), diag(n);
+  VectorXd x(n), fvec(m), diag(n), qtf;
   MatrixXd fjac;
   VectorXi ipvt;
 
@@ -469,7 +469,7 @@ void testLmstr()
   x.setConstant(n, 1.);
 
   // do the computation
-  info = ei_lmstr<lmstr_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
+  info = ei_lmstr<lmstr_functor, double>(x, fvec, nfev, njev, fjac, ipvt, qtf, diag);
   VectorXd wa(n);
 
   // check return values