DISCLAIMER!: I was asked to write this page a good while ago and it is possible that some details are out of date. Most of this page should still be useful.
numerical libraries FFTW, LAPACK, GSLThis page describes some of the freely available alternatives to using commercial C/Fortran compilers and numerical libraries. An increasing number of libraries are being written in C, but it is usually still possible to access them from Fortran (for GSL examples see the last section 'Calling C from Fortran'). The advantages of using free software include:
Linux is a free complete Unix-like operating system developed under the GNU General Public License. Most free software, especially that related to scientific research, is developed under Linux or Unix. Installing Linux on a home computer makes use of free software easier and provides a better software development environment. A few of the most widely available distributions are RedHat, SuSE, and Debian. The drawback is that administration of a Linux PC requires some computing experience.
Cygwin is an alternative to installing Linux. It's a Unix-like environment for Windows. Many GNU packages, including the GNU compilers and X-programs, are available under Cygwin. The full installation is large and in practice requires a fast network connection, but the installation procedure is relatively straight forward. The setup.exe may be used at any time to add/delete parts of the installation. The first time you open a cygwin terminal, type whoami and pwd. If you are not in a directory /home/username then create it, restart, and try pwd again. If you are still not in your home directory, run mkpasswd -l > /etc/passwd and mkgroup -l > /etc/group . Restart and type startx to get an X-terminal if the appropriate components are installed.
Many free numerical libraries can be found at the Netlib repository. They are generally of very high quality and well tested. GSL is a large collection of numerical numerical functions covering pure, applied and statistical functions. GSL is written in C, but related projects include wrappers for other languages including Perl and Python. Most functions can be called from Fortran through a hand-written wrapper, see below. Other widely used libraries include FFTW, a collection of very fast Fourier transform C routines (with built in wrappers for Fortran) and LAPACK, a linear algebra package written in Fortran 77 (versions are available for other languages).
The GNU Compiler Collection (GCC), contains front ends for C, C++, Objective-C, Fortran, Java and Ada. It is by far the most popular non-commercial compiler and has been ported to many platforms.
The GNU C/C++ compiler is part of the standard Linux installation and is also available under Cygwin.
The
GNU Fortran compiler
g77 is a program that calls the gcc compiler
with options to recognise programs written in Fortran. It supports ANSI
FORTRAN 77 conformance, plus popular extensions to Fortran including some
ANSI/ISO Fortran 90 features. The g95 project
is a full F95 compiler available for linux and cygwin.
The Intel Fortran compiler offers complete F90 support and generates fast binaries for Intel chips under Linux. The free unsupported version may be used for non-commercial purposes.
The C library GSL has been used for the following discussion, but the principles apply to any numerical C-Fortran integration.
Calling GSL from Fortran is not always straight forward and requires a good working knowledge of both Fortran and C. Despite difficulties outlined below, it may be desirable to integrate GSL into a large existing code, say. Most GSL functions can be accessed through a wrapper -- a routine written in C that is designed to be callable from Fortran. The wrapper accounts for the differences between C and Fortran, where the Fortran complier applies the following rules:
g77 and ifort compilers convert
subroutine names to lower case.
g77 compiler appends an underscore to
symbols, or two if the name contained an underscore. So
the Fortran subroutine WRAP_SUB appears as
wrap_sub__ in the compiled object file.
ifort compiler always appends only one underscore,
WRAP_SUB becomes wrap_sub_.
INTEGER N becomes in C
int* n.
The value N identifies with *n.
REAL A(10) becomes
float* a.
Note that the C array will be indexed
such that it starts at 0, here A(1) is identified
with a[0].
INTEGER A(3,2) are
A(1,1), A(2,1), A(3,1), A(1,2), A(2,2), A(3,2).
C arrays are stored in row major order, the elements of
int A[1][2] are
A[0][0], A[0][1], A[0][2], A[1][0], A[1][1], A[1][2].
void* .
char*, with the length of the string
appended to the end of the argument list as a
long int (note no *).
call wrap_sub(x,2)
g77 corresponds to the C function-prototype:
void wrap_sub__(double*,int*);
ifort as the Fortran complier, only one underscore
needs to be appended to the wrapper declaration in the C code.
Conversions for the numeric data types that appear in the declaration
may vary according to the compiler or compiler options. The standard
defaults are
float REAL 4 byte
double DOUBLE PRECISION 8
int INTEGER 4
long int INTEGER*8 8
float complex COMPLEX 4+4 (iso c99)
double complex DOUBLE COMPLEX 8+8
Fortran does not have equivalents for all the numeric data types (e.g. unsigned) but with care they can be cast to or from an available data type in the wrapper (see Example 2). The Fortran compilers deal with functions in different ways, but most compilers pass the address of subroutines. For example, the Fortran subroutine:
subroutine sub(x)
double precision x
EXTERNAL-type argument:
call use_sub(sub)
use_sub is:
void use_sub__( void(*s)(double*) );
sub through the function
pointer (also see Example 3):
double x;
...
s(&x);
Simple C structs can usually be passed by sending each item within the struct separately, but some GSL functions use pointers to types which are undefined in the documentation, however. Also GSL constants may be required for which the values are unknown without searching through the header files. In summary, GSL uses many C features which do not extend naturally to a wrapper routine for Fortran. Unfortunately, several calls to GSL functions are often required to achieve one logical goal. In this case it would require a lot of effort to write wrappers for each GSL call, and its more efficient to group these calls in one C-function with a simple function prototype (see Example 2).
In the following example the Fortran code gsl.f calls the GSL Bessel special function routine via the C wrapper in gslwr.c . The object files are linked via the Fortran compiler.
/*---------------------------------------- gslwr.c --*/
#include <gsl/gsl_sf_bessel.h>
void wrapper_j0__(double* y, double* x){
*y = gsl_sf_bessel_J0(*x);
}
c ---------------------------------------- gsl.f ------
double precision x, y
x = 1d0
call wrapper_j0(y,x)
print*, y
-------------------------------------------------------
> g77 -c gsl.f
> gcc -c gslwr.c
> g77 gsl.o gslwr.o -lgsl -lgslcblas
-------------------------------------------------------
In the next example all the variables actually passed between C and
Fortran are standard numeric types. Numeric items in C structs can
be passed one by one. The Fortran program does not
usually need to know about gsl_ types that are
undefined in the documentation.
The example is just one way to
get n samples from a Gaussian distribution, using the seed s.
Once the GSL random number generator is initalised by a call to
rng_initialise, the sampling function
rng_sample_gaussian can be called repeatedly.
/*---------------------------------------- gslwr.c --*/
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
static gsl_rng* r; /* r is (almost) never used explicitly and does
not need to be seen by the Fortran program. */
void rng_initialise__(int* s){
r = gsl_rng_alloc(gsl_rng_taus); /* constant gsl_rng_taus is unknown to Fortran */
gsl_rng_set(r, (unsigned long int)(*s)); /* s is cast from int to (unsigned long int) */
}
void rng_sample_gaussian__(double* x, int* n, double* sigma){
int i;
for(i=0; i<*n; i++)
x[i] = gsl_ran_gaussian(r,*sigma);
}
c ---------------------------------------- gsl.f ------
integer n, seed
parameter (n=50)
double precision x
dimension x(n)
...
call system_clock(seed)
call rng_initialise(seed)
call rng_sample_gaussian(x,n,1d0)
-------------------------------------------------------
EXTERNAL arguments:
Differentiation
In this example, the Fortran subroutine pwr which
sets y:=x^1.5, is passed to wrapper_diff_central__.
As functions are complied differently by different Fortran compliers,
pwr is written as a subroutine rather than a function.
It is then assigned to a gsl_function type via the
sub2fn routine.
The gsl_function data type is defined by
typedef struct{
double (*function)(double* x, void* params);
void* params;
} gsl_function;
gsl_diff_central evaluates the
derivative at x=2. Any subroutine of the form sub(y,x)
which sets y as a function of x could be passed to the wrapper.
/*---------------------------------------- gslwr.c --*/
#include <gsl/gsl_diff.h>
double sub2fn(double x, void* p){
void (*subp)(double*,double*) = p;
double y;
subp(&y,&x);
return y;
}
void wrapper_diff_central__(
void* sub, double* x, double* result, double* abserr)
{
gsl_function F;
F.function = &sub2fn;
F.params = sub;
gsl_diff_central(&F, *x, result, abserr);
}
c ---------------------------------------- gsl.f ------
subroutine pwr(y,x)
double precision y,x
y = x**1.5d0
end
c - - - - - - - - - - - - - - - - - - - - - - - - - - -
double precision result, abserr
external pwr
call wrapper_diff_central(pwr,2d0,result,abserr)
-------------------------------------------------------