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

use eigen stableNorm() instead of cminpack 'enorm'. The results are mostly

Browse Source 
slightly better in tests (one test needs 15 iterations intead of 16, for
the same result). Some numerical results have improved slightly, too.

If one uses blueNorm() instead, an assert for 'overflow' is raised from
blueNorm()
This commit is contained in:
Thomas Capricelli 2009-08-20 21:04:38 +02:00
parent 6953cad81d
commit 9a876806e1
7 changed files with 59 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
unsupported/Eigen/NonLinear
View File

@ -116,9 +116,16 @@ int lmstr ( minpack_funcderstr_mn fcn, void *p, int m,
void chkder ( int m, int n, const double *x, double *fvec, double *fjac,
int ldfjac, double *xp, double *fvecp, int mode,
double *err );
void covar(int n, double *r__, int ldr, const int *ipvt, double tol, double *wa);
#endif
template<typename Scalar>
Scalar ei_enorm ( int n, const Scalar *x ){
return Map< Matrix< Scalar, Dynamic, 1 > >(x,n).stableNorm();
// return Map< Matrix< Scalar, Dynamic, 1 > >(x,n).blueNorm();
}
#include "src/NonLinear/lmder1.h"
#include "src/NonLinear/lmder.h"

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

@ -77,7 +77,7 @@ L20:
if (iflag < 0) {
goto L300;
}
fnorm = enorm(n, &fvec[1]);
fnorm = ei_enorm<T>(n, &fvec[1]);
/* determine the number of calls to fcn needed to compute */
/* the jacobian matrix. */
@ -140,7 +140,7 @@ L50:
wa3[j] = diag[j] * x[j];
/* L60: */
}
xnorm = enorm(n, &wa3[1]);
xnorm = ei_enorm<T>(n, &wa3[1]);
delta = factor * xnorm;
if (delta == 0.) {
delta = factor;
@ -250,7 +250,7 @@ L190:
wa3[j] = diag[j] * wa1[j];
/* L200: */
}
pnorm = enorm(n, &wa3[1]);
pnorm = ei_enorm<T>(n, &wa3[1]);
/* on the first iteration, adjust the initial step bound. */
@ -265,7 +265,7 @@ L190:
if (iflag < 0) {
goto L300;
}
fnorm1 = enorm(n, &wa4[1]);
fnorm1 = ei_enorm<T>(n, &wa4[1]);
/* compute the scaled actual reduction. */
@ -291,7 +291,7 @@ L190:
wa3[i__] = qtf[i__] + sum;
/* L220: */
}
temp = enorm(n, &wa3[1]);
temp = ei_enorm<T>(n, &wa3[1]);
prered = 0.;
if (temp < fnorm) {
/* Computing 2nd power */
@ -344,7 +344,7 @@ L240:
fvec[j] = wa4[j];
/* L250: */
}
xnorm = enorm(n, &wa2[1]);
xnorm = ei_enorm<T>(n, &wa2[1]);
fnorm = fnorm1;
++iter;
L260:

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

@ -78,7 +78,7 @@ L20:
if (iflag < 0) {
goto L300;
}
fnorm = enorm(n, &fvec[1]);
fnorm = ei_enorm<T>(n, &fvec[1]);
/* initialize iteration counter and monitors. */
@ -133,7 +133,7 @@ L50:
wa3[j] = diag[j] * x[j];
/* L60: */
}
xnorm = enorm(n, &wa3[1]);
xnorm = ei_enorm<T>(n, &wa3[1]);
delta = factor * xnorm;
if (delta == 0.) {
delta = factor;
@ -243,7 +243,7 @@ L190:
wa3[j] = diag[j] * wa1[j];
/* L200: */
}
pnorm = enorm(n, &wa3[1]);
pnorm = ei_enorm<T>(n, &wa3[1]);
/* on the first iteration, adjust the initial step bound. */
@ -258,7 +258,7 @@ L190:
if (iflag < 0) {
goto L300;
}
fnorm1 = enorm(n, &wa4[1]);
fnorm1 = ei_enorm<T>(n, &wa4[1]);
/* compute the scaled actual reduction. */
@ -284,7 +284,7 @@ L190:
wa3[i__] = qtf[i__] + sum;
/* L220: */
}
temp = enorm(n, &wa3[1]);
temp = ei_enorm<T>(n, &wa3[1]);
prered = 0.;
if (temp < fnorm) {
/* Computing 2nd power */
@ -337,7 +337,7 @@ L240:
fvec[j] = wa4[j];
/* L250: */
}
xnorm = enorm(n, &wa2[1]);
xnorm = ei_enorm<T>(n, &wa2[1]);
fnorm = fnorm1;
++iter;
L260:

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

