Add floating point hardware checking for Intel on GNU compilers

When using the -check function (the default) it is enabled for Kinship
computation and LM/LMM up to individual SNP computation. This means
there can no longer be NaN values for matrices that are reused for every
SNP, but it is possible to have NaN for individual SNPs.

Fixes #161

1 files changed, 35 insertions, 0 deletions

diff --git a/src/debug.cpp b/src/debug.cpp
index 4e58d5d..0b9c19d 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -26,6 +26,7 @@
 #include <string>
 #include <sys/stat.h>
 #include <vector>
+#include <fenv.h>
 
 #include "gsl/gsl_blas.h"
 #include "gsl/gsl_cdf.h"
@@ -59,6 +60,40 @@ bool is_quiet_mode() { return debug_quiet; };
 bool is_issue(uint issue) { return issue == debug_issue; };
 bool is_legacy_mode() { return debug_legacy; };
 
+#include <stdio.h>
+#include <sys/types.h>
+#include <unistd.h>
+#include <signal.h>
+
+void sighandler(int signum)
+{
+  cout << R"(
+FATAL ERROR: GEMMA caused a floating point error which suggests machine boundaries were reached.
+
+You can disable floating point tests with the -no-check switch (use at your own risk!)
+)" << endl;
+  signal(signum, SIG_DFL);
+  kill(getpid(), signum); // should force a core dump
+}
+
+/*
+   Force the floating point processor to throw an exception when the result of
+   a double/float computation is overflow, underflow, NaN or inf. In principle
+   this is an Intel hardware feature that does not slow down computations.
+*/
+
+void enable_segfpe() {
+  #ifdef __GNUC__
+  signal(SIGFPE, sighandler);
+  feenableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW);
+  #endif
+}
+
+void disable_segfpe() {
+  #ifdef __GNUC__
+  fedisableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW);
+  #endif
+}
 
 /*
   Helper function to make sure gsl allocations do their job because