dogleg, lmpar : use more eigen features

unsupported/Eigen/src/NonLinear/dogleg.h

@@ -47,16 +47,12 @@ void ei_dogleg(
             }
         }
         x[j] = (qtb[j] - sum) / temp;
-        /* L50: */
     }
 
     /*     test whether the gauss-newton direction is acceptable. */
 
-    for (j = 0; j < n; ++j) {
+    wa1.fill(0.);
-        wa1[j] = 0.;
+    wa2 = diag.cwise() * x;
-        wa2[j] = diag[j] * x[j];
-        /* L60: */
-    }
     qnorm = wa2.stableNorm();
     if (qnorm <= delta)
         return;
@@ -70,10 +66,8 @@ void ei_dogleg(
         for (i = j; i < n; ++i) {
             wa1[i] += r[l] * temp;
             ++l;
-            /* L70: */
         }
         wa1[j] /= diag[j];
-        /* L80: */
     }
 
     /*     calculate the norm of the scaled gradient and test for */
@@ -82,17 +76,13 @@ void ei_dogleg(
     gnorm = wa1.stableNorm();
     sgnorm = 0.;
     alpha = delta / qnorm;
-    if (gnorm == 0.) {
+    if (gnorm == 0.)
         goto L120;
-    }
 
     /*     calculate the point along the scaled gradient */
     /*     at which the quadratic is minimized. */
 
-    for (j = 0; j < n; ++j) {
+    wa1.cwise() /= diag*gnorm;
-        wa1[j] = wa1[j] / gnorm / diag[j];
-        /* L90: */
-    }
     l = 0;
     for (j = 0; j < n; ++j) {
         sum = 0.;
@@ -129,10 +119,8 @@ L120:
     /*     form appropriate convex combination of the gauss-newton */
     /*     direction and the scaled gradient direction. */
 
-    temp = (1. - alpha) * std::min(sgnorm,delta);
+    temp = (1.-alpha) * std::min(sgnorm,delta);
-    for (j = 0; j < n; ++j) {
+    x = temp * wa1 + alpha * x;
-        x[j] = temp * wa1[j] + alpha * x[j];
-    }
     return;
 
 }

unsupported/Eigen/src/NonLinear/lmpar.h

@@ -59,23 +59,19 @@ void ei_lmpar(
     /*     for acceptance of the gauss-newton direction. */
 
     iter = 0;
-    for (j = 0; j < n; ++j) {
+    wa2 = diag.cwise() * x;
-        wa2[j] = diag[j] * x[j];
-    }
     dxnorm = wa2.blueNorm();
     fp = dxnorm - delta;
-    if (fp <= Scalar(0.1) * delta) {
+    if (fp <= Scalar(0.1) * delta)
         goto L220;
-    }
 
     /*     if the jacobian is not rank deficient, the newton */
     /*     step provides a lower bound, parl, for the zero of */
     /*     the function. otherwise set this bound to zero. */
 
     parl = 0.;
-    if (nsing < n-1) {
+    if (nsing < n-1)
         goto L120;
-    }
     for (j = 0; j < n; ++j) {
         l = ipvt[j]-1;
         wa1[j] = diag[l] * (wa2[l] / dxnorm);
@@ -94,13 +90,10 @@ L120:
 
     for (j = 0; j < n; ++j) {
         sum = 0.;
-        for (i = 0; i <= j; ++i) {
+        for (i = 0; i <= j; ++i)
             sum += r(i,j) * qtb[i];
-            /* L130: */
-        }
         l = ipvt[j]-1;
         wa1[j] = sum / diag[l];
-        /* L140: */
     }
     gnorm = wa1.stableNorm();
     paru = gnorm / delta;
@@ -113,9 +106,8 @@ L120:
 
     par = std::max(par,parl);
     par = std::min(par,paru);
-    if (par == 0.) {
+    if (par == 0.)
         par = gnorm / dxnorm;
-    }
 
     /*     beginning of an iteration. */
 
@@ -124,20 +116,15 @@ L150:
 
     /*        evaluate the function at the current value of par. */
 
-    if (par == 0.) {
+    if (par == 0.)
-        /* Computing MAX */
+        par = std::max(dwarf,Scalar(.001) * paru); /* Computing MAX */
-        par = std::max(dwarf,Scalar(.001) * paru);
+
-    }
     temp = ei_sqrt(par);
-    for (j = 0; j < n; ++j) {
+    wa1 = temp * diag;
-        wa1[j] = temp * diag[j];
+
-        /* L160: */
-    }
     ei_qrsolv<Scalar>(n, r.data(), r.rows(), ipvt.data(), wa1.data(), qtb.data(), x.data(), sdiag.data(), wa2.data());
-    for (j = 0; j < n; ++j) {
+
-        wa2[j] = diag[j] * x[j];
+    wa2 = diag.cwise() * x;
-        /* L170: */
-    }
     dxnorm = wa2.blueNorm();
     temp = fp;
     fp = dxnorm - delta;