cleaning f2c mess / trivial stuff

9 files changed, 156 insertions, 267 deletions

diff --git a/unsupported/Eigen/src/NonLinear/covar.h b/unsupported/Eigen/src/NonLinear/covar.h
index 69e649d63..7aa4976d5 100644
--- a/unsupported/Eigen/src/NonLinear/covar.h
+++ b/unsupported/Eigen/src/NonLinear/covar.h
@@ -4,10 +4,10 @@ void ei_covar(int n, Scalar *r__, int ldr,
 	const int *ipvt, Scalar tol, Scalar *wa)
 {
     /* System generated locals */
-    int r_dim1, r_offset, i__1, i__2, i__3;
+    int r_dim1, r_offset;
 
     /* Local variables */
-    int i__, j, k, l, ii, jj, km1;
+    int i, j, k, l, ii, jj, km1;
     int sing;
     Scalar temp, tolr;
 
@@ -24,8 +24,7 @@ void ei_covar(int n, Scalar *r__, int ldr,
 /*     form the inverse of r in the full upper triangle of r. */
 
     l = 0;
-    i__1 = n;
-    for (k = 1; k <= i__1; ++k) {
+    for (k = 1; k <= n; ++k) {
 	if (ei_abs(r__[k + k * r_dim1]) <= tolr) {
 	    goto L50;
 	}
@@ -34,13 +33,11 @@ void ei_covar(int n, Scalar *r__, int ldr,
 	if (km1 < 1) {
 	    goto L30;
 	}
-	i__2 = km1;
-	for (j = 1; j <= i__2; ++j) {
+	for (j = 1; j <= km1; ++j) {
 	    temp = r__[k + k * r_dim1] * r__[j + k * r_dim1];
 	    r__[j + k * r_dim1] = 0.;
-	    i__3 = j;
-	    for (i__ = 1; i__ <= i__3; ++i__) {
-		r__[i__ + k * r_dim1] -= temp * r__[i__ + j * r_dim1];
+	    for (i = 1; i <= j; ++i) {
+		r__[i + k * r_dim1] -= temp * r__[i + j * r_dim1];
 /* L10: */
 	    }
 /* L20: */
@@ -57,27 +54,23 @@ L50:
     if (l < 1) {
 	goto L110;
     }
-    i__1 = l;
-    for (k = 1; k <= i__1; ++k) {
+    for (k = 1; k <= l; ++k) {
 	km1 = k - 1;
 	if (km1 < 1) {
 	    goto L80;
 	}
-	i__2 = km1;
-	for (j = 1; j <= i__2; ++j) {
+	for (j = 1; j <= km1; ++j) {
 	    temp = r__[j + k * r_dim1];
-	    i__3 = j;
-	    for (i__ = 1; i__ <= i__3; ++i__) {
-		r__[i__ + j * r_dim1] += temp * r__[i__ + k * r_dim1];
+	    for (i = 1; i <= j; ++i) {
+		r__[i + j * r_dim1] += temp * r__[i + k * r_dim1];
 /* L60: */
 	    }
 /* L70: */
 	}
 L80:
 	temp = r__[k + k * r_dim1];
-	i__2 = k;
-	for (i__ = 1; i__ <= i__2; ++i__) {
-	    r__[i__ + k * r_dim1] = temp * r__[i__ + k * r_dim1];
+	for (i = 1; i <= k; ++i) {
+	    r__[i + k * r_dim1] = temp * r__[i + k * r_dim1];
 /* L90: */
 	}
 /* L100: */
@@ -87,21 +80,19 @@ L110:
 /*     form the full lower triangle of the covariance matrix */
 /*     in the strict lower triangle of r and in wa. */
 
-    i__1 = n;
-    for (j = 1; j <= i__1; ++j) {
+    for (j = 1; j <= n; ++j) {
 	jj = ipvt[j];
 	sing = j > l;
-	i__2 = j;
-	for (i__ = 1; i__ <= i__2; ++i__) {
+	for (i = 1; i <= j; ++i) {
 	    if (sing) {
-		r__[i__ + j * r_dim1] = 0.;
+		r__[i + j * r_dim1] = 0.;
 	    }
-	    ii = ipvt[i__];
+	    ii = ipvt[i];
 	    if (ii > jj) {
-		r__[ii + jj * r_dim1] = r__[i__ + j * r_dim1];
+		r__[ii + jj * r_dim1] = r__[i + j * r_dim1];
 	    }
 	    if (ii < jj) {
-		r__[jj + ii * r_dim1] = r__[i__ + j * r_dim1];
+		r__[jj + ii * r_dim1] = r__[i + j * r_dim1];
 	    }
 /* L120: */
 	}
@@ -111,11 +102,9 @@ L110:
 
 /*     symmetrize the covariance matrix in r. */
 
-    i__1 = n;
-    for (j = 1; j <= i__1; ++j) {
-	i__2 = j;
-	for (i__ = 1; i__ <= i__2; ++i__) {
-	    r__[i__ + j * r_dim1] = r__[j + i__ * r_dim1];
+    for (j = 1; j <= n; ++j) {
+	for (i = 1; i <= j; ++i) {
+	    r__[i + j * r_dim1] = r__[j + i * r_dim1];
 /* L140: */
 	}
 	r__[j + j * r_dim1] = wa[j];
diff --git a/unsupported/Eigen/src/NonLinear/fdjac1.h b/unsupported/Eigen/src/NonLinear/fdjac1.h
index e62bfe8d9..9d12b231f 100644
--- a/unsupported/Eigen/src/NonLinear/fdjac1.h
+++ b/unsupported/Eigen/src/NonLinear/fdjac1.h
@@ -5,11 +5,11 @@ int ei_fdjac1(minpack_func_nn fcn, void *p, int n, Scalar *x, const Scalar *
 	int mu, Scalar epsfcn, Scalar *wa1, Scalar *wa2)
 {
     /* System generated locals */
-    int fjac_dim1, fjac_offset, i__1, i__2, i__3, i__4;
+    int fjac_dim1, fjac_offset;
 
     /* Local variables */
     Scalar h__;
-    int i__, j, k;
+    int i, j, k;
     Scalar eps, temp;
     int msum;
     Scalar epsmch;
@@ -38,8 +38,7 @@ int ei_fdjac1(minpack_func_nn fcn, void *p, int n, Scalar *x, const Scalar *
 
 /*        computation of dense approximate jacobian. */
 
-    i__1 = n;
-    for (j = 1; j <= i__1; ++j) {
+    for (j = 1; j <= n; ++j) {
 	temp = x[j];
 	h__ = eps * ei_abs(temp);
 	if (h__ == 0.) {
@@ -51,9 +50,8 @@ int ei_fdjac1(minpack_func_nn fcn, void *p, int n, Scalar *x, const Scalar *
 	    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__;
+	for (i = 1; i <= n; ++i) {
+	    fjac[i + j * fjac_dim1] = (wa1[i] - fvec[i]) / h__;
 /* L10: */
 	}
 /* L20: */
@@ -65,11 +63,8 @@ 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) {
+    for (k = 1; k <= msum; ++k) {
+	for (j = k; msum< 0 ? j >= n: j <= n; j += msum) {
 	    wa2[j] = x[j];
 	    h__ = eps * ei_abs(wa2[j]);
 	    if (h__ == 0.) {
@@ -83,19 +78,16 @@ L40:
 	    /* 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) {
+	for (j = k; msum< 0 ? j >= n: j <= n; j += msum) {
 	    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__;
+	    for (i = 1; i <= n; ++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: */
 	    }
diff --git a/unsupported/Eigen/src/NonLinear/fdjac2.h b/unsupported/Eigen/src/NonLinear/fdjac2.h
index 324fbc0fe..b1c52bddf 100644
--- a/unsupported/Eigen/src/NonLinear/fdjac2.h
+++ b/unsupported/Eigen/src/NonLinear/fdjac2.h
@@ -5,11 +5,11 @@ int ei_fdjac2(minpack_func_mn fcn, void *p, int m, int n, Scalar *x,
 	Scalar epsfcn, Scalar *wa)
 {
     /* System generated locals */
-    int fjac_dim1, fjac_offset, i__1, i__2;
+    int fjac_dim1, fjac_offset;
 
     /* Local variables */
     Scalar h__;
-    int i__, j;
+    int i, j;
     Scalar eps, temp, epsmch;
     int iflag;
 
@@ -28,8 +28,7 @@ int ei_fdjac2(minpack_func_mn fcn, void *p, int m, int n, Scalar *x,
     epsmch = epsilon<Scalar>();
 
     eps = ei_sqrt((std::max(epsfcn,epsmch)));
-    i__1 = n;
-    for (j = 1; j <= i__1; ++j) {
+    for (j = 1; j <= n; ++j) {
 	temp = x[j];
 	h__ = eps * ei_abs(temp);
 	if (h__ == 0.) {
@@ -42,9 +41,8 @@ int ei_fdjac2(minpack_func_mn fcn, void *p, int m, int n, Scalar *x,
             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__;
+	for (i = 1; i <= m; ++i) {
+	    fjac[i + j * fjac_dim1] = (wa[i] - fvec[i]) / h__;
 /* L10: */
 	}
 /* L20: */
diff --git a/unsupported/Eigen/src/NonLinear/qform.h b/unsupported/Eigen/src/NonLinear/qform.h
index a9920e7b3..844c8c2f6 100644
--- a/unsupported/Eigen/src/NonLinear/qform.h
+++ b/unsupported/Eigen/src/NonLinear/qform.h
@@ -4,10 +4,10 @@ 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;
+    int q_dim1, q_offset;
 
     /* Local variables */
-    int i__, j, k, l, jm1, np1;
+    int i, j, k, l, jm1, np1;
     Scalar sum, temp;
     int minmn;
 
@@ -25,12 +25,10 @@ void ei_qform(int m, int n, Scalar *q, int
     if (minmn < 2) {
 	goto L30;
     }
-    i__1 = minmn;
-    for (j = 2; j <= i__1; ++j) {
+    for (j = 2; j <= minmn; ++j) {
 	jm1 = j - 1;
-	i__2 = jm1;
-	for (i__ = 1; i__ <= i__2; ++i__) {
-	    q[i__ + j * q_dim1] = 0.;
+	for (i = 1; i <= jm1; ++i) {
+	    q[i + j * q_dim1] = 0.;
 /* L10: */
 	}
 /* L20: */
@@ -43,11 +41,9 @@ L30:
     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.;
+    for (j = np1; j <= m; ++j) {
+	for (i = 1; i <= m; ++i) {
+	    q[i + j * q_dim1] = 0.;
 /* L40: */
 	}
 	q[j + j * q_dim1] = 1.;
@@ -57,31 +53,26 @@ L60:
 
 /*     accumulate q from its factored form. */
 
-    i__1 = minmn;
-    for (l = 1; l <= i__1; ++l) {
+    for (l = 1; l <= minmn; ++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.;
+	for (i = k; i <= m; ++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) {
+	for (j = k; j <= m; ++j) {
 	    sum = 0.;
-	    i__3 = m;
-	    for (i__ = k; i__ <= i__3; ++i__) {
-		sum += q[i__ + j * q_dim1] * wa[i__];
+	    for (i = k; i <= m; ++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__];
+	    for (i = k; i <= m; ++i) {
+		q[i + j * q_dim1] -= temp * wa[i];
 /* L90: */
 	    }
 /* L100: */
diff --git a/unsupported/Eigen/src/NonLinear/qrfac.h b/unsupported/Eigen/src/NonLinear/qrfac.h
index 245603567..3a669bbc4 100644
--- a/unsupported/Eigen/src/NonLinear/qrfac.h
+++ b/unsupported/Eigen/src/NonLinear/qrfac.h
@@ -4,16 +4,11 @@ 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;
+    int a_dim1, a_offset;
 
     /* Local variables */
-    int i__, j, k, jp1;
+    int i, j, k, jp1;
     Scalar sum;
     int kmax;
     Scalar temp;
@@ -38,8 +33,7 @@ void ei_qrfac(int m, int n, Scalar *a, int
 
 /*     compute the initial column norms and initialize several arrays. */
 
-    i__1 = n;
-    for (j = 1; j <= i__1; ++j) {
+    for (j = 1; j <= n; ++j) {
 	acnorm[j] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j * a_dim1 + 1],m).blueNorm();
 	rdiag[j] = acnorm[j];
 	wa[j] = rdiag[j];
@@ -52,8 +46,7 @@ void ei_qrfac(int m, int n, Scalar *a, int
 /*     reduce a to r with householder transformations. */
 
     minmn = std::min(m,n);
-    i__1 = minmn;
-    for (j = 1; j <= i__1; ++j) {
+    for (j = 1; j <= minmn; ++j) {
 	if (! (pivot)) {
 	    goto L40;
 	}
@@ -61,8 +54,7 @@ void ei_qrfac(int m, int n, Scalar *a, int
 /*        bring the column of largest norm into the pivot position. */
 
 	kmax = j;
-	i__2 = n;
-	for (k = j; k <= i__2; ++k) {
+	for (k = j; k <= n; ++k) {
 	    if (rdiag[k] > rdiag[kmax]) {
 		kmax = k;
 	    }
@@ -71,11 +63,10 @@ void ei_qrfac(int m, int n, Scalar *a, int
 	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;
+	for (i = 1; i <= m; ++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];
@@ -88,17 +79,15 @@ 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();
+	ajnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j + j * a_dim1],m-j+1).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;
+	for (i = j; i <= m; ++i) {
+	    a[i + j * a_dim1] /= ajnorm;
 /* L50: */
 	}
 	a[j + j * a_dim1] += 1.;
@@ -110,18 +99,15 @@ L40:
 	if (n < jp1) {
 	    goto L100;
 	}
-	i__2 = n;
-	for (k = jp1; k <= i__2; ++k) {
+	for (k = jp1; k <= n; ++k) {
 	    sum = 0.;
-	    i__3 = m;
-	    for (i__ = j; i__ <= i__3; ++i__) {
-		sum += a[i__ + j * a_dim1] * a[i__ + k * a_dim1];
+	    for (i = j; i <= m; ++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];
+	    for (i = j; i <= m; ++i) {
+		a[i + k * a_dim1] -= temp * a[i + j * a_dim1];
 /* L70: */
 	    }
 	    if (! (pivot) || rdiag[k] == 0.) {
@@ -130,16 +116,12 @@ L40:
 	    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)));
+	    rdiag[k] *= ei_sqrt((std::max(Scalar(0.), Scalar(1.)-ei_abs2(temp))));
 /* Computing 2nd power */
-	    d__1 = rdiag[k] / wa[k];
-	    if (p05 * (d__1 * d__1) > epsmch) {
+	    if (Scalar(.05) * ei_abs2(rdiag[k] / wa[k]) > epsmch) {
 		goto L80;
 	    }
-	    i__3 = m - j;
-	    rdiag[k] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[jp1 + k * a_dim1],i__3).blueNorm();
+	    rdiag[k] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[jp1 + k * a_dim1],m-j).blueNorm();
 	    wa[k] = rdiag[k];
 L80:
 /* L90: */
diff --git a/unsupported/Eigen/src/NonLinear/qrsolv.h b/unsupported/Eigen/src/NonLinear/qrsolv.h
index f4795ce95..778c66a82 100644
--- a/unsupported/Eigen/src/NonLinear/qrsolv.h
+++ b/unsupported/Eigen/src/NonLinear/qrsolv.h
@@ -4,17 +4,11 @@ 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;
+    int r_dim1, r_offset;
 
     /* Local variables */
-    int i__, j, k, l, jp1, kp1;
+    int i, j, k, l, jp1, kp1;
     Scalar tan__, cos__, sin__, sum, temp, cotan;
     int nsing;
     Scalar qtbpj;
@@ -35,11 +29,9 @@ void ei_qrsolv(int n, Scalar *r__, int ldr,
 /*     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];
+    for (j = 1; j <= n; ++j) {
+	for (i = j; i <= n; ++i) {
+	    r__[i + j * r_dim1] = r__[j + i * r_dim1];
 /* L10: */
 	}
 	x[j] = r__[j + j * r_dim1];
@@ -49,8 +41,7 @@ void ei_qrsolv(int n, Scalar *r__, int ldr,
 
 /*     eliminate the diagonal matrix d using a givens rotation. */
 
-    i__1 = n;
-    for (j = 1; j <= i__1; ++j) {
+    for (j = 1; j <= n; ++j) {
 
 /*        prepare the row of d to be eliminated, locating the */
 /*        diagonal element using p from the qr factorization. */
@@ -59,8 +50,7 @@ void ei_qrsolv(int n, Scalar *r__, int ldr,
 	if (diag[l] == 0.) {
 	    goto L90;
 	}
-	i__2 = n;
-	for (k = j; k <= i__2; ++k) {
+	for (k = j; k <= n; ++k) {
 	    sdiag[k] = 0.;
 /* L30: */
 	}
@@ -71,30 +61,24 @@ void ei_qrsolv(int n, Scalar *r__, int ldr,
 /*        beyond the first n, which is initially zero. */
 
 	qtbpj = 0.;
-	i__2 = n;
-	for (k = j; k <= i__2; ++k) {
+	for (k = j; k <= n; ++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;
-	    }
+	    if (sdiag[k] == 0.)
+            goto L70;
+	    if ( ei_abs(r__[k + k * r_dim1]) >= ei_abs(sdiag[k]))
+            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));
+	    sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan));
 	    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));
+	    cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__));
 	    sin__ = cos__ * tan__;
 L50:
 
@@ -113,12 +97,11 @@ L50:
 	    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;
+	    for (i = kp1; i <= n; ++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:
@@ -139,8 +122,7 @@ L90:
 /*     singular, then obtain a least squares solution. */
 
     nsing = n;
-    i__1 = n;
-    for (j = 1; j <= i__1; ++j) {
+    for (j = 1; j <= n; ++j) {
 	if (sdiag[j] == 0. && nsing == n) {
 	    nsing = j - 1;
 	}
@@ -152,17 +134,15 @@ L90:
     if (nsing < 1) {
 	goto L150;
     }
-    i__1 = nsing;
-    for (k = 1; k <= i__1; ++k) {
+    for (k = 1; k <= nsing; ++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__];
+	for (i = jp1; i <= nsing; ++i) {
+	    sum += r__[i + j * r_dim1] * wa[i];
 /* L120: */
 	}
 L130:
@@ -173,8 +153,7 @@ L150:
 
 /*     permute the components of z back to components of x. */
 
-    i__1 = n;
-    for (j = 1; j <= i__1; ++j) {
+    for (j = 1; j <= n; ++j) {
 	l = ipvt[j];
 	x[l] = wa[j];
 /* L160: */
diff --git a/unsupported/Eigen/src/NonLinear/r1mpyq.h b/unsupported/Eigen/src/NonLinear/r1mpyq.h
index 12448dabb..0ae7e0bbe 100644
--- a/unsupported/Eigen/src/NonLinear/r1mpyq.h
+++ b/unsupported/Eigen/src/NonLinear/r1mpyq.h
@@ -4,11 +4,10 @@ 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;
+    int a_dim1, a_offset;
 
     /* Local variables */
-    int i__, j, nm1, nmj;
+    int i, j, nm1, nmj;
     Scalar cos__, sin__, temp;
 
     /* Parameter adjustments */
@@ -27,31 +26,27 @@ void ei_r1mpyq(int m, int n, Scalar *a, int
 	/* goto L50; */
         return;
     }
-    i__1 = nm1;
-    for (nmj = 1; nmj <= i__1; ++nmj) {
+    for (nmj = 1; nmj <= nm1; ++nmj) {
 	j = n - nmj;
-	if ((d__1 = v[j], ei_abs(d__1)) > 1.) {
+	if (ei_abs(v[j]) > 1.) {
 	    cos__ = 1. / v[j];
 	}
-	if ((d__1 = v[j], ei_abs(d__1)) > 1.) {
+	if (ei_abs(v[j]) > 1.) {
 /* Computing 2nd power */
-	    d__2 = cos__;
-	    sin__ = ei_sqrt(1. - d__2 * d__2);
+	    sin__ = ei_sqrt(1. - ei_abs2(cos__));
 	}
-	if ((d__1 = v[j], ei_abs(d__1)) <= 1.) {
+	if (ei_abs(v[j]) <= 1.) {
 	    sin__ = v[j];
 	}
-	if ((d__1 = v[j], ei_abs(d__1)) <= 1.) {
+	if (ei_abs(v[j]) <= 1.) {
 /* Computing 2nd power */
-	    d__2 = sin__;
-	    cos__ = ei_sqrt(1. - d__2 * d__2);
+	    cos__ = ei_sqrt(1. - ei_abs2(sin__));
 	}
-	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;
+	for (i = 1; i <= m; ++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: */
@@ -59,30 +54,26 @@ void ei_r1mpyq(int m, int n, Scalar *a, int
 
 /*     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.) {
+    for (j = 1; j <= nm1; ++j) {
+	if (ei_abs(w[j]) > 1.) {
 	    cos__ = 1. / w[j];
 	}
-	if ((d__1 = w[j], ei_abs(d__1)) > 1.) {
+	if (ei_abs(w[j]) > 1.) {
 /* Computing 2nd power */
-	    d__2 = cos__;
-	    sin__ = ei_sqrt(1. - d__2 * d__2);
+	    sin__ = ei_sqrt(1. - ei_abs2(cos__));
 	}
-	if ((d__1 = w[j], ei_abs(d__1)) <= 1.) {
+	if (ei_abs(w[j]) <= 1.) {
 	    sin__ = w[j];
 	}
-	if ((d__1 = w[j], ei_abs(d__1)) <= 1.) {
+	if (ei_abs(w[j]) <= 1.) {
 /* Computing 2nd power */
-	    d__2 = sin__;
-	    cos__ = ei_sqrt(1. - d__2 * d__2);
+	    cos__ = ei_sqrt(1. - ei_abs2(sin__));
 	}
-	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;
+	for (i = 1; i <= m; ++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: */
diff --git a/unsupported/Eigen/src/NonLinear/r1updt.h b/unsupported/Eigen/src/NonLinear/r1updt.h
index 8df78321b..947a0e1e3 100644
--- a/unsupported/Eigen/src/NonLinear/r1updt.h
+++ b/unsupported/Eigen/src/NonLinear/r1updt.h
@@ -2,17 +2,8 @@
 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;
+    int i, j, l, jj, nm1;
     Scalar tan__;
     int nmj;
     Scalar cos__, sin__, tau, temp, giant, cotan;
@@ -36,9 +27,8 @@ void ei_r1updt(int m, int n, Scalar *s, int /* ls */, const Scalar *u, Scalar *v
 /*     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];
+    for (i = n; i <= m; ++i) {
+	w[i] = s[l];
 	++l;
 /* L10: */
     }
@@ -50,8 +40,7 @@ void ei_r1updt(int m, int n, Scalar *s, int /* ls */, const Scalar *u, Scalar *v
     if (nm1 < 1) {
 	goto L70;
     }
-    i__1 = nm1;
-    for (nmj = 1; nmj <= i__1; ++nmj) {
+    for (nmj = 1; nmj <= nm1; ++nmj) {
 	j = n - nmj;
 	jj -= m - j + 1;
 	w[j] = 0.;
@@ -62,13 +51,11 @@ void ei_r1updt(int m, int n, Scalar *s, int /* ls */, const Scalar *u, Scalar *v
 /*        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))) {
+	if (ei_abs(v[n]) >= ei_abs(v[j]))
 	    goto L20;
-	}
 	cotan = v[n] / v[j];
 /* Computing 2nd power */
-	d__1 = cotan;
-	sin__ = p5 / ei_sqrt(p25 + p25 * (d__1 * d__1));
+	sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan));
 	cos__ = sin__ * cotan;
 	tau = 1.;
 	if (ei_abs(cos__) * giant > 1.) {
@@ -78,8 +65,7 @@ void ei_r1updt(int m, int n, Scalar *s, int /* ls */, const Scalar *u, Scalar *v
 L20:
 	tan__ = v[j] / v[n];
 /* Computing 2nd power */
-	d__1 = tan__;
-	cos__ = p5 / ei_sqrt(p25 + p25 * (d__1 * d__1));
+	cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__));
 	sin__ = cos__ * tan__;
 	tau = sin__;
 L30:
@@ -93,10 +79,9 @@ L30:
 /*        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__];
+	for (i = j; i <= m; ++i) {
+	    temp = cos__ * s[l] - sin__ * w[i];
+	    w[i] = sin__ * s[l] + cos__ * w[i];
 	    s[l] = temp;
 	    ++l;
 /* L40: */
@@ -109,9 +94,8 @@ 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__];
+    for (i = 1; i <= m; ++i) {
+	w[i] += v[n] * u[i];
 /* L80: */
     }
 
@@ -121,8 +105,7 @@ L70:
     if (nm1 < 1) {
 	goto L140;
     }
-    i__1 = nm1;
-    for (j = 1; j <= i__1; ++j) {
+    for (j = 1; j <= nm1; ++j) {
 	if (w[j] == 0.) {
 	    goto L120;
 	}
@@ -130,13 +113,11 @@ L70:
 /*        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))) {
+	if (ei_abs(s[jj]) >= ei_abs(w[j]))
 	    goto L90;
-	}
 	cotan = s[jj] / w[j];
 /* Computing 2nd power */
-	d__1 = cotan;
-	sin__ = p5 / ei_sqrt(p25 + p25 * (d__1 * d__1));
+	sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan));
 	cos__ = sin__ * cotan;
 	tau = 1.;
 	if (ei_abs(cos__) * giant > 1.) {
@@ -146,8 +127,7 @@ L70:
 L90:
 	tan__ = w[j] / s[jj];
 /* Computing 2nd power */
-	d__1 = tan__;
-	cos__ = p5 / ei_sqrt(p25 + p25 * (d__1 * d__1));
+	cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__));
 	sin__ = cos__ * tan__;
 	tau = sin__;
 L100:
@@ -155,10 +135,9 @@ 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__];
+	for (i = j; i <= m; ++i) {
+	    temp = cos__ * s[l] + sin__ * w[i];
+	    w[i] = -sin__ * s[l] + cos__ * w[i];
 	    s[l] = temp;
 	    ++l;
 /* L110: */
@@ -183,9 +162,8 @@ 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__];
+    for (i = n; i <= m; ++i) {
+	s[l] = w[i];
 	++l;
 /* L150: */
     }
diff --git a/unsupported/Eigen/src/NonLinear/rwupdt.h b/unsupported/Eigen/src/NonLinear/rwupdt.h
index 4f80e490a..6be292f6a 100644
--- a/unsupported/Eigen/src/NonLinear/rwupdt.h
+++ b/unsupported/Eigen/src/NonLinear/rwupdt.h
@@ -4,17 +4,11 @@ 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;
+    int r_dim1, r_offset;
 
     /* Local variables */
-    int i__, j, jm1;
+    int i, j, jm1;
     Scalar tan__, temp, rowj, cotan;
 
     /* Parameter adjustments */
@@ -28,8 +22,7 @@ void ei_rwupdt(int n, Scalar *r__, int ldr,
 
     /* Function Body */
 
-    i__1 = n;
-    for (j = 1; j <= i__1; ++j) {
+    for (j = 1; j <= n; ++j) {
 	rowj = w[j];
 	jm1 = j - 1;
 
@@ -39,11 +32,10 @@ void ei_rwupdt(int n, Scalar *r__, int ldr,
 	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;
+	for (i = 1; i <= jm1; ++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:
@@ -55,20 +47,17 @@ L20:
 	if (rowj == 0.) {
 	    goto L50;
 	}
-	if ((d__1 = r__[j + j * r_dim1], ei_abs(d__1)) >= ei_abs(rowj)) {
+	if (ei_abs(r__[j + j * r_dim1]) >= 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));
+	sin__[j] = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan));
 	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));
+	cos__[j] = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__));
 	sin__[j] = cos__[j] * tan__;
 L40: