use machine precision from eigen instead of the one from cminpack. The test

pass though they would not if using a value of 2.220e-16 (the real value
for 'double' is 2.22044604926e-16). How sensitive those tests are :)
This commit is contained in:
Thomas Capricelli 2009-08-19 19:56:51 +02:00
parent 3093e92da5
commit 01622e9855
7 changed files with 91 additions and 46 deletions

View File

@ -48,7 +48,78 @@ namespace Eigen {
#define p001 .001
#define p0001 1e-4
#if 1
#include <cminpack.h>
#else
/* Declarations for minpack */
typedef int (*minpack_func_nn)(void *p, int n, const double *x, double *fvec, int iflag );
typedef int (*minpack_funcder_nn)(void *p, int n, const double *x, double *fvec, double *fjac, int ldfjac, int iflag );
typedef int (*minpack_func_mn)(void *p, int m, int n, const double *x, double *fvec, int iflag );
typedef int (*minpack_funcder_mn)(void *p, int m, int n, const double *x, double *fvec, double *fjac, int ldfjac, int iflag );
typedef int (*minpack_funcderstr_mn)(void *p, int m, int n, const double *x, double *fvec, double *fjrow, int iflag );
/* MINPACK functions: */
int hybrd1 ( minpack_func_nn fcn, void *p, int n, double *x, double *fvec, double tol,
double *wa, int lwa );
int hybrd ( minpack_func_nn fcn,
void *p, int n, double *x, double *fvec, double xtol, int maxfev,
int ml, int mu, double epsfcn, double *diag, int mode,
double factor, int nprint, int *nfev,
double *fjac, int ldfjac, double *r, int lr, double *qtf,
double *wa1, double *wa2, double *wa3, double *wa4);
int hybrj1 ( minpack_funcder_nn fcn, void *p, int n, double *x,
double *fvec, double *fjac, int ldfjac, double tol,
double *wa, int lwa );
int hybrj ( minpack_funcder_nn fcn, void *p, int n, double *x,
double *fvec, double *fjac, int ldfjac, double xtol,
int maxfev, double *diag, int mode, double factor,
int nprint, int *nfev, int *njev, double *r,
int lr, double *qtf, double *wa1, double *wa2,
double *wa3, double *wa4 );
int lmdif1 ( minpack_func_mn fcn,
void *p, int m, int n, double *x, double *fvec, double tol,
int *iwa, double *wa, int lwa );
int lmdif ( minpack_func_mn fcn,
void *p, int m, int n, double *x, double *fvec, double ftol,
double xtol, double gtol, int maxfev, double epsfcn,
double *diag, int mode, double factor, int nprint,
int *nfev, double *fjac, int ldfjac, int *ipvt,
double *qtf, double *wa1, double *wa2, double *wa3,
double *wa4 );
int lmder1 ( minpack_funcder_mn fcn,
void *p, int m, int n, double *x, double *fvec, double *fjac,
int ldfjac, double tol, int *ipvt,
double *wa, int lwa );
int lmder ( minpack_funcder_mn fcn,
void *p, int m, int n, double *x, double *fvec, double *fjac,
int ldfjac, double ftol, double xtol, double gtol,
int maxfev, double *diag, int mode, double factor,
int nprint, int *nfev, int *njev, int *ipvt,
double *qtf, double *wa1, double *wa2, double *wa3,
double *wa4 );
int lmstr1 ( minpack_funcderstr_mn fcn, void *p, int m, int n,
double *x, double *fvec, double *fjac, int ldfjac,
double tol, int *ipvt, double *wa, int lwa );
int lmstr ( minpack_funcderstr_mn fcn, void *p, int m,
int n, double *x, double *fvec, double *fjac,
int ldfjac, double ftol, double xtol, double gtol,
int maxfev, double *diag, int mode, double factor,
int nprint, int *nfev, int *njev, int *ipvt,
double *qtf, double *wa1, double *wa2, double *wa3,
double *wa4 );
void chkder ( int m, int n, const double *x, double *fvec, double *fjac,
int ldfjac, double *xp, double *fvecp, int mode,
double *err );
#endif
#include "src/NonLinear/lmder1.h"
#include "src/NonLinear/lmder.h"
#include "src/NonLinear/hybrd1.h"

View File

@ -14,7 +14,7 @@ void chkder_template(int m, int n, const T *x,
/* Local variables */
int i__, j;
T eps, epsf, temp, epsmch;
T eps, epsf, temp;
T epslog;
/* Parameter adjustments */
@ -29,11 +29,7 @@ void chkder_template(int m, int n, const T *x,
/* Function Body */
/* epsmch is the machine precision. */
epsmch = dpmpar(1);
eps = sqrt(epsmch);
eps = ei_sqrt(epsilon<T>());
if (mode == 2) {
goto L20;
@ -56,7 +52,7 @@ L20:
/* mode = 2. */
epsf = chkder_factor * epsmch;
epsf = chkder_factor * epsilon<T>();
epslog = chkder_log10e * log(eps);
i__1 = m;
for (i__ = 1; i__ <= i__1; ++i__) {
@ -87,7 +83,7 @@ L20:
fabs(fvecp[i__]));
}
err[i__] = 1.;
if (temp > epsmch && temp < eps) {
if (temp > epsilon<T>() && temp < eps) {
err[i__] = (chkder_log10e * log(temp) - epslog) / epslog;
}
if (temp >= eps) {

View File

@ -28,7 +28,7 @@ int hybrd_template(minpack_func_nn fcn, void *p, int n, T *x, T *
T pnorm, xnorm, fnorm1;
int nslow1, nslow2;
int ncfail;
T actred, epsmch, prered;
T actred, prered;
int info;
/* Parameter adjustments */
@ -47,10 +47,6 @@ int hybrd_template(minpack_func_nn fcn, void *p, int n, T *x, T *
/* Function Body */
/* epsmch is the machine precision. */
epsmch = dpmpar(1);
info = 0;
iflag = 0;
*nfev = 0;
@ -382,7 +378,7 @@ L260:
}
/* Computing MAX */
d__1 = p1 * delta;
if (p1 * max(d__1,pnorm) <= epsmch * xnorm) {
if (p1 * max(d__1,pnorm) <= epsilon<T>() * xnorm) {
info = 3;
}
if (nslow2 == 5) {

View File

@ -28,7 +28,7 @@ int hybrj_template(minpack_funcder_nn fcn, void *p, int n, T *x, T *
T pnorm, xnorm, fnorm1;
int nslow1, nslow2;
int ncfail;
T actred, epsmch, prered;
T actred, prered;
int info;
/* Parameter adjustments */
@ -47,10 +47,6 @@ int hybrj_template(minpack_funcder_nn fcn, void *p, int n, T *x, T *
/* Function Body */
/* epsmch is the machine precision. */
epsmch = dpmpar(1);
info = 0;
iflag = 0;
*nfev = 0;
@ -375,7 +371,7 @@ L260:
}
/* Computing MAX */
d__1 = p1 * delta;
if (p1 * max(d__1,pnorm) <= epsmch * xnorm) {
if (p1 * max(d__1,pnorm) <= epsilon<T>() * xnorm) {
info = 3;
}
if (nslow2 == 5) {

View File

@ -21,8 +21,7 @@ int lmder_template(minpack_funcder_mn fcn, void *p, int m, int n, T *x,
int iflag;
T delta;
T ratio;
T fnorm, gnorm, pnorm, xnorm, fnorm1, actred, dirder,
epsmch, prered;
T fnorm, gnorm, pnorm, xnorm, fnorm1, actred, dirder, prered;
int info;
/* Parameter adjustments */
@ -41,10 +40,6 @@ int lmder_template(minpack_funcder_mn fcn, void *p, int m, int n, T *x,
/* Function Body */
/* epsmch is the machine precision. */
epsmch = dpmpar(1);
info = 0;
iflag = 0;
*nfev = 0;
@ -382,13 +377,13 @@ L290:
if (*nfev >= maxfev) {
info = 5;
}
if (fabs(actred) <= epsmch && prered <= epsmch && p5 * ratio <= 1.) {
if (fabs(actred) <= epsilon<T>() && prered <= epsilon<T>() && p5 * ratio <= 1.) {
info = 6;
}
if (delta <= epsmch * xnorm) {
if (delta <= epsilon<T>() * xnorm) {
info = 7;
}
if (gnorm <= epsmch) {
if (gnorm <= epsilon<T>()) {
info = 8;
}
if (info != 0) {

View File

@ -24,7 +24,7 @@ int lmdif_template(minpack_func_mn fcn, void *p, int m, int n, T *x,
T delta;
T ratio;
T fnorm, gnorm;
T pnorm, xnorm, fnorm1, actred, dirder, epsmch, prered;
T pnorm, xnorm, fnorm1, actred, dirder, prered;
int info;
/* Parameter adjustments */
@ -43,10 +43,6 @@ int lmdif_template(minpack_func_mn fcn, void *p, int m, int n, T *x,
/* Function Body */
/* epsmch is the machine precision. */
epsmch = dpmpar(1);
info = 0;
iflag = 0;
*nfev = 0;
@ -384,13 +380,13 @@ L290:
if (*nfev >= maxfev) {
info = 5;
}
if (fabs(actred) <= epsmch && prered <= epsmch && p5 * ratio <= 1.) {
if (fabs(actred) <= epsilon<T>() && prered <= epsilon<T>() && p5 * ratio <= 1.) {
info = 6;
}
if (delta <= epsmch * xnorm) {
if (delta <= epsilon<T>() * xnorm) {
info = 7;
}
if (gnorm <= epsmch) {
if (gnorm <= epsilon<T>()) {
info = 8;
}
if (info != 0) {

View File

@ -22,8 +22,7 @@ int lmstr_template(minpack_funcderstr_mn fcn, void *p, int m, int n, T *x,
int iflag;
T delta;
T ratio;
T fnorm, gnorm, pnorm, xnorm, fnorm1, actred, dirder,
epsmch, prered;
T fnorm, gnorm, pnorm, xnorm, fnorm1, actred, dirder, prered;
int info;
/* Parameter adjustments */
@ -42,10 +41,6 @@ int lmstr_template(minpack_funcderstr_mn fcn, void *p, int m, int n, T *x,
/* Function Body */
/* epsmch is the machine precision. */
epsmch = dpmpar(1);
info = 0;
iflag = 0;
*nfev = 0;
@ -409,13 +404,13 @@ L330:
if (*nfev >= maxfev) {
info = 5;
}
if (abs(actred) <= epsmch && prered <= epsmch && p5 * ratio <= 1.) {
if (abs(actred) <= epsilon<T>() && prered <= epsilon<T>() && p5 * ratio <= 1.) {
info = 6;
}
if (delta <= epsmch * xnorm) {
if (delta <= epsilon<T>() * xnorm) {
info = 7;
}
if (gnorm <= epsmch) {
if (gnorm <= epsilon<T>()) {
info = 8;
}
if (info != 0) {