exp_repro: Reproducible implementation of exp()

.testing/Makefile

@@ -113,6 +113,7 @@ FCFLAGS_FMS ?= $(FCFLAGS)
 # - These flags can be configured outside of the Makefile, either with
 #   config.mk or as environment variables.
 
+LDFLAGS_OPT ?=
 LDFLAGS_COVERAGE ?= --coverage
 LDFLAGS_USER ?=
 
@@ -248,6 +249,7 @@ COV_FCFLAGS := FCFLAGS="$(FCFLAGS_COVERAGE) $(FCFLAGS_DEPS)"
 
 # Linker flags
 MOM_LDFLAGS := LDFLAGS="$(LDFLAGS_DEPS) $(LDFLAGS_USER)"
+OPT_LDFLAGS := LDFLAGS="$(LDFLAGS_OPT) $(LDFLAGS_DEPS) $(LDFLAGS_USER)"
 COV_LDFLAGS := LDFLAGS="$(LDFLAGS_COVERAGE) $(LDFLAGS_DEPS) $(LDFLAGS_USER)"
 
 # Environment variable configuration
@@ -258,13 +260,13 @@ $(BUILD)/asymmetric/Makefile: MOM_ENV += $(ASYMMETRIC_FCFLAGS) $(MOM_LDFLAGS) \
 $(BUILD)/repro/Makefile: MOM_ENV += $(REPRO_FCFLAGS) $(MOM_LDFLAGS)
 $(BUILD)/openmp/Makefile: MOM_ENV += $(OPENMP_FCFLAGS) $(MOM_LDFLAGS)
 $(BUILD)/target/Makefile: MOM_ENV += $(TARGET_FCFLAGS) $(MOM_LDFLAGS)
-$(BUILD)/opt/Makefile: MOM_ENV += $(OPT_FCFLAGS) $(MOM_LDFLAGS)
+$(BUILD)/opt/Makefile: MOM_ENV += $(OPT_FCFLAGS) $(OPT_LDFLAGS)
 $(BUILD)/opt_target/Makefile: MOM_ENV += $(OPT_FCFLAGS) $(MOM_LDFLAGS)
 $(BUILD)/coupled/Makefile: MOM_ENV += $(SYMMETRIC_FCFLAGS) $(MOM_LDFLAGS)
 $(BUILD)/nuopc/Makefile: MOM_ENV += $(SYMMETRIC_FCFLAGS) $(MOM_LDFLAGS)
 $(BUILD)/cov/Makefile: MOM_ENV += $(COV_FCFLAGS) $(COV_LDFLAGS)
 $(BUILD)/unit/Makefile: MOM_ENV += $(COV_FCFLAGS) $(COV_LDFLAGS)
-$(BUILD)/timing/Makefile: MOM_ENV += $(OPT_FCFLAGS) $(MOM_LDFLAGS)
+$(BUILD)/timing/Makefile: MOM_ENV += $(OPT_FCFLAGS) $(OPT_LDFLAGS)
 
 # Configure script flags
 MOM_ACFLAGS := --with-gsw --with-cvmix
@@ -326,7 +328,7 @@ $(BUILD)/%/config.status: \
 	cd $(@D) && $(MOM_ENV) ../../../configure -n $(MOM_ACFLAGS) \
 	 || (cat config.log && false)
 
-../configure:
+../configure: FORCE
 	$(MAKE) -C ../ configure
 
 ALL_EXECS = symmetric asymmetric repro openmp opt opt_target coupled nuopc \

ac/configure.ac

@@ -183,6 +183,13 @@ AS_IF(
 )
 
 
+# Detect if parentheses are honored in floating point arithmetic
+MOM6_FC_FAST_RINT
+AS_IF([test -n "$FAST_RINT_FCFLAGS"],
+  [FCFLAGS="$FCFLAGS $FAST_RINT_FCFLAGS"]
+)
+
+
 # OpenMP configuration
 
 # NOTE: AC_OPENMP fails on `Fortran` for Autoconf <2.69 due to a m4 bug.

ac/m4/mom6_fc_fast_rint.m4

