diff --git a/unsupported/Eigen/NonLinear b/unsupported/Eigen/NonLinear
index 3fa11d214..79c7e4c4c 100644
--- a/unsupported/Eigen/NonLinear
+++ b/unsupported/Eigen/NonLinear
@@ -50,6 +50,17 @@ typedef int (*minpack_funcderstr_mn)(void *p, int m, int n, const double *x, dou
 void covar(int n, double *r__, int ldr, const int *ipvt, double tol, double *wa);
 #endif
 
+#include "src/NonLinear/qrsolv.h"
+#include "src/NonLinear/r1updt.h"
+#include "src/NonLinear/r1mpyq.h"
+#include "src/NonLinear/rwupdt.h"
+#include "src/NonLinear/qrfac.h"
+#include "src/NonLinear/fdjac2.h"
+#include "src/NonLinear/fdjac1.h"
+#include "src/NonLinear/qform.h"
+#include "src/NonLinear/lmpar.h"
+#include "src/NonLinear/dogleg.h"
+#include "src/NonLinear/covar.h"
 
 #include "src/NonLinear/lmder.h"
 #include "src/NonLinear/hybrd.h"
diff --git a/unsupported/Eigen/src/NonLinear/dogleg.h b/unsupported/Eigen/src/NonLinear/dogleg.h
new file mode 100644
index 000000000..0909b9621
--- /dev/null
+++ b/unsupported/Eigen/src/NonLinear/dogleg.h
@@ -0,0 +1,181 @@
+
+template <typename Scalar>
+void ei_dogleg(int n, const Scalar *r__, int /* lr*/ , 
+	const Scalar *diag, const Scalar *qtb, Scalar delta, Scalar *x, 
+	Scalar *wa1, Scalar *wa2)
+{
+    /* System generated locals */
+    int i__1, i__2;
+    Scalar d__1, d__2, d__3, d__4;
+
+    /* Local variables */
+    int i__, j, k, l, jj, jp1;
+    Scalar sum, temp, alpha, bnorm;
+    Scalar gnorm, qnorm, epsmch;
+    Scalar sgnorm;
+
+    /* Parameter adjustments */
+    --wa2;
+    --wa1;
+    --x;
+    --qtb;
+    --diag;
+    --r__;
+
+    /* Function Body */
+
+/*     epsmch is the machine precision. */
+
+    epsmch = epsilon<Scalar>();
+
+/*     first, calculate the gauss-newton direction. */
+
+    jj = n * (n + 1) / 2 + 1;
+    i__1 = n;
+    for (k = 1; k <= i__1; ++k) {
+	j = n - k + 1;
+	jp1 = j + 1;
+	jj -= k;
+	l = jj + 1;
+	sum = 0.;
+	if (n < jp1) {
+	    goto L20;
+	}
+	i__2 = n;
+	for (i__ = jp1; i__ <= i__2; ++i__) {
+	    sum += r__[l] * x[i__];
+	    ++l;
+/* L10: */
+	}
+L20:
+	temp = r__[jj];
+	if (temp != 0.) {
+	    goto L40;
+	}
+	l = j;
+	i__2 = j;
+	for (i__ = 1; i__ <= i__2; ++i__) {
+/* Computing MAX */
+	    d__2 = temp, d__3 = fabs(r__[l]);
+	    temp = std::max(d__2,d__3);
+	    l = l + n - i__;
+/* L30: */
+	}
+	temp = epsmch * temp;
+	if (temp == 0.) {
+	    temp = epsmch;
+	}
+L40:
+	x[j] = (qtb[j] - sum) / temp;
+/* L50: */
+    }
+
+/*     test whether the gauss-newton direction is acceptable. */
+
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	wa1[j] = 0.;
+	wa2[j] = diag[j] * x[j];
+/* L60: */
+    }
+    qnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&wa2[1],n).stableNorm();
+    if (qnorm <= delta) {
+	/* goto L140; */
+        return;
+    }
+
+/*     the gauss-newton direction is not acceptable. */
+/*     next, calculate the scaled gradient direction. */
+
+    l = 1;
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	temp = qtb[j];
+	i__2 = n;
+	for (i__ = j; i__ <= i__2; ++i__) {
+	    wa1[i__] += r__[l] * temp;
+	    ++l;
+/* L70: */
+	}
+	wa1[j] /= diag[j];
+/* L80: */
+    }
+
+/*     calculate the norm of the scaled gradient and test for */
+/*     the special case in which the scaled gradient is zero. */
+
+    gnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&wa1[1],n).stableNorm();
+    sgnorm = 0.;
+    alpha = delta / qnorm;
+    if (gnorm == 0.) {
+	goto L120;
+    }
+
+/*     calculate the point along the scaled gradient */
+/*     at which the quadratic is minimized. */
+
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	wa1[j] = wa1[j] / gnorm / diag[j];
+/* L90: */
+    }
+    l = 1;
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	sum = 0.;
+	i__2 = n;
+	for (i__ = j; i__ <= i__2; ++i__) {
+	    sum += r__[l] * wa1[i__];
+	    ++l;
+/* L100: */
+	}
+	wa2[j] = sum;
+/* L110: */
+    }
+    temp = Map< Matrix< Scalar, Dynamic, 1 > >(&wa2[1],n).stableNorm();
+    sgnorm = gnorm / temp / temp;
+
+/*     test whether the scaled gradient direction is acceptable. */
+
+    alpha = 0.;
+    if (sgnorm >= delta) {
+	goto L120;
+    }
+
+/*     the scaled gradient direction is not acceptable. */
+/*     finally, calculate the point along the dogleg */
+/*     at which the quadratic is minimized. */
+
+    bnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&qtb[1],n).stableNorm();
+    temp = bnorm / gnorm * (bnorm / qnorm) * (sgnorm / delta);
+/* Computing 2nd power */
+    d__1 = sgnorm / delta;
+/* Computing 2nd power */
+    d__2 = temp - delta / qnorm;
+/* Computing 2nd power */
+    d__3 = delta / qnorm;
+/* Computing 2nd power */
+    d__4 = sgnorm / delta;
+    temp = temp - delta / qnorm * (d__1 * d__1) + sqrt(d__2 * d__2 + (1. - 
+	    d__3 * d__3) * (1. - d__4 * d__4));
+/* Computing 2nd power */
+    d__1 = sgnorm / delta;
+    alpha = delta / qnorm * (1. - d__1 * d__1) / temp;
+L120:
+
+/*     form appropriate convex combination of the gauss-newton */
+/*     direction and the scaled gradient direction. */
+
+    temp = (1. - alpha) * std::min(sgnorm,delta);
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	x[j] = temp * wa1[j] + alpha * x[j];
+/* L130: */
+    }
+/* L140: */
+    return;
+
+/*     last card of subroutine dogleg. */
+
+} /* dogleg_ */
+
diff --git a/unsupported/Eigen/src/NonLinear/fdjac1.h b/unsupported/Eigen/src/NonLinear/fdjac1.h
new file mode 100644
index 000000000..e62bfe8d9
--- /dev/null
+++ b/unsupported/Eigen/src/NonLinear/fdjac1.h
@@ -0,0 +1,113 @@
+
+template <typename Scalar>
+int ei_fdjac1(minpack_func_nn fcn, void *p, int n, Scalar *x, const Scalar *
+	fvec, Scalar *fjac, int ldfjac, int ml, 
+	int mu, Scalar epsfcn, Scalar *wa1, Scalar *wa2)
+{
+    /* System generated locals */
+    int fjac_dim1, fjac_offset, i__1, i__2, i__3, i__4;
+
+    /* Local variables */
+    Scalar h__;
+    int i__, j, k;
+    Scalar eps, temp;
+    int msum;
+    Scalar epsmch;
+    int iflag = 0;
+
+    /* Parameter adjustments */
+    --wa2;
+    --wa1;
+    --fvec;
+    --x;
+    fjac_dim1 = ldfjac;
+    fjac_offset = 1 + fjac_dim1 * 1;
+    fjac -= fjac_offset;
+
+    /* Function Body */
+
+/*     epsmch is the machine precision. */
+
+    epsmch = epsilon<Scalar>();
+
+    eps = ei_sqrt((std::max(epsfcn,epsmch)));
+    msum = ml + mu + 1;
+    if (msum < n) {
+	goto L40;
+    }
+
+/*        computation of dense approximate jacobian. */
+
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	temp = x[j];
+	h__ = eps * ei_abs(temp);
+	if (h__ == 0.) {
+	    h__ = eps;
+	}
+	x[j] = temp + h__;
+	iflag = (*fcn)(p, n, &x[1], &wa1[1], 1);
+	if (iflag < 0) {
+	    goto L30;
+	}
+	x[j] = temp;
+	i__2 = n;
+	for (i__ = 1; i__ <= i__2; ++i__) {
+	    fjac[i__ + j * fjac_dim1] = (wa1[i__] - fvec[i__]) / h__;
+/* L10: */
+	}
+/* L20: */
+    }
+L30:
+    /* goto L110; */
+    return iflag;
+L40:
+
+/*        computation of banded approximate jacobian. */
+
+    i__1 = msum;
+    for (k = 1; k <= i__1; ++k) {
+	i__2 = n;
+	i__3 = msum;
+	for (j = k; i__3 < 0 ? j >= i__2 : j <= i__2; j += i__3) {
+	    wa2[j] = x[j];
+	    h__ = eps * ei_abs(wa2[j]);
+	    if (h__ == 0.) {
+		h__ = eps;
+	    }
+	    x[j] = wa2[j] + h__;
+/* L60: */
+	}
+	iflag = (*fcn)(p, n, &x[1], &wa1[1], 1);
+	if (iflag < 0) {
+	    /* goto L100; */
+            return iflag;
+	}
+	i__3 = n;
+	i__2 = msum;
+	for (j = k; i__2 < 0 ? j >= i__3 : j <= i__3; j += i__2) {
+	    x[j] = wa2[j];
+	    h__ = eps * ei_abs(wa2[j]);
+	    if (h__ == 0.) {
+		h__ = eps;
+	    }
+	    i__4 = n;
+	    for (i__ = 1; i__ <= i__4; ++i__) {
+		fjac[i__ + j * fjac_dim1] = 0.;
+		if (i__ >= j - mu && i__ <= j + ml) {
+		    fjac[i__ + j * fjac_dim1] = (wa1[i__] - fvec[i__]) / h__;
+		}
+/* L70: */
+	    }
+/* L80: */
+	}
+/* L90: */
+    }
+/* L100: */
+/* L110: */
+    return iflag;
+
+/*     last card of subroutine fdjac1. */
+
+} /* fdjac1_ */
+
diff --git a/unsupported/Eigen/src/NonLinear/fdjac2.h b/unsupported/Eigen/src/NonLinear/fdjac2.h
new file mode 100644
index 000000000..324fbc0fe
--- /dev/null
+++ b/unsupported/Eigen/src/NonLinear/fdjac2.h
@@ -0,0 +1,58 @@
+
+template <typename Scalar>
+int ei_fdjac2(minpack_func_mn fcn, void *p, int m, int n, Scalar *x, 
+	const Scalar *fvec, Scalar *fjac, int ldfjac,
+	Scalar epsfcn, Scalar *wa)
+{
+    /* System generated locals */
+    int fjac_dim1, fjac_offset, i__1, i__2;
+
+    /* Local variables */
+    Scalar h__;
+    int i__, j;
+    Scalar eps, temp, epsmch;
+    int iflag;
+
+    /* Parameter adjustments */
+    --wa;
+    --fvec;
+    --x;
+    fjac_dim1 = ldfjac;
+    fjac_offset = 1 + fjac_dim1 * 1;
+    fjac -= fjac_offset;
+
+    /* Function Body */
+
+/*     epsmch is the machine precision. */
+
+    epsmch = epsilon<Scalar>();
+
+    eps = ei_sqrt((std::max(epsfcn,epsmch)));
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	temp = x[j];
+	h__ = eps * ei_abs(temp);
+	if (h__ == 0.) {
+	    h__ = eps;
+	}
+	x[j] = temp + h__;
+	iflag = (*fcn)(p, m, n, &x[1], &wa[1], 1);
+	if (iflag < 0) {
+	    /* goto L30; */
+            return iflag;
+	}
+	x[j] = temp;
+	i__2 = m;
+	for (i__ = 1; i__ <= i__2; ++i__) {
+	    fjac[i__ + j * fjac_dim1] = (wa[i__] - fvec[i__]) / h__;
+/* L10: */
+	}
+/* L20: */
+    }
+/* L30: */
+    return iflag;
+
+/*     last card of subroutine fdjac2. */
+
+} /* fdjac2_ */
+
diff --git a/unsupported/Eigen/src/NonLinear/hybrd.h b/unsupported/Eigen/src/NonLinear/hybrd.h
index dc4d9636c..445d2eb4e 100644
--- a/unsupported/Eigen/src/NonLinear/hybrd.h
+++ b/unsupported/Eigen/src/NonLinear/hybrd.h
@@ -103,7 +103,7 @@ L30:
 
     /*        calculate the jacobian matrix. */
 
