Fortran Best Practices
======================

.. highlight:: fortran

This page collects a modern canonical way of doing things in Fortran. It is meant to be short, and it is assumed that you already know how to program in other languages (like Python, C/C++, ...) and also know Fortran syntax a bit. Some things in Fortran are obsolete, so this guide only shows the "one correct/canonical modern way" how to do things.

Summary of the language: http://www.cs.umbc.edu/~squire/fortranclass/summary.shtml

Language features are explained at: http://en.wikipedia.org/wiki/Fortran_language_features

The best resource is a recent Fortran standard, for example the Fortran 2003 standard: http://www.j3-fortran.org/doc/year/04/04-007.pdf

Fortran Style Guide
-------------------

Naming Convention
~~~~~~~~~~~~~~~~~

Ultimately this is a matter of preference. Here is a style guide that we like
and that seems to be prevalent in most scientific codes (as well as the Fortran
standard library) and you are welcome to follow it.

1. Use lowercase for all Fortran constructs (``do``, ``subroutine``, ``module``,
   ...).
2. Follow short mathematical notation for mathematical variables/functions
   (``Ylm``, ``Gamma``, ``gamma``, ``Enl``, ``Rnl``, ...).
3. For other names use all lowercase: try to keep names to one or two
   syllables; if more are required, use underscores to clarify (``sortpair``,
   ``whitechar``, ``meshexp``, ``numstrings``, ``linspace``, ``meshgrid``,
   ``argsort``, ``spline``, ``spline_interp``, ``spline_interpolate``,
   ``stoperr``, ``stop_error``, ``meshexp_der``).

For example "spline interpolation" can be shortened to
``spline_interpolation``, ``spline_interpolate``, ``spline_interp``,
``spline``, but not to ``splineint`` ("int" could mean integration, integer,
etc. --- too much ambiguity, even in the clear context of a computational
code). This is in contrast to ``get_argument()`` where ``getarg()`` is
perfectly clean and clear.

The above are general guidelines.  In general, choosing the right name
certainly depends on the word being truncated as to whether the first syllable
is sufficient. Usually it is but clearly not always. Thus some thought should
go into step "try to keep names to 2 syllables or less" since it can really
affect the indicativeness and simplicity. Simple consistent naming rules are a
real help in this regard -- for both collaboration and for one's own sanity
when going back to some old code you haven't seen in while.

Indentation
~~~~~~~~~~~

Use 4 spaces indentation (this is again a matter of preference as
some people prefer 2 or 3 spaces).

Comparison to Other Languages
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

On the other hand, in most of the rest of the programming world, where the main
focus is, in one form or another, on defining and using large sets of complex
objects, with tons of properties and behaviors, known only in the code in which
they are defined (as opposed to defined by the same  notation throughout the
literature), it makes more sense to use longer, more descriptive names. The
naming conventions one sees used in more general-purpose languages such as C++
and Python, therefore, are perfectly consistent with their more general-purpose
missions. But Fortran has a different mission (numerical scientific computing).

.. _floating_point_numbers:

Floating Point Numbers
----------------------

Somewhere create and export a parameter ``dp``::

    integer, parameter:: dp=kind(0.d0)                   ! double precision

and declare floats as::

    real(dp) :: a, b, c

Always write all floating point constants with the ``_dp`` suffix::

    1.0_dp, 3.5_dp, 1.34e8_dp

and never any other way (see also the gotcha
:ref:`floating_point_numbers_gotcha`).
Omitting the dot in the literal constant is also incorrect.

To print floating point double precision
numbers without losing precision, use the ``(es23.16)`` format (see
http://stackoverflow.com/questions/6118231/why-do-i-need-17-significant-digits-and-not-16-to-represent-a-double/).

It is safe to assign integers to floating point numbers without losing
accuracy::

    real(dp) :: a
    a = 3

In order to impose floating point division (as opposed to integer division
``1/2`` equal to ``0``), one can convert the integer to a floating point number
by::

    real(dp) :: a
    a = real(1, dp) / 2       ! 'a' is equal to 0.5_dp

Integer Division
----------------

Just like in Python 2.x or C, when doing things like ``(N-1)/N`` where ``N`` is an integer and you want a floating point division, force the compiler to use floats at the right hand side, for example by::

    (N - 1.0_dp)/N

As long as one of the ``/`` operands is a float, a floating point division will be used.

.. _modules:

Modules and Programs
--------------------

Only use modules and programs. Always setup a module in the following way::

    module ode1d
    use types, only: dp, pi
    use utils, only: stop_error
    implicit none
    private
    public integrate, normalize, parsefunction, get_val, rk4step, eulerstep, &
            rk4step2, get_midpoints, rk4_integrate, rk4_integrate_inward, &
            rk4_integrate_inward2, rk4_integrate3, rk4_integrate4, &
            rk4_integrate_inward4

    contains

    subroutine get_val(...)
    ...
    end subroutine
    ...

    end module

The ``implicit none`` statement works for the whole module (so you don't need to worry about it). By keeping the ``private`` empty, all your subroutines/data types will be private to the module by default. Then you export things by putting it into the ``public`` clause.

Setup programs in the following way::

    program uranium
    use fmesh, only: mesh_exp
    use utils, only: stop_error, dp
    use dft, only: atom
    implicit none

    integer, parameter :: Z = 92
    real(dp), parameter :: r_min = 8e-9_dp, r_max = 50.0_dp, a = 1e7_dp
    ...
    print *, "I am running"
    end program

Notice the "explicit imports" (using Python terminology) in the ``use`` statements. You can also use "implicit imports" like::

    use fmesh

But just like in Python, this should be avoided ("explicit is better than implicit") in most cases.

Arrays
------

When passing arrays in and out of a subroutine/function, use
the following pattern for 1D arrays (it is called `assumed-shape`)::

    subroutine f(r)
    real(dp), intent(out) :: r(:)
    integer :: n, i
    n = size(r)
    do i = 1, n
        r(i) = 1.0_dp / i**2
    enddo
    end subroutine

2D arrays::

    subroutine g(A)
    real(dp), intent(in) :: A(:, :)
    ...
    end subroutine

and call it like this::

    real(dp) :: r(5)
    call f(r)

No array copying is done above. It has the following
advantages:

* the shape and size of the array is passed in automatically
* the shape is checked at compile time, the size optionally at runtime
* allows to use strides and all kinds of array
  arithmetic without actually copying any data.

This should always be your default
way of passing arrays in and out of subroutines. However
in the following cases one can (or has to) use `explicit-shape` arrays:

* returning an array from a function
* interfacing with C code or legacy Fortran (like Lapack)
* operating on arbitrary shape array with the given function (however there are
  also other ways to do that, see :ref:`elemental` for more information)

To use `explicit-shape` arrays, do::

    subroutine f(n, r)
    integer, intent(in) :: n
    real(dp), intent(out) :: r(n)
    integer :: i
    do i = 1, n
        r(i) = 1.0_dp / i**2
    enddo
    end subroutine

2D arrays::

    subroutine g(m, n, A)
    integer, intent(in) :: m, n
    real(dp), intent(in) :: A(m, n)
    ...
    end subroutine

and call it like this::

    real(dp) :: r(5)
    call f(size(r), r)

In order to return an array from a function, do::

    function f(n) result(r)
    integer, intent(in) :: n
    real(dp) :: r(n)
    integer :: i
    do i = 1, n
        r(i) = 1.0_dp / i**2
    enddo
    end function

If you want to enforce/check the size of the arrays, put at the beginning of
the function::

    if (size(r) /= 4) stop "Incorrect size of 'r'"

To initialize an array, do::

    integer :: r(5)
    r = [1, 2, 3, 4, 5]

This syntax is valid since the Fortran 2003 standard and it is the preferred
syntax (the old syntax ``r = (/ 1, 2, 3, 4, 5 /)`` should only be used if you
cannot use Fortran 2003).

In order for the array to start with different index than 1, do::

    subroutine print_eigenvalues(kappa_min, lam)
    integer, intent(in) :: kappa_min
    real(dp), intent(in) :: lam(kappa_min:)

    integer :: kappa
    do kappa = kappa_min, ubound(lam, 1)
        print *, kappa, lam(kappa)
    end do
    end subroutine

Multidimensional Arrays
-----------------------

Always access slices as ``V(:, 1)``, ``V(:, 2)``, or ``V(:, :, 1)``, e.g. the colons should be on the left. That way the stride is contiguous and it will be fast. So when you need some slice in your algorithm, always setup the array in a way, so that you call it as above. If you put the colon on the right, it will be slow.

Example::

    dydx = matmul(C(:, :, i), y) ! fast
    dydx = matmul(C(i, :, :), y) ! slow

In other words, the "fortran storage order" is: smallest/fastest changing/innermost-loop index first, largest/slowest/outermost-loop index last ("Inner-most are left-most."). So the elements of a 3D array ``A(N1,N2,N3)`` are stored, and thus most efficiently accessed, as::

    do i3 = 1, N3
        do i2 = 1, N2
            do i1 = 1, N1
                A(i1, i2, i3)
            end do
        end do
    end do

Associated array of vectors would then be most efficiently accessed as::

    do i3 = 1, N3
        do i2 = 1, N2
            A(:, i2, i3)
        end do
    end do

And associated set of matrices would be most efficiently accessed as::

    do i3 = 1, N3
        A(:, :, i3)
    end do

Storing/accessing as above then accesses always contiguous blocks of memory, directly adjacent to one another; no skips/strides.

When not sure, always rewrite (in your head) the algorithm to use strides, for example the first loop would become::

    do i3 = 1, N3
        Ai3 = A(:, :, i3)
        do i2 = 1, N2
            Ai2i3 = Ai3(:, i2)
            do i1 = 1, N1
                Ai2i3(i1)
            end do
        end do
    end do

the second loop would become::

    do i3 = 1, N3
        Ai3 = A(:, :, i3)
        do i2 = 1, N2
            Ai3(:, i2)
        end do
    end do

And then make sure that all the strides are always on the left. Then it will be fast.

.. _elemental:

Element-wise Operations on Arrays Using Subroutines/Functions
-------------------------------------------------------------

There are three approaches:

* ``elemental`` subroutines
* `explicit-shape` arrays
* implementing the operation for vectors and write simple wrapper subroutines
  (that use ``reshape`` internally) for each array shape

In the first approach,
one uses the ``elemental`` keyword to create a function like this::

    real(dp) elemental function nroot(n, x) result(y)
    integer, intent(in) :: n
    real(dp), intent(in) :: x
    y = x**(1._dp / n)
    end function

All arguments (in and out) must be scalars. You can then use
this function with arrays of any (compatible) shape, for example::

    print *, nroot(2, 9._dp)
    print *, nroot(2, [1._dp, 4._dp, 9._dp, 10._dp])
    print *, nroot(2, reshape([1._dp, 4._dp, 9._dp, 10._dp], [2, 2]))
    print *, nroot([2, 3, 4, 5], [1._dp, 4._dp, 9._dp, 10._dp])
    print *, nroot([2, 3, 4, 5], 4._dp)

The output will be::

   3.0000000000000000
   1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795
   1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795
   1.0000000000000000        1.5874010519681994        1.7320508075688772        1.5848931924611136
   2.0000000000000000        1.5874010519681994        1.4142135623730951        1.3195079107728942

In the above, typically ``n`` is a parameter and ``x`` is the array of an
arbitrary shape, but as you can see, Fortran does not care as long as the final
operation makes sense (if one argument is an array, then the other arguments
must be either arrays of the same shape or scalars). If it does not, you will
get a compiler error.

The ``elemental`` keyword implies the ``pure`` keyword, so the subroutine must
be pure (can only use ``pure`` subroutines and have no side effects).

If the elemental function's algorithm can be made faster using array operations
inside, or if for some reason the arguments must be arrays of incompatible
shapes,
then one should use the other two approaches.
One can make ``nroot`` operate
on a vector and write a simple wrappers for other array shapes::

    function nroot(n, x) result(y)
    integer, intent(in) :: n
    real(dp), intent(in) :: x(:)
    real(dp) :: y(size(x))
    y = x**(1._dp / n)
    end function

    function nroot_0d(n, x) result(y)
    integer, intent(in) :: n
    real(dp), intent(in) :: x
    real(dp) :: y
    real(dp) :: tmp(1)
    tmp = nroot(n, [x])
    y = tmp(1)
    end function

    function nroot_2d(n, x) result(y)
    integer, intent(in) :: n
    real(dp), intent(in) :: x(:, :)
    real(dp) :: y(size(x, 1), size(x, 2))
    y = reshape(nroot(n, reshape(x, [size(x)])), [size(x, 1), size(x, 2)])
    end function

And use as follows::

    print *, nroot_0d(2, 9._dp)
    print *, nroot(2, [1._dp, 4._dp, 9._dp, 10._dp])
    print *, nroot_2d(2, reshape([1._dp, 4._dp, 9._dp, 10._dp], [2, 2]))

This will print::

    3.0000000000000000
    1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795
    1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795


Or one can use `explicit-shape` arrays as follows::

    function nroot(n, k, x) result(y)
    integer, intent(in) :: n, k
    real(dp), intent(in) :: x(k)
    real(dp) :: y(k)
    y = x**(1._dp / n)
    end function

Use as follows::

    print *, nroot(2, 1, [9._dp])
    print *, nroot(2, 4, [1._dp, 4._dp, 9._dp, 10._dp])
    print *, nroot(2, 4, reshape([1._dp, 4._dp, 9._dp, 10._dp], [2, 2]))

The output is the same as before::

       3.0000000000000000
       1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795
       1.0000000000000000        2.0000000000000000        3.0000000000000000        3.1622776601683795


Allocatable Arrays
------------------

When using allocatable arrays (as opposed to pointers), Fortran manages the
memory automatically and it is not possible to create memory leaks.

For example you can allocate it inside a subroutine::

    subroutine foo(lam)
    real(dp), allocatable, intent(out) :: lam(:)
    allocate(lam(5))
    end subroutine

And use somewhere else::

    real(dp), allocatable :: lam(:)
    call foo(lam)

When the ``lam`` symbol goes out of scope, Fortran will deallocate it. If
``allocate`` is called twice on the same array, Fortran will abort with a
runtime error. One can check if ``lam`` is already allocated and deallocate it
if needed (before another allocation)::

    if (allocated(lam)) deallocate(lam)
    allocate(lam(10))


File Input/Output
-----------------

To read from a file::

    integer :: u
    open(newunit=u, file="log.txt", status="old")
    read(u, *) a, b
    close(u)

Write to a file::

    integer :: u
    open(newunit=u, file="log.txt", status="replace")
    write(u, *) a, b
    close(u)

To append to an existing file::

    integer :: u
    open(newunit=u, file="log.txt", position="append", status="old")
    write(u, *) N, V(N)
    close(u)

The ``newunit`` keyword argument to ``open`` is a Fortran 2008 standard, in older compilers, just replace
``open(newunit=u, ...)`` by::

    open(newunit(u), ...)

where the ``newunit`` function is defined by::

    integer function newunit(unit) result(n)
    ! returns lowest i/o unit number not in use
    integer, intent(out), optional :: unit
    logical inuse
    integer, parameter :: nmin=10   ! avoid lower numbers which are sometimes reserved
    integer, parameter :: nmax=999  ! may be system-dependent
    do n = nmin, nmax
        inquire(unit=n, opened=inuse)
        if (.not. inuse) then
            if (present(unit)) unit=n
            return
        end if
    end do
    call stop_error("newunit ERROR: available unit not found.")
    end function

.. _c_interface:

Interfacing with C
------------------

Write a C wrapper using the ``iso_c_binding`` module::

    module fmesh_wrapper

    use iso_c_binding, only: c_double, c_int
    use fmesh, only: mesh_exp

    implicit none

    contains

    subroutine c_mesh_exp(r_min, r_max, a, N, mesh) bind(c)
    real(c_double), intent(in) :: r_min
    real(c_double), intent(in) :: r_max
    real(c_double), intent(in) :: a
    integer(c_int), intent(in) :: N
    real(c_double), intent(out) :: mesh(N)
    call mesh_exp(r_min, r_max, a, N, mesh)
    end subroutine

    ! wrap more functions here
    ! ...

    end module

You need to declare the length of all arrays (``mesh(N)``) and pass it as a
parameter. The Fortran compiler will check that the C and Fortran types match.
If it compiles, you can then trust it, and call it from C using the following
declaration:

.. code-block:: c

    void c_mesh_exp(double *r_min, double *r_max, double *a, int *N,
            double *mesh);

use it as:

.. code-block:: c

    int N=5;
    double r_min, r_max, a, mesh[N];
    c_mesh_exp(&r_min, &r_max, &a, &N, mesh);

No matter if you are passing arrays in or out, always allocate them in C first, and you are (in C) responsible for the memory management. Use Fortran to fill (or use) your arrays (that you own in C).

If calling the Fortran ``exp_mesh`` subroutine from the ``c_exp_mesh`` subroutine is a problem (CPU efficiency), you can simply implement whatever the routine does directly in the ``c_exp_mesh`` subroutine. In other words, use the ``iso_c_binding`` module as a direct way to call Fortran code from C, and you can make it as fast as needed.

.. _python_interface:

Interfacing with Python
-----------------------

Using Cython
~~~~~~~~~~~~

To wrap Fortran code in Python, export it to C first (see above) and then write
this Cython code:

.. code-block:: cython

    from numpy cimport ndarray
    from numpy import empty

    cdef extern:
        void c_mesh_exp(double *r_min, double *r_max, double *a, int *N,
                double *mesh)

    def mesh_exp(double r_min, double r_max, double a, int N):
        cdef ndarray[double, mode="c"] mesh = empty(N, dtype=double)
        c_mesh_exp(&r_min, &r_max, &a, &N, &mesh[0])
        return mesh

The memory is allocated and owned (reference counted) by Python, and a pointer is given to the Fortran code. Use this approach for both "in" and "out" arrays.

Notice that we didn't write any C code --- we only told fortran to use the C
calling convention when producing the ".o" files, and then we pretended in
Cython, that the function is implemented in C, but in fact, it is linked in
from Fortran directly. So this is the most direct way of calling Fortran from
Python. There is no intermediate step, and no unnecessary processing/wrapping
involved.


Using ctypes
~~~~~~~~~~~~

Alternatively, you can assign C-callable names to your Fortran
routines like this::

    subroutine mesh_exp(r_min, r_max, a, N, mesh) bind(c, name='mesh_exp')
      real(c_double), intent(in), value :: r_min
      real(c_double), intent(in), value :: r_max
      real(c_double), intent(in), value :: a
      integer(c_int), intent(in), value :: N
      real(c_double), intent(out) :: mesh(N)

      ! ...

    end subroutine mesh_exp

and use the builtin `ctypes <http://docs.python.org/library/ctypes.html>`_
Python package to dynamically load
shared object files containing your C-callable Fortran routines and
call them directly:

.. code-block:: python

    from ctypes import CDLL, POINTER, c_int, c_double
    from numpy import empty

    fortran = CDLL('./libmyfortranroutines.so')

    mesh = empty(N, dtype="double")
    fortran.mesh_exp(c_double(r_min), c_double(r_max), c_double(a), c_int(N),
                     mesh.ctypes.data_as(POINTER(c_double)))



Callbacks
---------

There are two ways to implement callbacks to be used like this::

    subroutine foo(a, k)
    use integrals, only: simpson
    real(dp) :: a, k
    print *, simpson(f, 0._dp, pi)
    print *, simpson(f, 0._dp, 2*pi)

    contains

    real(dp) function f(x) result(y)
    real(dp), intent(in) :: x
    y = a*sin(k*x)
    end function f

    end subroutine foo

The traditional approach is to simply declare the ``f`` dummy variable as a
subroutine/function using::

    module integrals
    use types, only: dp
    implicit none
    private
    public simpson

    contains

    real(dp) function simpson(f, a, b) result(s)
    real(dp), intent(in) :: a, b
    interface
        real(dp) function f(x)
        use types, only: dp
        implicit none
        real(dp), intent(in) :: x
        end function
    end interface
    s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b))
    end function

    end module

The other approach since f2003 is to first define a new type for our callback,
and then use ``procedure(func)`` as the type of the dummy argument::

    module integrals
    use types, only: dp
    implicit none
    private
    public simpson

    contains

    real(dp) function simpson(f, a, b) result(s)
    real(dp), intent(in) :: a, b
    interface
        real(dp) function func(x)
        use types, only: dp
        implicit none
        real(dp), intent(in) :: x
        end function
    end interface
    procedure(func) :: f
    s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b))
    end function

    end module

The new type can also be defined outside of the function (and reused), like::

    module integrals
    use types, only: dp
    implicit none
    private
    public simpson

    interface
        real(dp) function func(x)
        use types, only: dp
        implicit none
        real(dp), intent(in) :: x
        end function
    end interface

    contains

    real(dp) function simpson(f, a, b) result(s)
    real(dp), intent(in) :: a, b
    procedure(func) :: f
    s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b))
    end function

    real(dp) function simpson2(f, a, b) result(s)
    real(dp), intent(in) :: a, b
    procedure(func) :: f
    real(dp) :: mid
    mid = (a + b)/2
    s = simpson(f, a, mid) + simpson(f, mid, b)
    end function

    end module

.. _callbacks:

Type Casting in Callbacks
-------------------------

There are essentially five different ways to do that, each
with its own advantages and disadvantages.

The methods I, II and V can be used both in C and Fortran.
The methods III and IV are only available in Fortran.
The method VI is obsolete and should not be used.

I: Work Arrays
~~~~~~~~~~~~~~~

Pass a "work array" or two which are packed with everything needed by the
caller and unpacked by the called routine. This is the old way -- e.g., how
LAPACK does it.

Integrator::

    module integrals
    use types, only: dp
    implicit none
    private
    public simpson

    contains

    real(dp) function simpson(f, a, b, data) result(s)
    real(dp), intent(in) :: a, b
    interface
        real(dp) function func(x, data)
        use types, only: dp
        implicit none
        real(dp), intent(in) :: x
        real(dp), intent(inout) :: data(:)
        end function
    end interface
    procedure(func) :: f
    real(dp), intent(inout) :: data(:)
    s = (b-a) / 6 * (f(a, data) + 4*f((a+b)/2, data) + f(b, data))
    end function

    end module

Usage::

    module test
    use types, only: dp
    use integrals, only: simpson
    implicit none
    private
    public foo

    contains

    real(dp) function f(x, data) result(y)
    real(dp), intent(in) :: x
    real(dp), intent(inout) :: data(:)
    real(dp) :: a, k
    a = data(1)
    k = data(2)
    y = a*sin(k*x)
    end function

    subroutine foo(a, k)
    real(dp) :: a, k
    real(dp) :: data(2)
    data(1) = a
    data(2) = k
    print *, simpson(f, 0._dp, pi, data)
    print *, simpson(f, 0._dp, 2*pi, data)
    end subroutine

    end module


