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

merging ei_lmder and lmder_template into ei_lmder() which takes eigen

Browse Source 
argument, but still uses f2c code inside.
This commit is contained in:
Thomas Capricelli 2009-08-21 02:34:40 +02:00
parent 6a8b52b3aa
commit 1715e2cb3b
3 changed files with 63 additions and 106 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
unsupported/Eigen/src/NonLinear/MathFunctions.h
View File

@ -159,49 +159,6 @@ int ei_lmstr(
); );
} }
template<typename Functor, typename Scalar>
int ei_lmder(
Matrix< Scalar, Dynamic, 1 > &x,
Matrix< Scalar, Dynamic, 1 > &fvec,
int &nfev,
int &njev,
Matrix< Scalar, Dynamic, Dynamic > &fjac,
VectorXi &ipvt,
Matrix< Scalar, Dynamic, 1 > &diag,
int mode=1,
Scalar factor = 100.,
int maxfev = 400,
Scalar ftol = ei_sqrt(epsilon<Scalar>()),
Scalar xtol = ei_sqrt(epsilon<Scalar>()),
Scalar gtol = Scalar(0.),
int nprint=0
)
{
Matrix< Scalar, Dynamic, 1 >
qtf(x.size()),
wa1(x.size()), wa2(x.size()), wa3(x.size()),
wa4(fvec.size());
int ldfjac = fvec.size();
ipvt.resize(x.size());
fjac.resize(ldfjac, x.size());
diag.resize(x.size());
return lmder_template<Scalar>(
Functor::f, 0,
fvec.size(), x.size(), x.data(), fvec.data(),
fjac.data() , ldfjac,
ftol, xtol, gtol,
maxfev,
diag.data(), mode,
factor,
nprint,
nfev, njev,
ipvt.data(),
qtf.data(),
wa1.data(), wa2.data(), wa3.data(), wa4.data()
);
}
template<typename Functor, typename Scalar> template<typename Functor, typename Scalar>
int ei_lmdif( int ei_lmdif(
Matrix< Scalar, Dynamic, 1 > &x, Matrix< Scalar, Dynamic, 1 > &x,

124
unsupported/Eigen/src/NonLinear/lmder.h
View File

@ -1,16 +1,29 @@
template<typename Scalar> template<typename Functor, typename Scalar>
int lmder_template(minpack_funcder_mn fcn, void *p, int m, int n, Scalar *x, int ei_lmder(
Scalar *fvec, Scalar *fjac, int ldfjac, Scalar ftol, Matrix< Scalar, Dynamic, 1 > &x,
Scalar xtol, Scalar gtol, int maxfev, Scalar * Matrix< Scalar, Dynamic, 1 > &fvec,
diag, int mode, Scalar factor, int nprint, int &nfev,
int &nfev, int &njev, int *ipvt, Scalar *qtf, int &njev,
Scalar *wa1, Scalar *wa2, Scalar *wa3, Scalar *wa4) Matrix< Scalar, Dynamic, Dynamic > &fjac,
VectorXi &ipvt,
Matrix< Scalar, Dynamic, 1 > &diag,
int mode=1,
Scalar factor = 100.,
int maxfev = 400,
Scalar ftol = ei_sqrt(epsilon<Scalar>()),
Scalar xtol = ei_sqrt(epsilon<Scalar>()),
Scalar gtol = Scalar(0.),
int nprint=0
)
{ {
/* Initialized data */ const int m = fvec.size(), n = x.size();
Matrix< Scalar, Dynamic, 1 > qtf(n), wa1(n), wa2(n), wa3(n), wa4(m);
int ldfjac = m;
/* System generated locals */ ipvt.resize(n);
int fjac_offset; fjac.resize(ldfjac, n);
diag.resize(n);
/* Local variables */ /* Local variables */
int i, j, l; int i, j, l;
@ -23,18 +36,6 @@ int lmder_template(minpack_funcder_mn fcn, void *p, int m, int n, Scalar *x,
Scalar fnorm, gnorm, pnorm, xnorm, fnorm1, actred, dirder, prered; Scalar fnorm, gnorm, pnorm, xnorm, fnorm1, actred, dirder, prered;
int info; int info;
/* Parameter adjustments */
--wa4;
--fvec;
--wa3;
--wa2;
--wa1;
--qtf;
--ipvt;
--diag;
--x;
fjac_offset = 1 + ldfjac;
fjac -= fjac_offset;
/* Function Body */ /* Function Body */
@ -52,7 +53,7 @@ int lmder_template(minpack_funcder_mn fcn, void *p, int m, int n, Scalar *x,
if (mode != 2) { if (mode != 2) {
goto L20; goto L20;
} }
for (j = 1; j <= n; ++j) { for (j = 0; j < n; ++j) {
if (diag[j] <= 0.) { if (diag[j] <= 0.) {
goto L300; goto L300;
} }
@ -63,12 +64,12 @@ L20:
/* evaluate the function at the starting point */ /* evaluate the function at the starting point */
/* and calculate its norm. */ /* and calculate its norm. */
iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 1); iflag = Functor::f(0, m, n, x.data(), fvec.data(), fjac.data(), ldfjac, 1);
nfev = 1; nfev = 1;
if (iflag < 0) { if (iflag < 0) {
goto L300; goto L300;
} }
fnorm = ei_enorm<Scalar>(m, &fvec[1]); fnorm = fvec.stableNorm();
/* initialize levenberg-marquardt parameter and iteration counter. */ /* initialize levenberg-marquardt parameter and iteration counter. */
@ -81,20 +82,20 @@ L30:
/* calculate the jacobian matrix. */ /* calculate the jacobian matrix. */
iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 2); iflag = Functor::f(0, m, n, x.data(), fvec.data(), fjac.data(), ldfjac, 2);
++njev; ++njev;
if (iflag < 0) { if (iflag < 0) {
goto L300; goto L300;
} }
/* if requested, call fcn to enable printing of iterates. */ /* if requested, call Functor::f to enable printing of iterates. */
if (nprint <= 0) { if (nprint <= 0) {
goto L40; goto L40;
} }
iflag = 0; iflag = 0;
if ((iter - 1) % nprint == 0) { if ((iter - 1) % nprint == 0) {
iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 0); iflag = Functor::f(0, m, n, x.data(), fvec.data(), fjac.data(), ldfjac, 0);
} }
if (iflag < 0) { if (iflag < 0) {
goto L300; goto L300;
@ -103,8 +104,7 @@ L40:
/* compute the qr factorization of the jacobian. */ /* compute the qr factorization of the jacobian. */
qrfac(m, n, &fjac[fjac_offset], ldfjac, TRUE_, &ipvt[1], n, &wa1[1], & qrfac(m, n, fjac.data(), ldfjac, TRUE_, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data());
wa2[1], &wa3[1]);
/* 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. */
@ -115,7 +115,7 @@ L40:
if (mode == 2) { if (mode == 2) {
goto L60; goto L60;
} }
for (j = 1; 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.;
@ -127,11 +127,11 @@ L60:
/* 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. */
for (j = 1; j <= n; ++j) { for (j = 0; j < n; ++j) {
wa3[j] = diag[j] * x[j]; wa3[j] = diag[j] * x[j];
/* L70: */ /* L70: */
} }
xnorm = ei_enorm<Scalar>(n, &wa3[1]); xnorm = wa3.stableNorm();
delta = factor * xnorm; delta = factor * xnorm;
if (delta == 0.) { if (delta == 0.) {
delta = factor; delta = factor;
@ -141,26 +141,26 @@ L80:
/* form (q transpose)*fvec and store the first n components in */ /* form (q transpose)*fvec and store the first n components in */
/* qtf. */ /* qtf. */
for (i = 1; i <= m; ++i) { for (i = 0; i < m; ++i) {
wa4[i] = fvec[i]; wa4[i] = fvec[i];
/* L90: */ /* L90: */
} }
for (j = 1; 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: */
} }
@ -171,14 +171,14 @@ L120:
if (fnorm == 0.) { if (fnorm == 0.) {
goto L170; goto L170;
} }
for (j = 1; j <= n; ++j) { for (j = 0; j < n; ++j) {
l = ipvt[j]; l = ipvt[j]-1;
if (wa2[l] == 0.) { if (wa2[l] == 0.) {
goto L150; goto L150;
} }
sum = 0.; sum = 0.;
for (i = 1; 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: */ /* L140: */
} }
/* Computing MAX */ /* Computing MAX */
@ -203,7 +203,7 @@ L170:
if (mode == 2) { if (mode == 2) {
goto L190; goto L190;
} }
for (j = 1; j <= n; ++j) for (j = 0; j < n; ++j)
diag[j] = max( diag[j], wa2[j]); diag[j] = max( diag[j], wa2[j]);
L190: L190:
@ -213,18 +213,18 @@ L200:
/* determine the levenberg-marquardt parameter. */ /* determine the levenberg-marquardt parameter. */
lmpar(n, &fjac[fjac_offset], ldfjac, &ipvt[1], &diag[1], &qtf[1], delta, lmpar(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta,
&par, &wa1[1], &wa2[1], &wa3[1], &wa4[1]); &par, wa1.data(), wa2.data(), wa3.data(), wa4.data());
/* store the direction p and x + p. calculate the norm of p. */ /* store the direction p and x + p. calculate the norm of p. */
for (j = 1; j <= n; ++j) { for (j = 0; j < n; ++j) {
wa1[j] = -wa1[j]; wa1[j] = -wa1[j];
wa2[j] = x[j] + wa1[j]; wa2[j] = x[j] + wa1[j];
wa3[j] = diag[j] * wa1[j]; wa3[j] = diag[j] * wa1[j];
/* L210: */ /* L210: */
} }
pnorm = ei_enorm<Scalar>(n, &wa3[1]); pnorm = wa3.stableNorm();
/* on the first iteration, adjust the initial step bound. */ /* on the first iteration, adjust the initial step bound. */
@ -234,12 +234,12 @@ L200:
/* evaluate the function at x + p and calculate its norm. */ /* evaluate the function at x + p and calculate its norm. */
iflag = (*fcn)(p, m, n, &wa2[1], &wa4[1], &fjac[fjac_offset], ldfjac, 1); iflag = Functor::f(0, m, n, wa2.data(), wa4.data(), fjac.data(), ldfjac, 1);
++nfev; ++nfev;
if (iflag < 0) { if (iflag < 0) {
goto L300; goto L300;
} }
fnorm1 = ei_enorm<Scalar>(m, &wa4[1]); fnorm1 = wa4.stableNorm();
/* compute the scaled actual reduction. */ /* compute the scaled actual reduction. */
@ -250,17 +250,17 @@ L200:
/* compute the scaled predicted reduction and */ /* compute the scaled predicted reduction and */
/* the scaled directional derivative. */ /* the scaled directional derivative. */
for (j = 1; j <= n; ++j) { for (j = 0; j < n; ++j) {
wa3[j] = 0.; wa3[j] = 0.;
l = ipvt[j]; l = ipvt[j]-1;
temp = wa1[l]; temp = wa1[l];
for (i = 1; 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: */
} }
temp1 = ei_abs2(ei_enorm<Scalar>(n, &wa3[1]) / fnorm); 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 */ /* Computing 2nd power */
prered = temp1 + temp2 / p5; prered = temp1 + temp2 / p5;
@ -309,16 +309,16 @@ L260:
/* successful iteration. update x, fvec, and their norms. */ /* successful iteration. update x, fvec, and their norms. */
for (j = 1; j <= n; ++j) { for (j = 0; j < n; ++j) {
x[j] = wa2[j]; x[j] = wa2[j];
wa2[j] = diag[j] * x[j]; wa2[j] = diag[j] * x[j];
/* L270: */ /* L270: */
} }
for (i = 1; i <= m; ++i) { for (i = 0; i < m; ++i) {
fvec[i] = wa4[i]; fvec[i] = wa4[i];
/* L280: */ /* L280: */
} }
xnorm = ei_enorm<Scalar>(n, &wa2[1]); xnorm = wa2.stableNorm();
fnorm = fnorm1; fnorm = fnorm1;
++iter; ++iter;
L290: L290:
@ -375,7 +375,7 @@ L300:
} }
iflag = 0; iflag = 0;
if (nprint > 0) { if (nprint > 0) {
iflag = (*fcn)(p, m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, 0); iflag = Functor::f(0, m, n, x.data(), fvec.data(), fjac.data(), ldfjac, 0);
} }
return info; return info;

2
unsupported/test/NonLinear.cpp
View File

@ -1925,7 +1925,7 @@ void test_NonLinear()
/* /*
* Can be useful for debugging... * Can be useful for debugging...
printf("info, nfev, jev : %d, %d, %d\n", info, nfev, njev); printf("info, nfev, njev : %d, %d, %d\n", info, nfev, njev);
printf("x[0] : %.32g\n", x[0]); printf("x[0] : %.32g\n", x[0]);
printf("x[1] : %.32g\n", x[1]); printf("x[1] : %.32g\n", x[1]);
printf("x[2] : %.32g\n", x[2]); printf("x[2] : %.32g\n", x[2]);