From acd757737a9a2e728ea701fb3b55d525dba3fc0c Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Sun, 23 Aug 2009 04:39:22 +0200 Subject: [PATCH] beautify Functor for lmder : we now have f,df,debug methods --- unsupported/Eigen/src/NonLinear/lmder.h | 10 +- unsupported/test/NonLinear.cpp | 591 ++++++++++++------------ 2 files changed, 310 insertions(+), 291 deletions(-) diff --git a/unsupported/Eigen/src/NonLinear/lmder.h b/unsupported/Eigen/src/NonLinear/lmder.h index 71b8c3ed2..a4422258e 100644 --- a/unsupported/Eigen/src/NonLinear/lmder.h +++ b/unsupported/Eigen/src/NonLinear/lmder.h @@ -64,7 +64,7 @@ L20: /* evaluate the function at the starting point */ /* and calculate its norm. */ - iflag = Functor::f(x, fvec, fjac, 1); + iflag = Functor::f(x, fvec); nfev = 1; if (iflag < 0) { goto L300; @@ -82,7 +82,7 @@ L30: /* calculate the jacobian matrix. */ - iflag = Functor::f(x, fvec, fjac, 2); + iflag = Functor::df(x, fjac); ++njev; if (iflag < 0) { goto L300; @@ -95,7 +95,7 @@ L30: } iflag = 0; if ((iter - 1) % nprint == 0) { - iflag = Functor::f(x, fvec, fjac, 0); + iflag = Functor::debug(x, fvec, fjac); } if (iflag < 0) { goto L300; @@ -237,7 +237,7 @@ L200: /* evaluate the function at x + p and calculate its norm. */ - iflag = Functor::f(wa2, wa4, fjac, 1); + iflag = Functor::f(wa2, wa4); ++nfev; if (iflag < 0) { goto L300; @@ -378,7 +378,7 @@ L300: } iflag = 0; if (nprint > 0) { - iflag = Functor::f(x, fvec, fjac, 0); + iflag = Functor::debug(x, fvec, fjac); } return info; diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp index 73d07b0bb..35dddba8b 100644 --- a/unsupported/test/NonLinear.cpp +++ b/unsupported/test/NonLinear.cpp @@ -102,38 +102,35 @@ void testChkder() struct lmder_functor { - static int f(const VectorXd &x, VectorXd &fvec, MatrixXd &fjac, int iflag) + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &x, VectorXd &fvec) { - - /* subroutine fcn for lmder1 example. */ - - int i; - double tmp1, tmp2, tmp3, tmp4; + double tmp1, tmp2, tmp3; double y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; - if (iflag != 2) + for (int i = 0; i < 15; i++) { - for (i = 0; i < 15; i++) - { - tmp1 = i+1; - tmp2 = 16 - i - 1; - tmp3 = (i>=8)? tmp2 : tmp1; - fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); - } + tmp1 = i+1; + tmp2 = 16 - i - 1; + tmp3 = (i>=8)? tmp2 : tmp1; + fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); } - else + return 0; + } + + static int df(const VectorXd &x, MatrixXd &fjac) + { + double tmp1, tmp2, tmp3, tmp4; + for (int i = 0; i < 15; i++) { - for ( i = 0; i < 15; i++) - { - tmp1 = i+1; - tmp2 = 16 - i - 1; - tmp3 = (i>=8)? tmp2 : tmp1; - tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; - fjac(i,0) = -1; - fjac(i,1) = tmp1*tmp2/tmp4; - fjac(i,2) = tmp1*tmp3/tmp4; - } + tmp1 = i+1; + tmp2 = 16 - i - 1; + tmp3 = (i>=8)? tmp2 : tmp1; + tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; + fjac(i,0) = -1; + fjac(i,1) = tmp1*tmp2/tmp4; + fjac(i,2) = tmp1*tmp3/tmp4; } return 0; } @@ -598,35 +595,39 @@ void testLmdif() } struct chwirut2_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double m_x[54]; + static const double m_y[54]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double _x[54] = { 0.500E0, 1.000E0, 1.750E0, 3.750E0, 5.750E0, 0.875E0, 2.250E0, 3.250E0, 5.250E0, 0.750E0, 1.750E0, 2.750E0, 4.750E0, 0.625E0, 1.250E0, 2.250E0, 4.250E0, .500E0, 3.000E0, .750E0, 3.000E0, 1.500E0, 6.000E0, 3.000E0, 6.000E0, 1.500E0, 3.000E0, .500E0, 2.000E0, 4.000E0, .750E0, 2.000E0, 5.000E0, .750E0, 2.250E0, 3.750E0, 5.750E0, 3.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .500E0, 6.000E0, 3.000E0, .500E0, 2.750E0, .500E0, 1.750E0}; - static const double y[54] = { 92.9000E0 ,57.1000E0 ,31.0500E0 ,11.5875E0 ,8.0250E0 ,63.6000E0 ,21.4000E0 ,14.2500E0 ,8.4750E0 ,63.8000E0 ,26.8000E0 ,16.4625E0 ,7.1250E0 ,67.3000E0 ,41.0000E0 ,21.1500E0 ,8.1750E0 ,81.5000E0 ,13.1200E0 ,59.9000E0 ,14.6200E0 ,32.9000E0 ,5.4400E0 ,12.5600E0 ,5.4400E0 ,32.0000E0 ,13.9500E0 ,75.8000E0 ,20.0000E0 ,10.4200E0 ,59.5000E0 ,21.6700E0 ,8.5500E0 ,62.0000E0 ,20.2000E0 ,7.7600E0 ,3.7500E0 ,11.8100E0 ,54.7000E0 ,23.7000E0 ,11.5500E0 ,61.3000E0 ,17.7000E0 ,8.7400E0 ,59.2000E0 ,16.3000E0 ,8.6200E0 ,81.0000E0 ,4.8700E0 ,14.6200E0 ,81.7000E0 ,17.1700E0 ,81.3000E0 ,28.9000E0 }; int i; assert(b.size()==3); assert(fvec.size()==54); + for(i=0; i<54; i++) { + double x = m_x[i]; + fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i]; + } + return 0; + } + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==3); assert(fjac.rows()==54); assert(fjac.cols()==3); - if (iflag != 2) {// compute fvec at b - for(i=0; i<54; i++) { - double x = _x[i]; - fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - y[i]; - } - } - else { // compute fjac at b - for(i=0; i<54; i++) { - double x = _x[i]; - double factor = 1./(b[1]+b[2]*x); - double e = exp(-b[0]*x); - fjac(i,0) = -x*e*factor; - fjac(i,1) = -e*factor*factor; - fjac(i,2) = -x*e*factor*factor; - } + for(int i=0; i<54; i++) { + double x = m_x[i]; + double factor = 1./(b[1]+b[2]*x); + double e = exp(-b[0]*x); + fjac(i,0) = -x*e*factor; + fjac(i,1) = -e*factor*factor; + fjac(i,2) = -x*e*factor*factor; } return 0; } }; +const double chwirut2_functor::m_x[54] = { 0.500E0, 1.000E0, 1.750E0, 3.750E0, 5.750E0, 0.875E0, 2.250E0, 3.250E0, 5.250E0, 0.750E0, 1.750E0, 2.750E0, 4.750E0, 0.625E0, 1.250E0, 2.250E0, 4.250E0, .500E0, 3.000E0, .750E0, 3.000E0, 1.500E0, 6.000E0, 3.000E0, 6.000E0, 1.500E0, 3.000E0, .500E0, 2.000E0, 4.000E0, .750E0, 2.000E0, 5.000E0, .750E0, 2.250E0, 3.750E0, 5.750E0, 3.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .500E0, 6.000E0, 3.000E0, .500E0, 2.750E0, .500E0, 1.750E0}; +const double chwirut2_functor::m_y[54] = { 92.9000E0 ,57.1000E0 ,31.0500E0 ,11.5875E0 ,8.0250E0 ,63.6000E0 ,21.4000E0 ,14.2500E0 ,8.4750E0 ,63.8000E0 ,26.8000E0 ,16.4625E0 ,7.1250E0 ,67.3000E0 ,41.0000E0 ,21.1500E0 ,8.1750E0 ,81.5000E0 ,13.1200E0 ,59.9000E0 ,14.6200E0 ,32.9000E0 ,5.4400E0 ,12.5600E0 ,5.4400E0 ,32.0000E0 ,13.9500E0 ,75.8000E0 ,20.0000E0 ,10.4200E0 ,59.5000E0 ,21.6700E0 ,8.5500E0 ,62.0000E0 ,20.2000E0 ,7.7600E0 ,3.7500E0 ,11.8100E0 ,54.7000E0 ,23.7000E0 ,11.5500E0 ,61.3000E0 ,17.7000E0 ,8.7400E0 ,59.2000E0 ,16.3000E0 ,8.6200E0 ,81.0000E0 ,4.8700E0 ,14.6200E0 ,81.7000E0 ,17.1700E0 ,81.3000E0 ,28.9000E0 }; // http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml void testNistChwirut2(void) @@ -678,30 +679,32 @@ void testNistChwirut2(void) struct misra1a_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double m_x[14]; + static const double m_y[14]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0}; - static const double y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0}; - int i; - assert(b.size()==2); assert(fvec.size()==14); + for(int i=0; i<14; i++) { + fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ; + } + return 0; + } + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==2); assert(fjac.rows()==14); assert(fjac.cols()==2); - if (iflag != 2) {// compute fvec at b - for(i=0; i<14; i++) { - fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i] ; - } - } - else { // compute fjac at b - for(i=0; i<14; i++) { - fjac(i,0) = (1.-exp(-b[1]*x[i])); - fjac(i,1) = (b[0]*x[i]*exp(-b[1]*x[i])); - } + for(int i=0; i<14; i++) { + fjac(i,0) = (1.-exp(-b[1]*m_x[i])); + fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i])); } return 0; } }; +const double misra1a_functor::m_x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0}; +const double misra1a_functor::m_y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0}; // http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml void testNistMisra1a(void) @@ -749,41 +752,44 @@ void testNistMisra1a(void) } struct hahn1_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double m_x[236]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double _x[236] = { 24.41E0 , 34.82E0 , 44.09E0 , 45.07E0 , 54.98E0 , 65.51E0 , 70.53E0 , 75.70E0 , 89.57E0 , 91.14E0 , 96.40E0 , 97.19E0 , 114.26E0 , 120.25E0 , 127.08E0 , 133.55E0 , 133.61E0 , 158.67E0 , 172.74E0 , 171.31E0 , 202.14E0 , 220.55E0 , 221.05E0 , 221.39E0 , 250.99E0 , 268.99E0 , 271.80E0 , 271.97E0 , 321.31E0 , 321.69E0 , 330.14E0 , 333.03E0 , 333.47E0 , 340.77E0 , 345.65E0 , 373.11E0 , 373.79E0 , 411.82E0 , 419.51E0 , 421.59E0 , 422.02E0 , 422.47E0 , 422.61E0 , 441.75E0 , 447.41E0 , 448.7E0 , 472.89E0 , 476.69E0 , 522.47E0 , 522.62E0 , 524.43E0 , 546.75E0 , 549.53E0 , 575.29E0 , 576.00E0 , 625.55E0 , 20.15E0 , 28.78E0 , 29.57E0 , 37.41E0 , 39.12E0 , 50.24E0 , 61.38E0 , 66.25E0 , 73.42E0 , 95.52E0 , 107.32E0 , 122.04E0 , 134.03E0 , 163.19E0 , 163.48E0 , 175.70E0 , 179.86E0 , 211.27E0 , 217.78E0 , 219.14E0 , 262.52E0 , 268.01E0 , 268.62E0 , 336.25E0 , 337.23E0 , 339.33E0 , 427.38E0 , 428.58E0 , 432.68E0 , 528.99E0 , 531.08E0 , 628.34E0 , 253.24E0 , 273.13E0 , 273.66E0 , 282.10E0 , 346.62E0 , 347.19E0 , 348.78E0 , 351.18E0 , 450.10E0 , 450.35E0 , 451.92E0 , 455.56E0 , 552.22E0 , 553.56E0 , 555.74E0 , 652.59E0 , 656.20E0 , 14.13E0 , 20.41E0 , 31.30E0 , 33.84E0 , 39.70E0 , 48.83E0 , 54.50E0 , 60.41E0 , 72.77E0 , 75.25E0 , 86.84E0 , 94.88E0 , 96.40E0 , 117.37E0 , 139.08E0 , 147.73E0 , 158.63E0 , 161.84E0 , 192.11E0 , 206.76E0 , 209.07E0 , 213.32E0 , 226.44E0 , 237.12E0 , 330.90E0 , 358.72E0 , 370.77E0 , 372.72E0 , 396.24E0 , 416.59E0 , 484.02E0 , 495.47E0 , 514.78E0 , 515.65E0 , 519.47E0 , 544.47E0 , 560.11E0 , 620.77E0 , 18.97E0 , 28.93E0 , 33.91E0 , 40.03E0 , 44.66E0 , 49.87E0 , 55.16E0 , 60.90E0 , 72.08E0 , 85.15E0 , 97.06E0 , 119.63E0 , 133.27E0 , 143.84E0 , 161.91E0 , 180.67E0 , 198.44E0 , 226.86E0 , 229.65E0 , 258.27E0 , 273.77E0 , 339.15E0 , 350.13E0 , 362.75E0 , 371.03E0 , 393.32E0 , 448.53E0 , 473.78E0 , 511.12E0 , 524.70E0 , 548.75E0 , 551.64E0 , 574.02E0 , 623.86E0 , 21.46E0 , 24.33E0 , 33.43E0 , 39.22E0 , 44.18E0 , 55.02E0 , 94.33E0 , 96.44E0 , 118.82E0 , 128.48E0 , 141.94E0 , 156.92E0 , 171.65E0 , 190.00E0 , 223.26E0 , 223.88E0 , 231.50E0 , 265.05E0 , 269.44E0 , 271.78E0 , 273.46E0 , 334.61E0 , 339.79E0 , 349.52E0 , 358.18E0 , 377.98E0 , 394.77E0 , 429.66E0 , 468.22E0 , 487.27E0 , 519.54E0 , 523.03E0 , 612.99E0 , 638.59E0 , 641.36E0 , 622.05E0 , 631.50E0 , 663.97E0 , 646.9E0 , 748.29E0 , 749.21E0 , 750.14E0 , 647.04E0 , 646.89E0 , 746.9E0 , 748.43E0 , 747.35E0 , 749.27E0 , 647.61E0 , 747.78E0 , 750.51E0 , 851.37E0 , 845.97E0 , 847.54E0 , 849.93E0 , 851.61E0 , 849.75E0 , 850.98E0 , 848.23E0}; - static const double _y[236] = { .591E0 , 1.547E0 , 2.902E0 , 2.894E0 , 4.703E0 , 6.307E0 , 7.03E0 , 7.898E0 , 9.470E0 , 9.484E0 , 10.072E0 , 10.163E0 , 11.615E0 , 12.005E0 , 12.478E0 , 12.982E0 , 12.970E0 , 13.926E0 , 14.452E0 , 14.404E0 , 15.190E0 , 15.550E0 , 15.528E0 , 15.499E0 , 16.131E0 , 16.438E0 , 16.387E0 , 16.549E0 , 16.872E0 , 16.830E0 , 16.926E0 , 16.907E0 , 16.966E0 , 17.060E0 , 17.122E0 , 17.311E0 , 17.355E0 , 17.668E0 , 17.767E0 , 17.803E0 , 17.765E0 , 17.768E0 , 17.736E0 , 17.858E0 , 17.877E0 , 17.912E0 , 18.046E0 , 18.085E0 , 18.291E0 , 18.357E0 , 18.426E0 , 18.584E0 , 18.610E0 , 18.870E0 , 18.795E0 , 19.111E0 , .367E0 , .796E0 , 0.892E0 , 1.903E0 , 2.150E0 , 3.697E0 , 5.870E0 , 6.421E0 , 7.422E0 , 9.944E0 , 11.023E0 , 11.87E0 , 12.786E0 , 14.067E0 , 13.974E0 , 14.462E0 , 14.464E0 , 15.381E0 , 15.483E0 , 15.59E0 , 16.075E0 , 16.347E0 , 16.181E0 , 16.915E0 , 17.003E0 , 16.978E0 , 17.756E0 , 17.808E0 , 17.868E0 , 18.481E0 , 18.486E0 , 19.090E0 , 16.062E0 , 16.337E0 , 16.345E0 , 16.388E0 , 17.159E0 , 17.116E0 , 17.164E0 , 17.123E0 , 17.979E0 , 17.974E0 , 18.007E0 , 17.993E0 , 18.523E0 , 18.669E0 , 18.617E0 , 19.371E0 , 19.330E0 , 0.080E0 , 0.248E0 , 1.089E0 , 1.418E0 , 2.278E0 , 3.624E0 , 4.574E0 , 5.556E0 , 7.267E0 , 7.695E0 , 9.136E0 , 9.959E0 , 9.957E0 , 11.600E0 , 13.138E0 , 13.564E0 , 13.871E0 , 13.994E0 , 14.947E0 , 15.473E0 , 15.379E0 , 15.455E0 , 15.908E0 , 16.114E0 , 17.071E0 , 17.135E0 , 17.282E0 , 17.368E0 , 17.483E0 , 17.764E0 , 18.185E0 , 18.271E0 , 18.236E0 , 18.237E0 , 18.523E0 , 18.627E0 , 18.665E0 , 19.086E0 , 0.214E0 , 0.943E0 , 1.429E0 , 2.241E0 , 2.951E0 , 3.782E0 , 4.757E0 , 5.602E0 , 7.169E0 , 8.920E0 , 10.055E0 , 12.035E0 , 12.861E0 , 13.436E0 , 14.167E0 , 14.755E0 , 15.168E0 , 15.651E0 , 15.746E0 , 16.216E0 , 16.445E0 , 16.965E0 , 17.121E0 , 17.206E0 , 17.250E0 , 17.339E0 , 17.793E0 , 18.123E0 , 18.49E0 , 18.566E0 , 18.645E0 , 18.706E0 , 18.924E0 , 19.1E0 , 0.375E0 , 0.471E0 , 1.504E0 , 2.204E0 , 2.813E0 , 4.765E0 , 9.835E0 , 10.040E0 , 11.946E0 , 12.596E0 , 13.303E0 , 13.922E0 , 14.440E0 , 14.951E0 , 15.627E0 , 15.639E0 , 15.814E0 , 16.315E0 , 16.334E0 , 16.430E0 , 16.423E0 , 17.024E0 , 17.009E0 , 17.165E0 , 17.134E0 , 17.349E0 , 17.576E0 , 17.848E0 , 18.090E0 , 18.276E0 , 18.404E0 , 18.519E0 , 19.133E0 , 19.074E0 , 19.239E0 , 19.280E0 , 19.101E0 , 19.398E0 , 19.252E0 , 19.89E0 , 20.007E0 , 19.929E0 , 19.268E0 , 19.324E0 , 20.049E0 , 20.107E0 , 20.062E0 , 20.065E0 , 19.286E0 , 19.972E0 , 20.088E0 , 20.743E0 , 20.83E0 , 20.935E0 , 21.035E0 , 20.93E0 , 21.074E0 , 21.085E0 , 20.935E0 }; - int i; + static const double m_y[236] = { .591E0 , 1.547E0 , 2.902E0 , 2.894E0 , 4.703E0 , 6.307E0 , 7.03E0 , 7.898E0 , 9.470E0 , 9.484E0 , 10.072E0 , 10.163E0 , 11.615E0 , 12.005E0 , 12.478E0 , 12.982E0 , 12.970E0 , 13.926E0 , 14.452E0 , 14.404E0 , 15.190E0 , 15.550E0 , 15.528E0 , 15.499E0 , 16.131E0 , 16.438E0 , 16.387E0 , 16.549E0 , 16.872E0 , 16.830E0 , 16.926E0 , 16.907E0 , 16.966E0 , 17.060E0 , 17.122E0 , 17.311E0 , 17.355E0 , 17.668E0 , 17.767E0 , 17.803E0 , 17.765E0 , 17.768E0 , 17.736E0 , 17.858E0 , 17.877E0 , 17.912E0 , 18.046E0 , 18.085E0 , 18.291E0 , 18.357E0 , 18.426E0 , 18.584E0 , 18.610E0 , 18.870E0 , 18.795E0 , 19.111E0 , .367E0 , .796E0 , 0.892E0 , 1.903E0 , 2.150E0 , 3.697E0 , 5.870E0 , 6.421E0 , 7.422E0 , 9.944E0 , 11.023E0 , 11.87E0 , 12.786E0 , 14.067E0 , 13.974E0 , 14.462E0 , 14.464E0 , 15.381E0 , 15.483E0 , 15.59E0 , 16.075E0 , 16.347E0 , 16.181E0 , 16.915E0 , 17.003E0 , 16.978E0 , 17.756E0 , 17.808E0 , 17.868E0 , 18.481E0 , 18.486E0 , 19.090E0 , 16.062E0 , 16.337E0 , 16.345E0 , 16.388E0 , 17.159E0 , 17.116E0 , 17.164E0 , 17.123E0 , 17.979E0 , 17.974E0 , 18.007E0 , 17.993E0 , 18.523E0 , 18.669E0 , 18.617E0 , 19.371E0 , 19.330E0 , 0.080E0 , 0.248E0 , 1.089E0 , 1.418E0 , 2.278E0 , 3.624E0 , 4.574E0 , 5.556E0 , 7.267E0 , 7.695E0 , 9.136E0 , 9.959E0 , 9.957E0 , 11.600E0 , 13.138E0 , 13.564E0 , 13.871E0 , 13.994E0 , 14.947E0 , 15.473E0 , 15.379E0 , 15.455E0 , 15.908E0 , 16.114E0 , 17.071E0 , 17.135E0 , 17.282E0 , 17.368E0 , 17.483E0 , 17.764E0 , 18.185E0 , 18.271E0 , 18.236E0 , 18.237E0 , 18.523E0 , 18.627E0 , 18.665E0 , 19.086E0 , 0.214E0 , 0.943E0 , 1.429E0 , 2.241E0 , 2.951E0 , 3.782E0 , 4.757E0 , 5.602E0 , 7.169E0 , 8.920E0 , 10.055E0 , 12.035E0 , 12.861E0 , 13.436E0 , 14.167E0 , 14.755E0 , 15.168E0 , 15.651E0 , 15.746E0 , 16.216E0 , 16.445E0 , 16.965E0 , 17.121E0 , 17.206E0 , 17.250E0 , 17.339E0 , 17.793E0 , 18.123E0 , 18.49E0 , 18.566E0 , 18.645E0 , 18.706E0 , 18.924E0 , 19.1E0 , 0.375E0 , 0.471E0 , 1.504E0 , 2.204E0 , 2.813E0 , 4.765E0 , 9.835E0 , 10.040E0 , 11.946E0 , 12.596E0 , 13.303E0 , 13.922E0 , 14.440E0 , 14.951E0 , 15.627E0 , 15.639E0 , 15.814E0 , 16.315E0 , 16.334E0 , 16.430E0 , 16.423E0 , 17.024E0 , 17.009E0 , 17.165E0 , 17.134E0 , 17.349E0 , 17.576E0 , 17.848E0 , 18.090E0 , 18.276E0 , 18.404E0 , 18.519E0 , 19.133E0 , 19.074E0 , 19.239E0 , 19.280E0 , 19.101E0 , 19.398E0 , 19.252E0 , 19.89E0 , 20.007E0 , 19.929E0 , 19.268E0 , 19.324E0 , 20.049E0 , 20.107E0 , 20.062E0 , 20.065E0 , 19.286E0 , 19.972E0 , 20.088E0 , 20.743E0 , 20.83E0 , 20.935E0 , 21.035E0 , 20.93E0 , 21.074E0 , 21.085E0 , 20.935E0 }; -// static int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; + // static int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; assert(b.size()==7); assert(fvec.size()==236); + for(int i=0; i<236; i++) { + double x=m_x[i], xx=x*x, xxx=xx*x; + fvec[i] = (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) / (1.+b[4]*x+b[5]*xx+b[6]*xxx) - m_y[i]; + } + return 0; + } + + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==7); assert(fjac.rows()==236); assert(fjac.cols()==7); - if (iflag != 2) {// compute fvec at x - for(i=0; i<236; i++) { - double x=_x[i], xx=x*x, xxx=xx*x; - fvec[i] = (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) / (1.+b[4]*x+b[5]*xx+b[6]*xxx) - _y[i]; - } - } - else { // compute fjac at x - for(i=0; i<236; i++) { - double x=_x[i], xx=x*x, xxx=xx*x; - double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); - fjac(i,0) = 1.*fact; - fjac(i,1) = x*fact; - fjac(i,2) = xx*fact; - fjac(i,3) = xxx*fact; - fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; - fjac(i,4) = x*fact; - fjac(i,5) = xx*fact; - fjac(i,6) = xxx*fact; - } + for(int i=0; i<236; i++) { + double x=m_x[i], xx=x*x, xxx=xx*x; + double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); + fjac(i,0) = 1.*fact; + fjac(i,1) = x*fact; + fjac(i,2) = xx*fact; + fjac(i,3) = xxx*fact; + fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; + fjac(i,4) = x*fact; + fjac(i,5) = xx*fact; + fjac(i,6) = xxx*fact; } return 0; } }; +const double hahn1_functor::m_x[236] = { 24.41E0 , 34.82E0 , 44.09E0 , 45.07E0 , 54.98E0 , 65.51E0 , 70.53E0 , 75.70E0 , 89.57E0 , 91.14E0 , 96.40E0 , 97.19E0 , 114.26E0 , 120.25E0 , 127.08E0 , 133.55E0 , 133.61E0 , 158.67E0 , 172.74E0 , 171.31E0 , 202.14E0 , 220.55E0 , 221.05E0 , 221.39E0 , 250.99E0 , 268.99E0 , 271.80E0 , 271.97E0 , 321.31E0 , 321.69E0 , 330.14E0 , 333.03E0 , 333.47E0 , 340.77E0 , 345.65E0 , 373.11E0 , 373.79E0 , 411.82E0 , 419.51E0 , 421.59E0 , 422.02E0 , 422.47E0 , 422.61E0 , 441.75E0 , 447.41E0 , 448.7E0 , 472.89E0 , 476.69E0 , 522.47E0 , 522.62E0 , 524.43E0 , 546.75E0 , 549.53E0 , 575.29E0 , 576.00E0 , 625.55E0 , 20.15E0 , 28.78E0 , 29.57E0 , 37.41E0 , 39.12E0 , 50.24E0 , 61.38E0 , 66.25E0 , 73.42E0 , 95.52E0 , 107.32E0 , 122.04E0 , 134.03E0 , 163.19E0 , 163.48E0 , 175.70E0 , 179.86E0 , 211.27E0 , 217.78E0 , 219.14E0 , 262.52E0 , 268.01E0 , 268.62E0 , 336.25E0 , 337.23E0 , 339.33E0 , 427.38E0 , 428.58E0 , 432.68E0 , 528.99E0 , 531.08E0 , 628.34E0 , 253.24E0 , 273.13E0 , 273.66E0 , 282.10E0 , 346.62E0 , 347.19E0 , 348.78E0 , 351.18E0 , 450.10E0 , 450.35E0 , 451.92E0 , 455.56E0 , 552.22E0 , 553.56E0 , 555.74E0 , 652.59E0 , 656.20E0 , 14.13E0 , 20.41E0 , 31.30E0 , 33.84E0 , 39.70E0 , 48.83E0 , 54.50E0 , 60.41E0 , 72.77E0 , 75.25E0 , 86.84E0 , 94.88E0 , 96.40E0 , 117.37E0 , 139.08E0 , 147.73E0 , 158.63E0 , 161.84E0 , 192.11E0 , 206.76E0 , 209.07E0 , 213.32E0 , 226.44E0 , 237.12E0 , 330.90E0 , 358.72E0 , 370.77E0 , 372.72E0 , 396.24E0 , 416.59E0 , 484.02E0 , 495.47E0 , 514.78E0 , 515.65E0 , 519.47E0 , 544.47E0 , 560.11E0 , 620.77E0 , 18.97E0 , 28.93E0 , 33.91E0 , 40.03E0 , 44.66E0 , 49.87E0 , 55.16E0 , 60.90E0 , 72.08E0 , 85.15E0 , 97.06E0 , 119.63E0 , 133.27E0 , 143.84E0 , 161.91E0 , 180.67E0 , 198.44E0 , 226.86E0 , 229.65E0 , 258.27E0 , 273.77E0 , 339.15E0 , 350.13E0 , 362.75E0 , 371.03E0 , 393.32E0 , 448.53E0 , 473.78E0 , 511.12E0 , 524.70E0 , 548.75E0 , 551.64E0 , 574.02E0 , 623.86E0 , 21.46E0 , 24.33E0 , 33.43E0 , 39.22E0 , 44.18E0 , 55.02E0 , 94.33E0 , 96.44E0 , 118.82E0 , 128.48E0 , 141.94E0 , 156.92E0 , 171.65E0 , 190.00E0 , 223.26E0 , 223.88E0 , 231.50E0 , 265.05E0 , 269.44E0 , 271.78E0 , 273.46E0 , 334.61E0 , 339.79E0 , 349.52E0 , 358.18E0 , 377.98E0 , 394.77E0 , 429.66E0 , 468.22E0 , 487.27E0 , 519.54E0 , 523.03E0 , 612.99E0 , 638.59E0 , 641.36E0 , 622.05E0 , 631.50E0 , 663.97E0 , 646.9E0 , 748.29E0 , 749.21E0 , 750.14E0 , 647.04E0 , 646.89E0 , 746.9E0 , 748.43E0 , 747.35E0 , 749.27E0 , 647.61E0 , 747.78E0 , 750.51E0 , 851.37E0 , 845.97E0 , 847.54E0 , 849.93E0 , 851.61E0 , 849.75E0 , 850.98E0 , 848.23E0}; // http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml void testNistHahn1(void) @@ -842,31 +848,33 @@ void testNistHahn1(void) } struct misra1d_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double x[14]; + static const double y[14]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0}; - static const double y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0}; - int i; - assert(b.size()==2); assert(fvec.size()==14); + for(int i=0; i<14; i++) { + fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i]; + } + return 0; + } + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==2); assert(fjac.rows()==14); assert(fjac.cols()==2); - if (iflag != 2) {// compute fvec at b - for(i=0; i<14; i++) { - fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i]; - } - } - else { // compute fjac at b - for(i=0; i<14; i++) { - double den = 1.+b[1]*x[i]; - fjac(i,0) = b[1]*x[i] / den; - fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den; - } + for(int i=0; i<14; i++) { + double den = 1.+b[1]*x[i]; + fjac(i,0) = b[1]*x[i] / den; + fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den; } return 0; } }; +const double misra1d_functor::x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0}; +const double misra1d_functor::y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0}; // http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml void testNistMisra1d(void) @@ -915,34 +923,35 @@ void testNistMisra1d(void) struct lanczos1_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double x[24]; + static const double y[24]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double x[24] = { 0.000000000000E+00, 5.000000000000E-02, 1.000000000000E-01, 1.500000000000E-01, 2.000000000000E-01, 2.500000000000E-01, 3.000000000000E-01, 3.500000000000E-01, 4.000000000000E-01, 4.500000000000E-01, 5.000000000000E-01, 5.500000000000E-01, 6.000000000000E-01, 6.500000000000E-01, 7.000000000000E-01, 7.500000000000E-01, 8.000000000000E-01, 8.500000000000E-01, 9.000000000000E-01, 9.500000000000E-01, 1.000000000000E+00, 1.050000000000E+00, 1.100000000000E+00, 1.150000000000E+00 }; - static const double y[24] = { 2.513400000000E+00 ,2.044333373291E+00 ,1.668404436564E+00 ,1.366418021208E+00 ,1.123232487372E+00 ,9.268897180037E-01 ,7.679338563728E-01 ,6.388775523106E-01 ,5.337835317402E-01 ,4.479363617347E-01 ,3.775847884350E-01 ,3.197393199326E-01 ,2.720130773746E-01 ,2.324965529032E-01 ,1.996589546065E-01 ,1.722704126914E-01 ,1.493405660168E-01 ,1.300700206922E-01 ,1.138119324644E-01 ,1.000415587559E-01 ,8.833209084540E-02 ,7.833544019350E-02 ,6.976693743449E-02 ,6.239312536719E-02 }; - int i; - assert(b.size()==6); assert(fvec.size()==24); + for(int i=0; i<24; i++) + fvec[i] = b[0]*exp(-b[1]*x[i]) + b[2]*exp(-b[3]*x[i]) + b[4]*exp(-b[5]*x[i]) - y[i]; + return 0; + } + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==6); assert(fjac.rows()==24); assert(fjac.cols()==6); - if (iflag != 2) {// compute fvec at b - for(i=0; i<24; i++) { - fvec[i] = b[0]*exp(-b[1]*x[i]) + b[2]*exp(-b[3]*x[i]) + b[4]*exp(-b[5]*x[i]) - y[i]; - } - } - else { // compute fjac at b - for(i=0; i<24; i++) { - fjac(i,0) = exp(-b[1]*x[i]); - fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]); - fjac(i,2) = exp(-b[3]*x[i]); - fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]); - fjac(i,4) = exp(-b[5]*x[i]); - fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]); - } + for(int i=0; i<24; i++) { + fjac(i,0) = exp(-b[1]*x[i]); + fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]); + fjac(i,2) = exp(-b[3]*x[i]); + fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]); + fjac(i,4) = exp(-b[5]*x[i]); + fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]); } return 0; } }; +const double lanczos1_functor::x[24] = { 0.000000000000E+00, 5.000000000000E-02, 1.000000000000E-01, 1.500000000000E-01, 2.000000000000E-01, 2.500000000000E-01, 3.000000000000E-01, 3.500000000000E-01, 4.000000000000E-01, 4.500000000000E-01, 5.000000000000E-01, 5.500000000000E-01, 6.000000000000E-01, 6.500000000000E-01, 7.000000000000E-01, 7.500000000000E-01, 8.000000000000E-01, 8.500000000000E-01, 9.000000000000E-01, 9.500000000000E-01, 1.000000000000E+00, 1.050000000000E+00, 1.100000000000E+00, 1.150000000000E+00 }; +const double lanczos1_functor::y[24] = { 2.513400000000E+00 ,2.044333373291E+00 ,1.668404436564E+00 ,1.366418021208E+00 ,1.123232487372E+00 ,9.268897180037E-01 ,7.679338563728E-01 ,6.388775523106E-01 ,5.337835317402E-01 ,4.479363617347E-01 ,3.775847884350E-01 ,3.197393199326E-01 ,2.720130773746E-01 ,2.324965529032E-01 ,1.996589546065E-01 ,1.722704126914E-01 ,1.493405660168E-01 ,1.300700206922E-01 ,1.138119324644E-01 ,1.000415587559E-01 ,8.833209084540E-02 ,7.833544019350E-02 ,6.976693743449E-02 ,6.239312536719E-02 }; // http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml void testNistLanczos1(void) @@ -999,32 +1008,35 @@ void testNistLanczos1(void) } struct rat42_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double x[9]; + static const double y[9]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 }; - static const double y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 }; - int i; - assert(b.size()==3); assert(fvec.size()==9); + for(int i=0; i<9; i++) { + fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i]; + } + return 0; + } + + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==3); assert(fjac.rows()==9); assert(fjac.cols()==3); - if (iflag != 2) {// compute fvec at b - for(i=0; i<9; i++) { - fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i]; - } - } - else { // compute fjac at b - for(i=0; i<9; i++) { - double e = exp(b[1]-b[2]*x[i]); - fjac(i,0) = 1./(1.+e); - fjac(i,1) = -b[0]*e/(1.+e)/(1.+e); - fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e); - } + for(int i=0; i<9; i++) { + double e = exp(b[1]-b[2]*x[i]); + fjac(i,0) = 1./(1.+e); + fjac(i,1) = -b[0]*e/(1.+e)/(1.+e); + fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e); } return 0; } }; +const double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 }; +const double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 }; // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml void testNistRat42(void) @@ -1074,33 +1086,34 @@ void testNistRat42(void) } struct MGH10_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double x[16]; + static const double y[16]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double x[16] = { 5.000000E+01, 5.500000E+01, 6.000000E+01, 6.500000E+01, 7.000000E+01, 7.500000E+01, 8.000000E+01, 8.500000E+01, 9.000000E+01, 9.500000E+01, 1.000000E+02, 1.050000E+02, 1.100000E+02, 1.150000E+02, 1.200000E+02, 1.250000E+02 }; - static const double y[16] = { 3.478000E+04, 2.861000E+04, 2.365000E+04, 1.963000E+04, 1.637000E+04, 1.372000E+04, 1.154000E+04, 9.744000E+03, 8.261000E+03, 7.030000E+03, 6.005000E+03, 5.147000E+03, 4.427000E+03, 3.820000E+03, 3.307000E+03, 2.872000E+03 }; - int i; - assert(b.size()==3); assert(fvec.size()==16); + for(int i=0; i<16; i++) + fvec[i] = b[0] * exp(b[1]/(x[i]+b[2])) - y[i]; + return 0; + } + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==3); assert(fjac.rows()==16); assert(fjac.cols()==3); - if (iflag != 2) {// compute fvec at b - for(i=0; i<16; i++) { - fvec[i] = b[0] * exp(b[1]/(x[i]+b[2])) - y[i]; - } - } - else { // compute fjac at b - for(i=0; i<16; i++) { - double factor = 1./(x[i]+b[2]); - double e = exp(b[1]*factor); - fjac(i,0) = e; - fjac(i,1) = b[0]*factor*e; - fjac(i,2) = -b[1]*b[0]*factor*factor*e; - } + for(int i=0; i<16; i++) { + double factor = 1./(x[i]+b[2]); + double e = exp(b[1]*factor); + fjac(i,0) = e; + fjac(i,1) = b[0]*factor*e; + fjac(i,2) = -b[1]*b[0]*factor*factor*e; } return 0; } }; +const double MGH10_functor::x[16] = { 5.000000E+01, 5.500000E+01, 6.000000E+01, 6.500000E+01, 7.000000E+01, 7.500000E+01, 8.000000E+01, 8.500000E+01, 9.000000E+01, 9.500000E+01, 1.000000E+02, 1.050000E+02, 1.100000E+02, 1.150000E+02, 1.200000E+02, 1.250000E+02 }; +const double MGH10_functor::y[16] = { 3.478000E+04, 2.861000E+04, 2.365000E+04, 1.963000E+04, 1.637000E+04, 1.372000E+04, 1.154000E+04, 9.744000E+03, 8.261000E+03, 7.030000E+03, 6.005000E+03, 5.147000E+03, 4.427000E+03, 3.820000E+03, 3.307000E+03, 2.872000E+03 }; // http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml void testNistMGH10(void) @@ -1151,31 +1164,31 @@ void testNistMGH10(void) struct BoxBOD_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double x[6]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double x[6] = { 1., 2., 3., 5., 7., 10. }; static const double y[6] = { 109., 149., 149., 191., 213., 224. }; - int i; - assert(b.size()==2); assert(fvec.size()==6); + for(int i=0; i<6; i++) + fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i]; + return 0; + } + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==2); assert(fjac.rows()==6); assert(fjac.cols()==2); - if (iflag != 2) {// compute fvec at b - for(i=0; i<6; i++) { - fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i]; - } - } - else { // compute fjac at b - for(i=0; i<6; i++) { - double e = exp(-b[1]*x[i]); - fjac(i,0) = 1.-e; - fjac(i,1) = b[0]*x[i]*e; - } + for(int i=0; i<6; i++) { + double e = exp(-b[1]*x[i]); + fjac(i,0) = 1.-e; + fjac(i,1) = b[0]*x[i]*e; } return 0; } }; +const double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. }; // http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml void testNistBoxBOD(void) @@ -1225,33 +1238,34 @@ void testNistBoxBOD(void) } struct MGH17_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double x[33]; + static const double y[33]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double x[33] = { 0.000000E+00, 1.000000E+01, 2.000000E+01, 3.000000E+01, 4.000000E+01, 5.000000E+01, 6.000000E+01, 7.000000E+01, 8.000000E+01, 9.000000E+01, 1.000000E+02, 1.100000E+02, 1.200000E+02, 1.300000E+02, 1.400000E+02, 1.500000E+02, 1.600000E+02, 1.700000E+02, 1.800000E+02, 1.900000E+02, 2.000000E+02, 2.100000E+02, 2.200000E+02, 2.300000E+02, 2.400000E+02, 2.500000E+02, 2.600000E+02, 2.700000E+02, 2.800000E+02, 2.900000E+02, 3.000000E+02, 3.100000E+02, 3.200000E+02 }; - static const double y[33] = { 8.440000E-01, 9.080000E-01, 9.320000E-01, 9.360000E-01, 9.250000E-01, 9.080000E-01, 8.810000E-01, 8.500000E-01, 8.180000E-01, 7.840000E-01, 7.510000E-01, 7.180000E-01, 6.850000E-01, 6.580000E-01, 6.280000E-01, 6.030000E-01, 5.800000E-01, 5.580000E-01, 5.380000E-01, 5.220000E-01, 5.060000E-01, 4.900000E-01, 4.780000E-01, 4.670000E-01, 4.570000E-01, 4.480000E-01, 4.380000E-01, 4.310000E-01, 4.240000E-01, 4.200000E-01, 4.140000E-01, 4.110000E-01, 4.060000E-01 }; - int i; - assert(b.size()==5); assert(fvec.size()==33); + for(int i=0; i<33; i++) + fvec[i] = b[0] + b[1]*exp(-b[3]*x[i]) + b[2]*exp(-b[4]*x[i]) - y[i]; + return 0; + } + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==5); assert(fjac.rows()==33); assert(fjac.cols()==5); - if (iflag != 2) {// compute fvec at b - for(i=0; i<33; i++) { - fvec[i] = b[0] + b[1]*exp(-b[3]*x[i]) + b[2]*exp(-b[4]*x[i]) - y[i]; - } - } - else { // compute fjac at b - for(i=0; i<33; i++) { - fjac(i,0) = 1.; - fjac(i,1) = exp(-b[3]*x[i]); - fjac(i,2) = exp(-b[4]*x[i]); - fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]); - fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]); - } + for(int i=0; i<33; i++) { + fjac(i,0) = 1.; + fjac(i,1) = exp(-b[3]*x[i]); + fjac(i,2) = exp(-b[4]*x[i]); + fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]); + fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]); } return 0; } }; +const double MGH17_functor::x[33] = { 0.000000E+00, 1.000000E+01, 2.000000E+01, 3.000000E+01, 4.000000E+01, 5.000000E+01, 6.000000E+01, 7.000000E+01, 8.000000E+01, 9.000000E+01, 1.000000E+02, 1.100000E+02, 1.200000E+02, 1.300000E+02, 1.400000E+02, 1.500000E+02, 1.600000E+02, 1.700000E+02, 1.800000E+02, 1.900000E+02, 2.000000E+02, 2.100000E+02, 2.200000E+02, 2.300000E+02, 2.400000E+02, 2.500000E+02, 2.600000E+02, 2.700000E+02, 2.800000E+02, 2.900000E+02, 3.000000E+02, 3.100000E+02, 3.200000E+02 }; +const double MGH17_functor::y[33] = { 8.440000E-01, 9.080000E-01, 9.320000E-01, 9.360000E-01, 9.250000E-01, 9.080000E-01, 8.810000E-01, 8.500000E-01, 8.180000E-01, 7.840000E-01, 7.510000E-01, 7.180000E-01, 6.850000E-01, 6.580000E-01, 6.280000E-01, 6.030000E-01, 5.800000E-01, 5.580000E-01, 5.380000E-01, 5.220000E-01, 5.060000E-01, 4.900000E-01, 4.780000E-01, 4.670000E-01, 4.570000E-01, 4.480000E-01, 4.380000E-01, 4.310000E-01, 4.240000E-01, 4.200000E-01, 4.140000E-01, 4.110000E-01, 4.060000E-01 }; // http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml void testNistMGH17(void) @@ -1309,35 +1323,37 @@ void testNistMGH17(void) } struct MGH09_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double _x[11]; + static const double y[11]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double _x[11] = { 4., 2., 1., 5.E-1 , 2.5E-01, 1.670000E-01, 1.250000E-01, 1.E-01, 8.330000E-02, 7.140000E-02, 6.250000E-02 }; - static const double y[11] = { 1.957000E-01, 1.947000E-01, 1.735000E-01, 1.600000E-01, 8.440000E-02, 6.270000E-02, 4.560000E-02, 3.420000E-02, 3.230000E-02, 2.350000E-02, 2.460000E-02 }; - int i; - assert(b.size()==4); assert(fvec.size()==11); + for(int i=0; i<11; i++) { + double x = _x[i], xx=x*x; + fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i]; + } + return 0; + } + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==4); assert(fjac.rows()==11); assert(fjac.cols()==4); - if (iflag != 2) {// compute fvec at b - for(i=0; i<11; i++) { - double x = _x[i], xx=x*x; - fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i]; - } - } - else { // compute fjac at b - for(i=0; i<11; i++) { - double x = _x[i], xx=x*x; - double factor = 1./(xx+x*b[2]+b[3]); - fjac(i,0) = (xx+x*b[1]) * factor; - fjac(i,1) = b[0]*x* factor; - fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor; - fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor; - } + for(int i=0; i<11; i++) { + double x = _x[i], xx=x*x; + double factor = 1./(xx+x*b[2]+b[3]); + fjac(i,0) = (xx+x*b[1]) * factor; + fjac(i,1) = b[0]*x* factor; + fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor; + fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor; } return 0; } }; +const double MGH09_functor::_x[11] = { 4., 2., 1., 5.E-1 , 2.5E-01, 1.670000E-01, 1.250000E-01, 1.E-01, 8.330000E-02, 7.140000E-02, 6.250000E-02 }; +const double MGH09_functor::y[11] = { 1.957000E-01, 1.947000E-01, 1.735000E-01, 1.600000E-01, 8.440000E-02, 6.270000E-02, 4.560000E-02, 3.420000E-02, 3.230000E-02, 2.350000E-02, 2.460000E-02 }; // http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml void testNistMGH09(void) @@ -1393,33 +1409,33 @@ void testNistMGH09(void) struct Bennett5_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double x[154]; + static const double y[154]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double x[154] = { 7.447168E0, 8.102586E0, 8.452547E0, 8.711278E0, 8.916774E0, 9.087155E0, 9.232590E0, 9.359535E0, 9.472166E0, 9.573384E0, 9.665293E0, 9.749461E0, 9.827092E0, 9.899128E0, 9.966321E0, 10.029280E0, 10.088510E0, 10.144430E0, 10.197380E0, 10.247670E0, 10.295560E0, 10.341250E0, 10.384950E0, 10.426820E0, 10.467000E0, 10.505640E0, 10.542830E0, 10.578690E0, 10.613310E0, 10.646780E0, 10.679150E0, 10.710520E0, 10.740920E0, 10.770440E0, 10.799100E0, 10.826970E0, 10.854080E0, 10.880470E0, 10.906190E0, 10.931260E0, 10.955720E0, 10.979590E0, 11.002910E0, 11.025700E0, 11.047980E0, 11.069770E0, 11.091100E0, 11.111980E0, 11.132440E0, 11.152480E0, 11.172130E0, 11.191410E0, 11.210310E0, 11.228870E0, 11.247090E0, 11.264980E0, 11.282560E0, 11.299840E0, 11.316820E0, 11.333520E0, 11.349940E0, 11.366100E0, 11.382000E0, 11.397660E0, 11.413070E0, 11.428240E0, 11.443200E0, 11.457930E0, 11.472440E0, 11.486750E0, 11.500860E0, 11.514770E0, 11.528490E0, 11.542020E0, 11.555380E0, 11.568550E0, 11.581560E0, 11.594420E0, 11.607121E0, 11.619640E0, 11.632000E0, 11.644210E0, 11.656280E0, 11.668200E0, 11.679980E0, 11.691620E0, 11.703130E0, 11.714510E0, 11.725760E0, 11.736880E0, 11.747890E0, 11.758780E0, 11.769550E0, 11.780200E0, 11.790730E0, 11.801160E0, 11.811480E0, 11.821700E0, 11.831810E0, 11.841820E0, 11.851730E0, 11.861550E0, 11.871270E0, 11.880890E0, 11.890420E0, 11.899870E0, 11.909220E0, 11.918490E0, 11.927680E0, 11.936780E0, 11.945790E0, 11.954730E0, 11.963590E0, 11.972370E0, 11.981070E0, 11.989700E0, 11.998260E0, 12.006740E0, 12.015150E0, 12.023490E0, 12.031760E0, 12.039970E0, 12.048100E0, 12.056170E0, 12.064180E0, 12.072120E0, 12.080010E0, 12.087820E0, 12.095580E0, 12.103280E0, 12.110920E0, 12.118500E0, 12.126030E0, 12.133500E0, 12.140910E0, 12.148270E0, 12.155570E0, 12.162830E0, 12.170030E0, 12.177170E0, 12.184270E0, 12.191320E0, 12.198320E0, 12.205270E0, 12.212170E0, 12.219030E0, 12.225840E0, 12.232600E0, 12.239320E0, 12.245990E0, 12.252620E0, 12.259200E0, 12.265750E0, 12.272240E0 }; - static const double y[154] = { -34.834702E0 ,-34.393200E0 ,-34.152901E0 ,-33.979099E0 ,-33.845901E0 ,-33.732899E0 ,-33.640301E0 ,-33.559200E0 ,-33.486801E0 ,-33.423100E0 ,-33.365101E0 ,-33.313000E0 ,-33.260899E0 ,-33.217400E0 ,-33.176899E0 ,-33.139198E0 ,-33.101601E0 ,-33.066799E0 ,-33.035000E0 ,-33.003101E0 ,-32.971298E0 ,-32.942299E0 ,-32.916302E0 ,-32.890202E0 ,-32.864101E0 ,-32.841000E0 ,-32.817799E0 ,-32.797501E0 ,-32.774300E0 ,-32.757000E0 ,-32.733799E0 ,-32.716400E0 ,-32.699100E0 ,-32.678799E0 ,-32.661400E0 ,-32.644001E0 ,-32.626701E0 ,-32.612202E0 ,-32.597698E0 ,-32.583199E0 ,-32.568699E0 ,-32.554298E0 ,-32.539799E0 ,-32.525299E0 ,-32.510799E0 ,-32.499199E0 ,-32.487598E0 ,-32.473202E0 ,-32.461601E0 ,-32.435501E0 ,-32.435501E0 ,-32.426800E0 ,-32.412300E0 ,-32.400799E0 ,-32.392101E0 ,-32.380501E0 ,-32.366001E0 ,-32.357300E0 ,-32.348598E0 ,-32.339901E0 ,-32.328400E0 ,-32.319698E0 ,-32.311001E0 ,-32.299400E0 ,-32.290699E0 ,-32.282001E0 ,-32.273300E0 ,-32.264599E0 ,-32.256001E0 ,-32.247299E0 ,-32.238602E0 ,-32.229900E0 ,-32.224098E0 ,-32.215401E0 ,-32.203800E0 ,-32.198002E0 ,-32.189400E0 ,-32.183601E0 ,-32.174900E0 ,-32.169102E0 ,-32.163300E0 ,-32.154598E0 ,-32.145901E0 ,-32.140099E0 ,-32.131401E0 ,-32.125599E0 ,-32.119801E0 ,-32.111198E0 ,-32.105400E0 ,-32.096699E0 ,-32.090900E0 ,-32.088001E0 ,-32.079300E0 ,-32.073502E0 ,-32.067699E0 ,-32.061901E0 ,-32.056099E0 ,-32.050301E0 ,-32.044498E0 ,-32.038799E0 ,-32.033001E0 ,-32.027199E0 ,-32.024300E0 ,-32.018501E0 ,-32.012699E0 ,-32.004002E0 ,-32.001099E0 ,-31.995300E0 ,-31.989500E0 ,-31.983700E0 ,-31.977900E0 ,-31.972099E0 ,-31.969299E0 ,-31.963501E0 ,-31.957701E0 ,-31.951900E0 ,-31.946100E0 ,-31.940300E0 ,-31.937401E0 ,-31.931601E0 ,-31.925800E0 ,-31.922899E0 ,-31.917101E0 ,-31.911301E0 ,-31.908400E0 ,-31.902599E0 ,-31.896900E0 ,-31.893999E0 ,-31.888201E0 ,-31.885300E0 ,-31.882401E0 ,-31.876600E0 ,-31.873699E0 ,-31.867901E0 ,-31.862101E0 ,-31.859200E0 ,-31.856300E0 ,-31.850500E0 ,-31.844700E0 ,-31.841801E0 ,-31.838900E0 ,-31.833099E0 ,-31.830200E0 ,-31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 }; - - int i; - assert(b.size()==3); assert(fvec.size()==154); + for(int i=0; i<154; i++) + fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i]; + return 0; + } + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==3); assert(fjac.rows()==154); assert(fjac.cols()==3); - if (iflag != 2) {// compute fvec at b - for(i=0; i<154; i++) { - fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i]; - } - } - else { // compute fjac at b - for(i=0; i<154; i++) { - double e = pow(b[1]+x[i],-1./b[2]); - fjac(i,0) = e; - fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]); - fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2]; - } + for(int i=0; i<154; i++) { + double e = pow(b[1]+x[i],-1./b[2]); + fjac(i,0) = e; + fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]); + fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2]; } return 0; } }; +const double Bennett5_functor::x[154] = { 7.447168E0, 8.102586E0, 8.452547E0, 8.711278E0, 8.916774E0, 9.087155E0, 9.232590E0, 9.359535E0, 9.472166E0, 9.573384E0, 9.665293E0, 9.749461E0, 9.827092E0, 9.899128E0, 9.966321E0, 10.029280E0, 10.088510E0, 10.144430E0, 10.197380E0, 10.247670E0, 10.295560E0, 10.341250E0, 10.384950E0, 10.426820E0, 10.467000E0, 10.505640E0, 10.542830E0, 10.578690E0, 10.613310E0, 10.646780E0, 10.679150E0, 10.710520E0, 10.740920E0, 10.770440E0, 10.799100E0, 10.826970E0, 10.854080E0, 10.880470E0, 10.906190E0, 10.931260E0, 10.955720E0, 10.979590E0, 11.002910E0, 11.025700E0, 11.047980E0, 11.069770E0, 11.091100E0, 11.111980E0, 11.132440E0, 11.152480E0, 11.172130E0, 11.191410E0, 11.210310E0, 11.228870E0, 11.247090E0, 11.264980E0, 11.282560E0, 11.299840E0, 11.316820E0, 11.333520E0, 11.349940E0, 11.366100E0, 11.382000E0, 11.397660E0, 11.413070E0, 11.428240E0, 11.443200E0, 11.457930E0, 11.472440E0, 11.486750E0, 11.500860E0, 11.514770E0, 11.528490E0, 11.542020E0, 11.555380E0, 11.568550E0, 11.581560E0, 11.594420E0, 11.607121E0, 11.619640E0, 11.632000E0, 11.644210E0, 11.656280E0, 11.668200E0, 11.679980E0, 11.691620E0, 11.703130E0, 11.714510E0, 11.725760E0, 11.736880E0, 11.747890E0, 11.758780E0, 11.769550E0, 11.780200E0, 11.790730E0, 11.801160E0, 11.811480E0, 11.821700E0, 11.831810E0, 11.841820E0, 11.851730E0, 11.861550E0, 11.871270E0, 11.880890E0, 11.890420E0, 11.899870E0, 11.909220E0, 11.918490E0, 11.927680E0, 11.936780E0, 11.945790E0, 11.954730E0, 11.963590E0, 11.972370E0, 11.981070E0, 11.989700E0, 11.998260E0, 12.006740E0, 12.015150E0, 12.023490E0, 12.031760E0, 12.039970E0, 12.048100E0, 12.056170E0, 12.064180E0, 12.072120E0, 12.080010E0, 12.087820E0, 12.095580E0, 12.103280E0, 12.110920E0, 12.118500E0, 12.126030E0, 12.133500E0, 12.140910E0, 12.148270E0, 12.155570E0, 12.162830E0, 12.170030E0, 12.177170E0, 12.184270E0, 12.191320E0, 12.198320E0, 12.205270E0, 12.212170E0, 12.219030E0, 12.225840E0, 12.232600E0, 12.239320E0, 12.245990E0, 12.252620E0, 12.259200E0, 12.265750E0, 12.272240E0 }; +const double Bennett5_functor::y[154] = { -34.834702E0 ,-34.393200E0 ,-34.152901E0 ,-33.979099E0 ,-33.845901E0 ,-33.732899E0 ,-33.640301E0 ,-33.559200E0 ,-33.486801E0 ,-33.423100E0 ,-33.365101E0 ,-33.313000E0 ,-33.260899E0 ,-33.217400E0 ,-33.176899E0 ,-33.139198E0 ,-33.101601E0 ,-33.066799E0 ,-33.035000E0 ,-33.003101E0 ,-32.971298E0 ,-32.942299E0 ,-32.916302E0 ,-32.890202E0 ,-32.864101E0 ,-32.841000E0 ,-32.817799E0 ,-32.797501E0 ,-32.774300E0 ,-32.757000E0 ,-32.733799E0 ,-32.716400E0 ,-32.699100E0 ,-32.678799E0 ,-32.661400E0 ,-32.644001E0 ,-32.626701E0 ,-32.612202E0 ,-32.597698E0 ,-32.583199E0 ,-32.568699E0 ,-32.554298E0 ,-32.539799E0 ,-32.525299E0 ,-32.510799E0 ,-32.499199E0 ,-32.487598E0 ,-32.473202E0 ,-32.461601E0 ,-32.435501E0 ,-32.435501E0 ,-32.426800E0 ,-32.412300E0 ,-32.400799E0 ,-32.392101E0 ,-32.380501E0 ,-32.366001E0 ,-32.357300E0 ,-32.348598E0 ,-32.339901E0 ,-32.328400E0 ,-32.319698E0 ,-32.311001E0 ,-32.299400E0 ,-32.290699E0 ,-32.282001E0 ,-32.273300E0 ,-32.264599E0 ,-32.256001E0 ,-32.247299E0 ,-32.238602E0 ,-32.229900E0 ,-32.224098E0 ,-32.215401E0 ,-32.203800E0 ,-32.198002E0 ,-32.189400E0 ,-32.183601E0 ,-32.174900E0 ,-32.169102E0 ,-32.163300E0 ,-32.154598E0 ,-32.145901E0 ,-32.140099E0 ,-32.131401E0 ,-32.125599E0 ,-32.119801E0 ,-32.111198E0 ,-32.105400E0 ,-32.096699E0 ,-32.090900E0 ,-32.088001E0 ,-32.079300E0 ,-32.073502E0 ,-32.067699E0 ,-32.061901E0 ,-32.056099E0 ,-32.050301E0 ,-32.044498E0 ,-32.038799E0 ,-32.033001E0 ,-32.027199E0 ,-32.024300E0 ,-32.018501E0 ,-32.012699E0 ,-32.004002E0 ,-32.001099E0 ,-31.995300E0 ,-31.989500E0 ,-31.983700E0 ,-31.977900E0 ,-31.972099E0 ,-31.969299E0 ,-31.963501E0 ,-31.957701E0 ,-31.951900E0 ,-31.946100E0 ,-31.940300E0 ,-31.937401E0 ,-31.931601E0 ,-31.925800E0 ,-31.922899E0 ,-31.917101E0 ,-31.911301E0 ,-31.908400E0 ,-31.902599E0 ,-31.896900E0 ,-31.893999E0 ,-31.888201E0 ,-31.885300E0 ,-31.882401E0 ,-31.876600E0 ,-31.873699E0 ,-31.867901E0 ,-31.862101E0 ,-31.859200E0 ,-31.856300E0 ,-31.850500E0 ,-31.844700E0 ,-31.841801E0 ,-31.838900E0 ,-31.833099E0 ,-31.830200E0 ,-31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 }; // http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml void testNistBennett5(void) @@ -1469,40 +1485,42 @@ void testNistBennett5(void) } struct thurber_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double _x[37]; + static const double _y[37]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double _x[37] = { -3.067E0, -2.981E0, -2.921E0, -2.912E0, -2.840E0, -2.797E0, -2.702E0, -2.699E0, -2.633E0, -2.481E0, -2.363E0, -2.322E0, -1.501E0, -1.460E0, -1.274E0, -1.212E0, -1.100E0, -1.046E0, -0.915E0, -0.714E0, -0.566E0, -0.545E0, -0.400E0, -0.309E0, -0.109E0, -0.103E0, 0.010E0, 0.119E0, 0.377E0, 0.790E0, 0.963E0, 1.006E0, 1.115E0, 1.572E0, 1.841E0, 2.047E0, 2.200E0 }; - static const double _y[37] = { 80.574E0, 84.248E0, 87.264E0, 87.195E0, 89.076E0, 89.608E0, 89.868E0, 90.101E0, 92.405E0, 95.854E0, 100.696E0, 101.060E0, 401.672E0, 390.724E0, 567.534E0, 635.316E0, 733.054E0, 759.087E0, 894.206E0, 990.785E0, 1090.109E0, 1080.914E0, 1122.643E0, 1178.351E0, 1260.531E0, 1273.514E0, 1288.339E0, 1327.543E0, 1353.863E0, 1414.509E0, 1425.208E0, 1421.384E0, 1442.962E0, 1464.350E0, 1468.705E0, 1447.894E0, 1457.628E0}; - int i; - -// static int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; + // static int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; assert(b.size()==7); assert(fvec.size()==37); + for(int i=0; i<37; i++) { + double x=_x[i], xx=x*x, xxx=xx*x; + fvec[i] = (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) / (1.+b[4]*x+b[5]*xx+b[6]*xxx) - _y[i]; + } + return 0; + } + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==7); assert(fjac.rows()==37); assert(fjac.cols()==7); - if (iflag != 2) {// compute fvec at x - for(i=0; i<37; i++) { - double x=_x[i], xx=x*x, xxx=xx*x; - fvec[i] = (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) / (1.+b[4]*x+b[5]*xx+b[6]*xxx) - _y[i]; - } - } - else { // compute fjac at x - for(i=0; i<37; i++) { - double x=_x[i], xx=x*x, xxx=xx*x; - double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); - fjac(i,0) = 1.*fact; - fjac(i,1) = x*fact; - fjac(i,2) = xx*fact; - fjac(i,3) = xxx*fact; - fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; - fjac(i,4) = x*fact; - fjac(i,5) = xx*fact; - fjac(i,6) = xxx*fact; - } + for(int i=0; i<37; i++) { + double x=_x[i], xx=x*x, xxx=xx*x; + double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); + fjac(i,0) = 1.*fact; + fjac(i,1) = x*fact; + fjac(i,2) = xx*fact; + fjac(i,3) = xxx*fact; + fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; + fjac(i,4) = x*fact; + fjac(i,5) = xx*fact; + fjac(i,6) = xxx*fact; } return 0; } }; +const double thurber_functor::_x[37] = { -3.067E0, -2.981E0, -2.921E0, -2.912E0, -2.840E0, -2.797E0, -2.702E0, -2.699E0, -2.633E0, -2.481E0, -2.363E0, -2.322E0, -1.501E0, -1.460E0, -1.274E0, -1.212E0, -1.100E0, -1.046E0, -0.915E0, -0.714E0, -0.566E0, -0.545E0, -0.400E0, -0.309E0, -0.109E0, -0.103E0, 0.010E0, 0.119E0, 0.377E0, 0.790E0, 0.963E0, 1.006E0, 1.115E0, 1.572E0, 1.841E0, 2.047E0, 2.200E0 }; +const double thurber_functor::_y[37] = { 80.574E0, 84.248E0, 87.264E0, 87.195E0, 89.076E0, 89.608E0, 89.868E0, 90.101E0, 92.405E0, 95.854E0, 100.696E0, 101.060E0, 401.672E0, 390.724E0, 567.534E0, 635.316E0, 733.054E0, 759.087E0, 894.206E0, 990.785E0, 1090.109E0, 1080.914E0, 1122.643E0, 1178.351E0, 1260.531E0, 1273.514E0, 1288.339E0, 1327.543E0, 1353.863E0, 1414.509E0, 1425.208E0, 1421.384E0, 1442.962E0, 1464.350E0, 1468.705E0, 1447.894E0, 1457.628E0}; // http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml void testNistThurber(void) @@ -1562,34 +1580,35 @@ void testNistThurber(void) } struct rat43_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double x[15]; + static const double y[15]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. }; - static const double y[15] = { 16.08, 33.83, 65.80, 97.20, 191.55, 326.20, 386.87, 520.53, 590.03, 651.92, 724.93, 699.56, 689.96, 637.56, 717.41 }; - int i; - assert(b.size()==4); assert(fvec.size()==15); + for(int i=0; i<15; i++) + fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i]; + return 0; + } + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==4); assert(fjac.rows()==15); assert(fjac.cols()==4); - if (iflag != 2) {// compute fvec at b - for(i=0; i<15; i++) { - fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i]; - } - } - else { // compute fjac at b - for(i=0; i<15; i++) { - double e = exp(b[1]-b[2]*x[i]); - double power = -1./b[3]; - fjac(i,0) = pow(1.+e, power); - fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.); - fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.); - fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power); - } + for(int i=0; i<15; i++) { + double e = exp(b[1]-b[2]*x[i]); + double power = -1./b[3]; + fjac(i,0) = pow(1.+e, power); + fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.); + fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.); + fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power); } return 0; } }; +const double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. }; +const double rat43_functor::y[15] = { 16.08, 33.83, 65.80, 97.20, 191.55, 326.20, 386.87, 520.53, 590.03, 651.92, 724.93, 699.56, 689.96, 637.56, 717.41 }; // http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml void testNistRat43(void) @@ -1645,34 +1664,34 @@ void testNistRat43(void) struct eckerle4_functor { - static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag) + static const double x[35]; + static const double y[35]; + static int debug(const VectorXd & /* x */, const VectorXd & /* fvec */, const MatrixXd & /* fjac */) { return 0;} + static int f(const VectorXd &b, VectorXd &fvec) { - static const double x[35] = { 400.0, 405.0, 410.0, 415.0, 420.0, 425.0, 430.0, 435.0, 436.5, 438.0, 439.5, 441.0, 442.5, 444.0, 445.5, 447.0, 448.5, 450.0, 451.5, 453.0, 454.5, 456.0, 457.5, 459.0, 460.5, 462.0, 463.5, 465.0, 470.0, 475.0, 480.0, 485.0, 490.0, 495.0, 500.0}; - static const double y[35] = { 0.0001575, 0.0001699, 0.0002350, 0.0003102, 0.0004917, 0.0008710, 0.0017418, 0.0046400, 0.0065895, 0.0097302, 0.0149002, 0.0237310, 0.0401683, 0.0712559, 0.1264458, 0.2073413, 0.2902366, 0.3445623, 0.3698049, 0.3668534, 0.3106727, 0.2078154, 0.1164354, 0.0616764, 0.0337200, 0.0194023, 0.0117831, 0.0074357, 0.0022732, 0.0008800, 0.0004579, 0.0002345, 0.0001586, 0.0001143, 0.0000710 }; - - int i; - assert(b.size()==3); assert(fvec.size()==35); + for(int i=0; i<35; i++) + fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i]; + return 0; + } + static int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==3); assert(fjac.rows()==35); assert(fjac.cols()==3); - if (iflag != 2) {// compute fvec at b - for(i=0; i<35; i++) { - fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i]; - } - } - else { // compute fjac at b - for(i=0; i<35; i++) { - double b12 = b[1]*b[1]; - double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12); - fjac(i,0) = e / b[1]; - fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12; - fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12; - } + for(int i=0; i<35; i++) { + double b12 = b[1]*b[1]; + double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12); + fjac(i,0) = e / b[1]; + fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12; + fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12; } return 0; } }; +const double eckerle4_functor::x[35] = { 400.0, 405.0, 410.0, 415.0, 420.0, 425.0, 430.0, 435.0, 436.5, 438.0, 439.5, 441.0, 442.5, 444.0, 445.5, 447.0, 448.5, 450.0, 451.5, 453.0, 454.5, 456.0, 457.5, 459.0, 460.5, 462.0, 463.5, 465.0, 470.0, 475.0, 480.0, 485.0, 490.0, 495.0, 500.0}; +const double eckerle4_functor::y[35] = { 0.0001575, 0.0001699, 0.0002350, 0.0003102, 0.0004917, 0.0008710, 0.0017418, 0.0046400, 0.0065895, 0.0097302, 0.0149002, 0.0237310, 0.0401683, 0.0712559, 0.1264458, 0.2073413, 0.2902366, 0.3445623, 0.3698049, 0.3668534, 0.3106727, 0.2078154, 0.1164354, 0.0616764, 0.0337200, 0.0194023, 0.0117831, 0.0074357, 0.0022732, 0.0008800, 0.0004579, 0.0002345, 0.0001586, 0.0001143, 0.0000710 }; // http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml void testNistEckerle4(void)