MPI with mpi_f08

MPI defines three ways to call it from Fortran (MPI 5.0 section 19.1):

The recommendation is simple: use mpi_f08 in new code. A plain use mpi_f08 is acceptable because every MPI name starts with MPI_ and a program must not declare names beginning with that prefix (section 2.6.2). Keep the use mpi_f08 statement in the one module or program unit that talks to MPI, as described in Intrinsic Modules and iso_fortran_env: it is a vendor-supplied module, not an intrinsic one, so it is not portable across installations and its .mod files are tied to the compiler that built the MPI library. mpi_f08 names are written with a capital letter after the prefix in the specification (MPI_Comm_rank), but Fortran is case insensitive and mpi_comm_rank works in all three bindings.

Two consequences of the strong typing matter from the start. First, an integer cannot stand for a handle, and the compiler says so. Example 7 passes an integer where MPI_Comm_rank wants a communicator:

Real gfortran 16.2 diagnostic through mpifort:

Second, the check stops at the handle. The buffer arguments are type(*), so the compiler cannot compare the buffer with the MPI datatype that describes it: describing a real(dp) array as MPI_REAL is a runtime error that nothing reports. Describing data: kinds and derived datatypes shows the consequence.

Source: MPI 5.0 report: 2.5.1 Opaque objects, 2.6.2 Fortran Binding Issues, 19.1 Support for Fortran, 19.1.2 Fortran Support Through the mpi_f08 Module and Appendix A.1.3 Prototype Definitions, J3/23-007r1: 7.3.2.2 TYPE type specifier (assumed type) and 8.5.8.7 Assumed-rank entity

Every MPI program begins with MPI_Init and ends with MPI_Finalize. Between them, a communicator (a group of ranks plus a private communication context) names who is talking; MPI_COMM_WORLD is the one containing all the ranks started together. MPI_Comm_rank returns the caller’s number, from 0, and MPI_Comm_size the number of ranks. With mpi_f08 the ierror argument is optional, so the calls are subroutines with the results in the last argument: there is no ierr variable to declare or forget to check.

Built with mpifort 01_hello.f90 and run with mpirun -n 4 ./a.out:

The four lines are written by four different processes to the same terminal, so their order is not defined; here they happened to come out sorted. The program also shows that this MPICH implements MPI 4.1 although MPI 5.0 is the current specification: MPI_Get_version reports the version the library implements, which can lag the published standard. MPI_Get_library_version returns a free-form string identifying the implementation and is the one to print in bug reports. Started without mpirun the same executable runs as a single rank (rank 0 of 1), which the fpm example below relies on.

Rules worth knowing (MPI 5.0 chapter 11):

Source: MPI 5.0 report: Chapter 11 Process Initialization, Creation, and Management, 19.1.2 Fortran Support Through the mpi_f08 Module and Chapter 20 Application Binary Interface

In mpi_f08 every MPI object is represented by a derived type whose only component is an integer, MPI_VAL, so a handle can be assigned, passed and compared but not confused with an integer or with a handle of another kind. The operators == and /= are overloaded for each handle type, which is how a communicator is tested against MPI_COMM_NULL. The next program splits MPI_COMM_WORLD into two groups, compares handles and frees one.

Real output with mpirun -n 4 (the lines of the four ranks arrive in any order):

Points to remember:

Source: MPI 5.0 report: 2.5.1 Opaque objects, 7.4 Communicator Management and 19.3.4 Transfer of Handles

MPI_Send and MPI_Recv move a buffer between two ranks of a communicator. The arguments are the buffer, the count of elements, the MPI datatype, the peer rank, a tag that lets a receiver select among messages, the communicator and, for a receive, a type(MPI_Status) that records where the message came from. The Fortran buffer is passed as written; because type(*), dimension(..) accepts any array or scalar, no loc, no flattening and no explicit shape is needed.

The messages are matched by the triple (source, tag, communicator), and the messages between a fixed pair of ranks do not overtake each other. A receive that names a smaller count than the message is an error (MPI_ERR_TRUNCATE); a larger count is allowed and MPI_Get_count tells how much arrived. That is the basis of the usual idiom for a message of unknown size: probe, ask the count, allocate, receive.

Real output with mpirun -n 4:

Source: MPI 5.0 report: Chapter 3 Point-to-Point Communication

Rank 0 sends a different number of values to each other rank, and the receivers size their allocatable array (see Allocatable Arrays and Dynamic Memory) from the probed count. MPI_STATUS_IGNORE says that the status is not wanted; it is a special constant, not a variable, and has the matching MPI_STATUSES_IGNORE for arrays of statuses. Useful members of type(MPI_Status) are MPI_SOURCE, MPI_TAG and MPI_ERROR, as the ring exchange shows:

