From 1d53ce8d48f487f48bca8ed5e44038555b16424d Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Mon, 10 Aug 2009 16:21:22 +0200 Subject: [PATCH] wrapper for hybrj --- .../Eigen/src/NonLinear/MathFunctions.h | 41 +++++ unsupported/test/NonLinear.cpp | 144 ++++++++---------- 2 files changed, 106 insertions(+), 79 deletions(-) diff --git a/unsupported/Eigen/src/NonLinear/MathFunctions.h b/unsupported/Eigen/src/NonLinear/MathFunctions.h index 3a967ebdd..e032368cb 100644 --- a/unsupported/Eigen/src/NonLinear/MathFunctions.h +++ b/unsupported/Eigen/src/NonLinear/MathFunctions.h @@ -89,6 +89,47 @@ int ei_hybrd( ); } +template +int ei_hybrj( + Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > &x, + Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > &fvec, + int &nfev, + int &njev, + Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > &fjac, + Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > &R, + Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > &qtf, + Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > &diag, + int mode=1, + int maxfev = 1000, + Scalar factor = Scalar(100.), + Scalar xtol = Eigen::ei_sqrt(Eigen::machine_epsilon()), + int nprint=0 + ) +{ + int n = x.size(); + int lr = (n*(n+1))/2; + Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > wa1(n), wa2(n), wa3(n), wa4(n); + + fvec.resize(n); + qtf.resize(n); + R.resize(lr); + int ldfjac = n; + fjac.resize(ldfjac, n); + return hybrj ( + Functor::f, 0, + n, x.data(), fvec.data(), + fjac.data(), ldfjac, + xtol, maxfev, + diag.data(), mode, + factor, + nprint, + &nfev, + &njev, + R.data(), lr, + qtf.data(), + wa1.data(), wa2.data(), wa3.data(), wa4.data() + ); +} template int ei_lmder1( diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp index 5ec89cc61..c2899d11a 100644 --- a/unsupported/test/NonLinear.cpp +++ b/unsupported/test/NonLinear.cpp @@ -349,100 +349,86 @@ void testHybrj1() for (j=1; j<=n; j++) VERIFY_IS_APPROX(x[j-1], x_ref[j-1]); } -int fcn_hybrj(void * /*p*/, int n, const double *x, double *fvec, double *fjac, int ldfjac, - int iflag) -{ - - /* subroutine fcn for hybrj example. */ - int j, k; - double one=1, temp, temp1, temp2, three=3, two=2, zero=0, four=4; +struct hybrj_functor { + static int f(void * /*p*/, int n, const double *x, double *fvec, double *fjac, int ldfjac, + int iflag) + { - if (iflag == 0) - { - /* insert print statements here when nprint is positive. */ - return 0; - } + /* subroutine fcn for hybrj example. */ - if (iflag != 2) - { - for (k=1; k <= n; k++) - { - temp = (three - two*x[k-1])*x[k-1]; - temp1 = zero; - if (k != 1) temp1 = x[k-1-1]; - temp2 = zero; - if (k != n) temp2 = x[k+1-1]; - fvec[k-1] = temp - temp1 - two*temp2 + one; - } + int j, k; + double one=1, temp, temp1, temp2, three=3, two=2, zero=0, four=4; + + if (iflag == 0) + { + /* insert print statements here when nprint is positive. */ + return 0; + } + + if (iflag != 2) + { + for (k=1; k <= n; k++) + { + temp = (three - two*x[k-1])*x[k-1]; + temp1 = zero; + if (k != 1) temp1 = x[k-1-1]; + temp2 = zero; + if (k != n) temp2 = x[k+1-1]; + fvec[k-1] = temp - temp1 - two*temp2 + one; + } + } + else + { + for (k = 1; k <= n; k++) + { + for (j=1; j <= n; j++) + { + fjac[k-1 + ldfjac*(j-1)] = zero; + } + fjac[k-1 + ldfjac*(k-1)] = three - four*x[k-1]; + if (k != 1) fjac[k-1 + ldfjac*(k-1-1)] = -one; + if (k != n) fjac[k-1 + ldfjac*(k+1-1)] = -two; + } + } + return 0; } - else - { - for (k = 1; k <= n; k++) - { - for (j=1; j <= n; j++) - { - fjac[k-1 + ldfjac*(j-1)] = zero; - } - fjac[k-1 + ldfjac*(k-1)] = three - four*x[k-1]; - if (k != 1) fjac[k-1 + ldfjac*(k-1-1)] = -one; - if (k != n) fjac[k-1 + ldfjac*(k+1-1)] = -two; - } - } - return 0; -} +}; void testHybrj() { + const int n=9; + int info, nfev, njev, mode; + Eigen::VectorXd x(n), fvec, diag(n), R, qtf; + Eigen::MatrixXd fjac; + VectorXi ipvt; - int j, n, ldfjac, maxfev, mode, nprint, info, nfev, njev, lr; - double xtol, factor, fnorm; - double x[9], fvec[9], fjac[9*9], diag[9], r[45], qtf[9], - wa1[9], wa2[9], wa3[9], wa4[9]; + /* the following starting values provide a rough fit. */ + x.setConstant(n, -1.); - n = 9; - -/* the following starting values provide a rough solution. */ - - for (j=1; j<=9; j++) - { - x[j-1] = -1.; - } - - ldfjac = 9; - lr = 45; - -/* set xtol to the square root of the machine precision. */ -/* unless high solutions are required, */ -/* this is the recommended setting. */ - - xtol = sqrt(dpmpar(1)); - - maxfev = 1000; mode = 2; - for (j=1; j<=9; j++) - { - diag[j-1] = 1.; - } - factor = 1.e2; - nprint = 0; + diag.setConstant(n, 1.); - info = hybrj(fcn_hybrj, 0, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, - mode, factor, nprint, &nfev, &njev, r, lr, qtf, - wa1, wa2, wa3, wa4); - fnorm = enorm(n, fvec); + // do the computation + info = ei_hybrj(x,fvec, nfev, njev, fjac, R, qtf, diag, mode); - VERIFY_IS_APPROX(fnorm, 1.192636e-08); - VERIFY(nfev==11); - VERIFY(njev==1); - VERIFY(info==1); + // check return value + VERIFY( 1 == info); + VERIFY(nfev==11); + VERIFY(njev==1); - double x_ref[] = { + // check norm + VERIFY_IS_APPROX(fvec.norm(), 1.192636e-08); + + +// check x + VectorXd x_ref(n); + x_ref << -0.5706545, -0.6816283, -0.7017325, -0.7042129, -0.701369, -0.6918656, - -0.665792, -0.5960342, -0.4164121 - }; - for (j=1; j<=n; j++) VERIFY_IS_APPROX(x[j-1], x_ref[j-1]); + -0.665792, -0.5960342, -0.4164121; + VERIFY_IS_APPROX(x, x_ref); + } struct hybrd1_functor {