-    iflag = fdjac1(Functor::f, 0, n, x.data(), fvec.data(), fjac.data(), ldfjac,
+    iflag = ei_fdjac1<Scalar>(Functor::f, 0, n, x.data(), fvec.data(), fjac.data(), ldfjac,
             nb_of_subdiagonals, nb_of_superdiagonals, epsfcn, wa1.data(), wa2.data());
     nfev += msum;
     if (iflag < 0) {
@@ -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());
+    ei_qrfac<Scalar>(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. */
@@ -192,7 +192,7 @@ L110:
 
     /*        accumulate the orthogonal factor in fjac. */
 
-    qform(n, n, fjac.data(), ldfjac, wa1.data());
+    ei_qform<Scalar>(n, n, fjac.data(), ldfjac, wa1.data());
 
     /*        rescale if necessary. */
 
@@ -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. */
 
@@ -383,9 +383,9 @@ L260:
 
     /*           compute the qr factorization of the updated jacobian. */
 
-    r1updt(n, n, R.data(), lr, wa1.data(), wa2.data(), wa3.data(), &sing);
-    r1mpyq(n, n, fjac.data(), ldfjac, wa2.data(), wa3.data());
-    r1mpyq(1, n, qtf.data(), 1, wa2.data(), wa3.data());
+    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>(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 92cd02efe..e2cb08257 100644
--- a/unsupported/Eigen/src/NonLinear/hybrj.h
+++ b/unsupported/Eigen/src/NonLinear/hybrj.h
@@ -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());
+    ei_qrfac<Scalar>(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. */
@@ -180,7 +180,7 @@ L110:
 
     /*        accumulate the orthogonal factor in fjac. */
 
-    qform(n, n, fjac.data(), ldfjac, wa1.data());
+    ei_qform<Scalar>(n, n, fjac.data(), ldfjac, wa1.data());
 
     /*        rescale if necessary. */
 
@@ -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. */
 
@@ -375,9 +375,9 @@ L260:
 
     /*           compute the qr factorization of the updated jacobian. */
 
-    r1updt(n, n, R.data(), lr, wa1.data(), wa2.data(), wa3.data(), &sing);
-    r1mpyq(n, n, fjac.data(), ldfjac, wa2.data(), wa3.data());
-    r1mpyq(1, n, qtf.data(), 1, wa2.data(), wa3.data());
+    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>(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 81b9b0745..1d02c3684 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());
+    ei_qrfac<Scalar>(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 */
@@ -215,8 +215,13 @@ L200:
     /*           determine the levenberg-marquardt parameter. */
 
     ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
+#if 1
+    ei_lmpar<Scalar>(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta,
+            &par, wa1.data(), wa2.data(), wa3.data(), wa4.data());
+#else
     lmpar(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta,
             &par, wa1.data(), wa2.data(), wa3.data(), wa4.data());
+#endif
     ipvt.cwise()-=1;
 
     /*           store the direction p and x + p. calculate the norm of p. */
diff --git a/unsupported/Eigen/src/NonLinear/lmdif.h b/unsupported/Eigen/src/NonLinear/lmdif.h
index 1c7e6cfc6..19af2005c 100644
--- a/unsupported/Eigen/src/NonLinear/lmdif.h
+++ b/unsupported/Eigen/src/NonLinear/lmdif.h
@@ -84,7 +84,7 @@ L30:
 
     /*        calculate the jacobian matrix. */
 
-    iflag = fdjac2(Functor::f, 0, m, n, x.data(), fvec.data(), fjac.data(), ldfjac,
+    iflag = ei_fdjac2<Scalar>(Functor::f, 0, m, n, x.data(), fvec.data(), fjac.data(), ldfjac,
             epsfcn, wa4.data());
     nfev += n;
     if (iflag < 0) {
@@ -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());
+    ei_qrfac<Scalar>(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 */
@@ -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;
 
diff --git a/unsupported/Eigen/src/NonLinear/lmpar.h b/unsupported/Eigen/src/NonLinear/lmpar.h
new file mode 100644
index 000000000..2d75b3b15
--- /dev/null
+++ b/unsupported/Eigen/src/NonLinear/lmpar.h
@@ -0,0 +1,264 @@
+
+template <typename Scalar>
+void ei_lmpar(int n, Scalar *r__, int ldr, 
+	const int *ipvt, const Scalar *diag, const Scalar *qtb, Scalar delta, 
+	Scalar *par, Scalar *x, Scalar *sdiag, Scalar *wa1, 
+	Scalar *wa2)
+{
+    /* Initialized data */
+
+#define p1 .1
+#define p001 .001
+
+    /* System generated locals */
+    int r_dim1, r_offset, i__1, i__2;
+    Scalar d__1, d__2;
+
+    /* Local variables */
+    int i__, j, k, l;
+    Scalar fp;
+    int jm1, jp1;
+    Scalar sum, parc, parl;
+    int iter;
+    Scalar temp, paru, dwarf;
+    int nsing;
+    Scalar gnorm;
+    Scalar dxnorm;
+
+    /* Parameter adjustments */
+    --wa2;
+    --wa1;
+    --sdiag;
+    --x;
+    --qtb;
+    --diag;
+    --ipvt;
+    r_dim1 = ldr;
+    r_offset = 1 + r_dim1 * 1;
+    r__ -= r_offset;
+
+    /* Function Body */
+
+/*     dwarf is the smallest positive magnitude. */
+
+    dwarf = std::numeric_limits<Scalar>::min();
+
+/*     compute and store in x the gauss-newton direction. if the */
+/*     jacobian is rank-deficient, obtain a least squares solution. */
+
+    nsing = n;
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	wa1[j] = qtb[j];
+	if (r__[j + j * r_dim1] == 0. && nsing == n) {
+	    nsing = j - 1;
+	}
+	if (nsing < n) {
+	    wa1[j] = 0.;
+	}
+/* L10: */
+    }
+    if (nsing < 1) {
+	goto L50;
+    }
+    i__1 = nsing;
+    for (k = 1; k <= i__1; ++k) {
+	j = nsing - k + 1;
+	wa1[j] /= r__[j + j * r_dim1];
+	temp = wa1[j];
+	jm1 = j - 1;
+	if (jm1 < 1) {
+	    goto L30;
+	}
+	i__2 = jm1;
+	for (i__ = 1; i__ <= i__2; ++i__) {
+	    wa1[i__] -= r__[i__ + j * r_dim1] * temp;
+/* L20: */
+	}
+L30:
+/* L40: */
+	;
+    }
+L50:
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	l = ipvt[j];
+	x[l] = wa1[j];
+/* L60: */
+    }
+
+/*     initialize the iteration counter. */
+/*     evaluate the function at the origin, and test */
+/*     for acceptance of the gauss-newton direction. */
+
+    iter = 0;
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	wa2[j] = diag[j] * x[j];
+/* L70: */
+    }
+    dxnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&wa2[1],n).blueNorm();
+    fp = dxnorm - delta;
+    if (fp <= p1 * delta) {
+	goto L220;
+    }
+
+/*     if the jacobian is not rank deficient, the newton */
+/*     step provides a lower bound, parl, for the zero of */
+/*     the function. otherwise set this bound to zero. */
+
+    parl = 0.;
+    if (nsing < n) {
+	goto L120;
+    }
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	l = ipvt[j];
+	wa1[j] = diag[l] * (wa2[l] / dxnorm);
+/* L80: */
+    }
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	sum = 0.;
+	jm1 = j - 1;
+	if (jm1 < 1) {
+	    goto L100;
+	}
+	i__2 = jm1;
+	for (i__ = 1; i__ <= i__2; ++i__) {
+	    sum += r__[i__ + j * r_dim1] * wa1[i__];
+/* L90: */
+	}
+L100:
+	wa1[j] = (wa1[j] - sum) / r__[j + j * r_dim1];
+/* L110: */
+    }
+    temp = Map< Matrix< Scalar, Dynamic, 1 > >(&wa1[1],n).blueNorm();
+    parl = fp / delta / temp / temp;
+L120:
+
+/*     calculate an upper bound, paru, for the zero of the function. */
+
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	sum = 0.;
+	i__2 = j;
+	for (i__ = 1; i__ <= i__2; ++i__) {
+	    sum += r__[i__ + j * r_dim1] * qtb[i__];
+/* L130: */
+	}
+	l = ipvt[j];
+	wa1[j] = sum / diag[l];
+/* L140: */
+    }
+    gnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&wa1[1],n).stableNorm();
+    paru = gnorm / delta;
+    if (paru == 0.) {
+	paru = dwarf / std::min(delta,p1);
+    }
+
+/*     if the input par lies outside of the interval (parl,paru), */
+/*     set par to the closer endpoint. */
+
+    *par = std::max(*par,parl);
+    *par = std::min(*par,paru);
+    if (*par == 0.) {
+	*par = gnorm / dxnorm;
+    }
+
+/*     beginning of an iteration. */
+
+L150:
+    ++iter;
+
+/*        evaluate the function at the current value of par. */
+
+    if (*par == 0.) {
+/* Computing MAX */
+	d__1 = dwarf, d__2 = p001 * paru;
+	*par = std::max(d__1,d__2);
+    }
+    temp = ei_sqrt(*par);
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	wa1[j] = temp * diag[j];
+/* L160: */
+    }
+    ei_qrsolv<Scalar>(n, &r__[r_offset], ldr, &ipvt[1], &wa1[1], &qtb[1], &x[1], &sdiag[1], &wa2[1]);
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	wa2[j] = diag[j] * x[j];
+/* L170: */
+    }
+    dxnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&wa2[1],n).blueNorm();
+    temp = fp;
+    fp = dxnorm - delta;
+
+/*        if the function is small enough, accept the current value */
+/*        of par. also test for the exceptional cases where parl */
+/*        is zero or the number of iterations has reached 10. */
+
+    if (ei_abs(fp) <= p1 * delta || (parl == 0. && fp <= temp && temp < 0.) ||
+	     iter == 10) {
+	goto L220;
+    }
+
+/*        compute the newton correction. */
+
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	l = ipvt[j];
+	wa1[j] = diag[l] * (wa2[l] / dxnorm);
+/* L180: */
+    }
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	wa1[j] /= sdiag[j];
+	temp = wa1[j];
+	jp1 = j + 1;
+	if (n < jp1) {
+	    goto L200;
+	}
+	i__2 = n;
+	for (i__ = jp1; i__ <= i__2; ++i__) {
+	    wa1[i__] -= r__[i__ + j * r_dim1] * temp;
+/* L190: */
+	}
+L200:
+/* L210: */
+	;
+    }
+    temp = Map< Matrix< Scalar, Dynamic, 1 > >(&wa1[1],n).blueNorm();
+    parc = fp / delta / temp / temp;
+
+/*        depending on the sign of the function, update parl or paru. */
+
+    if (fp > 0.) {
+	parl = std::max(parl,*par);
+    }
+    if (fp < 0.) {
+	paru = std::min(paru,*par);
+    }
+
+/*        compute an improved estimate for par. */
+
+/* Computing MAX */
+    d__1 = parl, d__2 = *par + parc;
+    *par = std::max(d__1,d__2);
+
+/*        end of an iteration. */
+
+    goto L150;
+L220:
+
+/*     termination. */
+
+    if (iter == 0) {
+	*par = 0.;
+    }
+    return;
+
+/*     last card of subroutine lmpar. */
+
+} /* lmpar_ */
+
diff --git a/unsupported/Eigen/src/NonLinear/lmstr.h b/unsupported/Eigen/src/NonLinear/lmstr.h
index 89a2d7002..d05dc1ece 100644
--- a/unsupported/Eigen/src/NonLinear/lmstr.h
+++ b/unsupported/Eigen/src/NonLinear/lmstr.h
@@ -115,7 +115,7 @@ L40:
             goto L340;
         }
         temp = fvec[i];
