3. OpenMP and MPI
FoBiS knows each compiler's switch for OpenMP and its MPI wrapper: a mode says openmp = True or mpi = True, and the right flags or command follow.
OpenMP
A new program, sweep, computes the errors of both rules for 12 grids. The cases are independent, so OpenMP threads compute them in parallel:
!< The errors of both rules for 2, 4, ... 4096 intervals: the cases computed in parallel by OpenMP threads.
program sweep
!< The errors of both rules for 2, 4, ... 4096 intervals: the cases computed in parallel by OpenMP threads.
use quad_rules, only : R8P, trapezoid, simpson
use integrands, only : sine
!$ use omp_lib, only : omp_get_max_threads
implicit none
integer, parameter :: cases = 12 !< Number of cases.
real(R8P), parameter :: pi = acos(-1._R8P) !< Pi.
real(R8P) :: error(2, cases) !< Errors, trapezoid and Simpson, of each case.
integer :: k !< Counter.
!$omp parallel do
do k = 1, cases
error(1, k) = abs(trapezoid(sine, 0._R8P, pi, 2 ** k) - 2._R8P)
error(2, k) = abs(simpson(sine, 0._R8P, pi, 2 ** k) - 2._R8P)
enddo
!$ print '(A, I0)', 'threads: ', omp_get_max_threads()
print '(A6, 2A18)', 'n', 'trapezoid error', 'simpson error'
do k = 1, cases
print '(I6, 2ES18.4)', 2 ** k, error(:, k)
enddo
endprogram sweepThe rules are pure: they can be called from parallel threads safely (and a parallel construct cannot appear inside a pure procedure, so the parallelism belongs to the caller). A mode for it:
[openmp]
help = The sweep, multi-threaded with OpenMP
template = template-gnu
build_profile = release
cflags = -c
openmp = True
build_dir = ./build/openmp/
target = src/sweep.F90$ fobis build --mode openmp
Builder options
Directories
Root of source files directory: "src .fobis_deps/penf"
Building directory: "build/openmp"
Compiled-objects .o directory: "build/openmp/obj"
Compiled-objects .mod directory: "build/openmp/mod"
Compiler options
Vendor: "gnu"
Compiler command: "gfortran"
Module directory switch: "-J"
Compiling flags: "-O3 -funroll-loops -c -fopenmp"
Linking flags: "-fopenmp"
Preprocessing flags: ""
Coverage: False
Profile: False
Archiver: ar -rcs
Ranlib: ranlib
Preprocessor used: None
Preprocessor args:
Preprocessor output directory: None
Preprocessor extensions processed: []
Building src/sweep.F90
Compiling src/quad_rules.F90 serially
Compiling src/integrands.F90 serially
Compiling src/sweep.F90 serially
Linking build/openmp/sweep
Target src/sweep.F90 has been successfully built-fopenmp was added to both the compiling and the linking flags.
$ OMP_NUM_THREADS=4 ./build/openmp/sweep
threads: 4
n trapezoid error simpson error
2 4.2920E-01 9.4395E-02
4 1.0388E-01 4.5598E-03
8 2.5768E-02 2.6917E-04
16 6.4297E-03 1.6591E-05
32 1.6066E-03 1.0334E-06
64 4.0161E-04 6.4530E-08
128 1.0040E-04 4.0323E-09
256 2.5100E-05 2.5200E-10
512 6.2749E-06 1.5752E-11
1024 1.5687E-06 9.8277E-13
2048 3.9218E-07 5.7732E-14
4096 9.8046E-08 3.9968E-15MPI
Another program splits [0, pi] among the MPI processes; each integrates its piece, and mpi_reduce adds them up:
!< Simpson's rule over [0, pi] split among the MPI processes: each integrates its piece, rank 0 adds them up.
program integrate_mpi
!< Simpson's rule over [0, pi] split among the MPI processes: each integrates its piece, rank 0 adds them up.
use mpi
use quad_rules, only : R8P, simpson
use integrands, only : sine
implicit none
real(R8P), parameter :: pi = acos(-1._R8P) !< Pi.
integer, parameter :: n = 64 !< Intervals of each piece.
real(R8P) :: a, b !< Bounds of this process's piece.
real(R8P) :: piece !< Integral over this process's piece.
real(R8P) :: total !< Sum of the pieces.
integer :: rank, size, ierr !< MPI rank, number of processes, error code.
call mpi_init(ierr)
call mpi_comm_rank(mpi_comm_world, rank, ierr)
call mpi_comm_size(mpi_comm_world, size, ierr)
a = pi * rank / size
b = pi * (rank + 1) / size
piece = simpson(sine, a, b, n)
call mpi_reduce(piece, total, 1, mpi_double_precision, mpi_sum, 0, mpi_comm_world, ierr)
if (rank == 0) print '(I0, A, I0, A, ES12.4)', size, ' processes x ', n, ' intervals, error ', abs(total - 2._R8P)
call mpi_finalize(ierr)
endprogram integrate_mpi[mpi]
help = Simpson's rule split among MPI processes
template = template-gnu
build_profile = release
cflags = -c
mpi = True
build_dir = ./build/mpi/
target = src/integrate_mpi.F90With mpi = True FoBiS compiles with the MPI wrapper of the compiler (mpif90 for gfortran):
$ fobis build --mode mpi --quiet
Building src/integrate_mpi.F90
Compiling src/quad_rules.F90 serially
Compiling src/integrands.F90 serially
Compiling src/integrate_mpi.F90 serially
Linking build/mpi/integrate_mpi
Target src/integrate_mpi.F90 has been successfully built$ mpirun -np 2 ./build/mpi/integrate_mpi && mpirun -np 4 ./build/mpi/integrate_mpi
2 processes x 64 intervals, error 4.0323E-09
4 processes x 64 intervals, error 2.5200E-10Two processes with 64 intervals each make 128 intervals over the whole range: the error is the one sweep printed for 128 intervals.
fobis run starts the mode's program through a launcher: {ranks} is replaced by --ranks. A mode can keep its launcher in the fobos (launcher = mpirun -np {ranks}), and fobis test takes the same options, so MPI tests run with one command:
$ fobis run --mode mpi --no-build --launcher "mpirun -np {ranks}" --ranks 4
4 processes x 64 intervals, error 2.5200E-10A timer in C
sweep could time each grid with clock_gettime, which Fortran reaches only through C. A C file among the sources is enough:
/* src/clock.c */
#define _POSIX_C_SOURCE 199309L
#include <time.h>
double wall_seconds(void) {
struct timespec now;
clock_gettime(CLOCK_MONOTONIC, &now);
return (double)now.tv_sec + 1.0e-9 * (double)now.tv_nsec;
}with an interface function wall_seconds() bind(c, name="wall_seconds") returning real(c_double) in a module. FoBiS compiles clock.c with the family's C compiler (gcc here, nvc for nvfortran, icx for ifx) and links it into every program; c_cflags holds its flags, since the Fortran cflags are not C's. Cookbook: mixed-language projects.
MPI outputs
The MPI outputs were recorded with Open MPI 5.0 and gfortran 15 (the documentation's CI has no MPI); the rest of the page, with gfortran 13.
What you learned
openmp = True, mpi = True. For libraries found through pkg-config or a wrapper (HDF5, NetCDF, FFTW, MPI too) see external libraries.