II: General Structure
~~~~~~~~~~~~~~~~~~~~~

Define general structure or two which encompass the variations you actually
need (or are even remotely likely to need going forward). This single structure
type or two can then change if needed as future needs/ideas permit but won't
likely need to change from passing, say, real numbers to, say, and
instantiation of a text editor.

Integrator::

    module integrals
    use types, only: dp
    implicit none
    private
    public simpson, context

    type context
        ! This would be adjusted according to the problem to be solved.
        ! For example:
        real(dp) :: a, b, c, d
        integer :: i, j, k, l
        real(dp), pointer :: x(:), y(:)
        integer, pointer :: z(:)
    end type

    contains

    real(dp) function simpson(f, a, b, data) result(s)
    real(dp), intent(in) :: a, b
    interface
        real(dp) function func(x, data)
        use types, only: dp
        implicit none
        real(dp), intent(in) :: x
        type(context), intent(inout) :: data
        end function
    end interface
    procedure(func) :: f
    type(context), intent(inout) :: data
    s = (b-a) / 6 * (f(a, data) + 4*f((a+b)/2, data) + f(b, data))
    end function

    end module

Usage::

    module test
    use types, only: dp
    use integrals, only: simpson, context
    implicit none
    private
    public foo

    contains

    real(dp) function f(x, data) result(y)
    real(dp), intent(in) :: x
    type(context), intent(inout) :: data
    real(dp) :: a, k
    a = data%a
    k = data%b
    y = a*sin(k*x)
    end function

    subroutine foo(a, k)
    real(dp) :: a, k
    type(context) :: data
    data%a = a
    data%b = k
    print *, simpson(f, 0._dp, pi, data)
    print *, simpson(f, 0._dp, 2*pi, data)
    end subroutine

    end module


