diff options
Diffstat (limited to 'unsupported/Eigen')
-rw-r--r-- | unsupported/Eigen/src/NonLinear/dogleg.h | 24 | ||||
-rw-r--r-- | unsupported/Eigen/src/NonLinear/lmpar.h | 37 |
2 files changed, 18 insertions, 43 deletions
diff --git a/unsupported/Eigen/src/NonLinear/dogleg.h b/unsupported/Eigen/src/NonLinear/dogleg.h index ea475041f..cc1283847 100644 --- a/unsupported/Eigen/src/NonLinear/dogleg.h +++ b/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[j] = 0.; - wa2[j] = diag[j] * x[j]; - /* L60: */ - } + wa1.fill(0.); + wa2 = diag.cwise() * x; 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[j] = wa1[j] / gnorm / diag[j]; - /* L90: */ - } + wa1.cwise() /= diag*gnorm; 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); - for (j = 0; j < n; ++j) { - x[j] = temp * wa1[j] + alpha * x[j]; - } + temp = (1.-alpha) * std::min(sgnorm,delta); + x = temp * wa1 + alpha * x; return; } diff --git a/unsupported/Eigen/src/NonLinear/lmpar.h b/unsupported/Eigen/src/NonLinear/lmpar.h index 663b0cd7b..28fcf390b 100644 --- a/unsupported/Eigen/src/NonLinear/lmpar.h +++ b/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[j] = diag[j] * x[j]; - } + wa2 = diag.cwise() * x; 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.) { - /* Computing MAX */ - par = std::max(dwarf,Scalar(.001) * paru); - } + if (par == 0.) + par = std::max(dwarf,Scalar(.001) * paru); /* Computing MAX */ + temp = ei_sqrt(par); - for (j = 0; j < n; ++j) { - wa1[j] = temp * diag[j]; - /* L160: */ - } + wa1 = temp * diag; + 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]; - /* L170: */ - } + + wa2 = diag.cwise() * x; dxnorm = wa2.blueNorm(); temp = fp; fp = dxnorm - delta; |