Description: expicitly compute Jacobians of fdfsolver structs
 as it is no more available as a J member in GSL 2.x
Author: Cédric Boutillier <boutil@debian.org>
Origin: vendor
Bug-Debian: https://bugs.debian.org/804499
Bug: https://github.com/SciRuby/rb-gsl/issues/24
Forwarded: https://github.com/SciRuby/rb-gsl/pull/31
Last-Update: 2016-02-29

--- a/ext/gsl_native/histogram.c
+++ b/ext/gsl_native/histogram.c
@@ -1169,6 +1169,7 @@
   size_t n, dof;      /* # of data points */
   size_t p = 3;  /* # of fitting parameters */
   gsl_multifit_function_fdf f;
+  gsl_matrix *J = NULL;
   gsl_matrix *covar = NULL;
   gsl_vector *x = NULL;
   double sigma, mean, height, errs, errm, errh, chi2;
@@ -1197,6 +1198,7 @@
   hh.binend = binend;
   n = binend - binstart + 1;
 
+  J = gsl_matrix_alloc(n, p);
   covar = gsl_matrix_alloc(p, p);
 
   f.f = Gaussian_f;
@@ -1219,7 +1221,8 @@
   sigma = sqrt(gsl_vector_get(s->x, 0));
   mean = gsl_vector_get(s->x, 1);
   height = gsl_vector_get(s->x, 2)*sigma*sqrt(2*M_PI);
-  gsl_multifit_covar(s->J, 0.0, covar);
+  gsl_multifit_fdfsolver_jac(s, J);
+  gsl_multifit_covar(J, 0.0, covar);
   chi2 = gsl_pow_2(gsl_blas_dnrm2(s->f));   /* not reduced chi-square */
   dof = n - p;
   errs = sqrt(chi2/dof*gsl_matrix_get(covar, 0, 0))/sigma/2;
@@ -1305,6 +1308,7 @@
   size_t n, dof;      /* # of data points */
   size_t p = 2;  /* # of fitting parameters */
   gsl_multifit_function_fdf f;
+  gsl_matrix *J = NULL;
   gsl_matrix *covar = NULL;
   gsl_vector *x = NULL;
   double sigma, height, errs, errh, chi2;
@@ -1332,6 +1336,7 @@
   hh.binend = binend;
   n = binend - binstart + 1;
 
+  J = gsl_matrix_alloc(n, p);
   covar = gsl_matrix_alloc(p, p);
 
   f.f = Rayleigh_f;
@@ -1353,7 +1358,8 @@
   } while (status == GSL_CONTINUE);
   sigma = sqrt(gsl_vector_get(s->x, 0));
   height = gsl_vector_get(s->x, 1)*sigma*sigma;
-  gsl_multifit_covar(s->J, 0.0, covar);
+  gsl_multifit_fdfsolver_jac(s, J);
+  gsl_multifit_covar(J, 0.0, covar);
   chi2 = gsl_pow_2(gsl_blas_dnrm2(s->f));   /* not reduced chi-square */
   dof = n - p;
   errs = sqrt(chi2/dof*gsl_matrix_get(covar, 0, 0))/sigma/2;
@@ -1441,6 +1447,7 @@
   size_t n, dof;      /* # of data points */
   size_t p = 2;  /* # of fitting parameters */
   gsl_multifit_function_fdf f;
+  gsl_matrix *J = NULL;
   gsl_matrix *covar = NULL;
   gsl_vector *x = NULL;
   double b, height, errs, errh, chi2;
@@ -1468,6 +1475,7 @@
   hh.binend = binend;
   n = binend - binstart + 1;
 
+  J = gsl_matrix_alloc(n, p);
   covar = gsl_matrix_alloc(p, p);
 
   f.f = xExponential_f;
@@ -1489,7 +1497,8 @@
   } while (status == GSL_CONTINUE);
   b = gsl_vector_get(s->x, 0);
   height = gsl_vector_get(s->x, 1);
-  gsl_multifit_covar(s->J, 0.0, covar);
+  gsl_multifit_fdfsolver_jac(s, J);
+  gsl_multifit_covar(J, 0.0, covar);
   chi2 = gsl_pow_2(gsl_blas_dnrm2(s->f));   /* not reduced chi-square */
   dof = n - p;
   errs = sqrt(chi2/dof*gsl_matrix_get(covar, 0, 0));
