various cleaning and homogeneization

This commit is contained in:
Thomas Capricelli
2009-08-24 08:28:31 +02:00
parent 930651ff9a
commit 950eb4a254
7 changed files with 102 additions and 160 deletions

View File

@@ -192,19 +192,16 @@ L170:
/* beginning of the inner loop. */
L180:
/* if requested, call fcn to enable printing of iterates. */
/* if requested, call Functor::f to enable printing of iterates. */
if (nprint <= 0) {
goto L190;
}
iflag = 0;
if ((iter - 1) % nprint == 0) {
if ((iter - 1) % nprint == 0)
iflag = Functor::debug(x, fvec);
}
if (iflag < 0) {
if (iflag < 0)
goto L300;
}
L190:
/* determine the direction p. */
@@ -220,17 +217,15 @@ L190:
/* on the first iteration, adjust the initial step bound. */
if (iter == 1) {
if (iter == 1)
delta = std::min(delta,pnorm);
}
/* evaluate the function at x + p and calculate its norm. */
iflag = Functor::f(wa2, wa4);
++nfev;
if (iflag < 0) {
if (iflag < 0)
goto L300;
}
fnorm1 = wa4.stableNorm();
/* compute the scaled actual reduction. */
@@ -295,7 +290,7 @@ L240:
x = wa2;
wa2 = diag.cwise() * x;
fvec = wa4;
temp = wa2.stableNorm();
xnorm = wa2.stableNorm();
fnorm = fnorm1;
++iter;
L260:
@@ -376,12 +371,8 @@ L300:
if (iflag < 0) {
info = iflag;
}
if (nprint > 0) {
if (nprint > 0)
iflag = Functor::debug(x, fvec);
}
return info;
/* last card of subroutine hybrd. */
} /* hybrd_ */
}

View File

@@ -57,8 +57,7 @@ int ei_hybrj(
}
if (mode == 2)
for (j = 0; j < n; ++j)
if (diag[j] <= 0.)
goto L300;
if (diag[j] <= 0.) goto L300;
/* evaluate the function at the starting point */
/* and calculate its norm. */
@@ -68,7 +67,7 @@ int ei_hybrj(
if (iflag < 0) {
goto L300;
}
fnorm = fvec.stableNorm();;
fnorm = fvec.stableNorm();
/* initialize iteration counter and monitors. */
@@ -93,7 +92,7 @@ L30:
/* compute the qr factorization of the jacobian. */
ei_qrfac<Scalar>(n, n,fjac.data(), fjac.rows(), false, iwa, 1, wa1.data(), wa2.data(), wa3.data());
ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, 1, wa1.data(), wa2.data(), wa3.data());
/* on the first iteration and if mode is 1, scale according */
/* to the norms of the columns of the initial jacobian. */
@@ -117,7 +116,7 @@ L50:
/* and initialize the step bound delta. */
wa3 = diag.cwise() * x;
xnorm = wa3.stableNorm();;
xnorm = wa3.stableNorm();
delta = factor * xnorm;
if (delta == 0.) {
delta = factor;
@@ -175,14 +174,12 @@ L110:
goto L170;
}
/* Computing MAX */
for (j = 0; j < n; ++j)
diag[j] = std::max(diag[j], wa2[j]);
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) {
@@ -191,9 +188,8 @@ L180:
iflag = 0;
if ((iter - 1) % nprint == 0)
iflag = Functor::debug(x, fvec, fjac);
if (iflag < 0) {
if (iflag < 0)
goto L300;
}
L190:
/* determine the direction p. */
@@ -202,26 +198,22 @@ L190:
/* store the direction p and x + p. calculate the norm of p. */
for (j = 0; j < n; ++j) {
wa1[j] = -wa1[j];
wa2[j] = x[j] + wa1[j];
wa3[j] = diag[j] * wa1[j];
}
wa1 = -wa1;
wa2 = x + wa1;
wa3 = diag.cwise() * wa1;
pnorm = wa3.stableNorm();
/* on the first iteration, adjust the initial step bound. */
if (iter == 1) {
if (iter == 1)
delta = std::min(delta,pnorm);
}
/* evaluate the function at x + p and calculate its norm. */
iflag = Functor::f(wa2, wa4);
++nfev;
if (iflag < 0) {
if (iflag < 0)
goto L300;
}
fnorm1 = wa4.stableNorm();
/* compute the scaled actual reduction. */
@@ -283,7 +275,7 @@ L240:
/* successful iteration. update x, fvec, and their norms. */
x =wa2;
x = wa2;
wa2 = diag.cwise() * x;
fvec = wa4;
xnorm = wa2.stableNorm();
@@ -319,24 +311,19 @@ L260:
info = 2;
}
/* 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;
}
if (nslow2 == 5) {
if (nslow2 == 5)
info = 4;
}
if (nslow1 == 10) {
if (nslow1 == 10)
info = 5;
}
if (info != 0) {
if (info != 0)
goto L300;
}
/* criterion for recalculating jacobian. */
if (ncfail == 2) {
if (ncfail == 2)
goto L290;
}
/* calculate the rank one modification to the jacobian */
/* and update qtf if necessary. */
@@ -345,10 +332,8 @@ L260:
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)) {
if (ratio >= Scalar(1e-4))
qtf[j] = sum;
}
/* L280: */
}
/* compute the qr factorization of the updated jacobian. */
@@ -376,8 +361,5 @@ L300:
if (nprint > 0)
iflag = Functor::debug(x, fvec, fjac);
return info;
/* last card of subroutine hybrj. */
} /* hybrj_ */
}

View File

@@ -7,6 +7,7 @@ int ei_lmder(
int &njev,
Matrix< Scalar, Dynamic, Dynamic > &fjac,
VectorXi &ipvt,
Matrix< Scalar, Dynamic, 1 > &qtf,
Matrix< Scalar, Dynamic, 1 > &diag,
int mode=1,
Scalar factor = 100.,
@@ -18,11 +19,12 @@ int ei_lmder(
)
{
const int m = fvec.size(), n = x.size();
Matrix< Scalar, Dynamic, 1 > qtf(n), wa1(n), wa2(n), wa3(n), wa4(m);
Matrix< Scalar, Dynamic, 1 > wa1(n), wa2(n), wa3(n), wa4(m);
ipvt.resize(n);
fjac.resize(m, n);
diag.resize(n);
qtf.resize(n);
/* Local variables */
int i, j, l;
@@ -32,12 +34,11 @@ int ei_lmder(
int iflag;
Scalar delta;
Scalar ratio;
Scalar fnorm, gnorm, pnorm, xnorm, fnorm1, actred, dirder, prered;
Scalar fnorm, gnorm;
Scalar pnorm, xnorm, fnorm1, actred, dirder, prered;
int info;
/* Function Body */
info = 0;
iflag = 0;
nfev = 0;
@@ -51,8 +52,7 @@ int ei_lmder(
}
if (mode == 2)
for (j = 0; j < n; ++j)
if (diag[j] <= 0.)
goto L300;
if (diag[j] <= 0.) goto L300;
/* evaluate the function at the starting point */
/* and calculate its norm. */
@@ -188,8 +188,8 @@ L170:
if (mode == 2) {
goto L190;
}
for (j = 0; j < n; ++j)
diag[j] = std::max( diag[j], wa2[j]);
/* Computing MAX */
diag = diag.cwise().max(wa2);
L190:
/* beginning of the inner loop. */
@@ -235,7 +235,6 @@ L200:
wa3.fill(0.);
for (j = 0; j < n; ++j) {
wa3[j] = 0.;
l = ipvt[j];
temp = wa1[l];
for (i = 0; i <= j; ++i) {
@@ -245,7 +244,7 @@ L200:
/* L230: */
}
temp1 = ei_abs2(wa3.stableNorm() / fnorm);
temp2 = ei_abs2( ei_sqrt(par) * pnorm / fnorm);
temp2 = ei_abs2(ei_sqrt(par) * pnorm / fnorm);
/* Computing 2nd power */
prered = temp1 + temp2 / Scalar(.5);
dirder = -(temp1 + temp2);
@@ -269,11 +268,10 @@ L200:
if (actred < 0.) {
temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
}
if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) {
if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1))
temp = Scalar(.1);
}
/* Computing MIN */
delta = temp * std::min(delta, pnorm/Scalar(.1));
delta = temp * std::min(delta, pnorm / Scalar(.1));
par /= temp;
goto L260;
L240:
@@ -356,8 +354,5 @@ L300:
iflag = Functor::debug(x, fvec, fjac);
}
return info;
/* last card of subroutine lmder. */
} /* lmder_ */
}

View File

