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

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