diff --git a/unsupported/Eigen/src/NonLinear/hybrd.h b/unsupported/Eigen/src/NonLinear/hybrd.h
index e66a0f918..5c56c0576 100644
--- a/unsupported/Eigen/src/NonLinear/hybrd.h
+++ b/unsupported/Eigen/src/NonLinear/hybrd.h
@@ -1,10 +1,10 @@
 
 template<typename Scalar>
 int hybrd_template(minpack_func_nn fcn, void *p, int n, Scalar *x, Scalar *
-	fvec, Scalar xtol, int maxfev, int ml, int mu, 
-	Scalar epsfcn, Scalar *diag, int mode, Scalar factor, int nprint, int *nfev, Scalar *
-	fjac, int ldfjac, Scalar *r__, int lr, Scalar *qtf, 
-	Scalar *wa1, Scalar *wa2, Scalar *wa3, Scalar *wa4)
+        fvec, Scalar xtol, int maxfev, int ml, int mu, 
+        Scalar epsfcn, Scalar *diag, int mode, Scalar factor, int nprint, int *nfev, Scalar *
+        fjac, int ldfjac, Scalar *r__, int lr, Scalar *qtf, 
+        Scalar *wa1, Scalar *wa2, Scalar *wa3, Scalar *wa4)
 {
     /* Initialized data */
 
@@ -50,40 +50,40 @@ int hybrd_template(minpack_func_nn fcn, void *p, int n, Scalar *x, Scalar *
     iflag = 0;
     *nfev = 0;
 
-/*     check the input parameters for errors. */
+    /*     check the input parameters for errors. */
 
     if (n <= 0 || xtol < 0. || maxfev <= 0 || ml < 0 || mu < 0 ||
-	    factor <= 0. || ldfjac < n || lr < n * (n + 1) / 2) {
-	goto L300;
+            factor <= 0. || ldfjac < n || lr < n * (n + 1) / 2) {
+        goto L300;
     }
     if (mode != 2) {
-	goto L20;
+        goto L20;
     }
     for (j = 1; j <= n; ++j) {
-	if (diag[j] <= 0.) {
-	    goto L300;
-	}
-/* L10: */
+        if (diag[j] <= 0.) {
+            goto L300;
+        }
+        /* L10: */
     }
 L20:
 
-/*     evaluate the function at the starting point */
-/*     and calculate its norm. */
+    /*     evaluate the function at the starting point */
+    /*     and calculate its norm. */
 
     iflag = (*fcn)(p, n, &x[1], &fvec[1], 1);
     *nfev = 1;
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
     fnorm = ei_enorm<Scalar>(n, &fvec[1]);
 
-/*     determine the number of calls to fcn needed to compute */
-/*     the jacobian matrix. */
+    /*     determine the number of calls to fcn needed to compute */
+    /*     the jacobian matrix. */
 
-/* Computing MIN */
+    /* Computing MIN */
     msum = min(ml + mu + 1, n);
 
-/*     initialize iteration counter and monitors. */
+    /*     initialize iteration counter and monitors. */
 
     iter = 1;
     ncsuc = 0;
@@ -91,211 +91,211 @@ L20:
     nslow1 = 0;
     nslow2 = 0;
 
-/*     beginning of the outer loop. */
+    /*     beginning of the outer loop. */
 
 L30:
     jeval = TRUE_;
 
-/*        calculate the jacobian matrix. */
+    /*        calculate the jacobian matrix. */
 
     iflag = fdjac1(fcn, p, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac,
-	     ml, mu, epsfcn, &wa1[1], &wa2[1]);
+            ml, mu, epsfcn, &wa1[1], &wa2[1]);
     *nfev += msum;
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
 
-/*        compute the qr factorization of the jacobian. */
+    /*        compute the qr factorization of the jacobian. */
 
     qrfac(n, n, &fjac[fjac_offset], ldfjac, FALSE_, iwa, 1, &wa1[1], &
-	    wa2[1], &wa3[1]);
+            wa2[1], &wa3[1]);
 
-/*        on the first iteration and if mode is 1, scale according */
-/*        to the norms of the columns of the initial jacobian. */
+    /*        on the first iteration and if mode is 1, scale according */
+    /*        to the norms of the columns of the initial jacobian. */
 
     if (iter != 1) {
-	goto L70;
+        goto L70;
     }
     if (mode == 2) {
-	goto L50;
+        goto L50;
     }
     for (j = 1; j <= n; ++j) {
-	diag[j] = wa2[j];
-	if (wa2[j] == 0.) {
-	    diag[j] = 1.;
-	}
-/* L40: */
+        diag[j] = wa2[j];
+        if (wa2[j] == 0.) {
+            diag[j] = 1.;
+        }
+        /* L40: */
     }
 L50:
 
-/*        on the first iteration, calculate the norm of the scaled x */
-/*        and initialize the step bound delta. */
+    /*        on the first iteration, calculate the norm of the scaled x */
+    /*        and initialize the step bound delta. */
 
     for (j = 1; j <= n; ++j) {
-	wa3[j] = diag[j] * x[j];
-/* L60: */
+        wa3[j] = diag[j] * x[j];
+        /* L60: */
     }
     xnorm = ei_enorm<Scalar>(n, &wa3[1]);
     delta = factor * xnorm;
     if (delta == 0.) {
-	delta = factor;
+        delta = factor;
     }
 L70:
 
-/*        form (q transpose)*fvec and store in qtf. */
+    /*        form (q transpose)*fvec and store in qtf. */
 
     for (i__ = 1; i__ <= n; ++i__) {
-	qtf[i__] = fvec[i__];
-/* L80: */
+        qtf[i__] = fvec[i__];
+        /* L80: */
     }
     for (j = 1; j <= n; ++j) {
-	if (fjac[j + j * fjac_dim1] == 0.) {
-	    goto L110;
-	}
-	sum = 0.;
-	for (i__ = j; i__ <= n; ++i__) {
-	    sum += fjac[i__ + j * fjac_dim1] * qtf[i__];
-/* L90: */
-	}
-	temp = -sum / fjac[j + j * fjac_dim1];
-	for (i__ = j; i__ <= n; ++i__) {
-	    qtf[i__] += fjac[i__ + j * fjac_dim1] * temp;
-/* L100: */
-	}
+        if (fjac[j + j * fjac_dim1] == 0.) {
+            goto L110;
+        }
+        sum = 0.;
+        for (i__ = j; i__ <= n; ++i__) {
+            sum += fjac[i__ + j * fjac_dim1] * qtf[i__];
+            /* L90: */
+        }
+        temp = -sum / fjac[j + j * fjac_dim1];
+        for (i__ = j; i__ <= n; ++i__) {
+            qtf[i__] += fjac[i__ + j * fjac_dim1] * temp;
+            /* L100: */
+        }
 L110:
-/* L120: */
-	;
+        /* L120: */
+        ;
     }
 
-/*        copy the triangular factor of the qr factorization into r. */
+    /*        copy the triangular factor of the qr factorization into r. */
 
     sing = FALSE_;
     for (j = 1; j <= n; ++j) {
-	l = j;
-	jm1 = j - 1;
-	if (jm1 < 1) {
-	    goto L140;
-	}
-	for (i__ = 1; i__ <= jm1; ++i__) {
-	    r__[l] = fjac[i__ + j * fjac_dim1];
-	    l = l + n - i__;
-/* L130: */
-	}
+        l = j;
+        jm1 = j - 1;
+        if (jm1 < 1) {
+            goto L140;
+        }
+        for (i__ = 1; i__ <= jm1; ++i__) {
+            r__[l] = fjac[i__ + j * fjac_dim1];
+            l = l + n - i__;
+            /* L130: */
+        }
 L140:
-	r__[l] = wa1[j];
-	if (wa1[j] == 0.) {
-	    sing = TRUE_;
-	}
-/* L150: */
+        r__[l] = wa1[j];
+        if (wa1[j] == 0.) {
+            sing = TRUE_;
+        }
+        /* L150: */
     }
 
-/*        accumulate the orthogonal factor in fjac. */
+    /*        accumulate the orthogonal factor in fjac. */
 
     qform(n, n, &fjac[fjac_offset], ldfjac, &wa1[1]);
 
-/*        rescale if necessary. */
+    /*        rescale if necessary. */
 
     if (mode == 2) {
-	goto L170;
+        goto L170;
     }
     for (j = 1; j <= n; ++j) {
-/* Computing MAX */
-	d__1 = diag[j], d__2 = wa2[j];
-	diag[j] = max(d__1,d__2);
-/* L160: */
+        /* Computing MAX */
+        d__1 = diag[j], d__2 = wa2[j];
+        diag[j] = max(d__1,d__2);
+        /* L160: */
     }
 L170:
 
-/*        beginning of the inner loop. */
+    /*        beginning of the inner loop. */
 
 L180:
 
-/*           if requested, call fcn to enable printing of iterates. */
+    /*           if requested, call fcn to enable printing of iterates. */
 
     if (nprint <= 0) {
-	goto L190;
+        goto L190;
     }
     iflag = 0;
     if ((iter - 1) % nprint == 0) {
-	iflag = (*fcn)(p, n, &x[1], &fvec[1], 0);
+        iflag = (*fcn)(p, n, &x[1], &fvec[1], 0);
     }
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
 L190:
 
-/*           determine the direction p. */
+    /*           determine the direction p. */
 
     dogleg(n, &r__[1], lr, &diag[1], &qtf[1], delta, &wa1[1], &wa2[1], &wa3[
-	    1]);
+            1]);
 
-/*           store the direction p and x + p. calculate the norm of p. */
+    /*           store the direction p and x + p. calculate the norm of p. */
 
     for (j = 1; j <= n; ++j) {
-	wa1[j] = -wa1[j];
-	wa2[j] = x[j] + wa1[j];
-	wa3[j] = diag[j] * wa1[j];
-/* L200: */
+        wa1[j] = -wa1[j];
+        wa2[j] = x[j] + wa1[j];
+        wa3[j] = diag[j] * wa1[j];
+        /* L200: */
     }
     pnorm = ei_enorm<Scalar>(n, &wa3[1]);
 
-/*           on the first iteration, adjust the initial step bound. */
+    /*           on the first iteration, adjust the initial step bound. */
 
     if (iter == 1) {
-	delta = min(delta,pnorm);
+        delta = min(delta,pnorm);
     }
 
-/*           evaluate the function at x + p and calculate its norm. */
+    /*           evaluate the function at x + p and calculate its norm. */
 
     iflag = (*fcn)(p, n, &wa2[1], &wa4[1], 1);
     ++(*nfev);
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
     fnorm1 = ei_enorm<Scalar>(n, &wa4[1]);
 