-        rwupdt(n, fjac.data(), ldfjac, wa3.data(), qtf.data(), &temp, wa1.data(), wa2.data());
+        ei_rwupdt<Scalar>(n, fjac.data(), ldfjac, wa3.data(), qtf.data(), &temp, wa1.data(), wa2.data());
         ++iflag;
         /* L70: */
     }
@@ -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());
+    ei_qrfac<Scalar>(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.) {
@@ -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;
 
diff --git a/unsupported/Eigen/src/NonLinear/qform.h b/unsupported/Eigen/src/NonLinear/qform.h
new file mode 100644
index 000000000..a9920e7b3
--- /dev/null
+++ b/unsupported/Eigen/src/NonLinear/qform.h
@@ -0,0 +1,98 @@
+
+template <typename Scalar>
+void ei_qform(int m, int n, Scalar *q, int
+	ldq, Scalar *wa)
+{
+    /* System generated locals */
+    int q_dim1, q_offset, i__1, i__2, i__3;
+
+    /* Local variables */
+    int i__, j, k, l, jm1, np1;
+    Scalar sum, temp;
+    int minmn;
+
+    /* Parameter adjustments */
+    --wa;
+    q_dim1 = ldq;
+    q_offset = 1 + q_dim1 * 1;
+    q -= q_offset;
+
+    /* Function Body */
+
+/*     zero out upper triangle of q in the first min(m,n) columns. */
+
+    minmn = std::min(m,n);
+    if (minmn < 2) {
+	goto L30;
+    }
+    i__1 = minmn;
+    for (j = 2; j <= i__1; ++j) {
+	jm1 = j - 1;
+	i__2 = jm1;
+	for (i__ = 1; i__ <= i__2; ++i__) {
+	    q[i__ + j * q_dim1] = 0.;
+/* L10: */
+	}
+/* L20: */
+    }
+L30:
+
+/*     initialize remaining columns to those of the identity matrix. */
+
+    np1 = n + 1;
+    if (m < np1) {
+	goto L60;
+    }
+    i__1 = m;
+    for (j = np1; j <= i__1; ++j) {
+	i__2 = m;
+	for (i__ = 1; i__ <= i__2; ++i__) {
+	    q[i__ + j * q_dim1] = 0.;
+/* L40: */
+	}
+	q[j + j * q_dim1] = 1.;
+/* L50: */
+    }
+L60:
+
+/*     accumulate q from its factored form. */
+
+    i__1 = minmn;
+    for (l = 1; l <= i__1; ++l) {
+	k = minmn - l + 1;
+	i__2 = m;
+	for (i__ = k; i__ <= i__2; ++i__) {
+	    wa[i__] = q[i__ + k * q_dim1];
+	    q[i__ + k * q_dim1] = 0.;
+/* L70: */
+	}
+	q[k + k * q_dim1] = 1.;
+	if (wa[k] == 0.) {
+	    goto L110;
+	}
+	i__2 = m;
+	for (j = k; j <= i__2; ++j) {
+	    sum = 0.;
+	    i__3 = m;
+	    for (i__ = k; i__ <= i__3; ++i__) {
+		sum += q[i__ + j * q_dim1] * wa[i__];
+/* L80: */
+	    }
+	    temp = sum / wa[k];
+	    i__3 = m;
+	    for (i__ = k; i__ <= i__3; ++i__) {
+		q[i__ + j * q_dim1] -= temp * wa[i__];
+/* L90: */
+	    }
+/* L100: */
+	}
+L110:
+/* L120: */
+	;
+    }
+    return;
+
+/*     last card of subroutine qform. */
+
+} /* qform_ */
+
diff --git a/unsupported/Eigen/src/NonLinear/qrfac.h b/unsupported/Eigen/src/NonLinear/qrfac.h
new file mode 100644
index 000000000..245603567
--- /dev/null
+++ b/unsupported/Eigen/src/NonLinear/qrfac.h
@@ -0,0 +1,157 @@
+
+template <typename Scalar>
+void ei_qrfac(int m, int n, Scalar *a, int
+	lda, int pivot, int *ipvt, int /* lipvt */, Scalar *rdiag,
+	 Scalar *acnorm, Scalar *wa)
+{
+    /* Initialized data */
+
+#define p05 .05
+
+    /* System generated locals */
+    int a_dim1, a_offset, i__1, i__2, i__3;
+    Scalar d__1, d__2, d__3;
+
+    /* Local variables */
+    int i__, j, k, jp1;
+    Scalar sum;
+    int kmax;
+    Scalar temp;
+    int minmn;
+    Scalar epsmch;
+    Scalar ajnorm;
+
+    /* Parameter adjustments */
+    --wa;
+    --acnorm;
+    --rdiag;
+    a_dim1 = lda;
+    a_offset = 1 + a_dim1 * 1;
+    a -= a_offset;
+    --ipvt;
+
+    /* Function Body */
+
+/*     epsmch is the machine precision. */
+
+    epsmch = epsilon<Scalar>();
+
+/*     compute the initial column norms and initialize several arrays. */
+
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	acnorm[j] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j * a_dim1 + 1],m).blueNorm();
+	rdiag[j] = acnorm[j];
+	wa[j] = rdiag[j];
+	if (pivot) {
+	    ipvt[j] = j;
+	}
+/* L10: */
+    }
+
+/*     reduce a to r with householder transformations. */
+
+    minmn = std::min(m,n);
+    i__1 = minmn;
+    for (j = 1; j <= i__1; ++j) {
+	if (! (pivot)) {
+	    goto L40;
+	}
+
+/*        bring the column of largest norm into the pivot position. */
+
+	kmax = j;
+	i__2 = n;
+	for (k = j; k <= i__2; ++k) {
+	    if (rdiag[k] > rdiag[kmax]) {
+		kmax = k;
+	    }
+/* L20: */
+	}
+	if (kmax == j) {
+	    goto L40;
+	}
+	i__2 = m;
+	for (i__ = 1; i__ <= i__2; ++i__) {
+	    temp = a[i__ + j * a_dim1];
+	    a[i__ + j * a_dim1] = a[i__ + kmax * a_dim1];
+	    a[i__ + kmax * a_dim1] = temp;
+/* L30: */
+	}
+	rdiag[kmax] = rdiag[j];
+	wa[kmax] = wa[j];
+	k = ipvt[j];
+	ipvt[j] = ipvt[kmax];
+	ipvt[kmax] = k;
+L40:
+
+/*        compute the householder transformation to reduce the */
+/*        j-th column of a to a multiple of the j-th unit vector. */
+
+	i__2 = m - j + 1;
+	ajnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j + j * a_dim1],i__2).blueNorm();
+	if (ajnorm == 0.) {
+	    goto L100;
+	}
+	if (a[j + j * a_dim1] < 0.) {
+	    ajnorm = -ajnorm;
+	}
+	i__2 = m;
+	for (i__ = j; i__ <= i__2; ++i__) {
+	    a[i__ + j * a_dim1] /= ajnorm;
+/* L50: */
+	}
+	a[j + j * a_dim1] += 1.;
+
+/*        apply the transformation to the remaining columns */
+/*        and update the norms. */
+
+	jp1 = j + 1;
+	if (n < jp1) {
+	    goto L100;
+	}
+	i__2 = n;
+	for (k = jp1; k <= i__2; ++k) {
+	    sum = 0.;
+	    i__3 = m;
+	    for (i__ = j; i__ <= i__3; ++i__) {
+		sum += a[i__ + j * a_dim1] * a[i__ + k * a_dim1];
+/* L60: */
+	    }
+	    temp = sum / a[j + j * a_dim1];
+	    i__3 = m;
+	    for (i__ = j; i__ <= i__3; ++i__) {
+		a[i__ + k * a_dim1] -= temp * a[i__ + j * a_dim1];
+/* L70: */
+	    }
+	    if (! (pivot) || rdiag[k] == 0.) {
+		goto L80;
+	    }
+	    temp = a[j + k * a_dim1] / rdiag[k];
+/* Computing MAX */
+/* Computing 2nd power */
+	    d__3 = temp;
+	    d__1 = 0., d__2 = 1. - d__3 * d__3;
+	    rdiag[k] *= ei_sqrt((std::max(d__1,d__2)));
+/* Computing 2nd power */
+	    d__1 = rdiag[k] / wa[k];
+	    if (p05 * (d__1 * d__1) > epsmch) {
+		goto L80;
+	    }
+	    i__3 = m - j;
+	    rdiag[k] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[jp1 + k * a_dim1],i__3).blueNorm();
+	    wa[k] = rdiag[k];
+L80:
+/* L90: */
+	    ;
+	}
+L100:
+	rdiag[j] = -ajnorm;
+/* L110: */
+    }
+    return;
+
+/*     last card of subroutine qrfac. */
+
+} /* qrfac_ */
+
diff --git a/unsupported/Eigen/src/NonLinear/qrsolv.h b/unsupported/Eigen/src/NonLinear/qrsolv.h
new file mode 100644
index 000000000..f4795ce95
--- /dev/null
+++ b/unsupported/Eigen/src/NonLinear/qrsolv.h
@@ -0,0 +1,187 @@
+
+template <typename Scalar>
+void ei_qrsolv(int n, Scalar *r__, int ldr, 
+	const int *ipvt, const Scalar *diag, const Scalar *qtb, Scalar *x, 
+	Scalar *sdiag, Scalar *wa)
+{
+    /* Initialized data */
+
+#define p5 .5
+#define p25 .25
+
+    /* System generated locals */
+    int r_dim1, r_offset, i__1, i__2, i__3;
+    Scalar d__1, d__2;
+
+    /* Local variables */
+    int i__, j, k, l, jp1, kp1;
+    Scalar tan__, cos__, sin__, sum, temp, cotan;
+    int nsing;
+    Scalar qtbpj;
+
+    /* Parameter adjustments */
+    --wa;
+    --sdiag;
+    --x;
+    --qtb;
+    --diag;
+    --ipvt;
+    r_dim1 = ldr;
+    r_offset = 1 + r_dim1 * 1;
+    r__ -= r_offset;
+
+    /* Function Body */
+
+/*     copy r and (q transpose)*b to preserve input and initialize s. */
+/*     in particular, save the diagonal elements of r in x. */
+
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	i__2 = n;
+	for (i__ = j; i__ <= i__2; ++i__) {
+	    r__[i__ + j * r_dim1] = r__[j + i__ * r_dim1];
+/* L10: */
+	}
+	x[j] = r__[j + j * r_dim1];
+	wa[j] = qtb[j];
+/* L20: */
+    }
+
+/*     eliminate the diagonal matrix d using a givens rotation. */
+
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+
+/*        prepare the row of d to be eliminated, locating the */
+/*        diagonal element using p from the qr factorization. */
+
+	l = ipvt[j];
+	if (diag[l] == 0.) {
+	    goto L90;
+	}
+	i__2 = n;
+	for (k = j; k <= i__2; ++k) {
+	    sdiag[k] = 0.;
+/* L30: */
+	}
+	sdiag[j] = diag[l];
+
+/*        the transformations to eliminate the row of d */
+/*        modify only a single element of (q transpose)*b */
+/*        beyond the first n, which is initially zero. */
+
+	qtbpj = 0.;
+	i__2 = n;
+	for (k = j; k <= i__2; ++k) {
+
+/*           determine a givens rotation which eliminates the */
+/*           appropriate element in the current row of d. */
+
+	    if (sdiag[k] == 0.) {
+		goto L70;
+	    }
+	    if ((d__1 = r__[k + k * r_dim1], ei_abs(d__1)) >= (d__2 = sdiag[k], 
+		    ei_abs(d__2))) {
+		goto L40;
+	    }
+	    cotan = r__[k + k * r_dim1] / sdiag[k];
+/* Computing 2nd power */
+	    d__1 = cotan;
+	    sin__ = p5 / ei_sqrt(p25 + p25 * (d__1 * d__1));
+	    cos__ = sin__ * cotan;
+	    goto L50;
+L40:
+	    tan__ = sdiag[k] / r__[k + k * r_dim1];
+/* Computing 2nd power */
+	    d__1 = tan__;
+	    cos__ = p5 / ei_sqrt(p25 + p25 * (d__1 * d__1));
+	    sin__ = cos__ * tan__;
+L50:
+
+/*           compute the modified diagonal element of r and */
+/*           the modified element of ((q transpose)*b,0). */
+
+	    r__[k + k * r_dim1] = cos__ * r__[k + k * r_dim1] + sin__ * sdiag[
+		    k];
+	    temp = cos__ * wa[k] + sin__ * qtbpj;
+	    qtbpj = -sin__ * wa[k] + cos__ * qtbpj;
+	    wa[k] = temp;
+
+/*           accumulate the tranformation in the row of s. */
+
+	    kp1 = k + 1;
+	    if (n < kp1) {
+		goto L70;
+	    }
+	    i__3 = n;
+	    for (i__ = kp1; i__ <= i__3; ++i__) {
+		temp = cos__ * r__[i__ + k * r_dim1] + sin__ * sdiag[i__];
+		sdiag[i__] = -sin__ * r__[i__ + k * r_dim1] + cos__ * sdiag[
+			i__];
+		r__[i__ + k * r_dim1] = temp;
+/* L60: */
+	    }
+L70:
+/* L80: */
+	    ;
+	}
+L90:
+
+/*        store the diagonal element of s and restore */
+/*        the corresponding diagonal element of r. */
+
+	sdiag[j] = r__[j + j * r_dim1];
+	r__[j + j * r_dim1] = x[j];
+/* L100: */
+    }
+
+/*     solve the triangular system for z. if the system is */
+/*     singular, then obtain a least squares solution. */
+
+    nsing = n;
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	if (sdiag[j] == 0. && nsing == n) {
+	    nsing = j - 1;
+	}
+	if (nsing < n) {
+	    wa[j] = 0.;
+	}
+/* L110: */
+    }
+    if (nsing < 1) {
+	goto L150;
+    }
+    i__1 = nsing;
+    for (k = 1; k <= i__1; ++k) {
+	j = nsing - k + 1;
+	sum = 0.;
+	jp1 = j + 1;
+	if (nsing < jp1) {
+	    goto L130;
+	}
+	i__2 = nsing;
+	for (i__ = jp1; i__ <= i__2; ++i__) {
+	    sum += r__[i__ + j * r_dim1] * wa[i__];
+/* L120: */
+	}
+L130:
+	wa[j] = (wa[j] - sum) / sdiag[j];
+/* L140: */
+    }
+L150:
+
+/*     permute the components of z back to components of x. */
+
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	l = ipvt[j];
+	x[l] = wa[j];
+/* L160: */
+    }
+    return;
+
+/*     last card of subroutine qrsolv. */
+
+} /* qrsolv_ */
+
diff --git a/unsupported/Eigen/src/NonLinear/r1mpyq.h b/unsupported/Eigen/src/NonLinear/r1mpyq.h
new file mode 100644
index 000000000..12448dabb
--- /dev/null
+++ b/unsupported/Eigen/src/NonLinear/r1mpyq.h
@@ -0,0 +1,96 @@
+
+template <typename Scalar>
+void ei_r1mpyq(int m, int n, Scalar *a, int
+	lda, const Scalar *v, const Scalar *w)
+{
+    /* System generated locals */
+    int a_dim1, a_offset, i__1, i__2;
+    Scalar d__1, d__2;
+
+    /* Local variables */
+    int i__, j, nm1, nmj;
+    Scalar cos__, sin__, temp;
+
+    /* Parameter adjustments */
+    --w;
+    --v;
+    a_dim1 = lda;
+    a_offset = 1 + a_dim1 * 1;
+    a -= a_offset;
+
+    /* Function Body */
+
+/*     apply the first set of givens rotations to a. */
+
+    nm1 = n - 1;
+    if (nm1 < 1) {
+	/* goto L50; */
+        return;
+    }
+    i__1 = nm1;
+    for (nmj = 1; nmj <= i__1; ++nmj) {
+	j = n - nmj;
+	if ((d__1 = v[j], ei_abs(d__1)) > 1.) {
+	    cos__ = 1. / v[j];
+	}
+	if ((d__1 = v[j], ei_abs(d__1)) > 1.) {
+/* Computing 2nd power */
+	    d__2 = cos__;
+	    sin__ = ei_sqrt(1. - d__2 * d__2);
+	}
+	if ((d__1 = v[j], ei_abs(d__1)) <= 1.) {
+	    sin__ = v[j];
+	}
+	if ((d__1 = v[j], ei_abs(d__1)) <= 1.) {
+/* Computing 2nd power */
+	    d__2 = sin__;
+	    cos__ = ei_sqrt(1. - d__2 * d__2);
+	}
+	i__2 = m;
+	for (i__ = 1; i__ <= i__2; ++i__) {
+	    temp = cos__ * a[i__ + j * a_dim1] - sin__ * a[i__ + n * a_dim1];
+	    a[i__ + n * a_dim1] = sin__ * a[i__ + j * a_dim1] + cos__ * a[
+		    i__ + n * a_dim1];
+	    a[i__ + j * a_dim1] = temp;
+/* L10: */
+	}
+/* L20: */
+    }
+
+/*     apply the second set of givens rotations to a. */
+
+    i__1 = nm1;
+    for (j = 1; j <= i__1; ++j) {
+	if ((d__1 = w[j], ei_abs(d__1)) > 1.) {
+	    cos__ = 1. / w[j];
+	}
+	if ((d__1 = w[j], ei_abs(d__1)) > 1.) {
+/* Computing 2nd power */
+	    d__2 = cos__;
+	    sin__ = ei_sqrt(1. - d__2 * d__2);
+	}
+	if ((d__1 = w[j], ei_abs(d__1)) <= 1.) {
+	    sin__ = w[j];
+	}
+	if ((d__1 = w[j], ei_abs(d__1)) <= 1.) {
+/* Computing 2nd power */
+	    d__2 = sin__;
+	    cos__ = ei_sqrt(1. - d__2 * d__2);
+	}
+	i__2 = m;
+	for (i__ = 1; i__ <= i__2; ++i__) {
+	    temp = cos__ * a[i__ + j * a_dim1] + sin__ * a[i__ + n * a_dim1];
+	    a[i__ + n * a_dim1] = -sin__ * a[i__ + j * a_dim1] + cos__ * a[
+		    i__ + n * a_dim1];
+	    a[i__ + j * a_dim1] = temp;
+/* L30: */
+	}
+/* L40: */
+    }
+/* L50: */
+    return;
+
+/*     last card of subroutine r1mpyq. */
+
+} /* r1mpyq_ */
+
diff --git a/unsupported/Eigen/src/NonLinear/r1updt.h b/unsupported/Eigen/src/NonLinear/r1updt.h
new file mode 100644
index 000000000..8df78321b
--- /dev/null
+++ b/unsupported/Eigen/src/NonLinear/r1updt.h
@@ -0,0 +1,200 @@
+
+template <typename Scalar>
+void ei_r1updt(int m, int n, Scalar *s, int /* ls */, const Scalar *u, Scalar *v, Scalar *w, int *sing)
+{
+    /* Initialized data */
+
+#define p5 .5
+#define p25 .25
+
+    /* System generated locals */
+    int i__1, i__2;
+    Scalar d__1, d__2;
+
+    /* Local variables */
+    int i__, j, l, jj, nm1;
+    Scalar tan__;
+    int nmj;
+    Scalar cos__, sin__, tau, temp, giant, cotan;
+
+    /* Parameter adjustments */
+    --w;
+    --u;
+    --v;
+    --s;
+
+    /* Function Body */
+
+/*     giant is the largest magnitude. */
+
+    giant = std::numeric_limits<Scalar>::max();
+
+/*     initialize the diagonal element pointer. */
+
+    jj = n * ((m << 1) - n + 1) / 2 - (m - n);
+
+/*     move the nontrivial part of the last column of s into w. */
+
+    l = jj;
+    i__1 = m;
+    for (i__ = n; i__ <= i__1; ++i__) {
+	w[i__] = s[l];
+	++l;
+/* L10: */
+    }
+
+/*     rotate the vector v into a multiple of the n-th unit vector */
+/*     in such a way that a spike is introduced into w. */
+
+    nm1 = n - 1;
+    if (nm1 < 1) {
+	goto L70;
+    }
+    i__1 = nm1;
+    for (nmj = 1; nmj <= i__1; ++nmj) {
+	j = n - nmj;
+	jj -= m - j + 1;
+	w[j] = 0.;
+	if (v[j] == 0.) {
+	    goto L50;
+	}
+
+/*        determine a givens rotation which eliminates the */
+/*        j-th element of v. */
+
+	if ((d__1 = v[n], ei_abs(d__1)) >= (d__2 = v[j], ei_abs(d__2))) {
+	    goto L20;
+	}
+	cotan = v[n] / v[j];
+/* Computing 2nd power */
+	d__1 = cotan;
+	sin__ = p5 / ei_sqrt(p25 + p25 * (d__1 * d__1));
+	cos__ = sin__ * cotan;
+	tau = 1.;
+	if (ei_abs(cos__) * giant > 1.) {
+	    tau = 1. / cos__;
+	}
+	goto L30;
+L20:
+	tan__ = v[j] / v[n];
+/* Computing 2nd power */
+	d__1 = tan__;
+	cos__ = p5 / ei_sqrt(p25 + p25 * (d__1 * d__1));
+	sin__ = cos__ * tan__;
+	tau = sin__;
+L30:
+
+/*        apply the transformation to v and store the information */
+/*        necessary to recover the givens rotation. */
+
+	v[n] = sin__ * v[j] + cos__ * v[n];
+	v[j] = tau;
+
+/*        apply the transformation to s and extend the spike in w. */
+
+	l = jj;
+	i__2 = m;
+	for (i__ = j; i__ <= i__2; ++i__) {
+	    temp = cos__ * s[l] - sin__ * w[i__];
+	    w[i__] = sin__ * s[l] + cos__ * w[i__];
+	    s[l] = temp;
+	    ++l;
+/* L40: */
+	}
+L50:
+/* L60: */
+	;
+    }
+L70:
+
+/*     add the spike from the rank 1 update to w. */
+
+    i__1 = m;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+	w[i__] += v[n] * u[i__];
+/* L80: */
+    }
+
+/*     eliminate the spike. */
+
+    *sing = false;
+    if (nm1 < 1) {
+	goto L140;
+    }
+    i__1 = nm1;
+    for (j = 1; j <= i__1; ++j) {
+	if (w[j] == 0.) {
+	    goto L120;
+	}
+
+/*        determine a givens rotation which eliminates the */
+/*        j-th element of the spike. */
+
+	if ((d__1 = s[jj], ei_abs(d__1)) >= (d__2 = w[j], ei_abs(d__2))) {
+	    goto L90;
+	}
+	cotan = s[jj] / w[j];
+/* Computing 2nd power */
+	d__1 = cotan;
+	sin__ = p5 / ei_sqrt(p25 + p25 * (d__1 * d__1));
+	cos__ = sin__ * cotan;
+	tau = 1.;
+	if (ei_abs(cos__) * giant > 1.) {
+	    tau = 1. / cos__;
+	}
+	goto L100;
+L90:
+	tan__ = w[j] / s[jj];
+/* Computing 2nd power */
+	d__1 = tan__;
+	cos__ = p5 / ei_sqrt(p25 + p25 * (d__1 * d__1));
+	sin__ = cos__ * tan__;
+	tau = sin__;
+L100:
+
+/*        apply the transformation to s and reduce the spike in w. */
+
+	l = jj;
+	i__2 = m;
+	for (i__ = j; i__ <= i__2; ++i__) {
+	    temp = cos__ * s[l] + sin__ * w[i__];
+	    w[i__] = -sin__ * s[l] + cos__ * w[i__];
+	    s[l] = temp;
+	    ++l;
+/* L110: */
+	}
+
+/*        store the information necessary to recover the */
+/*        givens rotation. */
+
+	w[j] = tau;
+L120:
+
+/*        test for zero diagonal elements in the output s. */
+
+	if (s[jj] == 0.) {
+	    *sing = true;
+	}
+	jj += m - j + 1;
+/* L130: */
+    }
+L140:
+
+/*     move w back into the last column of the output s. */
+
+    l = jj;
+    i__1 = m;
+    for (i__ = n; i__ <= i__1; ++i__) {
+	s[l] = w[i__];
+	++l;
+/* L150: */
+    }
+    if (s[jj] == 0.) {
+	*sing = true;
+    }
+    return;
+
+/*     last card of subroutine r1updt. */
+
+} /* r1updt_ */
+
diff --git a/unsupported/Eigen/src/NonLinear/rwupdt.h b/unsupported/Eigen/src/NonLinear/rwupdt.h
new file mode 100644
index 000000000..4f80e490a
--- /dev/null
+++ b/unsupported/Eigen/src/NonLinear/rwupdt.h
@@ -0,0 +1,91 @@
+
+template <typename Scalar>
+void ei_rwupdt(int n, Scalar *r__, int ldr, 
+	const Scalar *w, Scalar *b, Scalar *alpha, Scalar *cos__, 
+	Scalar *sin__)
+{
+    /* Initialized data */
+
+#define p5 .5
+#define p25 .25
+
+    /* System generated locals */
+    int r_dim1, r_offset, i__1, i__2;
+    Scalar d__1;
+
+    /* Local variables */
+    int i__, j, jm1;
+    Scalar tan__, temp, rowj, cotan;
+
+    /* Parameter adjustments */
+    --sin__;
+    --cos__;
+    --b;
+    --w;
+    r_dim1 = ldr;
+    r_offset = 1 + r_dim1 * 1;
+    r__ -= r_offset;
+
+    /* Function Body */
+
+    i__1 = n;
+    for (j = 1; j <= i__1; ++j) {
+	rowj = w[j];
+	jm1 = j - 1;
+
+/*        apply the previous transformations to */
+/*        r(i,j), i=1,2,...,j-1, and to w(j). */
+
+	if (jm1 < 1) {
+	    goto L20;
+	}
+	i__2 = jm1;
+	for (i__ = 1; i__ <= i__2; ++i__) {
+	    temp = cos__[i__] * r__[i__ + j * r_dim1] + sin__[i__] * rowj;
+	    rowj = -sin__[i__] * r__[i__ + j * r_dim1] + cos__[i__] * rowj;
+	    r__[i__ + j * r_dim1] = temp;
+/* L10: */
+	}
+L20:
+
+/*        determine a givens rotation which eliminates w(j). */
+
+	cos__[j] = 1.;
+	sin__[j] = 0.;
+	if (rowj == 0.) {
+	    goto L50;
+	}
+	if ((d__1 = r__[j + j * r_dim1], ei_abs(d__1)) >= ei_abs(rowj)) {
+	    goto L30;
+	}
+	cotan = r__[j + j * r_dim1] / rowj;
+/* Computing 2nd power */
+	d__1 = cotan;
+	sin__[j] = p5 / ei_sqrt(p25 + p25 * (d__1 * d__1));
+	cos__[j] = sin__[j] * cotan;
+	goto L40;
+L30:
+	tan__ = rowj / r__[j + j * r_dim1];
+/* Computing 2nd power */
+	d__1 = tan__;
+	cos__[j] = p5 / ei_sqrt(p25 + p25 * (d__1 * d__1));
+	sin__[j] = cos__[j] * tan__;
+L40:
+
+/*        apply the current transformation to r(j,j), b(j), and alpha. */
+
+	r__[j + j * r_dim1] = cos__[j] * r__[j + j * r_dim1] + sin__[j] * 
+		rowj;
+	temp = cos__[j] * b[j] + sin__[j] * *alpha;
+	*alpha = -sin__[j] * b[j] + cos__[j] * *alpha;
+	b[j] = temp;
+L50:
+/* L60: */
+	;
+    }
+    return;
+
+/*     last card of subroutine rwupdt. */
+
+} /* rwupdt_ */
+
diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp
index 11abd7b47..6092ad59b 100644
--- a/unsupported/test/NonLinear.cpp
+++ b/unsupported/test/NonLinear.cpp
@@ -1447,9 +1447,9 @@ void testNistMGH17(void)
           1, 100., 5000, epsilon<double>(), epsilon<double>());
 
   // check return value
-  VERIFY( 2 == info); 
-  VERIFY( 604 == nfev); 
-  VERIFY( 545 == njev); 
+  VERIFY( 1 == info); 
+  VERIFY( 599 == nfev); 
+  VERIFY( 544 == njev); 
   // check norm^2
   VERIFY_IS_APPROX(fvec.squaredNorm(), 5.4648946975E-05);
   // check x