@@ -10,7 +10,7 @@ int ei_lmder1(
const int n = x.size(), m=fvec.size();
int info, nfev=0, njev=0;
Matrix< Scalar, Dynamic, Dynamic > fjac(m, n);
Matrix< Scalar, Dynamic, 1> diag;
Matrix< Scalar, Dynamic, 1> diag, qtf;
/* check the input parameters for errors. */
if (n <= 0 || m < n || tol < 0.) {
@@ -22,7 +22,7 @@ int ei_lmder1(
info = ei_lmder<Functor,Scalar>(
x, fvec,
nfev, njev,
fjac, ipvt, diag,
fjac, ipvt, qtf, diag,
1,
100.,
(n+1)*100,

View File

@@ -19,9 +19,7 @@ int ei_lmdif(
)
{
const int m = fvec.size(), n = x.size();
Matrix< Scalar, Dynamic, 1 >
wa1(n), wa2(n), wa3(n),
wa4(m);
Matrix< Scalar, Dynamic, 1 > wa1(n), wa2(n), wa3(n), wa4(m);
ipvt.resize(n);
fjac.resize(m, n);
@@ -53,8 +51,7 @@ int ei_lmdif(
}
if (mode == 2)
for (j = 0; j < n; ++j)
if (diag[j] <= 0.)
goto L300;
if (diag[j] <= 0.) goto L300;
/* evaluate the function at the starting point */
/* and calculate its norm. */
@@ -187,8 +184,8 @@ L170:
if (mode == 2) {
goto L190;
}
for (j = 0; j < n; ++j) /* Computing MAX */
diag[j] = std::max(diag[j], wa2[j]);
/* Computing MAX */
diag = diag.cwise().max(wa2);
L190:
/* beginning of the inner loop. */
@@ -232,8 +229,8 @@ L200:
/* compute the scaled predicted reduction and */
/* the scaled directional derivative. */
wa3.fill(0.);
for (j = 0; j < n; ++j) {
wa3[j] = 0.;
l = ipvt[j];
temp = wa1[l];
for (i = 0; i <= j; ++i) {
@@ -267,9 +264,8 @@ L200:
if (actred < 0.) {
temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
}
if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) {
if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1))
temp = Scalar(.1);
}
/* Computing MIN */
delta = temp * std::min(delta, pnorm / Scalar(.1));
par /= temp;
@@ -354,8 +350,5 @@ L300:
iflag = Functor::debug(x, fvec);
}
return info;
/* last card of subroutine lmdif. */
} /* lmdif_ */
}

View File

@@ -54,8 +54,7 @@ int ei_lmstr(
if (mode == 2)
for (j = 0; j < n; ++j)
if (diag[j] <= 0.)
goto L300;
if (diag[j] <= 0.) goto L300;
/* evaluate the function at the starting point */
/* and calculate its norm. */
@@ -95,23 +94,15 @@ L40:
/* forming (q transpose)*fvec and storing the first */
/* n components in qtf. */
for (j = 0; j < n; ++j) {
qtf[j] = 0.;
for (i = 0; i < n; ++i) {
fjac(i,j) = 0.;
/* L50: */
}
/* L60: */
}
qtf.fill(0.);
fjac.fill(0.);
iflag = 2;
for (i = 0; i < m; ++i) {
if (Functor::df(x, wa3, iflag) < 0) {
if (Functor::df(x, wa3, iflag) < 0)
goto L340;
}
temp = fvec[i];
ei_rwupdt<Scalar>(n, fjac.data(), fjac.rows(), wa3.data(), qtf.data(), &temp, wa1.data(), wa2.data());
++iflag;
/* L70: */
}
++njev;
@@ -126,25 +117,21 @@ L40:
ipvt[j] = j;
wa2[j] = fjac.col(j).start(j).stableNorm();
}
if (! sing) {
if (! sing)
goto L130;
}
ipvt.cwise()+=1;
ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), true, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data());
ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1)
for (j = 0; j < n; ++j) {
if (fjac(j,j) == 0.) {
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: */
}
L110:
fjac(j,j) = wa1[j];
@@ -173,10 +160,7 @@ L150:
/* on the first iteration, calculate the norm of the scaled x */
/* and initialize the step bound delta. */
for (j = 0; j < n; ++j) {
wa3[j] = diag[j] * x[j];
/* L160: */
}
wa3 = diag.cwise() * x;
xnorm = wa3.stableNorm();
delta = factor * xnorm;
if (delta == 0.) {
@@ -222,8 +206,8 @@ L210:
if (mode == 2) {
goto L230;
}
for (j = 0; j < n; ++j) /* Computing MAX */
diag[j] = std::max(diag[j], wa2[j]);
/* Computing MAX */
diag = diag.cwise().max(wa2);
L230:
/* beginning of the inner loop. */
@@ -390,8 +374,5 @@ L340:
iflag = Functor::debug(x, fvec, wa3);
}
return info;
/* last card of subroutine lmstr. */
} /* lmstr_ */
}