From a9d24f91057441bbd2e3ed9e7536b071121526cb Mon Sep 17 00:00:00 2001
From: "D. V. Wiebe" <dvw@ketiltrout.net>
Date: Thu, 10 Mar 2016 14:09:26 -0800
Subject: [PATCH] GSL-2.x support.

---
 src/plugins/fits/non_linear.h          | 67 +++++++++++++++++++++-------------
 src/plugins/fits/non_linear_weighted.h | 66 ++++++++++++++++++++-------------
 2 files changed, 82 insertions(+), 51 deletions(-)

diff --git a/src/plugins/fits/non_linear.h b/src/plugins/fits/non_linear.h
index 4506704..74e82e7 100644
--- a/src/plugins/fits/non_linear.h
+++ b/src/plugins/fits/non_linear.h
@@ -18,6 +18,7 @@
 #include <gsl/gsl_blas.h>
 #include <gsl/gsl_multifit_nlin.h>
 #include <gsl/gsl_statistics.h>
+#include <gsl/gsl_version.h>
 #include "common.h"
 
 struct data {
@@ -100,6 +101,7 @@ bool kstfit_nonlinear(
   gsl_multifit_function_fdf	function;
   gsl_vector_view vectorViewInitial;
   gsl_matrix* pMatrixCovariance;
+  gsl_matrix *pMatrixJacobian;
   struct data d;  
   double dXInitial[NUM_PARAMS];
   double* pInputX;
@@ -177,37 +179,50 @@ bool kstfit_nonlinear(
             }
             iIterations++;
           } while( iStatus == GSL_CONTINUE && iIterations < MAX_NUM_ITERATIONS );
-          gsl_multifit_covar( pSolver->J, 0.0, pMatrixCovariance );
-
-          //
-          // determine the fitted values...
-          //
-          for( i=0; i<NUM_PARAMS; i++ ) {
-            dXInitial[i] = gsl_vector_get( pSolver->x, i );
-          }
-
-          for( i=0; i<iLength; i++ ) {
-            vectorOutYFitted->value()[i] = function_calculate( pInputX[i], dXInitial );
-            vectorOutYResiduals->value()[i] = pInputY[i] - vectorOutYFitted->value()[i];
-          }
+#if GSL_MAJOR_VERSION >= 2
+          pMatrixJacobian = gsl_matrix_alloc( iLength, NUM_PARAMS );
+#else
+          pMatrixJacobian = pSolver->J;
+#endif
+          if ( pMatrixJacobian != NULL) {
+#if GSL_MAJOR_VERSION >= 2
+            gsl_multifit_fdfsolver_jac( pSolver, pMatrixJacobian );
+#endif
+            gsl_multifit_covar( pMatrixJacobian, 0.0, pMatrixCovariance );
+
+            //
+            // determine the fitted values...
+            //
+            for( i=0; i<NUM_PARAMS; i++ ) {
+              dXInitial[i] = gsl_vector_get( pSolver->x, i );
+            }
 
-          //
-          // fill in the parameter values and covariance matrix...
-          //
-          for( i=0; i<NUM_PARAMS; i++ ) {
-            vectorOutYParameters->value()[i] = gsl_vector_get( pSolver->x, i );
-            for( j=0; j<NUM_PARAMS; j++ ) {
-              vectorOutYCovariance->value()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
+            for( i=0; i<iLength; i++ ) {
+              vectorOutYFitted->value()[i] = function_calculate( pInputX[i], dXInitial );
+              vectorOutYResiduals->value()[i] = pInputY[i] - vectorOutYFitted->value()[i];
             }
-          }
 
-          //
-          // determine the value of chi^2/nu
-          //
-          scalarOutChi->setValue(gsl_blas_dnrm2( pSolver->f ));
+            //
+            // fill in the parameter values and covariance matrix...
+            //
+            for( i=0; i<NUM_PARAMS; i++ ) {
+              vectorOutYParameters->value()[i] = gsl_vector_get( pSolver->x, i );
+              for( j=0; j<NUM_PARAMS; j++ ) {
+                vectorOutYCovariance->value()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
+              }
+            }
 
-          bReturn = true;
+            //
+            // determine the value of chi^2/nu
+            //
+            scalarOutChi->setValue(gsl_blas_dnrm2( pSolver->f ));
 
+            bReturn = true;
+            
+#if GSL_MAJOR_VERSION >= 2
+            gsl_matrix_free( pMatrixJacobian );
+#endif
+          }
           gsl_matrix_free( pMatrixCovariance );
         }
         gsl_multifit_fdfsolver_free( pSolver );
diff --git a/src/plugins/fits/non_linear_weighted.h b/src/plugins/fits/non_linear_weighted.h
index 6ca7d6f..347ae9d 100644
--- a/src/plugins/fits/non_linear_weighted.h
+++ b/src/plugins/fits/non_linear_weighted.h
@@ -18,6 +18,7 @@
 #include <gsl/gsl_blas.h>
 #include <gsl/gsl_multifit_nlin.h>
 #include <gsl/gsl_statistics.h>
+#include <gsl/gsl_version.h>
 #include "common.h"
 
 struct data {
@@ -101,6 +102,7 @@ bool kstfit_nonlinear_weighted(
   gsl_multifit_function_fdf	function;
   gsl_vector_view vectorViewInitial;
   gsl_matrix* pMatrixCovariance;
+  gsl_matrix *pMatrixJacobian;
   struct data d;  
   double dXInitial[NUM_PARAMS];
   double* pInputs[3];
@@ -193,37 +195,51 @@ bool kstfit_nonlinear_weighted(
           }
           while( iStatus == GSL_CONTINUE && iIterations < MAX_NUM_ITERATIONS );
 
-          gsl_multifit_covar( pSolver->J, 0.0, pMatrixCovariance );
-
-          //
-          // determine the fitted values...
-          //
-          for( i=0; i<NUM_PARAMS; i++ ) {
-            dXInitial[i] = gsl_vector_get( pSolver->x, i );
-          }
+#if GSL_MAJOR_VERSION >= 2
+          pMatrixJacobian = gsl_matrix_alloc( iLength, NUM_PARAMS );
+#else
+          pMatrixJacobian = pSolver->J;
+#endif
+
+          if ( pMatrixJacobian != NULL) {
+#if GSL_MAJOR_VERSION >= 2
+            gsl_multifit_fdfsolver_jac( pSolver, pMatrixJacobian );
+#endif
+            gsl_multifit_covar( pMatrixJacobian, 0.0, pMatrixCovariance );
+
+            //
+            // determine the fitted values...
+            //
+            for( i=0; i<NUM_PARAMS; i++ ) {
+              dXInitial[i] = gsl_vector_get( pSolver->x, i );
+            }
 
-          for( i=0; i<iLength; i++ ) {
-            vectorOutYFitted->value()[i] = function_calculate( pInputs[XVALUES][i], dXInitial );
-            vectorOutYResiduals->value()[i] = pInputs[YVALUES][i] - vectorOutYFitted->value()[i];
-          }
+            for( i=0; i<iLength; i++ ) {
+              vectorOutYFitted->value()[i] = function_calculate( pInputs[XVALUES][i], dXInitial );
+              vectorOutYResiduals->value()[i] = pInputs[YVALUES][i] - vectorOutYFitted->value()[i];
+            }
 
-          //
-          // fill in the parameter values and covariance matrix...
-          //
-          for( i=0; i<NUM_PARAMS; i++ ) {
-            vectorOutYParameters->value()[i] = gsl_vector_get( pSolver->x, i );
-            for( j=0; j<NUM_PARAMS; j++ ) {
-              vectorOutYCovariance->value()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
+            //
+            // fill in the parameter values and covariance matrix...
+            //
+            for( i=0; i<NUM_PARAMS; i++ ) {
+              vectorOutYParameters->value()[i] = gsl_vector_get( pSolver->x, i );
+              for( j=0; j<NUM_PARAMS; j++ ) {
+                vectorOutYCovariance->value()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
+              }
             }
-          }
 
-          //
-          // determine the value of chi^2/nu
-          //
-          scalarOutChi->setValue(gsl_blas_dnrm2( pSolver->f ));
+            //
+            // determine the value of chi^2/nu
+            //
+            scalarOutChi->setValue(gsl_blas_dnrm2( pSolver->f ));
 
-          bReturn = true;
+            bReturn = true;
 
+#if GSL_MAJOR_VERSION >= 2
+            gsl_matrix_free( pMatrixJacobian );
+#endif
+          }
           gsl_matrix_free( pMatrixCovariance );
         }
         gsl_multifit_fdfsolver_free( pSolver );