FPE proposal

applications/clawpack/shallow/2d/radialdam/radialdam.cpp

@@ -112,6 +112,8 @@ void run_program(fclaw2d_global_t* glob)
     }
 
     radialdam_link_solvers(glob);
+
+    fclaw_fpe_handling(glob);
 
     /* ---------------------------------------------------------------
        Run

applications/clawpack/shallow/2d/radialdam/radialdam_user.cpp

@@ -25,6 +25,8 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include "radialdam_user.h"
 
+#include <fclaw_fpe.h>
+
 #include <fclaw2d_clawpatch.h>
 #include <fclaw2d_clawpatch46_fort.h>
 #include <fclaw2d_clawpatch5_fort.h>
@@ -115,6 +117,7 @@ void radialdam_link_solvers(fclaw2d_global_t *glob)
     fclaw2d_vtable_t *vt = fclaw2d_vt();
 
     vt->problem_setup = &radialdam_problem_setup;  /* Version-independent */
+    vt->fpe_enable = fclaw_fpe_handling_macos;
 
     const user_options_t* user = radialdam_get_options(glob);
     if (user->claw_version == 4)

snippets/tst_fpe/Makefile

@@ -0,0 +1,8 @@
+TARGET=tst_fpe
+
+all:
+	gcc -g -c main_fpe.c 
+	gcc -o ${TARGET} main_fpe.o 
+
+clean:
+	rm -f *.o *~ ${TARGET}

snippets/tst_fpe/README.txt

@@ -0,0 +1,32 @@
+This example implements DYI exception handling.  It can be tested with
+
+    $ tst_fpe <fpe_id>
+
+where <fpe_id> is an FPE type to trap. 
+
+Valid types are indicated by integer 1-5 : 
+
+1 : INVALID, e.g. sqrt(-1), log(-1), acos(2) etc. 
+
+2 : OVERFLOW, e.g. 2^1024, exp(1000)
+
+3 : UNDERFLOW, e.g. denormalized numbers such as 2^(-1074), exp(-710)
+
+4 : DIVBYZERO, e.g. 1/0
+
+5 : Quiet NAN : This is useful for detecting use of uninitialized data
+
+6 : Signaling NAN : This NAN will get trapped. 
+
+Ideas for code taken from 
+
+    https://developer.apple.com/forums/thread/689159
+
+and
+
+    https://stackoverflow.com/questions/69059981/how-to-trap-floating-point-exceptions-on-m1-macs
+
+
+
+
+

snippets/tst_fpe/main_fpe.c

@@ -0,0 +1,133 @@
+#include <signal.h>      // for signal()                                         
+#include <fenv.h>       // for fegetenv(), fesetenv()        
+#include <math.h>        // for sqrt()                                           
+
+#include <stdio.h>
+#include <stdlib.h>
+
+/* Needed for print_hexval, below */
+#include <string.h>    /* Needed for memcpy */
+#include <inttypes.h>  /* Needed for PRIx64 */
+
+/* This doesn't seem to do anything */
+void fpe_signal_handler(int sig) {
+  printf("Floating point exception!\n");
+  return;
+}
+
+void enable_floating_point_exceptions() 
+{
+ fenv_t env;
+ fegetenv(&env);
+ env.__fpcr = env.__fpcr | __fpcr_trap_invalid | __fpcr_trap_overflow 
+     | __fpcr_trap_divbyzero | __fpcr_trap_underflow;
+ fesetenv(&env);
+ signal(SIGFPE, fpe_signal_handler);
+}
+
+
+static
+void set_qnan(double *f)
+{
+    /*
+     The quiet nan from math.h
+    */
+    *f = NAN;
+}
+
+static
+void set_snan(double* f)
+{
+    /* From :
+      "NaNs, Uninitialized Variables, and C++"
+      http://codingcastles.blogspot.fr/2008/12/nans-in-c.html
+    */
+    *((long long*)f) = 0x7ff0000000000001LL;
+}
+
+static
+void print_hexval(double x)
+{
+    /* Kind of cool - didn't know it was possible to do this */
+    uint64_t xn; 
+    memcpy(&xn, &x, sizeof x);
+    printf("x = %16f (%" PRIx64 ")\n", x, xn);
+}
+
+static
+double create_exception(int choice)
+{
+    double x,y;
+    switch(choice)
+    {
+        case 1:
+        x = -1;
+        y = sqrt(x);  // Or acos(x-1), log(x)
+        print_hexval(y);
+        printf("0: INVALID : y = sqrt(-1)\n");
+        break;
+
+        case 2:
+        x = 1000;
+        y = exp(x);
+        print_hexval(y);
+        printf("1: OVERFLOW : y = exp(1000)\n");
+        break;
+
+        case 3:
+        y = pow(2,-1074);
+        print_hexval(y);
+        printf("1: UNDERFLOW : y = 2^(-1074)\n");
+        break;
+
+        case 4:
+        x = 0;
+        y = 1/x;
+        print_hexval(y);
+        printf("2: DIVBYZERO : y = 1/0\n");
+        break;
+
+        case 5:
+        set_qnan(&x);
+        print_hexval(x);
+        y = x + 1;
+        printf("3: Quiet NAN (not trapped) : x = NAN; y = x + 1;\n");
+        break;
+
+        case 6:
+        set_snan(&x);
+        print_hexval(x);
+        y = x + 1;
+        printf("3: Signaling NAN (trapped) : x = snan; y = x + 1;\n");
+        break;
+
+        default:
+        y = -1;
+        printf("Usage : tst_fpe <fpe_id>\n");
+        printf("        fpe_id in [1-5]\n\n");
+        exit(0);
+    }
+    return y;
+}
+
+int main(int argc, char** argv) 
+{
+    int fpe_id = 1;
+
+    if (argc > 1)        
+        fpe_id = atoi(argv[1]);
+
+    printf("Exception handling disabled : \n");
+
+    double y = create_exception(fpe_id);
+
+    printf("Result : y = %g\n",y);
+
+    printf("\n");
+    printf("Enabling exception handling\n");  
+    enable_floating_point_exceptions();
+
+    y = create_exception(fpe_id);
+    printf("Result : y = %g\n",y);
+}
+

