GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 19:59:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

clean, remove goto's

Browse Source
This commit is contained in:
Thomas Capricelli 2009-08-24 15:32:06 +02:00
parent d4968cd059
commit 91a2145cb3
2 changed files with 436 additions and 542 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

133
unsupported/Eigen/src/NonLinear/hybrd.h
View File

@ -19,7 +19,6 @@ int ei_hybrd(
) )
{ {
const int n = x.size(); const int n = x.size();
int lr = (n*(n+1))/2;
Matrix< Scalar, Dynamic, 1 > wa1(n), wa2(n), wa3(n), wa4(n); 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; if (nb_of_superdiagonals<0) nb_of_superdiagonals = n-1;
fvec.resize(n); fvec.resize(n);
qtf.resize(n); qtf.resize(n);
R.resize(lr); R.resize( (n*(n+1))/2);
fjac.resize(n, n); fjac.resize(n, n);
/* Local variables */ /* Local variables */
@ -56,10 +55,8 @@ int ei_hybrd(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || xtol < 0. || maxfev <= 0 || nb_of_subdiagonals < 0 || nb_of_superdiagonals < 0 || if (n <= 0 || xtol < 0. || maxfev <= 0 || nb_of_subdiagonals < 0 || nb_of_superdiagonals < 0 || factor <= 0. )
factor <= 0. || lr < n * (n + 1) / 2) {
goto L300; goto L300;
}
if (mode == 2) if (mode == 2)
for (j = 0; j < n; ++j) for (j = 0; j < n; ++j)
if (diag[j] <= 0.) goto L300; if (diag[j] <= 0.) goto L300;
@ -69,9 +66,8 @@ int ei_hybrd(
iflag = Functor::f(x, fvec); iflag = Functor::f(x, fvec);
nfev = 1; nfev = 1;
if (iflag < 0) { if (iflag < 0)
goto L300; goto L300;
}
fnorm = fvec.stableNorm(); fnorm = fvec.stableNorm();
/* determine the number of calls to fcn needed to compute */ /* determine the number of calls to fcn needed to compute */
@ -90,7 +86,7 @@ int ei_hybrd(
/* beginning of the outer loop. */ /* beginning of the outer loop. */
L30: while (true) {
jeval = true; jeval = true;
/* calculate the jacobian matrix. */ /* calculate the jacobian matrix. */
@ -98,9 +94,8 @@ L30:
iflag = ei_fdjac1<Functor,Scalar>(x, fvec, fjac, iflag = ei_fdjac1<Functor,Scalar>(x, fvec, fjac,
nb_of_subdiagonals, nb_of_superdiagonals, epsfcn, wa1, wa2); nb_of_subdiagonals, nb_of_superdiagonals, epsfcn, wa1, wa2);
nfev += msum; nfev += msum;
if (iflag < 0) { if (iflag < 0)
goto L300; break;
}
/* compute the qr factorization of the jacobian. */ /* compute the qr factorization of the jacobian. */
@ -109,20 +104,13 @@ L30:
/* on the first iteration and if mode is 1, scale according */ /* on the first iteration and if mode is 1, scale according */
/* to the norms of the columns of the initial jacobian. */ /* to the norms of the columns of the initial jacobian. */
if (iter != 1) { if (iter == 1) {
goto L70; if (mode != 2)
}
if (mode == 2) {
goto L50;
}
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
diag[j] = wa2[j]; diag[j] = wa2[j];
if (wa2[j] == 0.) { if (wa2[j] == 0.)
diag[j] = 1.; diag[j] = 1.;
} }
/* L40: */
}
L50:
/* on the first iteration, calculate the norm of the scaled x */ /* on the first iteration, calculate the norm of the scaled x */
/* and initialize the step bound delta. */ /* and initialize the step bound delta. */
@ -130,31 +118,21 @@ L50:
wa3 = diag.cwise() * x; wa3 = diag.cwise() * x;
xnorm = wa3.stableNorm(); xnorm = wa3.stableNorm();
delta = factor * xnorm; delta = factor * xnorm;
if (delta == 0.) { if (delta == 0.)
delta = factor; delta = factor;
} }
L70:
/* form (q transpose)*fvec and store in qtf. */ /* form (q transpose)*fvec and store in qtf. */
qtf = fvec; qtf = fvec;
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j)
if (fjac(j,j) == 0.) { if (fjac(j,j) != 0.) {
goto L110;
}
sum = 0.; sum = 0.;
for (i = j; i < n; ++i) { for (i = j; i < n; ++i)
sum += fjac(i,j) * qtf[i]; sum += fjac(i,j) * qtf[i];
/* L90: */
}
temp = -sum / fjac(j,j); temp = -sum / fjac(j,j);
for (i = j; i < n; ++i) { for (i = j; i < n; ++i)
qtf[i] += fjac(i,j) * temp; qtf[i] += fjac(i,j) * temp;
/* L100: */
}
L110:
/* L120: */
;
} }
/* copy the triangular factor of the qr factorization into r. */ /* copy the triangular factor of the qr factorization into r. */
@ -162,19 +140,15 @@ L110:
sing = false; sing = false;
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
l = j; l = j;
if (j) { if (j)
for (i = 0; i < j; ++i) { for (i = 0; i < j; ++i) {
R[l] = fjac(i,j); R[l] = fjac(i,j);
l = l + n - i -1; l = l + n - i -1;
/* L130: */
}
} }
R[l] = wa1[j]; R[l] = wa1[j];
if (wa1[j] == 0.) { if (wa1[j] == 0.)
sing = true; sing = true;
} }
/* L150: */
}
/* accumulate the orthogonal factor in fjac. */ /* accumulate the orthogonal factor in fjac. */
@ -182,27 +156,22 @@ L110:
/* rescale if necessary. */ /* rescale if necessary. */
if (mode == 2) {
goto L170;
}
/* Computing MAX */ /* Computing MAX */
if (mode != 2)
diag = diag.cwise().max(wa2); diag = diag.cwise().max(wa2);
L170:
/* beginning of the inner loop. */ /* beginning of the inner loop. */
L180: while (true) {
/* if requested, call Functor::f to enable printing of iterates. */ /* if requested, call Functor::f to enable printing of iterates. */
if (nprint <= 0) { if (nprint > 0) {
goto L190;
}
iflag = 0; iflag = 0;
if ((iter - 1) % nprint == 0) if ((iter - 1) % nprint == 0)
iflag = Functor::debug(x, fvec); iflag = Functor::debug(x, fvec);
if (iflag < 0) if (iflag < 0)
goto L300; goto L300;
L190: }
/* determine the direction p. */ /* determine the direction p. */
@ -242,10 +211,8 @@ L190:
for (j = i; j < n; ++j) { for (j = i; j < n; ++j) {
sum += R[l] * wa1[j]; sum += R[l] * wa1[j];
++l; ++l;
/* L210: */
} }
wa3[i] = qtf[i] + sum; wa3[i] = qtf[i] + sum;
/* L220: */
} }
temp = wa3.stableNorm(); temp = wa3.stableNorm();
prered = 0.; prered = 0.;
@ -256,20 +223,16 @@ L190:
/* reduction. */ /* reduction. */
ratio = 0.; ratio = 0.;
if (prered > 0.) { if (prered > 0.)
ratio = actred / prered; ratio = actred / prered;
}
/* update the step bound. */ /* update the step bound. */
if (ratio >= Scalar(.1)) { if (ratio < Scalar(.1)) {
goto L230;
}
ncsuc = 0; ncsuc = 0;
++ncfail; ++ncfail;
delta = Scalar(.5) * delta; delta = Scalar(.5) * delta;
goto L240; } else {
L230:
ncfail = 0; ncfail = 0;
++ncsuc; ++ncsuc;
if (ratio >= Scalar(.5) || ncsuc > 1) /* Computing MAX */ if (ratio >= Scalar(.5) || ncsuc > 1) /* Computing MAX */
@ -277,51 +240,41 @@ L230:
if (ei_abs(ratio - 1.) <= Scalar(.1)) { if (ei_abs(ratio - 1.) <= Scalar(.1)) {
delta = pnorm / Scalar(.5); delta = pnorm / Scalar(.5);
} }
L240: }
/* test for successful iteration. */ /* test for successful iteration. */
if (ratio < Scalar(1e-4)) { if (ratio >= Scalar(1e-4)) {
goto L260;
}
/* successful iteration. update x, fvec, and their norms. */ /* successful iteration. update x, fvec, and their norms. */
x = wa2; x = wa2;
wa2 = diag.cwise() * x; wa2 = diag.cwise() * x;
fvec = wa4; fvec = wa4;
xnorm = wa2.stableNorm(); xnorm = wa2.stableNorm();
fnorm = fnorm1; fnorm = fnorm1;
++iter; ++iter;
L260: }
/* determine the progress of the iteration. */ /* determine the progress of the iteration. */
++nslow1; ++nslow1;
if (actred >= Scalar(.001)) { if (actred >= Scalar(.001))
nslow1 = 0; nslow1 = 0;
} if (jeval)
if (jeval) {
++nslow2; ++nslow2;
} if (actred >= Scalar(.1))
if (actred >= Scalar(.1)) {
nslow2 = 0; nslow2 = 0;
}
/* test for convergence. */ /* test for convergence. */
if (delta <= xtol * xnorm || fnorm == 0.) { if (delta <= xtol * xnorm || fnorm == 0.)
info = 1; info = 1;
} if (info != 0)
if (info != 0) {
goto L300; goto L300;
}
/* tests for termination and stringent tolerances. */ /* tests for termination and stringent tolerances. */
if (nfev >= maxfev) { if (nfev >= maxfev)
info = 2; info = 2;
}
/* Computing MAX */ /* Computing MAX */
if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm) if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm)
info = 3; info = 3;
@ -336,7 +289,7 @@ L260:
/* by forward differences. */ /* by forward differences. */
if (ncfail == 2) if (ncfail == 2)
goto L290; break;
/* calculate the rank one modification to the jacobian */ /* calculate the rank one modification to the jacobian */
/* and update qtf if necessary. */ /* and update qtf if necessary. */
@ -351,26 +304,20 @@ L260:
/* compute the qr factorization of the updated jacobian. */ /* compute the qr factorization of the updated jacobian. */
ei_r1updt<Scalar>(n, n, R.data(), lr, wa1.data(), wa2.data(), wa3.data(), &sing); 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>(n, n, fjac.data(), fjac.rows(), wa2.data(), wa3.data());
ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data()); ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data());
/* end of the inner loop. */ /* end of the inner loop. */
jeval = false; jeval = false;
goto L180;
L290:
/* end of the outer loop. */
goto L30;
L300:
/* termination, either normal or user imposed. */
if (iflag < 0) {
info = iflag;
} }
/* end of the outer loop. */
}
L300:
/* termination, either normal or user imposed. */
if (iflag < 0)
info = iflag;
if (nprint > 0) if (nprint > 0)
iflag = Functor::debug(x, fvec); iflag = Functor::debug(x, fvec);
return info; return info;