-/*           compute the scaled actual reduction. */
+    /*           compute the scaled actual reduction. */
 
     actred = -1.;
     if (fnorm1 < fnorm) {
-/* Computing 2nd power */
-	d__1 = fnorm1 / fnorm;
-	actred = 1. - d__1 * d__1;
+        /* Computing 2nd power */
+        d__1 = fnorm1 / fnorm;
+        actred = 1. - d__1 * d__1;
     }
 
-/*           compute the scaled predicted reduction. */
+    /*           compute the scaled predicted reduction. */
 
     l = 1;
     for (i__ = 1; i__ <= n; ++i__) {
-	sum = 0.;
-	for (j = i__; j <= n; ++j) {
-	    sum += r__[l] * wa1[j];
-	    ++l;
-/* L210: */
-	}
-	wa3[i__] = qtf[i__] + sum;
-/* L220: */
+        sum = 0.;
+        for (j = i__; j <= n; ++j) {
+            sum += r__[l] * wa1[j];
+            ++l;
+            /* L210: */
+        }
+        wa3[i__] = qtf[i__] + sum;
+        /* L220: */
     }
     temp = ei_enorm<Scalar>(n, &wa3[1]);
     prered = 0.;
     if (temp < fnorm) {
-/* Computing 2nd power */
-	d__1 = temp / fnorm;
-	prered = 1. - d__1 * d__1;
+        /* Computing 2nd power */
+        d__1 = temp / fnorm;
+        prered = 1. - d__1 * d__1;
     }
 
-/*           compute the ratio of the actual to the predicted */
-/*           reduction. */
+    /*           compute the ratio of the actual to the predicted */
+    /*           reduction. */
 
     ratio = 0.;
     if (prered > 0.) {
-	ratio = actred / prered;
+        ratio = actred / prered;
     }
 
-/*           update the step bound. */
+    /*           update the step bound. */
 
     if (ratio >= p1) {
-	goto L230;
+        goto L230;
     }
     ncsuc = 0;
     ++ncfail;
@@ -305,128 +305,128 @@ L230:
     ncfail = 0;
     ++ncsuc;
     if (ratio >= p5 || ncsuc > 1) {
-/* Computing MAX */
-	d__1 = delta, d__2 = pnorm / p5;
-	delta = max(d__1,d__2);
+        /* Computing MAX */
+        d__1 = delta, d__2 = pnorm / p5;
+        delta = max(d__1,d__2);
     }
     if (fabs(ratio - 1.) <= p1) {
-	delta = pnorm / p5;
+        delta = pnorm / p5;
     }
 L240:
 
-/*           test for successful iteration. */
+    /*           test for successful iteration. */
 
     if (ratio < p0001) {
-	goto L260;
+        goto L260;
     }
 
-/*           successful iteration. update x, fvec, and their norms. */
+    /*           successful iteration. update x, fvec, and their norms. */
 
     for (j = 1; j <= n; ++j) {
-	x[j] = wa2[j];
-	wa2[j] = diag[j] * x[j];
-	fvec[j] = wa4[j];
-/* L250: */
+        x[j] = wa2[j];
+        wa2[j] = diag[j] * x[j];
+        fvec[j] = wa4[j];
+        /* L250: */
     }
     xnorm = ei_enorm<Scalar>(n, &wa2[1]);
     fnorm = fnorm1;
     ++iter;
 L260:
 
-/*           determine the progress of the iteration. */
+    /*           determine the progress of the iteration. */
 
     ++nslow1;
     if (actred >= p001) {
-	nslow1 = 0;
+        nslow1 = 0;
     }
     if (jeval) {
-	++nslow2;
+        ++nslow2;
     }
     if (actred >= p1) {
-	nslow2 = 0;
+        nslow2 = 0;
     }
 
-/*           test for convergence. */
+    /*           test for convergence. */
 
     if (delta <= xtol * xnorm || fnorm == 0.) {
-	info = 1;
+        info = 1;
     }
     if (info != 0) {
-	goto L300;
+        goto L300;
     }
 
-/*           tests for termination and stringent tolerances. */
+    /*           tests for termination and stringent tolerances. */
 
     if (*nfev >= maxfev) {
-	info = 2;
+        info = 2;
     }
-/* Computing MAX */
+    /* Computing MAX */
     d__1 = p1 * delta;
     if (p1 * max(d__1,pnorm) <= epsilon<Scalar>() * xnorm) {
-	info = 3;
+        info = 3;
     }
     if (nslow2 == 5) {
-	info = 4;
+        info = 4;
     }
     if (nslow1 == 10) {
-	info = 5;
+        info = 5;
     }
     if (info != 0) {
-	goto L300;
+        goto L300;
     }
 
-/*           criterion for recalculating jacobian approximation */
-/*           by forward differences. */
+    /*           criterion for recalculating jacobian approximation */
+    /*           by forward differences. */
 
     if (ncfail == 2) {
-	goto L290;
+        goto L290;
     }
 
-/*           calculate the rank one modification to the jacobian */
-/*           and update qtf if necessary. */
+    /*           calculate the rank one modification to the jacobian */
+    /*           and update qtf if necessary. */
 
     for (j = 1; j <= n; ++j) {
-	sum = 0.;
-	for (i__ = 1; i__ <= n; ++i__) {
-	    sum += fjac[i__ + j * fjac_dim1] * wa4[i__];
-/* L270: */
-	}
-	wa2[j] = (sum - wa3[j]) / pnorm;
-	wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
-	if (ratio >= p0001) {
-	    qtf[j] = sum;
-	}
-/* L280: */
+        sum = 0.;
+        for (i__ = 1; i__ <= n; ++i__) {
+            sum += fjac[i__ + j * fjac_dim1] * wa4[i__];
+            /* L270: */
+        }
+        wa2[j] = (sum - wa3[j]) / pnorm;
+        wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
+        if (ratio >= p0001) {
+            qtf[j] = sum;
+        }
+        /* L280: */
     }
 
-/*           compute the qr factorization of the updated jacobian. */
+    /*           compute the qr factorization of the updated jacobian. */
 
     r1updt(n, n, &r__[1], lr, &wa1[1], &wa2[1], &wa3[1], &sing);
     r1mpyq(n, n, &fjac[fjac_offset], ldfjac, &wa2[1], &wa3[1]);
     r1mpyq(1, n, &qtf[1], 1, &wa2[1], &wa3[1]);
 
-/*           end of the inner loop. */
+    /*           end of the inner loop. */
 
     jeval = FALSE_;
     goto L180;
 L290:
 
-/*        end of the outer loop. */
+    /*        end of the outer loop. */
 
     goto L30;
 L300:
 
-/*     termination, either normal or user imposed. */
+    /*     termination, either normal or user imposed. */
 
     if (iflag < 0) {
-	info = iflag;
+        info = iflag;
     }
     if (nprint > 0) {
-	(*fcn)(p, n, &x[1], &fvec[1], 0);
+        (*fcn)(p, n, &x[1], &fvec[1], 0);
     }
     return info;
 
-/*     last card of subroutine hybrd. */
+    /*     last card of subroutine hybrd. */
 
 } /* hybrd_ */
 
