From b09ebe01da1381208a6007afbc8501edfb9f91d6 Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Thu, 20 Aug 2009 22:59:09 +0200 Subject: [PATCH] * porting lmdif1 to eigen * qtf was missing in lmdif signature (this is an output of the method) --- unsupported/Eigen/NonLinear | 2 +- .../Eigen/src/NonLinear/MathFunctions.h | 27 +------ unsupported/Eigen/src/NonLinear/lmdif1.h | 76 +++++++------------ unsupported/test/NonLinear.cpp | 4 +- 4 files changed, 31 insertions(+), 78 deletions(-) diff --git a/unsupported/Eigen/NonLinear b/unsupported/Eigen/NonLinear index fda61e231..7a4006aef 100644 --- a/unsupported/Eigen/NonLinear +++ b/unsupported/Eigen/NonLinear @@ -130,7 +130,6 @@ Scalar ei_enorm ( int n, const Scalar *x ){ #include "src/NonLinear/lmder.h" #include "src/NonLinear/hybrd.h" #include "src/NonLinear/lmstr.h" -#include "src/NonLinear/lmdif1.h" #include "src/NonLinear/lmdif.h" #include "src/NonLinear/hybrj.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/hybrd1.h" #include "src/NonLinear/hybrj1.h" +#include "src/NonLinear/lmdif1.h" //@} diff --git a/unsupported/Eigen/src/NonLinear/MathFunctions.h b/unsupported/Eigen/src/NonLinear/MathFunctions.h index d4efa5e01..ece132fd9 100644 --- a/unsupported/Eigen/src/NonLinear/MathFunctions.h +++ b/unsupported/Eigen/src/NonLinear/MathFunctions.h @@ -209,6 +209,7 @@ int ei_lmdif( int &nfev, Matrix< Scalar, Dynamic, Dynamic > &fjac, VectorXi &ipvt, + Matrix< Scalar, Dynamic, 1 > &qtf, Matrix< Scalar, Dynamic, 1 > &diag, int mode=1, Scalar factor = 100., @@ -221,7 +222,6 @@ int ei_lmdif( ) { Matrix< Scalar, Dynamic, 1 > - qtf(x.size()), wa1(x.size()), wa2(x.size()), wa3(x.size()), wa4(fvec.size()); int ldfjac = fvec.size(); @@ -229,6 +229,7 @@ int ei_lmdif( ipvt.resize(x.size()); fjac.resize(ldfjac, x.size()); diag.resize(x.size()); + qtf.resize(x.size()); return lmdif_template ( Functor::f, 0, fvec.size(), x.size(), x.data(), fvec.data(), @@ -246,30 +247,6 @@ int ei_lmdif( ); } -template -int ei_lmdif1( - Matrix< Scalar, Dynamic, 1 > &x, - Matrix< Scalar, Dynamic, 1 > &fvec, - Scalar tol = ei_sqrt(epsilon()) - ) -{ - 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 ( - Functor::f, 0, - fvec.size(), n, x.data(), fvec.data(), - tol, - iwa.data(), - wa.data(), lwa - ); -} - template void ei_chkder( Matrix< Scalar, Dynamic, 1 > &x, diff --git a/unsupported/Eigen/src/NonLinear/lmdif1.h b/unsupported/Eigen/src/NonLinear/lmdif1.h index 8dfd59db7..7a82c551f 100644 --- a/unsupported/Eigen/src/NonLinear/lmdif1.h +++ b/unsupported/Eigen/src/NonLinear/lmdif1.h @@ -1,56 +1,32 @@ -template -int lmdif1_template(minpack_func_mn fcn, void *p, int m, int n, Scalar *x, - Scalar *fvec, Scalar tol, int *iwa, - Scalar *wa, int lwa) +template +int ei_lmdif1( + Matrix< Scalar, Dynamic, 1 > &x, + Matrix< Scalar, Dynamic, 1 > &fvec, + Scalar tol = ei_sqrt(epsilon()) + ) { - /* 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.; - - int mp5n, mode, nfev; - Scalar ftol, gtol, xtol; - 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; + /* check the input parameters for errors. */ + if (n <= 0 || m < n || tol < 0.) { + printf("ei_lmder1 bad args : m,n,tol,..."); + return 0; } -/* call lmdif. */ - - maxfev = (n + 1) * 200; - ftol = tol; - xtol = tol; - gtol = 0.; - epsfcn = 0.; - mode = 1; - nprint = 0; - mp5n = m + n * 5; - 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_ */ + info = ei_lmdif( + x, fvec, + nfev, + fjac, ipvt, qtf, diag, + 1, + 100., + (n+1)*200, + tol, tol, Scalar(0.), Scalar(0.) + ); + return (info==8)?4:info; +} diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp index 7591634ff..26d1bdf6b 100644 --- a/unsupported/test/NonLinear.cpp +++ b/unsupported/test/NonLinear.cpp @@ -734,7 +734,7 @@ void testLmdif() const int m=15, n=3; int info, nfev; double fnorm, covfac, covar_ftol; - VectorXd x(n), fvec(m), diag(n); + VectorXd x(n), fvec(m), diag(n), qtf; MatrixXd fjac; VectorXi ipvt; @@ -742,7 +742,7 @@ void testLmdif() x.setConstant(n, 1.); // do the computation - info = ei_lmdif(x, fvec, nfev, fjac, ipvt, diag); + info = ei_lmdif(x, fvec, nfev, fjac, ipvt, qtf, diag); // check return values VERIFY( 1 == info);