There is only so much flexibility really
needed. For example, you could define two structure types for this purpose, one
for Schroedinger and one for Dirac. Each would then be sufficiently general and
contain all the needed pieces with all the right labels.

Point is: it needn't
be "one abstract type to encompass all" or bust. There are natural and viable
options between "all" and "none".

III: Private Module Variables
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Hide the variable arguments completely by passing in module variables.

Integrator::

    module integrals
    use types, only: dp
    implicit none
    private
    public simpson

    contains

    real(dp) function simpson(f, a, b) result(s)
    real(dp), intent(in) :: a, b
    interface
        real(dp) function func(x)
        use types, only: dp
        implicit none
        real(dp), intent(in) :: x
        end function
    end interface
    procedure(func) :: f
    s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b))
    end function

    end module

Usage::

    module test
    use types, only: dp
    use integrals, only: simpson
    implicit none
    private
    public foo

    real(dp) :: global_a, global_k

    contains

    real(dp) function f(x) result(y)
    real(dp), intent(in) :: x
    y = global_a*sin(global_k*x)
    end function

    subroutine foo(a, k)
    real(dp) :: a, k
    global_a = a
    global_k = k
    print *, simpson(f, 0._dp, pi)
    print *, simpson(f, 0._dp, 2*pi)
    end subroutine

    end module