diff --git a/unsupported/Eigen/src/NonLinear/hybrj.h b/unsupported/Eigen/src/NonLinear/hybrj.h
index c82e196a4..ee207e75d 100644
--- a/unsupported/Eigen/src/NonLinear/hybrj.h
+++ b/unsupported/Eigen/src/NonLinear/hybrj.h
@@ -1,11 +1,11 @@
 
 template<typename Scalar>
 int hybrj_template(minpack_funcder_nn fcn, void *p, int n, Scalar *x, Scalar *
-	fvec, Scalar *fjac, int ldfjac, Scalar xtol, int
-	maxfev, Scalar *diag, int mode, Scalar factor, int
-	nprint, int *nfev, int *njev, Scalar *r__, 
-	int lr, Scalar *qtf, Scalar *wa1, Scalar *wa2, 
-	Scalar *wa3, Scalar *wa4)
+        fvec, Scalar *fjac, int ldfjac, Scalar xtol, int
+        maxfev, Scalar *diag, int mode, Scalar factor, int
+        nprint, int *nfev, int *njev, Scalar *r__, 
+        int lr, Scalar *qtf, Scalar *wa1, Scalar *wa2, 
+        Scalar *wa3, Scalar *wa4)
 {
     /* Initialized data */
 
@@ -52,34 +52,34 @@ int hybrj_template(minpack_funcder_nn fcn, void *p, int n, Scalar *x, Scalar *
     *nfev = 0;
     *njev = 0;
 
-/*     check the input parameters for errors. */
+    /*     check the input parameters for errors. */
 
     if (n <= 0 || ldfjac < n || xtol < 0. || maxfev <= 0 || factor <= 
-	    0. || lr < n * (n + 1) / 2) {
-	goto L300;
+            0. || lr < n * (n + 1) / 2) {
+        goto L300;
     }
     if (mode != 2) {
-	goto L20;
+        goto L20;
     }
     for (j = 1; j <= n; ++j) {
-	if (diag[j] <= 0.) {
-	    goto L300;
-	}
-/* L10: */
+        if (diag[j] <= 0.) {
+            goto L300;
+        }
+        /* L10: */
     }
 L20:
 
-/*     evaluate the function at the starting point */
-/*     and calculate its norm. */
+    /*     evaluate the function at the starting point */
+    /*     and calculate its norm. */
 
     iflag = (*fcn)(p, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 1);
     *nfev = 1;
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
     fnorm = ei_enorm<Scalar>(n, &fvec[1]);
 
-/*     initialize iteration counter and monitors. */
+    /*     initialize iteration counter and monitors. */
 
     iter = 1;
     ncsuc = 0;
@@ -87,210 +87,210 @@ L20:
     nslow1 = 0;
     nslow2 = 0;
 
-/*     beginning of the outer loop. */
+    /*     beginning of the outer loop. */
 
 L30:
     jeval = TRUE_;
 
-/*        calculate the jacobian matrix. */
+    /*        calculate the jacobian matrix. */
 
     iflag = (*fcn)(p, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 2);
     ++(*njev);
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
 
-/*        compute the qr factorization of the jacobian. */
+    /*        compute the qr factorization of the jacobian. */
 
     qrfac(n, n, &fjac[fjac_offset], ldfjac, FALSE_, iwa, 1, &wa1[1], &
-	    wa2[1], &wa3[1]);
+            wa2[1], &wa3[1]);
 
-/*        on the first iteration and if mode is 1, scale according */
-/*        to the norms of the columns of the initial jacobian. */
+    /*        on the first iteration and if mode is 1, scale according */
+    /*        to the norms of the columns of the initial jacobian. */
 
     if (iter != 1) {
-	goto L70;
+        goto L70;
     }
     if (mode == 2) {
-	goto L50;
+        goto L50;
     }
     for (j = 1; j <= n; ++j) {
-	diag[j] = wa2[j];
-	if (wa2[j] == 0.) {
-	    diag[j] = 1.;
-	}
-/* L40: */
+        diag[j] = wa2[j];
+        if (wa2[j] == 0.) {
+            diag[j] = 1.;
+        }
+        /* L40: */
     }
 L50:
 
-/*        on the first iteration, calculate the norm of the scaled x */
-/*        and initialize the step bound delta. */
+    /*        on the first iteration, calculate the norm of the scaled x */
+    /*        and initialize the step bound delta. */
 
     for (j = 1; j <= n; ++j) {
-	wa3[j] = diag[j] * x[j];
-/* L60: */
+        wa3[j] = diag[j] * x[j];
+        /* L60: */
     }
     xnorm = ei_enorm<Scalar>(n, &wa3[1]);
     delta = factor * xnorm;
     if (delta == 0.) {
-	delta = factor;
+        delta = factor;
     }
 L70:
 
-/*        form (q transpose)*fvec and store in qtf. */
+    /*        form (q transpose)*fvec and store in qtf. */
 
     for (i__ = 1; i__ <= n; ++i__) {
-	qtf[i__] = fvec[i__];
-/* L80: */
+        qtf[i__] = fvec[i__];
+        /* L80: */
     }
     for (j = 1; j <= n; ++j) {
-	if (fjac[j + j * fjac_dim1] == 0.) {
-	    goto L110;
-	}
-	sum = 0.;
-	for (i__ = j; i__ <= n; ++i__) {
-	    sum += fjac[i__ + j * fjac_dim1] * qtf[i__];
-/* L90: */
-	}
-	temp = -sum / fjac[j + j * fjac_dim1];
-	for (i__ = j; i__ <= n; ++i__) {
-	    qtf[i__] += fjac[i__ + j * fjac_dim1] * temp;
-/* L100: */
-	}
+        if (fjac[j + j * fjac_dim1] == 0.) {
+            goto L110;
+        }
+        sum = 0.;
+        for (i__ = j; i__ <= n; ++i__) {
+            sum += fjac[i__ + j * fjac_dim1] * qtf[i__];
+            /* L90: */
+        }
+        temp = -sum / fjac[j + j * fjac_dim1];
+        for (i__ = j; i__ <= n; ++i__) {
+            qtf[i__] += fjac[i__ + j * fjac_dim1] * temp;
+            /* L100: */
+        }
 L110:
-/* L120: */
-	;
+        /* L120: */
+        ;
     }
 
-/*        copy the triangular factor of the qr factorization into r. */
+    /*        copy the triangular factor of the qr factorization into r. */
 
     sing = FALSE_;
     for (j = 1; j <= n; ++j) {
-	l = j;
-	jm1 = j - 1;
-	if (jm1 < 1) {
-	    goto L140;
-	}
-	for (i__ = 1; i__ <= jm1; ++i__) {
-	    r__[l] = fjac[i__ + j * fjac_dim1];
-	    l = l + n - i__;
-/* L130: */
-	}
+        l = j;
+        jm1 = j - 1;
+        if (jm1 < 1) {
+            goto L140;
+        }
+        for (i__ = 1; i__ <= jm1; ++i__) {
+            r__[l] = fjac[i__ + j * fjac_dim1];
+            l = l + n - i__;
+            /* L130: */
+        }
 L140:
-	r__[l] = wa1[j];
-	if (wa1[j] == 0.) {
-	    sing = TRUE_;
-	}
-/* L150: */
+        r__[l] = wa1[j];
+        if (wa1[j] == 0.) {
+            sing = TRUE_;
+        }
+        /* L150: */
     }
 
-/*        accumulate the orthogonal factor in fjac. */
+    /*        accumulate the orthogonal factor in fjac. */
 
     qform(n, n, &fjac[fjac_offset], ldfjac, &wa1[1]);
 
-/*        rescale if necessary. */
+    /*        rescale if necessary. */
 
     if (mode == 2) {
-	goto L170;
+        goto L170;
     }
     for (j = 1; j <= n; ++j) {
-/* Computing MAX */
-	d__1 = diag[j], d__2 = wa2[j];
-	diag[j] = max(d__1,d__2);
-/* L160: */
+        /* Computing MAX */
+        d__1 = diag[j], d__2 = wa2[j];
+        diag[j] = max(d__1,d__2);
+        /* L160: */
     }
 L170:
 
-/*        beginning of the inner loop. */
+    /*        beginning of the inner loop. */
 
 L180:
 
-/*           if requested, call fcn to enable printing of iterates. */
+    /*           if requested, call fcn to enable printing of iterates. */
 
     if (nprint <= 0) {
-	goto L190;
+        goto L190;
     }
     iflag = 0;
     if ((iter - 1) % nprint == 0) {
-	iflag = (*fcn)(p, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 0);
+        iflag = (*fcn)(p, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 0);
     }
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
 L190:
 
-/*           determine the direction p. */
+    /*           determine the direction p. */
 
     dogleg(n, &r__[1], lr, &diag[1], &qtf[1], delta, &wa1[1], &wa2[1], &wa3[
-	    1]);
+            1]);
 
-/*           store the direction p and x + p. calculate the norm of p. */
+    /*           store the direction p and x + p. calculate the norm of p. */
 
     for (j = 1; j <= n; ++j) {
-	wa1[j] = -wa1[j];
-	wa2[j] = x[j] + wa1[j];
-	wa3[j] = diag[j] * wa1[j];
-/* L200: */
+        wa1[j] = -wa1[j];
+        wa2[j] = x[j] + wa1[j];
+        wa3[j] = diag[j] * wa1[j];
+        /* L200: */
     }
     pnorm = ei_enorm<Scalar>(n, &wa3[1]);
 
-/*           on the first iteration, adjust the initial step bound. */
+    /*           on the first iteration, adjust the initial step bound. */
 
     if (iter == 1) {
-	delta = min(delta,pnorm);
+        delta = min(delta,pnorm);
     }
 
-/*           evaluate the function at x + p and calculate its norm. */
+    /*           evaluate the function at x + p and calculate its norm. */
 
     iflag = (*fcn)(p, n, &wa2[1], &wa4[1], &fjac[fjac_offset], ldfjac, 1);
     ++(*nfev);
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
     fnorm1 = ei_enorm<Scalar>(n, &wa4[1]);
 
-/*           compute the scaled actual reduction. */
+    /*           compute the scaled actual reduction. */
 
     actred = -1.;
     if (fnorm1 < fnorm) {
-/* Computing 2nd power */
-	d__1 = fnorm1 / fnorm;
-	actred = 1. - d__1 * d__1;
+        /* Computing 2nd power */
+        d__1 = fnorm1 / fnorm;
+        actred = 1. - d__1 * d__1;
     }
 
-/*           compute the scaled predicted reduction. */
+    /*           compute the scaled predicted reduction. */
 
     l = 1;
     for (i__ = 1; i__ <= n; ++i__) {
-	sum = 0.;
-	for (j = i__; j <= n; ++j) {
-	    sum += r__[l] * wa1[j];
-	    ++l;
-/* L210: */
-	}
-	wa3[i__] = qtf[i__] + sum;
-/* L220: */
+        sum = 0.;
+        for (j = i__; j <= n; ++j) {
+            sum += r__[l] * wa1[j];
+            ++l;
+            /* L210: */
+        }
+        wa3[i__] = qtf[i__] + sum;
+        /* L220: */
     }
     temp = ei_enorm<Scalar>(n, &wa3[1]);
     prered = 0.;
     if (temp < fnorm) {
-/* Computing 2nd power */
-	d__1 = temp / fnorm;
-	prered = 1. - d__1 * d__1;
+        /* Computing 2nd power */
+        d__1 = temp / fnorm;
+        prered = 1. - d__1 * d__1;
     }
 
-/*           compute the ratio of the actual to the predicted */
-/*           reduction. */
+    /*           compute the ratio of the actual to the predicted */
+    /*           reduction. */
 
     ratio = 0.;
     if (prered > 0.) {
-	ratio = actred / prered;
+        ratio = actred / prered;
     }
 
-/*           update the step bound. */
+    /*           update the step bound. */
 
     if (ratio >= p1) {
-	goto L230;
+        goto L230;
     }
     ncsuc = 0;
     ++ncfail;
@@ -300,127 +300,127 @@ L230:
     ncfail = 0;
     ++ncsuc;
     if (ratio >= p5 || ncsuc > 1) {
-/* Computing MAX */
-	d__1 = delta, d__2 = pnorm / p5;
-	delta = max(d__1,d__2);
+        /* Computing MAX */
+        d__1 = delta, d__2 = pnorm / p5;
+        delta = max(d__1,d__2);
     }
     if (fabs(ratio - 1.) <= p1) {
-	delta = pnorm / p5;
+        delta = pnorm / p5;
     }
 L240:
 
-/*           test for successful iteration. */
+    /*           test for successful iteration. */
 
     if (ratio < p0001) {
-	goto L260;
+        goto L260;
     }
 
-/*           successful iteration. update x, fvec, and their norms. */
+    /*           successful iteration. update x, fvec, and their norms. */
 
     for (j = 1; j <= n; ++j) {
-	x[j] = wa2[j];
-	wa2[j] = diag[j] * x[j];
-	fvec[j] = wa4[j];
-/* L250: */
+        x[j] = wa2[j];
+        wa2[j] = diag[j] * x[j];
+        fvec[j] = wa4[j];
+        /* L250: */
     }
     xnorm = ei_enorm<Scalar>(n, &wa2[1]);
     fnorm = fnorm1;
     ++iter;
 L260:
 
-/*           determine the progress of the iteration. */
+    /*           determine the progress of the iteration. */
 
     ++nslow1;
     if (actred >= p001) {
-	nslow1 = 0;
+        nslow1 = 0;
     }
     if (jeval) {
-	++nslow2;
+        ++nslow2;
     }
     if (actred >= p1) {
-	nslow2 = 0;
+        nslow2 = 0;
     }
 
-/*           test for convergence. */
+    /*           test for convergence. */
 
     if (delta <= xtol * xnorm || fnorm == 0.) {
-	info = 1;
+        info = 1;
     }
     if (info != 0) {
-	goto L300;
+        goto L300;
     }
 
-/*           tests for termination and stringent tolerances. */
+    /*           tests for termination and stringent tolerances. */
 
     if (*nfev >= maxfev) {
-	info = 2;
+        info = 2;
     }
-/* Computing MAX */
+    /* Computing MAX */
     d__1 = p1 * delta;
     if (p1 * max(d__1,pnorm) <= epsilon<Scalar>() * xnorm) {
-	info = 3;
+        info = 3;
     }
     if (nslow2 == 5) {
-	info = 4;
+        info = 4;
     }
     if (nslow1 == 10) {
-	info = 5;
+        info = 5;
     }
     if (info != 0) {
-	goto L300;
+        goto L300;
     }
 
-/*           criterion for recalculating jacobian. */
+    /*           criterion for recalculating jacobian. */
 
     if (ncfail == 2) {
-	goto L290;
+        goto L290;
     }
 
