Physics 6720: Using GSL with C++ or C


The Gnu Scientific Library (GSL) is a public domain package with a wide variety of codes for doing scientific calculations, such as Fourier transforms, wavelet transforms, root finding, and least squares fitting among others.

Documentation for the library is available from the Gnu software site. Local documentation is available in the Gnu info system.


Example: Linear Regression


Suppose you want to fit a straight line to a set of ndata points xi, yi, with errors ei. First, look at the example code. This code fits a set of built-in data. The key routine is gsl_fit_wlinear. Look it up in the function index. You should find the prototype
int gsl_fit_wlinear (const double * x, const size_t xstride, const double * w, const size_t wstride, const double * y, const size_t ystride, size_t n, double * c0, double * c1, double * cov00, double * cov01, double * cov11, double * chisq)
and the explanation of the parameters. The stride parameter is usually 1. Setting it to 2 causes the fitter to use every other value in the list, beginning with the first (??).

Armed with this information we write our code as follows:

#include <iostream>
#include <gsl/gsl_fit.h>
#define MAX 100

int main(){

   int i, n;
   double x[MAX], y[MAX], e[MAX];
   double c0, c1, cov00, cov01, cov11, chisq;

   // Read the values of the matrix
   std::cout << "Enter n \n";
   std::cin >> n;
   std::cout << "On each line type x y sigma:\n";
   for(i = 0; i < n; i++){
     for(i = 0; i < n; i++)std::cin >> x[i] >> y[i] >> e[i];
     if(std::cin.fail())
     {
        std::cout << "Failed to read data\n";
	return 1;
     }
   }

   // Fit the straight line
   gsl_fit_wlinear (x, 1, w, 1, y, 1, n, 
                    &c0, &c1, &cov00, &cov01, &cov11, 
                    &chisq);

   // Write the answer
   std::cout << "Fit to y = m*x + b gives\n";
   std::cout << "m = " << c1 << " +/- " << sqrt(cov11) 
             << "b = " << c0 << " +/- " << sqrt)cov00) << "\n";
   std::cout << "chisq = " << chisq << " for " << n-2 << "df\n";
   return 0;
}
Note that we are not required to use the same names as given in the prototype.
Physics Department Home Page

This page is maintained by:

Carleton DeTar Mail Form Last modified 25 October 2007