use eigen objects for the lmder callback

1 files changed, 129 insertions, 114 deletions

diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp
index 5a2482163..d9b878c94 100644
--- a/unsupported/test/NonLinear.cpp
+++ b/unsupported/test/NonLinear.cpp
@@ -102,7 +102,7 @@ void testChkder()
 
 
 struct lmder_functor {
-    static int f(int /*m*/, int /*n*/, const double *x, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &x, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
 
         /*      subroutine fcn for lmder1 example. */
@@ -130,9 +130,9 @@ struct lmder_functor {
                 tmp2 = 16 - i - 1;
                 tmp3 = (i>=8)? tmp2 : tmp1;
                 tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4;
-                fjac[i + ldfjac*(0)] = -1.;
-                fjac[i + ldfjac*(1)] = tmp1*tmp2/tmp4;
-                fjac[i + ldfjac*(2)] = tmp1*tmp3/tmp4;
+                fjac(i,0) = -1;
+                fjac(i,1) = tmp1*tmp2/tmp4;
+                fjac(i,2) = tmp1*tmp3/tmp4;
             }
         }
         return 0;
@@ -587,15 +587,16 @@ void testLmdif()
 }
 
 struct chwirut2_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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(m==54);
-        assert(n==3);
-        assert(ldfjac==54);
+        assert(b.size()==3);
+        assert(fvec.size()==54);
+        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];
@@ -607,9 +608,9 @@ struct chwirut2_functor {
                 double x = _x[i];
                 double factor = 1./(b[1]+b[2]*x);
                 double e = exp(-b[0]*x);
-                fjac[i+ldfjac*0] = -x*e*factor;
-                fjac[i+ldfjac*1] = -e*factor*factor;
-                fjac[i+ldfjac*2] = -x*e*factor*factor;
+                fjac(i,0) = -x*e*factor;
+                fjac(i,1) = -e*factor*factor;
+                fjac(i,2) = -x*e*factor*factor;
             }
         }
         return 0;
@@ -666,15 +667,16 @@ void testNistChwirut2(void)
 
 
 struct misra1a_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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(m==14);
-        assert(n==2);
-        assert(ldfjac==14);
+        assert(b.size()==2);
+        assert(fvec.size()==14);
+        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] ;
@@ -682,8 +684,8 @@ struct misra1a_functor {
         }
         else { // compute fjac at b
             for(i=0; i<14; i++) {
-                fjac[i+ldfjac*0] = (1.-exp(-b[1]*x[i]));
-                fjac[i+ldfjac*1] = (b[0]*x[i]*exp(-b[1]*x[i]));
+                fjac(i,0) = (1.-exp(-b[1]*x[i]));
+                fjac(i,1) = (b[0]*x[i]*exp(-b[1]*x[i]));
             }
         }
         return 0;
@@ -736,16 +738,18 @@ void testNistMisra1a(void)
 }
 
 struct hahn1_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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 int called=0; printf("call hahn1_functor with  iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
-        assert(m==236);
-        assert(n==7);
-        assert(ldfjac==236);
+
+        assert(b.size()==7);
+        assert(fvec.size()==236);
+        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;
@@ -756,14 +760,14 @@ struct hahn1_functor {
             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+ldfjac*0] = 1.*fact;
-                fjac[i+ldfjac*1] = x*fact;
-                fjac[i+ldfjac*2] = xx*fact;
-                fjac[i+ldfjac*3] = xxx*fact;
+                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+ldfjac*4] = x*fact;
-                fjac[i+ldfjac*5] = xx*fact;
-                fjac[i+ldfjac*6] = xxx*fact;
+                fjac(i,4) = x*fact;
+                fjac(i,5) = xx*fact;
+                fjac(i,6) = xxx*fact;
             }
         }
         return 0;
@@ -827,15 +831,16 @@ void testNistHahn1(void)
 }
 
 struct misra1d_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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(m==14);
-        assert(n==2);
-        assert(ldfjac==14);
+        assert(b.size()==2);
+        assert(fvec.size()==14);
+        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];
@@ -844,8 +849,8 @@ struct misra1d_functor {
         else { // compute fjac at b
             for(i=0; i<14; i++) {
                 double den = 1.+b[1]*x[i];
-                fjac[i+ldfjac*0] = b[1]*x[i] / den;
-                fjac[i+ldfjac*1] = b[0]*x[i]*(den-b[1]*x[i])/den/den;
+                fjac(i,0) = b[1]*x[i] / den;
+                fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den;
             }
         }
         return 0;
@@ -899,15 +904,16 @@ void testNistMisra1d(void)
 
 
 struct lanczos1_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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(m==24);
-        assert(n==6);
-        assert(ldfjac==24);
+        assert(b.size()==6);
+        assert(fvec.size()==24);
+        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];
@@ -915,12 +921,12 @@ struct lanczos1_functor {
         }
         else { // compute fjac at b
             for(i=0; i<24; i++) {
-                fjac[i+ldfjac*0] = exp(-b[1]*x[i]);
-                fjac[i+ldfjac*1] = -b[0]*x[i]*exp(-b[1]*x[i]);
-                fjac[i+ldfjac*2] = exp(-b[3]*x[i]);
-                fjac[i+ldfjac*3] = -b[2]*x[i]*exp(-b[3]*x[i]);
-                fjac[i+ldfjac*4] = exp(-b[5]*x[i]);
-                fjac[i+ldfjac*5] = -b[4]*x[i]*exp(-b[5]*x[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;
@@ -982,15 +988,16 @@ void testNistLanczos1(void)
 }
 
 struct rat42_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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(m==9);
-        assert(n==3);
-        assert(ldfjac==9);
+        assert(b.size()==3);
+        assert(fvec.size()==9);
+        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];
@@ -999,9 +1006,9 @@ struct rat42_functor {
         else { // compute fjac at b
             for(i=0; i<9; i++) {
                 double e = exp(b[1]-b[2]*x[i]);
-                fjac[i+ldfjac*0] = 1./(1.+e);
-                fjac[i+ldfjac*1] = -b[0]*e/(1.+e)/(1.+e);
-                fjac[i+ldfjac*2] = +b[0]*e*x[i]/(1.+e)/(1.+e);
+                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;
@@ -1056,15 +1063,16 @@ void testNistRat42(void)
 }
 
 struct MGH10_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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(m==16);
-        assert(n==3);
-        assert(ldfjac==16);
+        assert(b.size()==3);
+        assert(fvec.size()==16);
+        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];
@@ -1074,9 +1082,9 @@ struct MGH10_functor {
             for(i=0; i<16; i++) {
                 double factor = 1./(x[i]+b[2]);
                 double e = exp(b[1]*factor);
-                fjac[i+ldfjac*0] = e;
-                fjac[i+ldfjac*1] = b[0]*factor*e;
-                fjac[i+ldfjac*2] = -b[1]*b[0]*factor*factor*e;
+                fjac(i,0) = e;
+                fjac(i,1) = b[0]*factor*e;
+                fjac(i,2) = -b[1]*b[0]*factor*factor*e;
             }
         }
         return 0;
@@ -1132,15 +1140,16 @@ void testNistMGH10(void)
 
 
 struct BoxBOD_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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(m==6);
-        assert(n==2);
-        assert(ldfjac==6);
+        assert(b.size()==2);
+        assert(fvec.size()==6);
+        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];
@@ -1149,8 +1158,8 @@ struct BoxBOD_functor {
         else { // compute fjac at b
             for(i=0; i<6; i++) {
                 double e = exp(-b[1]*x[i]);
-                fjac[i+ldfjac*0] = 1.-e;
-                fjac[i+ldfjac*1] = b[0]*x[i]*e;
+                fjac(i,0) = 1.-e;
+                fjac(i,1) = b[0]*x[i]*e;
             }
         }
         return 0;
@@ -1205,15 +1214,16 @@ void testNistBoxBOD(void)
 }
 
 struct MGH17_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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(m==33);
-        assert(n==5);
-        assert(ldfjac==33);
+        assert(b.size()==5);
+        assert(fvec.size()==33);
+        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];
@@ -1221,11 +1231,11 @@ struct MGH17_functor {
         }
         else { // compute fjac at b
             for(i=0; i<33; i++) {
-                fjac[i+ldfjac*0] = 1.;
-                fjac[i+ldfjac*1] = exp(-b[3]*x[i]);
-                fjac[i+ldfjac*2] = exp(-b[4]*x[i]);
-                fjac[i+ldfjac*3] = -x[i]*b[1]*exp(-b[3]*x[i]);
-                fjac[i+ldfjac*4] = -x[i]*b[2]*exp(-b[4]*x[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;
@@ -1288,15 +1298,16 @@ void testNistMGH17(void)
 }
 
 struct MGH09_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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(m==11);
-        assert(n==4);
-        assert(ldfjac==11);
+        assert(b.size()==4);
+        assert(fvec.size()==11);
+        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;
@@ -1307,10 +1318,10 @@ struct MGH09_functor {
             for(i=0; i<11; i++) {
                 double x = _x[i], xx=x*x;
                 double factor = 1./(xx+x*b[2]+b[3]);
-                fjac[i+ldfjac*0] = (xx+x*b[1]) * factor;
-                fjac[i+ldfjac*1] = b[0]*x* factor;
-                fjac[i+ldfjac*2] = - b[0]*(xx+x*b[1]) * x * factor * factor;
-                fjac[i+ldfjac*3] = - b[0]*(xx+x*b[1]) * factor * factor;
+                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;
@@ -1371,16 +1382,17 @@ void testNistMGH09(void)
 
 
 struct Bennett5_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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(m==154);
-        assert(n==3);
-        assert(ldfjac==154);
+        assert(b.size()==3);
+        assert(fvec.size()==154);
+        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];
@@ -1389,9 +1401,9 @@ struct Bennett5_functor {
         else { // compute fjac at b
             for(i=0; i<154; i++) {
                 double e = pow(b[1]+x[i],-1./b[2]);
-                fjac[i+ldfjac*0] = e;
-                fjac[i+ldfjac*1] = - b[0]*e/b[2]/(b[1]+x[i]);
-                fjac[i+ldfjac*2] = b[0]*e*log(b[1]+x[i])/b[2]/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;
@@ -1446,16 +1458,17 @@ void testNistBennett5(void)
 }
 
 struct thurber_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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++;
-        assert(m==37);
-        assert(n==7);
-        assert(ldfjac==37);
+        assert(b.size()==7);
+        assert(fvec.size()==37);
+        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;
@@ -1466,14 +1479,14 @@ struct thurber_functor {
             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+ldfjac*0] = 1.*fact;
-                fjac[i+ldfjac*1] = x*fact;
-                fjac[i+ldfjac*2] = xx*fact;
-                fjac[i+ldfjac*3] = xxx*fact;
+                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+ldfjac*4] = x*fact;
-                fjac[i+ldfjac*5] = xx*fact;
-                fjac[i+ldfjac*6] = xxx*fact;
+                fjac(i,4) = x*fact;
+                fjac(i,5) = xx*fact;
+                fjac(i,6) = xxx*fact;
             }
         }
         return 0;
@@ -1538,15 +1551,16 @@ void testNistThurber(void)
 }
 
 struct rat43_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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(m==15);
-        assert(n==4);
-        assert(ldfjac==15);
+        assert(b.size()==4);
+        assert(fvec.size()==15);
+        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];
@@ -1556,10 +1570,10 @@ struct rat43_functor {
             for(i=0; i<15; i++) {
                 double e = exp(b[1]-b[2]*x[i]);
                 double power = -1./b[3];
-                fjac[i+ldfjac*0] = pow(1.+e, power);
-                fjac[i+ldfjac*1] = power*b[0]*e*pow(1.+e, power-1.);
-                fjac[i+ldfjac*2] = -power*b[0]*e*x[i]*pow(1.+e, power-1.);
-                fjac[i+ldfjac*3] = b[0]*power*power*log(1.+e)*pow(1.+e, power);
+                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;
@@ -1620,16 +1634,17 @@ void testNistRat43(void)
 
 
 struct eckerle4_functor {
-    static int f(int m, int n, const double *b, double *fvec, double *fjac, int ldfjac, int iflag)
+    static int f(const VectorXd &b, VectorXd &fvec, MatrixXd &fjac, int iflag)
     {
         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(m==35);
-        assert(n==3);
-        assert(ldfjac==35);
+        assert(b.size()==3);
+        assert(fvec.size()==35);
+        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];
@@ -1639,9 +1654,9 @@ struct eckerle4_functor {
             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+ldfjac*0] = e / b[1];
-                fjac[i+ldfjac*1] = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12;
-                fjac[i+ldfjac*2] = (x[i]-b[2])*e*b[0]/b[1]/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;