-/*           calculate the rank one modification to the jacobian */
-/*           and update qtf if necessary. */
+    /*           calculate the rank one modification to the jacobian */
+    /*           and update qtf if necessary. */
 
     for (j = 1; j <= n; ++j) {
-	sum = 0.;
-	for (i__ = 1; i__ <= n; ++i__) {
-	    sum += fjac[i__ + j * fjac_dim1] * wa4[i__];
-/* L270: */
-	}
-	wa2[j] = (sum - wa3[j]) / pnorm;
-	wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
-	if (ratio >= p0001) {
-	    qtf[j] = sum;
-	}
-/* L280: */
+        sum = 0.;
+        for (i__ = 1; i__ <= n; ++i__) {
+            sum += fjac[i__ + j * fjac_dim1] * wa4[i__];
+            /* L270: */
+        }
+        wa2[j] = (sum - wa3[j]) / pnorm;
+        wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
+        if (ratio >= p0001) {
+            qtf[j] = sum;
+        }
+        /* L280: */
     }
 
-/*           compute the qr factorization of the updated jacobian. */
+    /*           compute the qr factorization of the updated jacobian. */
 
     r1updt(n, n, &r__[1], lr, &wa1[1], &wa2[1], &wa3[1], &sing);
     r1mpyq(n, n, &fjac[fjac_offset], ldfjac, &wa2[1], &wa3[1]);
     r1mpyq(1, n, &qtf[1], 1, &wa2[1], &wa3[1]);
 
-/*           end of the inner loop. */
+    /*           end of the inner loop. */
 
     jeval = FALSE_;
     goto L180;
 L290:
 
-/*        end of the outer loop. */
+    /*        end of the outer loop. */
 
     goto L30;
 L300:
 
-/*     termination, either normal or user imposed. */
+    /*     termination, either normal or user imposed. */
 
     if (iflag < 0) {
-	info = iflag;
+        info = iflag;
     }
     if (nprint > 0) {
-	iflag = (*fcn)(p, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 0);
+        iflag = (*fcn)(p, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 0);
     }
     return info;
 
-/*     last card of subroutine hybrj. */
+    /*     last card of subroutine hybrj. */
 
 } /* hybrj_ */
 
diff --git a/unsupported/Eigen/src/NonLinear/lmder.h b/unsupported/Eigen/src/NonLinear/lmder.h
index e66b758d1..1e8289557 100644
--- a/unsupported/Eigen/src/NonLinear/lmder.h
+++ b/unsupported/Eigen/src/NonLinear/lmder.h
@@ -1,11 +1,11 @@
 
 template<typename Scalar>
 int lmder_template(minpack_funcder_mn fcn, void *p, int m, int n, Scalar *x, 
-	Scalar *fvec, Scalar *fjac, int ldfjac, Scalar ftol,
-	Scalar xtol, Scalar gtol, int maxfev, Scalar *
-	diag, int mode, Scalar factor, int nprint,
-	int *nfev, int *njev, int *ipvt, Scalar *qtf, 
-	Scalar *wa1, Scalar *wa2, Scalar *wa3, Scalar *wa4)
+        Scalar *fvec, Scalar *fjac, int ldfjac, Scalar ftol,
+        Scalar xtol, Scalar gtol, int maxfev, Scalar *
+        diag, int mode, Scalar factor, int nprint,
+        int *nfev, int *njev, int *ipvt, Scalar *qtf, 
+        Scalar *wa1, Scalar *wa2, Scalar *wa3, Scalar *wa4)
 {
     /* Initialized data */
 
@@ -45,359 +45,359 @@ int lmder_template(minpack_funcder_mn fcn, void *p, int m, int n, Scalar *x,
     *nfev = 0;
     *njev = 0;
 
-/*     check the input parameters for errors. */
+    /*     check the input parameters for errors. */
 
     if (n <= 0 || m < n || ldfjac < m || ftol < 0. || xtol < 0. || 
-	    gtol < 0. || maxfev <= 0 || factor <= 0.) {
-	goto L300;
+            gtol < 0. || maxfev <= 0 || factor <= 0.) {
+        goto L300;
     }
     if (mode != 2) {
-	goto L20;
+        goto L20;
     }
     for (j = 1; j <= n; ++j) {
-	if (diag[j] <= 0.) {
-	    goto L300;
-	}
-/* L10: */
+        if (diag[j] <= 0.) {
+            goto L300;
+        }
+        /* L10: */
     }
 L20:
 
-/*     evaluate the function at the starting point */
-/*     and calculate its norm. */
+    /*     evaluate the function at the starting point */
+    /*     and calculate its norm. */
 
     iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 1);
     *nfev = 1;
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
     fnorm = ei_enorm<Scalar>(m, &fvec[1]);
 
-/*     initialize levenberg-marquardt parameter and iteration counter. */
+    /*     initialize levenberg-marquardt parameter and iteration counter. */
 
     par = 0.;
     iter = 1;
 
-/*     beginning of the outer loop. */
+    /*     beginning of the outer loop. */
 
 L30:
 
-/*        calculate the jacobian matrix. */
+    /*        calculate the jacobian matrix. */
 
     iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 2);
     ++(*njev);
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
 
-/*        if requested, call fcn to enable printing of iterates. */
+    /*        if requested, call fcn to enable printing of iterates. */
 
     if (nprint <= 0) {
-	goto L40;
+        goto L40;
     }
     iflag = 0;
     if ((iter - 1) % nprint == 0) {
-	iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 0);
+        iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 0);
     }
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
 L40:
 
-/*        compute the qr factorization of the jacobian. */
+    /*        compute the qr factorization of the jacobian. */
 
     qrfac(m, n, &fjac[fjac_offset], ldfjac, TRUE_, &ipvt[1], n, &wa1[1], &
-	    wa2[1], &wa3[1]);
+            wa2[1], &wa3[1]);
 
-/*        on the first iteration and if mode is 1, scale according */
-/*        to the norms of the columns of the initial jacobian. */
+    /*        on the first iteration and if mode is 1, scale according */
+    /*        to the norms of the columns of the initial jacobian. */
 
     if (iter != 1) {
-	goto L80;
+        goto L80;
     }
     if (mode == 2) {
-	goto L60;
+        goto L60;
     }
     for (j = 1; j <= n; ++j) {
-	diag[j] = wa2[j];
-	if (wa2[j] == 0.) {
-	    diag[j] = 1.;
-	}
-/* L50: */
+        diag[j] = wa2[j];
+        if (wa2[j] == 0.) {
+            diag[j] = 1.;
+        }
+        /* L50: */
     }
 L60:
 
-/*        on the first iteration, calculate the norm of the scaled x */
-/*        and initialize the step bound delta. */
+    /*        on the first iteration, calculate the norm of the scaled x */
+    /*        and initialize the step bound delta. */
 
     for (j = 1; j <= n; ++j) {
-	wa3[j] = diag[j] * x[j];
-/* L70: */
+        wa3[j] = diag[j] * x[j];
+        /* L70: */
     }
     xnorm = ei_enorm<Scalar>(n, &wa3[1]);
     delta = factor * xnorm;
     if (delta == 0.) {
-	delta = factor;
+        delta = factor;
     }
 L80:
 
-/*        form (q transpose)*fvec and store the first n components in */
-/*        qtf. */
+    /*        form (q transpose)*fvec and store the first n components in */
+    /*        qtf. */
 
     for (i__ = 1; i__ <= m; ++i__) {
-	wa4[i__] = fvec[i__];
-/* L90: */
+        wa4[i__] = fvec[i__];
+        /* L90: */
     }
     for (j = 1; j <= n; ++j) {
-	if (fjac[j + j * fjac_dim1] == 0.) {
-	    goto L120;
-	}
-	sum = 0.;
-	for (i__ = j; i__ <= m; ++i__) {
-	    sum += fjac[i__ + j * fjac_dim1] * wa4[i__];
-/* L100: */
-	}
-	temp = -sum / fjac[j + j * fjac_dim1];
-	for (i__ = j; i__ <= m; ++i__) {
-	    wa4[i__] += fjac[i__ + j * fjac_dim1] * temp;
-/* L110: */
-	}
+        if (fjac[j + j * fjac_dim1] == 0.) {
+            goto L120;
+        }
+        sum = 0.;
+        for (i__ = j; i__ <= m; ++i__) {
+            sum += fjac[i__ + j * fjac_dim1] * wa4[i__];
+            /* L100: */
+        }
+        temp = -sum / fjac[j + j * fjac_dim1];
+        for (i__ = j; i__ <= m; ++i__) {
+            wa4[i__] += fjac[i__ + j * fjac_dim1] * temp;
+            /* L110: */
+        }
 L120:
-	fjac[j + j * fjac_dim1] = wa1[j];
-	qtf[j] = wa4[j];
-/* L130: */
+        fjac[j + j * fjac_dim1] = wa1[j];
+        qtf[j] = wa4[j];
+        /* L130: */
     }
 
-/*        compute the norm of the scaled gradient. */
+    /*        compute the norm of the scaled gradient. */
 
     gnorm = 0.;
     if (fnorm == 0.) {
-	goto L170;
+        goto L170;
     }
     for (j = 1; j <= n; ++j) {
-	l = ipvt[j];
-	if (wa2[l] == 0.) {
-	    goto L150;
-	}
-	sum = 0.;
-	for (i__ = 1; i__ <= j; ++i__) {
-	    sum += fjac[i__ + j * fjac_dim1] * (qtf[i__] / fnorm);
-/* L140: */
-	}
-/* Computing MAX */
-	d__2 = gnorm, d__3 = fabs(sum / wa2[l]);
-	gnorm = max(d__2,d__3);
+        l = ipvt[j];
+        if (wa2[l] == 0.) {
+            goto L150;
+        }
+        sum = 0.;
+        for (i__ = 1; i__ <= j; ++i__) {
+            sum += fjac[i__ + j * fjac_dim1] * (qtf[i__] / fnorm);
+            /* L140: */
+        }
+        /* Computing MAX */
+        d__2 = gnorm, d__3 = fabs(sum / wa2[l]);
+        gnorm = max(d__2,d__3);
 L150:
-/* L160: */
-	;
+        /* L160: */
+        ;
     }
 L170:
 
-/*        test for convergence of the gradient norm. */
+    /*        test for convergence of the gradient norm. */
 
     if (gnorm <= gtol) {
-	info = 4;
+        info = 4;
     }
     if (info != 0) {
-	goto L300;
+        goto L300;
     }
 
-/*        rescale if necessary. */
+    /*        rescale if necessary. */
 
     if (mode == 2) {
-	goto L190;
+        goto L190;
     }
     for (j = 1; j <= n; ++j) {
-/* Computing MAX */
-	d__1 = diag[j], d__2 = wa2[j];
-	diag[j] = max(d__1,d__2);
-/* L180: */
+        /* Computing MAX */
+        d__1 = diag[j], d__2 = wa2[j];
+        diag[j] = max(d__1,d__2);
+        /* L180: */
     }
 L190:
 