src/CMakeLists.txt

@@ -9,7 +9,10 @@ add_library(forestclaw_c OBJECT
   fclaw_math.c 
   fclaw_timer.c 
   fclaw_mpi.c 
-  fclaw2d_block.c 
+  fclaw_fpe_default.c
+  fclaw_fpe_linux.c
+  fclaw_fpe_macos.c
+  fclaw2d_block.c   
   fclaw2d_options.c 
   fclaw2d_global.c
   fclaw2d_forestclaw.c 

src/fclaw2d_vtable.c

@@ -27,6 +27,8 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include <fclaw2d_output.h>
 #include <fclaw2d_global.h>
 
+#include <fclaw_fpe.h>
+
 #include <forestclaw2d.h>
 
 static fclaw2d_vtable_t s_vt;
@@ -64,6 +66,16 @@ void fclaw2d_after_regrid(fclaw2d_global_t *glob)
     }
 }
 
+void fclaw_fpe_handling(fclaw2d_global_t *glob)
+{
+    fclaw2d_vtable_t *fclaw_vt = fclaw2d_vt();
+    if (fclaw_vt->fpe_enable != NULL)
+    {
+        fclaw_vt->fpe_enable(glob);
+    }
+}
+
+
 /* Initialize any settings that can be set here */
 void fclaw2d_vtable_initialize()
 {
@@ -86,5 +98,8 @@ void fclaw2d_vtable_initialize()
     /* Defaults for output */
     vt->output_frame               = NULL;
 
+    vt->fpe_enable                 = fclaw_fpe_handling_default;
+    vt->fpe_signal_handler         = fclaw_fpe_signal_handler_default;
+
     vt->is_set = 1;
 }

src/fclaw2d_vtable.h

@@ -29,9 +29,6 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #ifdef __cplusplus
 extern "C"
 {
-#if 0
-}                               /* need this because indent is dumb */
-#endif
 #endif
 
 struct fclaw2d_global;
@@ -42,6 +39,8 @@ void fclaw2d_after_regrid(struct fclaw2d_global *glob);
 
 void fclaw2d_after_init(struct fclaw2d_global *glob);
 
+void fclaw_fpe_handling(struct fclaw2d_global *glob);
+
 
 /* ---------------------------------- Typedefs ---------------------------------------- */  
 typedef void (*fclaw2d_vtable_initialize_t)();
@@ -54,19 +53,26 @@ typedef void (*fclaw2d_after_regrid_t)(struct fclaw2d_global *glob);
 
 typedef void (*fclaw2d_after_initialization_t)(struct fclaw2d_global *glob);
 
+typedef void (*fclaw_fpe_handling_t)(struct fclaw2d_global *glob);
+
+typedef void (*fclaw_fpe_signal_handler_t)(int sig);
+
 /* ------------------------------------ vtable ---------------------------------------- */  
 typedef struct fclaw2d_vtable
 {
 
-	fclaw2d_vtable_initialize_t          vtable_init;
-	fclaw2d_problem_setup_t              problem_setup;
-	fclaw2d_after_initialization_t       after_init;
+	fclaw2d_vtable_initialize_t      vtable_init;
+	fclaw2d_problem_setup_t          problem_setup;
+	fclaw2d_after_initialization_t   after_init;
 
 	/* regridding functions */
-	fclaw2d_after_regrid_t               after_regrid;
+	fclaw2d_after_regrid_t           after_regrid;
 
 	/* Output functions */
-	fclaw2d_output_frame_t               output_frame;
+	fclaw2d_output_frame_t           output_frame;
+
+    fclaw_fpe_handling_t             fpe_enable;
+    fclaw_fpe_signal_handler_t       fpe_signal_handler;
 
 	int is_set;
 
@@ -80,9 +86,6 @@ void fclaw2d_vtable_initialize();
 
 
 #ifdef __cplusplus
-#if 0
-{
-#endif
 }
 #endif
 

src/fclaw_fpe.h

@@ -0,0 +1,57 @@
+/*
+Copyright (c) 2012-2022 Carsten Burstedde, Donna Calhoun, Scott Aiton
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+ * Redistributions of source code must retain the above copyright notice, this
+list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above copyright notice,
+this list of conditions and the following disclaimer in the documentation
+and/or other materials provided with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+/** \file
+ *
+ * Routines for virtualized exception handling
+ *
+ */
+
+#ifndef FCLAW_FPE_H
+#define FCLAW_FPE_H
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
+struct fclaw2d_global;
+
+void fclaw_fpe_handling_default(struct fclaw2d_global *glob);
+
+void fclaw_fpe_signal_handler_default(int sig);
+
+
+/* Thse files are only compiled in optionally */
+void fclaw_fpe_handling_macos(struct fclaw2d_global* glob);
+
+void fclaw_fpe_handling_linux(struct fclaw2d_global* glob);
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif