From 6f3648dbd784079826a7ca1183e77929007c8489 Mon Sep 17 00:00:00 2001 From: vmagnin Date: Wed, 21 Apr 2021 11:10:19 +0200 Subject: [PATCH 01/24] Transfer fortran90.org Best Practices content (#79) --- _data/learning.yml | 20 + learn/best_practices/allocatable_arrays.md | 34 ++ learn/best_practices/arrays.md | 128 +++++ learn/best_practices/c_interfacing.md | 60 +++ learn/best_practices/callbacks.md | 117 +++++ learn/best_practices/element_operations.md | 128 +++++ learn/best_practices/file_io.md | 59 +++ learn/best_practices/floating_point.md | 48 ++ learn/best_practices/index.md | 21 + learn/best_practices/integer_division.md | 16 + learn/best_practices/modules_programs.md | 60 +++ learn/best_practices/multidim_arrays.md | 83 ++++ learn/best_practices/parallel_programming.md | 51 ++ learn/best_practices/python_interfacing.md | 70 +++ learn/best_practices/style_guide.md | 57 +++ learn/best_practices/type_casting.md | 494 +++++++++++++++++++ 16 files changed, 1446 insertions(+) create mode 100644 learn/best_practices/allocatable_arrays.md create mode 100644 learn/best_practices/arrays.md create mode 100644 learn/best_practices/c_interfacing.md create mode 100644 learn/best_practices/callbacks.md create mode 100644 learn/best_practices/element_operations.md create mode 100644 learn/best_practices/file_io.md create mode 100644 learn/best_practices/floating_point.md create mode 100644 learn/best_practices/index.md create mode 100644 learn/best_practices/integer_division.md create mode 100644 learn/best_practices/modules_programs.md create mode 100644 learn/best_practices/multidim_arrays.md create mode 100644 learn/best_practices/parallel_programming.md create mode 100644 learn/best_practices/python_interfacing.md create mode 100644 learn/best_practices/style_guide.md create mode 100644 learn/best_practices/type_casting.md diff --git a/_data/learning.yml b/_data/learning.yml index def5dd289..51a9df001 100644 --- a/_data/learning.yml +++ b/_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: diff --git a/learn/best_practices/allocatable_arrays.md b/learn/best_practices/allocatable_arrays.md new file mode 100644 index 000000000..2869d2aa6 --- /dev/null +++ b/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)) +``` diff --git a/learn/best_practices/arrays.md b/learn/best_practices/arrays.md new file mode 100644 index 000000000..55c02e9d3 --- /dev/null +++ b/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 assumed-shape): + +``` 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 +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 `elemental` for more + information) + +To use explicit-shape 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 +``` diff --git a/learn/best_practices/c_interfacing.md b/learn/best_practices/c_interfacing.md new file mode 100644 index 000000000..6e41ef70a --- /dev/null +++ b/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. diff --git a/learn/best_practices/callbacks.md b/learn/best_practices/callbacks.md new file mode 100644 index 000000000..a52740353 --- /dev/null +++ b/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 +``` diff --git a/learn/best_practices/element_operations.md b/learn/best_practices/element_operations.md new file mode 100644 index 000000000..f8b9ea791 --- /dev/null +++ b/learn/best_practices/element_operations.md @@ -0,0 +1,128 @@ +--- +layout: book +title: Element-wise Operations on Arrays Using Subroutines/Functions +permalink: /learn/best_practices/element_operations +--- + +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: + +``` fortran +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: + +``` fortran +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: + +``` fortran +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: + +``` fortran +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: + +``` fortran +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: + +``` fortran +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: + +``` fortran +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: + +``` fortran +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: + +``` fortran +3.0000000000000000 +1.0000000000000000 2.0000000000000000 3.0000000000000000 3.1622776601683795 +1.0000000000000000 2.0000000000000000 3.0000000000000000 3.1622776601683795 +``` diff --git a/learn/best_practices/file_io.md b/learn/best_practices/file_io.md new file mode 100644 index 000000000..a396444d0 --- /dev/null +++ b/learn/best_practices/file_io.md @@ -0,0 +1,59 @@ +--- +layout: book +title: File Input/Output +permalink: /learn/best_practices/file_io +--- + +To read from a file: + +``` fortran +integer :: u +open(newunit=u, file="log.txt", status="old") +read(u, *) a, b +close(u) +``` + +Write to a file: + +``` fortran +integer :: u +open(newunit=u, file="log.txt", status="replace") +write(u, *) a, b +close(u) +``` + +To append to an existing file: + +``` fortran +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: + +``` fortran +open(newunit(u), ...) +``` + +where the `newunit` function is defined by: + +``` fortran +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 +``` diff --git a/learn/best_practices/floating_point.md b/learn/best_practices/floating_point.md new file mode 100644 index 000000000..ee555a4b3 --- /dev/null +++ b/learn/best_practices/floating_point.md @@ -0,0 +1,48 @@ +--- +layout: book +title: Floating Point Numbers +permalink: /learn/best_practices/floating_point +--- + +Somewhere create and export a parameter `dp`: + +``` fortran +integer, parameter:: dp=kind(0.d0) ! double precision +``` + +and declare floats as: + +``` fortran +real(dp) :: a, b, c +``` + +Always write all floating point constants with the `_dp` suffix: + +``` fortran +1.0_dp, 3.5_dp, 1.34e8_dp +``` + +and never any other way (see also the gotcha +`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 +). + +It is safe to assign integers to floating point numbers without losing +accuracy: + +``` fortran +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: + +``` fortran +real(dp) :: a +a = real(1, dp) / 2 ! 'a' is equal to 0.5_dp +``` diff --git a/learn/best_practices/index.md b/learn/best_practices/index.md new file mode 100644 index 000000000..40f2b8055 --- /dev/null +++ b/learn/best_practices/index.md @@ -0,0 +1,21 @@ +--- +layout: book +title: Fortran Best Practices +permalink: /learn/best_practices +author: Ondřej Čertík, John Pask, Jed Brown, Matthew Emmett, Juan Luis Cano Rodríguez, Neil Carlson, Andrea Vigliotti, Pierre Haessig +--- + +This mini-book 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: + + +Language features are explained at: + + +The best resource is a recent Fortran standard, for example the Fortran +2003 standard: diff --git a/learn/best_practices/integer_division.md b/learn/best_practices/integer_division.md new file mode 100644 index 000000000..4377b976f --- /dev/null +++ b/learn/best_practices/integer_division.md @@ -0,0 +1,16 @@ +--- +layout: book +title: Integer Division +permalink: /learn/best_practices/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: + +``` fortran +(N - 1.0_dp)/N +``` + +As long as one of the `/` operands is a float, a floating point division +will be used. diff --git a/learn/best_practices/modules_programs.md b/learn/best_practices/modules_programs.md new file mode 100644 index 000000000..f6e21331a --- /dev/null +++ b/learn/best_practices/modules_programs.md @@ -0,0 +1,60 @@ +--- +layout: book +title: Modules and Programs +permalink: /learn/best_practices/modules_programs +--- + +Only use modules and programs. Always setup a module in the following +way: + +``` fortran +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: + +``` fortran +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: + +``` fortran +use fmesh +``` + +But just like in Python, this should be avoided ("explicit is better +than implicit") in most cases. diff --git a/learn/best_practices/multidim_arrays.md b/learn/best_practices/multidim_arrays.md new file mode 100644 index 000000000..82a336c64 --- /dev/null +++ b/learn/best_practices/multidim_arrays.md @@ -0,0 +1,83 @@ +--- +layout: book +title: Multidimensional Arrays +permalink: /learn/best_practices/multidim_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: + +``` fortran +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: + +``` fortran +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: + +``` fortran +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: + +``` fortran +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: + +``` fortran +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: + +``` fortran +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. diff --git a/learn/best_practices/parallel_programming.md b/learn/best_practices/parallel_programming.md new file mode 100644 index 000000000..9e8625134 --- /dev/null +++ b/learn/best_practices/parallel_programming.md @@ -0,0 +1,51 @@ +--- +layout: book +title: Parallel programming +permalink: /learn/best_practices/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 : + +``` fortran +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. diff --git a/learn/best_practices/python_interfacing.md b/learn/best_practices/python_interfacing.md new file mode 100644 index 000000000..8abe8f0f8 --- /dev/null +++ b/learn/best_practices/python_interfacing.md @@ -0,0 +1,70 @@ +--- +layout: book +title: Interfacing with Python +permalink: /learn/best_practices/python_interfacing +--- + +Using Cython +------------ + +To wrap Fortran code in Python, export it to C first (see above) and +then write this Cython code: + +``` 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: + +``` fortran +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: + +``` 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))) +``` diff --git a/learn/best_practices/style_guide.md b/learn/best_practices/style_guide.md new file mode 100644 index 000000000..c2b4ba2e5 --- /dev/null +++ b/learn/best_practices/style_guide.md @@ -0,0 +1,57 @@ +--- +layout: book +title: Fortran Style Guide +permalink: /learn/best_practices/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). diff --git a/learn/best_practices/type_casting.md b/learn/best_practices/type_casting.md new file mode 100644 index 000000000..3aa1088f4 --- /dev/null +++ b/learn/best_practices/type_casting.md @@ -0,0 +1,494 @@ +--- +layout: book +title: Type Casting in Callbacks +permalink: /learn/best_practices/type_casting +--- + +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: + +``` fortran +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: + +``` fortran +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: + +``` fortran +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: + +``` fortran +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: + +``` 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 +``` + +Usage: + +``` fortran +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: + +``` 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 +``` + +Usage: + +``` 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 +``` + +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: + +``` fortran +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: + +``` fortran +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: + + + + + + + +VII: Object Oriented Approach +----------------------------- + +The module: + +``` fortran +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: + +``` fortran +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 *` | | +| Fortran | `type(c_ptr)`   | | +| Fortran | `transfer()` | | + +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. From b713a9bf2a4223d5a67d7182e2290d83d81705fc Mon Sep 17 00:00:00 2001 From: vmagnin Date: Wed, 21 Apr 2021 14:26:56 +0200 Subject: [PATCH 02/24] Removed vestigial best_practices.md (2020-04-15) --- learn/best_practices.md | 5 ----- 1 file changed, 5 deletions(-) delete mode 100644 learn/best_practices.md diff --git a/learn/best_practices.md b/learn/best_practices.md deleted file mode 100644 index 7ecfc63ef..000000000 --- a/learn/best_practices.md +++ /dev/null @@ -1,5 +0,0 @@ ---- -layout: page -title: Fortran best practices -permalink: /learn/best_practices ---- \ No newline at end of file From 16b65a5dcbe3c08baf7d30d5d402d38e9bf4ac5f Mon Sep 17 00:00:00 2001 From: vmagnin Date: Sat, 24 Apr 2021 14:25:12 +0200 Subject: [PATCH 03/24] Replaced HTML code by Markdown italics *...* --- learn/best_practices/arrays.md | 8 +++----- learn/best_practices/element_operations.md | 4 ++-- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/learn/best_practices/arrays.md b/learn/best_practices/arrays.md index 55c02e9d3..08bfeecac 100644 --- a/learn/best_practices/arrays.md +++ b/learn/best_practices/arrays.md @@ -5,8 +5,7 @@ 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 assumed-shape): +following pattern for 1D arrays (it is called *assumed-shape*): ``` fortran subroutine f(r) @@ -43,8 +42,7 @@ No array copying is done above. It has the following advantages: 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: +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) @@ -52,7 +50,7 @@ subroutines. However in the following cases one can (or has to) use there are also other ways to do that, see `elemental` for more information) -To use explicit-shape arrays, do: +To use *explicit-shape* arrays, do: ``` fortran subroutine f(n, r) diff --git a/learn/best_practices/element_operations.md b/learn/best_practices/element_operations.md index f8b9ea791..831e07d78 100644 --- a/learn/best_practices/element_operations.md +++ b/learn/best_practices/element_operations.md @@ -7,7 +7,7 @@ permalink: /learn/best_practices/element_operations There are three approaches: - `elemental` subroutines -- explicit-shape arrays +- *explicit-shape* arrays - implementing the operation for vectors and write simple wrapper subroutines (that use `reshape` internally) for each array shape @@ -99,7 +99,7 @@ This will print: 1.0000000000000000 2.0000000000000000 3.0000000000000000 3.1622776601683795 ``` -Or one can use explicit-shape arrays as +Or one can use *explicit-shape* arrays as follows: ``` fortran From 1a256d33d7970606b236a0e5b093258fca7b77b1 Mon Sep 17 00:00:00 2001 From: vmagnin Date: Sat, 24 Apr 2021 14:27:09 +0200 Subject: [PATCH 04/24] Typo: replaced 3 'enddo' by 'end do' --- learn/best_practices/arrays.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/learn/best_practices/arrays.md b/learn/best_practices/arrays.md index 08bfeecac..2c5d856b8 100644 --- a/learn/best_practices/arrays.md +++ b/learn/best_practices/arrays.md @@ -14,7 +14,7 @@ integer :: n, i n = size(r) do i = 1, n r(i) = 1.0_dp / i**2 -enddo +end do end subroutine ``` @@ -59,7 +59,7 @@ real(dp), intent(out) :: r(n) integer :: i do i = 1, n r(i) = 1.0_dp / i**2 -enddo +end do end subroutine ``` @@ -89,7 +89,7 @@ real(dp) :: r(n) integer :: i do i = 1, n r(i) = 1.0_dp / i**2 -enddo +end do end function ``` From 6c15deac7decab03067e8e910d9462d5e683510d Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Sat, 8 May 2021 11:50:03 +0200 Subject: [PATCH 05/24] Consistent indentation (2 spaces) with quickstart minibook --- learn/best_practices/allocatable_arrays.md | 4 +- learn/best_practices/arrays.md | 50 +-- learn/best_practices/c_interfacing.md | 18 +- learn/best_practices/callbacks.md | 74 ++--- learn/best_practices/element_operations.md | 42 +-- learn/best_practices/file_io.md | 20 +- learn/best_practices/modules_programs.md | 36 +-- learn/best_practices/multidim_arrays.md | 36 +-- learn/best_practices/parallel_programming.md | 14 +- learn/best_practices/type_casting.md | 306 +++++++++---------- 10 files changed, 300 insertions(+), 300 deletions(-) diff --git a/learn/best_practices/allocatable_arrays.md b/learn/best_practices/allocatable_arrays.md index 2869d2aa6..3a8d8cf53 100644 --- a/learn/best_practices/allocatable_arrays.md +++ b/learn/best_practices/allocatable_arrays.md @@ -11,8 +11,8 @@ For example you can allocate it inside a subroutine: ``` fortran subroutine foo(lam) -real(dp), allocatable, intent(out) :: lam(:) -allocate(lam(5)) + real(dp), allocatable, intent(out) :: lam(:) + allocate(lam(5)) end subroutine ``` diff --git a/learn/best_practices/arrays.md b/learn/best_practices/arrays.md index 2c5d856b8..5330d2147 100644 --- a/learn/best_practices/arrays.md +++ b/learn/best_practices/arrays.md @@ -9,12 +9,12 @@ following pattern for 1D arrays (it is called *assumed-shape*): ``` fortran subroutine f(r) -real(dp), intent(out) :: r(:) -integer :: n, i -n = size(r) -do i = 1, n + real(dp), intent(out) :: r(:) + integer :: n, i + n = size(r) + do i = 1, n r(i) = 1.0_dp / i**2 -end do + end do end subroutine ``` @@ -22,8 +22,8 @@ end subroutine ``` fortran subroutine g(A) -real(dp), intent(in) :: A(:, :) -... + real(dp), intent(in) :: A(:, :) + ... end subroutine ``` @@ -54,12 +54,12 @@ To use *explicit-shape* arrays, do: ``` fortran subroutine f(n, r) -integer, intent(in) :: n -real(dp), intent(out) :: r(n) -integer :: i -do i = 1, n + integer, intent(in) :: n + real(dp), intent(out) :: r(n) + integer :: i + do i = 1, n r(i) = 1.0_dp / i**2 -end do + end do end subroutine ``` @@ -67,9 +67,9 @@ end subroutine ``` fortran subroutine g(m, n, A) -integer, intent(in) :: m, n -real(dp), intent(in) :: A(m, n) -... + integer, intent(in) :: m, n + real(dp), intent(in) :: A(m, n) + ... end subroutine ``` @@ -84,12 +84,12 @@ 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 + integer, intent(in) :: n + real(dp) :: r(n) + integer :: i + do i = 1, n r(i) = 1.0_dp / i**2 -end do + end do end function ``` @@ -115,12 +115,12 @@ 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, intent(in) :: kappa_min + real(dp), intent(in) :: lam(kappa_min:) -integer :: kappa -do kappa = kappa_min, ubound(lam, 1) + integer :: kappa + do kappa = kappa_min, ubound(lam, 1) print *, kappa, lam(kappa) -end do + end do end subroutine ``` diff --git a/learn/best_practices/c_interfacing.md b/learn/best_practices/c_interfacing.md index 6e41ef70a..6617eec19 100644 --- a/learn/best_practices/c_interfacing.md +++ b/learn/best_practices/c_interfacing.md @@ -9,20 +9,20 @@ 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 + use iso_c_binding, only: c_double, c_int + use fmesh, only: mesh_exp -implicit none + 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) + 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 diff --git a/learn/best_practices/callbacks.md b/learn/best_practices/callbacks.md index a52740353..1a81f971d 100644 --- a/learn/best_practices/callbacks.md +++ b/learn/best_practices/callbacks.md @@ -8,16 +8,16 @@ 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) + 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) + real(dp), intent(in) :: x + y = a*sin(k*x) end function f end subroutine foo @@ -28,23 +28,23 @@ a subroutine/function using: ``` fortran module integrals -use types, only: dp -implicit none -private -public simpson + 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), 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 interface + s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b)) end function end module @@ -56,24 +56,24 @@ argument: ``` fortran module integrals -use types, only: dp -implicit none -private -public simpson + 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), 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 interface + procedure(func) :: f + s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b)) end function end module @@ -84,33 +84,33 @@ like: ``` fortran module integrals -use types, only: dp -implicit none -private -public simpson + use types, only: dp + implicit none + private + public simpson -interface + interface real(dp) function func(x) use types, only: dp implicit none real(dp), intent(in) :: x end function -end interface + 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)) + 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) + 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 diff --git a/learn/best_practices/element_operations.md b/learn/best_practices/element_operations.md index 831e07d78..13b12de74 100644 --- a/learn/best_practices/element_operations.md +++ b/learn/best_practices/element_operations.md @@ -16,9 +16,9 @@ function like this: ``` fortran real(dp) elemental function nroot(n, x) result(y) -integer, intent(in) :: n -real(dp), intent(in) :: x -y = x**(1._dp / n) + integer, intent(in) :: n + real(dp), intent(in) :: x + y = x**(1._dp / n) end function ``` @@ -60,26 +60,26 @@ other array shapes: ``` fortran 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) + 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) + 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)]) + 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 ``` @@ -104,10 +104,10 @@ follows: ``` fortran 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) + integer, intent(in) :: n, k + real(dp), intent(in) :: x(k) + real(dp) :: y(k) + y = x**(1._dp / n) end function ``` diff --git a/learn/best_practices/file_io.md b/learn/best_practices/file_io.md index a396444d0..1e15877ff 100644 --- a/learn/best_practices/file_io.md +++ b/learn/best_practices/file_io.md @@ -42,18 +42,18 @@ where the `newunit` function is defined by: ``` fortran 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 + ! 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 + if (present(unit)) unit=n + return end if -end do -call stop_error("newunit ERROR: available unit not found.") + end do + call stop_error("newunit ERROR: available unit not found.") end function ``` diff --git a/learn/best_practices/modules_programs.md b/learn/best_practices/modules_programs.md index f6e21331a..cab56cb1b 100644 --- a/learn/best_practices/modules_programs.md +++ b/learn/best_practices/modules_programs.md @@ -9,19 +9,19 @@ way: ``` fortran 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 + 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 ... @@ -37,15 +37,15 @@ Setup programs in the following way: ``` fortran 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" + 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 ``` diff --git a/learn/best_practices/multidim_arrays.md b/learn/best_practices/multidim_arrays.md index 82a336c64..a981898db 100644 --- a/learn/best_practices/multidim_arrays.md +++ b/learn/best_practices/multidim_arrays.md @@ -24,11 +24,11 @@ index last ("Inner-most are left-most."). So the elements of a 3D array ``` fortran do i3 = 1, N3 - do i2 = 1, N2 - do i1 = 1, N1 - A(i1, i2, i3) - end do + do i2 = 1, N2 + do i1 = 1, N1 + A(i1, i2, i3) end do + end do end do ``` @@ -36,9 +36,9 @@ Associated array of vectors would then be most efficiently accessed as: ``` fortran do i3 = 1, N3 - do i2 = 1, N2 - A(:, i2, i3) - end do + do i2 = 1, N2 + A(:, i2, i3) + end do end do ``` @@ -46,7 +46,7 @@ And associated set of matrices would be most efficiently accessed as: ``` fortran do i3 = 1, N3 - A(:, :, i3) + A(:, :, i3) end do ``` @@ -58,13 +58,13 @@ strides, for example the first loop would become: ``` fortran do i3 = 1, N3 - Ai3 = A(:, :, i3) - do i2 = 1, N2 - Ai2i3 = Ai3(:, i2) - do i1 = 1, N1 - Ai2i3(i1) - end do + Ai3 = A(:, :, i3) + do i2 = 1, N2 + Ai2i3 = Ai3(:, i2) + do i1 = 1, N1 + Ai2i3(i1) end do + end do end do ``` @@ -72,10 +72,10 @@ the second loop would become: ``` fortran do i3 = 1, N3 - Ai3 = A(:, :, i3) - do i2 = 1, N2 - Ai3(:, i2) - end do + Ai3 = A(:, :, i3) + do i2 = 1, N2 + Ai3(:, i2) + end do end do ``` diff --git a/learn/best_practices/parallel_programming.md b/learn/best_practices/parallel_programming.md index 9e8625134..ef3ebe96a 100644 --- a/learn/best_practices/parallel_programming.md +++ b/learn/best_practices/parallel_programming.md @@ -15,16 +15,16 @@ regular Fortran code. The following code : ``` fortran program test_openmpi - !$ use omp_lib - implicit none + !$ use omp_lib + implicit none - integer :: nthreads + integer :: nthreads - nthreads = -1 - !$ nthreads = omp_get_num_threads() + 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 + ! will print the number of running threads when compiled with OpenMP, else will print -1 + print*, "nthreads=", nthreads end program ``` diff --git a/learn/best_practices/type_casting.md b/learn/best_practices/type_casting.md index 3aa1088f4..d56de6a08 100644 --- a/learn/best_practices/type_casting.md +++ b/learn/best_practices/type_casting.md @@ -22,26 +22,26 @@ Integrator: ``` fortran module integrals -use types, only: dp -implicit none -private -public simpson + 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), 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 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 @@ -51,30 +51,30 @@ Usage: ``` fortran module test -use types, only: dp -use integrals, only: simpson -implicit none -private -public foo + 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) + 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) + 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 @@ -93,35 +93,35 @@ Integrator: ``` fortran module integrals -use types, only: dp -implicit none -private -public simpson, context + use types, only: dp + implicit none + private + public simpson, context -type 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 + end type contains real(dp) function simpson(f, a, b, data) result(s) -real(dp), intent(in) :: a, b -interface + 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 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 @@ -131,30 +131,30 @@ Usage: ``` fortran module test -use types, only: dp -use integrals, only: simpson, context -implicit none -private -public foo + 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) + 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) + 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 @@ -177,24 +177,24 @@ Integrator: ``` fortran module integrals -use types, only: dp -implicit none -private -public simpson + 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), 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 interface + procedure(func) :: f + s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b)) end function end module @@ -204,27 +204,27 @@ Usage: ``` fortran module test -use types, only: dp -use integrals, only: simpson -implicit none -private -public foo + use types, only: dp + use integrals, only: simpson + implicit none + private + public foo -real(dp) :: global_a, global_k + 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) + 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) + 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 @@ -242,24 +242,24 @@ Integrator: ``` fortran module integrals -use types, only: dp -implicit none -private -public simpson + 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), 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 interface + procedure(func) :: f + s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b)) end function end module @@ -294,27 +294,27 @@ Integrator: ``` fortran module integrals -use types, only: dp -use iso_c_binding, only: c_ptr -implicit none -private -public simpson + 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), 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 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 @@ -324,35 +324,35 @@ Usage: ``` fortran 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 + 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 + 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) + 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)) + 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 @@ -394,33 +394,33 @@ The module: ``` fortran module integrals -use types, only: dp -implicit none -private + use types, only: dp + implicit none + private -public :: integrand, simpson + public :: integrand, simpson -! User extends this type -type, abstract :: integrand -contains + ! User extends this type + type, abstract :: integrand + contains procedure(func), deferred :: eval -end type + end type -abstract interface + abstract interface function func(this, x) result(fx) - import :: integrand, dp - class(integrand) :: this - real(dp), intent(in) :: x - real(dp) :: fx + import :: integrand, dp + class(integrand) :: this + real(dp), intent(in) :: x + real(dp) :: fx end function -end interface + 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)) + 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 @@ -437,35 +437,35 @@ Usage: ``` fortran module example_usage -use types, only: dp -use integrals, only: integrand, simpson -implicit none -private + use types, only: dp + use integrals, only: integrand, simpson + implicit none + private -public :: foo + public :: foo -type, extends(integrand) :: my_integrand + type, extends(integrand) :: my_integrand real(dp) :: a, k -contains + contains procedure :: eval => f -end type + 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) + 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) + 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 From c84e3f001fb39ef8f4594c3a5a09d88f6ad30c76 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Sat, 8 May 2021 12:01:28 +0200 Subject: [PATCH 06/24] Minor adjustments, wording, stop -> error stop --- learn/best_practices/arrays.md | 5 ++-- learn/best_practices/parallel_programming.md | 2 +- learn/best_practices/style_guide.md | 5 ++-- learn/best_practices/type_casting.md | 28 ++++++++++---------- 4 files changed, 21 insertions(+), 19 deletions(-) diff --git a/learn/best_practices/arrays.md b/learn/best_practices/arrays.md index 5330d2147..e1813c29c 100644 --- a/learn/best_practices/arrays.md +++ b/learn/best_practices/arrays.md @@ -76,8 +76,9 @@ end subroutine and call it like this: ``` fortran -real(dp) :: r(5) +real(dp) :: r(5), s(3, 4) call f(size(r), r) +call f(size(s), s) ``` In order to return an array from a function, do: @@ -97,7 +98,7 @@ 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'" +if (size(r) /= 4) error stop "Incorrect size of 'r'" ``` To initialize an array, do: diff --git a/learn/best_practices/parallel_programming.md b/learn/best_practices/parallel_programming.md index ef3ebe96a..ffce72588 100644 --- a/learn/best_practices/parallel_programming.md +++ b/learn/best_practices/parallel_programming.md @@ -14,7 +14,7 @@ ignore them. For OpenMP compilers, these lines will be considered as regular Fortran code. The following code : ``` fortran -program test_openmpi +program test_openmp !$ use omp_lib implicit none diff --git a/learn/best_practices/style_guide.md b/learn/best_practices/style_guide.md index c2b4ba2e5..b871714e7 100644 --- a/learn/best_practices/style_guide.md +++ b/learn/best_practices/style_guide.md @@ -40,8 +40,9 @@ 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). +Use a consistent indentation to make your code readable. +The amount of indentation is a matter of preference, the most common choices +are two, three or four spaces. Comparison to Other Languages ----------------------------- diff --git a/learn/best_practices/type_casting.md b/learn/best_practices/type_casting.md index d56de6a08..80be54185 100644 --- a/learn/best_practices/type_casting.md +++ b/learn/best_practices/type_casting.md @@ -11,8 +11,8 @@ 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 --------------- +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 -- @@ -80,8 +80,8 @@ end subroutine end module ``` -II: General Structure ---------------------- +General Structure +----------------- Define general structure or two which encompass the variations you actually need (or are even remotely likely to need going forward). This @@ -168,8 +168,8 @@ 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 ------------------------------ +Private Module Variables +------------------------ Hide the variable arguments completely by passing in module variables. @@ -235,8 +235,8 @@ 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 --------------------- +Nested functions +---------------- Integrator: @@ -284,8 +284,8 @@ end function f end subroutine foo ``` -V: Using type(c\_ptr) Pointer ------------------------------ +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. @@ -370,8 +370,8 @@ 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 ---------------------------------- +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 @@ -386,8 +386,8 @@ Examples: -VII: Object Oriented Approach ------------------------------ +Object Oriented Approach +------------------------ The module: From 33830397cacdf6278ca46d7c5df70f2be3155e81 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Tue, 1 Jun 2021 11:39:48 +0200 Subject: [PATCH 07/24] Extend the alloctable arrays section --- learn/best_practices/allocatable_arrays.md | 120 ++++++++++++++++++--- 1 file changed, 108 insertions(+), 12 deletions(-) diff --git a/learn/best_practices/allocatable_arrays.md b/learn/best_practices/allocatable_arrays.md index 3a8d8cf53..26591289f 100644 --- a/learn/best_practices/allocatable_arrays.md +++ b/learn/best_practices/allocatable_arrays.md @@ -4,31 +4,127 @@ 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. +The ``allocatable`` attribute provides safe way to for memory handling. +In comparison to variables with ``pointer`` attribute the memory is managed +automatically and will go deallocated automatically once the variable goes +out-of-scope. Using ``allocatable`` variables removes the possibility to +create memory leaks in an application. -For example you can allocate it inside a subroutine: +They can be used in subroutines to create scratch or work arrays, where +automatic arrays would become too large to fit on the stack. -``` fortran +```fortran +real(dp), allocatable :: temp(:) +allocate(temp(10)) +``` + +The allocation status can be checked using the ``allocated`` intrinsic +to avoid uninitialized access + +```fortran +subroutine show_arr(arr) + integer, allocatable, intent(in) :: arr(:) + + if (allocated(arr)) then + print *, arr + end if +end subroutine show_arr +``` + +To allocate variables inside a procedure the dummy argument has to carry +the ``allocatable`` attribute. Using it in combination with ``intent(out)`` +will deallocate previous allocations before entring the procedure: + +```fortran subroutine foo(lam) real(dp), allocatable, intent(out) :: lam(:) allocate(lam(5)) -end subroutine +end subroutine foo ``` -And use somewhere else: +The allocated array can be used afterwards like a normal array -``` fortran +```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): +An already allocated array cannot be allocated again without prior deallocation. +Similarly, deallocation can only be invoked for allocated arrays. To reallocate +an array use -``` fortran +```fortran if (allocated(lam)) deallocate(lam) allocate(lam(10)) ``` + +Passing allocated arrays to procedures does not require the ``allocatable`` attribute +for the dummy arguments anymore. + +```fortran +subroutine show_arr(arr) + integer, intent(in) :: arr(:) + + print *, arr +end subroutine show_arr + +subroutine proc + integer :: i + integer, allocatable :: arr + + allocate(arr(5)) + + do i = 1, size(arr) + arr(i) = 2*i + 1 + end do + call show_arr(arr) +end subroutine proc +``` + +Passing an unallocated array in this context will lead to an invalid memory access. +Allocatable arrays can be passed to ``optional`` dummy arguments, if they are unallocated +the argument will not be present. The ``allocatable`` attribute is not limited to +arrays and can also be associated with scalars, which can be useful in combination +with ``optional`` dummy arguments. + +Allocations can be moved between different arrays with ``allocatable`` attribute. + +```fortran +subroutine resize(var, n) + real(wp), allocatable, intent(inout) :: var(:) + integer, intent(in), optional :: n + integer :: this_size, new_size + integer, parameter :: inital_size = 16 + + if (allocated(var)) then + this_size = size(var, 1) + call move_alloc(var, tmp) + else + this_size = initial_size + end if + + if (present(n)) then + new_size = n + else + new_size = this_size + this_size/2 + 1 + end if + + allocate(var(new_size)) + + if (allocated(tmp)) then + this_size = min(size(tmp, 1), size(var, 1)) + var(:this_size) = tmp(:this_size) + end if +end subroutine resize +``` + +Finally, allocations do not initialize the array, the content of the uninitialized +array is most likely just the bytes of whatever was previously at the respective address. +The allocation supports initialization using the source attribute: + +```fortran +real(dp), allocatable :: arr(:) +allocate(arr(10), source=0.0_dp) +``` + +The ``source`` keyword supports scalar and array valued variables and constants. From 976fd62142093a88c65fa2f3c66c3f9f6acad190 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Tue, 1 Jun 2021 12:27:52 +0200 Subject: [PATCH 08/24] Expand chapter on arrays - add discussion of assumed-rank and assumed-size arrays --- learn/best_practices/arrays.md | 145 +++++++++++++++++++++++++-------- 1 file changed, 109 insertions(+), 36 deletions(-) diff --git a/learn/best_practices/arrays.md b/learn/best_practices/arrays.md index e1813c29c..8b5eb28d1 100644 --- a/learn/best_practices/arrays.md +++ b/learn/best_practices/arrays.md @@ -4,10 +4,19 @@ 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 *assumed-shape*): +Arrays are a central object in Fortran. The creation of dynamic sized arrays +is discussed in the [allocatable arrays arrays](./allocatable_arrays.html). -``` fortran +To pass arrays to procedures four ways are available + +1. *assumed-shape* arrays +4. *assumed-rank* arrays +2. *explicit-shape* arrays +3. *assumed-size* arrays + +The preferred way to pass arrays to procedures is as *assumed-shape* arrays + +```fortran subroutine f(r) real(dp), intent(out) :: r(:) integer :: n, i @@ -15,42 +24,71 @@ subroutine f(r) do i = 1, n r(i) = 1.0_dp / i**2 end do -end subroutine +end subroutine f ``` -2D arrays: +Higher dimensional arrays can be passed in a similar way. -``` fortran +```fortran subroutine g(A) real(dp), intent(in) :: A(:, :) ... -end subroutine +end subroutine g ``` -and call it like this: +The array is simply passed by -``` fortran +```fortran real(dp) :: r(5) call f(r) ``` -No array copying is done above. It has the following advantages: +In this case no array copy is done, which has the advantage that the shape and size +information is automatically passed along and checked at compile and optionally at +runtime. +Similarly, array strides can be passed without requiring a copy of the array but as +*assumed-shape* discriptor: -- 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. +```fortran +real(dp) :: r(10) +call f(r(1:10:2)) +call f(r(2:10:2)) +``` + +This should always be your default way of passing arrays in and out of subroutines. +Avoid passing arrays as whole slices, as it obfuscates the actual intent of the code: + +```fortran +real(dp) :: r(10) +call f(r(:)) +``` -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: +In case more general arrays should be passed to a procedure the *assumed-rank* +functionality introduced in the Fortran 2018 standard can be used + +```fortran +subroutine h(r) + real(dp), intent(in) :: r(..) + select rank(r) + rank(1) + ! ... + rank(2) + ! ... + end select +end subroutine h +``` + +The actual rank can be queried at runtime using the ``select rank`` construct. +This easily allows to create more generic functions that have to deal with +differnet array ranks. -- 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) +*Explicit-shape* arrays can be useful for returning data from functions. +Most of their functionality can be provided by *assumed-shape* and *assumed-rank* +arrays but they find frequent use for interfacing with C or in legacy Fortran +procedures, therefore they will be discussed briefly here. -To use *explicit-shape* arrays, do: +To use *explicit-shape* arrays, the dimension has to be passed explicitly as dummy +argument like in the example below ``` fortran subroutine f(n, r) @@ -63,7 +101,7 @@ subroutine f(n, r) end subroutine ``` -2D arrays: +For high-dimensional arrays additional indices have to be passed. ``` fortran subroutine g(m, n, A) @@ -73,15 +111,28 @@ subroutine g(m, n, A) end subroutine ``` -and call it like this: +The routines can be invoked by ``` fortran real(dp) :: r(5), s(3, 4) call f(size(r), r) -call f(size(s), s) +call g(size(s, 1), size(s, 2), s) +``` + +Note that the shape is not checked, therefore the following would be valid code +with will potentially yield incorrect results: + +``` +real(dp) :: s(3, 4) +call g(size(s), 1, s) ! s(12, 1) in g +call g(size(s, 2), size(s, 1), s) ! s(4, 3) in g ``` -In order to return an array from a function, do: +In this case the memory layout is preserved but the shape is changed. +Also, *explicit-shape* arrays require contiguous memory and will create temporary +arrays in case non-contiguous array strides are passed. + +To return an array from a function with *explicit-shape* use ``` fortran function f(n) result(r) @@ -94,27 +145,49 @@ function f(n) result(r) end function ``` -If you want to enforce/check the size of the arrays, put at the -beginning of the function: +Finally, there are *assumed-size* arrays, which provide the least compile and runtime +checking and can be found be found frequently in legacy code, they should be avoided +in favour of *assumed-shape* or *assumed-rank* arrays. +An *assumed-size* array dummy argument is identified by an asterisk as the last dimension, +this disables the usage of this array with many intrinsic functions, like ``size`` or +``shape``. -``` fortran +To check for the correct size and shape of an *assumed-shape* array the ``size`` and +``shape`` intrinsic functions can be used to query for those properties + +```fortran if (size(r) /= 4) error stop "Incorrect size of 'r'" +if (any(shape(r) /= [2, 2])) error stop "Incorrect shape of 'r'" ``` -To initialize an array, do: +Note that ``size`` returns the total size of all dimensions, to obtain the shape of +a specific dimension add it as second argument to the function. -``` fortran +Arrays can be initialized by using an array constructor + +```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). +The array constructor can be annoted with the type of the constructed array + +```fortran +real(dp) :: r(5) +r = [real(dp) :: 1, 2, 3, 4, 5] +``` + +Implicit do loops can be used inside an array constructor as well + +```fortran +integer :: i +real(dp) :: r(5) +r = [(real(i**2, dp), i = 1, size(r))] +``` In order for the array to start with different index than 1, do: -``` fortran +```fortran subroutine print_eigenvalues(kappa_min, lam) integer, intent(in) :: kappa_min real(dp), intent(in) :: lam(kappa_min:) @@ -123,5 +196,5 @@ subroutine print_eigenvalues(kappa_min, lam) do kappa = kappa_min, ubound(lam, 1) print *, kappa, lam(kappa) end do -end subroutine +end subroutine print_eigenvalues ``` From c358acea9d698af9ffefce9d42a33fe256db27b6 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Tue, 1 Jun 2021 12:44:32 +0200 Subject: [PATCH 09/24] Expand chapter on integer division --- learn/best_practices/integer_division.md | 30 ++++++++++++++++++------ 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/learn/best_practices/integer_division.md b/learn/best_practices/integer_division.md index 4377b976f..b3fb26f6c 100644 --- a/learn/best_practices/integer_division.md +++ b/learn/best_practices/integer_division.md @@ -4,13 +4,29 @@ title: Integer Division permalink: /learn/best_practices/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: +Fortran distinguishs between floating point and integer arithmetic. +It is important to note that division for integers is always using +integer arithmetic. +Consider the following example for integer division of an odd number: -``` fortran -(N - 1.0_dp)/N +```fortran +integer :: n +n = 3 +print *, n / 2 ! prints 1 +print *, n*(n + 1)/2 ! prints 6 +print *, n/2*(n + 1) ! prints 4 +n = -3 +print *, n / 2 ! prints -1 ``` -As long as one of the `/` operands is a float, a floating point division -will be used. +Be careful about whether you want to actually use integer arithmetic +in this context. If you want to use floating point arithmetic instead +make sure to cast to reals before using the division operator: + +```fortran +integer :: n +n = 3 +print *, real(n, dp) / 2 ! prints 1.5 +n = -3 +print *, real(n, dp) / 2 ! prints -1.5 +``` From 4114262cef63e097d0e8bed9a375377859e2c4c2 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Tue, 1 Jun 2021 12:53:38 +0200 Subject: [PATCH 10/24] Add definition of callback Co-authored-by: Jeremie Vandenplas --- learn/best_practices/callbacks.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/learn/best_practices/callbacks.md b/learn/best_practices/callbacks.md index 1a81f971d..2d4189ba0 100644 --- a/learn/best_practices/callbacks.md +++ b/learn/best_practices/callbacks.md @@ -3,7 +3,7 @@ layout: book title: Callbacks permalink: /learn/best_practices/callbacks --- - +A callback is a function that is passed as an argument to another fucntion. There are two ways to implement callbacks to be used like this: ``` fortran From 4b3de58b4156cb64c4ed656cf8dbc90d7f41e765 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Tue, 1 Jun 2021 12:54:17 +0200 Subject: [PATCH 11/24] Add intent in callback snippet Co-authored-by: Jeremie Vandenplas --- learn/best_practices/callbacks.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/learn/best_practices/callbacks.md b/learn/best_practices/callbacks.md index 2d4189ba0..1d10f3aad 100644 --- a/learn/best_practices/callbacks.md +++ b/learn/best_practices/callbacks.md @@ -9,7 +9,7 @@ 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 + real(dp), intent(in) :: a, k print *, simpson(f, 0._dp, pi) print *, simpson(f, 0._dp, 2*pi) From f14634277d34ea40ca5e928b561b96c541565ad7 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Wed, 2 Jun 2021 15:29:34 +0200 Subject: [PATCH 12/24] Rework multidimensional array section --- learn/best_practices/multidim_arrays.md | 105 +++++++++++------------- 1 file changed, 46 insertions(+), 59 deletions(-) diff --git a/learn/best_practices/multidim_arrays.md b/learn/best_practices/multidim_arrays.md index a981898db..d8f39049f 100644 --- a/learn/best_practices/multidim_arrays.md +++ b/learn/best_practices/multidim_arrays.md @@ -4,80 +4,67 @@ title: Multidimensional Arrays permalink: /learn/best_practices/multidim_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. +Multidimensional arrays are stored in column-major order. This means the +left-most (inner-most) index addresses elements contiguously. +From a practical point this means that the array slice ``V(:, 1)`` is +contiguous, while the stride between elements in the slice ``V(1, :)`` +is the dimension of the columns. This is important when passing array +slices to procedures which expect to work on contiguous data. -Example: +The locality of the memory is important to consider depending on +your application, usually when performing operations on a multidimensional +the sequential access should always advance in unity strides. -``` fortran -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: +In the following example the inverse distance between two sets of points +is evaluated, note that the points are stored contiguously in the arrays +``xyz1``/``xyz2``, while the inner-most loop is advancing the left-most +index of the matrix ``a``. -``` fortran -do i3 = 1, N3 - do i2 = 1, N2 - do i1 = 1, N1 - A(i1, i2, i3) +``` +subroutine coulomb_matrix(xyz1, xyz2, a) + real(dp), intent(in) :: xyz1(:, :) + real(dp), intent(in) :: xyz2(:, :) + real(dp), intent(out) :: a(:, :) + integer :: i, j + do i = 1, size(a, 2) + do j = 1, size(a, 1) + a(j, i) = 1.0_dp/norm2(xyz1(:, j) - xyz2(:, i)) end do end do -end do +end subroutine coulomb_matrix ``` -Associated array of vectors would then be most efficiently accessed as: +Another example would be the contraction of the third dimension of a rank +three array: -``` fortran -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: - -``` fortran -do i3 = 1, N3 - A(:, :, i3) +do i = 1, size(amat, 3) + do j = 1, size(amat, 2) + do k = 1, size(amat, 1) + cmat(k, j) = cmat(k, j) + amat(k, j, i) * bvec(i) + end do + end do end do ``` -Storing/accessing as above then accesses always contiguous blocks of -memory, directly adjacent to one another; no skips/strides. +Contiguous array slices can be used in array-bound remapping to allow usage +of higher rank arrays as lower rank arrays without requiring to reshape +and potentially create a temporary copy of the array. -When not sure, always rewrite (in your head) the algorithm to use -strides, for example the first loop would become: +For example this can be used to contract the third dimension of a rank +three array using a matrix-vector operation: -``` fortran -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 ``` +subroutine matmul312(amat, bvec, cmat) + real(dp), contiguous, intent(in) :: amat(:, :, :) + real(dp), intent(in) :: bvec(:) + real(dp), contiguous, intent(out) :: cmat(:, :) + real(dp), pointer :: aptr(:, :) + real(dp), pointer :: cptr(:) -the second loop would become: + aptr(1:size(amat, 1)*size(amat, 2), 1:size(amat, 3)) => amat + cptr(1:size(cmat)) => cmat -``` fortran -do i3 = 1, N3 - Ai3 = A(:, :, i3) - do i2 = 1, N2 - Ai3(:, i2) - end do -end do + cptr = matmul(aptr, bvec) +end subroutine matmul312 ``` - -And then make sure that all the strides are always on the left. Then it -will be fast. From b88565402a055554225b7f998bfbcef91f2f122c Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Wed, 2 Jun 2021 16:37:32 +0200 Subject: [PATCH 13/24] Rework floating point chapter --- learn/best_practices/floating_point.md | 79 +++++++++++++++++++------- 1 file changed, 59 insertions(+), 20 deletions(-) diff --git a/learn/best_practices/floating_point.md b/learn/best_practices/floating_point.md index ee555a4b3..2ecd9dbc0 100644 --- a/learn/best_practices/floating_point.md +++ b/learn/best_practices/floating_point.md @@ -4,36 +4,71 @@ title: Floating Point Numbers permalink: /learn/best_practices/floating_point --- -Somewhere create and export a parameter `dp`: +The default representation of floating point numbers is using single precision +(usually 32 bits / 4 bytes). For most application a higher precision is required. +For this purpose a custom kind parameter can be defined. +The recommended way of defining kind parameters is to use -``` fortran -integer, parameter:: dp=kind(0.d0) ! double precision +```fortran +integer, parameter :: dp = selected_real_kind(15) ``` -and declare floats as: +For many purposes it also suffices to directly infer the kind parameter from +a literal like here + +```fortran +integer, parameter :: dp = kind(0.0d0) +``` + +or to rename the imported kind parameter from the ``iso_fortran_env`` module + +```fortran +use, intrinsic :: iso_fortran_env, only : dp => real64 +``` + +For some insightful thoughts on kind parameters see +Doctor Fortran in it takes all KINDs. + +It is recommended to have a central module to define kind parameters and include +them with use as necessary. An example for such a module is given with -``` fortran -real(dp) :: a, b, c ``` +!> Numerical storage size parameters for real and integer values +module kind_parameter + implicit none + public + + !> Single precision real numbers, 6 digits, range 10⁻³⁷ to 10³⁷-1; 32 bits + integer, parameter :: sp = selected_real_kind(6, 37) + !> Double precision real numbers, 15 digits, range 10⁻³⁰⁷ to 10³⁰⁷-1; 64 bits + integer, parameter :: dp = selected_real_kind(15, 307) + !> Quadruple precision real numbers, 33 digits, range 10⁻⁴⁹³¹ to 10⁴⁹³¹-1; 128 bits + integer, parameter :: qp = selected_real_kind(33, 4931) -Always write all floating point constants with the `_dp` suffix: + !> Char length for integers, range -2⁷ to 2⁷-1; 8 bits + integer, parameter :: i1 = selected_int_kind(2) + !> Short length for integers, range -2¹⁵ to 2¹⁵-1; 16 bits + integer, parameter :: i2 = selected_int_kind(4) + !> Length of default integers, range -2³¹ to 2³¹-1; 32 bits + integer, parameter :: i4 = selected_int_kind(9) + !> Long length for integers, range -2⁶³ to 2⁶³-1; 64 bits + integer, parameter :: i8 = selected_int_kind(18) -``` fortran -1.0_dp, 3.5_dp, 1.34e8_dp +end module kind_parameter ``` -and never any other way (see also the gotcha -`floating_point_numbers_gotcha`). Omitting the dot in the literal -constant is also incorrect. +Floating point constants should always be declared including a kind parameter suffix: -To print floating point double precision numbers without losing -precision, use the `(es23.16)` format (see -). +```fortran +real(dp) :: a, b, c +a = 1.0_dp +b = 3.5_dp +c = 1.34e8_dp +``` -It is safe to assign integers to floating point numbers without losing -accuracy: +It is safe to assign integers to floating point numbers without losing accuracy: -``` fortran +```fortran real(dp) :: a a = 3 ``` @@ -42,7 +77,11 @@ 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: -``` fortran +```fortran real(dp) :: a -a = real(1, dp) / 2 ! 'a' is equal to 0.5_dp +a = real(1, dp) / 2 ! 'a' is equal to 0.5_dp ``` + +To print floating point numbers without losing precision use the unlimited +format specificer ``(g0)`` or the exponential representation ``(es24.16e3)``, +which will give you 17 significant digits of printout. From 1fb3b4455dcde6c0da3973061ad0731bf42ad270 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Wed, 2 Jun 2021 18:22:40 +0200 Subject: [PATCH 14/24] Remove parallel_programming from best practices -> deserves a separate mini-book with OpenMP, OpenACC, MPI, Coarrays --- _data/learning.yml | 1 - learn/best_practices/parallel_programming.md | 51 -------------------- 2 files changed, 52 deletions(-) delete mode 100644 learn/best_practices/parallel_programming.md diff --git a/_data/learning.yml b/_data/learning.yml index 51a9df001..d46fbe914 100644 --- a/_data/learning.yml +++ b/_data/learning.yml @@ -71,7 +71,6 @@ books: - 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 # diff --git a/learn/best_practices/parallel_programming.md b/learn/best_practices/parallel_programming.md deleted file mode 100644 index ffce72588..000000000 --- a/learn/best_practices/parallel_programming.md +++ /dev/null @@ -1,51 +0,0 @@ ---- -layout: book -title: Parallel programming -permalink: /learn/best_practices/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 : - -``` fortran -program test_openmp - !$ 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. From 7a9fcbc17b96ffabd1c4a391fda520386160f951 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Wed, 2 Jun 2021 18:24:12 +0200 Subject: [PATCH 15/24] Remove Python interfacing chapter from best practices -> deserves its own mini-book with more detailed tutorials --- _data/learning.yml | 1 - learn/best_practices/python_interfacing.md | 70 ---------------------- 2 files changed, 71 deletions(-) delete mode 100644 learn/best_practices/python_interfacing.md diff --git a/_data/learning.yml b/_data/learning.yml index d46fbe914..94851c185 100644 --- a/_data/learning.yml +++ b/_data/learning.yml @@ -68,7 +68,6 @@ books: - 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 diff --git a/learn/best_practices/python_interfacing.md b/learn/best_practices/python_interfacing.md deleted file mode 100644 index 8abe8f0f8..000000000 --- a/learn/best_practices/python_interfacing.md +++ /dev/null @@ -1,70 +0,0 @@ ---- -layout: book -title: Interfacing with Python -permalink: /learn/best_practices/python_interfacing ---- - -Using Cython ------------- - -To wrap Fortran code in Python, export it to C first (see above) and -then write this Cython code: - -``` 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: - -``` fortran -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: - -``` 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))) -``` From f5978fbe7cfe6807bad68da441a8e9fb6a3c0c8d Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Wed, 2 Jun 2021 18:28:44 +0200 Subject: [PATCH 16/24] Remove Fortran-C interfacing from best practice -> deserves its own minibook (iso_c_binding + iso_fortran_binding) --- _data/learning.yml | 1 - learn/best_practices/c_interfacing.md | 60 --------------------------- 2 files changed, 61 deletions(-) delete mode 100644 learn/best_practices/c_interfacing.md diff --git a/_data/learning.yml b/_data/learning.yml index 94851c185..46abe662f 100644 --- a/_data/learning.yml +++ b/_data/learning.yml @@ -67,7 +67,6 @@ books: - 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/callbacks - link: /learn/best_practices/type_casting diff --git a/learn/best_practices/c_interfacing.md b/learn/best_practices/c_interfacing.md deleted file mode 100644 index 6617eec19..000000000 --- a/learn/best_practices/c_interfacing.md +++ /dev/null @@ -1,60 +0,0 @@ ---- -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. From 81aa4d5d4476005a9ca36dc9bca513ac0e798595 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Wed, 2 Jun 2021 18:35:13 +0200 Subject: [PATCH 17/24] Apply suggestions from code review Co-authored-by: Jeremie Vandenplas --- learn/best_practices/callbacks.md | 2 +- learn/best_practices/element_operations.md | 16 ++++++++-------- learn/best_practices/file_io.md | 16 ++++++++-------- learn/best_practices/modules_programs.md | 4 ++-- learn/best_practices/type_casting.md | 8 ++++---- 5 files changed, 23 insertions(+), 23 deletions(-) diff --git a/learn/best_practices/callbacks.md b/learn/best_practices/callbacks.md index 1d10f3aad..063a15330 100644 --- a/learn/best_practices/callbacks.md +++ b/learn/best_practices/callbacks.md @@ -50,7 +50,7 @@ end function end module ``` -The other approach since f2003 is to first define a new type for our +Since the Fortran 2003 Standard, the other approach is to first define a new type for our callback, and then use `procedure(func)` as the type of the dummy argument: diff --git a/learn/best_practices/element_operations.md b/learn/best_practices/element_operations.md index 13b12de74..395181068 100644 --- a/learn/best_practices/element_operations.md +++ b/learn/best_practices/element_operations.md @@ -4,9 +4,9 @@ title: Element-wise Operations on Arrays Using Subroutines/Functions permalink: /learn/best_practices/element_operations --- -There are three approaches: +There are three approaches to perform element-wise operations on arrays when using subroutines and functions: -- `elemental` subroutines +- `elemental` procedures - *explicit-shape* arrays - implementing the operation for vectors and write simple wrapper subroutines (that use `reshape` internally) for each array shape @@ -49,14 +49,14 @@ 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). +The `elemental` keyword implies the `pure` keyword, so the procedure +must be pure. It results that `elemental procedures` can only use `pure` procedures 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 +If the elemental procedure algorithm can be made faster using array +operations inside, or if for some reasons 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: +can make `nroot` operate on a vector and write a simple wrapper for +other array shapes, e.g.: ``` fortran function nroot(n, x) result(y) diff --git a/learn/best_practices/file_io.md b/learn/best_practices/file_io.md index 1e15877ff..cd098f79c 100644 --- a/learn/best_practices/file_io.md +++ b/learn/best_practices/file_io.md @@ -8,21 +8,21 @@ To read from a file: ``` fortran integer :: u -open(newunit=u, file="log.txt", status="old") +open(newunit=u, file="log.txt", status="old", action='read') read(u, *) a, b close(u) ``` -Write to a file: +Write to a file as follows: ``` fortran integer :: u -open(newunit=u, file="log.txt", status="replace") +open(newunit=u, file="log.txt", status="replace", action='write') write(u, *) a, b close(u) ``` -To append to an existing file: +It is possible to append to an existing file as follows: ``` fortran integer :: u @@ -31,8 +31,8 @@ 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: +The `newunit` keyword argument to `open` is a Fortran 2008 standard feature. Therefore for +older compilers that do not suport it, just replace `open(newunit=u, ...)` by: ``` fortran open(newunit(u), ...) @@ -44,7 +44,7 @@ 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 + 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 @@ -54,6 +54,6 @@ integer function newunit(unit) result(n) return end if end do - call stop_error("newunit ERROR: available unit not found.") + error stop 'newunit ERROR: available unit not found.' end function ``` diff --git a/learn/best_practices/modules_programs.md b/learn/best_practices/modules_programs.md index cab56cb1b..ae59d863c 100644 --- a/learn/best_practices/modules_programs.md +++ b/learn/best_practices/modules_programs.md @@ -29,8 +29,8 @@ 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 +need to repeat it within each procedure). By keeping the `private` empty, all your +procedures/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: diff --git a/learn/best_practices/type_casting.md b/learn/best_practices/type_casting.md index 80be54185..c53544240 100644 --- a/learn/best_practices/type_casting.md +++ b/learn/best_practices/type_casting.md @@ -4,7 +4,7 @@ title: Type Casting in Callbacks permalink: /learn/best_practices/type_casting --- -There are essentially five different ways to do that, each with its own +There are essentially five different ways to do type casting, each with its own advantages and disadvantages. The methods I, II and V can be used both in C and Fortran. The methods @@ -14,7 +14,7 @@ should not be used. Work Arrays ----------- -Pass a "work array" or two which are packed with everything needed by +Pass a "work array" which is packed with everything needed by the caller and unpacked by the called routine. This is the old way -- e.g., how LAPACK does it. @@ -83,9 +83,9 @@ end module General Structure ----------------- -Define general structure or two which encompass the variations you +Define a general structure 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 +single structure type 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. From e17cba7a825f5d4830818cc62d75567e16fb1a56 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Wed, 2 Jun 2021 18:44:58 +0200 Subject: [PATCH 18/24] Fix some typos in allocatable arrays section --- learn/best_practices/allocatable_arrays.md | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/learn/best_practices/allocatable_arrays.md b/learn/best_practices/allocatable_arrays.md index 26591289f..7bf5e9296 100644 --- a/learn/best_practices/allocatable_arrays.md +++ b/learn/best_practices/allocatable_arrays.md @@ -4,9 +4,9 @@ title: Allocatable Arrays permalink: /learn/best_practices/allocatable_arrays --- -The ``allocatable`` attribute provides safe way to for memory handling. +The ``allocatable`` attribute provides safe way for memory handling. In comparison to variables with ``pointer`` attribute the memory is managed -automatically and will go deallocated automatically once the variable goes +automatically and will be deallocated automatically once the variable goes out-of-scope. Using ``allocatable`` variables removes the possibility to create memory leaks in an application. @@ -33,7 +33,7 @@ end subroutine show_arr To allocate variables inside a procedure the dummy argument has to carry the ``allocatable`` attribute. Using it in combination with ``intent(out)`` -will deallocate previous allocations before entring the procedure: +will deallocate previous allocations before entering the procedure: ```fortran subroutine foo(lam) @@ -87,7 +87,8 @@ the argument will not be present. The ``allocatable`` attribute is not limited t arrays and can also be associated with scalars, which can be useful in combination with ``optional`` dummy arguments. -Allocations can be moved between different arrays with ``allocatable`` attribute. +Allocations can be moved between different arrays with ``allocatable`` attribute +using the ``move_alloc`` intrinsic subroutine. ```fortran subroutine resize(var, n) From e73677a1d9bd6f41d92164562708cf6ad5191ff2 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Wed, 2 Jun 2021 18:51:13 +0200 Subject: [PATCH 19/24] Add correct syntax highlighting --- learn/best_practices/arrays.md | 2 +- learn/best_practices/multidim_arrays.md | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/learn/best_practices/arrays.md b/learn/best_practices/arrays.md index 8b5eb28d1..05d302154 100644 --- a/learn/best_practices/arrays.md +++ b/learn/best_practices/arrays.md @@ -122,7 +122,7 @@ call g(size(s, 1), size(s, 2), s) Note that the shape is not checked, therefore the following would be valid code with will potentially yield incorrect results: -``` +```fortran real(dp) :: s(3, 4) call g(size(s), 1, s) ! s(12, 1) in g call g(size(s, 2), size(s, 1), s) ! s(4, 3) in g diff --git a/learn/best_practices/multidim_arrays.md b/learn/best_practices/multidim_arrays.md index d8f39049f..0d64ee0cf 100644 --- a/learn/best_practices/multidim_arrays.md +++ b/learn/best_practices/multidim_arrays.md @@ -20,7 +20,7 @@ is evaluated, note that the points are stored contiguously in the arrays ``xyz1``/``xyz2``, while the inner-most loop is advancing the left-most index of the matrix ``a``. -``` +```fortran subroutine coulomb_matrix(xyz1, xyz2, a) real(dp), intent(in) :: xyz1(:, :) real(dp), intent(in) :: xyz2(:, :) @@ -37,7 +37,7 @@ end subroutine coulomb_matrix Another example would be the contraction of the third dimension of a rank three array: -``` +```fortran do i = 1, size(amat, 3) do j = 1, size(amat, 2) do k = 1, size(amat, 1) @@ -54,7 +54,7 @@ and potentially create a temporary copy of the array. For example this can be used to contract the third dimension of a rank three array using a matrix-vector operation: -``` +```fortran subroutine matmul312(amat, bvec, cmat) real(dp), contiguous, intent(in) :: amat(:, :, :) real(dp), intent(in) :: bvec(:) From 158a970304371007e5f3e10c956809e83dbfead9 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Wed, 2 Jun 2021 19:26:12 +0200 Subject: [PATCH 20/24] Rework file IO section - remove mention of impure newunit function hack --- learn/best_practices/file_io.md | 122 +++++++++++++++++++++----------- 1 file changed, 82 insertions(+), 40 deletions(-) diff --git a/learn/best_practices/file_io.md b/learn/best_practices/file_io.md index cd098f79c..00c893160 100644 --- a/learn/best_practices/file_io.md +++ b/learn/best_practices/file_io.md @@ -4,56 +4,98 @@ title: File Input/Output permalink: /learn/best_practices/file_io --- -To read from a file: +In Fortran files are managed by unit identifiers. Interaction with the filesystem +mainly happens through the ``open`` and ``inquire`` built-in procedures. +Generally, the workflow is to open a file to a unit identifier, read and/or write +to it and close it again. -``` fortran -integer :: u -open(newunit=u, file="log.txt", status="old", action='read') -read(u, *) a, b -close(u) +```fortran +integer :: io +open(newunit=io, file="log.txt") +! ... +close(io) ``` -Write to a file as follows: +By default the file will be created if it is not existing already and opened for +both reading and writing. Writing to an existing file will start in the first +record (line) and therefore overwrite the file by default. -``` fortran -integer :: u -open(newunit=u, file="log.txt", status="replace", action='write') -write(u, *) a, b -close(u) +To create a read-only access to a file the ``status`` and ``action`` have to be +specified with + +```fortran +integer :: io +open(newunit=io, file="log.txt", status="old", action="read") +read(io, *) a, b +close(io) +``` + +In case the file is not present a runtime error will occur. To check for the existence +of a file prior to opening it the ``inquire`` function can be used + +```fortran +logical :: exists +inquire(file="log.txt", exist=exists) +if (exists) then + ! ... +end if ``` -It is possible to append to an existing file as follows: +Alternatively, the ``open`` procedure can return an optional *iostat* and *iomsg*: -``` fortran -integer :: u -open(newunit=u, file="log.txt", position="append", status="old") -write(u, *) N, V(N) -close(u) +```fortran +integer :: io, stat +character(len=512) :: msg +open(newunit=io, file="log.txt", status="old", action="read", & + iostat=stat, iomsg=msg) +if (stat /= 0) then + print *, trim(msg) +end if ``` -The `newunit` keyword argument to `open` is a Fortran 2008 standard feature. Therefore for -older compilers that do not suport it, just replace `open(newunit=u, ...)` by: +Note that *iomsg* requires a fixed-length character variable with sufficent storage +size to hold the error message. + +Similarly, writing to a file happens by using the *status* and *action* keyword. +To create a new file use -``` fortran -open(newunit(u), ...) +```fortran +integer :: io +open(newunit=io, file="log.txt", status="new", action="write") +write(io, *) a, b +close(io) ``` -where the `newunit` function is defined by: - -``` fortran -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 - error stop 'newunit ERROR: available unit not found.' -end function +Alternatively, ``status="replace"`` can be used to overwrite an existing file, +it is highly recommended to first check for the existence of a file before deciding +on the *status* to use. +To append to an output file the *position* keyword can be specified explicitly with + +```fortran +integer :: io +open(newunit=io, file="log.txt", position="append", & + & status="old", action="write") +write(io, *) size(v) +write(io, *) v(:) +close(io) ``` + +To reset the position in a file the built-in procedures ``rewind`` and ``backspace`` +can be used. ``rewind`` will reset to the first record (line), while ``backspace`` will +return to the previous record (line). + +Finally, to delete a file the file has to be opened and can be deleted after closing +with + +```fortran +logical :: exists +integer :: io, stat +inquire(file="log.txt", exist=exists) +if (exists) then + open(file="log.txt", newunit=io, iostat=stat) + if (stat == 0) close(io, status="delete", iostat=stat) +end if +``` + +A useful IO feature are scratch files, which can be opened with ``status="scratch"``. +They are automatically deleted after closing the unit identifier. From a8da64c7e06650525415e649b5eef5afe5b77117 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Thu, 3 Jun 2021 11:16:45 +0200 Subject: [PATCH 21/24] Rework module and program chapter --- learn/best_practices/modules_programs.md | 160 ++++++++++++++++++----- 1 file changed, 124 insertions(+), 36 deletions(-) diff --git a/learn/best_practices/modules_programs.md b/learn/best_practices/modules_programs.md index ae59d863c..a8a702ee5 100644 --- a/learn/best_practices/modules_programs.md +++ b/learn/best_practices/modules_programs.md @@ -4,57 +4,145 @@ title: Modules and Programs permalink: /learn/best_practices/modules_programs --- -Only use modules and programs. Always setup a module in the following -way: +Modules are the preferred way create modern Fortran libraries and applications. +As a convention, one source file should always contain only one module, while +the module name should match the filepath to allow easy navigation in larger +projects. It is also recommended to prefix module names with the library name +to avoid nameclashes when used as dependency in other projects. -``` fortran -module ode1d - use types, only: dp, pi - use utils, only: stop_error +An example for such a module file is given here + +``` +!> Interface to TOML processing library. +!> +!> ... +module fpm_toml + use fpm_error, only : error_t, fatal_error, file_not_found_error + use fpm_strings, only : string_t + use tomlf, only : toml_table, toml_array, toml_key, toml_stat, get_value, & + & set_value, toml_parse, toml_error, new_table, add_table, add_array, & + & toml_serializer, len 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 + + public :: read_package_file + public :: toml_table, toml_array, toml_key, toml_stat, get_value, set_value + public :: new_table, add_table, add_array, len + public :: toml_error, toml_serializer, toml_parse contains -subroutine get_val(...) - ... -end subroutine -... + !> Process the configuration file to a TOML data structure + subroutine read_package_file(table, manifest, error) + !> TOML data structure + type(toml_table), allocatable, intent(out) :: table + !> Name of the package configuration file + character(len=*), intent(in) :: manifest + !> Error status of the operation + type(error_t), allocatable, intent(out) :: error + ! ... + end subroutine read_package_file -end module +end module fpm_toml ``` -The `implicit none` statement works for the whole module (so you don't -need to repeat it within each procedure). By keeping the `private` empty, all your -procedures/data types will be private to the module by default. Then -you export things by putting it into the `public` clause. +There are a few things in this example module to highlight. First, every module +starts with comments documenting the purpose and content of the module, +similarly, every procedure starts with a comment briefly describing its +purse and the intent of the dummy arguments. Documentation is one of the most +important parts of creating long-living software, regardless of language. + +Second, imports (*use*) and exports (*public*) are explicitly given, this +allows on a glance at the module source to check the used and available +procedures, constants and derived types. The imports are usually limited +to the module scope rather than reimported in every procedure or interface +scope. Similarly, exports are made explicitly by adding a *private* statement +one a single line and explicitly listing all exported symbols in *public* +statements. + +Finally, the `implicit none` statement works for the whole module and there +is no need to repeat it within each procedure. -Setup programs in the following way: +Variables inside a module are static (*implicitly saved*), it is highly +recommended to limit the usage of module variables to constant expressions, +like parameters or enumerators only or export them as *protected* rather +than *public*. -``` fortran -program uranium - use fmesh, only: mesh_exp - use utils, only: stop_error, dp - use dft, only: atom +Submodules can be used to break long dependency chains and shorten recompilation +cascades in Fortran programs. They also offer the possibility to provide specialized +and optimized implementations without requiring the use of preprocessor. + +An example from the Fortran standard library is the quadrature module, which +only defines interfaces to module procedures, but no implementations + +```fortran +!> Numerical integration +!> +!> ... +module stdlib_quadrature + use stdlib_kinds, only: sp, dp, qp implicit none + private + + public :: trapz + ! ... - 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 + !> Integrates sampled values using trapezoidal rule + interface trapz + pure module function trapz_dx_dp(y, dx) result(integral) + real(dp), intent(in) :: y(:) + real(dp), intent(in) :: dx + real(dp) :: integral + end function trapz_dx_dp + module function trapz_x_dp(y, x) result(integral) + real(dp), intent(in) :: y(:) + real(dp), intent(in) :: x(:) + real(dp) :: integral + end function trapz_x_dp + end interface trapz + + ! ... +end module stdlib_quadrature ``` -Notice the "explicit imports" (using Python terminology) in the `use` -statements. You can also use "implicit imports" like: +While the implementation is provided in separate submodules like the one for the +trapezoidal integration rule given here. + +```fortran +!> Actual implementation of the trapezoidal integration rule +!> +!> ... +submodule (stdlib_quadrature) stdlib_quadrature_trapz + use stdlib_error, only: check + implicit none + +contains + + pure module function trapz_dx_dp(y, dx) result(integral) + real(dp), intent(in) :: y(:) + real(dp), intent(in) :: dx + real(dp) :: integral + integer :: n -``` fortran -use fmesh + n = size(y) + select case (n) + case (0:1) + integral = 0.0_dp + case (2) + integral = 0.5_dp*dx*(y(1) + y(2)) + case default + integral = dx*(sum(y(2:n-1)) + 0.5_dp*(y(1) + y(n))) + end select + end function trapz_dx_dp + + ! ... +end submodule stdlib_quadrature_trapz ``` -But just like in Python, this should be avoided ("explicit is better -than implicit") in most cases. +Note that the module procedures do not have to be implemented in the same submodule. +Several submodules can be used to reduce the compilation load for huge modules. + +Finally, when setting up a program, it is recommended to keep the actual implementations +in the program body at minimum. Reusing implementations from modules allows to write +reusable code and focus in the program unit only on translating user input to the +respective library functions and objects. From 5cf80cc2ecb95001247d47e89b53910068aafbeb Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Thu, 3 Jun 2021 11:27:35 +0200 Subject: [PATCH 22/24] Update authors and introduction for best practice minibook --- learn/best_practices/element_operations.md | 2 +- learn/best_practices/index.md | 21 ++++++--------------- 2 files changed, 7 insertions(+), 16 deletions(-) diff --git a/learn/best_practices/element_operations.md b/learn/best_practices/element_operations.md index 395181068..241657412 100644 --- a/learn/best_practices/element_operations.md +++ b/learn/best_practices/element_operations.md @@ -1,6 +1,6 @@ --- layout: book -title: Element-wise Operations on Arrays Using Subroutines/Functions +title: Element-wise Operations on Arrays permalink: /learn/best_practices/element_operations --- diff --git a/learn/best_practices/index.md b/learn/best_practices/index.md index 40f2b8055..e087f567a 100644 --- a/learn/best_practices/index.md +++ b/learn/best_practices/index.md @@ -2,20 +2,11 @@ layout: book title: Fortran Best Practices permalink: /learn/best_practices -author: Ondřej Čertík, John Pask, Jed Brown, Matthew Emmett, Juan Luis Cano Rodríguez, Neil Carlson, Andrea Vigliotti, Pierre Haessig +author: Ondřej Čertík, John Pask, Jed Brown, Matthew Emmett, Juan Luis Cano Rodríguez, Neil Carlson, Andrea Vigliotti, Pierre Haessig, Vincent Magnin, Sebastian Ehlert, Jeremie Vandenplas --- -This mini-book 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: - - -Language features are explained at: - - -The best resource is a recent Fortran standard, for example the Fortran -2003 standard: +This mini-book collects a modern canonical way of doing things in Fortran. +It serves as a style guide and best practice recommendation for popular topics +and common tasks. Generally, a cononical solution or pattern is presented and +discussed. It is meant for programmers basic familarity of the Fortran syntax +and programming in general. From 5a86e10ad831fb5c8c7fe5810e30b6ae518da163 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Thu, 3 Jun 2021 11:30:54 +0200 Subject: [PATCH 23/24] Fix typo --- learn/best_practices/index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/learn/best_practices/index.md b/learn/best_practices/index.md index e087f567a..2abed3bce 100644 --- a/learn/best_practices/index.md +++ b/learn/best_practices/index.md @@ -8,5 +8,5 @@ author: Ondřej Čertík, John Pask, Jed Brown, Matthew Emmett, Juan Luis Cano R This mini-book collects a modern canonical way of doing things in Fortran. It serves as a style guide and best practice recommendation for popular topics and common tasks. Generally, a cononical solution or pattern is presented and -discussed. It is meant for programmers basic familarity of the Fortran syntax +discussed. It is meant for programmers with basic familarity of the Fortran syntax and programming in general. From c1e492a5b307bda9dbb3dbedf3bd19914c003e9e Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Thu, 3 Jun 2021 11:48:10 +0200 Subject: [PATCH 24/24] Rework callback section in best practise minibook --- learn/best_practices/callbacks.md | 138 ++++++++++++------------------ 1 file changed, 57 insertions(+), 81 deletions(-) diff --git a/learn/best_practices/callbacks.md b/learn/best_practices/callbacks.md index 063a15330..be7121e35 100644 --- a/learn/best_practices/callbacks.md +++ b/learn/best_practices/callbacks.md @@ -3,115 +3,91 @@ layout: book title: Callbacks permalink: /learn/best_practices/callbacks --- -A callback is a function that is passed as an argument to another fucntion. -There are two ways to implement callbacks to be used like this: -``` fortran -subroutine foo(a, k) - use integrals, only: simpson - real(dp), intent(in) :: a, k - print *, simpson(f, 0._dp, pi) - print *, simpson(f, 0._dp, 2*pi) +A callback is a function that is passed as an argument to another function. -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: +The preferred way of creating such a callback is to provide an *abstract interface* +declaring the signature of the callback. This allows to use compile time checks +for the passed callback. -``` fortran +```fortran module integrals use types, only: dp implicit none private - public simpson + public :: simpson, integratable_function -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 + abstract interface + function integratable_function(x) result(func) + import :: dp + real(dp), intent(in) :: x + real(dp) :: func end function end interface - s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b)) -end function -end module +contains + + function simpson(f, a, b) result(s) + real(dp), intent(in) :: a, b + procedure(integratable_function) :: f + real(dp) :: s + + s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b)) + end function simpson + +end module integrals ``` -Since the Fortran 2003 Standard, the other approach is to first define a new type for our -callback, and then use `procedure(func)` as the type of the dummy -argument: +The function can than be used with a callback by importing the module +as shown in the following example -``` fortran -module integrals +```fortran +module demo_functions use types, only: dp implicit none private - public simpson + public :: test_integral 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 + subroutine test_integral(a, k) + real(dp), intent(in) :: a, k + + print *, simpson(f, 0._dp, pi) + print *, simpson(f, 0._dp, 2*pi) + contains + + function f(x) result(y) + real(dp), intent(in) :: x + real(dp) :: y + y = a*sin(k*x) + end function f + end subroutine test_integral -end module +end module demo_functions ``` -The new type can also be defined outside of the function (and reused), -like: +Exporting the abstract interface allows to create procedure pointers with the +correct signature and also to extend the callback further like shown here -``` fortran -module integrals +```fortran +module demo_integrals use types, only: dp + use integrals, only: simpson, integratable_function 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 + public :: simpson2, integratable_function 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 + function simpson2(f, a, b) result(s) + real(dp), intent(in) :: a, b + procedure(integratable_function) :: f + real(dp) :: s + real(dp) :: mid + mid = (a + b)/2 + s = simpson(f, a, mid) + simpson(f, mid, b) + end function simpson2 + +end module demo_integrals ```