However it is best to avoid such global variables -- even though really just
semi-global -- if possible. But sometimes it may be the simplest cleanest way.
However, with a bit of thought, usually there is a better, safer, more explicit
way along the lines of II or IV.

IV: Nested functions
~~~~~~~~~~~~~~~~~~~~

Integrator::

    module integrals
    use types, only: dp
    implicit none
    private
    public simpson

    contains

    real(dp) function simpson(f, a, b) result(s)
    real(dp), intent(in) :: a, b
    interface
        real(dp) function func(x)
        use types, only: dp
        implicit none
        real(dp), intent(in) :: x
        end function
    end interface
    procedure(func) :: f
    s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b))
    end function

    end module

Usage::

    subroutine foo(a, k)
    use integrals, only: simpson
    real(dp) :: a, k
    print *, simpson(f, 0._dp, pi)
    print *, simpson(f, 0._dp, 2*pi)

    contains

    real(dp) function f(x) result(y)
    real(dp), intent(in) :: x
    y = a*sin(k*x)
    end function f

    end subroutine foo



V: Using type(c_ptr) Pointer
~~~~~~~~~~~~~~~~~~~~~~~~~~~~

In C, one would use the ``void *`` pointer. In Fortran, one
can use ``type(c_ptr)`` for exactly the same purpose.

Integrator::

    module integrals
    use types, only: dp
    use iso_c_binding, only: c_ptr
    implicit none
    private
    public simpson

    contains

    real(dp) function simpson(f, a, b, data) result(s)
    real(dp), intent(in) :: a, b
    interface
        real(dp) function func(x, data)
        use types, only: dp
        implicit none
        real(dp), intent(in) :: x
        type(c_ptr), intent(in) :: data
        end function
    end interface
    procedure(func) :: f
    type(c_ptr), intent(in) :: data
    s = (b-a) / 6 * (f(a, data) + 4*f((a+b)/2, data) + f(b, data))
    end function

    end module

Usage::

    module test
    use types, only: dp
    use integrals, only: simpson
    use iso_c_binding, only: c_ptr, c_loc, c_f_pointer
    implicit none
    private
    public foo

    type f_data
        ! Only contains data that we need for our particular callback.
        real(dp) :: a, k
    end type

    contains

    real(dp) function f(x, data) result(y)
    real(dp), intent(in) :: x
    type(c_ptr), intent(in) :: data
    type(f_data), pointer :: d
    call c_f_pointer(data, d)
    y = d%a * sin(d%k * x)
    end function

    subroutine foo(a, k)
    real(dp) :: a, k
    type(f_data), target :: data
    data%a = a
    data%k = k
    print *, simpson(f, 0._dp, pi, c_loc(data))
    print *, simpson(f, 0._dp, 2*pi, c_loc(data))
    end subroutine

    end module

As always, with the advantages of such re-casting, as Fortran lets you
do if you really want to, come also the disadvantages that fewer compile- and
run-time checks are possible to catch errors; and with that, inevitably more
leaky, bug-prone code. So one always has to balance the costs and benefits.

Usually, in the context of scientific programming, where the main thrust
is to represent and solve precise mathematical formulations (as opposed to
create a GUI with some untold number of buttons, drop-downs, and other
interface elements), simplest, least bug-prone, and fastest is to use one of
the previous approaches.

VI: transfer() Intrinsic Function
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Before Fortran 2003, the only way to do type casting
was using the ``transfer`` intrinsic function. It is functionally equivalent to
the method V, but more verbose and more error prone.
It is now obsolete and one should use the method V instead.

Examples:

http://jblevins.org/log/transfer

http://jblevins.org/research/generic-list.pdf

http://www.macresearch.org/advanced_fortran_90_callbacks_with_the_transfer_function