--- a/ext/gsl_native/multifit.c
+++ b/ext/gsl_native/multifit.c
@@ -324,6 +324,7 @@
 {
   gsl_multifit_fdfsolver *solver = NULL;
   gsl_vector *g = NULL;
+  gsl_matrix *J = NULL;
   int status;
   double epsabs;
   Data_Get_Struct(obj, gsl_multifit_fdfsolver, solver);
@@ -331,7 +332,9 @@
   case 1:
     Need_Float(argv[0]);
     g = gsl_vector_alloc(solver->x->size);
-    gsl_multifit_gradient(solver->J, solver->f, g);
+    J = gsl_matrix_alloc(solver->f->size, solver->x->size);
+    gsl_multifit_fdfsolver_jac(solver, J);
+    gsl_multifit_gradient(J, solver->f, g);
     epsabs = NUM2DBL(argv[0]);
     status = gsl_multifit_test_gradient(g, epsabs);
     gsl_vector_free(g);
@@ -353,15 +356,17 @@
 {
   gsl_multifit_fdfsolver *solver = NULL;
   gsl_vector *g = NULL;
+  gsl_matrix *J = gsl_matrix_alloc(solver->f->size, solver->x->size);
   // local variable "status" declared and set, but never used
   //int status;
   Data_Get_Struct(obj, gsl_multifit_fdfsolver, solver);
+  gsl_multifit_fdfsolver_jac(solver, J);
   if (argc == 1) {
     Data_Get_Vector(argv[0], g);
-    return INT2FIX(gsl_multifit_gradient(solver->J, solver->f, g));
+    return INT2FIX(gsl_multifit_gradient(J, solver->f, g));
   } else {
     g = gsl_vector_alloc(solver->x->size);
-    /*status =*/ gsl_multifit_gradient(solver->J, solver->f, g);
+    /*status =*/ gsl_multifit_gradient(J, solver->f, g);
     return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, g);
   }
 }
@@ -377,15 +382,17 @@
   Need_Float(argv[0]);
   Data_Get_Struct(obj, gsl_multifit_fdfsolver, solver);
   epsrel = NUM2DBL(argv[0]);
+  gsl_matrix *J = gsl_matrix_alloc(solver->f->size, solver->x->size);
+  gsl_multifit_fdfsolver_jac(solver, J);
   switch (argc) {
   case 1:
     covar = gsl_matrix_alloc(solver->x->size, solver->x->size);
-    /*status =*/ gsl_multifit_covar(solver->J, epsrel, covar);
+    /*status =*/ gsl_multifit_covar(J, epsrel, covar);
     return Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, covar);
     break;
   case 2:
     Data_Get_Matrix(argv[1], covar);
-    return INT2FIX(gsl_multifit_covar(solver->J, epsrel, covar));
+    return INT2FIX(gsl_multifit_covar(J, epsrel, covar));
     break;
   default:
     rb_raise(rb_eArgError, "wrong number of arguments");
@@ -418,7 +425,9 @@
 {
   gsl_multifit_fdfsolver *solver = NULL;
   Data_Get_Struct(obj, gsl_multifit_fdfsolver, solver);
-  return Data_Wrap_Struct(cgsl_matrix_view_ro, 0, NULL, solver->J);
+  gsl_matrix *J = gsl_matrix_alloc(solver->f->size, solver->x->size);
+  gsl_multifit_fdfsolver_jac(solver, J);
+  return Data_Wrap_Struct(cgsl_matrix_view_ro, 0, NULL, J);
 }
 
 /* singleton */
@@ -1699,7 +1708,9 @@
   covar = gsl_matrix_alloc(p, p);
   chi2 = gsl_pow_2(gsl_blas_dnrm2(solver->f));   /* not reduced chi-square */
   dof = n - p;
-  gsl_multifit_covar(solver->J, 0.0, covar);
+  gsl_matrix *J = gsl_matrix_alloc(n,p);
+  gsl_multifit_fdfsolver_jac(solver, J);
+  gsl_multifit_covar(J, 0.0, covar);
   for (i = 0; i < p; i++)
     gsl_vector_set(verr, i, sqrt(chi2/dof*gsl_matrix_get(covar, i, i)));
   gsl_matrix_free(covar);