3 files changed, 55 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
diff --git a/src/debug.h b/src/debug.h
index 67764df..ced8e72 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -45,6 +45,9 @@ bool is_quiet_mode();
 bool is_issue(uint issue);
 bool is_legacy_mode();
 
+void enable_segfpe();
+void disable_segfpe();
+
 #define check_int_mult_overflow(m,n) \
   { auto x = m * n;                                      \
     enforce_msg(x / m == n, "multiply integer overflow"); }
diff --git a/src/gemma.cpp b/src/gemma.cpp
index 809cd8e..a183edc 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -1624,6 +1624,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
   clock_t time_begin, time_start;
   time_begin = clock();
 
+  if (is_check_mode()) enable_segfpe(); // fast NaN checking
+
   // Read Files.
   cout << "Reading Files ... " << endl;
   cPar.ReadFiles();
@@ -1882,6 +1884,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     enforce_msg(G, "allocate G"); // just to be sure
 
     time_start = clock();
+
     cPar.CalcKin(G);
 
     cPar.time_G = (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0);
@@ -1909,6 +1912,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
     VARCOV cVarcov;
     cVarcov.CopyFromParam(cPar);
 
+    if (is_check_mode()) disable_segfpe(); // fast NaN checking
+
     if (!cPar.file_bfile.empty()) {
       cVarcov.AnalyzePlink();
     } else {
@@ -2022,6 +2027,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
     VARCOV cVarcov;
     cVarcov.CopyFromParam(cPar);
 
+    if (is_check_mode()) disable_segfpe(); // fast NaN checking
+
     if (!cPar.file_bfile.empty()) {
       cVarcov.AnalyzePlink();
     } else {
@@ -2048,6 +2055,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
       gsl_vector_view Y_col = gsl_matrix_column(Y, 0);
 
+      if (is_check_mode()) disable_segfpe(); // fast NaN checking
+
       if (!cPar.file_gene.empty()) {
         cLm.AnalyzeGene(W,
                         &Y_col.vector); // y is the predictor, not the phenotype
@@ -2754,6 +2763,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
           LMM cLmm;
           cLmm.CopyFromParam(cPar);
 
+          if (is_check_mode()) disable_segfpe(); // fast NaN checking
+
           gsl_vector_view Y_col = gsl_matrix_column(Y, 0);
           gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0);
 
@@ -2770,6 +2781,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
           }
           else {
             // BIMBAM analysis
+
             if (cPar.file_gxe.empty()) {
               cLmm.AnalyzeBimbam(U, eval, UtW, &UtY_col.vector, W,
                                  &Y_col.vector, cPar.setGWASnps);
@@ -2785,6 +2797,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
           MVLMM cMvlmm;
           cMvlmm.CopyFromParam(cPar);
 
+          if (is_check_mode()) disable_segfpe(); // fast NaN checking
+
           if (!cPar.file_bfile.empty()) {
             if (cPar.file_gxe.empty()) {
               cMvlmm.AnalyzePlink(U, eval, UtW, UtY);
@@ -3107,6 +3121,9 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) {
   outfile << "##" << endl;
   outfile << "## GEMMA Version    = " << version << " (" << date << ")" << endl;
   outfile << "## Build profile    = " << GEMMA_PROFILE << endl ;
+  #ifdef __GNUC__
+  outfile << "## GCC version      = " << __GNUC__ << "." << __GNUC_MINOR__ << "." << __GNUC_PATCHLEVEL__ << endl;
+  #endif
   outfile << "## GSL Version      = " << GSL_VERSION << endl;
   outfile << "## Eigen Version    = " << EIGEN_WORLD_VERSION << "." << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION << endl;
 #ifdef OPENBLAS