MPI_Sendrecv posts the send and the receive together and is the safe way to write a ring or a halo swap: the pattern "everybody MPI_Sends to the right, then everybody MPI_Recvs from the left" works for small messages because they are buffered, and deadlocks for large ones, because MPI_Send may not return until the receiver has started the matching receive. Which of the two happens depends on the message size and the implementation, so a program that is correct only for small messages is not correct. Wildcards MPI_ANY_SOURCE and MPI_ANY_TAG exist for receives; MPI_PROC_NULL as a peer turns a send or receive into a no-operation, which removes the special cases at the ends of a decomposition (the heatflow example below uses it). The send modes are blocking (MPI_Send), synchronous (MPI_Ssend), buffered (MPI_Bsend) and ready (MPI_Rsend); MPI_Send is the one to start from.

Source: MPI 5.0 report: Chapter 3 Point-to-Point Communication

A collective operation is called by every rank of the communicator, with consistent arguments, and the library chooses the algorithm. They replace long hand-written send and receive loops and are usually faster. The common ones are:

The example distributes a root’s array, reduces it two ways and gathers the partial sums.

Real output with mpirun -n 4 (lines from different ranks interleave in any order):

Things that surprise newcomers:

Source: MPI 5.0 report: Chapter 6 Collective Communication

MPI_Isend and MPI_Irecv return immediately with an MPI_Request; the transfer proceeds while the program does other work, and MPI_Wait, MPI_Waitall or MPI_Test completes it. Until completion the buffer belongs to MPI: it must not be modified (send) or read (receive). The cost of the freedom is that the Fortran compiler does not know any of this.

Source: MPI 5.0 report: Chapter 3.7 Nonblocking Communication

The buffers are declared with the asynchronous attribute (J3/23-007r1 8.5.4). It is the same attribute that marks variables involved in asynchronous input/output (see Derived-Type and Asynchronous I/O; the book chapter on it is MRCB section 8.11), and it tells the compiler that the variable "might be associated with a pending input/output storage sequence" and may change at any time, which disables code motion across calls. The problem it solves is the one in the next two programs. They are identical except for the attribute, and both send three values to the right neighbour and receive three from the left, using persistent requests (MPI_Recv_init, MPI_Send_init, then MPI_Startall and MPI_Waitall on every step). The counts are of kind MPI_COUNT_KIND, the argument kind that selects the large-count variants of the routines through the generic names (MPI 5.0 section 19.2).

Without optimization (-O0) this prints what one expects. Compiled the way mpifort does it by default (-O2) and run with mpirun -n 4:

Source: MPI 5.0 report: Chapter 3.7 Nonblocking Communication, 19.1.8 Additional Support for Fortran Register-Memory-Synchronization and 19.1.17 Problems with Code Movement and Register Optimization, J3/23-007r1: 8.5.4 ASYNCHRONOUS attribute

The values never arrive in the program’s view of in: the compiler saw that in is not modified between its initialization and the print, because the call that completes the receive has no visible connection to it, and kept it in a register or folded the value. The program is wrong, not the library, and no compiler warning exists. Adding the attribute to both buffers is the whole fix:

With that single change (file 14b_persistent_async.f90) the same build prints:

The rules, from MPI 5.0 section 19.1.17 and J3/23-007r1 8.5.4:

Source: MPI 5.0 report: Chapter 3.7 Nonblocking Communication, 19.1.8 Additional Support for Fortran Register-Memory-Synchronization and 19.1.17 Problems with Code Movement and Register Optimization, J3/23-007r1: 8.5.4 ASYNCHRONOUS attribute

An array section such as a(2, :) is not contiguous in column-major order (see Arrays Fundamentals and Array Arguments: Assumed Shape, Size and Rank). With mpi_f08 and a library that sets MPI_SUBARRAYS_SUPPORTED, it can be passed directly, and the compiler passes a descriptor so that the library gathers the elements itself, with no copy that could be freed too early.

Source: MPI 5.0 report: Chapter 3 Point-to-Point Communication and 19.1.2 Fortran Support Through the mpi_f08 Module

With the mpi module or mpif.h a section is passed by sequence association, which makes the compiler copy it into a contiguous temporary; that is harmless for a blocking call and a bug for a non-blocking one, because the temporary may be released before the transfer completes (MPI 5.0 section 19.1.12). The traditional remedy is a derived datatype, shown next, and it remains the better choice for strided data on the performance-sensitive path.

