segfpe handling

4 files changed, 27 insertions, 11 deletions

diff --git a/src/gemma.cpp b/src/gemma.cpp
index a183edc..78633e8 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -1624,7 +1624,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
   clock_t time_begin, time_start;
   time_begin = clock();
 
-  if (is_check_mode()) enable_segfpe(); // fast NaN checking
+  if (is_check_mode()) enable_segfpe(); // fast NaN checking by default
 
   // Read Files.
   cout << "Reading Files ... " << endl;
@@ -1912,7 +1912,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     VARCOV cVarcov;
     cVarcov.CopyFromParam(cPar);
 
-    if (is_check_mode()) disable_segfpe(); // fast NaN checking
+    if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now
 
     if (!cPar.file_bfile.empty()) {
       cVarcov.AnalyzePlink();
@@ -2027,7 +2027,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     VARCOV cVarcov;
     cVarcov.CopyFromParam(cPar);
 
-    if (is_check_mode()) disable_segfpe(); // fast NaN checking
+    if (is_check_mode()) disable_segfpe(); // fast NaN checking for now
 
     if (!cPar.file_bfile.empty()) {
       cVarcov.AnalyzePlink();
@@ -2055,7 +2055,7 @@ 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 (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now
 
       if (!cPar.file_gene.empty()) {
         cLm.AnalyzeGene(W,
@@ -2763,7 +2763,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
           LMM cLmm;
           cLmm.CopyFromParam(cPar);
 
-          if (is_check_mode()) disable_segfpe(); // fast NaN checking
+          if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now
 
           gsl_vector_view Y_col = gsl_matrix_column(Y, 0);
           gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0);
@@ -2790,14 +2790,13 @@ void GEMMA::BatchRun(PARAM &cPar) {
                                     &Y_col.vector, env);
             }
           }
-
           cLmm.WriteFiles();
           cLmm.CopyToParam(cPar);
         } else {
           MVLMM cMvlmm;
           cMvlmm.CopyFromParam(cPar);
 
-          if (is_check_mode()) disable_segfpe(); // fast NaN checking
+          if (is_check_mode()) disable_segfpe(); // disable fast NaN checking
 
           if (!cPar.file_bfile.empty()) {
             if (cPar.file_gxe.empty()) {
diff --git a/src/gemma_io.cpp b/src/gemma_io.cpp
index 818c5e8..fad3e47 100644
--- a/src/gemma_io.cpp
+++ b/src/gemma_io.cpp
@@ -1483,12 +1483,14 @@ bool BimbamKin(const string file_geno, const set<string> ksnps,
   // and transpose
   // FIXME: the following is very slow
 
+  debug_msg("begin transpose");
   for (size_t i = 0; i < ni_total; ++i) {
     for (size_t j = 0; j < i; ++j) {
       d = gsl_matrix_get(matrix_kin, j, i);
       gsl_matrix_set(matrix_kin, i, j, d);
     }
   }
+  debug_msg("end transpose");
   // GSL is faster - and there are even faster methods
   // enforce_gsl(gsl_matrix_transpose(matrix_kin));
 
diff --git a/src/lapack.cpp b/src/lapack.cpp
index 187daa8..b309a1b 100644
--- a/src/lapack.cpp
+++ b/src/lapack.cpp
@@ -193,6 +193,8 @@ void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
     double WORK_temp[1];
     int IWORK_temp[1];
 
+    if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now
+
     // DSYEVR - computes selected eigenvalues and, optionally,
     // eigenvectors of a real symmetric matrix
     // Here compute both (JOBZ is V), all eigenvalues (RANGE is A)
@@ -200,7 +202,12 @@ void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
     dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU,
             &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK_temp,
             &LWORK, IWORK_temp, &LIWORK, &INFO);
-    enforce_msg(INFO != 0, "lapack_eigen_symmv failed");
+    // If info = 0, the execution is successful.
+    // If info = -i, the i-th parameter had an illegal value.
+    // If info = i, an internal error has occurred.
+
+    if (INFO != 0) cerr << "ERROR: value of INFO is " << INFO;
+    enforce_msg(INFO == 0, "lapack_eigen_symmv failed");
     LWORK = (int)WORK_temp[0];
     LIWORK = (int)IWORK_temp[0];
 
@@ -210,7 +217,10 @@ void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
     dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU,
             &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK, &LWORK,
             IWORK, &LIWORK, &INFO);
-    enforce_msg(INFO != 0, "lapack_eigen_symmv failed");
+    if (INFO != 0) cerr << "ERROR: value of INFO is " << INFO;
+    enforce_msg(INFO == 0, "lapack_eigen_symmv failed");
+
+    if (is_check_mode()) disable_segfpe(); // reinstate fast NaN checking
 
     gsl_matrix_transpose(evec);
 
diff --git a/src/lmm.cpp b/src/lmm.cpp
index acd9667..7e4c3e1 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -97,8 +97,9 @@ void LMM::CopyToParam(PARAM &cPar) {
 }
 
 void LMM::WriteFiles() {
-
   string file_str;
+  debug_msg("LMM::WriteFiles");
+
   file_str = path_out + "/" + file_out;
   file_str += ".assoc.txt";
 
@@ -1319,6 +1320,7 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
 
   auto batch_compute = [&](size_t l) { // using a C++ closure
     // Compute SNPs in batch, note the computations are independent per SNP
+    debug_msg("enter batch_compute");
     gsl_matrix_view Xlarge_sub = gsl_matrix_submatrix(Xlarge, 0, 0, inds, l);
     gsl_matrix_view UtXlarge_sub =
         gsl_matrix_submatrix(UtXlarge, 0, 0, inds, l);
@@ -1366,6 +1368,7 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
                       p_wald, p_lrt, p_score, logl_H1};
       sumStat.push_back(SNPs);
     }
+    debug_msg("exit batch_compute");
   };
 
   const auto num_snps = indicator_snp.size();
@@ -1441,9 +1444,11 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
     gsl_vector_safe_memcpy(&Xlarge_col.vector, x);
     c++; // count SNPs going in
 
-    if (c % msize == 0)
+    if (c % msize == 0) {
       batch_compute(msize);
+    }
   }
+
   batch_compute(c % msize);
   ProgressBar("Reading SNPs", num_snps - 1, num_snps - 1);
   // cout << "Counted SNPs " << c << " sumStat " << sumStat.size() << endl;