diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/debug.cpp | 35 | ||||
| -rw-r--r-- | src/debug.h | 3 | ||||
| -rw-r--r-- | src/gemma.cpp | 17 |
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 |