An MPI datatype tells the library what the bytes of a buffer are. For intrinsic types the names are MPI_INTEGER, MPI_REAL, MPI_DOUBLE_PRECISION, MPI_COMPLEX, MPI_LOGICAL and MPI_CHARACTER, for the default kinds, plus optional size-specific names (MPI_INTEGER8 and MPI_REAL8, for example). The Fortran standard does not say that real64 is the same as the default double precision, so code that uses kinds from Intrinsic Modules and iso_fortran_env should ask MPI for the matching datatype with MPI_Type_create_f90_real, MPI_Type_create_f90_integer or MPI_Type_create_f90_complex, passing the result of precision/range or selected_int_kind (MPI 5.0 section 19.1.9). The heatflow kinds are used below.

Source: MPI 5.0 report: Chapter 5 Datatypes and 19.1.9 Additional Support for Fortran Numeric Intrinsic Types

The datatypes are the right size, and a reduction with them gives the expected 10.0 and 4 * 2147483647 + 6 for the 64-bit sum (MPI_INTEGER8 stands for integer(lp) here, which is correct only because lp is int64). The two F lines are the point: the datatype that MPI_Type_create_f90_real returns is a distinct handle and so compares unequal to MPI_DOUBLE_PRECISION, even though it describes the same layout. Compare sizes with MPI_Type_size, not handles.

Mistakes in the datatype are not diagnosed. This program sums integers while describing them as MPI_REAL, and the result is not an integer sum:

Real output with mpirun -n 4, which compiles without a warning:

Source: MPI 5.0 report: Chapter 5 Datatypes and Chapter 6 Collective Communication

(With small integers the same mistake can look correct, because their bit patterns are tiny subnormal reals that add like integers, which is why this example uses 10000000 plus the rank.)

For data that is not contiguous, build a derived datatype and commit it. MPI_Type_vector(count, blocklength, stride, oldtype, newtype) describes count blocks of blocklength elements, with the start of consecutive blocks stride elements apart; MPI_Type_create_subarray describes a rectangular piece of a larger array, and MPI_Type_create_struct a derived type (which should be sequence or bind(c) so its layout is defined, see Derived Types). The program below sends one grid row of a column-major array, whose elements are nx apart, and stores it into another row of the receiver’s array.

Real output on rank 1 (mpirun -n 2):

The starting address passed to MPI_Send and MPI_Recv is the first element of the row (a(2, 1), a(4, 1)), and the count is one instance of the type. A datatype must be committed with MPI_Type_commit before use and released with MPI_Type_free afterwards. A derived type with allocatable or pointer components cannot be described this way, since its components are not at fixed offsets: send the components separately.

Source: MPI 5.0 report: Chapter 5 Datatypes and 19.1.9 Additional Support for Fortran Numeric Intrinsic Types

The decomposition arithmetic (who is my neighbour?) can be left to MPI. MPI_Dims_create chooses a balanced grid of ranks for a given count, MPI_Cart_create builds a communicator with that structure (and may reorder the ranks to match the network when reorder is true), MPI_Cart_coords converts a rank to coordinates, and MPI_Cart_shift returns the source and destination of a shift along a dimension, with MPI_PROC_NULL at a non-periodic edge (MPI 5.0 section 8.5).

Real output with mpirun -n 4:

Source: MPI 5.0 report: 8.5 Topology Constructors

Coordinates and dimension numbers are counted from zero, unlike Fortran array indices, and the value of MPI_PROC_NULL (-1 here) is implementation defined, so always compare with the constant. The .true. argument of MPI_Cart_create is reorder.

By default an MPI error is fatal: the default handler of a communicator is MPI_ERRORS_ARE_FATAL and aborts the job (MPI_ERRORS_RETURN is the other predefined one, and MPI_Comm_create_errhandler builds your own). Programs that want to react set MPI_ERRORS_RETURN and pass the optional ierror argument, then compare with MPI_SUCCESS. The error code is implementation specific; MPI_Error_class maps it to a portable class such as MPI_ERR_RANK, MPI_ERR_TRUNCATE, MPI_ERR_COUNT or MPI_ERR_TYPE, and MPI_Error_string produces text.

Real output with mpirun -n 2:

Printing the string from MPI_Error_string gave a multi-line stack with addresses on this MPICH, which is why the example prints the class. After an error with MPI_ERRORS_RETURN the standard does not promise that the communicator is still usable; for most programs the right reaction is to report and MPI_Abort. Errors in Fortran itself (an out-of-bounds access, an unallocated array) are not MPI errors and are handled as in Error Handling and Program Termination and Testing, Debugging and Profiling.

