clean, remove goto's

2 files changed, 382 insertions, 488 deletions

diff --git a/unsupported/Eigen/src/NonLinear/hybrd.h b/unsupported/Eigen/src/NonLinear/hybrd.h
index 143f88832..4d0416917 100644
--- a/unsupported/Eigen/src/NonLinear/hybrd.h
+++ b/unsupported/Eigen/src/NonLinear/hybrd.h
@@ -19,7 +19,6 @@ int ei_hybrd(
         )
 {
     const int n = x.size();
-    int lr = (n*(n+1))/2;
     Matrix< Scalar, Dynamic, 1 > wa1(n), wa2(n), wa3(n), wa4(n);
 
 
@@ -27,7 +26,7 @@ int ei_hybrd(
     if (nb_of_superdiagonals<0) nb_of_superdiagonals = n-1;
     fvec.resize(n);
     qtf.resize(n);
-    R.resize(lr);
+    R.resize( (n*(n+1))/2);
     fjac.resize(n, n);
 
     /* Local variables */
@@ -56,10 +55,8 @@ int ei_hybrd(
 
     /*     check the input parameters for errors. */
 
-    if (n <= 0 || xtol < 0. || maxfev <= 0 || nb_of_subdiagonals < 0 || nb_of_superdiagonals < 0 ||
-            factor <= 0. || lr < n * (n + 1) / 2) {
+    if (n <= 0 || xtol < 0. || maxfev <= 0 || nb_of_subdiagonals < 0 || nb_of_superdiagonals < 0 || factor <= 0. )
         goto L300;
-    }
     if (mode == 2)
         for (j = 0; j < n; ++j)
             if (diag[j] <= 0.) goto L300;
@@ -69,9 +66,8 @@ int ei_hybrd(
 
     iflag = Functor::f(x, fvec);
     nfev = 1;
-    if (iflag < 0) {
+    if (iflag < 0)
         goto L300;
-    }
     fnorm = fvec.stableNorm();
 
     /*     determine the number of calls to fcn needed to compute */
@@ -90,287 +86,238 @@ int ei_hybrd(
 
     /*     beginning of the outer loop. */
 
-L30:
-    jeval = true;
+    while (true) {
+        jeval = true;
 
-    /*        calculate the jacobian matrix. */
+        /*        calculate the jacobian matrix. */
 
-    iflag = ei_fdjac1<Functor,Scalar>(x, fvec, fjac,
-            nb_of_subdiagonals, nb_of_superdiagonals, epsfcn, wa1, wa2);
-    nfev += msum;
-    if (iflag < 0) {
-        goto L300;
-    }
+        iflag = ei_fdjac1<Functor,Scalar>(x, fvec, fjac,
+                nb_of_subdiagonals, nb_of_superdiagonals, epsfcn, wa1, wa2);
+        nfev += msum;
+        if (iflag < 0)
+            break;
 
-    /*        compute the qr factorization of the jacobian. */
+        /*        compute the qr factorization of the jacobian. */
 
-    ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, 1, wa1.data(), wa2.data());
-
-    /*        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;
-    }
-    if (mode == 2) {
-        goto L50;
-    }
-    for (j = 0; j < n; ++j) {
-        diag[j] = wa2[j];
-        if (wa2[j] == 0.) {
-            diag[j] = 1.;
-        }
-        /* L40: */
-    }
-L50:
+        ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, 1, wa1.data(), wa2.data());
 
-    /*        on the first iteration, calculate the norm of the scaled x */
-    /*        and initialize the step bound delta. */
+        /*        on the first iteration and if mode is 1, scale according */
+        /*        to the norms of the columns of the initial jacobian. */
 
-    wa3 = diag.cwise() * x;
-    xnorm = wa3.stableNorm();
-    delta = factor * xnorm;
-    if (delta == 0.) {
-        delta = factor;
-    }
-L70:
+        if (iter == 1) {
+            if (mode != 2)
+                for (j = 0; j < n; ++j) {
+                    diag[j] = wa2[j];
+                    if (wa2[j] == 0.)
+                        diag[j] = 1.;
+                }
 
-    /*        form (q transpose)*fvec and store in qtf. */
+            /*        on the first iteration, calculate the norm of the scaled x */
+            /*        and initialize the step bound delta. */
 
-    qtf = fvec;
-    for (j = 0; j < n; ++j) {
-        if (fjac(j,j) == 0.) {
-            goto L110;
-        }
-        sum = 0.;
-        for (i = j; i < n; ++i) {
-            sum += fjac(i,j) * qtf[i];
-            /* L90: */
-        }
-        temp = -sum / fjac(j,j);
-        for (i = j; i < n; ++i) {
-            qtf[i] += fjac(i,j) * temp;
-            /* L100: */
+            wa3 = diag.cwise() * x;
+            xnorm = wa3.stableNorm();
+            delta = factor * xnorm;
+            if (delta == 0.)
+                delta = factor;
         }
-L110:
-        /* L120: */
-        ;
-    }
 
-    /*        copy the triangular factor of the qr factorization into r. */
+        /*        form (q transpose)*fvec and store in qtf. */
 
-    sing = false;
-    for (j = 0; j < n; ++j) {
-        l = j;
-        if (j) {
-            for (i = 0; i < j; ++i) {
-                R[l] = fjac(i,j);
-                l = l + n - i -1;
-                /* L130: */
+        qtf = fvec;
+        for (j = 0; j < n; ++j)
+            if (fjac(j,j) != 0.) {
+                sum = 0.;
+                for (i = j; i < n; ++i)
+                    sum += fjac(i,j) * qtf[i];
+                temp = -sum / fjac(j,j);
+                for (i = j; i < n; ++i)
+                    qtf[i] += fjac(i,j) * temp;
             }
-        }
-        R[l] = wa1[j];
-        if (wa1[j] == 0.) {
-            sing = true;
-        }
-        /* L150: */
-    }
-
-    /*        accumulate the orthogonal factor in fjac. */
-
-    ei_qform<Scalar>(n, n, fjac.data(), fjac.rows(), wa1.data());
-
-    /*        rescale if necessary. */
-
-    if (mode == 2) {
-        goto L170;
-    }
-    /* Computing MAX */
-    diag = diag.cwise().max(wa2);
-L170:
-
-    /*        beginning of the inner loop. */
-
-L180:
-    /*           if requested, call Functor::f to enable printing of iterates. */
-
-    if (nprint <= 0) {
-        goto L190;
-    }
-    iflag = 0;
-    if ((iter - 1) % nprint == 0)
-        iflag = Functor::debug(x, fvec);
-    if (iflag < 0)
-        goto L300;
-L190:
-
-    /*           determine the direction p. */
-
-    ei_dogleg<Scalar>(R, diag, qtf, delta, wa1);
-
-    /*           store the direction p and x + p. calculate the norm of p. */
 
-    wa1 = -wa1;
-    wa2 = x + wa1;
-    wa3 = diag.cwise() * wa1;
-    pnorm = wa3.stableNorm();
+        /*        copy the triangular factor of the qr factorization into r. */
+
+        sing = false;
+        for (j = 0; j < n; ++j) {
+            l = j;
+            if (j)
+                for (i = 0; i < j; ++i) {
+                    R[l] = fjac(i,j);
+                    l = l + n - i -1;
+                }
+            R[l] = wa1[j];
+            if (wa1[j] == 0.)
+                sing = true;
+        }
 
-    /*           on the first iteration, adjust the initial step bound. */
+        /*        accumulate the orthogonal factor in fjac. */
 
-    if (iter == 1)
-        delta = std::min(delta,pnorm);
+        ei_qform<Scalar>(n, n, fjac.data(), fjac.rows(), wa1.data());
 
-    /*           evaluate the function at x + p and calculate its norm. */
+        /*        rescale if necessary. */
 
-    iflag = Functor::f(wa2, wa4);
-    ++nfev;
-    if (iflag < 0)
-        goto L300;
-    fnorm1 = wa4.stableNorm();
+        /* Computing MAX */
+        if (mode != 2)
+            diag = diag.cwise().max(wa2);
 
-    /*           compute the scaled actual reduction. */
+        /*        beginning of the inner loop. */
 
-    actred = -1.;
-    if (fnorm1 < fnorm) /* Computing 2nd power */
-        actred = 1. - ei_abs2(fnorm1 / fnorm);
+        while (true) {
+            /*           if requested, call Functor::f to enable printing of iterates. */
 
-    /*           compute the scaled predicted reduction. */
+            if (nprint > 0) {
+                iflag = 0;
+                if ((iter - 1) % nprint == 0)
+                    iflag = Functor::debug(x, fvec);
+                if (iflag < 0)
+                    goto L300;
+            }
 
-    l = 0;
-    for (i = 0; i < n; ++i) {
-        sum = 0.;
-        for (j = i; j < n; ++j) {
-            sum += R[l] * wa1[j];
-            ++l;
-            /* L210: */
-        }
-        wa3[i] = qtf[i] + sum;
-        /* L220: */
-    }
-    temp = wa3.stableNorm();
-    prered = 0.;
-    if (temp < fnorm) /* Computing 2nd power */
-        prered = 1. - ei_abs2(temp / fnorm);
+            /*           determine the direction p. */
 
-    /*           compute the ratio of the actual to the predicted */
-    /*           reduction. */
+            ei_dogleg<Scalar>(R, diag, qtf, delta, wa1);
 
-    ratio = 0.;
-    if (prered > 0.) {
-        ratio = actred / prered;
-    }
+            /*           store the direction p and x + p. calculate the norm of p. */
 
-    /*           update the step bound. */
+            wa1 = -wa1;
+            wa2 = x + wa1;
+            wa3 = diag.cwise() * wa1;
+            pnorm = wa3.stableNorm();
 
-    if (ratio >= Scalar(.1)) {
-        goto L230;
-    }
-    ncsuc = 0;
-    ++ncfail;
-    delta = Scalar(.5) * delta;
-    goto L240;
-L230:
-    ncfail = 0;
-    ++ncsuc;
-    if (ratio >= Scalar(.5) || ncsuc > 1) /* Computing MAX */
-        delta = std::max(delta, pnorm / Scalar(.5));
-    if (ei_abs(ratio - 1.) <= Scalar(.1)) {
-        delta = pnorm / Scalar(.5);
-    }
-L240:
+            /*           on the first iteration, adjust the initial step bound. */
 
-    /*           test for successful iteration. */
+            if (iter == 1)
+                delta = std::min(delta,pnorm);
 
-    if (ratio < Scalar(1e-4)) {
-        goto L260;
-    }
+            /*           evaluate the function at x + p and calculate its norm. */
 
-    /*           successful iteration. update x, fvec, and their norms. */
+            iflag = Functor::f(wa2, wa4);
+            ++nfev;
+            if (iflag < 0)
+                goto L300;
+            fnorm1 = wa4.stableNorm();
 
-    x = wa2;
-    wa2 = diag.cwise() * x;
-    fvec = wa4;
-    xnorm = wa2.stableNorm();
-    fnorm = fnorm1;
-    ++iter;
-L260:
+            /*           compute the scaled actual reduction. */
 
-    /*           determine the progress of the iteration. */
+            actred = -1.;
+            if (fnorm1 < fnorm) /* Computing 2nd power */
+                actred = 1. - ei_abs2(fnorm1 / fnorm);
 
-    ++nslow1;
-    if (actred >= Scalar(.001)) {
-        nslow1 = 0;
-    }
-    if (jeval) {
-        ++nslow2;
-    }
-    if (actred >= Scalar(.1)) {
-        nslow2 = 0;
-    }
+            /*           compute the scaled predicted reduction. */
 
-    /*           test for convergence. */
+            l = 0;
+            for (i = 0; i < n; ++i) {
+                sum = 0.;
+                for (j = i; j < n; ++j) {
+                    sum += R[l] * wa1[j];
+                    ++l;
+                }
+                wa3[i] = qtf[i] + sum;
+            }
+            temp = wa3.stableNorm();
+            prered = 0.;
+            if (temp < fnorm) /* Computing 2nd power */
+                prered = 1. - ei_abs2(temp / fnorm);
+
+            /*           compute the ratio of the actual to the predicted */
+            /*           reduction. */
+
+            ratio = 0.;
+            if (prered > 0.)
+                ratio = actred / prered;
+
+            /*           update the step bound. */
+
+            if (ratio < Scalar(.1)) {
+                ncsuc = 0;
+                ++ncfail;
+                delta = Scalar(.5) * delta;
+            } else {
+                ncfail = 0;
+                ++ncsuc;
+                if (ratio >= Scalar(.5) || ncsuc > 1) /* Computing MAX */
+                    delta = std::max(delta, pnorm / Scalar(.5));
+                if (ei_abs(ratio - 1.) <= Scalar(.1)) {
+                    delta = pnorm / Scalar(.5);
+                }
+            }
 
-    if (delta <= xtol * xnorm || fnorm == 0.) {
-        info = 1;
-    }
-    if (info != 0) {
-        goto L300;
-    }
+            /*           test for successful iteration. */
 
-    /*           tests for termination and stringent tolerances. */
+            if (ratio >= Scalar(1e-4)) {
+                /*           successful iteration. update x, fvec, and their norms. */
+                x = wa2;
+                wa2 = diag.cwise() * x;
+                fvec = wa4;
+                xnorm = wa2.stableNorm();
+                fnorm = fnorm1;
+                ++iter;
+            }
 
-    if (nfev >= maxfev) {
-        info = 2;
-    }
-    /* Computing MAX */
-    if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm)
-        info = 3;
-    if (nslow2 == 5)
-        info = 4;
-    if (nslow1 == 10)
-        info = 5;
-    if (info != 0)
-        goto L300;
+            /*           determine the progress of the iteration. */
+
+            ++nslow1;
+            if (actred >= Scalar(.001))
+                nslow1 = 0;
+            if (jeval)
+                ++nslow2;
+            if (actred >= Scalar(.1))
+                nslow2 = 0;
+
+            /*           test for convergence. */
+
+            if (delta <= xtol * xnorm || fnorm == 0.)
+                info = 1;
+            if (info != 0)
+                goto L300;
+
+            /*           tests for termination and stringent tolerances. */
+
+            if (nfev >= maxfev)
+                info = 2;
+            /* Computing MAX */
+            if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm)
+                info = 3;
+            if (nslow2 == 5)
+                info = 4;
+            if (nslow1 == 10)
+                info = 5;
+            if (info != 0)
+                goto L300;
+
+            /*           criterion for recalculating jacobian approximation */
+            /*           by forward differences. */
+
+            if (ncfail == 2)
+                break;
+
+            /*           calculate the rank one modification to the jacobian */
+            /*           and update qtf if necessary. */
+
+            for (j = 0; j < n; ++j) {
+                sum = wa4.dot(fjac.col(j));
+                wa2[j] = (sum - wa3[j]) / pnorm;
+                wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
+                if (ratio >= Scalar(1e-4))
+                    qtf[j] = sum;
+            }
 
-    /*           criterion for recalculating jacobian approximation */
-    /*           by forward differences. */
+            /*           compute the qr factorization of the updated jacobian. */
 
-    if (ncfail == 2)
-        goto L290;
+            ei_r1updt<Scalar>(n, n, R.data(), R.size(), wa1.data(), wa2.data(), wa3.data(), &sing);
+            ei_r1mpyq<Scalar>(n, n, fjac.data(), fjac.rows(), wa2.data(), wa3.data());
+            ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data());
 
-    /*           calculate the rank one modification to the jacobian */
-    /*           and update qtf if necessary. */
+            /*           end of the inner loop. */
 
-    for (j = 0; j < n; ++j) {
-        sum = wa4.dot(fjac.col(j));
-        wa2[j] = (sum - wa3[j]) / pnorm;
-        wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
-        if (ratio >= Scalar(1e-4))
-            qtf[j] = sum;
+            jeval = false;
+        }
+        /*        end of the outer loop. */
     }
-
-    /*           compute the qr factorization of the updated jacobian. */
-
-    ei_r1updt<Scalar>(n, n, R.data(), lr, wa1.data(), wa2.data(), wa3.data(), &sing);
-    ei_r1mpyq<Scalar>(n, n, fjac.data(), fjac.rows(), wa2.data(), wa3.data());
-    ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data());
-
-    /*           end of the inner loop. */
-
-    jeval = false;
-    goto L180;
-L290:
-
-    /*        end of the outer loop. */
-
-    goto L30;
 L300:
-
     /*     termination, either normal or user imposed. */
-
-    if (iflag < 0) {
+    if (iflag < 0)
         info = iflag;
-    }
     if (nprint > 0)
         iflag = Functor::debug(x, fvec);
     return info;
diff --git a/unsupported/Eigen/src/NonLinear/hybrj.h b/unsupported/Eigen/src/NonLinear/hybrj.h
index 8f78c1868..4276d804b 100644
--- a/unsupported/Eigen/src/NonLinear/hybrj.h
+++ b/unsupported/Eigen/src/NonLinear/hybrj.h
@@ -17,12 +17,11 @@ int ei_hybrj(
         )
 {
     const int n = x.size();
-    const int lr = (n*(n+1))/2;
     Matrix< Scalar, Dynamic, 1 > wa1(n), wa2(n), wa3(n), wa4(n);
 
     fvec.resize(n);
     qtf.resize(n);
-    R.resize(lr);
+    R.resize( (n*(n+1))/2);
     fjac.resize(n, n);
 
     /* Local variables */
@@ -51,10 +50,8 @@ int ei_hybrj(
 
     /*     check the input parameters for errors. */
 
-    if (n <= 0 || xtol < 0. || maxfev <= 0 || factor <= 
-            0. || lr < n * (n + 1) / 2) {
+    if (n <= 0 || xtol < 0. || maxfev <= 0 || factor <= 0. )
         goto L300;
-    }
     if (mode == 2)
         for (j = 0; j < n; ++j)
             if (diag[j] <= 0.) goto L300;
@@ -64,9 +61,8 @@ int ei_hybrj(
 
     iflag = Functor::f(x, fvec);
     nfev = 1;
-    if (iflag < 0) {
+    if (iflag < 0)
         goto L300;
-    }
     fnorm = fvec.stableNorm();
 
     /*     initialize iteration counter and monitors. */
@@ -79,285 +75,236 @@ int ei_hybrj(
 
     /*     beginning of the outer loop. */
 
-L30:
-    jeval = true;
+    while (true) {
+        jeval = true;
 
-    /*        calculate the jacobian matrix. */
+        /*        calculate the jacobian matrix. */
 
-    iflag = Functor::df(x, fjac);
-    ++njev;
-    if (iflag < 0) {
-        goto L300;
-    }
+        iflag = Functor::df(x, fjac);
+        ++njev;
+        if (iflag < 0)
+            break;
 
-    /*        compute the qr factorization of the jacobian. */
+        /*        compute the qr factorization of the jacobian. */
 
-    ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, 1, wa1.data(), wa2.data());
-
-    /*        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;
-    }
-    if (mode == 2) {
-        goto L50;
-    }
-    for (j = 0; j < n; ++j) {
-        diag[j] = wa2[j];
-        if (wa2[j] == 0.) {
-            diag[j] = 1.;
-        }
-        /* L40: */
-    }
-L50:
+        ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, 1, wa1.data(), wa2.data());
 
-    /*        on the first iteration, calculate the norm of the scaled x */
-    /*        and initialize the step bound delta. */
+        /*        on the first iteration and if mode is 1, scale according */
+        /*        to the norms of the columns of the initial jacobian. */
 
-    wa3 = diag.cwise() * x;
-    xnorm = wa3.stableNorm();
-    delta = factor * xnorm;
-    if (delta == 0.) {
-        delta = factor;
-    }
-L70:
+        if (iter == 1) {
+            if (mode != 2)
+                for (j = 0; j < n; ++j) {
+                    diag[j] = wa2[j];
+                    if (wa2[j] == 0.)
+                        diag[j] = 1.;
+                }
 
-    /*        form (q transpose)*fvec and store in qtf. */
+            /*        on the first iteration, calculate the norm of the scaled x */
+            /*        and initialize the step bound delta. */
 
-    qtf = fvec;
-    for (j = 0; j < n; ++j) {
-        if (fjac(j,j) == 0.) {
-            goto L110;
-        }
-        sum = 0.;
-        for (i = j; i < n; ++i) {
-            sum += fjac(i,j) * qtf[i];
-            /* L90: */
-        }
-        temp = -sum / fjac(j,j);
-        for (i = j; i < n; ++i) {
-            qtf[i] += fjac(i,j) * temp;
-            /* L100: */
+            wa3 = diag.cwise() * x;
+            xnorm = wa3.stableNorm();
+            delta = factor * xnorm;
+            if (delta == 0.)
+                delta = factor;
         }
-L110:
-        /* L120: */
-        ;
-    }
 
-    /*        copy the triangular factor of the qr factorization into r. */
+        /*        form (q transpose)*fvec and store in qtf. */
 
-    sing = false;
-    for (j = 0; j < n; ++j) {
-        l = j;
-        if (j) {
-            for (i = 0; i < j; ++i) {
-                R[l] = fjac(i,j);
-                l = l + n - i -1;
-                /* L130: */
+        qtf = fvec;
+        for (j = 0; j < n; ++j)
+            if (fjac(j,j) != 0.) {
+                sum = 0.;
+                for (i = j; i < n; ++i)
+                    sum += fjac(i,j) * qtf[i];
+                temp = -sum / fjac(j,j);
+                for (i = j; i < n; ++i)
+                    qtf[i] += fjac(i,j) * temp;
             }
-        }
-        R[l] = wa1[j];
-        if (wa1[j] == 0.) {
-            sing = true;
-        }
-        /* L150: */
-    }
-
-    /*        accumulate the orthogonal factor in fjac. */
-
-    ei_qform<Scalar>(n, n, fjac.data(), fjac.rows(), wa1.data());
-
-    /*        rescale if necessary. */
-
-    if (mode == 2) {
-        goto L170;
-    }
-    /* Computing MAX */
-    diag = diag.cwise().max(wa2);
-L170:
-
-    /*        beginning of the inner loop. */
-
-L180:
-    /*           if requested, call Functor::f to enable printing of iterates. */
-
-    if (nprint <= 0) {
-        goto L190;
-    }
-    iflag = 0;
-    if ((iter - 1) % nprint == 0)
-        iflag = Functor::debug(x, fvec, fjac);
-    if (iflag < 0)
-        goto L300;
-L190:
-
-    /*           determine the direction p. */
-
-    ei_dogleg<Scalar>(R, diag, qtf, delta, wa1);
-
-    /*           store the direction p and x + p. calculate the norm of p. */
 
-    wa1 = -wa1;
-    wa2 = x + wa1;
-    wa3 = diag.cwise() * wa1;
-    pnorm = wa3.stableNorm();
+        /*        copy the triangular factor of the qr factorization into r. */
+
+        sing = false;
+        for (j = 0; j < n; ++j) {
+            l = j;
+            if (j)
+                for (i = 0; i < j; ++i) {
+                    R[l] = fjac(i,j);
+                    l = l + n - i -1;
+                }
+            R[l] = wa1[j];
+            if (wa1[j] == 0.)
+                sing = true;
+        }
 
-    /*           on the first iteration, adjust the initial step bound. */
+        /*        accumulate the orthogonal factor in fjac. */
 
-    if (iter == 1)
-        delta = std::min(delta,pnorm);
+        ei_qform<Scalar>(n, n, fjac.data(), fjac.rows(), wa1.data());
 
-    /*           evaluate the function at x + p and calculate its norm. */
+        /*        rescale if necessary. */
 
-    iflag = Functor::f(wa2, wa4);
-    ++nfev;
-    if (iflag < 0)
-        goto L300;
-    fnorm1 = wa4.stableNorm();
+        /* Computing MAX */
+        if (mode != 2)
+            diag = diag.cwise().max(wa2);
 
-    /*           compute the scaled actual reduction. */
+        /*        beginning of the inner loop. */
 
-    actred = -1.;
-    if (fnorm1 < fnorm) /* Computing 2nd power */
-        actred = 1. - ei_abs2(fnorm1 / fnorm);
+        while (true) {
+            /*           if requested, call Functor::f to enable printing of iterates. */
 
-    /*           compute the scaled predicted reduction. */
+            if (nprint > 0) {
+                iflag = 0;
+                if ((iter - 1) % nprint == 0)
+                    iflag = Functor::debug(x, fvec, fjac);
+                if (iflag < 0)
+                    goto L300;
+            }
 
-    l = 0;
-    for (i = 0; i < n; ++i) {
-        sum = 0.;
-        for (j = i; j < n; ++j) {
-            sum += R[l] * wa1[j];
-            ++l;
-            /* L210: */
-        }
-        wa3[i] = qtf[i] + sum;
-        /* L220: */
-    }
-    temp = wa3.stableNorm();
-    prered = 0.;
-    if (temp < fnorm) /* Computing 2nd power */
-        prered = 1. - ei_abs2(temp / fnorm);
+            /*           determine the direction p. */
 
-    /*           compute the ratio of the actual to the predicted */
-    /*           reduction. */
+            ei_dogleg<Scalar>(R, diag, qtf, delta, wa1);
 
-    ratio = 0.;
-    if (prered > 0.) {
-        ratio = actred / prered;
-    }
+            /*           store the direction p and x + p. calculate the norm of p. */
 
-    /*           update the step bound. */
+            wa1 = -wa1;
+            wa2 = x + wa1;
+            wa3 = diag.cwise() * wa1;
+            pnorm = wa3.stableNorm();
 
-    if (ratio >= Scalar(.1)) {
-        goto L230;
-    }
-    ncsuc = 0;
-    ++ncfail;
-    delta = Scalar(.5) * delta;
-    goto L240;
-L230:
-    ncfail = 0;
-    ++ncsuc;
-    if (ratio >= Scalar(.5) || ncsuc > 1) /* Computing MAX */
-        delta = std::max(delta, pnorm / Scalar(.5));
-    if (ei_abs(ratio - 1.) <= Scalar(.1)) {
-        delta = pnorm / Scalar(.5);
-    }
-L240:
+            /*           on the first iteration, adjust the initial step bound. */
 
-    /*           test for successful iteration. */
+            if (iter == 1)
+                delta = std::min(delta,pnorm);
 
-    if (ratio < Scalar(1e-4)) {
-        goto L260;
-    }
+            /*           evaluate the function at x + p and calculate its norm. */
 
-    /*           successful iteration. update x, fvec, and their norms. */
+            iflag = Functor::f(wa2, wa4);
+            ++nfev;
+            if (iflag < 0)
+                goto L300;
+            fnorm1 = wa4.stableNorm();
 
-    x = wa2;
-    wa2 = diag.cwise() * x;
-    fvec = wa4;
-    xnorm = wa2.stableNorm();
-    fnorm = fnorm1;
-    ++iter;
-L260:
+            /*           compute the scaled actual reduction. */
 
-    /*           determine the progress of the iteration. */
+            actred = -1.;
+            if (fnorm1 < fnorm) /* Computing 2nd power */
+                actred = 1. - ei_abs2(fnorm1 / fnorm);
 
-    ++nslow1;
-    if (actred >= Scalar(.001)) {
-        nslow1 = 0;
-    }
-    if (jeval) {
-        ++nslow2;
-    }
-    if (actred >= Scalar(.1)) {
-        nslow2 = 0;
-    }
+            /*           compute the scaled predicted reduction. */
 
-    /*           test for convergence. */
+            l = 0;
+            for (i = 0; i < n; ++i) {
+                sum = 0.;
+                for (j = i; j < n; ++j) {
+                    sum += R[l] * wa1[j];
+                    ++l;
+                }
+                wa3[i] = qtf[i] + sum;
+            }
+            temp = wa3.stableNorm();
+            prered = 0.;
+            if (temp < fnorm) /* Computing 2nd power */
+                prered = 1. - ei_abs2(temp / fnorm);
+
+            /*           compute the ratio of the actual to the predicted */
+            /*           reduction. */
+
+            ratio = 0.;
+            if (prered > 0.)
+                ratio = actred / prered;
+
+            /*           update the step bound. */
+
+            if (ratio < Scalar(.1)) {
+                ncsuc = 0;
+                ++ncfail;
+                delta = Scalar(.5) * delta;
+            } else {
+                ncfail = 0;
+                ++ncsuc;
+                if (ratio >= Scalar(.5) || ncsuc > 1) /* Computing MAX */
+                    delta = std::max(delta, pnorm / Scalar(.5));
+                if (ei_abs(ratio - 1.) <= Scalar(.1)) {
+                    delta = pnorm / Scalar(.5);
+                }
+            }
 
-    if (delta <= xtol * xnorm || fnorm == 0.) {
-        info = 1;
-    }
-    if (info != 0) {
-        goto L300;
-    }
+            /*           test for successful iteration. */
 
-    /*           tests for termination and stringent tolerances. */
+            if (ratio >= Scalar(1e-4)) {
+                /*           successful iteration. update x, fvec, and their norms. */
+                x = wa2;
+                wa2 = diag.cwise() * x;
+                fvec = wa4;
+                xnorm = wa2.stableNorm();
+                fnorm = fnorm1;
+                ++iter;
+            }
 
-    if (nfev >= maxfev) {
-        info = 2;
-    }
-    /* Computing MAX */
-    if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm)
-        info = 3;
-    if (nslow2 == 5)
-        info = 4;
-    if (nslow1 == 10)
-        info = 5;
-    if (info != 0)
-        goto L300;
+            /*           determine the progress of the iteration. */
+
+            ++nslow1;
+            if (actred >= Scalar(.001))
+                nslow1 = 0;
+            if (jeval)
+                ++nslow2;
+            if (actred >= Scalar(.1))
+                nslow2 = 0;
+
+            /*           test for convergence. */
+
+            if (delta <= xtol * xnorm || fnorm == 0.)
+                info = 1;
+            if (info != 0)
+                goto L300;
+
+            /*           tests for termination and stringent tolerances. */
+
+            if (nfev >= maxfev)
+                info = 2;
+            /* Computing MAX */
+            if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm)
+                info = 3;
+            if (nslow2 == 5)
+                info = 4;
+            if (nslow1 == 10)
+                info = 5;
+            if (info != 0)
+                goto L300;
+
+            /*           criterion for recalculating jacobian. */
+
+            if (ncfail == 2)
+                break;
+
+            /*           calculate the rank one modification to the jacobian */
+            /*           and update qtf if necessary. */
+
+            for (j = 0; j < n; ++j) {
+                sum = wa4.dot(fjac.col(j));
+                wa2[j] = (sum - wa3[j]) / pnorm;
+                wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
+                if (ratio >= Scalar(1e-4))
+                    qtf[j] = sum;
+            }
 
-    /*           criterion for recalculating jacobian. */
+            /*           compute the qr factorization of the updated jacobian. */
 
-    if (ncfail == 2)
-        goto L290;
+            ei_r1updt<Scalar>(n, n, R.data(), R.size(), wa1.data(), wa2.data(), wa3.data(), &sing);
+            ei_r1mpyq<Scalar>(n, n, fjac.data(), fjac.rows(), wa2.data(), wa3.data());
+            ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data());
 
-    /*           calculate the rank one modification to the jacobian */
-    /*           and update qtf if necessary. */
+            /*           end of the inner loop. */
 
-    for (j = 0; j < n; ++j) {
-        sum = wa4.dot(fjac.col(j));
-        wa2[j] = (sum - wa3[j]) / pnorm;
-        wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
-        if (ratio >= Scalar(1e-4))
-            qtf[j] = sum;
+            jeval = false;
+        }
+        /*        end of the outer loop. */
     }
-
-    /*           compute the qr factorization of the updated jacobian. */
-
-    ei_r1updt<Scalar>(n, n, R.data(), lr, wa1.data(), wa2.data(), wa3.data(), &sing);
-    ei_r1mpyq<Scalar>(n, n, fjac.data(), fjac.rows(), wa2.data(), wa3.data());
-    ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data());
-
-    /*           end of the inner loop. */
-
-    jeval = false;
-    goto L180;
-L290:
-
-    /*        end of the outer loop. */
-
-    goto L30;
 L300:
-
     /*     termination, either normal or user imposed. */
-
-    if (iflag < 0) {
+    if (iflag < 0)
         info = iflag;
-    }
     if (nprint > 0)
         iflag = Functor::debug(x, fvec, fjac);
     return info;