diff --git a/unsupported/Eigen/NonLinear b/unsupported/Eigen/NonLinear
index cc3e2d147..3fa11d214 100644
--- a/unsupported/Eigen/NonLinear
+++ b/unsupported/Eigen/NonLinear
@@ -36,18 +36,6 @@ namespace Eigen {
   */
 //@{
 
-#define min(a,b) ((a) <= (b) ? (a) : (b))
-#define max(a,b) ((a) >= (b) ? (a) : (b))
-#define abs(x) ((x) >= 0 ? (x) : -(x))
-#define TRUE_ (1)
-#define FALSE_ (0)
-#define p1 .1
-#define p5 .5
-#define p25 .25
-#define p75 .75
-#define p001 .001
-#define p0001 1e-4
-
 #if 1
 #include <cminpack.h>
 #else
diff --git a/unsupported/Eigen/src/NonLinear/hybrd.h b/unsupported/Eigen/src/NonLinear/hybrd.h
index 626130f38..7b286c8d5 100644
--- a/unsupported/Eigen/src/NonLinear/hybrd.h
+++ b/unsupported/Eigen/src/NonLinear/hybrd.h
@@ -86,7 +86,7 @@ L20:
     /*     the jacobian matrix. */
 
     /* Computing MIN */
-    msum = min(nb_of_subdiagonals + nb_of_superdiagonals + 1, n);
+    msum = std::min(nb_of_subdiagonals + nb_of_superdiagonals + 1, n);
 
     /*     initialize iteration counter and monitors. */
 
@@ -99,7 +99,7 @@ L20:
     /*     beginning of the outer loop. */
 
 L30:
-    jeval = TRUE_;
+    jeval = true;
 
     /*        calculate the jacobian matrix. */
 
@@ -112,7 +112,7 @@ L30:
 
     /*        compute the qr factorization of the jacobian. */
 
-    qrfac(n, n, fjac.data(), ldfjac, FALSE_, iwa, 1, wa1.data(), wa2.data(), wa3.data());
+    qrfac(n, n, fjac.data(), ldfjac, 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. */
@@ -173,7 +173,7 @@ L110:
 
     /*        copy the triangular factor of the qr factorization into r. */
 
-    sing = FALSE_;
+    sing = false;
     for (j = 0; j < n; ++j) {
         l = j;
         if (j) {
@@ -185,7 +185,7 @@ L110:
         }
         R[l] = wa1[j];
         if (wa1[j] == 0.) {
-            sing = TRUE_;
+            sing = true;
         }
         /* L150: */
     }
@@ -201,7 +201,7 @@ L110:
     }
     /* Computing MAX */
     for (j = 0; j < n; ++j)
-        diag[j] = max(diag[j], wa2[j]);
+        diag[j] = std::max(diag[j], wa2[j]);
 L170:
 
     /*        beginning of the inner loop. */
@@ -224,7 +224,7 @@ L190:
 
     /*           determine the direction p. */
 
-    dogleg(n, R.data(), lr, diag.data(), qtf.data(), delta, wa1.data(), wa2.data(), wa3.data());
+    ei_dogleg<Scalar>(n, R.data(), lr, diag.data(), qtf.data(), delta, wa1.data(), wa2.data(), wa3.data());
 
     /*           store the direction p and x + p. calculate the norm of p. */
 
@@ -239,7 +239,7 @@ L190:
     /*           on the first iteration, adjust the initial step bound. */
 
     if (iter == 1) {
-        delta = min(delta,pnorm);
+        delta = std::min(delta,pnorm);
     }
 
     /*           evaluate the function at x + p and calculate its norm. */
@@ -285,26 +285,26 @@ L190:
 
     /*           update the step bound. */
 
-    if (ratio >= p1) {
+    if (ratio >= Scalar(.1)) {
         goto L230;
     }
     ncsuc = 0;
     ++ncfail;
-    delta = p5 * delta;
+    delta = Scalar(.5) * delta;
     goto L240;
 L230:
     ncfail = 0;
     ++ncsuc;
-    if (ratio >= p5 || ncsuc > 1) /* Computing MAX */
-        delta = max(delta, pnorm / p5);
-    if (ei_abs(ratio - 1.) <= p1) {
-        delta = pnorm / p5;
+    if (ratio >= Scalar(.5) || ncsuc > 1) /* Computing MAX */
+        delta = std::max(delta, pnorm / Scalar(.5));
+    if (ei_abs(ratio - 1.) <= Scalar(.1)) {
+        delta = pnorm / Scalar(.5);
     }
 L240:
 
     /*           test for successful iteration. */
 
-    if (ratio < p0001) {
+    if (ratio < Scalar(1e-4)) {
         goto L260;
     }
 
@@ -324,13 +324,13 @@ L260:
     /*           determine the progress of the iteration. */
 
     ++nslow1;
-    if (actred >= p001) {
+    if (actred >= Scalar(.001)) {
         nslow1 = 0;
     }
     if (jeval) {
         ++nslow2;
     }
-    if (actred >= p1) {
+    if (actred >= Scalar(.1)) {
         nslow2 = 0;
     }
 
@@ -349,7 +349,7 @@ L260:
         info = 2;
     }
     /* Computing MAX */
-    if (p1 * max(p1 * delta, pnorm) <= epsilon<Scalar>() * xnorm)
+    if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm)
         info = 3;
     if (nslow2 == 5)
         info = 4;
@@ -375,7 +375,7 @@ L260:
         }
         wa2[j] = (sum - wa3[j]) / pnorm;
         wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
-        if (ratio >= p0001) {
+        if (ratio >= Scalar(1e-4)) {
             qtf[j] = sum;
         }
         /* L280: */
@@ -389,7 +389,7 @@ L260:
 
     /*           end of the inner loop. */
 
-    jeval = FALSE_;
+    jeval = false;
     goto L180;
 L290:
 
diff --git a/unsupported/Eigen/src/NonLinear/hybrj.h b/unsupported/Eigen/src/NonLinear/hybrj.h
index 8e0945894..ceb8e8f62 100644
--- a/unsupported/Eigen/src/NonLinear/hybrj.h
+++ b/unsupported/Eigen/src/NonLinear/hybrj.h
@@ -88,7 +88,7 @@ L20:
     /*     beginning of the outer loop. */
 
 L30:
-    jeval = TRUE_;
+    jeval = true;
 
     /*        calculate the jacobian matrix. */
 
@@ -100,7 +100,7 @@ L30:
 
     /*        compute the qr factorization of the jacobian. */
 
-    qrfac(n, n,fjac.data(), ldfjac, FALSE_, iwa, 1, wa1.data(), wa2.data(), wa3.data());
+    qrfac(n, n,fjac.data(), ldfjac, 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. */
@@ -161,7 +161,7 @@ L110:
 
     /*        copy the triangular factor of the qr factorization into r. */
 
-    sing = FALSE_;
+    sing = false;
     for (j = 0; j < n; ++j) {
         l = j;
         if (j) {
@@ -173,7 +173,7 @@ L110:
         }
         R[l] = wa1[j];
         if (wa1[j] == 0.) {
-            sing = TRUE_;
+            sing = true;
         }
         /* L150: */
     }
@@ -189,7 +189,7 @@ L110:
     }
     /* Computing MAX */
     for (j = 0; j < n; ++j)
-        diag[j] = max(diag[j], wa2[j]);
+        diag[j] = std::max(diag[j], wa2[j]);
 L170:
 
     /*        beginning of the inner loop. */
@@ -212,7 +212,7 @@ L190:
 
     /*           determine the direction p. */
 
-    dogleg(n, R.data(), lr, diag.data(), qtf.data(), delta, wa1.data(), wa2.data(), wa3.data());
+    ei_dogleg<Scalar>(n, R.data(), lr, diag.data(), qtf.data(), delta, wa1.data(), wa2.data(), wa3.data());
 
     /*           store the direction p and x + p. calculate the norm of p. */
 
@@ -227,7 +227,7 @@ L190:
     /*           on the first iteration, adjust the initial step bound. */
 
     if (iter == 1) {
-        delta = min(delta,pnorm);
+        delta = std::min(delta,pnorm);
     }
 
     /*           evaluate the function at x + p and calculate its norm. */
@@ -273,26 +273,26 @@ L190:
 
     /*           update the step bound. */
 
-    if (ratio >= p1) {
+    if (ratio >= Scalar(.1)) {
         goto L230;
     }
     ncsuc = 0;
     ++ncfail;
-    delta = p5 * delta;
+    delta = Scalar(.5) * delta;
     goto L240;
 L230:
     ncfail = 0;
     ++ncsuc;
-    if (ratio >= p5 || ncsuc > 1) /* Computing MAX */
-        delta = max(delta, pnorm / p5);
-    if (ei_abs(ratio - 1.) <= p1) {
-        delta = pnorm / p5;
+    if (ratio >= Scalar(.5) || ncsuc > 1) /* Computing MAX */
+        delta = std::max(delta, pnorm / Scalar(.5));
+    if (ei_abs(ratio - 1.) <= Scalar(.1)) {
+        delta = pnorm / Scalar(.5);
     }
 L240:
 
     /*           test for successful iteration. */
 
-    if (ratio < p0001) {
+    if (ratio < Scalar(1e-4)) {
         goto L260;
     }
 
@@ -312,13 +312,13 @@ L260:
     /*           determine the progress of the iteration. */
 
     ++nslow1;
-    if (actred >= p001) {
+    if (actred >= Scalar(.001)) {
         nslow1 = 0;
     }
     if (jeval) {
         ++nslow2;
     }
-    if (actred >= p1) {
+    if (actred >= Scalar(.1)) {
         nslow2 = 0;
     }
 
@@ -337,7 +337,7 @@ L260:
         info = 2;
     }
     /* Computing MAX */
-    if (p1 * max(p1 * delta, pnorm) <= epsilon<Scalar>() * xnorm) {
+    if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm) {
         info = 3;
     }
     if (nslow2 == 5) {
@@ -367,7 +367,7 @@ L260:
         }
         wa2[j] = (sum - wa3[j]) / pnorm;
         wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
-        if (ratio >= p0001) {
+        if (ratio >= Scalar(1e-4)) {
             qtf[j] = sum;
         }
         /* L280: */
@@ -381,7 +381,7 @@ L260:
 
     /*           end of the inner loop. */
 
-    jeval = FALSE_;
+    jeval = false;
     goto L180;
 L290:
 
diff --git a/unsupported/Eigen/src/NonLinear/lmder.h b/unsupported/Eigen/src/NonLinear/lmder.h
index 4499b381e..2a8635661 100644
--- a/unsupported/Eigen/src/NonLinear/lmder.h
+++ b/unsupported/Eigen/src/NonLinear/lmder.h
@@ -104,7 +104,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());
+    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 */
@@ -183,7 +183,7 @@ L120:
             /* L140: */
         }
         /* Computing MAX */
-        gnorm = max(gnorm, ei_abs(sum / wa2[l]));
+        gnorm = std::max(gnorm, ei_abs(sum / wa2[l]));
 L150:
         /* L160: */
         ;
@@ -205,7 +205,7 @@ L170:
         goto L190;
     }
     for (j = 0; j < n; ++j)
-        diag[j] = max( diag[j], wa2[j]);
+        diag[j] = std::max( diag[j], wa2[j]);
 L190:
 
     /*        beginning of the inner loop. */
@@ -215,7 +215,7 @@ 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,
+    ei_lmpar<Scalar>(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta,
             &par, wa1.data(), wa2.data(), wa3.data(), wa4.data());
     ipvt.cwise()-=1;
 
@@ -232,7 +232,7 @@ L200:
     /*           on the first iteration, adjust the initial step bound. */
 
     if (iter == 1) {
-        delta = min(delta,pnorm);
+        delta = std::min(delta,pnorm);
     }
 
     /*           evaluate the function at x + p and calculate its norm. */
@@ -247,7 +247,7 @@ L200:
     /*           compute the scaled actual reduction. */
 
     actred = -1.;
-    if (p1 * fnorm1 < fnorm) /* Computing 2nd power */
+    if (Scalar(.1) * fnorm1 < fnorm) /* Computing 2nd power */
         actred = 1. - ei_abs2(fnorm1 / fnorm);
 
     /*           compute the scaled predicted reduction and */
@@ -266,7 +266,7 @@ L200:
     temp1 = ei_abs2(wa3.stableNorm() / fnorm);
     temp2 = ei_abs2( ei_sqrt(par) * pnorm / fnorm);
     /* Computing 2nd power */
-    prered = temp1 + temp2 / p5;
+    prered = temp1 + temp2 / Scalar(.5);
     dirder = -(temp1 + temp2);
 
     /*           compute the ratio of the actual to the predicted */
@@ -279,34 +279,34 @@ L200:
 
     /*           update the step bound. */
 
-    if (ratio > p25) {
+    if (ratio > Scalar(.25)) {
         goto L240;
     }
     if (actred >= 0.) {
-        temp = p5;
+        temp = Scalar(.5);
     }
     if (actred < 0.) {
-        temp = p5 * dirder / (dirder + p5 * actred);
+        temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
     }
-    if (p1 * fnorm1 >= fnorm || temp < p1) {
-        temp = p1;
+    if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) {
+        temp = Scalar(.1);
     }
     /* Computing MIN */
-    delta = temp * min(delta, pnorm/p1);
+    delta = temp * std::min(delta, pnorm/Scalar(.1));
     par /= temp;
     goto L260;
 L240:
-    if (par != 0. && ratio < p75) {
+    if (par != 0. && ratio < Scalar(.75)) {
         goto L250;
     }
-    delta = pnorm / p5;
-    par = p5 * par;
+    delta = pnorm / Scalar(.5);
+    par = Scalar(.5) * par;
 L250:
 L260:
 
     /*           test for successful iteration. */
 
-    if (ratio < p0001) {
+    if (ratio < Scalar(1e-4)) {
         goto L290;
     }
 
@@ -328,13 +328,13 @@ L290:
 
     /*           tests for convergence. */
 
-    if (ei_abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1.) {
+    if (ei_abs(actred) <= ftol && prered <= ftol && Scalar(.5) * ratio <= 1.) {
         info = 1;
     }
     if (delta <= xtol * xnorm) {
         info = 2;
     }
-    if (ei_abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1. && info 
+    if (ei_abs(actred) <= ftol && prered <= ftol && Scalar(.5) * ratio <= 1. && info 
             == 2) {
         info = 3;
     }
@@ -347,7 +347,7 @@ L290:
     if (nfev >= maxfev) {
         info = 5;
     }
-    if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && p5 * ratio <= 1.) {
+    if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && Scalar(.5) * ratio <= 1.) {
         info = 6;
     }
     if (delta <= epsilon<Scalar>() * xnorm) {
@@ -362,7 +362,7 @@ L290:
 
     /*           end of the inner loop. repeat if iteration unsuccessful. */
 
-    if (ratio < p0001) {
+    if (ratio < Scalar(1e-4)) {
         goto L200;
     }
 
diff --git a/unsupported/Eigen/src/NonLinear/lmdif.h b/unsupported/Eigen/src/NonLinear/lmdif.h
index ab97794ca..c12e3efa7 100644
--- a/unsupported/Eigen/src/NonLinear/lmdif.h
+++ b/unsupported/Eigen/src/NonLinear/lmdif.h
@@ -107,7 +107,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());
+    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 */
@@ -186,7 +186,7 @@ L120:
             /* L140: */
         }
         /* Computing MAX */
-        gnorm = max(gnorm, ei_abs(sum / wa2[l]));
+        gnorm = std::max(gnorm, ei_abs(sum / wa2[l]));
 L150:
         /* L160: */
         ;
@@ -208,7 +208,7 @@ L170:
         goto L190;
     }
     for (j = 0; j < n; ++j) /* Computing MAX */
-        diag[j] = max(diag[j], wa2[j]);
+        diag[j] = std::max(diag[j], wa2[j]);
 L190:
 
     /*        beginning of the inner loop. */
@@ -218,7 +218,7 @@ 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,
+    ei_lmpar<Scalar>(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta,
             &par, wa1.data(), wa2.data(), wa3.data(), wa4.data());
     ipvt.cwise()-=1;
 
@@ -235,7 +235,7 @@ L200:
     /*           on the first iteration, adjust the initial step bound. */
 
     if (iter == 1) {
-        delta = min(delta,pnorm);
+        delta = std::min(delta,pnorm);
     }
 
     /*           evaluate the function at x + p and calculate its norm. */
@@ -250,7 +250,7 @@ L200:
     /*           compute the scaled actual reduction. */
 
     actred = -1.;
-    if (p1 * fnorm1 < fnorm) /* Computing 2nd power */
+    if (Scalar(.1) * fnorm1 < fnorm) /* Computing 2nd power */
         actred = 1. - ei_abs2(fnorm1 / fnorm);
 
     /*           compute the scaled predicted reduction and */
@@ -269,7 +269,7 @@ L200:
     temp1 = ei_abs2(wa3.stableNorm() / fnorm);
     temp2 = ei_abs2(ei_sqrt(par) * pnorm / fnorm);
     /* Computing 2nd power */
-    prered = temp1 + temp2 / p5;
+    prered = temp1 + temp2 / Scalar(.5);
     dirder = -(temp1 + temp2);
 
     /*           compute the ratio of the actual to the predicted */
@@ -282,34 +282,34 @@ L200:
 
     /*           update the step bound. */
 
-    if (ratio > p25) {
+    if (ratio > Scalar(.25)) {
         goto L240;
     }
     if (actred >= 0.) {
-        temp = p5;
+        temp = Scalar(.5);
     }
     if (actred < 0.) {
-        temp = p5 * dirder / (dirder + p5 * actred);
+        temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
     }
-    if (p1 * fnorm1 >= fnorm || temp < p1) {
-        temp = p1;
+    if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) {
+        temp = Scalar(.1);
     }
     /* Computing MIN */
-    delta = temp * min(delta, pnorm / p1);
+    delta = temp * std::min(delta, pnorm / Scalar(.1));
     par /= temp;
     goto L260;
 L240:
-    if (par != 0. && ratio < p75) {
+    if (par != 0. && ratio < Scalar(.75)) {
         goto L250;
     }
-    delta = pnorm / p5;
-    par = p5 * par;
+    delta = pnorm / Scalar(.5);
+    par = Scalar(.5) * par;
 L250:
 L260:
 
     /*           test for successful iteration. */
 
-    if (ratio < p0001) {
+    if (ratio < Scalar(1e-4)) {
         goto L290;
     }
 
@@ -331,13 +331,13 @@ L290:
 
     /*           tests for convergence. */
 
-    if (ei_abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1.) {
+    if (ei_abs(actred) <= ftol && prered <= ftol && Scalar(.5) * ratio <= 1.) {
         info = 1;
     }
     if (delta <= xtol * xnorm) {
         info = 2;
     }
-    if (ei_abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1. && info 
+    if (ei_abs(actred) <= ftol && prered <= ftol && Scalar(.5) * ratio <= 1. && info 
             == 2) {
         info = 3;
     }
@@ -350,7 +350,7 @@ L290:
     if (nfev >= maxfev) {
         info = 5;
     }
-    if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && p5 * ratio <= 1.) {
+    if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && Scalar(.5) * ratio <= 1.) {
         info = 6;
     }
     if (delta <= epsilon<Scalar>() * xnorm) {
@@ -365,7 +365,7 @@ L290:
 
     /*           end of the inner loop. repeat if iteration unsuccessful. */
 
-    if (ratio < p0001) {
+    if (ratio < Scalar(1e-4)) {
         goto L200;
     }
 
diff --git a/unsupported/Eigen/src/NonLinear/lmstr.h b/unsupported/Eigen/src/NonLinear/lmstr.h
index a221b897c..db8ec1562 100644
--- a/unsupported/Eigen/src/NonLinear/lmstr.h
+++ b/unsupported/Eigen/src/NonLinear/lmstr.h
@@ -124,10 +124,10 @@ L40:
     /*        if the jacobian is rank deficient, call qrfac to */
     /*        reorder its columns and update the components of qtf. */
 
-    sing = FALSE_;
+    sing = false;
     for (j = 0; j < n; ++j) {
         if (fjac(j,j) == 0.) {
-            sing = TRUE_;
+            sing = true;
         }
         ipvt[j] = j;
         wa2[j] = fjac.col(j).start(j).stableNorm();
@@ -139,7 +139,7 @@ L40:
         goto L130;
     }
     ipvt.cwise()+=1;
-    qrfac(n, n, fjac.data(), ldfjac, TRUE_, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data());
+    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.) {
@@ -210,7 +210,7 @@ L170:
             /* L180: */
         }
         /* Computing MAX */
-        gnorm = max(gnorm, ei_abs(sum/wa2[l]));
+        gnorm = std::max(gnorm, ei_abs(sum/wa2[l]));
 L190:
         /* L200: */
         ;
@@ -232,7 +232,7 @@ L210:
         goto L230;
     }
     for (j = 0; j < n; ++j) /* Computing MAX */
-        diag[j] = max(diag[j], wa2[j]);
+        diag[j] = std::max(diag[j], wa2[j]);
 L230:
 
     /*        beginning of the inner loop. */
@@ -242,7 +242,7 @@ 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,
+    ei_lmpar<Scalar>(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta, &par,
             wa1.data(), wa2.data(), wa3.data(), wa4.data());
     ipvt.cwise()-=1;
 
@@ -259,7 +259,7 @@ L240:
     /*           on the first iteration, adjust the initial step bound. */
 
     if (iter == 1) {
-        delta = min(delta,pnorm);
+        delta = std::min(delta,pnorm);
     }
 
     /*           evaluate the function at x + p and calculate its norm. */
@@ -274,7 +274,7 @@ L240:
     /*           compute the scaled actual reduction. */
 
     actred = -1.;
-    if (p1 * fnorm1 < fnorm) /* Computing 2nd power */
+    if (Scalar(.1) * fnorm1 < fnorm) /* Computing 2nd power */
         actred = 1. - ei_abs2(fnorm1 / fnorm);
 
     /*           compute the scaled predicted reduction and */
@@ -293,7 +293,7 @@ L240:
     temp1 = ei_abs2(wa3.stableNorm() / fnorm);
     temp2 = ei_abs2( ei_sqrt(par) * pnorm / fnorm);
     /* Computing 2nd power */
-    prered = temp1 + temp2 / p5;
+    prered = temp1 + temp2 / Scalar(.5);
     dirder = -(temp1 + temp2);
 
     /*           compute the ratio of the actual to the predicted */
@@ -306,35 +306,35 @@ L240:
 
     /*           update the step bound. */
 
-    if (ratio > p25) {
+    if (ratio > Scalar(.25)) {
         goto L280;
     }
     if (actred >= 0.) {
-        temp = p5;
+        temp = Scalar(.5);
     }
     if (actred < 0.) {
-        temp = p5 * dirder / (dirder + p5 * actred);
+        temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
     }
-    if (p1 * fnorm1 >= fnorm || temp < p1) {
-        temp = p1;
+    if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) {
+        temp = Scalar(.1);
     }
     /* Computing MIN */
-    delta = temp * min(delta, pnorm / p1);
+    delta = temp * std::min(delta, pnorm / Scalar(.1));
 
     par /= temp;
     goto L300;
 L280:
-    if (par != 0. && ratio < p75) {
+    if (par != 0. && ratio < Scalar(.75)) {
         goto L290;
     }
-    delta = pnorm / p5;
-    par = p5 * par;
+    delta = pnorm / Scalar(.5);
+    par = Scalar(.5) * par;
 L290:
 L300:
 
     /*           test for successful iteration. */
 
-    if (ratio < p0001) {
+    if (ratio < Scalar(1e-4)) {
         goto L330;
     }
 
@@ -356,13 +356,13 @@ L330:
 
     /*           tests for convergence. */
 
-    if (ei_abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1.) {
+    if (ei_abs(actred) <= ftol && prered <= ftol && Scalar(.5) * ratio <= 1.) {
         info = 1;
     }
     if (delta <= xtol * xnorm) {
         info = 2;
     }
-    if (ei_abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1. && info 
+    if (ei_abs(actred) <= ftol && prered <= ftol && Scalar(.5) * ratio <= 1. && info 
             == 2) {
         info = 3;
     }
@@ -375,7 +375,7 @@ L330:
     if (nfev >= maxfev) {
         info = 5;
     }
-    if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && p5 * ratio <= 1.) {
+    if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && Scalar(.5) * ratio <= 1.) {
         info = 6;
     }
     if (delta <= epsilon<Scalar>() * xnorm) {
@@ -390,7 +390,7 @@ L330:
 
     /*           end of the inner loop. repeat if iteration unsuccessful. */
 
-    if (ratio < p0001) {
+    if (ratio < Scalar(1e-4)) {
         goto L240;
     }