port ei_rwupdt to c++, and misc cleaning

This commit is contained in:
Thomas Capricelli 2010-01-27 09:01:13 +01:00
parent e97529c2e3
commit 7ba9dc07ed
2 changed files with 30 additions and 36 deletions

View File

@ -289,7 +289,6 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
if (mode != 2) if (mode != 2)
diag = diag.cwiseMax(wa2); diag = diag.cwiseMax(wa2);
/* beginning of the inner loop. */
do { do {
/* determine the levenberg-marquardt parameter. */ /* determine the levenberg-marquardt parameter. */
@ -374,9 +373,9 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
return XtolTooSmall; return XtolTooSmall;
if (gnorm <= epsilon<Scalar>()) if (gnorm <= epsilon<Scalar>())
return GtolTooSmall; return GtolTooSmall;
/* end of the inner loop. repeat if iteration unsuccessful. */
} while (ratio < Scalar(1e-4)); } while (ratio < Scalar(1e-4));
/* end of the outer loop. */
return Running; return Running;
} }
@ -468,7 +467,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
int rownb = 2; int rownb = 2;
for (i = 0; i < m; ++i) { for (i = 0; i < m; ++i) {
if (functor.df(x, wa3, rownb) < 0) return UserAsked; if (functor.df(x, wa3, rownb) < 0) return UserAsked;
ei_rwupdt<Scalar>(n, fjac.data(), fjac.rows(), wa3.data(), qtf.data(), fvec[i]); ei_rwupdt<Scalar>(fjac, wa3, qtf, fvec[i]);
++rownb; ++rownb;
} }
++njev; ++njev;
@ -485,7 +484,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
if (sing) { if (sing) {
wa2 = fjac.colwise().blueNorm(); wa2 = fjac.colwise().blueNorm();
// TODO We have no unit test covering this code path, do not modify // TODO We have no unit test covering this code path, do not modify
// before it is carefully tested // until it is carefully tested
ColPivHouseholderQR<JacobianType> qrfac(fjac); ColPivHouseholderQR<JacobianType> qrfac(fjac);
fjac = qrfac.matrixQR(); fjac = qrfac.matrixQR();
wa1 = fjac.diagonal(); wa1 = fjac.diagonal();
@ -538,7 +537,6 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
if (mode != 2) if (mode != 2)
diag = diag.cwiseMax(wa2); diag = diag.cwiseMax(wa2);
/* beginning of the inner loop. */
do { do {
/* determine the levenberg-marquardt parameter. */ /* determine the levenberg-marquardt parameter. */
@ -623,9 +621,9 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
return XtolTooSmall; return XtolTooSmall;
if (gnorm <= epsilon<Scalar>()) if (gnorm <= epsilon<Scalar>())
return GtolTooSmall; return GtolTooSmall;
/* end of the inner loop. repeat if iteration unsuccessful. */
} while (ratio < Scalar(1e-4)); } while (ratio < Scalar(1e-4));
/* end of the outer loop. */
return Running; return Running;
} }

View File

@ -1,45 +1,41 @@
template <typename Scalar> template <typename Scalar>
void ei_rwupdt(int n, Scalar *r__, int ldr, const Scalar *w, Scalar *b, Scalar alpha) void ei_rwupdt(
Matrix< Scalar, Dynamic, Dynamic > &r,
const Matrix< Scalar, Dynamic, 1> &w,
Matrix< Scalar, Dynamic, 1> &b,
Scalar alpha)
{ {
const int n = r.cols();
assert(r.rows()>=n);
std::vector<PlanarRotation<Scalar> > givens(n); std::vector<PlanarRotation<Scalar> > givens(n);
/* System generated locals */
int r_dim1, r_offset;
/* Local variables */ /* Local variables */
Scalar temp, rowj; Scalar temp, rowj;
/* Parameter adjustments */
--b;
--w;
r_dim1 = ldr;
r_offset = 1 + r_dim1 * 1;
r__ -= r_offset;
/* Function Body */ /* Function Body */
for (int j = 1; j <= n; ++j) { for (int j = 0; j < n; ++j) {
rowj = w[j]; rowj = w[j];
/* apply the previous transformations to */ /* apply the previous transformations to */
/* r(i,j), i=1,2,...,j-1, and to w(j). */ /* r(i,j), i=0,1,...,j-1, and to w(j). */
if (j-1>=1) for (int i = 0; i < j; ++i) {
for (int i = 1; i <= j-1; ++i) { temp = givens[i].c() * r(i,j) + givens[i].s() * rowj;
temp = givens[i-1].c() * r__[i + j * r_dim1] + givens[i-1].s() * rowj; rowj = -givens[i].s() * r(i,j) + givens[i].c() * rowj;
rowj = -givens[i-1].s() * r__[i + j * r_dim1] + givens[i-1].c() * rowj; r(i,j) = temp;
r__[i + j * r_dim1] = temp; }
}
if (rowj == 0.)
continue;
/* determine a givens rotation which eliminates w(j). */ /* determine a givens rotation which eliminates w(j). */
if (rowj != 0.) { givens[j].makeGivens(-r(j,j), rowj);
givens[j-1].makeGivens(-r__[j + j * r_dim1], rowj);
/* apply the current transformation to r(j,j), b(j), and alpha. */ /* apply the current transformation to r(j,j), b(j), and alpha. */
r__[j + j * r_dim1] = givens[j-1].c() * r__[j + j * r_dim1] + givens[j-1].s() * rowj; r(j,j) = givens[j].c() * r(j,j) + givens[j].s() * rowj;
temp = givens[j-1].c() * b[j] + givens[j-1].s() * alpha; temp = givens[j].c() * b[j] + givens[j].s() * alpha;
alpha = -givens[j-1].s() * b[j] + givens[j-1].c() * alpha; alpha = -givens[j].s() * b[j] + givens[j].c() * alpha;
b[j] = temp; b[j] = temp;
}
} }
return;
} }