Source: MPI 5.0 report: 9.3 Error Handling and 9.4 Error Codes and Classes

The heatflow grid has nx points in x and ny rows in y. A one-dimensional decomposition gives each rank ny / nprocs rows stored as u(nx, 0:nloc + 1): rows 1:nloc are owned, and rows 0 and nloc + 1 are the halo that mirrors the neighbours' edge rows. Because the second index is the slow one in column-major order, a whole row u(:, j) is contiguous, so no derived datatype is needed. Every step the halo is refreshed, the interior is updated and the boundary rows of the first and last rank are held at zero.

The version in the parallel/ directory of the heatflow project posts one MPI_Sendrecv per direction and is the simplest correct one; the listing below is the non-blocking variant that overlaps communication with computation. Both give the same result.

Built with mpifort heatflow_kinds.f90 09_halo_nonblocking.f90 and run with mpirun -n 4:

Source: MPI 5.0 report: Chapter 3 Point-to-Point Communication and Chapter 6 Collective Communication

This is the same value the MPI_Sendrecv version in the project, the coarray version and the OpenMP version (Coarrays, Images and Synchronization and OpenMP and OpenACC) report for four ranks or images. The hot spot is placed at row nloc / 2 of rank 0, so the number depends on how many rows a rank owns: the same program printed 0.512 with one rank and 0.509 with two when it was tried, which is a property of the test setup and not an error. The important details are:

An MPI installation provides compiler wrappers that add the include and library paths and the MPI library to a normal Fortran compiler: mpifort (also installed as mpif90) for Fortran, plus mpicc and mpicxx. mpifort -show prints the real command line. On the MPICH used here it was:

The wrapper is a distribution’s choice (this one adds -O2, which is why the optimization example above misbehaves without asynchronous, and -fallow-argument-mismatch, which hides errors for the mpi module and mpif.h, though not for mpi_f08). Programs are started with mpirun -n 4 ./a.out or mpiexec -n 4 ./a.out; both are the launcher of the installed implementation, and the options beyond -n differ between MPICH and Open MPI. The same -std=f2023 -Wall -Wextra -pedantic -fcheck=all options used in the rest of the section can be passed to mpifort, and every example compiled with them (the deliberately wrong ones excepted). A run with 16 ranks failed in this container during MPI_Init (a UCX out-of-memory error), so the examples use at most four; use no more ranks than the machine has cores unless the launcher documentation says oversubscription is supported.

Two build-system integrations were tried and worked:

CMake has a FindMPI module whose Fortran component reports whether mpi_f08 is usable:

Source: CMake: FindMPI

and mpirun -n 2 build/hello printed the two ranks. The version reported is that of the MPI standard the library implements, 4.1.

fpm has an mpi metapackage (0.13.0 was used):

Source: Fortran Package Manager: metapackages

The two [fortran] entries are not optional: without them fpm refused to build with "metapackage does require implicit-typing" (and the same for implicit-external), because the metapackage also has to support mpif.h. They loosen the checks for the whole package, so a package that wants strict implicit none (type, external) enforcement in its own sources must enforce it with its own compiler flags. fpm run --runner "mpirun -n 2" started two ranks, and plain fpm run ran one. See Building with Make, CMake and Meson and Fortran Package Manager (fpm) for the surrounding build structure, and Compilers and Compiler Options for the options mentioned above.

A last practical constraint is that the mpi_f08 module files are produced by one specific Fortran compiler. Open MPI and MPICH each ship modules built for the compiler they were configured with, and using a different compiler (or even a different major version of the same one) typically gives an error when the module is read. For that reason, and for the ABI, build the whole program, MPI library and mpifort included, with a matching toolchain, and treat the MPI 5.0 ABI (chapter 20) as a possible future way around it.

Both are SPMD models in which every process runs the same program, but the abstraction differs:

The same strip decomposition is a single statement per halo with coarrays and about a dozen with MPI, and the result of the heatflow test was identical (max = 0.380 with four ranks or images). MPI remains the choice when the program must run on a cluster, interoperate with C or Python MPI code, or use the MPI parallel I/O and one-sided interfaces. Coarray runtimes are often built on MPI, and OpenCoarrays used the MPICH installed in this container. Combine MPI with OpenMP for the shared-memory level inside a node (the thread-support level is requested with MPI_Init_thread).

See also Performance and Memory Layout for why a contiguous halo row matters, and Common Gotchas and Pitfalls.