-/*        beginning of the inner loop. */
+    /*        beginning of the inner loop. */
 
 L200:
 
-/*           determine the levenberg-marquardt parameter. */
+    /*           determine the levenberg-marquardt parameter. */
 
     lmpar(n, &fjac[fjac_offset], ldfjac, &ipvt[1], &diag[1], &qtf[1], delta,
-	     &par, &wa1[1], &wa2[1], &wa3[1], &wa4[1]);
+            &par, &wa1[1], &wa2[1], &wa3[1], &wa4[1]);
 
-/*           store the direction p and x + p. calculate the norm of p. */
+    /*           store the direction p and x + p. calculate the norm of p. */
 
     for (j = 1; j <= n; ++j) {
-	wa1[j] = -wa1[j];
-	wa2[j] = x[j] + wa1[j];
-	wa3[j] = diag[j] * wa1[j];
-/* L210: */
+        wa1[j] = -wa1[j];
+        wa2[j] = x[j] + wa1[j];
+        wa3[j] = diag[j] * wa1[j];
+        /* L210: */
     }
     pnorm = ei_enorm<Scalar>(n, &wa3[1]);
 
-/*           on the first iteration, adjust the initial step bound. */
+    /*           on the first iteration, adjust the initial step bound. */
 
     if (iter == 1) {
-	delta = min(delta,pnorm);
+        delta = min(delta,pnorm);
     }
 
-/*           evaluate the function at x + p and calculate its norm. */
+    /*           evaluate the function at x + p and calculate its norm. */
 
     iflag = (*fcn)(p, m, n, &wa2[1], &wa4[1], &fjac[fjac_offset], ldfjac, 1);
     ++(*nfev);
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
     fnorm1 = ei_enorm<Scalar>(m, &wa4[1]);
 
-/*           compute the scaled actual reduction. */
+    /*           compute the scaled actual reduction. */
 
     actred = -1.;
     if (p1 * fnorm1 < fnorm) {
-/* Computing 2nd power */
-	d__1 = fnorm1 / fnorm;
-	actred = 1. - d__1 * d__1;
+        /* Computing 2nd power */
+        d__1 = fnorm1 / fnorm;
+        actred = 1. - d__1 * d__1;
     }
 
-/*           compute the scaled predicted reduction and */
-/*           the scaled directional derivative. */
+    /*           compute the scaled predicted reduction and */
+    /*           the scaled directional derivative. */
 
     for (j = 1; j <= n; ++j) {
-	wa3[j] = 0.;
-	l = ipvt[j];
-	temp = wa1[l];
-	for (i__ = 1; i__ <= j; ++i__) {
-	    wa3[i__] += fjac[i__ + j * fjac_dim1] * temp;
-/* L220: */
-	}
-/* L230: */
+        wa3[j] = 0.;
+        l = ipvt[j];
+        temp = wa1[l];
+        for (i__ = 1; i__ <= j; ++i__) {
+            wa3[i__] += fjac[i__ + j * fjac_dim1] * temp;
+            /* L220: */
+        }
+        /* L230: */
     }
     temp1 = ei_enorm<Scalar>(n, &wa3[1]) / fnorm;
     temp2 = sqrt(par) * pnorm / fnorm;
-/* Computing 2nd power */
+    /* Computing 2nd power */
     d__1 = temp1;
-/* Computing 2nd power */
+    /* Computing 2nd power */
     d__2 = temp2;
     prered = d__1 * d__1 + d__2 * d__2 / p5;
-/* Computing 2nd power */
+    /* Computing 2nd power */
     d__1 = temp1;
-/* Computing 2nd power */
+    /* Computing 2nd power */
     d__2 = temp2;
     dirder = -(d__1 * d__1 + d__2 * d__2);
 
-/*           compute the ratio of the actual to the predicted */
-/*           reduction. */
+    /*           compute the ratio of the actual to the predicted */
+    /*           reduction. */
 
     ratio = 0.;
     if (prered != 0.) {
-	ratio = actred / prered;
+        ratio = actred / prered;
     }
 
-/*           update the step bound. */
+    /*           update the step bound. */
 
     if (ratio > p25) {
-	goto L240;
+        goto L240;
     }
     if (actred >= 0.) {
-	temp = p5;
+        temp = p5;
     }
     if (actred < 0.) {
-	temp = p5 * dirder / (dirder + p5 * actred);
+        temp = p5 * dirder / (dirder + p5 * actred);
     }
     if (p1 * fnorm1 >= fnorm || temp < p1) {
-	temp = p1;
+        temp = p1;
     }
-/* Computing MIN */
+    /* Computing MIN */
     d__1 = delta, d__2 = pnorm / p1;
     delta = temp * min(d__1,d__2);
     par /= temp;
     goto L260;
 L240:
     if (par != 0. && ratio < p75) {
-	goto L250;
+        goto L250;
     }
     delta = pnorm / p5;
     par = p5 * par;
 L250:
 L260:
 
-/*           test for successful iteration. */
+    /*           test for successful iteration. */
 
     if (ratio < p0001) {
-	goto L290;
+        goto L290;
     }
 
-/*           successful iteration. update x, fvec, and their norms. */
+    /*           successful iteration. update x, fvec, and their norms. */
 
     for (j = 1; j <= n; ++j) {
-	x[j] = wa2[j];
-	wa2[j] = diag[j] * x[j];
-/* L270: */
+        x[j] = wa2[j];
+        wa2[j] = diag[j] * x[j];
+        /* L270: */
     }
     for (i__ = 1; i__ <= m; ++i__) {
-	fvec[i__] = wa4[i__];
-/* L280: */
+        fvec[i__] = wa4[i__];
+        /* L280: */
     }
     xnorm = ei_enorm<Scalar>(n, &wa2[1]);
     fnorm = fnorm1;
     ++iter;
 L290:
 
-/*           tests for convergence. */
+    /*           tests for convergence. */
 
     if (fabs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1.) {
-	info = 1;
+        info = 1;
     }
     if (delta <= xtol * xnorm) {
-	info = 2;
+        info = 2;
     }
     if (fabs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1. && info 
-	    == 2) {
-	info = 3;
+            == 2) {
+        info = 3;
     }
     if (info != 0) {
-	goto L300;
+        goto L300;
     }
 
-/*           tests for termination and stringent tolerances. */
+    /*           tests for termination and stringent tolerances. */
 
     if (*nfev >= maxfev) {
-	info = 5;
+        info = 5;
     }
     if (fabs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && p5 * ratio <= 1.) {
-	info = 6;
+        info = 6;
     }
     if (delta <= epsilon<Scalar>() * xnorm) {
-	info = 7;
+        info = 7;
     }
     if (gnorm <= epsilon<Scalar>()) {
-	info = 8;
+        info = 8;
     }
     if (info != 0) {
-	goto L300;
+        goto L300;
     }
 
-/*           end of the inner loop. repeat if iteration unsuccessful. */
+    /*           end of the inner loop. repeat if iteration unsuccessful. */
 
     if (ratio < p0001) {
-	goto L200;
+        goto L200;
     }
 
-/*        end of the outer loop. */
+    /*        end of the outer loop. */
 
     goto L30;
 L300:
 
-/*     termination, either normal or user imposed. */
+    /*     termination, either normal or user imposed. */
 
     if (iflag < 0) {
-	info = iflag;
+        info = iflag;
     }
     iflag = 0;
     if (nprint > 0) {
-	iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 0);
+        iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 0);
     }
     return info;
 
-/*     last card of subroutine lmder. */
+    /*     last card of subroutine lmder. */
 
 } /* lmder_ */
 
diff --git a/unsupported/Eigen/src/NonLinear/lmdif.h b/unsupported/Eigen/src/NonLinear/lmdif.h
index 0ed91f2d1..26296b631 100644
--- a/unsupported/Eigen/src/NonLinear/lmdif.h
+++ b/unsupported/Eigen/src/NonLinear/lmdif.h
@@ -1,12 +1,12 @@
 
 template<typename Scalar>
 int lmdif_template(minpack_func_mn fcn, void *p, int m, int n, Scalar *x, 
-	Scalar *fvec, Scalar ftol, Scalar xtol, Scalar gtol,
-    int maxfev, Scalar epsfcn, Scalar *diag, int
-	mode, Scalar factor, int nprint, int *
-	nfev, Scalar *fjac, int ldfjac, int *ipvt, Scalar *
-	qtf, Scalar *wa1, Scalar *wa2, Scalar *wa3, Scalar *
-	wa4)
+        Scalar *fvec, Scalar ftol, Scalar xtol, Scalar gtol,
+        int maxfev, Scalar epsfcn, Scalar *diag, int
+        mode, Scalar factor, int nprint, int *
+        nfev, Scalar *fjac, int ldfjac, int *ipvt, Scalar *
+        qtf, Scalar *wa1, Scalar *wa2, Scalar *wa3, Scalar *
+        wa4)
 {
     /* Initialized data */
 
@@ -46,360 +46,360 @@ int lmdif_template(minpack_func_mn fcn, void *p, int m, int n, Scalar *x,
     iflag = 0;
     *nfev = 0;
 
-/*     check the input parameters for errors. */
+    /*     check the input parameters for errors. */
 
     if (n <= 0 || m < n || ldfjac < m || ftol < 0. || xtol < 0. || 
-	    gtol < 0. || maxfev <= 0 || factor <= 0.) {
-	goto L300;
+            gtol < 0. || maxfev <= 0 || factor <= 0.) {
+        goto L300;
     }
     if (mode != 2) {
-	goto L20;
+        goto L20;
     }
     for (j = 1; j <= n; ++j) {
-	if (diag[j] <= 0.) {
-	    goto L300;
-	}
-/* L10: */
+        if (diag[j] <= 0.) {
+            goto L300;
+        }
+        /* L10: */
     }
 L20:
 
-/*     evaluate the function at the starting point */
-/*     and calculate its norm. */
+    /*     evaluate the function at the starting point */
+    /*     and calculate its norm. */
 
     iflag = (*fcn)(p, m, n, &x[1], &fvec[1], 1);
     *nfev = 1;
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
     fnorm = ei_enorm<Scalar>(m, &fvec[1]);
 
-/*     initialize levenberg-marquardt parameter and iteration counter. */
+    /*     initialize levenberg-marquardt parameter and iteration counter. */
 
     par = 0.;
     iter = 1;
 
-/*     beginning of the outer loop. */
+    /*     beginning of the outer loop. */
 
 L30:
 
-/*        calculate the jacobian matrix. */
+    /*        calculate the jacobian matrix. */
 
     iflag = fdjac2(fcn, p, m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac,
-	    epsfcn, &wa4[1]);
+            epsfcn, &wa4[1]);
     *nfev += n;
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
 
-/*        if requested, call fcn to enable printing of iterates. */
+    /*        if requested, call fcn to enable printing of iterates. */
 
     if (nprint <= 0) {
-	goto L40;
+        goto L40;
     }
     iflag = 0;
     if ((iter - 1) % nprint == 0) {
-	iflag = (*fcn)(p, m, n, &x[1], &fvec[1], 0);
+        iflag = (*fcn)(p, m, n, &x[1], &fvec[1], 0);
     }
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
 L40:
 
-/*        compute the qr factorization of the jacobian. */
+    /*        compute the qr factorization of the jacobian. */
 
     qrfac(m, n, &fjac[fjac_offset], ldfjac, TRUE_, &ipvt[1], n, &wa1[1], &
-	    wa2[1], &wa3[1]);
+            wa2[1], &wa3[1]);
 
