* porting lmdif1 to eigen

* qtf was missing in lmdif signature (this is an output of the method)
This commit is contained in:
Thomas Capricelli 2009-08-20 22:59:09 +02:00
parent 8d2f6ad7e1
commit b09ebe01da
4 changed files with 31 additions and 78 deletions

View File

@ -130,7 +130,6 @@ Scalar ei_enorm ( int n, const Scalar *x ){
#include "src/NonLinear/lmder.h" #include "src/NonLinear/lmder.h"
#include "src/NonLinear/hybrd.h" #include "src/NonLinear/hybrd.h"
#include "src/NonLinear/lmstr.h" #include "src/NonLinear/lmstr.h"
#include "src/NonLinear/lmdif1.h"
#include "src/NonLinear/lmdif.h" #include "src/NonLinear/lmdif.h"
#include "src/NonLinear/hybrj.h" #include "src/NonLinear/hybrj.h"
#include "src/NonLinear/chkder.h" #include "src/NonLinear/chkder.h"
@ -139,6 +138,7 @@ Scalar ei_enorm ( int n, const Scalar *x ){
#include "src/NonLinear/lmstr1.h" #include "src/NonLinear/lmstr1.h"
#include "src/NonLinear/hybrd1.h" #include "src/NonLinear/hybrd1.h"
#include "src/NonLinear/hybrj1.h" #include "src/NonLinear/hybrj1.h"
#include "src/NonLinear/lmdif1.h"
//@} //@}

View File

@ -209,6 +209,7 @@ int ei_lmdif(
int &nfev, int &nfev,
Matrix< Scalar, Dynamic, Dynamic > &fjac, Matrix< Scalar, Dynamic, Dynamic > &fjac,
VectorXi &ipvt, VectorXi &ipvt,
Matrix< Scalar, Dynamic, 1 > &qtf,
Matrix< Scalar, Dynamic, 1 > &diag, Matrix< Scalar, Dynamic, 1 > &diag,
int mode=1, int mode=1,
Scalar factor = 100., Scalar factor = 100.,
@ -221,7 +222,6 @@ int ei_lmdif(
) )
{ {
Matrix< Scalar, Dynamic, 1 > Matrix< Scalar, Dynamic, 1 >
qtf(x.size()),
wa1(x.size()), wa2(x.size()), wa3(x.size()), wa1(x.size()), wa2(x.size()), wa3(x.size()),
wa4(fvec.size()); wa4(fvec.size());
int ldfjac = fvec.size(); int ldfjac = fvec.size();
@ -229,6 +229,7 @@ int ei_lmdif(
ipvt.resize(x.size()); ipvt.resize(x.size());
fjac.resize(ldfjac, x.size()); fjac.resize(ldfjac, x.size());
diag.resize(x.size()); diag.resize(x.size());
qtf.resize(x.size());
return lmdif_template<Scalar> ( return lmdif_template<Scalar> (
Functor::f, 0, Functor::f, 0,
fvec.size(), x.size(), x.data(), fvec.data(), fvec.size(), x.size(), x.data(), fvec.data(),
@ -246,30 +247,6 @@ int ei_lmdif(
); );
} }
template<typename Functor, typename Scalar>
int ei_lmdif1(
Matrix< Scalar, Dynamic, 1 > &x,
Matrix< Scalar, Dynamic, 1 > &fvec,
Scalar tol = ei_sqrt(epsilon<Scalar>())
)
{
int n = x.size();
int ldfjac = fvec.size();
int lwa = ldfjac*n+5*n+ldfjac;
VectorXi iwa(n);
Matrix< Scalar, Dynamic, 1 > wa(lwa);
Matrix< Scalar, Dynamic, Dynamic > fjac(ldfjac, n);
wa.resize(lwa);
return lmdif1_template<Scalar> (
Functor::f, 0,
fvec.size(), n, x.data(), fvec.data(),
tol,
iwa.data(),
wa.data(), lwa
);
}
template<typename Scalar> template<typename Scalar>
void ei_chkder( void ei_chkder(
Matrix< Scalar, Dynamic, 1 > &x, Matrix< Scalar, Dynamic, 1 > &x,

View File

@ -1,56 +1,32 @@
template<typename Scalar> template<typename Functor, typename Scalar>
int lmdif1_template(minpack_func_mn fcn, void *p, int m, int n, Scalar *x, int ei_lmdif1(
Scalar *fvec, Scalar tol, int *iwa, Matrix< Scalar, Dynamic, 1 > &x,
Scalar *wa, int lwa) Matrix< Scalar, Dynamic, 1 > &fvec,
Scalar tol = ei_sqrt(epsilon<Scalar>())
)
{ {
/* Initialized data */ const int n = x.size(), m=fvec.size();
int info, nfev;
Matrix< Scalar, Dynamic, Dynamic > fjac(m, n);
Matrix< Scalar, Dynamic, 1> diag, qtf;
VectorXi ipvt;
const Scalar factor = 100.; /* check the input parameters for errors. */
if (n <= 0 || m < n || tol < 0.) {
int mp5n, mode, nfev; printf("ei_lmder1 bad args : m,n,tol,...");
Scalar ftol, gtol, xtol; return 0;
Scalar epsfcn;
int maxfev, nprint;
int info;
/* Parameter adjustments */
--fvec;
--iwa;
--x;
--wa;
/* Function Body */
info = 0;
/* check the input parameters for errors. */
if (n <= 0 || m < n || tol < 0. || lwa < m * n + n * 5 + m) {
/* goto L10; */
return info;
} }
/* call lmdif. */ info = ei_lmdif<Functor,Scalar>(
x, fvec,
maxfev = (n + 1) * 200; nfev,
ftol = tol; fjac, ipvt, qtf, diag,
xtol = tol; 1,
gtol = 0.; 100.,
epsfcn = 0.; (n+1)*200,
mode = 1; tol, tol, Scalar(0.), Scalar(0.)
nprint = 0; );
mp5n = m + n * 5; return (info==8)?4:info;
info = lmdif(fcn, p, m, n, &x[1], &fvec[1], ftol, xtol, gtol, maxfev, }
epsfcn, &wa[1], mode, factor, nprint, &nfev, &wa[mp5n +
1], m, &iwa[1], &wa[n + 1], &wa[(n << 1) + 1], &wa[n * 3 + 1],
&wa[(n << 2) + 1], &wa[n * 5 + 1]);
if (info == 8) {
info = 4;
}
/* L10: */
return info;
/* last card of subroutine lmdif1. */
} /* lmdif1_ */

View File

@ -734,7 +734,7 @@ void testLmdif()
const int m=15, n=3; const int m=15, n=3;
int info, nfev; int info, nfev;
double fnorm, covfac, covar_ftol; double fnorm, covfac, covar_ftol;
VectorXd x(n), fvec(m), diag(n); VectorXd x(n), fvec(m), diag(n), qtf;
MatrixXd fjac; MatrixXd fjac;
VectorXi ipvt; VectorXi ipvt;
@ -742,7 +742,7 @@ void testLmdif()
x.setConstant(n, 1.); x.setConstant(n, 1.);
// do the computation // do the computation
info = ei_lmdif<lmdif_functor, double>(x, fvec, nfev, fjac, ipvt, diag); info = ei_lmdif<lmdif_functor, double>(x, fvec, nfev, fjac, ipvt, qtf, diag);
// check return values // check return values
VERIFY( 1 == info); VERIFY( 1 == info);