misc cleaning

This commit is contained in:
Thomas Capricelli 2009-08-23 06:16:05 +02:00
parent 264e61932c
commit 9a8c5cbd2c
6 changed files with 41 additions and 63 deletions

View File

@ -20,10 +20,9 @@ void ei_chkder(
int i,j; int i,j;
const int m = fvec.size(), n = x.size(); const int m = fvec.size(), n = x.size();
int ldfjac = m;
if (mode != 2) { if (mode != 2) {
xp.resize(ldfjac); xp.resize(m);
/* mode = 1. */ /* mode = 1. */
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
temp = eps * ei_abs(x[j]); temp = eps * ei_abs(x[j]);
@ -34,7 +33,7 @@ void ei_chkder(
} }
else { else {
/* mode = 2. */ /* mode = 2. */
err.setZero(ldfjac); err.setZero(m);
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
temp = ei_abs(x[j]); temp = ei_abs(x[j]);
if (temp == 0.) if (temp == 0.)
@ -43,21 +42,14 @@ void ei_chkder(
} }
for (i = 0; i < m; ++i) { for (i = 0; i < m; ++i) {
temp = 1.; temp = 1.;
if (fvec[i] != 0. && fvecp[i] != 0. && ei_abs(fvecp[i] - if (fvec[i] != 0. && fvecp[i] != 0. && ei_abs(fvecp[i] - fvec[i]) >= epsf * ei_abs(fvec[i]))
fvec[i]) >= epsf * ei_abs(fvec[i])) temp = eps * ei_abs((fvecp[i] - fvec[i]) / eps - err[i]) / (ei_abs(fvec[i]) + ei_abs(fvecp[i]));
{
temp = eps * ei_abs((fvecp[i] - fvec[i]) / eps - err[i])
/ (ei_abs(fvec[i]) +
ei_abs(fvecp[i]));
}
err[i] = 1.; err[i] = 1.;
if (temp > epsilon<Scalar>() && temp < eps) { if (temp > epsilon<Scalar>() && temp < eps)
err[i] = (chkder_log10e * ei_log(temp) - epslog) / epslog; err[i] = (chkder_log10e * ei_log(temp) - epslog) / epslog;
} if (temp >= eps)
if (temp >= eps) {
err[i] = 0.; err[i] = 0.;
} }
} }
}
} /* chkder_ */ } /* chkder_ */

View File

