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++

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

 

Example. Compute arctanh(0.11) using GSL library function gsl_atanh. Here C++ language is used.
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;
}
Compiling and execution of the program
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.


For more  information on the GSL functions see
http://www.gnu.org/software/gsl/manual/html_node/Numerical-Differentiation-functions.html
http://www.gnu.org/software/gsl/manual/html_node/Numerical-Differentiation-Examples.html

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;
   }
}
Test fun. Only part of the output is listed here.
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

http://www.lrz-muenchen.de/services/software/mathematik/gsl/fortran/index.html