GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-13 04:09:10 +08:00
Code Issues Packages Projects Releases Wiki Activity

further cleaning/ goto removing

Browse Source
This commit is contained in:
Thomas Capricelli 2009-08-24 16:39:49 +02:00
parent 92a5bb4539
commit 312ab1abb3
2 changed files with 101 additions and 110 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
unsupported/Eigen/src/NonLinear/dogleg.h
View File

@ -20,7 +20,7 @@ void ei_dogleg(
assert(n==qtb.size()); assert(n==qtb.size());
assert(n==x.size()); assert(n==x.size());
/* first, calculate the gauss-newton direction. */ /* first, calculate the gauss-newton direction. */
jj = n * (n + 1) / 2; jj = n * (n + 1) / 2;
for (k = 0; k < n; ++k) { for (k = 0; k < n; ++k) {
@ -39,17 +39,15 @@ void ei_dogleg(
/* Computing MAX */ /* Computing MAX */
temp = std::max(temp,ei_abs(r[l])); temp = std::max(temp,ei_abs(r[l]));
l = l + n - i; l = l + n - i;
/* L30: */
} }
temp = epsmch * temp; temp = epsmch * temp;
if (temp == 0.) { if (temp == 0.)
temp = epsmch; temp = epsmch;
}
} }
x[j] = (qtb[j] - sum) / temp; x[j] = (qtb[j] - sum) / temp;
} }
/* test whether the gauss-newton direction is acceptable. */ /* test whether the gauss-newton direction is acceptable. */
wa1.fill(0.); wa1.fill(0.);
wa2 = diag.cwise() * x; wa2 = diag.cwise() * x;
@ -57,8 +55,8 @@ void ei_dogleg(
if (qnorm <= delta) if (qnorm <= delta)
return; return;
/* the gauss-newton direction is not acceptable. */ /* the gauss-newton direction is not acceptable. */
/* next, calculate the scaled gradient direction. */ /* next, calculate the scaled gradient direction. */
l = 0; l = 0;
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
@ -70,17 +68,17 @@ void ei_dogleg(
wa1[j] /= diag[j]; wa1[j] /= diag[j];
} }
/* calculate the norm of the scaled gradient and test for */ /* calculate the norm of the scaled gradient and test for */
/* the special case in which the scaled gradient is zero. */ /* the special case in which the scaled gradient is zero. */
gnorm = wa1.stableNorm(); gnorm = wa1.stableNorm();
sgnorm = 0.; sgnorm = 0.;
alpha = delta / qnorm; alpha = delta / qnorm;
if (gnorm == 0.) if (gnorm == 0.)
goto L120; goto algo_end;
/* calculate the point along the scaled gradient */ /* calculate the point along the scaled gradient */
/* at which the quadratic is minimized. */ /* at which the quadratic is minimized. */
wa1.cwise() /= diag*gnorm; wa1.cwise() /= diag*gnorm;
l = 0; l = 0;
@ -97,16 +95,15 @@ void ei_dogleg(
temp = wa2.stableNorm(); temp = wa2.stableNorm();
sgnorm = gnorm / temp / temp; sgnorm = gnorm / temp / temp;
/* test whether the scaled gradient direction is acceptable. */ /* test whether the scaled gradient direction is acceptable. */
alpha = 0.; alpha = 0.;
if (sgnorm >= delta) { if (sgnorm >= delta)
goto L120; goto algo_end;
}
/* the scaled gradient direction is not acceptable. */ /* the scaled gradient direction is not acceptable. */
/* finally, calculate the point along the dogleg */ /* finally, calculate the point along the dogleg */
/* at which the quadratic is minimized. */ /* at which the quadratic is minimized. */
bnorm = qtb.stableNorm(); bnorm = qtb.stableNorm();
temp = bnorm / gnorm * (bnorm / qnorm) * (sgnorm / delta); temp = bnorm / gnorm * (bnorm / qnorm) * (sgnorm / delta);
@ -114,10 +111,10 @@ void ei_dogleg(
temp = temp - delta / qnorm * ei_abs2(sgnorm / delta) + ei_sqrt(ei_abs2(temp - delta / qnorm) + (1.-ei_abs2(delta / qnorm)) * (1.-ei_abs2(sgnorm / delta))); temp = temp - delta / qnorm * ei_abs2(sgnorm / delta) + ei_sqrt(ei_abs2(temp - delta / qnorm) + (1.-ei_abs2(delta / qnorm)) * (1.-ei_abs2(sgnorm / delta)));
/* Computing 2nd power */ /* Computing 2nd power */
alpha = delta / qnorm * (1. - ei_abs2(sgnorm / delta)) / temp; alpha = delta / qnorm * (1. - ei_abs2(sgnorm / delta)) / temp;
L120: algo_end:
/* form appropriate convex combination of the gauss-newton */ /* form appropriate convex combination of the gauss-newton */
/* direction and the scaled gradient direction. */ /* direction and the scaled gradient direction. */
temp = (1.-alpha) * std::min(sgnorm,delta); temp = (1.-alpha) * std::min(sgnorm,delta);
x = temp * wa1 + alpha * x; x = temp * wa1 + alpha * x;

170
unsupported/Eigen/src/NonLinear/lmpar.h
View File

@ -30,8 +30,8 @@ void ei_lmpar(
Matrix< Scalar, Dynamic, 1 > wa1(n), wa2(n); Matrix< Scalar, Dynamic, 1 > wa1(n), wa2(n);
/* compute and store in x the gauss-newton direction. if the */ /* compute and store in x the gauss-newton direction. if the */
/* jacobian is rank-deficient, obtain a least squares solution. */ /* jacobian is rank-deficient, obtain a least squares solution. */
nsing = n-1; nsing = n-1;
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
@ -54,39 +54,40 @@ void ei_lmpar(
x[l] = wa1[j]; x[l] = wa1[j];
} }
/* initialize the iteration counter. */ /* initialize the iteration counter. */
/* evaluate the function at the origin, and test */ /* evaluate the function at the origin, and test */
/* for acceptance of the gauss-newton direction. */ /* for acceptance of the gauss-newton direction. */
iter = 0; iter = 0;
wa2 = diag.cwise() * x; wa2 = diag.cwise() * x;
dxnorm = wa2.blueNorm(); dxnorm = wa2.blueNorm();
fp = dxnorm - delta; fp = dxnorm - delta;
if (fp <= Scalar(0.1) * delta) if (fp <= Scalar(0.1) * delta) {
goto L220; par = 0;
return;
}
/* if the jacobian is not rank deficient, the newton */ /* if the jacobian is not rank deficient, the newton */
/* step provides a lower bound, parl, for the zero of */ /* step provides a lower bound, parl, for the zero of */
/* the function. otherwise set this bound to zero. */ /* the function. otherwise set this bound to zero. */
parl = 0.; parl = 0.;
if (nsing < n-1) if (nsing >= n-1) {
goto L120; for (j = 0; j < n; ++j) {
for (j = 0; j < n; ++j) { l = ipvt[j];
l = ipvt[j]; wa1[j] = diag[l] * (wa2[l] / dxnorm);
wa1[j] = diag[l] * (wa2[l] / dxnorm); }
for (j = 0; j < n; ++j) {
sum = 0.;
for (i = 0; i < j; ++i)
sum += r(i,j) * wa1[i];
wa1[j] = (wa1[j] - sum) / r(j,j);
}
temp = wa1.blueNorm();
parl = fp / delta / temp / temp;
} }
for (j = 0; j < n; ++j) {
sum = 0.;
for (i = 0; i < j; ++i)
sum += r(i,j) * wa1[i];
wa1[j] = (wa1[j] - sum) / r(j,j);
}
temp = wa1.blueNorm();
parl = fp / delta / temp / temp;
L120:
/* calculate an upper bound, paru, for the zero of the function. */ /* calculate an upper bound, paru, for the zero of the function. */
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
sum = 0.; sum = 0.;
@ -97,89 +98,82 @@ L120:
} }
gnorm = wa1.stableNorm(); gnorm = wa1.stableNorm();
paru = gnorm / delta; paru = gnorm / delta;
if (paru == 0.) { if (paru == 0.)
paru = dwarf / std::min(delta,Scalar(0.1)); paru = dwarf / std::min(delta,Scalar(0.1));
}
/* if the input par lies outside of the interval (parl,paru), */ /* if the input par lies outside of the interval (parl,paru), */
/* set par to the closer endpoint. */ /* set par to the closer endpoint. */
par = std::max(par,parl); par = std::max(par,parl);
par = std::min(par,paru); par = std::min(par,paru);
if (par == 0.) if (par == 0.)
par = gnorm / dxnorm; par = gnorm / dxnorm;
/* beginning of an iteration. */ /* beginning of an iteration. */
L150: while (true) {
++iter; ++iter;
/* evaluate the function at the current value of par. */ /* evaluate the function at the current value of par. */
if (par == 0.) if (par == 0.)
par = std::max(dwarf,Scalar(.001) * paru); /* Computing MAX */ par = std::max(dwarf,Scalar(.001) * paru); /* Computing MAX */
temp = ei_sqrt(par); temp = ei_sqrt(par);
wa1 = temp * diag; wa1 = temp * diag;
ipvt.cwise()+=1; // qrsolv() expects the fortran convention (as qrfac provides) ipvt.cwise()+=1; // qrsolv() expects the fortran convention (as qrfac provides)
ei_qrsolv<Scalar>(n, r.data(), r.rows(), ipvt.data(), wa1.data(), qtb.data(), x.data(), sdiag.data(), wa2.data()); ei_qrsolv<Scalar>(n, r.data(), r.rows(), ipvt.data(), wa1.data(), qtb.data(), x.data(), sdiag.data(), wa2.data());
ipvt.cwise()-=1; ipvt.cwise()-=1;
wa2 = diag.cwise() * x; wa2 = diag.cwise() * x;
dxnorm = wa2.blueNorm(); dxnorm = wa2.blueNorm();
temp = fp; temp = fp;
fp = dxnorm - delta; fp = dxnorm - delta;
/* if the function is small enough, accept the current value */ /* if the function is small enough, accept the current value */
/* of par. also test for the exceptional cases where parl */ /* of par. also test for the exceptional cases where parl */
/* is zero or the number of iterations has reached 10. */ /* is zero or the number of iterations has reached 10. */
if (ei_abs(fp) <= Scalar(0.1) * delta || (parl == 0. && fp <= temp && temp < 0.) || iter == 10)
break;
/* compute the newton correction. */
for (j = 0; j < n; ++j) {
l = ipvt[j];
wa1[j] = diag[l] * (wa2[l] / dxnorm);
/* L180: */
}
for (j = 0; j < n; ++j) {
wa1[j] /= sdiag[j];
temp = wa1[j];
for (i = j+1; i < n; ++i)
wa1[i] -= r(i,j) * temp;
}
temp = wa1.blueNorm();
parc = fp / delta / temp / temp;
/* depending on the sign of the function, update parl or paru. */
if (fp > 0.)
parl = std::max(parl,par);
if (fp < 0.)
paru = std::min(paru,par);
/* compute an improved estimate for par. */
/* Computing MAX */
par = std::max(parl,par+parc);
/* end of an iteration. */
if (ei_abs(fp) <= Scalar(0.1) * delta || (parl == 0. && fp <= temp && temp < 0.) ||
iter == 10) {
goto L220;
} }
/* compute the newton correction. */ /* termination. */
for (j = 0; j < n; ++j) { if (iter == 0)
l = ipvt[j];
wa1[j] = diag[l] * (wa2[l] / dxnorm);
/* L180: */
}
for (j = 0; j < n; ++j) {
wa1[j] /= sdiag[j];
temp = wa1[j];
for (i = j+1; i < n; ++i)
wa1[i] -= r(i,j) * temp;
}
temp = wa1.blueNorm();
parc = fp / delta / temp / temp;
/* depending on the sign of the function, update parl or paru. */
if (fp > 0.) {
parl = std::max(parl,par);
}
if (fp < 0.) {
paru = std::min(paru,par);
}
/* compute an improved estimate for par. */
/* Computing MAX */
par = std::max(parl,par+parc);
/* end of an iteration. */
goto L150;
L220:
/* termination. */
if (iter == 0) {
par = 0.; par = 0.;
}
return; return;
} }