VII: Object Oriented Approach
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The module::

    module integrals

    use types, only: dp
    implicit none
    private

    public :: integrand, simpson

    ! User extends this type
    type, abstract :: integrand
    contains
        procedure(func), deferred :: eval
    end type

    abstract interface
        function func(this, x) result(fx)
        import :: integrand, dp
        class(integrand) :: this
        real(dp), intent(in) :: x
        real(dp) :: fx
        end function
    end interface

    contains

    real(dp) function simpson(f, a, b) result(s)
    class(integrand) :: f
    real(dp), intent(in) :: a, b
    s = ((b-a)/6) * (f%eval(a) + 4*f%eval((a+b)/2) + f%eval(b))
    end function

    end module

The abstract type prescribes exactly what the integration routine
needs, namely a method to evaluate the function, but imposes nothing
else on the user.  The user extends this type, providing a concrete
implementation of the eval type bound procedure and adding necessary
context data as components of the extended type.

Usage::

    module example_usage

    use types, only: dp
    use integrals, only: integrand, simpson
    implicit none
    private

    public :: foo

    type, extends(integrand) :: my_integrand
        real(dp) :: a, k
    contains
        procedure :: eval => f
    end type

    contains

    function f(this, x) result(fx)
    class(my_integrand) :: this
    real(dp), intent(in) :: x
    real(dp) :: fx
    fx = this%a*sin(this%k*x)
    end function

    subroutine foo(a, k)
    real(dp) :: a, k
    type(my_integrand) :: my_f
    my_f%a = a
    my_f%k = k
    print *, simpson(my_f, 0.0_dp, 1.0_dp)
    print *, simpson(my_f, 0.0_dp, 2.0_dp)
    end subroutine

    end module

Complete Example of void * vs type(c_ptr) and transfer()
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Here are three equivalent codes: one in C using ``void *`` and two codes in
Fortran using ``type(c_ptr)`` and ``transfer()``:

========  ===============   ===============================
Language  Method            Link
========  ===============   ===============================
C         ``void *``        https://gist.github.com/1665641
Fortran   ``type(c_ptr)``   https://gist.github.com/1665626
Fortran   ``transfer()``    https://gist.github.com/1665630
========  ===============   ===============================

The C code uses the standard C approach for writing extensible libraries that
accept callbacks and contexts. The two Fortran codes show how to do the same.
The ``type(c_ptr)`` method is equivalent to the C version and that is the
approach that should be used.

The ``transfer()`` method is here for completeness only (before Fortran 2003,
it was the only way) and it is a little cumbersome, because the user needs to
create auxiliary conversion functions for each of his types.
As such, the ``type(c_ptr)`` method should be used instead.

.. _parallel:

Parallel programming
--------------------

OpenMP
~~~~~~

`OpenMP <http://www.openmp.org/>`_ should be compatible with
non-openMP compilers. This can be enforced by prepending all
OpenMP-specific calls by ``!\$``. Regular compilers will consider these
lines as comments and ignore them. For OpenMP compilers, these lines
will be considered as regular Fortran code. The following code ::

    program test_openmpi
       !\$ use omp_lib
       implicit none

       integer :: nthreads

       nthreads = -1
       !\$ nthreads = omp_get_num_threads()

       ! will print the number of running threads when compiled with OpenMP, else will print -1
       print*, "nthreads=", nthreads
    end program

will print the number of threads used when compiled with OpenMP. It will print by default -1 if compiled without OpenMP.

MPI
~~~

There are three ways of including MPI in a fortran program:

=============== ==================== ====================================================================
Fortran version  Method              Comments
=============== ==================== ====================================================================
Fortran 08      ``use mpi_f08``      Consistent with F08 standard, good type-checking; highly recommended
Fortran 90      ``use mpi``          Not consistent with standard, so-so type-checking; not recommended
Fortran 77      ``include "mpif.h"`` Not consistent with standard, no type-checking; strongly discouraged
=============== ==================== ====================================================================

On infrastructures where ``use mpi_f08`` is not available, one should
fallback to ``use mpi``. The use of ``include "mpif.h"`` is strongly
discouraged, as it does not check at all the types of the argument or
that the function calls provide the good arguments. For example, you
don’t get any compiler warnings if you call a subroutine and forget a
parameter, add an extra parameter, or pass a parameter of the wrong
type. It may also lead to silent data corruption.