@@ -0,0 +1,75 @@
+dnl Determine the flag required to honor parentheses in floating-point
+dnl expressions (i.e., prevent reassociation of (a + b) - b into a).
+dnl
+dnl This is required for the fast round-to-nearest-integer trick:
+dnl   K = (x + round_bias) - round_bias
+dnl where round_bias = 1.5 * 2^52.
+dnl
+dnl Compiler flags that enable this:
+dnl   GCC gfortran:    -fprotect-parens (default with -std=f2008+)
+dnl   Intel ifx/ifort: -assume protect_parens
+dnl   NVHPC nvfortran: (appears to be default)
+dnl
+AC_DEFUN([MOM6_FC_FAST_RINT], [
+  AC_ARG_ENABLE([fast-rint],
+    [AS_HELP_STRING([--disable-fast-rint],
+      [do not use fast round-to-nearest-integer])])
+  AS_IF([test "x$enable_fast_rint" != xno], [
+    AC_CACHE_CHECK([for $FC option to honor parentheses],
+        [mom6_cv_fc_fast_rint], [
+          mom6_cv_fc_fast_rint="unsupported"
+          ac_fc_fr_FCFLAGS_save=${FCFLAGS}
+          AC_LANG_PUSH([Fortran])
+          for ac_flag in none \
+            -fprotect-parens \
+            "-assume protect_parens"
+          do
+            test "$ac_flag" != none \
+              && FCFLAGS="$ac_fc_fr_FCFLAGS_save $ac_flag"
+            AC_RUN_IFELSE([
+              AC_LANG_PROGRAM([], [
+dnl ---
+      real, parameter :: bias = 1.5 * 2.**(digits(1.) - 1)
+      if (fast_rint(1.23) == 1.) stop 0
+      stop 1
+      contains
+      function fast_rint(x) result(y)
+        real, intent(in) :: x
+        real :: y
+        y = (x + bias) - bias
+      end function
+dnl ---
+              ])
+            ], [
+              mom6_cv_fc_fast_rint="$ac_flag"
+              break
+            ], [], [
+              mom6_cv_fc_fast_rint="unsupported"
+              break
+            ])
+          done
+          AC_LANG_POP([Fortran])
+          FCFLAGS=$ac_fc_fr_FCFLAGS_save
+        ]
+      )
+      AS_CASE([$mom6_cv_fc_fast_rint],
+        [none], [
+          mom6_cv_fc_fast_rint="none needed"
+          FAST_RINT_FCFLAGS=""
+          AC_DEFINE([ENABLE_FAST_RINT], [1],
+            [Define to 1 if parentheses are protected in FP expressions])
+        ],
+        [unsupported], [
+          AC_MSG_WARN(
+            [No known flag found to protect parentheses; using ieee_rint fallback])
+          FAST_RINT_FCFLAGS=""
+        ],
+        [
+          FAST_RINT_FCFLAGS=$mom6_cv_fc_fast_rint
+          AC_DEFINE([ENABLE_FAST_RINT], [1],
+            [Define to 1 if parentheses are protected in FP expressions])
+        ]
+      )
+  ])
+  AC_SUBST([FAST_RINT_FCFLAGS])
+])

ac/makedep

@@ -27,13 +27,19 @@ re_cpp_endif = re.compile(r"^ *# *endif")
 re_cpp_include = re.compile(r"^ *# *include *[<\"']([a-zA-Z_0-9\.]+)[>\"']")
 re_f90_include = re.compile(r"^ *include +[\"']([a-zA-Z_0-9\.]+)[\"']")
 re_program = re.compile(r"^ *program +([a-z_0-9]+)", re.IGNORECASE)
-re_end = re.compile(r"^ *end *(module|procedure) ", re.IGNORECASE)
-# NOTE: This excludes comments and tokens with substrings containing `function`
-# or `subroutine`, but will fail if the keywords appear in other contexts.
+
+re_end = re.compile(r"^ *end *(module|submodule|program) ", re.IGNORECASE)
+# NOTE: `end procedure` is intentionally excluded from re_end.  it should
+# not reset external namespace since we remain inside the containing
+# submodule or module.
+
 re_procedure = re.compile(
     r"^[^!]*(?<![a-z_])(function|subroutine)(?![a-z_])",
     re.IGNORECASE
 )
+# NOTE: re_procedure excludes comments and tokens with substrings
+# containing `function` or `subroutine`, but will fail if the keywords
+# appear in other contexts.
 
 
 # Preprocessor expression tokenization
@@ -366,16 +372,11 @@ def link_obj(obj, o2uses, mod2o, all_modules, parent2subobjs, o2mods):
                 if m in parent2subobjs:
                     for subobj in parent2subobjs[m]:
                         recur(subobj, depth=depth+1)
-        else:
-            return
-        uses = [m for m in o2uses[obj] if m in all_modules]
-        if len(uses) > 0:
-            ouses = [mod2o[m] for m in uses]
+            # Follow the modules used by this object (including submodules)
+            uses = [m for m in o2uses.get(obj, []) if m in all_modules]
             for m in uses:
                 o = mod2o[m]
                 recur(o, depth=depth+1)