-/*        on the first iteration and if mode is 1, scale according */
-/*        to the norms of the columns of the initial jacobian. */
+    /*        on the first iteration and if mode is 1, scale according */
+    /*        to the norms of the columns of the initial jacobian. */
 
     if (iter != 1) {
-	goto L80;
+        goto L80;
     }
     if (mode == 2) {
-	goto L60;
+        goto L60;
     }
     for (j = 1; j <= n; ++j) {
-	diag[j] = wa2[j];
-	if (wa2[j] == 0.) {
-	    diag[j] = 1.;
-	}
-/* L50: */
+        diag[j] = wa2[j];
+        if (wa2[j] == 0.) {
+            diag[j] = 1.;
+        }
+        /* L50: */
     }
 L60:
 
-/*        on the first iteration, calculate the norm of the scaled x */
-/*        and initialize the step bound delta. */
+    /*        on the first iteration, calculate the norm of the scaled x */
+    /*        and initialize the step bound delta. */
 
     for (j = 1; j <= n; ++j) {
-	wa3[j] = diag[j] * x[j];
-/* L70: */
+        wa3[j] = diag[j] * x[j];
+        /* L70: */
     }
     xnorm = ei_enorm<Scalar>(n, &wa3[1]);
     delta = factor * xnorm;
     if (delta == 0.) {
-	delta = factor;
+        delta = factor;
     }
 L80:
 
-/*        form (q transpose)*fvec and store the first n components in */
-/*        qtf. */
+    /*        form (q transpose)*fvec and store the first n components in */
+    /*        qtf. */
 
     for (i__ = 1; i__ <= m; ++i__) {
-	wa4[i__] = fvec[i__];
-/* L90: */
+        wa4[i__] = fvec[i__];
+        /* L90: */
     }
     for (j = 1; j <= n; ++j) {
-	if (fjac[j + j * fjac_dim1] == 0.) {
-	    goto L120;
-	}
-	sum = 0.;
-	for (i__ = j; i__ <= m; ++i__) {
-	    sum += fjac[i__ + j * fjac_dim1] * wa4[i__];
-/* L100: */
-	}
-	temp = -sum / fjac[j + j * fjac_dim1];
-	for (i__ = j; i__ <= m; ++i__) {
-	    wa4[i__] += fjac[i__ + j * fjac_dim1] * temp;
-/* L110: */
-	}
+        if (fjac[j + j * fjac_dim1] == 0.) {
+            goto L120;
+        }
+        sum = 0.;
+        for (i__ = j; i__ <= m; ++i__) {
+            sum += fjac[i__ + j * fjac_dim1] * wa4[i__];
+            /* L100: */
+        }
+        temp = -sum / fjac[j + j * fjac_dim1];
+        for (i__ = j; i__ <= m; ++i__) {
+            wa4[i__] += fjac[i__ + j * fjac_dim1] * temp;
+            /* L110: */
+        }
 L120:
-	fjac[j + j * fjac_dim1] = wa1[j];
-	qtf[j] = wa4[j];
-/* L130: */
+        fjac[j + j * fjac_dim1] = wa1[j];
+        qtf[j] = wa4[j];
+        /* L130: */
     }
 
-/*        compute the norm of the scaled gradient. */
+    /*        compute the norm of the scaled gradient. */
 
     gnorm = 0.;
     if (fnorm == 0.) {
-	goto L170;
+        goto L170;
     }
     for (j = 1; j <= n; ++j) {
-	l = ipvt[j];
-	if (wa2[l] == 0.) {
-	    goto L150;
-	}
-	sum = 0.;
-	for (i__ = 1; i__ <= j; ++i__) {
-	    sum += fjac[i__ + j * fjac_dim1] * (qtf[i__] / fnorm);
-/* L140: */
-	}
-/* Computing MAX */
-	d__2 = gnorm, d__3 = fabs(sum / wa2[l]);
-	gnorm = max(d__2,d__3);
+        l = ipvt[j];
+        if (wa2[l] == 0.) {
+            goto L150;
+        }
+        sum = 0.;
+        for (i__ = 1; i__ <= j; ++i__) {
+            sum += fjac[i__ + j * fjac_dim1] * (qtf[i__] / fnorm);
+            /* L140: */
+        }
+        /* Computing MAX */
+        d__2 = gnorm, d__3 = fabs(sum / wa2[l]);
+        gnorm = max(d__2,d__3);
 L150:
-/* L160: */
-	;
+        /* L160: */
+        ;
     }
 L170:
 
-/*        test for convergence of the gradient norm. */
+    /*        test for convergence of the gradient norm. */
 
     if (gnorm <= gtol) {
-	info = 4;
+        info = 4;
     }
     if (info != 0) {
-	goto L300;
+        goto L300;
     }
 
-/*        rescale if necessary. */
+    /*        rescale if necessary. */
 
     if (mode == 2) {
-	goto L190;
+        goto L190;
     }
     for (j = 1; j <= n; ++j) {
-/* Computing MAX */
-	d__1 = diag[j], d__2 = wa2[j];
-	diag[j] = max(d__1,d__2);
-/* L180: */
+        /* Computing MAX */
+        d__1 = diag[j], d__2 = wa2[j];
+        diag[j] = max(d__1,d__2);
+        /* L180: */
     }
 L190:
 
-/*        beginning of the inner loop. */
+    /*        beginning of the inner loop. */
 
 L200:
 
-/*           determine the levenberg-marquardt parameter. */
+    /*           determine the levenberg-marquardt parameter. */
 
     lmpar(n, &fjac[fjac_offset], ldfjac, &ipvt[1], &diag[1], &qtf[1], delta,
-	     &par, &wa1[1], &wa2[1], &wa3[1], &wa4[1]);
+            &par, &wa1[1], &wa2[1], &wa3[1], &wa4[1]);
 
-/*           store the direction p and x + p. calculate the norm of p. */
+    /*           store the direction p and x + p. calculate the norm of p. */
 
     for (j = 1; j <= n; ++j) {
-	wa1[j] = -wa1[j];
-	wa2[j] = x[j] + wa1[j];
-	wa3[j] = diag[j] * wa1[j];
-/* L210: */
+        wa1[j] = -wa1[j];
+        wa2[j] = x[j] + wa1[j];
+        wa3[j] = diag[j] * wa1[j];
+        /* L210: */
     }
     pnorm = ei_enorm<Scalar>(n, &wa3[1]);
 
-/*           on the first iteration, adjust the initial step bound. */
+    /*           on the first iteration, adjust the initial step bound. */
 
     if (iter == 1) {
-	delta = min(delta,pnorm);
+        delta = min(delta,pnorm);
     }
 
-/*           evaluate the function at x + p and calculate its norm. */
+    /*           evaluate the function at x + p and calculate its norm. */
 
     iflag = (*fcn)(p, m, n, &wa2[1], &wa4[1], 1);
     ++(*nfev);
     if (iflag < 0) {
-	goto L300;
+        goto L300;
     }
     fnorm1 = ei_enorm<Scalar>(m, &wa4[1]);
 
-/*           compute the scaled actual reduction. */
+    /*           compute the scaled actual reduction. */
 
     actred = -1.;
     if (p1 * fnorm1 < fnorm) {
-/* Computing 2nd power */
-	d__1 = fnorm1 / fnorm;
-	actred = 1. - d__1 * d__1;
+        /* Computing 2nd power */
+        d__1 = fnorm1 / fnorm;
+        actred = 1. - d__1 * d__1;
     }
 
-/*           compute the scaled predicted reduction and */
-/*           the scaled directional derivative. */
+    /*           compute the scaled predicted reduction and */
+    /*           the scaled directional derivative. */
 
     for (j = 1; j <= n; ++j) {
-	wa3[j] = 0.;
-	l = ipvt[j];
-	temp = wa1[l];
-	for (i__ = 1; i__ <= j; ++i__) {
-	    wa3[i__] += fjac[i__ + j * fjac_dim1] * temp;
-/* L220: */
-	}
-/* L230: */
+        wa3[j] = 0.;
+        l = ipvt[j];
+        temp = wa1[l];
+        for (i__ = 1; i__ <= j; ++i__) {
+            wa3[i__] += fjac[i__ + j * fjac_dim1] * temp;
+            /* L220: */
+        }
+        /* L230: */
     }
     temp1 = ei_enorm<Scalar>(n, &wa3[1]) / fnorm;
     temp2 = sqrt(par) * pnorm / fnorm;
-/* Computing 2nd power */
+    /* Computing 2nd power */
     d__1 = temp1;
-/* Computing 2nd power */
+    /* Computing 2nd power */
     d__2 = temp2;
     prered = d__1 * d__1 + d__2 * d__2 / p5;
-/* Computing 2nd power */
+    /* Computing 2nd power */
     d__1 = temp1;
-/* Computing 2nd power */
+    /* Computing 2nd power */
     d__2 = temp2;
     dirder = -(d__1 * d__1 + d__2 * d__2);
 
-/*           compute the ratio of the actual to the predicted */
-/*           reduction. */
+    /*           compute the ratio of the actual to the predicted */
+    /*           reduction. */
 
     ratio = 0.;
     if (prered != 0.) {
-	ratio = actred / prered;
+        ratio = actred / prered;
     }
 
-/*           update the step bound. */
+    /*           update the step bound. */
 
     if (ratio > p25) {
-	goto L240;
+        goto L240;
     }
     if (actred >= 0.) {
-	temp = p5;
+        temp = p5;
     }
     if (actred < 0.) {
-	temp = p5 * dirder / (dirder + p5 * actred);
+        temp = p5 * dirder / (dirder + p5 * actred);
     }
     if (p1 * fnorm1 >= fnorm || temp < p1) {
-	temp = p1;
+        temp = p1;
     }
-/* Computing MIN */
+    /* Computing MIN */
     d__1 = delta, d__2 = pnorm / p1;
     delta = temp * min(d__1,d__2);
     par /= temp;
     goto L260;
 L240:
     if (par != 0. && ratio < p75) {
-	goto L250;
+        goto L250;
     }
     delta = pnorm / p5;
     par = p5 * par;
 L250:
 L260:
 
