about summary refs log tree commit diff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile3
-rw-r--r--src/debug.cpp19
-rw-r--r--src/debug.h5
-rw-r--r--src/lapack.cpp2
-rw-r--r--src/lmm.cpp19
-rw-r--r--src/param.cpp11
6 files changed, 46 insertions, 13 deletions
diff --git a/Makefile b/Makefile
index 775a40e..e247a87 100644
--- a/Makefile
+++ b/Makefile
@@ -64,7 +64,8 @@ WITH_LAPACK            =                  # Force linking LAPACK (if OpenBlas la
 WITH_GSLCBLAS          =                  # Force linking gslcblas (if OpenBlas lacks it)
 OPENBLAS_LEGACY        =                  # Using older OpenBlas
 FORCE_STATIC           =                  # Static linking of libraries
-GCC_FLAGS              = -Wall -O3 -std=gnu++11 # extra flags -Wl,--allow-multiple-definition
+# GCC_FLAGS              = -Wall -O3 -std=gnu++11 # extra flags -Wl,--allow-multiple-definition
+GCC_FLAGS              = -Wall -Og -std=gnu++11 # extra flags -Wl,--allow-multiple-definition
 TRAVIS_CI              =                  # used by TRAVIS for testing
 
 GSL_INCLUDE_PATH =
diff --git a/src/debug.cpp b/src/debug.cpp
index 7728d92..37e4e8a 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -142,6 +142,7 @@ void disable_segfpe() {
   if (is_legacy_mode()) return;
   #ifdef __GNUC__
     #if defined(__x86_64__)
+      debug_msg("disable segfpe");
       fedisableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW);
     #endif
   #endif
@@ -228,6 +229,24 @@ gsl_vector *gsl_vector_safe_alloc(size_t n) {
   return v;
 }
 
