Fortran: FGSL

---

Vorlage:Wikipedia

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 108 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	...	109

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=\underset{0}{\overset{1}{\int }}\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 Vorlage:W

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

• GSL - GNU Scientific Library
• FGSL: A Fortran interface to the GNU Scientific Library

---