133
unsupported/Eigen/src/NonLinear/hybrj.h
View File

@ -17,12 +17,11 @@ int ei_hybrj(
) )
{ {
const int n = x.size(); const int n = x.size();
const int lr = (n*(n+1))/2;
Matrix< Scalar, Dynamic, 1 > wa1(n), wa2(n), wa3(n), wa4(n); Matrix< Scalar, Dynamic, 1 > wa1(n), wa2(n), wa3(n), wa4(n);
fvec.resize(n); fvec.resize(n);
qtf.resize(n); qtf.resize(n);
R.resize(lr); R.resize( (n*(n+1))/2);
fjac.resize(n, n); fjac.resize(n, n);
/* Local variables */ /* Local variables */
@ -51,10 +50,8 @@ int ei_hybrj(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || xtol < 0. || maxfev <= 0 || factor <= if (n <= 0 || xtol < 0. || maxfev <= 0 || factor <= 0. )
0. || lr < n * (n + 1) / 2) {
goto L300; goto L300;
}
if (mode == 2) if (mode == 2)
for (j = 0; j < n; ++j) for (j = 0; j < n; ++j)
if (diag[j] <= 0.) goto L300; if (diag[j] <= 0.) goto L300;
@ -64,9 +61,8 @@ int ei_hybrj(
iflag = Functor::f(x, fvec); iflag = Functor::f(x, fvec);
nfev = 1; nfev = 1;
if (iflag < 0) { if (iflag < 0)
goto L300; goto L300;
}
fnorm = fvec.stableNorm(); fnorm = fvec.stableNorm();
/* initialize iteration counter and monitors. */ /* initialize iteration counter and monitors. */
@ -79,16 +75,15 @@ int ei_hybrj(
/* beginning of the outer loop. */ /* beginning of the outer loop. */
L30: while (true) {
jeval = true; jeval = true;
/* calculate the jacobian matrix. */ /* calculate the jacobian matrix. */
iflag = Functor::df(x, fjac); iflag = Functor::df(x, fjac);
++njev; ++njev;
if (iflag < 0) { if (iflag < 0)
goto L300; break;
}
/* compute the qr factorization of the jacobian. */ /* compute the qr factorization of the jacobian. */
@ -97,20 +92,13 @@ L30:
/* on the first iteration and if mode is 1, scale according */ /* on the first iteration and if mode is 1, scale according */
/* to the norms of the columns of the initial jacobian. */ /* to the norms of the columns of the initial jacobian. */
if (iter != 1) { if (iter == 1) {
goto L70; if (mode != 2)
}
if (mode == 2) {
goto L50;
}
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
diag[j] = wa2[j]; diag[j] = wa2[j];
if (wa2[j] == 0.) { if (wa2[j] == 0.)
diag[j] = 1.; diag[j] = 1.;
} }
/* L40: */
}
L50:
/* on the first iteration, calculate the norm of the scaled x */ /* on the first iteration, calculate the norm of the scaled x */
/* and initialize the step bound delta. */ /* and initialize the step bound delta. */
@ -118,31 +106,21 @@ L50:
wa3 = diag.cwise() * x; wa3 = diag.cwise() * x;
xnorm = wa3.stableNorm(); xnorm = wa3.stableNorm();
delta = factor * xnorm; delta = factor * xnorm;
if (delta == 0.) { if (delta == 0.)
delta = factor; delta = factor;
} }
L70:
/* form (q transpose)*fvec and store in qtf. */ /* form (q transpose)*fvec and store in qtf. */
qtf = fvec; qtf = fvec;
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j)
if (fjac(j,j) == 0.) { if (fjac(j,j) != 0.) {
goto L110;
}
sum = 0.; sum = 0.;
for (i = j; i < n; ++i) { for (i = j; i < n; ++i)
sum += fjac(i,j) * qtf[i]; sum += fjac(i,j) * qtf[i];
/* L90: */
}
temp = -sum / fjac(j,j); temp = -sum / fjac(j,j);
for (i = j; i < n; ++i) { for (i = j; i < n; ++i)
qtf[i] += fjac(i,j) * temp; qtf[i] += fjac(i,j) * temp;
/* L100: */
}
L110:
/* L120: */
;
} }
/* copy the triangular factor of the qr factorization into r. */ /* copy the triangular factor of the qr factorization into r. */
@ -150,19 +128,15 @@ L110:
sing = false; sing = false;
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
l = j; l = j;
if (j) { if (j)
for (i = 0; i < j; ++i) { for (i = 0; i < j; ++i) {
R[l] = fjac(i,j); R[l] = fjac(i,j);
l = l + n - i -1; l = l + n - i -1;
/* L130: */
}
} }
R[l] = wa1[j]; R[l] = wa1[j];
if (wa1[j] == 0.) { if (wa1[j] == 0.)
sing = true; sing = true;
} }
/* L150: */
}
/* accumulate the orthogonal factor in fjac. */ /* accumulate the orthogonal factor in fjac. */
@ -170,27 +144,22 @@ L110:
/* rescale if necessary. */ /* rescale if necessary. */
if (mode == 2) {
goto L170;
}
/* Computing MAX */ /* Computing MAX */
if (mode != 2)
diag = diag.cwise().max(wa2); diag = diag.cwise().max(wa2);
L170:
/* beginning of the inner loop. */ /* beginning of the inner loop. */
L180: while (true) {
/* if requested, call Functor::f to enable printing of iterates. */ /* if requested, call Functor::f to enable printing of iterates. */
if (nprint <= 0) { if (nprint > 0) {
goto L190;
}
iflag = 0; iflag = 0;
if ((iter - 1) % nprint == 0) if ((iter - 1) % nprint == 0)
iflag = Functor::debug(x, fvec, fjac); iflag = Functor::debug(x, fvec, fjac);
if (iflag < 0) if (iflag < 0)
goto L300; goto L300;
L190: }
/* determine the direction p. */ /* determine the direction p. */
@ -230,10 +199,8 @@ L190:
for (j = i; j < n; ++j) { for (j = i; j < n; ++j) {
sum += R[l] * wa1[j]; sum += R[l] * wa1[j];
++l; ++l;
/* L210: */
} }
wa3[i] = qtf[i] + sum; wa3[i] = qtf[i] + sum;
/* L220: */
} }
temp = wa3.stableNorm(); temp = wa3.stableNorm();
prered = 0.; prered = 0.;
@ -244,20 +211,16 @@ L190:
/* reduction. */ /* reduction. */
ratio = 0.; ratio = 0.;
if (prered > 0.) { if (prered > 0.)
ratio = actred / prered; ratio = actred / prered;
}
/* update the step bound. */ /* update the step bound. */
if (ratio >= Scalar(.1)) { if (ratio < Scalar(.1)) {
goto L230;
}
ncsuc = 0; ncsuc = 0;
++ncfail; ++ncfail;
delta = Scalar(.5) * delta; delta = Scalar(.5) * delta;
goto L240; } else {
L230:
ncfail = 0; ncfail = 0;
++ncsuc; ++ncsuc;
if (ratio >= Scalar(.5) || ncsuc > 1) /* Computing MAX */ if (ratio >= Scalar(.5) || ncsuc > 1) /* Computing MAX */
@ -265,51 +228,41 @@ L230:
if (ei_abs(ratio - 1.) <= Scalar(.1)) { if (ei_abs(ratio - 1.) <= Scalar(.1)) {
delta = pnorm / Scalar(.5); delta = pnorm / Scalar(.5);
} }
L240: }
/* test for successful iteration. */ /* test for successful iteration. */
if (ratio < Scalar(1e-4)) { if (ratio >= Scalar(1e-4)) {
goto L260;
}
/* successful iteration. update x, fvec, and their norms. */ /* successful iteration. update x, fvec, and their norms. */
x = wa2; x = wa2;
wa2 = diag.cwise() * x; wa2 = diag.cwise() * x;
fvec = wa4; fvec = wa4;
xnorm = wa2.stableNorm(); xnorm = wa2.stableNorm();
fnorm = fnorm1; fnorm = fnorm1;
++iter; ++iter;
L260: }
/* determine the progress of the iteration. */ /* determine the progress of the iteration. */
++nslow1; ++nslow1;
if (actred >= Scalar(.001)) { if (actred >= Scalar(.001))
nslow1 = 0; nslow1 = 0;
} if (jeval)
if (jeval) {
++nslow2; ++nslow2;
} if (actred >= Scalar(.1))
if (actred >= Scalar(.1)) {
nslow2 = 0; nslow2 = 0;
}
/* test for convergence. */ /* test for convergence. */
if (delta <= xtol * xnorm || fnorm == 0.) { if (delta <= xtol * xnorm || fnorm == 0.)
info = 1; info = 1;
} if (info != 0)
if (info != 0) {
goto L300; goto L300;
}
/* tests for termination and stringent tolerances. */ /* tests for termination and stringent tolerances. */
if (nfev >= maxfev) { if (nfev >= maxfev)
info = 2; info = 2;
}
/* Computing MAX */ /* Computing MAX */
if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm) if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm)
info = 3; info = 3;
@ -323,7 +276,7 @@ L260:
/* criterion for recalculating jacobian. */ /* criterion for recalculating jacobian. */
if (ncfail == 2) if (ncfail == 2)
goto L290; break;
/* calculate the rank one modification to the jacobian */ /* calculate the rank one modification to the jacobian */
/* and update qtf if necessary. */ /* and update qtf if necessary. */
@ -338,26 +291,20 @@ L260:
/* compute the qr factorization of the updated jacobian. */ /* compute the qr factorization of the updated jacobian. */
ei_r1updt<Scalar>(n, n, R.data(), lr, wa1.data(), wa2.data(), wa3.data(), &sing); 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>(n, n, fjac.data(), fjac.rows(), wa2.data(), wa3.data());
ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data()); ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data());
/* end of the inner loop. */ /* end of the inner loop. */
jeval = false; jeval = false;
goto L180;
L290:
/* end of the outer loop. */
goto L30;
L300:
/* termination, either normal or user imposed. */
if (iflag < 0) {
info = iflag;
} }
/* end of the outer loop. */
}
L300:
/* termination, either normal or user imposed. */
if (iflag < 0)
info = iflag;
if (nprint > 0) if (nprint > 0)
iflag = Functor::debug(x, fvec, fjac); iflag = Functor::debug(x, fvec, fjac);
return info; return info;