+double do_gsl_matrix_safe_get(const gsl_matrix * m, const size_t i, const size_t j,
+                              const char *__pretty_function, const char *__file, int __line) {
+  enforce(m);
+  if (!is_legacy_mode() && (is_debug_mode() || is_check_mode() || is_strict_mode())) {
+    auto rows = m->size1;
+    auto cols = m->size2;
+    if (i >= cols || j >= rows) {
+      std::string msg = "Matrix out of bounds (" + std::to_string(cols) + "," + std::to_string(rows) + ") ";
+      msg += std::to_string(i);
+      msg += ",";
+      msg += std::to_string(j);
+      fail_at_msg(__file,__line,msg.c_str());
+    }
+  }
+  return gsl_matrix_get(m,i,j);
+}
+
+
 char *do_strtok_safe(char *tokenize, const char *delimiters, const char *__pretty_function, const char *__file, int __line,
                      const char *infile) {
   auto token = strtok(tokenize,delimiters);
diff --git a/src/debug.h b/src/debug.h
index b08d75d..64ab721 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -57,6 +57,9 @@ int gsl_matrix_safe_memcpy (gsl_matrix *dest, const gsl_matrix *src);
 void gsl_matrix_safe_free (gsl_matrix *v);
 void do_gsl_matrix_safe_free (gsl_matrix *m, const char *__pretty_function, const char *__file, int __line);
 
+double do_gsl_matrix_safe_get (const gsl_matrix * m, const size_t i, const size_t j, const char *__pretty_function, const char *__file, int __line);
+#define gsl_matrix_safe_get(m,i,j) do_gsl_matrix_safe_get(m, i, j,__SHOW_FUNC,__FILE__,__LINE__);
+
 gsl_vector *gsl_vector_safe_alloc(size_t n);
 int gsl_vector_safe_memcpy (gsl_vector *dest, const gsl_vector *src);
 void gsl_vector_safe_free (gsl_vector *v);
@@ -85,7 +88,7 @@ inline void warnfail_at_msg(bool strict, const char *__function, const char *__f
 }
 
 inline void fail_at_msg(const char *__file, int __line, std::string msg) {
-  std::cerr << msg << " in " << __file << " at line " << __line << std::endl;
+  std::cerr << "**** FAILED: " << msg << " in " << __file << " at line " << __line << std::endl;
   exit(1);
 }
 
diff --git a/src/lapack.cpp b/src/lapack.cpp
index 0fc2de4..6e43bd9 100644
--- a/src/lapack.cpp
+++ b/src/lapack.cpp
@@ -222,7 +222,7 @@ void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
     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
+    if (is_check_mode()) enable_segfpe(); // reinstate fast NaN checking
 
     gsl_matrix_transpose(evec);
 
diff --git a/src/lmm.cpp b/src/lmm.cpp
index ed1366d..cdddf40 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -41,8 +41,6 @@
 #include "gsl/gsl_roots.h"
 #include "gsl/gsl_vector.h"
 
-// #include "eigenlib.h"
-
 #include "gzstream.h"
 #include "gemma_io.h"
 #include "fastblas.h"
@@ -226,9 +224,9 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval,
   double p_ab;
   double ps_ab, ps_aw, ps_bw, ps_ww;
 
-  for (size_t p = 0; p <= n_cvt + 1; ++p) {
-    for (size_t a = p + 1; a <= n_cvt + 2; ++a) {
-      for (size_t b = a; b <= n_cvt + 2; ++b) {
+  for (size_t p = 0; p <= n_cvt + 1; ++p) {        // p walks phenotypes + covariates
+    for (size_t a = p + 1; a <= n_cvt + 2; ++a) {  // a in p+1..rest
+      for (size_t b = a; b <= n_cvt + 2; ++b) {    // b in a..rest
         index_ab = GetabIndex(a, b, n_cvt);
         if (p == 0) {
           gsl_vector_const_view Uab_col =
@@ -244,11 +242,12 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval,
           index_bw = GetabIndex(b, p, n_cvt);
           index_ww = GetabIndex(p, p, n_cvt);
 
-          ps_ab = gsl_matrix_get(Pab, p - 1, index_ab);
-          ps_aw = gsl_matrix_get(Pab, p - 1, index_aw);
-          ps_bw = gsl_matrix_get(Pab, p - 1, index_bw);
-          ps_ww = gsl_matrix_get(Pab, p - 1, index_ww);
+          ps_ab = gsl_matrix_safe_get(Pab, p - 1, index_ab);
+          ps_aw = gsl_matrix_safe_get(Pab, p - 1, index_aw);
+          ps_bw = gsl_matrix_safe_get(Pab, p - 1, index_bw);
+          ps_ww = gsl_matrix_safe_get(Pab, p - 1, index_ww);
 
+          assert(ps_ww != 0.0);
           p_ab = ps_ab - ps_aw * ps_bw / ps_ww;
           gsl_matrix_set(Pab, p, index_ab, p_ab);
         }
@@ -291,6 +290,7 @@ void CalcPPab(const size_t n_cvt, const size_t e_mode,
           ps2_bw = gsl_matrix_get(PPab, p - 1, index_bw);
           ps2_ww = gsl_matrix_get(PPab, p - 1, index_ww);
 
+          assert(ps_ww != 0.0);
           p2_ab = ps2_ab + ps_aw * ps_bw * ps2_ww / (ps_ww * ps_ww);
           p2_ab -= (ps_aw * ps2_bw + ps_bw * ps2_aw) / ps_ww;
           gsl_matrix_set(PPab, p, index_ab, p2_ab);
@@ -340,6 +340,7 @@ void CalcPPPab(const size_t n_cvt, const size_t e_mode,
           ps3_bw = gsl_matrix_get(PPPab, p - 1, index_bw);
           ps3_ww = gsl_matrix_get(PPPab, p - 1, index_ww);
 
+          assert(ps_ww != 0.0);
           p3_ab = ps3_ab -
                   ps_aw * ps_bw * ps2_ww * ps2_ww / (ps_ww * ps_ww * ps_ww);
           p3_ab -= (ps_aw * ps3_bw + ps_bw * ps3_aw + ps2_aw * ps2_bw) / ps_ww;
diff --git a/src/param.cpp b/src/param.cpp
index edee79d..1a2b57c 100644
--- a/src/param.cpp
+++ b/src/param.cpp
@@ -1341,6 +1341,8 @@ void compAKtoS(const gsl_matrix *A, const gsl_matrix *K, const size_t n_cvt,
       if (tr_A == 0 || tr_K == 0) {
         d = 0;
       } else {
+        assert((tr_A * tr_K) - 1 != 0);
+        assert(ni_test - n_cvt != 0);
         d = d / (tr_A * tr_K) - 1 / (double)(ni_test - n_cvt);
       }
 
@@ -1369,9 +1371,12 @@ size_t GetabIndex(const size_t a, const size_t b, const size_t n_cvt) {
   }
 
   size_t n = n_cvt + 2;
-  index = (2 * n - l + 2) * (l - 1) / 2 + h - l;
 
+  assert(2+h-1 != 0);
+  index = (2 * n - l + 2) * (l - 1) / 2 + h - l;
   return index;
+
+  // return ( b < a ?  ((2 * n - b + 2) * (b - 1) / 2 + a - b ): ((2 * n - a + 2) * (a - 1) / 2 + b - a) );
 }
 
 // From an existing n by nd (centered) G matrix, compute the d+1 by
@@ -1386,6 +1391,7 @@ void compKtoV(const gsl_matrix *G, gsl_matrix *V) {
   gsl_vector *trKiKj = gsl_vector_alloc(n_vc * (n_vc + 1) / 2);
   gsl_vector *trKi = gsl_vector_alloc(n_vc);
 
+  assert(ni_test != 1);
   double d, tr, r = (double)ni_test / (double)(ni_test - 1);
   size_t t, t_il, t_jm, t_lm, t_im, t_jl, t_ij;
 
@@ -1541,6 +1547,7 @@ void compKtoV(const gsl_matrix *G, gsl_matrix *V) {
     }
   }
 
+  assert(ni_test != 0);
   gsl_matrix_scale(V, 1.0 / pow((double)ni_test, 2));
 
   gsl_matrix_free(KiKj);
@@ -1603,6 +1610,7 @@ void JackknifeAKtoS(const gsl_matrix *W, const gsl_matrix *A,
     }
 
     for (size_t t = 0; t < ni_test; t++) {
+      assert(ni_test != 1);
       sumA[i][t] /= (double)(ni_test - 1);
       sumK[i][t] /= (double)(ni_test - 1);
     }
@@ -1636,6 +1644,7 @@ void JackknifeAKtoS(const gsl_matrix *W, const gsl_matrix *A,
       }
 
       for (size_t t = 0; t < ni_test; t++) {
+        assert(ni_test != 1);
         sumAK[i][j][t] /= (double)(ni_test - 1);
       }
generated by cgit-pink 1.4.1 (git 2.36.1) at 2025-10-07 20:31:43 +0000