Floating-point exceptions

components/omega/OmegaBuild.cmake

@@ -28,6 +28,7 @@ macro(common)
   option(OMEGA_LOG_FLUSH "Turn on unbuffered logging (default OFF)." OFF)
   option(OMEGA_TEST_CDASH "Turn on CDash support (default ON)." ON)
   option(OMEGA_EXTERNAL_PROF "Integration of Omega timers with external profiling tools (default OFF)." OFF)
+  option(OMEGA_FPEXCEPTS_IN_TESTS "Enable floating-point exceptions in Omega unit tests (default OFF)." OFF)
 
   if("${OMEGA_BUILD_TYPE}" STREQUAL "Debug" OR "${OMEGA_BUILD_TYPE}" STREQUAL "DEBUG")
     set(OMEGA_DEBUG ON)

components/omega/doc/devGuide/FloatExcept.md

@@ -0,0 +1,46 @@
+(omega-dev-float-except)=
+
+## Floating-point Exceptions
+
+Omega provides functions that enable or disable selected
+floating-point exceptions. Floating-point exception handling depends
+on the target architecture, and is not always supported.
+
+To enable floating-point exceptions `Exceptions` call
+```c++
+FloatExceptStatus Status = enableFloatExceptions(Exceptions);
+```
+The returned `FloatExceptStatus` struct contains two fields.
+`Status.Err` holds an error code indicating success or failure.
+On success, `Status.OldExceptions` contains exceptions that were in effect
+before the call to the enable function.
+`Exceptions` is an integer mask, which can be composed of exception
+types defined in `<cfenv>`. They are implementation defined, but
+typically include
+- `FE_DIVBYZERO`
+- `FE_INEXACT`
+- `FE_INVALID`
+- `FE_OVERFLOW`
+- `FE_UNDERFLOW`
+
+For example, to enable divide-by-zero and underflow exceptions do
+```c++
+FloatExceptStatus Status = enableFloatExceptions(FE_DIVBYZERO | FE_UNDERFLOW);
+```
+
+You can call `enableFloatExceptions` without an argument, which
+enables
+- `FE_DIVBYZERO`
+- `FE_INVALID`
+- `FE_OVERFLOW`
+
+To disable floating-point exceptions `Exceptions` call
+```c++
+FloatExceptStatus Status = disableFloatExceptions(Exceptions);
+```
+Calling this function without an argument disables all exceptions.
+
+Omega also provides a helper function `enableFloatExceptionsInTests(Exceptions)`.
+This function is used in unit tests to conditionally call
+`enableFloatExceptions(Exceptions)` based on the compiler and specified CMake options.
+See [User's Guide](#omega-user-float-except) for details.

components/omega/doc/index.md

@@ -58,6 +58,7 @@ userGuide/VertAdv
 userGuide/Forcing
 userGuide/SfcCoupling
 userGuide/Analysis
+userGuide/FloatExcept
 ```
 
 ```{toctree}
@@ -108,6 +109,7 @@ devGuide/VertAdv
 devGuide/Forcing
 devGuide/SfcCoupling
 devGuide/Analysis
+devGuide/FloatExcept
 ```
 
 ```{toctree}

components/omega/doc/userGuide/FloatExcept.md

@@ -0,0 +1,20 @@
+(omega-user-float-except)=
+
+## Floating-point exceptions
+
+Omega has support for enabling floating-point exceptions. In
+certain conditions floating-point exceptions may be enabled by
+default in Omega unit tests. Currently, this happens only when
+the GNU compiler is used. A CMake option `OMEGA_FPEXCEPTS_IN_TESTS`
+can be set to enable floating-point exceptions in unit tests with
+any compiler. This enables the following exceptions
+- divide by zero
+- invalid operation
+- overflow
+
+Note that enabling floating-point exceptions in other compilers
+may require adjustment to compiler flags. For example, `clang`
+has the command-line option  `-ffp-exception-behavior` which
+defaults to `ignore`. Enabling floating-point exceptions without
+changing this flag to either `maytrap` or `strict` can cause
+false positives.

components/omega/src/CMakeLists.txt

@@ -35,6 +35,25 @@ if(OMEGA_SINGLE_PRECISION)
   )
 endif()
 