-/*           test for successful iteration. */
+    /*           test for successful iteration. */
 
     if (ratio < p0001) {
-	goto L290;
+        goto L290;
     }
 
-/*           successful iteration. update x, fvec, and their norms. */
+    /*           successful iteration. update x, fvec, and their norms. */
 
     for (j = 1; j <= n; ++j) {
-	x[j] = wa2[j];
-	wa2[j] = diag[j] * x[j];
-/* L270: */
+        x[j] = wa2[j];
+        wa2[j] = diag[j] * x[j];
+        /* L270: */
     }
     for (i__ = 1; i__ <= m; ++i__) {
-	fvec[i__] = wa4[i__];
-/* L280: */
+        fvec[i__] = wa4[i__];
+        /* L280: */
     }
     xnorm = ei_enorm<Scalar>(n, &wa2[1]);
     fnorm = fnorm1;
     ++iter;
 L290:
 
-/*           tests for convergence. */
+    /*           tests for convergence. */
 
     if (fabs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1.) {
-	info = 1;
+        info = 1;
     }
     if (delta <= xtol * xnorm) {
-	info = 2;
+        info = 2;
     }
     if (fabs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1. && info 
-	    == 2) {
-	info = 3;
+            == 2) {
+        info = 3;
     }
     if (info != 0) {
-	goto L300;
+        goto L300;
     }
 
-/*           tests for termination and stringent tolerances. */
+    /*           tests for termination and stringent tolerances. */
 
     if (*nfev >= maxfev) {
-	info = 5;
+        info = 5;
     }
     if (fabs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && p5 * ratio <= 1.) {
-	info = 6;
+        info = 6;
     }
     if (delta <= epsilon<Scalar>() * xnorm) {
-	info = 7;
+        info = 7;
     }
     if (gnorm <= epsilon<Scalar>()) {
-	info = 8;
+        info = 8;
     }
     if (info != 0) {
-	goto L300;
+        goto L300;
     }
 
-/*           end of the inner loop. repeat if iteration unsuccessful. */
+    /*           end of the inner loop. repeat if iteration unsuccessful. */
 
     if (ratio < p0001) {
-	goto L200;
+        goto L200;
     }
 
-/*        end of the outer loop. */
+    /*        end of the outer loop. */
 
     goto L30;
 L300:
 
-/*     termination, either normal or user imposed. */
+    /*     termination, either normal or user imposed. */
 
     if (iflag < 0) {
-	info = iflag;
+        info = iflag;
     }
     iflag = 0;
     if (nprint > 0) {
-	iflag = (*fcn)(p, m, n, &x[1], &fvec[1], 0);
+        iflag = (*fcn)(p, m, n, &x[1], &fvec[1], 0);
     }
     return info;
 
-/*     last card of subroutine lmdif. */
+    /*     last card of subroutine lmdif. */
 
 } /* lmdif_ */
 
diff --git a/unsupported/Eigen/src/NonLinear/lmstr.h b/unsupported/Eigen/src/NonLinear/lmstr.h
index 3b939e197..6be32e094 100644
--- a/unsupported/Eigen/src/NonLinear/lmstr.h
+++ b/unsupported/Eigen/src/NonLinear/lmstr.h
@@ -1,11 +1,11 @@
 
 template<typename Scalar>
 int lmstr_template(minpack_funcderstr_mn fcn, void *p, int m, int n, Scalar *x, 
-	Scalar *fvec, Scalar *fjac, int ldfjac, Scalar ftol,
-	Scalar xtol, Scalar gtol, int maxfev, Scalar *
-	diag, int mode, Scalar factor, int nprint,
-	int *nfev, int *njev, int *ipvt, Scalar *qtf, 
-	Scalar *wa1, Scalar *wa2, Scalar *wa3, Scalar *wa4)
+        Scalar *fvec, Scalar *fjac, int ldfjac, Scalar ftol,
+        Scalar xtol, Scalar gtol, int maxfev, Scalar *
+        diag, int mode, Scalar factor, int nprint,
+        int *nfev, int *njev, int *ipvt, Scalar *qtf, 
+        Scalar *wa1, Scalar *wa2, Scalar *wa3, Scalar *wa4)
 {
     /* Initialized data */
 
@@ -46,382 +46,382 @@ int lmstr_template(minpack_funcderstr_mn fcn, void *p, int m, int n, Scalar *x,
     *nfev = 0;
     *njev = 0;
 
-/*     check the input parameters for errors. */
+    /*     check the input parameters for errors. */
 
     if (n <= 0 || m < n || ldfjac < n || ftol < 0. || xtol < 0. || 
-	    gtol < 0. || maxfev <= 0 || factor <= 0.) {
-	goto L340;
+            gtol < 0. || maxfev <= 0 || factor <= 0.) {
+        goto L340;
     }
     if (mode != 2) {
-	goto L20;
+        goto L20;
     }
     for (j = 1; j <= n; ++j) {
-	if (diag[j] <= 0.) {
-	    goto L340;
-	}
-/* L10: */
+        if (diag[j] <= 0.) {
+            goto L340;
+        }
+        /* L10: */
     }
 L20:
 
-/*     evaluate the function at the starting point */
-/*     and calculate its norm. */
+    /*     evaluate the function at the starting point */
+    /*     and calculate its norm. */
 
     iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &wa3[1], 1);
     *nfev = 1;
     if (iflag < 0) {
-	goto L340;
+        goto L340;
     }
     fnorm = ei_enorm<Scalar>(m, &fvec[1]);
 
-/*     initialize levenberg-marquardt parameter and iteration counter. */
+    /*     initialize levenberg-marquardt parameter and iteration counter. */
 
     par = 0.;
     iter = 1;
 
-/*     beginning of the outer loop. */
+    /*     beginning of the outer loop. */
 
 L30:
 
-/*        if requested, call fcn to enable printing of iterates. */
+    /*        if requested, call fcn to enable printing of iterates. */
 
     if (nprint <= 0) {
-	goto L40;
+        goto L40;
     }
     iflag = 0;
     if ((iter - 1) % nprint == 0) {
-	iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &wa3[1], 0);
+        iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &wa3[1], 0);
     }
     if (iflag < 0) {
-	goto L340;
+        goto L340;
     }
 L40:
 
-/*        compute the qr factorization of the jacobian matrix */
-/*        calculated one row at a time, while simultaneously */
-/*        forming (q transpose)*fvec and storing the first */
-/*        n components in qtf. */
+    /*        compute the qr factorization of the jacobian matrix */
+    /*        calculated one row at a time, while simultaneously */
+    /*        forming (q transpose)*fvec and storing the first */
+    /*        n components in qtf. */
 
     for (j = 1; j <= n; ++j) {
-	qtf[j] = 0.;
-	for (i__ = 1; i__ <= n; ++i__) {
-	    fjac[i__ + j * fjac_dim1] = 0.;
-/* L50: */
-	}
-/* L60: */
+        qtf[j] = 0.;
+        for (i__ = 1; i__ <= n; ++i__) {
+            fjac[i__ + j * fjac_dim1] = 0.;
+            /* L50: */
+        }
+        /* L60: */
     }
     iflag = 2;
     for (i__ = 1; i__ <= m; ++i__) {
-	if ((*fcn)(p, m, n, &x[1], &fvec[1], &wa3[1], iflag) < 0) {
-	    goto L340;
-	}
-	temp = fvec[i__];
-	rwupdt(n, &fjac[fjac_offset], ldfjac, &wa3[1], &qtf[1], &temp, &wa1[
-		1], &wa2[1]);
-	++iflag;
-/* L70: */
+        if ((*fcn)(p, m, n, &x[1], &fvec[1], &wa3[1], iflag) < 0) {
+            goto L340;
+        }
+        temp = fvec[i__];
+        rwupdt(n, &fjac[fjac_offset], ldfjac, &wa3[1], &qtf[1], &temp, &wa1[
+                1], &wa2[1]);
+        ++iflag;
+        /* L70: */
     }
     ++(*njev);
 
-/*        if the jacobian is rank deficient, call qrfac to */
-/*        reorder its columns and update the components of qtf. */
+    /*        if the jacobian is rank deficient, call qrfac to */
+    /*        reorder its columns and update the components of qtf. */
 
     sing = FALSE_;
     for (j = 1; j <= n; ++j) {
-	if (fjac[j + j * fjac_dim1] == 0.) {
-	    sing = TRUE_;
-	}
-	ipvt[j] = j;
-	wa2[j] = ei_enorm<Scalar>(j, &fjac[j * fjac_dim1 + 1]);
-/* L80: */
+        if (fjac[j + j * fjac_dim1] == 0.) {
+            sing = TRUE_;
+        }
+        ipvt[j] = j;
+        wa2[j] = ei_enorm<Scalar>(j, &fjac[j * fjac_dim1 + 1]);
+        /* L80: */
     }
     if (! sing) {
-	goto L130;
+        goto L130;
     }
     qrfac(n, n, &fjac[fjac_offset], ldfjac, TRUE_, &ipvt[1], n, &wa1[1], &
-	    wa2[1], &wa3[1]);
+            wa2[1], &wa3[1]);
     for (j = 1; j <= n; ++j) {
-	if (fjac[j + j * fjac_dim1] == 0.) {
-	    goto L110;
-	}
-	sum = 0.;
-	for (i__ = j; i__ <= n; ++i__) {
-	    sum += fjac[i__ + j * fjac_dim1] * qtf[i__];
-/* L90: */
-	}
-	temp = -sum / fjac[j + j * fjac_dim1];
-	for (i__ = j; i__ <= n; ++i__) {
-	    qtf[i__] += fjac[i__ + j * fjac_dim1] * temp;
-/* L100: */
-	}
+        if (fjac[j + j * fjac_dim1] == 0.) {
+            goto L110;
+        }
+        sum = 0.;
+        for (i__ = j; i__ <= n; ++i__) {
+            sum += fjac[i__ + j * fjac_dim1] * qtf[i__];
+            /* L90: */
+        }
+        temp = -sum / fjac[j + j * fjac_dim1];
+        for (i__ = j; i__ <= n; ++i__) {
+            qtf[i__] += fjac[i__ + j * fjac_dim1] * temp;
+            /* L100: */
+        }
 L110:
-	fjac[j + j * fjac_dim1] = wa1[j];
-/* L120: */
+        fjac[j + j * fjac_dim1] = wa1[j];
+        /* L120: */
     }
 L130:
 