@ -71,7 +71,7 @@ L20:
if (iflag < 0) {
goto L300;
}
fnorm = enorm(m, &fvec[1]);
fnorm = ei_enorm<T>(m, &fvec[1]);
/* initialize levenberg-marquardt parameter and iteration counter. */
@ -136,7 +136,7 @@ L60:
wa3[j] = diag[j] * x[j];
/* L70: */
}
xnorm = enorm(n, &wa3[1]);
xnorm = ei_enorm<T>(n, &wa3[1]);
delta = factor * xnorm;
if (delta == 0.) {
delta = factor;
@ -242,7 +242,7 @@ L200:
wa3[j] = diag[j] * wa1[j];
/* L210: */
}
pnorm = enorm(n, &wa3[1]);
pnorm = ei_enorm<T>(n, &wa3[1]);
/* on the first iteration, adjust the initial step bound. */
@ -257,7 +257,7 @@ L200:
if (iflag < 0) {
goto L300;
}
fnorm1 = enorm(m, &wa4[1]);
fnorm1 = ei_enorm<T>(m, &wa4[1]);
/* compute the scaled actual reduction. */
@ -283,7 +283,7 @@ L200:
}
/* L230: */
}
temp1 = enorm(n, &wa3[1]) / fnorm;
temp1 = ei_enorm<T>(n, &wa3[1]) / fnorm;
temp2 = sqrt(par) * pnorm / fnorm;
/* Computing 2nd power */
d__1 = temp1;
@ -351,7 +351,7 @@ L260:
fvec[i__] = wa4[i__];
/* L280: */
}
xnorm = enorm(n, &wa2[1]);
xnorm = ei_enorm<T>(n, &wa2[1]);
fnorm = fnorm1;
++iter;
L290:

12
unsupported/Eigen/src/NonLinear/lmdif.h
View File

@ -73,7 +73,7 @@ L20:
if (iflag < 0) {
goto L300;
}
fnorm = enorm(m, &fvec[1]);
fnorm = ei_enorm<T>(m, &fvec[1]);
/* initialize levenberg-marquardt parameter and iteration counter. */
@ -139,7 +139,7 @@ L60:
wa3[j] = diag[j] * x[j];
/* L70: */
}
xnorm = enorm(n, &wa3[1]);
xnorm = ei_enorm<T>(n, &wa3[1]);
delta = factor * xnorm;
if (delta == 0.) {
delta = factor;
@ -245,7 +245,7 @@ L200:
wa3[j] = diag[j] * wa1[j];
/* L210: */
}
pnorm = enorm(n, &wa3[1]);
pnorm = ei_enorm<T>(n, &wa3[1]);
/* on the first iteration, adjust the initial step bound. */
@ -260,7 +260,7 @@ L200:
if (iflag < 0) {
goto L300;
}
fnorm1 = enorm(m, &wa4[1]);
fnorm1 = ei_enorm<T>(m, &wa4[1]);
/* compute the scaled actual reduction. */
@ -286,7 +286,7 @@ L200:
}
/* L230: */
}
temp1 = enorm(n, &wa3[1]) / fnorm;
temp1 = ei_enorm<T>(n, &wa3[1]) / fnorm;
temp2 = sqrt(par) * pnorm / fnorm;
/* Computing 2nd power */
d__1 = temp1;
@ -354,7 +354,7 @@ L260:
fvec[i__] = wa4[i__];
/* L280: */
}
xnorm = enorm(n, &wa2[1]);
xnorm = ei_enorm<T>(n, &wa2[1]);
fnorm = fnorm1;
++iter;
L290:

14
unsupported/Eigen/src/NonLinear/lmstr.h
View File

