Fortran: FGSL

<< zur Fortran-Startseite
< BLAS und ATLAS	LAPACK >

Wikipedia hat einen Artikel zum Thema:

Gnu Scientific Library

Allgemeines

Bei der "GNU Scientific Library" (GSL) handelt es sich um eine in C geschriebene Bibliothek. Diese bietet Funktionen für ein weites Spektrum der Mathematik. Beispielhaft seien folgende Bereiche genannt:

Komplexe Zahlen
Lineare Algebra
Polynome
Statistik
Fast Fourier Transformation (FFT)
Numerische Differentiation
Numerische Integration
Gewöhnliche Differentialgleichungen
IEEE Floating-Point-Arithmetik

FGSL ist ein Fortran-Interface für diese GSL-Bibliothek. FGSL selbst deckt nicht die komplette Funktionspalette von GSL ab, inzwischen sind aber auch lineare Algebra und die FFT-Funktionen Bestandteil von FGSL, auch wenn zu diesem Zweck der Einsatz der optimierten Bibliotheken LAPACK und FFTW empfohlen wird. Für einige Teilbereiche werden nicht alle in C implementierten Datentypen unterstützt.

FGSL wurde unter Verwendung einiger Fortran 2003-Sprachmerkmale erstellt. Die Einbindung von FGSL in eigene Programme setzt aus diesem Grunde einen entsprechenden Fortran-Compiler voraus. Der g95-Compiler erfüllt z.B. diese Voraussetzungen.

Derzeit ist die FGSL-Version 0.9.4 vom 31. Mai 2011 aktuell.

Beispiele

Beispiel: Datentypen, Potenzierung und mathematische Konstanten

C Fortran

#include <stdio.h>
#include <gsl/gsl_math.h>
     
int main (void)
{
   double a;
   double d = 5.0;

   a = gsl_pow_2 (d) * M_PI_4;

   printf ("Kreisflaeche = %f\n", a);
   return 0;
}

program bsp
  use fgsl
  implicit none
  
  real( kind = fgsl_double ) :: a
  real( kind = fgsl_double ) :: d = 5.0_fgsl_double

  a = d ** 2 * m_pi_4

  write( *, * ) "Kreisflaeche = ", a
end program bsp

Kompilieren, Linken:

gcc bsp.c -I/usr/local/include -L/usr/local/lib -lgsl -lgslcblas -lm

Kompilieren, Linken:

g95 bsp.f03 -I/usr/local/include/g95 -L/usr/local/lib -lgsl -lfgsl_g95 -lgslcblas

Der kind-Wert fgsl_double entspricht dem c_double aus dem iso_c_binding-Modul. FGSL kennt die speziellen pow-Funktionen aus GSL nicht, da Fortran ohnehin über einen eigenen Potenzierungsoperator verfügt. Neben der m_pi_4-Konstante ( = $\pi /4$ ) kennt FGSL noch eine ganze Reihe anderer mathematischer Konstanten, z.B.:

`m_e`	...	e, Eulersche Zahl, 2,714...
`m_euler`	...	Eulersche Konstante, 0,577...
`m_pi`	...	$\pi$
`m_pi_2`	...	$\pi /2$
`m_sqrt2`	...	${\sqrt {2}}$

Auch jede Menge physikalische Konstanten kennt FGSL, z.B.:

`fgsl_const_mksa_speed_of_light`	...	Lichtgeschwindigkeit im Vakuum, 2.9979... x 10⁸ m / s
`fgsl_const_mksa_molar_gas`	...	Allgemeine Gaskonstante, 8.314472 J / (K · mol)
`fgsl_const_mksa_inch`	...	0.0254 m
`fgsl_const_mksa_torr`	...	133.322... Pa
`fgsl_const_num_giga`	...	10⁹

Beispiel: Lösen einer quadratischen Gleichung

Gesucht ist die Lösung der quadratischen Gleichung $x^{2}+12x+37=0$

C

Fortran

#include <stdio.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_poly.h>

int main (void)
{
  double a, b, c;
  gsl_complex z1, z2;
  int i;
   
  a = 1.0;
  b = 12.0;
  c = 37.0;
  
  i = gsl_poly_complex_solve_quadratic(a, b, c,
                                       &z1, &z2);

  if(i == 1)
  {
    /* nur eine Lsg.*/
    printf("z = (%f, %f)\n", z1.dat[0], 
                             z1.dat[1]);   
  }
  else    
  {
    /* 2 Lsg., reell oder komplex */
    printf( "z1 = (%f, %f)\n", z1.dat[0], 
                               z1.dat[1]); 
    printf( "z2 = (%f, %f)\n", z2.dat[0], 
                               z2.dat[1]); 
  }
     
  return 0;
  
/* Ausgabe:
    z1 = (-6.000000, -1.000000)
    z2 = (-6.000000, 1.000000) 
*/  
}

program bsp
  use fgsl
  implicit none
  
  real( kind = fgsl_double )                 :: a, b, c
  complex( fgsl_double ), dimension( 2 )     :: z
  integer( kind = fgsl_int )                 :: i
  
  a = 1.0_fgsl_double
  b = 12.0_fgsl_double
  c = 37.0_fgsl_double
  
  i = fgsl_poly_complex_solve_quadratic( a, b, c, z(1), z(2) )

  if( i == 1 ) then
    ! nur eine Lsg.
    write( *, * ) "   z =", z(1)   
  else    
    ! 2 Lsg., reell oder komplex
    write( *, * ) "   z1,2 =", z   
  end if  

! Ausgabe:
!   z1,2 = (-6.,-1.) (-6.,1.)
end program bsp

Beispiel: Numerische Integration