-/*        on the first iteration and if mode is 1, scale according */
-/*        to the norms of the columns of the initial jacobian. */
+    /*        on the first iteration and if mode is 1, scale according */
+    /*        to the norms of the columns of the initial jacobian. */
 
     if (iter != 1) {
-	goto L170;
+        goto L170;
     }
     if (mode == 2) {
-	goto L150;
+        goto L150;
     }
     for (j = 1; j <= n; ++j) {
-	diag[j] = wa2[j];
-	if (wa2[j] == 0.) {
-	    diag[j] = 1.;
-	}
-/* L140: */
+        diag[j] = wa2[j];
+        if (wa2[j] == 0.) {
+            diag[j] = 1.;
+        }
+        /* L140: */
     }
 L150:
 
-/*        on the first iteration, calculate the norm of the scaled x */
-/*        and initialize the step bound delta. */
+    /*        on the first iteration, calculate the norm of the scaled x */
+    /*        and initialize the step bound delta. */
 
     for (j = 1; j <= n; ++j) {
-	wa3[j] = diag[j] * x[j];
-/* L160: */
+        wa3[j] = diag[j] * x[j];
+        /* L160: */
     }
     xnorm = ei_enorm<Scalar>(n, &wa3[1]);
     delta = factor * xnorm;
     if (delta == 0.) {
-	delta = factor;
+        delta = factor;
     }
 L170:
 
-/*        compute the norm of the scaled gradient. */
+    /*        compute the norm of the scaled gradient. */
 
     gnorm = 0.;
     if (fnorm == 0.) {
-	goto L210;
+        goto L210;
     }
     for (j = 1; j <= n; ++j) {
-	l = ipvt[j];
-	if (wa2[l] == 0.) {
-	    goto L190;
-	}
-	sum = 0.;
-	for (i__ = 1; i__ <= j; ++i__) {
-	    sum += fjac[i__ + j * fjac_dim1] * (qtf[i__] / fnorm);
-/* L180: */
-	}
-/* Computing MAX */
-	d__2 = gnorm, d__3 = (d__1 = sum / wa2[l], abs(d__1));
-	gnorm = max(d__2,d__3);
+        l = ipvt[j];
+        if (wa2[l] == 0.) {
+            goto L190;
+        }
+        sum = 0.;
+        for (i__ = 1; i__ <= j; ++i__) {
+            sum += fjac[i__ + j * fjac_dim1] * (qtf[i__] / fnorm);
+            /* L180: */
+        }
+        /* Computing MAX */
+        d__2 = gnorm, d__3 = (d__1 = sum / wa2[l], abs(d__1));
+        gnorm = max(d__2,d__3);
 L190:
-/* L200: */
-	;
+        /* L200: */
+        ;
     }
 L210:
 
-/*        test for convergence of the gradient norm. */
+    /*        test for convergence of the gradient norm. */
 
     if (gnorm <= gtol) {
-	info = 4;
+        info = 4;
     }
     if (info != 0) {
-	goto L340;
+        goto L340;
     }
 
-/*        rescale if necessary. */
+    /*        rescale if necessary. */
 
     if (mode == 2) {
-	goto L230;
+        goto L230;
     }
     for (j = 1; j <= n; ++j) {
-/* Computing MAX */
-	d__1 = diag[j], d__2 = wa2[j];
-	diag[j] = max(d__1,d__2);
-/* L220: */
+        /* Computing MAX */
+        d__1 = diag[j], d__2 = wa2[j];
+        diag[j] = max(d__1,d__2);
+        /* L220: */
     }
 L230:
 
-/*        beginning of the inner loop. */
+    /*        beginning of the inner loop. */
 
 L240:
 
-/*           determine the levenberg-marquardt parameter. */
+    /*           determine the levenberg-marquardt parameter. */
 
     lmpar(n, &fjac[fjac_offset], ldfjac, &ipvt[1], &diag[1], &qtf[1], delta,
-	     &par, &wa1[1], &wa2[1], &wa3[1], &wa4[1]);
+            &par, &wa1[1], &wa2[1], &wa3[1], &wa4[1]);
 
-/*           store the direction p and x + p. calculate the norm of p. */
+    /*           store the direction p and x + p. calculate the norm of p. */
 
     for (j = 1; j <= n; ++j) {
-	wa1[j] = -wa1[j];
-	wa2[j] = x[j] + wa1[j];
-	wa3[j] = diag[j] * wa1[j];
-/* L250: */
+        wa1[j] = -wa1[j];
+        wa2[j] = x[j] + wa1[j];
+        wa3[j] = diag[j] * wa1[j];
+        /* L250: */
     }
     pnorm = ei_enorm<Scalar>(n, &wa3[1]);
 
-/*           on the first iteration, adjust the initial step bound. */
+    /*           on the first iteration, adjust the initial step bound. */
 
     if (iter == 1) {
-	delta = min(delta,pnorm);
+        delta = min(delta,pnorm);
     }
 
-/*           evaluate the function at x + p and calculate its norm. */
+    /*           evaluate the function at x + p and calculate its norm. */
 
     iflag = (*fcn)(p, m, n, &wa2[1], &wa4[1], &wa3[1], 1);
     ++(*nfev);
     if (iflag < 0) {
-	goto L340;
+        goto L340;
     }
     fnorm1 = ei_enorm<Scalar>(m, &wa4[1]);
 
-/*           compute the scaled actual reduction. */
+    /*           compute the scaled actual reduction. */
 
     actred = -1.;
     if (p1 * fnorm1 < fnorm) {
-/* Computing 2nd power */
-	d__1 = fnorm1 / fnorm;
-	actred = 1. - d__1 * d__1;
+        /* Computing 2nd power */
+        d__1 = fnorm1 / fnorm;
+        actred = 1. - d__1 * d__1;
     }
 
-/*           compute the scaled predicted reduction and */
-/*           the scaled directional derivative. */
+    /*           compute the scaled predicted reduction and */
+    /*           the scaled directional derivative. */
 
     for (j = 1; j <= n; ++j) {
-	wa3[j] = 0.;
-	l = ipvt[j];
-	temp = wa1[l];
-	for (i__ = 1; i__ <= j; ++i__) {
-	    wa3[i__] += fjac[i__ + j * fjac_dim1] * temp;
-/* L260: */
-	}
-/* L270: */
+        wa3[j] = 0.;
+        l = ipvt[j];
+        temp = wa1[l];
+        for (i__ = 1; i__ <= j; ++i__) {
+            wa3[i__] += fjac[i__ + j * fjac_dim1] * temp;
+            /* L260: */
+        }
+        /* L270: */
     }
     temp1 = ei_enorm<Scalar>(n, &wa3[1]) / fnorm;
     temp2 = sqrt(par) * pnorm / fnorm;
-/* Computing 2nd power */
+    /* Computing 2nd power */
     d__1 = temp1;
-/* Computing 2nd power */
+    /* Computing 2nd power */
     d__2 = temp2;
     prered = d__1 * d__1 + d__2 * d__2 / p5;
-/* Computing 2nd power */
+    /* Computing 2nd power */
     d__1 = temp1;
-/* Computing 2nd power */
+    /* Computing 2nd power */
     d__2 = temp2;
     dirder = -(d__1 * d__1 + d__2 * d__2);
 
-/*           compute the ratio of the actual to the predicted */
-/*           reduction. */
+    /*           compute the ratio of the actual to the predicted */
+    /*           reduction. */
 
     ratio = 0.;
     if (prered != 0.) {
-	ratio = actred / prered;
+        ratio = actred / prered;
     }
 
-/*           update the step bound. */
+    /*           update the step bound. */
 
     if (ratio > p25) {
-	goto L280;
+        goto L280;
     }
     if (actred >= 0.) {
-	temp = p5;
+        temp = p5;
     }
     if (actred < 0.) {
-	temp = p5 * dirder / (dirder + p5 * actred);
+        temp = p5 * dirder / (dirder + p5 * actred);
     }
     if (p1 * fnorm1 >= fnorm || temp < p1) {
-	temp = p1;
+        temp = p1;
     }
-/* Computing MIN */
+    /* Computing MIN */
     d__1 = delta, d__2 = pnorm / p1;
     delta = temp * min(d__1,d__2);
     par /= temp;
     goto L300;
 L280:
     if (par != 0. && ratio < p75) {
-	goto L290;
+        goto L290;
     }
     delta = pnorm / p5;
     par = p5 * par;
 L290:
 L300:
 
-/*           test for successful iteration. */
+    /*           test for successful iteration. */
 
     if (ratio < p0001) {
-	goto L330;
+        goto L330;
     }
 
-/*           successful iteration. update x, fvec, and their norms. */
+    /*           successful iteration. update x, fvec, and their norms. */
 
     for (j = 1; j <= n; ++j) {
-	x[j] = wa2[j];
-	wa2[j] = diag[j] * x[j];
-/* L310: */
+        x[j] = wa2[j];
+        wa2[j] = diag[j] * x[j];
+        /* L310: */
     }
     for (i__ = 1; i__ <= m; ++i__) {
-	fvec[i__] = wa4[i__];
-/* L320: */
+        fvec[i__] = wa4[i__];
+        /* L320: */
     }
     xnorm = ei_enorm<Scalar>(n, &wa2[1]);
     fnorm = fnorm1;
     ++iter;
 L330:
 
-/*           tests for convergence. */
+    /*           tests for convergence. */
 
     if (abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1.) {
-	info = 1;
+        info = 1;
     }
     if (delta <= xtol * xnorm) {
-	info = 2;
+        info = 2;
     }
     if (abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1. && info 
-	    == 2) {
-	info = 3;
+            == 2) {
+        info = 3;
     }
     if (info != 0) {
-	goto L340;
+        goto L340;
     }
 
-/*           tests for termination and stringent tolerances. */
+    /*           tests for termination and stringent tolerances. */
 
     if (*nfev >= maxfev) {
-	info = 5;
+        info = 5;
     }
     if (abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && p5 * ratio <= 1.) {
-	info = 6;
+        info = 6;
     }
     if (delta <= epsilon<Scalar>() * xnorm) {
-	info = 7;
+        info = 7;
     }
     if (gnorm <= epsilon<Scalar>()) {
-	info = 8;
+        info = 8;
     }
     if (info != 0) {
-	goto L340;
+        goto L340;
     }
 
-/*           end of the inner loop. repeat if iteration unsuccessful. */
+    /*           end of the inner loop. repeat if iteration unsuccessful. */
 
     if (ratio < p0001) {
-	goto L240;
+        goto L240;
     }
 
-/*        end of the outer loop. */
+    /*        end of the outer loop. */
 
     goto L30;
 L340:
 
-/*     termination, either normal or user imposed. */
+    /*     termination, either normal or user imposed. */
 
     if (iflag < 0) {
-	info = iflag;
+        info = iflag;
     }
     iflag = 0;
     if (nprint > 0) {
-	iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &wa3[1], 0);
+        iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &wa3[1], 0);
     }
     return info;
 
-/*     last card of subroutine lmstr. */
+    /*     last card of subroutine lmstr. */
 
 } /* lmstr_ */