9 #include <unsupported/Eigen/NonLinearOptimization>
16 #define LM_EVAL_COUNT_TOL 4/3
18 #define LM_CHECK_N_ITERS(SOLVER,NFEV,NJEV) { \
20 VERIFY_IS_EQUAL(SOLVER.nfev, NFEV); \
21 VERIFY_IS_EQUAL(SOLVER.njev, NJEV); \
23 VERIFY(SOLVER.nfev <= NFEV * LM_EVAL_COUNT_TOL); \
24 VERIFY(SOLVER.njev <= NJEV * LM_EVAL_COUNT_TOL); \
27 int fcn_chkder(
const VectorXd &
x, VectorXd &fvec, MatrixXd &fjac,
int iflag)
32 assert(15 == fvec.size());
33 assert(3 ==
x.size());
34 double tmp1, tmp2, tmp3, tmp4;
35 static const double y[15]={1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1,
36 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
43 for (
i=0;
i<15;
i++) {
47 if (
i >= 8) tmp3 = tmp2;
48 fvec[
i] =
y[
i] - (
x[0] + tmp1/(
x[1]*tmp2 +
x[2]*tmp3));
51 for (
i = 0;
i < 15;
i++) {
59 if (
i >= 8) tmp3 = tmp2;
60 tmp4 = (
x[1]*tmp2 +
x[2]*tmp3); tmp4=tmp4*tmp4;
62 fjac(
i,1) = tmp1*tmp2/tmp4;
63 fjac(
i,2) = tmp1*tmp3/tmp4;
73 VectorXd
x(
n), fvec(
m), xp, fvecp(
m), err;
79 x << 9.2e-1, 1.3e-1, 5.4e-1;
90 VectorXd fvec_ref(
m), fvecp_ref(
m), err_ref(
m);
92 -1.181606, -1.429655, -1.606344,
93 -1.745269, -1.840654, -1.921586,
94 -1.984141, -2.022537, -2.468977,
95 -2.827562, -3.473582, -4.437612,
96 -6.047662, -9.267761, -18.91806;
98 -7.724666e-09, -3.432406e-09, -2.034843e-10,
99 2.313685e-09, 4.331078e-09, 5.984096e-09,
100 7.363281e-09, 8.53147e-09, 1.488591e-08,
101 2.33585e-08, 3.522012e-08, 5.301255e-08,
102 8.26666e-08, 1.419747e-07, 3.19899e-07;
104 0.1141397, 0.09943516, 0.09674474,
105 0.09980447, 0.1073116, 0.1220445,
116 template<
typename _Scalar,
int NX=Dynamic,
int NY=Dynamic>
145 double tmp1, tmp2, tmp3;
146 static const double y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1,
147 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
153 tmp3 = (
i>=8)? tmp2 : tmp1;
154 fvec[
i] =
y[
i] - (
x[0] + tmp1/(
x[1]*tmp2 +
x[2]*tmp3));
159 int df(
const VectorXd &
x, MatrixXd &fjac)
const
161 double tmp1, tmp2, tmp3, tmp4;
166 tmp3 = (
i>=8)? tmp2 : tmp1;
167 tmp4 = (
x[1]*tmp2 +
x[2]*tmp3); tmp4 = tmp4*tmp4;
169 fjac(
i,1) = tmp1*tmp2/tmp4;
170 fjac(
i,2) = tmp1*tmp3/tmp4;
183 x.setConstant(
n, 1.);
199 x_ref << 0.08241058, 1.133037, 2.343695;
207 double fnorm, covfac;
211 x.setConstant(
n, 1.);
223 fnorm = lm.
fvec.blueNorm();
228 x_ref << 0.08241058, 1.133037, 2.343695;
232 covfac = fnorm*fnorm/(
m-
n);
235 MatrixXd cov_ref(
n,
n);
237 0.0001531202, 0.002869941, -0.002656662,
238 0.002869941, 0.09480935, -0.09098995,
239 -0.002656662, -0.09098995, 0.08778727;
244 cov = covfac*lm.
fjac.topLeftCorner<
n,
n>();
256 double temp, temp1, temp2;
258 assert(fvec.size()==
n);
261 temp = (3. - 2.*
x[k])*
x[k];
263 if (k) temp1 =
x[k-1];
265 if (k !=
n-1) temp2 =
x[k+1];
266 fvec[k] = temp - temp1 - 2.*temp2 + 1.;
270 int df(
const VectorXd &
x, MatrixXd &fjac)
273 assert(fjac.rows()==
n);
274 assert(fjac.cols()==
n);
279 fjac(k,k) = 3.- 4.*
x[k];
280 if (k) fjac(k,k-1) = -1.;
281 if (k !=
n-1) fjac(k,k+1) = -2.;
295 x.setConstant(
n, -1.);
313 -0.5706545, -0.6816283, -0.7017325,
314 -0.7042129, -0.701369, -0.6918656,
315 -0.665792, -0.5960342, -0.4164121;
326 x.setConstant(
n, -1.);
332 solver.diag.setConstant(
n, 1.);
333 solver.useExternalScaling =
true;
347 -0.5706545, -0.6816283, -0.7017325,
348 -0.7042129, -0.701369, -0.6918656,
349 -0.665792, -0.5960342, -0.4164121;
359 double temp, temp1, temp2;
362 assert(fvec.size()==
n);
365 temp = (3. - 2.*
x[k])*
x[k];
367 if (k) temp1 =
x[k-1];
369 if (k !=
n-1) temp2 =
x[k+1];
370 fvec[k] = temp - temp1 - 2.*temp2 + 1.;
382 x.setConstant(
n, -1.);
398 x_ref << -0.5706545, -0.6816283, -0.7017325, -0.7042129, -0.701369, -0.6918656, -0.665792, -0.5960342, -0.4164121;
409 x.setConstant(
n, -1.);
414 solver.parameters.nb_of_subdiagonals = 1;
415 solver.parameters.nb_of_superdiagonals = 1;
416 solver.diag.setConstant(
n, 1.);
417 solver.useExternalScaling =
true;
430 -0.5706545, -0.6816283, -0.7017325,
431 -0.7042129, -0.701369, -0.6918656,
432 -0.665792, -0.5960342, -0.4164121;
442 double tmp1, tmp2, tmp3;
443 static const double y[15]={1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1,
444 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
446 assert(15==fvec.size());
449 for (
int i=0;
i<15;
i++)
453 tmp3 = (
i>=8)? tmp2 : tmp1;
454 fvec[
i] =
y[
i] - (
x[0] + tmp1/(
x[1]*tmp2 +
x[2]*tmp3));
461 assert(jac_row.size()==
x.size());
462 double tmp1, tmp2, tmp3, tmp4;
467 tmp3 = (
i>=8)? tmp2 : tmp1;
468 tmp4 = (
x[1]*tmp2 +
x[2]*tmp3); tmp4 = tmp4*tmp4;
470 jac_row[1] = tmp1*tmp2/tmp4;
471 jac_row[2] = tmp1*tmp3/tmp4;
484 x.setConstant(
n, 1.);
500 x_ref << 0.08241058, 1.133037, 2.343695 ;
512 x.setConstant(
n, 1.);
524 fnorm = lm.
fvec.blueNorm();
529 x_ref << 0.08241058, 1.133037, 2.343695;
540 double tmp1,tmp2,tmp3;
541 static const 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,
542 3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0};
545 assert(fvec.size()==15);
552 if (
i >= 8) tmp3 = tmp2;
553 fvec[
i] =
y[
i] - (
x[0] + tmp1/(
x[1]*tmp2 +
x[2]*tmp3));
564 VectorXd
x(
n), fvec(15);
567 x.setConstant(
n, 1.);
584 x_ref << 0.0824106, 1.1330366, 2.3436947;
593 double fnorm, covfac;
597 x.setConstant(
n, 1.);
610 fnorm = lm.
fvec.blueNorm();
615 x_ref << 0.08241058, 1.133037, 2.343695;
619 covfac = fnorm*fnorm/(
m-
n);
622 MatrixXd cov_ref(
n,
n);
624 0.0001531202, 0.002869942, -0.002656662,
625 0.002869942, 0.09480937, -0.09098997,
626 -0.002656662, -0.09098997, 0.08778729;
631 cov = covfac*lm.
fjac.topLeftCorner<
n,
n>();
640 static const double m_x[54];
641 static const double m_y[54];
647 assert(fvec.size()==54);
648 for(
i=0;
i<54;
i++) {
654 int df(
const VectorXd &
b, MatrixXd &fjac)
657 assert(fjac.rows()==54);
658 assert(fjac.cols()==3);
659 for(
int i=0;
i<54;
i++) {
661 double factor = 1./(
b[1]+
b[2]*
x);
663 fjac(
i,0) = -
x*
e*factor;
664 fjac(
i,1) = -
e*factor*factor;
665 fjac(
i,2) = -
x*
e*factor*factor;
670 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};
671 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 };
703 x<< 0.15, 0.008, 0.010;
725 static const double m_x[14];
726 static const double m_y[14];
730 assert(fvec.size()==14);
731 for(
int i=0;
i<14;
i++) {
736 int df(
const VectorXd &
b, MatrixXd &fjac)
739 assert(fjac.rows()==14);
740 assert(fjac.cols()==2);
741 for(
int i=0;
i<14;
i++) {
748 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};
749 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};
797 static const double m_x[236];
800 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 ,
801 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 ,
802 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 };
807 assert(fvec.size()==236);
808 for(
int i=0;
i<236;
i++) {
810 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];
815 int df(
const VectorXd &
b, MatrixXd &fjac)
818 assert(fjac.rows()==236);
819 assert(fjac.cols()==7);
820 for(
int i=0;
i<236;
i++) {
822 double fact = 1./(1.+
b[4]*
x+
b[5]*xx+
b[6]*xxx);
826 fjac(
i,3) = xxx*fact;
827 fact = - (
b[0]+
b[1]*
x+
b[2]*xx+
b[3]*xxx) * fact * fact;
830 fjac(
i,6) = xxx*fact;
835 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 ,
836 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 ,
837 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};
850 x<< 10., -1., .05, -.00001, -.05, .001, -.000001;
873 x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001;
896 static const double x[14];
897 static const double y[14];
901 assert(fvec.size()==14);
902 for(
int i=0;
i<14;
i++) {
903 fvec[
i] =
b[0]*
b[1]*
x[
i]/(1.+
b[1]*
x[
i]) -
y[
i];
907 int df(
const VectorXd &
b, MatrixXd &fjac)
910 assert(fjac.rows()==14);
911 assert(fjac.cols()==2);
912 for(
int i=0;
i<14;
i++) {
913 double den = 1.+
b[1]*
x[
i];
914 fjac(
i,0) =
b[1]*
x[
i] / den;
915 fjac(
i,1) =
b[0]*
x[
i]*(den-
b[1]*
x[
i])/den/den;
920 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};
921 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};
970 static const double x[24];
971 static const double y[24];
975 assert(fvec.size()==24);
976 for(
int i=0;
i<24;
i++)
980 int df(
const VectorXd &
b, MatrixXd &fjac)
983 assert(fjac.rows()==24);
984 assert(fjac.cols()==6);
985 for(
int i=0;
i<24;
i++) {
996 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 };
997 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 };
1010 x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6;
1020 std::cout.precision(30);
1021 std::cout << lm.
fvec.squaredNorm() <<
"\n";
1022 VERIFY(lm.
fvec.squaredNorm() <= 1.4307867721E-25);
1034 x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3;
1042 VERIFY(lm.
fvec.squaredNorm() <= 1.4307867721E-25);
1056 static const double x[9];
1057 static const double y[9];
1060 assert(
b.size()==3);
1061 assert(fvec.size()==9);
1062 for(
int i=0;
i<9;
i++) {
1068 int df(
const VectorXd &
b, MatrixXd &fjac)
1070 assert(
b.size()==3);
1071 assert(fjac.rows()==9);
1072 assert(fjac.cols()==3);
1073 for(
int i=0;
i<9;
i++) {
1075 fjac(
i,0) = 1./(1.+
e);
1076 fjac(
i,1) = -
b[0]*
e/(1.+
e)/(1.+
e);
1077 fjac(
i,2) = +
b[0]*
e*
x[
i]/(1.+
e)/(1.+
e);
1082 const double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 };
1083 const double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 };
1133 static const double x[16];
1134 static const double y[16];
1137 assert(
b.size()==3);
1138 assert(fvec.size()==16);
1139 for(
int i=0;
i<16;
i++)
1143 int df(
const VectorXd &
b, MatrixXd &fjac)
1145 assert(
b.size()==3);
1146 assert(fjac.rows()==16);
1147 assert(fjac.cols()==3);
1148 for(
int i=0;
i<16;
i++) {
1149 double factor = 1./(
x[
i]+
b[2]);
1150 double e =
exp(
b[1]*factor);
1152 fjac(
i,1) =
b[0]*factor*
e;
1153 fjac(
i,2) = -
b[1]*
b[0]*factor*factor*
e;
1158 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 };
1159 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 };
1172 x<< 2., 400000., 25000.;
1191 x<< 0.02, 4000., 250.;
1210 static const double x[6];
1213 static const double y[6] = { 109., 149., 149., 191., 213., 224. };
1214 assert(
b.size()==2);
1215 assert(fvec.size()==6);
1216 for(
int i=0;
i<6;
i++)
1220 int df(
const VectorXd &
b, MatrixXd &fjac)
1222 assert(
b.size()==2);
1223 assert(fjac.rows()==6);
1224 assert(fjac.cols()==2);
1225 for(
int i=0;
i<6;
i++) {
1228 fjac(
i,1) =
b[0]*
x[
i]*
e;
1287 static const double x[33];
1288 static const double y[33];
1291 assert(
b.size()==5);
1292 assert(fvec.size()==33);
1293 for(
int i=0;
i<33;
i++)
1297 int df(
const VectorXd &
b, MatrixXd &fjac)
1299 assert(
b.size()==5);
1300 assert(fjac.rows()==33);
1301 assert(fjac.cols()==5);
1302 for(
int i=0;
i<33;
i++) {
1312 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 };
1313 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 };
1326 x<< 50., 150., -100., 1., 2.;
1351 x<< 0.5 ,1.5 ,-1 ,0.01 ,0.02;
1372 static const double _x[11];
1373 static const double y[11];
1376 assert(
b.size()==4);
1377 assert(fvec.size()==11);
1378 for(
int i=0;
i<11;
i++) {
1379 double x =
_x[
i], xx=
x*
x;
1380 fvec[
i] =
b[0]*(xx+
x*
b[1])/(xx+
x*
b[2]+
b[3]) -
y[
i];
1384 int df(
const VectorXd &
b, MatrixXd &fjac)
1386 assert(
b.size()==4);
1387 assert(fjac.rows()==11);
1388 assert(fjac.cols()==4);
1389 for(
int i=0;
i<11;
i++) {
1390 double x =
_x[
i], xx=
x*
x;
1391 double factor = 1./(xx+
x*
b[2]+
b[3]);
1392 fjac(
i,0) = (xx+
x*
b[1]) * factor;
1393 fjac(
i,1) =
b[0]*
x* factor;
1394 fjac(
i,2) = -
b[0]*(xx+
x*
b[1]) *
x * factor * factor;
1395 fjac(
i,3) = -
b[0]*(xx+
x*
b[1]) * factor * factor;
1400 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 };
1401 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 };
1414 x<< 25., 39, 41.5, 39.;
1435 x<< 0.25, 0.39, 0.415, 0.39;
1457 static const double x[154];
1458 static const double y[154];
1461 assert(
b.size()==3);
1462 assert(fvec.size()==154);
1463 for(
int i=0;
i<154;
i++)
1467 int df(
const VectorXd &
b, MatrixXd &fjac)
1469 assert(
b.size()==3);
1470 assert(fjac.rows()==154);
1471 assert(fjac.cols()==3);
1472 for(
int i=0;
i<154;
i++) {
1473 double e =
pow(
b[1]+
x[
i],-1./
b[2]);
1475 fjac(
i,1) = -
b[0]*
e/
b[2]/(
b[1]+
x[
i]);
1481 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,
1482 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 };
1483 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
1484 ,-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 ,
1485 -31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 };
1498 x<< -2000., 50., 0.8;
1517 x<< -1500., 45., 0.85;
1536 static const double _x[37];
1537 static const double _y[37];
1541 assert(
b.size()==7);
1542 assert(fvec.size()==37);
1543 for(
int i=0;
i<37;
i++) {
1544 double x=
_x[
i], xx=
x*
x, xxx=xx*
x;
1545 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];
1549 int df(
const VectorXd &
b, MatrixXd &fjac)
1551 assert(
b.size()==7);
1552 assert(fjac.rows()==37);
1553 assert(fjac.cols()==7);
1554 for(
int i=0;
i<37;
i++) {
1555 double x=
_x[
i], xx=
x*
x, xxx=xx*
x;
1556 double fact = 1./(1.+
b[4]*
x+
b[5]*xx+
b[6]*xxx);
1557 fjac(
i,0) = 1.*fact;
1559 fjac(
i,2) = xx*fact;
1560 fjac(
i,3) = xxx*fact;
1561 fact = - (
b[0]+
b[1]*
x+
b[2]*xx+
b[3]*xxx) * fact * fact;
1563 fjac(
i,5) = xx*fact;
1564 fjac(
i,6) = xxx*fact;
1569 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 };
1570 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};
1583 x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ;
1608 x<< 1300 ,1500 ,500 ,75 ,1 ,0.4 ,0.05 ;
1633 static const double x[15];
1634 static const double y[15];
1637 assert(
b.size()==4);
1638 assert(fvec.size()==15);
1639 for(
int i=0;
i<15;
i++)
1643 int df(
const VectorXd &
b, MatrixXd &fjac)
1645 assert(
b.size()==4);
1646 assert(fjac.rows()==15);
1647 assert(fjac.cols()==4);
1648 for(
int i=0;
i<15;
i++) {
1650 double power = -1./
b[3];
1651 fjac(
i,0) =
pow(1.+
e, power);
1652 fjac(
i,1) = power*
b[0]*
e*
pow(1.+
e, power-1.);
1653 fjac(
i,2) = -power*
b[0]*
e*
x[
i]*
pow(1.+
e, power-1.);
1654 fjac(
i,3) =
b[0]*power*power*
log(1.+
e)*
pow(1.+
e, power);
1659 const double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. };
1660 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 };
1673 x<< 100., 10., 1., 1.;
1695 x<< 700., 5., 0.75, 1.3;
1719 static const double x[35];
1720 static const double y[35];
1723 assert(
b.size()==3);
1724 assert(fvec.size()==35);
1725 for(
int i=0;
i<35;
i++)
1726 fvec[
i] =
b[0]/
b[1] *
exp(-0.5*(
x[
i]-
b[2])*(
x[
i]-
b[2])/(
b[1]*
b[1])) -
y[
i];
1729 int df(
const VectorXd &
b, MatrixXd &fjac)
1731 assert(
b.size()==3);
1732 assert(fjac.rows()==35);
1733 assert(fjac.cols()==3);
1734 for(
int i=0;
i<35;
i++) {
1735 double b12 =
b[1]*
b[1];
1736 double e =
exp(-0.5*(
x[
i]-
b[2])*(
x[
i]-
b[2])/b12);
1737 fjac(
i,0) =
e /
b[1];
1738 fjac(
i,1) = ((
x[
i]-
b[2])*(
x[
i]-
b[2])/b12-1.) *
b[0]*
e/b12;
1739 fjac(
i,2) = (
x[
i]-
b[2])*
e*
b[0]/
b[1]/b12;
1744 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};
1745 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 };