Now that the main algorithms are imported into eigen, we import subroutines

used by those algorithms (aka "second level").

This is a row import : we copy/paste the files from cminpack and make
very few changes :
* template<Scalar> them (replace double)
* dpmpar() replaced by c++ standard equivalent
* abs/fabs/sqrt/min/max replaced by ei_* or std::*
* use eigen norms instead of enorm()

Important Notes:
* The use of stableNorm() was not enough in some cases, but using
  blueNorm() instead fixed the problems (some tests gave bad results,
  either in number of iterations or precision of the results)
* As a whole, the only test that changed is testNistMGH17() : it now takes
  some few steps less to get the same result. So this is a small improvement.

After this commit, the only remaining dependency from the cminpack
static library is 'covar', only used from the tests.

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"

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_ */
+

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_ */
+

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_ */
+

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);
+    ei_r1updt<Scalar>(n, n, R.data(), lr, wa1.data(), wa2.data(), wa3.data(), &sing);
-    r1mpyq(n, n, fjac.data(), ldfjac, wa2.data(), wa3.data());
+    ei_r1mpyq<Scalar>(n, n, fjac.data(), ldfjac, wa2.data(), wa3.data());
-    r1mpyq(1, n, qtf.data(), 1, wa2.data(), wa3.data());
+    ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data());
 
     /*           end of the inner loop. */
 

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);
+    ei_r1updt<Scalar>(n, n, R.data(), lr, wa1.data(), wa2.data(), wa3.data(), &sing);
-    r1mpyq(n, n, fjac.data(), ldfjac, wa2.data(), wa3.data());
+    ei_r1mpyq<Scalar>(n, n, fjac.data(), ldfjac, wa2.data(), wa3.data());
-    r1mpyq(1, n, qtf.data(), 1, wa2.data(), wa3.data());
+    ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data());
 
     /*           end of the inner loop. */
 

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. */

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;
 

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_ */
+

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;
 

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_ */
+

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_ */
+

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_ */
+

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_ */
+

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_ */
+

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_ */
+

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( 1 == info); 
-  VERIFY( 604 == nfev); 
+  VERIFY( 599 == nfev); 
-  VERIFY( 545 == njev); 
+  VERIFY( 544 == njev); 
   // check norm^2
   VERIFY_IS_APPROX(fvec.squaredNorm(), 5.4648946975E-05);
   // check x