-            return
-        return
     olst = []
     recur(obj)
     return sorted(set(olst))

config_src/drivers/timing_tests/time_MOM_intrinsic_functions.F90

@@ -0,0 +1,172 @@
+! This file is part of MOM6, the Modular Ocean Model version 6.
+! See the LICENSE file for licensing information.
+! SPDX-License-Identifier: Apache-2.0
+
+!> Timing tests for MOM_intrinsic_functions (exp_repro)
+program time_MOM_intrinsic_functions
+
+use, intrinsic :: iso_fortran_env, only : int64, real64
+use MOM_intrinsic_functions, only : exp_repro
+
+implicit none
+
+integer, parameter :: npts = 100000
+  !< Number of test points
+integer, parameter :: niter = 200
+  !< Number of timing iterations
+real, parameter :: xmin = -10.
+  !< Minimum x value [nondim]
+real, parameter :: xmax = 10.
+  !< Maximum x value [nondim]
+
+real, allocatable :: x(:), val(:)
+integer(kind=int64) :: count_rate, c1, c2
+real(kind=real64) :: clock_rate, time_intrinsic, time_repro
+real(kind=real64) :: time_scalar_baseline, time_scalar_intrinsic, time_scalar_repro
+real :: I_npts
+real, volatile :: scalar_x, scalar_val
+integer :: i, j
+
+call system_clock(count_rate=count_rate)
+clock_rate = real(count_rate, real64)
+
+allocate(x(npts), val(npts))
+I_npts = 1. / (npts - 1)
+do i = 1, npts
+  x(i) = xmin + (i - 1) * ((xmax - xmin) * I_npts)
+enddo
+
+print '("=== MOM_intrinsic_functions timing ===")'
+print '("npts = ", i0, ", niter = ", i0)', npts, niter
+print '("x range: [", f6.1, ", ", f6.1, "]")', xmin, xmax
+print *
+
+! Warm-up and time intrinsic exp()
+do j = 1, 3
+  val = exp(x)
+  if (val(1) < 0.) val(1) = 0.
+enddo
+call system_clock(count=c1)
+do j = 1, niter
+  val = exp(x)
+  if (val(1) < 0.) val(1) = 0.
+enddo
+call system_clock(count=c2)
+time_intrinsic = real(c2 - c1, real64) / clock_rate / niter / npts * 1e9
+
+print '("exp() time/elem:", t30, f8.2, " ns")', time_intrinsic
+print '("  sum (to prevent elision):", t30, ES12.5)', sum(val)
+
+! Warm-up and time exp_repro()
+do j = 1, 3
+  val = exp_repro(x)
+  if (val(1) < 0.) val(1) = 0.
+enddo
+call system_clock(count=c1)
+do j = 1, niter
+  val = exp_repro(x)
+  if (val(1) < 0.) val(1) = 0.
+enddo
+call system_clock(count=c2)
+time_repro = real(c2 - c1, real64) / clock_rate / niter / npts * 1e9
+
+print '("exp_repro() time/elem:", t30, f8.2, " ns")', time_repro
+print '("  sum (to prevent elision):", t30, ES12.5)', sum(val)
+
+print *
+print '("slowdown factor:", t30, f8.2, "x")', time_repro / time_intrinsic
+
+print *
+print '("=== scalar loop-carried timing ===")'
+
+! Warm-up and time a scalar dependency chain without exp().
+! This gives an approximate loop overhead for interpreting the scalar timings.
+do j = 1, 3
+  scalar_x = 0.125
+  do i = 1, npts
+    scalar_val = 1. + scalar_x
+    scalar_x = scalar_x + (1.e-7 * (scalar_val - 1.))
+    if (scalar_x > 0.5) scalar_x = scalar_x - 1.
+  enddo
+enddo
+
+call system_clock(count=c1)
+do j = 1, niter
+  scalar_x = 0.125
+  do i = 1, npts
+    scalar_val = 1. + scalar_x
+    scalar_x = scalar_x + (1.e-7 * (scalar_val - 1.))
+    if (scalar_x > 0.5) scalar_x = scalar_x - 1.
+  enddo
+enddo
+call system_clock(count=c2)
+time_scalar_baseline = real(c2 - c1, real64) / clock_rate / niter / npts * 1e9
+
+print '("baseline scalar time/call:", t30, f8.2, " ns")', time_scalar_baseline
+print '("  final scalar value:", t30, ES12.5)', scalar_val
+
+! Warm-up and time intrinsic exp() in a scalar loop with a carried dependency.
+! This prevents the loop from being converted into independent vector lanes.
+do j = 1, 3
+  scalar_x = 0.125
+  do i = 1, npts
+    scalar_val = exp(scalar_x)
+    scalar_x = scalar_x + (1.e-7 * (scalar_val - 1.))
+    if (scalar_x > 0.5) scalar_x = scalar_x - 1.
+  enddo
+enddo
+
+call system_clock(count=c1)
+do j = 1, niter
+  scalar_x = 0.125
+  do i = 1, npts
+    scalar_val = exp(scalar_x)
+    scalar_x = scalar_x + (1.e-7 * (scalar_val - 1.))
+    if (scalar_x > 0.5) scalar_x = scalar_x - 1.
+  enddo
+enddo
+call system_clock(count=c2)
+time_scalar_intrinsic = real(c2 - c1, real64) / clock_rate / niter / npts * 1e9
+
+print '("exp() scalar time/call:", t30, f8.2, " ns")', time_scalar_intrinsic
+print '("  final scalar value:", t30, ES12.5)', scalar_val
+
+! Warm-up and time exp_repro() in the same scalar loop pattern.
+do j = 1, 3
+  scalar_x = 0.125
+  do i = 1, npts
+    scalar_val = exp_repro(scalar_x)
+    scalar_x = scalar_x + (1.e-7 * (scalar_val - 1.))
+    if (scalar_x > 0.5) scalar_x = scalar_x - 1.
+  enddo
+enddo
+
+call system_clock(count=c1)
+do j = 1, niter
+  scalar_x = 0.125
+  do i = 1, npts
+    scalar_val = exp_repro(scalar_x)
+    scalar_x = scalar_x + (1.e-7 * (scalar_val - 1.))
+    if (scalar_x > 0.5) scalar_x = scalar_x - 1.
+  enddo
+enddo
+call system_clock(count=c2)
+time_scalar_repro = real(c2 - c1, real64) / clock_rate / niter / npts * 1e9
+
+print '("exp_repro() scalar time/call:", t30, f8.2, " ns")', time_scalar_repro
+print '("  final scalar value:", t30, ES12.5)', scalar_val
+
+print *
+print '("scalar slowdown factor:", t30, f8.2, "x")', &
+    time_scalar_repro / time_scalar_intrinsic
+print '("exp() minus baseline:", t30, f8.2, " ns")', &
+    time_scalar_intrinsic - time_scalar_baseline
+print '("exp_repro() minus baseline:", t30, f8.2, " ns")', &
+    time_scalar_repro - time_scalar_baseline
+print '("adjusted slowdown factor:", t30, f8.2, "x")', &
+    (time_scalar_repro - time_scalar_baseline) &
+      / (time_scalar_intrinsic - time_scalar_baseline)
+
+deallocate(x, val)
+
+end program time_MOM_intrinsic_functions

