Gnu Scientific Library GSL

Gnu Scientific Library (GSL) is C library of numerical routines for scientific computing ( http://www.gnu.org/software/gsl/  ).

The manual of the library is at www-pages  http://www.gnu.org/software/gsl/manual/html_node/

The development of the GSL library was started 1996 and the library has become more and more comprehensive each year. At the moment, it includes routines for instance for computing values of special functions, one and multidimensional integration, fast fourier transforms, and least squares fitting.

The GSL library is available in the server mutteri (and ruuvi).  Its functions may be linked to C, C++, and Fortran programs.

GSL and C and C++

• The functions of the GSL library are used according to the instructions given at the library www-page.
• The header file of GSL function to be used is included in the program.
• Furthermore, when compiling the program libraries gsl and gslblas need to be linked to the program.

Each group of GSL functions has its own header file. For instance, special functions like erf, Bessel, spherical harmonics need the header file  gsl_sf.h.  The names of the header files may be found at the www-pages. Read carefully the manual to find the proper header file.  

Example. Compute arctanh(0.11) using GSL library function gsl_atanh. Here C language is used.

#include <stdio.h>
#include <gsl/gsl_math.h>

int main(void) {

   double a = 0.11;

   printf("%lf \n", gsl_atanh(a));

   return 0;
}

Linking

gcc atanh.c -lgsl -lgslcblas 

Execution

heavy> ./a.out
0.110447

ruuvi.it.helsinki.fi> cat koe.cpp
#include <iostream>
#include <gsl/gsl_math.h>

using namespace std;

int main(void) {
   double a = 0.11, tmp;

   tmp = gsl_atanh(a);
   cout << "atanh(0.11) " << tmp << " \n";
   return 0;

}

ruuvi.it.helsinki.fi> c++ koe.cpp -lgsl -lgslcblas
ruuvi.it.helsinki.fi> ./a.out
atanh(0.11) 0.110447

Example. This program computes the numerical derivative of Lorenz function. The main program forms the x-scale and calls GSL routine gsl_deriv_central. This function calls a user written function lorenz, which computes the values of this function. The name of this user written  function is passed to the GSL routine in a structure 'F' of type 'gsl_function' as a field 'F.function'.

Here a  part of the description of the function gls_deriv_central is listed:

— Function: int gsl_deriv_central (const gsl_function * f, double x, double h, double * result, double * abserr)

This function computes the numerical derivative of the function f at the point x using an adaptive central difference algorithm with a step-size of h. The derivative is returned in result and an estimate of its absolute error is returned in abserr.

The initial value of h is used to estimate an optimal step-size, based on the scaling of the truncation error and round-off error in the derivative calculation. The derivative is computed using a 5-point rule for equally spaced abscissae at x-h, x-h/2, x, x+h/2, x+h, with an error estimate taken from the difference between the 5-point rule and the corresponding 3-point rule x-h, x, x+h. Note that the value of the function at x does not contribute to the derivative calculation, so only 4-points are actually used.

ruuvi.it.helsinki.fi> cat gsl-der.c
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>

double lorenz(double x, void * params);
void xscale(double *x);

#define n 101

int main (void) {
   double x[n];
   int i;
   gsl_function F;
   double result, abserr;

/* form x scale */
   xscale(x);

   F.function = &lorenz;
   F.params = 0;
   printf ("Derivative of Lorenz function\n");

   for (i = 0; i < n; i++) {
      gsl_deriv_central (&F, x[i], 1e-8, &result, &abserr);
      printf ("x = %f, f'(x) = %.10f +/- %.10f\n", x[i], result, abserr);
   }

   return 0;
}

double lorenz(double x, void * params){
/* Lorenz function */

    return 1/(1 + pow(x,2));

}

void xscale(double *x) {
/* Compute x-scale */
   int i;
   double xa = -5.0, xstep = 0.1;

   for (i = 0; i < n; i++) {
       x[i] = xa + i*xstep;
   }
}

ruuvi.it.helsinki.fi> pico gsl-der.c
ruuvi.it.helsinki.fi> gcc gsl-der.c -lgsl -lgslcblas
ruuvi.it.helsinki.fi> ./a.out
Lorenz function
...
x = -0.600000, f'(x) = 0.6487889189 +/- 0.0000001214
x = -0.500000, f'(x) = 0.6399999956 +/- 0.0000001211
x = -0.400000, f'(x) = 0.5945303268 +/- 0.0000001275
x = -0.300000, f'(x) = 0.5050080058 +/- 0.0000001404
x = -0.200000, f'(x) = 0.3698224949 +/- 0.0000001297
x = -0.100000, f'(x) = 0.1960592170 +/- 0.0000001471
x = 0.000000, f'(x) = 0.0000000000 +/- 0.0000001332
x = 0.100000, f'(x) = -0.1960592170 +/- 0.0000001471
x = 0.200000, f'(x) = -0.3698224986 +/- 0.0000001445
x = 0.300000, f'(x) = -0.5050080021 +/- 0.0000001256
x = 0.400000, f'(x) = -0.5945303120 +/- 0.0000001275
x = 0.500000, f'(x) = -0.6399999956 +/- 0.0000001211
x = 0.600000, f'(x) = -0.6487889189 +/- 0.0000001214
x = 0.700000, f'(x) = -0.6306022241 +/- 0.0000001066
x = 0.800000, f'(x) = -0.5948839958 +/- 0.0000000992

GSL and Fortran

GSL library may also be used in a Fortran program. One needs to write a 'gateway' function using C language, which calls the GSl function and includes the proper GSL header file in the program. The gateway function is declared 'external' in the Fortran program.

Of cource, all aspects related to calling a C function from a Fortran program need to be taken into account, for instance the matching of C and Fortran data types.

The linking phase is also more complicated. The gateway function is compiled separately (with option -c). The object file is linked to the Fortran main program as well as the GSl libraries gsl and glscblas.

Example.  Compute arctanh(0.11) using GSL library function gsl_atanh. Here the main program was written in Fortran.

ruuvi.it.helsinki.fi> cat gsltest.f90
program gsltest
implicit none

 real(kind=selected_real_kind(12)) :: a = 0.11, res

 external :: glsgateway

 call gslgateway(a,res)

 write(*,*) 'x', a, 'atanh(x)', res

end program gsltest

C-function

#include <gsl/gsl_math.h>

void gslgateway_(double *x, double *res){

   *res = gsl_atanh(*x);

}

Compiling and linking

gcc -c gslgateway.c
gfortran gsltest.f90 gslgateway.o -lgsl -lgslcblas


ruuvi.it.helsinki.fi> ./a.out
 x 0.109999999403954 atanh(x) 0.110446915186750

See also www-pages:

http://www.maths.qmul.ac.uk/~dhruba/tips_and_tricks/node9.html

http://www.ladhyx.polytechnique.fr/people/willis/aco/freesoftware.html#CfromF