diff --git a/unsupported/Eigen/src/NonLinear/fdjac1.h b/unsupported/Eigen/src/NonLinear/fdjac1.h
index 2fd0696f7..9e9baf4a3 100644
--- a/unsupported/Eigen/src/NonLinear/fdjac1.h
+++ b/unsupported/Eigen/src/NonLinear/fdjac1.h
@@ -34,7 +34,7 @@ int ei_fdjac1(
 	if (h == 0.)
 	    h = eps;
 	x[j] = temp + h;
-	iflag = Functor::f(x, wa1, 1);
+	iflag = Functor::f(x, wa1);
 	if (iflag < 0)
 	    goto L30;
 	x[j] = temp;
@@ -59,7 +59,7 @@ L40:
 	    x[j] = wa2[j] + h;
 /* L60: */
 	}
-	iflag = Functor::f(x, wa1, 1);
+	iflag = Functor::f(x, wa1);
 	if (iflag < 0) {
 	    /* goto L100; */
             return iflag;
diff --git a/unsupported/Eigen/src/NonLinear/fdjac2.h b/unsupported/Eigen/src/NonLinear/fdjac2.h
index 2a4b9734b..fcd986a95 100644
--- a/unsupported/Eigen/src/NonLinear/fdjac2.h
+++ b/unsupported/Eigen/src/NonLinear/fdjac2.h
@@ -26,7 +26,7 @@ int ei_fdjac2(
 	    h = eps;
 	}
 	x[j] = temp + h;
-	iflag = Functor::f(x, wa, 1);
+	iflag = Functor::f(x, wa);
 	if (iflag < 0) {
 	    /* goto L30; */
             return iflag;
diff --git a/unsupported/Eigen/src/NonLinear/hybrd.h b/unsupported/Eigen/src/NonLinear/hybrd.h
index 772a57b67..774088516 100644
--- a/unsupported/Eigen/src/NonLinear/hybrd.h
+++ b/unsupported/Eigen/src/NonLinear/hybrd.h
@@ -75,7 +75,7 @@ L20:
     /*     evaluate the function at the starting point */
     /*     and calculate its norm. */
 
-    iflag = Functor::f(x, fvec, 1);
+    iflag = Functor::f(x, fvec);
     nfev = 1;
     if (iflag < 0) {
         goto L300;
@@ -215,7 +215,7 @@ L180:
     }
     iflag = 0;
     if ((iter - 1) % nprint == 0) {
-        iflag = Functor::f(x, fvec, 0);
+        iflag = Functor::debug(x, fvec);
     }
     if (iflag < 0) {
         goto L300;
@@ -244,7 +244,7 @@ L190:
 
     /*           evaluate the function at x + p and calculate its norm. */
 
-    iflag = Functor::f(wa2, wa4, 0);
+    iflag = Functor::f(wa2, wa4);
     ++nfev;
     if (iflag < 0) {
         goto L300;
@@ -404,7 +404,7 @@ L300:
         info = iflag;
     }
     if (nprint > 0) {
-        iflag = Functor::f(x, fvec, 0);
+        iflag = Functor::debug(x, fvec);
     }
     return info;
 
diff --git a/unsupported/Eigen/src/NonLinear/hybrj.h b/unsupported/Eigen/src/NonLinear/hybrj.h
index 5533196bc..1b5fec348 100644
--- a/unsupported/Eigen/src/NonLinear/hybrj.h
+++ b/unsupported/Eigen/src/NonLinear/hybrj.h
@@ -70,7 +70,7 @@ L20:
     /*     evaluate the function at the starting point */
     /*     and calculate its norm. */
 
-    iflag = Functor::f(x, fvec, fjac, 1);
+    iflag = Functor::f(x, fvec);
     nfev = 1;
     if (iflag < 0) {
         goto L300;
@@ -92,7 +92,7 @@ L30:
 
     /*        calculate the jacobian matrix. */
 
-    iflag = Functor::f(x, fvec, fjac, 2);
+    iflag = Functor::df(x, fjac);
     ++njev;
     if (iflag < 0) {
         goto L300;
@@ -203,7 +203,7 @@ L180:
     }
     iflag = 0;
     if ((iter - 1) % nprint == 0)
-        iflag = Functor::f(x, fvec, fjac, 0);
+        iflag = Functor::debug(x, fvec, fjac);
     if (iflag < 0) {
         goto L300;
     }
@@ -231,7 +231,7 @@ L190:
 
     /*           evaluate the function at x + p and calculate its norm. */
 
-    iflag = Functor::f(wa2, wa4, fjac, 1);
+    iflag = Functor::f(wa2, wa4);
     ++nfev;
     if (iflag < 0) {
         goto L300;
@@ -395,7 +395,7 @@ L300:
         info = iflag;
     }
     if (nprint > 0)
-        iflag = Functor::f(x, fvec, fjac, 0);
+        iflag = Functor::debug(x, fvec, fjac);
     return info;
 
     /*     last card of subroutine hybrj. */
diff --git a/unsupported/Eigen/src/NonLinear/lmdif.h b/unsupported/Eigen/src/NonLinear/lmdif.h
index afcc57fdf..267b1d6e3 100644
--- a/unsupported/Eigen/src/NonLinear/lmdif.h
+++ b/unsupported/Eigen/src/NonLinear/lmdif.h
@@ -66,7 +66,7 @@ L20:
     /*     evaluate the function at the starting point */
     /*     and calculate its norm. */
 
-    iflag = Functor::f(x, fvec, 1);
+    iflag = Functor::f(x, fvec);
     nfev = 1;
     if (iflag < 0) {
         goto L300;
@@ -97,7 +97,7 @@ L30:
     }
     iflag = 0;
     if ((iter - 1) % nprint == 0) {
-        iflag = Functor::f(x, fvec, 0);
+        iflag = Functor::debug(x, fvec);
     }
     if (iflag < 0) {
         goto L300;
@@ -239,7 +239,7 @@ L200:
 
     /*           evaluate the function at x + p and calculate its norm. */
 
-    iflag = Functor::f(wa2, wa4, 1);
+    iflag = Functor::f(wa2, wa4);
     ++nfev;
     if (iflag < 0) {
         goto L300;
@@ -380,7 +380,7 @@ L300:
     }
     iflag = 0;
     if (nprint > 0) {
-        iflag = Functor::f(x, fvec, 0);
+        iflag = Functor::debug(x, fvec);
     }
     return info;
 
diff --git a/unsupported/Eigen/src/NonLinear/lmstr.h b/unsupported/Eigen/src/NonLinear/lmstr.h
index 94d01f32c..44fcce623 100644
--- a/unsupported/Eigen/src/NonLinear/lmstr.h
+++ b/unsupported/Eigen/src/NonLinear/lmstr.h
@@ -66,7 +66,7 @@ L20:
     /*     evaluate the function at the starting point */
     /*     and calculate its norm. */
 
-    iflag = Functor::f(x, fvec, wa3, 1);
+    iflag = Functor::f(x, fvec);
     nfev = 1;
     if (iflag < 0) {
         goto L340;
@@ -89,7 +89,7 @@ L30:
     }
     iflag = 0;
     if ((iter - 1) % nprint == 0) {
-        iflag = Functor::f(x, fvec, wa3, 0);
+        iflag = Functor::debug(x, fvec, wa3);
     }
     if (iflag < 0) {
         goto L340;
@@ -111,7 +111,7 @@ L40:
     }
     iflag = 2;
     for (i = 0; i < m; ++i) {
-        if (Functor::f(x, fvec, wa3, iflag) < 0) {
+        if (Functor::df(x, wa3, iflag) < 0) {
             goto L340;
         }
         temp = fvec[i];
@@ -264,7 +264,7 @@ L240:
 
     /*           evaluate the function at x + p and calculate its norm. */
 
-    iflag = Functor::f(wa2, wa4, wa3, 1);
+    iflag = Functor::f(wa2, wa4);
     ++nfev;
     if (iflag < 0) {
         goto L340;
@@ -406,7 +406,7 @@ L340:
     }
     iflag = 0;
     if (nprint > 0) {
-        iflag = Functor::f(x, fvec, wa3, 0);
+        iflag = Functor::debug(x, fvec, wa3);
     }
     return info;
 
diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp
index 35dddba8b..354b3bf42 100644
--- a/unsupported/test/NonLinear.cpp
+++ b/unsupported/test/NonLinear.cpp
@@ -214,40 +214,36 @@ void testLmder()
 }
 
 struct hybrj_functor {
-    static int f(const VectorXd &x, VectorXd &fvec, MatrixXd &fjac, int iflag)
+
+    static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;}
+    static int f(const VectorXd &x, VectorXd &fvec)
     {
-        /*      subroutine fcn for hybrj1 example. */
-
-        int j, k;
-        double one=1, temp, temp1, temp2, three=3, two=2, zero=0, four=4;
+        double temp, temp1, temp2;
         const int n = x.size();
-
         assert(fvec.size()==n);
+        for (int k = 0; k < n; k++)
+        {
+            temp = (3. - 2.*x[k])*x[k];
+            temp1 = 0.;
+            if (k) temp1 = x[k-1];
+            temp2 = 0.;
+            if (k != n-1) temp2 = x[k+1];
+            fvec[k] = temp - temp1 - 2.*temp2 + 1.;
+        }
+        return 0;
+    }
+    static int df(const VectorXd &x, MatrixXd &fjac)
+    {
+        const int n = x.size();
         assert(fjac.rows()==n);
         assert(fjac.cols()==n);
-
-        if (iflag != 2)
+        for (int k = 0; k < n; k++)
         {
-            for (k = 0; k < n; k++)
-            {
-                temp = (three - two*x[k])*x[k];
-                temp1 = zero;
-                if (k) temp1 = x[k-1];
-                temp2 = zero;
-                if (k != n-1) temp2 = x[k+1];
-                fvec[k] = temp - temp1 - two*temp2 + one;
-            }
-        }
-        else
-        {
-            for (k = 0; k < n; k++)
-            {
-                for (j = 0; j < n; j++)
-                    fjac(k,j) = zero;
-                fjac(k,k) = three - four*x[k];
-                if (k) fjac(k,k-1) = -one;
-                if (k != n-1) fjac(k,k+1) = -two;
-            }
+            for (int j = 0; j < n; j++)
+                fjac(k,j) = 0.;
+            fjac(k,k) = 3.- 4.*x[k];
+            if (k) fjac(k,k-1) = -1.;
+            if (k != n-1) fjac(k,k+1) = -2.;
         }
         return 0;
     }
@@ -319,24 +315,21 @@ void testHybrj()
 }
 
 struct hybrd_functor {
-    static int f(const VectorXd &x, VectorXd &fvec, int /*iflag*/)
+    static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */) { return 0;}
+    static int f(const VectorXd &x, VectorXd &fvec)
     {
-        /*      subroutine fcn for hybrd1 example. */
-
-        int k;
-        double one=1, temp, temp1, temp2, three=3, two=2, zero=0;
+        double temp, temp1, temp2;
         const int n = x.size();
 
         assert(fvec.size()==n);
-
-        for (k=0; k < n; k++)
+        for (int k=0; k < n; k++)
         {
-            temp = (three - two*x[k])*x[k];
-            temp1 = zero;
+            temp = (3. - 2.*x[k])*x[k];
+            temp1 = 0.;
             if (k) temp1 = x[k-1];
-            temp2 = zero;
+            temp2 = 0.;
             if (k != n-1) temp2 = x[k+1];
-            fvec[k] = temp - temp1 - two*temp2 + one;
+            fvec[k] = temp - temp1 - 2.*temp2 + 1.;
         }
         return 0;
     }
@@ -400,41 +393,42 @@ void testHybrd()
 }
 
 struct lmstr_functor {
-    static int f(const VectorXd &x, VectorXd &fvec, VectorXd &fjrow, int iflag)
+    static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const VectorXd & /* fjac */) { return 0;}
+    static int f(const VectorXd &x, VectorXd &fvec)
     {
         /*  subroutine fcn for lmstr1 example. */
-        int i;
-        double tmp1, tmp2, tmp3, tmp4;
+        double tmp1, tmp2, tmp3;
         double y[15]={1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1,
             3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
 
         assert(15==fvec.size());
         assert(3==x.size());
-        assert(fjrow.size()==x.size());
 
-        if (iflag < 2)
+        for (int i=0; i<15; i++)
         {
-            for (i=0; i<15; i++)
-            {
-                tmp1 = i+1;
-                tmp2 = 16 - i - 1;
-                tmp3 = (i>=8)? tmp2 : tmp1;
-                fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
-            }
-        }
-        else
-        {
-            i = iflag-2;
             tmp1 = i+1;
             tmp2 = 16 - i - 1;
             tmp3 = (i>=8)? tmp2 : tmp1;
-            tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4;
-            fjrow[0] = -1;
-            fjrow[1] = tmp1*tmp2/tmp4;
-            fjrow[2] = tmp1*tmp3/tmp4;
+            fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
         }
         return 0;
     }
+    static int df(const VectorXd &x, VectorXd &jac_row, int rownb)
+    {
+        assert(x.size()==3);
+        assert(jac_row.size()==x.size());
+        double tmp1, tmp2, tmp3, tmp4;
+
+        int i = rownb-2;
+        tmp1 = i+1;
+        tmp2 = 16 - i - 1;
+        tmp3 = (i>=8)? tmp2 : tmp1;
+        tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4;
+        jac_row[0] = -1;
+        jac_row[1] = tmp1*tmp2/tmp4;
+        jac_row[2] = tmp1*tmp3/tmp4;
+        return 0;
+    }
 };
 
 void testLmstr1()
@@ -495,7 +489,8 @@ void testLmstr()
 }
 
 struct lmdif_functor {
-    static int f(const VectorXd &x, VectorXd &fvec, int /*iflag*/)
+    static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */) { return 0;}
+    static int f(const VectorXd &x, VectorXd &fvec)
     {
         /* function fcn for lmdif1 example */