@ -72,7 +72,7 @@ L20:
if (iflag < 0) {
goto L340;
}
fnorm = enorm(m, &fvec[1]);
fnorm = ei_enorm<T>(m, &fvec[1]);
/* initialize levenberg-marquardt parameter and iteration counter. */
@ -136,7 +136,7 @@ L40:
sing = TRUE_;
}
ipvt[j] = j;
wa2[j] = enorm(j, &fjac[j * fjac_dim1 + 1]);
wa2[j] = ei_enorm<T>(j, &fjac[j * fjac_dim1 + 1]);
/* L80: */
}
if (! sing) {
@ -194,7 +194,7 @@ L150:
wa3[j] = diag[j] * x[j];
/* L160: */
}
xnorm = enorm(n, &wa3[1]);
xnorm = ei_enorm<T>(n, &wa3[1]);
delta = factor * xnorm;
if (delta == 0.) {
delta = factor;
@ -269,7 +269,7 @@ L240:
wa3[j] = diag[j] * wa1[j];
/* L250: */
}
pnorm = enorm(n, &wa3[1]);
pnorm = ei_enorm<T>(n, &wa3[1]);
/* on the first iteration, adjust the initial step bound. */
@ -284,7 +284,7 @@ L240:
if (iflag < 0) {
goto L340;
}
fnorm1 = enorm(m, &wa4[1]);
fnorm1 = ei_enorm<T>(m, &wa4[1]);
/* compute the scaled actual reduction. */
@ -310,7 +310,7 @@ L240:
}
/* L270: */
}
temp1 = enorm(n, &wa3[1]) / fnorm;
temp1 = ei_enorm<T>(n, &wa3[1]) / fnorm;
temp2 = sqrt(par) * pnorm / fnorm;
/* Computing 2nd power */
d__1 = temp1;
@ -378,7 +378,7 @@ L300:
fvec[i__] = wa4[i__];
/* L320: */
}
xnorm = enorm(n, &wa2[1]);
xnorm = ei_enorm<T>(n, &wa2[1]);
fnorm = fnorm1;
++iter;
L330:

32
unsupported/test/NonLinear.cpp
View File

@ -1142,7 +1142,7 @@ void testNistLanczos1(void)
VERIFY( 79 == nfev);
VERIFY( 72 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 1.4292388868910E-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats
VERIFY_IS_APPROX(fvec.squaredNorm(), 1.429604433690E-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats
// check x
VERIFY_IS_APPROX(x[0], 9.5100000027E-02 );
VERIFY_IS_APPROX(x[1], 1.0000000001E+00 );
@ -1294,7 +1294,7 @@ void testNistMGH10(void)
info = ei_lmder<MGH10_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 3 == info);
VERIFY( 2 == info);
VERIFY( 285 == nfev);
VERIFY( 250 == njev);
// check norm^2
@ -1312,7 +1312,7 @@ void testNistMGH10(void)
info = ei_lmder<MGH10_functor, double>(x, fvec, nfev, njev, fjac, ipvt, diag);
// check return value
VERIFY( 3 == info);
VERIFY( 2 == info);
VERIFY( 126 == nfev);
VERIFY( 116 == njev);
// check norm^2
@ -1388,8 +1388,8 @@ void testNistBoxBOD(void)
// check return value
VERIFY( 1 == info);
VERIFY( 16 == nfev);
VERIFY( 15 == njev);
VERIFY( 15 == nfev);
VERIFY( 14 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 1.1680088766E+03);
// check x
@ -1531,15 +1531,15 @@ void testNistMGH09(void)
// check return value
VERIFY( 1 == info);
VERIFY( 487== nfev);
VERIFY( 378 == njev);
VERIFY( 503== nfev);
VERIFY( 385 == njev);
// check norm^2
VERIFY_IS_APPROX(fvec.squaredNorm(), 3.0750560385E-04);
// check x
VERIFY_IS_APPROX(x[0], 0.19280590); // should be 1.9280693458E-01
VERIFY_IS_APPROX(x[1], 0.19130543); // should be 1.9128232873E-01
VERIFY_IS_APPROX(x[2], 0.12306085); // should be 1.2305650693E-01
VERIFY_IS_APPROX(x[3], 0.13607303); // should be 1.3606233068E-01
VERIFY_IS_APPROX(x[0], 0.19280624); // should be 1.9280693458E-01
VERIFY_IS_APPROX(x[1], 0.19129774); // should be 1.9128232873E-01
VERIFY_IS_APPROX(x[2], 0.12305940); // should be 1.2305650693E-01
VERIFY_IS_APPROX(x[3], 0.13606946); // should be 1.3606233068E-01
/*
* Second try
@ -1924,3 +1924,13 @@ void test_NonLinear()
CALL_SUBTEST(testNistEckerle4());
}
/*
* Can be useful for debugging...
printf("info, nfev, jev : %d, %d, %d\n", info, nfev, njev);
printf("x[0] : %.32g\n", x[0]);
printf("x[1] : %.32g\n", x[1]);
printf("x[2] : %.32g\n", x[2]);
printf("x[3] : %.32g\n", x[3]);
printf("fvec.squaredNorm() : %.32g\n", fvec.squaredNorm());
*/