config_src/drivers/unit_tests/test_MOM_intrinsic_functions.F90

@@ -0,0 +1,18 @@
+! This file is part of MOM6, the Modular Ocean Model version 6.
+! See the LICENSE file for licensing information.
+! SPDX-License-Identifier: Apache-2.0
+
+!> Driver for MOM_intrinsic_functions unit tests
+program test_MOM_intrinsic_functions
+
+use MOM_error_handler, only : set_skip_mpi
+use MOM_intrinsic_functions, only : intrinsic_functions_unit_tests
+use MOM_intrinsic_functions_tests, only : run_intrinsic_functions_tests
+
+call set_skip_mpi(.true.)
+
+if (intrinsic_functions_unit_tests(.true.)) stop 1
+
+call run_intrinsic_functions_tests
+
+end program test_MOM_intrinsic_functions

src/framework/MOM_error_handler.F90

@@ -27,7 +27,7 @@ module MOM_error_handler
 public :: is_root_pe, stdlog, stdout
 !> Integer parameters encoding the severity of an error message
 public :: NOTE, WARNING, FATAL
-public :: disable_fatal_errors, enable_fatal_errors, set_skip_mpi
+public :: disable_fatal_errors, enable_fatal_errors, set_skip_mpi, query_skip_mpi
 
 integer :: verbosity = 6
 !< Verbosity level:
@@ -147,6 +147,14 @@ subroutine set_skip_mpi(skip)
 
 end subroutine set_skip_mpi
 
+!> Query whether MPI-dependent behaviors should be skipped
+function query_skip_mpi() result(skip)
+  logical :: skip !< True if MPI should be skipped
+
+  skip = skip_mpi_dep
+
+end function query_skip_mpi
+
 !> This provides a convenient interface for writing an error message
 !! with run-time filter based on a verbosity and the severity of the error.
 subroutine MOM_error(level, message, all_print)