Gesucht ist die Lösung des Integrals $I=\int _{0}^{1}{\frac {\sin x}{x}}dx$

C

Fortran

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
     

double f(double x, void *params) 
{  
  return(sin(x) / x);
}
     
int main ()
{
  double result, error;
  size_t neval;
  gsl_function func;

  func.function = &f;
  func.params = 0;
     
  gsl_integration_qng (&func, 0.0, 1.0, 1e-9, 1e-9,
                       &result, &error, &neval); 
     
  printf ("Ergebnis = %f\n", result);
     
  return 0;
  
/* Ausgabe:
     Ergebnis = 0.946083
*/
}

module integral
  implicit none

  contains
    function f( x, params ) bind(c)
      use, intrinsic :: iso_c_binding
      implicit none
  
      real( kind = c_double )        :: f
      real( kind = c_double ), value :: x
      type( c_ptr ), value           :: params
 
      f = sin( x ) / x
    end function f 
end module integral


program bsp
  use fgsl
  use integral
  use, intrinsic :: iso_c_binding
  implicit none
  
  
  real( kind = fgsl_double )   :: result, error
  integer( kind = fgsl_size_t) :: neval
  integer( kind = fgsl_int)    :: i
  type( fgsl_function )        :: func


  func = fgsl_function_init( f, c_null_ptr )
     
  i = fgsl_integration_qng ( func,              &
                         0.0_fgsl_double,       &
                         1.0_fgsl_double,       &
                         1e-9_fgsl_double,      & 
                         1e-9_fgsl_double,      &   
                         result, error, neval ) 
     
  write( *, * ) "Ergebnis =", result
  
! Ausgabe:
!   Ergebnis = 0.946083070367183
end program bsp

(F)GSL stellt verschiedene Möglichkeiten der numerischen Integration zur Verfügung. Hier wurde die Funktion für den QNG-Algorithmus (non-adaptive Gauss-Kronrod) gewählt.

Beispiel: IEEE-Floating-Point-Arithmetik

Darstellung einer Fließkommazahl nach IEEE 754-Standard:

(-1)^{s}(1.ffff...)2^{E}\

Beispielsweise wird eine 32-bit-Fließkommazahl binär so aufgegliedert:

seeeeeeeefffffffffffffffffffffff

s ... sign, 1 bit
e, E ... exponent, 8 bit (2⁸ = 256; E_min = -127; E_max = 128)
f ... fraction, 23 bit

Näheres zum IEEE 754-Standard findet sich z.B. bei IEEE 754

FGSL bietet Subroutinen, um Fießkommazahlen anschaulich entsprechend dem IEEE 754-Standard auszugeben:

fgsl_ieee_printf( x ) ... Ausgabe der Zahl x im IEEE-Format auf stdout
fgsl_ieee_fprintf( str, x ) ... Ausgabe der Zahl x im IEEE-Format. str ist ein C-Zeiger (C: FILE *, Fortran: type( c_ptr )).

C

Fortran

#include <gsl/gsl_ieee_utils.h>

int main()
{
  float  a = 1.5;
  double b = -2.6666666666666666;

  
  gsl_ieee_printf_float( &a );
  gsl_ieee_printf_double( &b );
}

program bsp
  use fgsl
  implicit none
  
  real( kind = fgsl_float )  :: a = 1.5
  real( kind = fgsl_double ) :: b = -2.6666666666666666d0
  
  call fgsl_ieee_printf( a )
  call fgsl_ieee_printf( b )
end program bsp

Ausgabe:

 1.10000000000000000000000*2^0
-1.0101010101010101010101010101010101010101010101010101*2^1

Mit der Subroutine fgsl_ieee_env_setup() lassen sich über die Environment-Variable GSL_IEEE_MODE einige nützliche Attribute (Rundungsmodus etc.) festlegen.

C

Fortran

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

int main()
{
  float a = 1.5;
  int i;

  gsl_ieee_env_setup();

  for(i = 0; i < 5; i++)
  {
    a /= 0.11;
    printf("%f\n", a);
  } 
}

program bsp
  use fgsl
  implicit none
  
  real    :: a = 1.5
  integer :: i
  
  call fgsl_ieee_env_setup()

  do i = 0, 4
    a = a / 0.11
    write( *, * ) a 
  end do  
end program bsp

Programmaufruf:

GSL_IEEE_MODE="round-to-nearest" ./a.out

Ausgabe:

GSL_IEEE_MODE="round-to-nearest,trap-common"
13.636364
123.966942
1126.972168
10245.201172
93138.195312

Programmaufruf:

GSL_IEEE_MODE="round-to-nearest" ./a.out

Ausgabe:

GSL_IEEE_MODE="round-to-nearest,trap-common"
13.636364
123.96695
1126.9723
10245.203
93138.21

Programmaufruf:

GSL_IEEE_MODE="round-up" ./a.out

Ausgabe:

GSL_IEEE_MODE="round-up,trap-common"
13.636364
123.966949
1126.972290
10245.203125
93138.210938

Programmaufruf:

GSL_IEEE_MODE="round-up" ./a.out

Ausgabe:

GSL_IEEE_MODE="round-up,trap-common"
13.636364
123.96695
1126.9723
10245.203
93138.21

Programmaufruf:

GSL_IEEE_MODE="round-down" ./a.out

Ausgabe:

GSL_IEEE_MODE="round-down,trap-common"
13.636363
123.966934
1126.972046
10245.200195
93138.179688

Programmaufruf:

GSL_IEEE_MODE="round-down" ./a.out

Ausgabe:

GSL_IEEE_MODE="round-down,trap-common"
13.636363
123.966934
1126.972
10245.2
93138.18

Weblinks

<< zur Fortran-Startseite
< BLAS und ATLAS	LAPACK >