@ -28,8 +28,7 @@ int ei_hybrd(
fvec.resize(n); fvec.resize(n);
qtf.resize(n); qtf.resize(n);
R.resize(lr); R.resize(lr);
int ldfjac = n; fjac.resize(n, n);
fjac.resize(ldfjac, n);
/* Local variables */ /* Local variables */
int i, j, l, iwa[1]; int i, j, l, iwa[1];
@ -58,7 +57,7 @@ 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. || ldfjac < n || lr < n * (n + 1) / 2) { factor <= 0. || lr < n * (n + 1) / 2) {
goto L300; goto L300;
} }
if (mode == 2) if (mode == 2)
@ -105,7 +104,7 @@ L30:
/* compute the qr factorization of the jacobian. */ /* compute the qr factorization of the jacobian. */
ei_qrfac<Scalar>(n, n, fjac.data(), ldfjac, 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 */ /* 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. */
@ -179,7 +178,7 @@ L110:
/* accumulate the orthogonal factor in fjac. */ /* accumulate the orthogonal factor in fjac. */
ei_qform<Scalar>(n, n, fjac.data(), ldfjac, wa1.data()); ei_qform<Scalar>(n, n, fjac.data(), fjac.rows(), wa1.data());
/* rescale if necessary. */ /* rescale if necessary. */
@ -358,7 +357,7 @@ 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(), lr, wa1.data(), wa2.data(), wa3.data(), &sing);
ei_r1mpyq<Scalar>(n, n, fjac.data(), ldfjac, 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. */

View File

@ -23,8 +23,7 @@ int ei_hybrj(
fvec.resize(n); fvec.resize(n);
qtf.resize(n); qtf.resize(n);
R.resize(lr); R.resize(lr);
int ldfjac = n; fjac.resize(n, n);
fjac.resize(ldfjac, n);
/* Local variables */ /* Local variables */
int i, j, l, iwa[1]; int i, j, l, iwa[1];
@ -52,7 +51,7 @@ int ei_hybrj(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || ldfjac < n || xtol < 0. || maxfev <= 0 || factor <= if (n <= 0 || xtol < 0. || maxfev <= 0 || factor <=
0. || lr < n * (n + 1) / 2) { 0. || lr < n * (n + 1) / 2) {
goto L300; goto L300;
} }
@ -94,7 +93,7 @@ L30:
/* compute the qr factorization of the jacobian. */ /* compute the qr factorization of the jacobian. */
ei_qrfac<Scalar>(n, n,fjac.data(), ldfjac, 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 */ /* 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. */
@ -168,7 +167,7 @@ L110:
/* accumulate the orthogonal factor in fjac. */ /* accumulate the orthogonal factor in fjac. */
ei_qform<Scalar>(n, n, fjac.data(), ldfjac, wa1.data()); ei_qform<Scalar>(n, n, fjac.data(), fjac.rows(), wa1.data());
/* rescale if necessary. */ /* rescale if necessary. */
@ -355,7 +354,7 @@ 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(), lr, wa1.data(), wa2.data(), wa3.data(), &sing);
ei_r1mpyq<Scalar>(n, n, fjac.data(), ldfjac, 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. */

View File

@ -19,10 +19,9 @@ int ei_lmder(
{ {
const int m = fvec.size(), n = x.size(); 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 > qtf(n), wa1(n), wa2(n), wa3(n), wa4(m);
int ldfjac = m;
ipvt.resize(n); ipvt.resize(n);
fjac.resize(ldfjac, n); fjac.resize(m, n);
diag.resize(n); diag.resize(n);
/* Local variables */ /* Local variables */
@ -46,7 +45,7 @@ int ei_lmder(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || m < n || ldfjac < m || ftol < 0. || xtol < 0. || if (n <= 0 || m < n || ftol < 0. || xtol < 0. ||
gtol < 0. || maxfev <= 0 || factor <= 0.) { gtol < 0. || maxfev <= 0 || factor <= 0.) {
goto L300; goto L300;
} }
@ -98,7 +97,7 @@ L40:
/* compute the qr factorization of the jacobian. */ /* compute the qr factorization of the jacobian. */
ei_qrfac<Scalar>(m, n, fjac.data(), ldfjac, true, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data()); ei_qrfac<Scalar>(m, 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) ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1)
/* on the first iteration and if mode is 1, scale according */ /* on the first iteration and if mode is 1, scale according */
@ -200,7 +199,7 @@ L200:
/* determine the levenberg-marquardt parameter. */ /* determine the levenberg-marquardt parameter. */
ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides) ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
ei_lmpar<Scalar>(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta, ei_lmpar<Scalar>(n, fjac.data(), fjac.rows(), ipvt.data(), diag.data(), qtf.data(), delta,
par, wa1.data(), wa2.data(), wa3.data(), wa4.data()); par, wa1.data(), wa2.data(), wa3.data(), wa4.data());
ipvt.cwise()-=1; ipvt.cwise()-=1;

View File

@ -22,10 +22,9 @@ int ei_lmdif(
Matrix< Scalar, Dynamic, 1 > Matrix< Scalar, Dynamic, 1 >
wa1(n), wa2(n), wa3(n), wa1(n), wa2(n), wa3(n),
wa4(m); wa4(m);
int ldfjac = m;
ipvt.resize(n); ipvt.resize(n);
fjac.resize(ldfjac, n); fjac.resize(m, n);
diag.resize(n); diag.resize(n);
qtf.resize(n); qtf.resize(n);
@ -48,7 +47,7 @@ int ei_lmdif(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || m < n || ldfjac < m || ftol < 0. || xtol < 0. || if (n <= 0 || m < n || ftol < 0. || xtol < 0. ||
gtol < 0. || maxfev <= 0 || factor <= 0.) { gtol < 0. || maxfev <= 0 || factor <= 0.) {
goto L300; goto L300;
} }
@ -100,7 +99,7 @@ L40:
/* compute the qr factorization of the jacobian. */ /* compute the qr factorization of the jacobian. */
ei_qrfac<Scalar>(m, n, fjac.data(), ldfjac, true, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data()); ei_qrfac<Scalar>(m, 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) ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1)
/* on the first iteration and if mode is 1, scale according */ /* on the first iteration and if mode is 1, scale according */
@ -137,21 +136,21 @@ L80:
wa4 = fvec; wa4 = fvec;
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
if (fjac[j + j * ldfjac] == 0.) { if (fjac(j,j) == 0.) {
goto L120; goto L120;
} }
sum = 0.; sum = 0.;
for (i = j; i < m; ++i) { for (i = j; i < m; ++i) {
sum += fjac[i + j * ldfjac] * wa4[i]; sum += fjac(i,j) * wa4[i];
/* L100: */ /* L100: */
} }
temp = -sum / fjac[j + j * ldfjac]; temp = -sum / fjac(j,j);
for (i = j; i < m; ++i) { for (i = j; i < m; ++i) {
wa4[i] += fjac[i + j * ldfjac] * temp; wa4[i] += fjac(i,j) * temp;
/* L110: */ /* L110: */
} }
L120: L120:
fjac[j + j * ldfjac] = wa1[j]; fjac(j,j) = wa1[j];
qtf[j] = wa4[j]; qtf[j] = wa4[j];
/* L130: */ /* L130: */
} }
@ -164,19 +163,13 @@ L120:
} }
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
l = ipvt[j]; l = ipvt[j];
if (wa2[l] == 0.) { if (wa2[l] != 0.) {
goto L150;
}
sum = 0.; sum = 0.;
for (i = 0; i <= j; ++i) { for (i = 0; i <= j; ++i)
sum += fjac[i + j * ldfjac] * (qtf[i] / fnorm); sum += fjac(i,j) * (qtf[i] / fnorm);
/* L140: */
}
/* Computing MAX */ /* Computing MAX */
gnorm = std::max(gnorm, ei_abs(sum / wa2[l])); gnorm = std::max(gnorm, ei_abs(sum / wa2[l]));
L150: }
/* L160: */
;
} }
L170: L170:
@ -205,7 +198,7 @@ L200:
/* determine the levenberg-marquardt parameter. */ /* determine the levenberg-marquardt parameter. */
ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides) ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
ei_lmpar<Scalar>(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta, ei_lmpar<Scalar>(n, fjac.data(), fjac.rows(), ipvt.data(), diag.data(), qtf.data(), delta,
par, wa1.data(), wa2.data(), wa3.data(), wa4.data()); par, wa1.data(), wa2.data(), wa3.data(), wa4.data());
ipvt.cwise()-=1; ipvt.cwise()-=1;
@ -245,7 +238,7 @@ L200:
l = ipvt[j]; l = ipvt[j];
temp = wa1[l]; temp = wa1[l];
for (i = 0; i <= j; ++i) { for (i = 0; i <= j; ++i) {
wa3[i] += fjac[i + j * ldfjac] * temp; wa3[i] += fjac(i,j) * temp;
/* L220: */ /* L220: */
} }
/* L230: */ /* L230: */