+include(CheckCXXSymbolExists)
+# check if cfenv provides feenableexcept which is a gnu extension for enabling floating-point exceptions
+check_cxx_symbol_exists(feenableexcept "cfenv" OMEGA_HAVE_FEENABLEEXCEPT)
+if (OMEGA_HAVE_FEENABLEEXCEPT)
+  target_compile_definitions(
+      OmegaLibFlags
+      INTERFACE
+      OMEGA_HAVE_FEENABLEEXCEPT=1
+  )
+endif()
+
+if(OMEGA_FPEXCEPTS_IN_TESTS)
+  target_compile_definitions(
+      OmegaLibFlags
+      INTERFACE
+      OMEGA_FPEXCEPTS_IN_TESTS=1
+  )
+endif()
+
 
 target_link_options(
     OmegaLibFlags

components/omega/src/base/IO.cpp

@@ -255,6 +255,10 @@ void openFileRead(
     const std::string &Filename, // [in] name (incl path) of file to open
     FileFmt InFormat             // [in] (optional) file format
 ) {
+   // Workaround for HDF5 v1.14.3 raising floating-point exceptions
+   // Disable fpes in this function and re-enable before exiting
+   // Remove this workaround after HDF5 is upgraded on every machine
+   FloatExceptStatus Status = disableFloatExceptions();
 
    int PIOErr = 0;        // internal SCORPIO/PIO return call
    int Format = InFormat; // coerce to integer for PIO calls
@@ -264,6 +268,10 @@ void openFileRead(
    if (PIOErr != PIO_NOERR)
       ABORT_ERROR("IO::openFile: PIO error opening file {} for read", Filename);
 
+   if (Status.Err.isSuccess()) {
+      enableFloatExceptions(Status.OldExceptions);
+   }
+
    return;
 
 } // End openFileRead

components/omega/src/base/IO.h

@@ -40,6 +40,7 @@
 
 #include "DataTypes.h"
 #include "Error.h"
+#include "FloatExcept.h"
 #include "mpi.h"
 #include "pio.h"
 

components/omega/src/infra/FloatExcept.cpp

@@ -0,0 +1,207 @@
+//===-- FloatExcept.cpp - Omega floating-point exception handling implementation
+//-*- C++ -*-===//
+//
+/// \file
+/// \brief Utilities for handling floating-point exceptions
+///
+/// This file contains the implementation of the floating-point exception
+/// handling functions for Omega.
+//
+//===---------------------------------------------------------------------------------------===//
+
+#include "FloatExcept.h"
+#include "OmegaKokkos.h"
+#include <Kokkos_Core.hpp>
+#include <cfenv>
+
+namespace OMEGA {
+
+// Floating-point exception handling using the feenableexcept GNU extension
+#ifdef OMEGA_HAVE_FEENABLEEXCEPT
+static FloatExceptStatus enableFloatExceptionsGnuExtension(int Exceptions) {
+   FloatExceptStatus Status;
+
+   // Clear previous exceptions
+   int Clear = std::feclearexcept(FE_ALL_EXCEPT);
+   if (Clear != 0) {
+      Status.Err +=
+          Error(ErrorCode::Fail, "Failed to clear previous exceptions");
+      return Status;
+   }
+
+   int EnableVal = feenableexcept(Exceptions);
+   if (EnableVal == -1) {
+      Status.Err +=
+          Error(ErrorCode::Fail, "Failed to enable floating-point exceptions");
+   } else {
+      Status.OldExceptions = EnableVal;
+   }
+   return Status;
+}
+
+static FloatExceptStatus disableFloatExceptionsGnuExtension(int Exceptions) {
+   FloatExceptStatus Status;
+
+   int DisableVal = fedisableexcept(Exceptions);
+
+   if (DisableVal == -1) {
+      Status.Err +=
+          Error(ErrorCode::Fail, "Failed to disable floating-point exceptions");
+   } else {
+      Status.OldExceptions = DisableVal;
+   }
+
+   return Status;
+}
+#endif
+
+// Floating-point exception handling for ARM64
+#ifdef __aarch64__
+
+// Exception bits in fpcr are shifted by this amount
+static constexpr int FPCRShift = 8;
+
+static FloatExceptStatus enableFloatExceptionsArm64(int Exceptions) {
+   FloatExceptStatus Status;
+
+   std::fenv_t env;
+   if (fegetenv(&env) != 0) {
+      Status.Err += Error(ErrorCode::Fail, "Failed to get floating-point env");
+      return Status;
+   }
+
+   // Get currently enabled exceptions
+   Status.OldExceptions = (env.__fpcr >> FPCRShift) & FE_ALL_EXCEPT;
+
+   // Clear all exceptions
+   env.__fpsr &= ~FE_ALL_EXCEPT;
+
+   // Enable selected exceptions
+   env.__fpcr |= (Exceptions << FPCRShift);
+
+   if (fesetenv(&env) != 0) {
+      Status.Err += Error(ErrorCode::Fail, "Failed to set floating-point env");
+   }
+
+   return Status;
+}
+
+static FloatExceptStatus disableFloatExceptionsArm64(int Exceptions) {
+   FloatExceptStatus Status;
+
+   std::fenv_t env;
+   if (fegetenv(&env) != 0) {
+      Status.Err += Error(ErrorCode::Fail, "Failed to get floating-point env");
+      return Status;
+   }
+
+   // Get currently enabled exceptions
+   Status.OldExceptions = (env.__fpcr >> FPCRShift) & FE_ALL_EXCEPT;
+
+   // Disable selected exceptions
+   env.__fpcr &= ~(Exceptions << FPCRShift);
+
+   if (fesetenv(&env) != 0) {
+      Status.Err += Error(ErrorCode::Fail, "Failed to set floating-point env");
+   }
+
+   return Status;
+}
+#endif
+
+// Enable selected floating-point exceptions
+// If feenableexcept exists then use it, else
+// dispatch based on architecture
+FloatExceptStatus enableFloatExceptions(int Exceptions) {
+   auto EnableLambda = [=]() {
+#ifdef OMEGA_HAVE_FEENABLEEXCEPT
+      return enableFloatExceptionsGnuExtension(Exceptions);
+#elif defined __aarch64__
+      return enableFloatExceptionsArm64(Exceptions);
+#else // unsupported arch
+      return FloatExceptStatus{Error(ErrorCode::Fail,
+                                     "Omega doesn't support enabling floating "
+                                     "point exceptions on this architecture")};
+#endif
+   };
+
+   // Enable exceptions on the main thread
+   FloatExceptStatus Status = EnableLambda();
+
+   // Enable exceptions on threads spawned by Kokkos
+#if defined(KOKKOS_ENABLE_OPENMP) || defined(KOKKOS_ENABLE_THREADS)
+   if (Kokkos::is_initialized()) {
+      const auto NumThreads = ExecSpace().concurrency();
+      // Not using Omege wrappers because we want to always use a static
+      // schedule
+      const auto Policy =
+          Kokkos::RangePolicy<ExecSpace, Kokkos::Schedule<Kokkos::Static>>(
+              0, NumThreads);
+      Kokkos::parallel_reduce(
+          Policy, KOKKOS_LAMBDA(int, Error &Err) { Err += EnableLambda().Err; },
+          Status.Err);
+   }
+#endif
+
+   return Status;
+}
+
+// Disable selected floating-point exceptions
+// If feenableexcept exists then use it, else
+// dispatch based on architecture
+FloatExceptStatus disableFloatExceptions(int Exceptions) {
+   auto DisableLambda = [=]() {
+#ifdef OMEGA_HAVE_FEENABLEEXCEPT
+      return disableFloatExceptionsGnuExtension(Exceptions);
+#elif defined __aarch64__
+      return disableFloatExceptionsArm64(Exceptions);
+#else // unsupported arch
+      return FloatExceptStatus{Error(ErrorCode::Fail,
+                                     "Omega doesn't support disabling floating "
+                                     "point exceptions on this architecture")};
+#endif
+   };
+
+   // Disable exceptions on the main thread
+   FloatExceptStatus Status = DisableLambda();
+
+   // Disable exceptions on threads spawned by Kokkos
+#if defined(KOKKOS_ENABLE_OPENMP) || defined(KOKKOS_ENABLE_THREADS)
+
+   if (Kokkos::is_initialized()) {
+      const auto NumThreads = ExecSpace().concurrency();
+      // Not using Omega wrappers because we want to always use a static
+      // schedule
+      const auto Policy =
+          Kokkos::RangePolicy<ExecSpace, Kokkos::Schedule<Kokkos::Static>>(
+              0, NumThreads);
+      Kokkos::parallel_reduce(
+          Policy,
+          KOKKOS_LAMBDA(int, Error &Err) { Err += DisableLambda().Err; },
+          Status.Err);
+   }
+#endif
+
+   return Status;
+}
+
+// A helper function to conditionally enable floating-point exceptions in Omega
+// unit tests Floating-point exceptions are enabled when explicitly requested
+// using the OMEGA_FPEXCEPTS_IN_TESTS CMake option or when using the GNU
+// compiler
+FloatExceptStatus enableFloatExceptionsInTests(int Exceptions) {
+   FloatExceptStatus Status;
+// Enable if explicitly requested
+#ifdef OMEGA_FPEXCEPTS_IN_TESTS
+   Status = enableFloatExceptions(Exceptions);
+#else
+   // Enable if using the GNU compiler
+#if defined(KOKKOS_COMPILER_GNU)
+   Status = enableFloatExceptions(Exceptions);
+#endif
+#endif
+   return Status;
+}
+
+} // namespace OMEGA
+//===----------------------------------------------------------------------===//