Transferring fortran90.org "Fortran Best Practices" into a mini-book

_data/learning.yml

@@ -53,6 +53,26 @@ books:
       - link: /learn/os_setup/ides
       - link: /learn/os_setup/tips
 
+  - title: Fortran Best Practices
+    description: This tutorial collects a modern canonical way of doing things in Fortran.
+    category: Getting started
+    link: /learn/best_practices
+    pages:
+        - link: /learn/best_practices/style_guide
+        - link: /learn/best_practices/floating_point
+        - link: /learn/best_practices/integer_division
+        - link: /learn/best_practices/modules_programs
+        - link: /learn/best_practices/arrays
+        - link: /learn/best_practices/multidim_arrays
+        - link: /learn/best_practices/element_operations
+        - link: /learn/best_practices/allocatable_arrays
+        - link: /learn/best_practices/file_io
+        - link: /learn/best_practices/c_interfacing
+        - link: /learn/best_practices/python_interfacing
+        - link: /learn/best_practices/callbacks
+        - link: /learn/best_practices/type_casting
+        - link: /learn/best_practices/parallel_programming
+
 # Web links listed at the bottom of the 'Learn' landing page
 #
 reference-links:

learn/best_practices/allocatable_arrays.md

@@ -0,0 +1,34 @@
+---
+layout: book
+title: Allocatable Arrays
+permalink: /learn/best_practices/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:
+
+``` fortran
+subroutine foo(lam)
+real(dp), allocatable, intent(out) :: lam(:)
+allocate(lam(5))
+end subroutine
+```
+
+And use somewhere else:
+
+``` fortran
+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):
+
+``` fortran
+if (allocated(lam)) deallocate(lam)
+allocate(lam(10))
+```

learn/best_practices/arrays.md

@@ -0,0 +1,128 @@
+---
+layout: book
+title: Arrays
+permalink: /learn/best_practices/arrays
+---
+
+When passing arrays in and out of a subroutine/function, use the
+following pattern for 1D arrays (it is called <span
+class="title-ref">assumed-shape</span>):
+
+``` fortran
+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:
+
+``` fortran
+subroutine g(A)
+real(dp), intent(in) :: A(:, :)
+...
+end subroutine
+```
+
+and call it like this:
+
+``` fortran
+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
+<span class="title-ref">explicit-shape</span> 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 `elemental` for more
+    information)
+
+To use <span class="title-ref">explicit-shape</span> arrays, do:
+
+``` fortran
+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:
+
+``` fortran
+subroutine g(m, n, A)
+integer, intent(in) :: m, n
+real(dp), intent(in) :: A(m, n)
+...
+end subroutine
+```
+
+and call it like this:
+
+``` fortran
+real(dp) :: r(5)
+call f(size(r), r)
+```
+
+In order to return an array from a function, do:
+
+``` fortran
+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:
+
+``` fortran
+if (size(r) /= 4) stop "Incorrect size of 'r'"
+```
+
+To initialize an array, do:
+
+``` fortran
+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:
+
+``` fortran
+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
+```

learn/best_practices/c_interfacing.md

@@ -0,0 +1,60 @@
+---
+layout: book
+title: Interfacing with C
+permalink: /learn/best_practices/c_interfacing
+---
+
+Write a C wrapper using the `iso_c_binding` module:
+
+``` fortran
+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:
+
+``` c
+void c_mesh_exp(double *r_min, double *r_max, double *a, int *N,
+        double *mesh);
+```
+
+use it as:
+
+``` 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.

learn/best_practices/callbacks.md

@@ -0,0 +1,117 @@
+---
+layout: book
+title: Callbacks
+permalink: /learn/best_practices/callbacks
+---
+
+There are two ways to implement callbacks to be used like this:
+
+``` fortran
+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:
+
+``` fortran
+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:
+
+``` fortran
+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:
+
+``` fortran
+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
+```