View File

@ -22,10 +22,9 @@ int ei_lmstr(
qtf(n), qtf(n),
wa1(n), wa2(n), wa3(n), wa1(n), wa2(n), wa3(n),
wa4(m); wa4(m);
int ldfjac = m;
ipvt.resize(n); ipvt.resize(n);
fjac.resize(ldfjac, n); fjac.resize(m, n);
diag.resize(n); diag.resize(n);
/* Local variables */ /* Local variables */
@ -48,7 +47,7 @@ int ei_lmstr(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || m < n || ldfjac < n || ftol < 0. || xtol < 0. || if (n <= 0 || m < n || ftol < 0. || xtol < 0. ||
gtol < 0. || maxfev <= 0 || factor <= 0.) { gtol < 0. || maxfev <= 0 || factor <= 0.) {
goto L340; goto L340;
} }
@ -110,7 +109,7 @@ L40:
goto L340; goto L340;
} }
temp = fvec[i]; temp = fvec[i];
ei_rwupdt<Scalar>(n, fjac.data(), ldfjac, wa3.data(), qtf.data(), &temp, wa1.data(), wa2.data()); ei_rwupdt<Scalar>(n, fjac.data(), fjac.rows(), wa3.data(), qtf.data(), &temp, wa1.data(), wa2.data());
++iflag; ++iflag;
/* L70: */ /* L70: */
} }
@ -126,15 +125,12 @@ L40:
} }
ipvt[j] = j; ipvt[j] = j;
wa2[j] = fjac.col(j).start(j).stableNorm(); wa2[j] = fjac.col(j).start(j).stableNorm();
// wa2[j] = ei_enorm<Scalar>(j, &fjac[j * ldfjac + 1]);
// sum += fjac[i + j * ldfjac] * (qtf[i] / fnorm);
/* L80: */
} }
if (! sing) { if (! sing) {
goto L130; goto L130;
} }
ipvt.cwise()+=1; ipvt.cwise()+=1;
ei_qrfac<Scalar>(n, n, fjac.data(), ldfjac, true, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data()); 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) ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1)
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
if (fjac(j,j) == 0.) { if (fjac(j,j) == 0.) {
@ -237,7 +233,7 @@ L240:
/* determine the levenberg-marquardt parameter. */ /* determine the levenberg-marquardt parameter. */
ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides) ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
ei_lmpar<Scalar>(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta, par, ei_lmpar<Scalar>(n, fjac.data(), fjac.rows(), ipvt.data(), diag.data(), qtf.data(), delta, par,
wa1.data(), wa2.data(), wa3.data(), wa4.data()); wa1.data(), wa2.data(), wa3.data(), wa4.data());
ipvt.cwise()-=1; ipvt.cwise()-=1;