Nans: introducing checking on mem free

6 files changed, 357 insertions, 266 deletions

diff --git a/src/debug.cpp b/src/debug.cpp
index d133344..a00190d 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -56,6 +56,55 @@ gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols) {
   return m;
 }
 
+int gsl_matrix_safe_memcpy (gsl_matrix *dest, const gsl_matrix *src) {
+  enforce(dest->size1 == src->size1);
+  enforce(dest->size2 == src->size2);
+  return gsl_matrix_memcpy(dest,src);
+}
+
+void do_gsl_matrix_safe_free (gsl_matrix *m, const char *__pretty_function, const char *__file, int __line) {
+  enforce(m);
+  if (is_check_mode() && is_debug_mode()) {
+    bool has_NaN = has_nan(m);
+    bool has_Inf = has_inf(m);
+    if (has_NaN || has_Inf) {
+      std::string msg = "Matrix (size ";
+      msg += std::to_string(m->size1);
+      msg += "x";
+      msg += std::to_string(m->size2);
+      msg += ")";
+      if (has_Inf)
+        warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,(msg+" contains Infinite on free!").c_str());
+      if (has_NaN)
+        warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,(msg+" contains NaN on free!").c_str());
+    }
+  }
+  return gsl_matrix_free(m);
+}
+
+int gsl_vector_safe_memcpy (gsl_vector *dest, const gsl_vector *src) {
+  enforce(dest->size == src->size);
+  return gsl_vector_memcpy(dest,src);
+}
+
+void do_gsl_vector_safe_free (gsl_vector *v, const char *__pretty_function, const char *__file, int __line) {
+  enforce(v);
+  if (is_check_mode() && is_debug_mode()) {
+    bool has_NaN = has_nan(v);
+    bool has_Inf = has_inf(v);
+    if (has_NaN || has_Inf) {
+      std::string msg = "Vector (size ";
+      msg += std::to_string(v->size);
+      msg += ")";
+      if (has_Inf)
+        warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,(msg+" contains Infinite on free!").c_str());
+      if (has_NaN)
+        warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,(msg+" contains NaN on free!").c_str());
+    }
+  }
+  return gsl_vector_free(v);
+}
+
 /*
   Helper function to make sure gsl allocations do their job because
   gsl_vector_alloc does not initiatize values (behaviour that changed
@@ -81,7 +130,7 @@ char *do_strtok_safe(char *tokenize, const char *delimiters, const char *__prett
 
 // Helper function called by macro validate_K(K, check). K is validated
 // unless -no-check option is used.
-void do_validate_K(const gsl_matrix *K, const char *__file, int __line) {
+void do_validate_K(const gsl_matrix *K, const char *__pretty_function, const char *__file, int __line) {
   if (is_check_mode()) {
     // debug_msg("Validating K");
     auto eigenvalues = getEigenValues(K);
@@ -95,13 +144,13 @@ void do_validate_K(const gsl_matrix *K, const char *__file, int __line) {
     if (!isMatrixIllConditioned(eigenvalues))
       warning_at_msg(__file,__line,"K is ill conditioned!");
     if (!isMatrixSymmetric(K))
-      warnfail_at_msg(is_strict_mode(),__file,__line,"K is not symmetric!" );
+      warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,"K is not symmetric!" );
     const bool negative_values = has_negative_values_but_one(eigenvalues);
     if (negative_values) {
       warning_at_msg(__file,__line,"K has more than one negative eigenvalues!");
     }
     if (count_small>1 && negative_values && !isMatrixPositiveDefinite(K))
-      warnfail_at_msg(is_strict_mode(),__file,__line,"K is not positive definite!");
+      warnfail_at_msg(is_strict_mode(),__pretty_function,__file,__line,"K is not positive definite!");
     gsl_vector_free(eigenvalues);
   }
 }
diff --git a/src/debug.h b/src/debug.h
index 6b82b3b..69b0a7c 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -26,25 +26,32 @@ bool is_issue(uint issue);
 bool is_legacy_mode();
 
 gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols);
+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);
+
 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);
+void do_gsl_vector_safe_free (gsl_vector *v, const char *__pretty_function, const char *__file, int __line);
 
 char *do_strtok_safe(char *tokenize, const char *delimiters, const char *__pretty_function, const char *__file, int __line);
 #define strtok_safe(string,delimiters) do_strtok_safe(string,delimiters,__PRETTY_FUNCTION__,__FILE__,__LINE__)
 
 // Validation routines
-void do_validate_K(const gsl_matrix *K, const char *__file, int __line);
+void do_validate_K(const gsl_matrix *K, const char*__pretty_func, const char *__file, int __line);
 
 #define ROUND(f) round(f * 10000.)/10000
-#define validate_K(K) do_validate_K(K,__FILE__,__LINE__)
+#define validate_K(K) do_validate_K(K,__PRETTY_FUNCTION__,__FILE__,__LINE__)
 
 #define warning_at_msg(__file,__line,msg) cerr << "**** WARNING: " << msg << " in " << __file << " at line " << __line << endl;
 
-inline void warnfail_at_msg(bool strict, const char *__file, int __line, const char *msg) {
+inline void warnfail_at_msg(bool strict, const char *__function, const char *__file, int __line, const char *msg) {
   if (strict)
     std::cerr << "**** STRICT FAIL: ";
   else
     std::cerr << "**** WARNING: ";
-  std::cerr << msg << " in " << __file << " at line " << __line << std::endl;
+  std::cerr << msg << " in " << __file << " at line " << __line << " in " << __function << std::endl;
   if (strict)
     exit(1);
 }
@@ -127,5 +134,9 @@ inline void __enforce_fail(const char *__assertion, const char *__file,
     enforce_msg(stat(fn.c_str(), &fileInfo) == 0,                              \
                 ((std::string(__STRING(fn)) + " " + fn + ": " + msg).c_str()));
 
+#define gsl_matrix_safe_free(m) \
+  do_gsl_matrix_safe_free(m,__ASSERT_FUNCTION,__FILE__,__LINE__);
+#define gsl_vector_safe_free(v) \
+  do_gsl_vector_safe_free(v,__ASSERT_FUNCTION,__FILE__,__LINE__);
 
 #endif
diff --git a/src/gemma.cpp b/src/gemma.cpp
index 2c54672..8977255 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -1673,8 +1673,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
       // read u
       cPRDT.AddBV(G, u_hat, y_prdt);
 
-      gsl_matrix_free(G);
-      gsl_vector_free(u_hat);
+      gsl_matrix_safe_free(G);
+      gsl_vector_safe_free(u_hat);
     }
 
     // add beta
@@ -1701,7 +1701,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
     cPRDT.WriteFiles(y_prdt);
 
-    gsl_vector_free(y_prdt);
+    gsl_vector_safe_free(y_prdt);
   }
 
   // Prediction with kinship matrix only; for one or more phenotypes
@@ -1790,8 +1790,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
       gsl_matrix_add(H_full, G_full);
 
       // free matrices
-      gsl_vector_free(beta);
-      gsl_vector_free(se_beta);
+      gsl_vector_safe_free(beta);
+      gsl_vector_safe_free(se_beta);
     } else {
       gsl_matrix *Vg = gsl_matrix_safe_alloc(cPar.n_ph, cPar.n_ph);
       gsl_matrix *Ve = gsl_matrix_safe_alloc(cPar.n_ph, cPar.n_ph);
@@ -1838,10 +1838,10 @@ void GEMMA::BatchRun(PARAM &cPar) {
       }
 
       // free matrices
-      gsl_matrix_free(Vg);
-      gsl_matrix_free(Ve);
-      gsl_matrix_free(B);
-      gsl_matrix_free(se_B);
+      gsl_matrix_safe_free(Vg);
+      gsl_matrix_safe_free(Ve);
+      gsl_matrix_safe_free(B);
+      gsl_matrix_safe_free(se_B);
     }
 
     PRDT cPRDT;
@@ -1855,19 +1855,19 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
     cPRDT.WriteFiles(Y_full);
 
-    gsl_matrix_free(Y);
-    gsl_matrix_free(W);
-    gsl_matrix_free(G);
-    gsl_matrix_free(U);
-    gsl_matrix_free(UtW);
-    gsl_matrix_free(UtY);
-    gsl_vector_free(eval);
-
-    gsl_matrix_free(Y_full);
-    gsl_matrix_free(Y_hat);
-    gsl_matrix_free(W_full);
-    gsl_matrix_free(G_full);
-    gsl_matrix_free(H_full);
+    gsl_matrix_safe_free(Y);
+    gsl_matrix_safe_free(W);
+    gsl_matrix_safe_free(G);
+    gsl_matrix_safe_free(U);
+    gsl_matrix_safe_free(UtW);
+    gsl_matrix_safe_free(UtY);
+    gsl_vector_safe_free(eval);
+
+    gsl_matrix_safe_free(Y_full);
+    gsl_matrix_safe_free(Y_hat);
+    gsl_matrix_safe_free(W_full);
+    gsl_matrix_safe_free(G_full);
+    gsl_matrix_safe_free(H_full);
   }
 
   // Generate Kinship matrix (optionally using LOCO)
@@ -1895,7 +1895,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
       cPar.WriteMatrix(G, "sXX");
     }
 
-    gsl_matrix_free(G);
+    gsl_matrix_safe_free(G);
   }
 
   // Compute the LDSC weights (not implemented yet)
@@ -1959,14 +1959,14 @@ void GEMMA::BatchRun(PARAM &cPar) {
     cPar.WriteVector(ns, "size");
     cPar.WriteVar("snps");
 
-    gsl_matrix_free(S);
-    gsl_vector_free(ns);
+    gsl_matrix_safe_free(S);
+    gsl_vector_safe_free(ns);
 
-    gsl_matrix_free(A);
-    gsl_matrix_free(K);
+    gsl_matrix_safe_free(A);
+    gsl_matrix_safe_free(K);
 
-    gsl_vector_free(y);
-    gsl_matrix_free(K);
+    gsl_vector_safe_free(y);
+    gsl_matrix_safe_free(K);
   }
 
   // Compute the q vector, that is used for variance component estimation using
@@ -2008,9 +2008,9 @@ void GEMMA::BatchRun(PARAM &cPar) {
     cPar.WriteMatrix(Vq, "Vq");
     cPar.WriteVector(q, "q");
     cPar.WriteVector(s, "size");
-    gsl_matrix_free(Vq);
-    gsl_vector_free(q);
-    gsl_vector_free(s);
+    gsl_matrix_safe_free(Vq);
+    gsl_vector_safe_free(q);
+    gsl_vector_safe_free(s);
   }
 
   // Calculate SNP covariance.
@@ -2057,8 +2057,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
       cLm.CopyToParam(cPar);
     }
     // release all matrices and vectors
-    gsl_matrix_free(Y);
-    gsl_matrix_free(W);
+    gsl_matrix_safe_free(Y);
+    gsl_matrix_safe_free(W);
   }
 
   // VC estimation with one or multiple kinship matrices
@@ -2186,15 +2186,15 @@ void GEMMA::BatchRun(PARAM &cPar) {
       cPar.WriteVector(q, "q");
       cPar.WriteVector(s, "size");
 
-      gsl_matrix_free(S);
-      gsl_matrix_free(Vq);
-      gsl_vector_free(q);
-      gsl_vector_free(s);
+      gsl_matrix_safe_free(S);
+      gsl_matrix_safe_free(Vq);
+      gsl_vector_safe_free(q);
+      gsl_vector_safe_free(s);
 
-      gsl_matrix_free(A);
-      gsl_matrix_free(K);
-      gsl_vector_free(y);
-      gsl_matrix_free(W);
+      gsl_matrix_safe_free(A);
+      gsl_matrix_safe_free(K);
+      gsl_vector_safe_free(y);
+      gsl_matrix_safe_free(W);
     } else if (!cPar.file_study.empty() || !cPar.file_mstudy.empty()) {
       if (!cPar.file_study.empty()) {
         string sfile = cPar.file_study + ".size.txt";
@@ -2270,7 +2270,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
       assert(!has_nan(cPar.v_se_pve));
 
       gsl_vector_view s_sub = gsl_vector_subvector(s, 0, cPar.n_vc);
-      gsl_vector_memcpy(&s_sub.vector, s_ref);
+      gsl_vector_safe_memcpy(&s_sub.vector, s_ref);
       gsl_vector_set(s, cPar.n_vc, cPar.ni_ref);
 
       cPar.WriteMatrix(S, "S");
@@ -2278,14 +2278,14 @@ void GEMMA::BatchRun(PARAM &cPar) {
       cPar.WriteVector(q, "q");
       cPar.WriteVector(s, "size");
 
-      gsl_matrix_free(S);
-      gsl_matrix_free(Vq);
-      // gsl_matrix_free (V);
-      // gsl_matrix_free (Vslope);
-      gsl_vector_free(q);
-      gsl_vector_free(s_study);
-      gsl_vector_free(s_ref);
-      gsl_vector_free(s);
+      gsl_matrix_safe_free(S);
+      gsl_matrix_safe_free(Vq);
+      // gsl_matrix_safe_free (V);
+      // gsl_matrix_safe_free (Vslope);
+      gsl_vector_safe_free(q);
+      gsl_vector_safe_free(s_study);
+      gsl_vector_safe_free(s_ref);
+      gsl_vector_safe_free(s);
     } else {
       gsl_matrix *Y = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_ph);
       gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt);
@@ -2462,7 +2462,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
                cPar.mindicator_snp, vec_cat, w1, z, Xz);
     }
     if (cPar.a_mode == 66) {
-      gsl_matrix_memcpy(XWz, Xz);
+      gsl_matrix_safe_memcpy(XWz, Xz);
     } else if (cPar.a_mode == 67) {
       cout << "Calculating XWz ... " << endl;
 
@@ -2507,17 +2507,17 @@ void GEMMA::BatchRun(PARAM &cPar) {
              cPar.v_se_sigma2, cPar.v_enrich, cPar.v_se_enrich);
     assert(!has_nan(cPar.v_se_pve));
 
-    gsl_matrix_free(S);
-    gsl_matrix_free(Svar);
-    gsl_vector_free(s_ref);
-
-    gsl_matrix_free(Xz);
-    gsl_matrix_free(XWz);
-    gsl_matrix_free(XtXWz);
-    gsl_vector_free(w);
-    gsl_vector_free(w1);
-    gsl_vector_free(z);
-    gsl_vector_free(s_vec);
+    gsl_matrix_safe_free(S);
+    gsl_matrix_safe_free(Svar);
+    gsl_vector_safe_free(s_ref);
+
+    gsl_matrix_safe_free(Xz);
+    gsl_matrix_safe_free(XWz);
+    gsl_matrix_safe_free(XtXWz);
+    gsl_vector_safe_free(w);
+    gsl_vector_safe_free(w1);
+    gsl_vector_safe_free(z);
+    gsl_vector_safe_free(s_vec);
   }
 
   // LMM or mvLMM or Eigen-Decomposition
@@ -2716,7 +2716,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
           gsl_vector *y_hat = gsl_vector_safe_alloc(Y->size1);
 
           // obtain Utu and Ute
-          gsl_vector_memcpy(y_hat, &UtY_col.vector);
+          gsl_vector_safe_memcpy(y_hat, &UtY_col.vector);
           gsl_blas_dgemv(CblasNoTrans, -1.0, UtW, &beta.vector, 1.0, y_hat);
 
           double d, u, e;
@@ -2737,9 +2737,9 @@ void GEMMA::BatchRun(PARAM &cPar) {
           cPar.WriteVector(u_hat, "residU");
           cPar.WriteVector(e_hat, "residE");
 
-          gsl_vector_free(u_hat);
-          gsl_vector_free(e_hat);
-          gsl_vector_free(y_hat);
+          gsl_vector_safe_free(u_hat);
+          gsl_vector_safe_free(e_hat);
+          gsl_vector_safe_free(y_hat);
         }
       }
 
@@ -2802,16 +2802,16 @@ void GEMMA::BatchRun(PARAM &cPar) {
     }
 
     // release all matrices and vectors
-    gsl_matrix_free(Y);
-    gsl_matrix_free(W);
-    gsl_matrix_free(B);
-    gsl_matrix_free(se_B);
-    gsl_matrix_free(G);
-    gsl_matrix_free(U);
-    gsl_matrix_free(UtW);
-    gsl_matrix_free(UtY);
-    gsl_vector_free(eval);
-    gsl_vector_free(env);
+    gsl_matrix_safe_free(Y);
+    gsl_matrix_safe_free(W);
+    gsl_matrix_safe_free(B);
+    gsl_matrix_safe_free(se_B);
+    gsl_matrix_safe_free(G);
+    gsl_matrix_safe_free(U);
+    gsl_matrix_safe_free(UtW);
+    gsl_matrix_safe_free(UtY);
+    gsl_vector_safe_free(eval);
+    gsl_vector_safe_free(env);
   }
 
   // BSLMM
@@ -2912,15 +2912,15 @@ void GEMMA::BatchRun(PARAM &cPar) {
       }
 
       // release all matrices and vectors
-      gsl_matrix_free(G);
-      gsl_matrix_free(U);
-      gsl_matrix_free(UtW);
-      gsl_vector_free(eval);
-      gsl_vector_free(Uty);
+      gsl_matrix_safe_free(G);
+      gsl_matrix_safe_free(U);
+      gsl_matrix_safe_free(UtW);
+      gsl_vector_safe_free(eval);
+      gsl_vector_safe_free(Uty);
     }
-    gsl_matrix_free(W);
-    gsl_vector_free(y);
-    gsl_matrix_free(UtX);
+    gsl_matrix_safe_free(W);
+    gsl_vector_safe_free(y);
+    gsl_matrix_safe_free(UtX);
   }
 
   // BSLMM-DAP
@@ -3016,16 +3016,16 @@ void GEMMA::BatchRun(PARAM &cPar) {
         cBslmmDap.CopyToParam(cPar);
 
         // release all matrices and vectors
-        gsl_matrix_free(G);
-        gsl_matrix_free(U);
-        gsl_matrix_free(UtW);
-        gsl_vector_free(eval);
-        gsl_vector_free(Uty);
+        gsl_matrix_safe_free(G);
+        gsl_matrix_safe_free(U);
+        gsl_matrix_safe_free(UtW);
+        gsl_vector_safe_free(eval);
+        gsl_vector_safe_free(Uty);
       }
 
-      gsl_matrix_free(W);
-      gsl_vector_free(y);
-      gsl_matrix_free(UtX);
+      gsl_matrix_safe_free(W);
+      gsl_vector_safe_free(y);
+      gsl_matrix_safe_free(UtX);
     } else if (cPar.a_mode == 15) {
       // perform EM algorithm and estimate parameters
       vector<string> vec_rs;
@@ -3074,7 +3074,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
       cPar.time_opt = (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0);
       cBslmmDap.CopyToParam(cPar);
 
-      gsl_matrix_free(Ac);
+      gsl_matrix_safe_free(Ac);
       gsl_matrix_int_free(Ad);
       gsl_vector_int_free(dlevel);
     } else {
diff --git a/src/lmm.cpp b/src/lmm.cpp
index 683858d..1a2614b 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -367,12 +367,12 @@ double LogL_f(double l, void *params) {
   gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size);
   gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
 
-  gsl_vector_memcpy(v_temp, p->eval);
+  gsl_vector_safe_memcpy(v_temp, p->eval);
   gsl_vector_scale(v_temp, l);
   if (p->e_mode == 0) {
     gsl_vector_set_all(Hi_eval, 1.0);
   } else {
-    gsl_vector_memcpy(Hi_eval, v_temp);
+    gsl_vector_safe_memcpy(Hi_eval, v_temp);
   }
   gsl_vector_add_constant(v_temp, 1.0);
   gsl_vector_div(Hi_eval, v_temp);
@@ -391,9 +391,9 @@ double LogL_f(double l, void *params) {
   double P_yy = gsl_matrix_get(Pab, nc_total, index_yy);
   f = c - 0.5 * logdet_h - 0.5 * (double)ni_test * log(P_yy);
 
-  gsl_matrix_free(Pab);
-  gsl_vector_free(Hi_eval);
-  gsl_vector_free(v_temp);
+  gsl_matrix_safe_free(Pab); // FIXME
+  gsl_vector_safe_free(Hi_eval);
+  gsl_vector_safe_free(v_temp);
   return f;
 }
 
@@ -419,17 +419,17 @@ double LogL_dev1(double l, void *params) {
   gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
   gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
 
-  gsl_vector_memcpy(v_temp, p->eval);
+  gsl_vector_safe_memcpy(v_temp, p->eval);
   gsl_vector_scale(v_temp, l);
   if (p->e_mode == 0) {
     gsl_vector_set_all(Hi_eval, 1.0);
   } else {
-    gsl_vector_memcpy(Hi_eval, v_temp);
+    gsl_vector_safe_memcpy(Hi_eval, v_temp);
   }
   gsl_vector_add_constant(v_temp, 1.0);
   gsl_vector_div(Hi_eval, v_temp);
 
-  gsl_vector_memcpy(HiHi_eval, Hi_eval);
+  gsl_vector_safe_memcpy(HiHi_eval, Hi_eval);
   gsl_vector_mul(HiHi_eval, Hi_eval);
 
   gsl_vector_set_all(v_temp, 1.0);
@@ -451,11 +451,11 @@ double LogL_dev1(double l, void *params) {
   double yPKPy = (P_yy - PP_yy) / l;
   dev1 = -0.5 * trace_HiK + 0.5 * (double)ni_test * yPKPy / P_yy;
 
-  gsl_matrix_free(Pab);
-  gsl_matrix_free(PPab);
-  gsl_vector_free(Hi_eval);
-  gsl_vector_free(HiHi_eval);
-  gsl_vector_free(v_temp);
+  gsl_matrix_safe_free(Pab);   // FIXME: may contain NaN
+  gsl_matrix_safe_free(PPab);  // FIXME: may contain NaN
+  gsl_vector_safe_free(Hi_eval);
+  gsl_vector_safe_free(HiHi_eval);
+  gsl_vector_safe_free(v_temp);
 
   return dev1;
 }
@@ -484,19 +484,19 @@ double LogL_dev2(double l, void *params) {
   gsl_vector *HiHiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
   gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
 
-  gsl_vector_memcpy(v_temp, p->eval);
+  gsl_vector_safe_memcpy(v_temp, p->eval);
   gsl_vector_scale(v_temp, l);
   if (p->e_mode == 0) {
     gsl_vector_set_all(Hi_eval, 1.0);
   } else {
-    gsl_vector_memcpy(Hi_eval, v_temp);
+    gsl_vector_safe_memcpy(Hi_eval, v_temp);
   }
   gsl_vector_add_constant(v_temp, 1.0);
   gsl_vector_div(Hi_eval, v_temp);
 
-  gsl_vector_memcpy(HiHi_eval, Hi_eval);
+  gsl_vector_safe_memcpy(HiHi_eval, Hi_eval);
   gsl_vector_mul(HiHi_eval, Hi_eval);
-  gsl_vector_memcpy(HiHiHi_eval, HiHi_eval);
+  gsl_vector_safe_memcpy(HiHiHi_eval, HiHi_eval);
   gsl_vector_mul(HiHiHi_eval, Hi_eval);
 
   gsl_vector_set_all(v_temp, 1.0);
@@ -526,13 +526,13 @@ double LogL_dev2(double l, void *params) {
          0.5 * (double)ni_test * (2.0 * yPKPKPy * P_yy - yPKPy * yPKPy) /
              (P_yy * P_yy);
 
-  gsl_matrix_free(Pab);
-  gsl_matrix_free(PPab);
-  gsl_matrix_free(PPPab);
-  gsl_vector_free(Hi_eval);
-  gsl_vector_free(HiHi_eval);
-  gsl_vector_free(HiHiHi_eval);
-  gsl_vector_free(v_temp);
+  gsl_matrix_safe_free(Pab);  // FIXME
+  gsl_matrix_safe_free(PPab);
+  gsl_matrix_safe_free(PPPab);
+  gsl_vector_safe_free(Hi_eval);
+  gsl_vector_safe_free(HiHi_eval);
+  gsl_vector_safe_free(HiHiHi_eval);
+  gsl_vector_safe_free(v_temp);
 
   return dev2;
 }
@@ -561,19 +561,19 @@ void LogL_dev12(double l, void *params, double *dev1, double *dev2) {
   gsl_vector *HiHiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
   gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
 
-  gsl_vector_memcpy(v_temp, p->eval);
+  gsl_vector_safe_memcpy(v_temp, p->eval);
   gsl_vector_scale(v_temp, l);
   if (p->e_mode == 0) {
     gsl_vector_set_all(Hi_eval, 1.0);
   } else {
-    gsl_vector_memcpy(Hi_eval, v_temp);
+    gsl_vector_safe_memcpy(Hi_eval, v_temp);
   }
   gsl_vector_add_constant(v_temp, 1.0);
   gsl_vector_div(Hi_eval, v_temp);
 
-  gsl_vector_memcpy(HiHi_eval, Hi_eval);
+  gsl_vector_safe_memcpy(HiHi_eval, Hi_eval);
   gsl_vector_mul(HiHi_eval, Hi_eval);
-  gsl_vector_memcpy(HiHiHi_eval, HiHi_eval);
+  gsl_vector_safe_memcpy(HiHiHi_eval, HiHi_eval);
   gsl_vector_mul(HiHiHi_eval, Hi_eval);
 
   gsl_vector_set_all(v_temp, 1.0);
@@ -606,13 +606,13 @@ void LogL_dev12(double l, void *params, double *dev1, double *dev2) {
           0.5 * (double)ni_test * (2.0 * yPKPKPy * P_yy - yPKPy * yPKPy) /
               (P_yy * P_yy);
 
-  gsl_matrix_free(Pab);
-  gsl_matrix_free(PPab);
-  gsl_matrix_free(PPPab);
-  gsl_vector_free(Hi_eval);
-  gsl_vector_free(HiHi_eval);
-  gsl_vector_free(HiHiHi_eval);
-  gsl_vector_free(v_temp);
+  gsl_matrix_safe_free(Pab);   // FIXME: may contain NaN
+  gsl_matrix_safe_free(PPab);  // FIXME: may contain NaN
+  gsl_matrix_safe_free(PPPab); // FIXME: may contain NaN
+  gsl_vector_safe_free(Hi_eval);
+  gsl_vector_safe_free(HiHi_eval);
+  gsl_vector_safe_free(HiHiHi_eval);
+  gsl_vector_safe_free(v_temp);
 
   return;
 }
@@ -641,12 +641,12 @@ double LogRL_f(double l, void *params) {
   gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size);
   gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
 
-  gsl_vector_memcpy(v_temp, p->eval);
+  gsl_vector_safe_memcpy(v_temp, p->eval);
   gsl_vector_scale(v_temp, l);
   if (p->e_mode == 0) {
     gsl_vector_set_all(Hi_eval, 1.0);
   } else {
-    gsl_vector_memcpy(Hi_eval, v_temp);
+    gsl_vector_safe_memcpy(Hi_eval, v_temp);
   }
   gsl_vector_add_constant(v_temp, 1.0);
   gsl_vector_div(Hi_eval, v_temp);
@@ -675,10 +675,10 @@ double LogRL_f(double l, void *params) {
   double c = 0.5 * df * (log(df) - log(2 * M_PI) - 1.0);
   f = c - 0.5 * logdet_h - 0.5 * logdet_hiw - 0.5 * df * log(P_yy);
 
-  gsl_matrix_free(Pab);
-  gsl_matrix_free(Iab);
-  gsl_vector_free(Hi_eval);
-  gsl_vector_free(v_temp);
+  gsl_matrix_safe_free(Pab);
+  gsl_matrix_safe_free(Iab);
+  gsl_vector_safe_free(Hi_eval);
+  gsl_vector_safe_free(v_temp);
   return f;
 }
 
@@ -707,17 +707,17 @@ double LogRL_dev1(double l, void *params) {
   gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
   gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
 
-  gsl_vector_memcpy(v_temp, p->eval);
+  gsl_vector_safe_memcpy(v_temp, p->eval);
   gsl_vector_scale(v_temp, l);
   if (p->e_mode == 0) {
     gsl_vector_set_all(Hi_eval, 1.0);
   } else {
-    gsl_vector_memcpy(Hi_eval, v_temp);
+    gsl_vector_safe_memcpy(Hi_eval, v_temp);
   }
   gsl_vector_add_constant(v_temp, 1.0);
   gsl_vector_div(Hi_eval, v_temp);
 
-  gsl_vector_memcpy(HiHi_eval, Hi_eval);
+  gsl_vector_safe_memcpy(HiHi_eval, Hi_eval);
   gsl_vector_mul(HiHi_eval, Hi_eval);
 
   gsl_vector_set_all(v_temp, 1.0);
@@ -749,11 +749,11 @@ double LogRL_dev1(double l, void *params) {
 
   dev1 = -0.5 * trace_PK + 0.5 * df * yPKPy / P_yy;
 
-  gsl_matrix_free(Pab);
-  gsl_matrix_free(PPab);
-  gsl_vector_free(Hi_eval);
-  gsl_vector_free(HiHi_eval);
-  gsl_vector_free(v_temp);
+  gsl_matrix_safe_free(Pab);  // FIXME: may contain NaN
+  gsl_matrix_safe_free(PPab); // FIXME: may contain NaN
+  gsl_vector_safe_free(Hi_eval);
+  gsl_vector_safe_free(HiHi_eval);
+  gsl_vector_safe_free(v_temp);
 
   return dev1;
 }
@@ -785,19 +785,19 @@ double LogRL_dev2(double l, void *params) {
   gsl_vector *HiHiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
   gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
 
-  gsl_vector_memcpy(v_temp, p->eval);
+  gsl_vector_safe_memcpy(v_temp, p->eval);
   gsl_vector_scale(v_temp, l);
   if (p->e_mode == 0) {
     gsl_vector_set_all(Hi_eval, 1.0);
   } else {
-    gsl_vector_memcpy(Hi_eval, v_temp);
+    gsl_vector_safe_memcpy(Hi_eval, v_temp);
   }
   gsl_vector_add_constant(v_temp, 1.0);
   gsl_vector_div(Hi_eval, v_temp);
 
-  gsl_vector_memcpy(HiHi_eval, Hi_eval);
+  gsl_vector_safe_memcpy(HiHi_eval, Hi_eval);
   gsl_vector_mul(HiHi_eval, Hi_eval);
-  gsl_vector_memcpy(HiHiHi_eval, HiHi_eval);
+  gsl_vector_safe_memcpy(HiHiHi_eval, HiHi_eval);
   gsl_vector_mul(HiHiHi_eval, Hi_eval);
 
   gsl_vector_set_all(v_temp, 1.0);
@@ -837,13 +837,13 @@ double LogRL_dev2(double l, void *params) {
   dev2 = 0.5 * trace_PKPK -
          0.5 * df * (2.0 * yPKPKPy * P_yy - yPKPy * yPKPy) / (P_yy * P_yy);
 
-  gsl_matrix_free(Pab);
-  gsl_matrix_free(PPab);
-  gsl_matrix_free(PPPab);
-  gsl_vector_free(Hi_eval);
-  gsl_vector_free(HiHi_eval);
-  gsl_vector_free(HiHiHi_eval);
-  gsl_vector_free(v_temp);
+  gsl_matrix_safe_free(Pab);  // FIXME
+  gsl_matrix_safe_free(PPab);
+  gsl_matrix_safe_free(PPPab);
+  gsl_vector_safe_free(Hi_eval);
+  gsl_vector_safe_free(HiHi_eval);
+  gsl_vector_safe_free(HiHiHi_eval);
+  gsl_vector_safe_free(v_temp);
 
   return dev2;
 }
@@ -875,19 +875,19 @@ void LogRL_dev12(double l, void *params, double *dev1, double *dev2) {
   gsl_vector *HiHiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
   gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
 
-  gsl_vector_memcpy(v_temp, p->eval);
+  gsl_vector_safe_memcpy(v_temp, p->eval);
   gsl_vector_scale(v_temp, l);
   if (p->e_mode == 0) {
     gsl_vector_set_all(Hi_eval, 1.0);
   } else {
-    gsl_vector_memcpy(Hi_eval, v_temp);
+    gsl_vector_safe_memcpy(Hi_eval, v_temp);
   }
   gsl_vector_add_constant(v_temp, 1.0);
   gsl_vector_div(Hi_eval, v_temp);
 
-  gsl_vector_memcpy(HiHi_eval, Hi_eval);
+  gsl_vector_safe_memcpy(HiHi_eval, Hi_eval);
   gsl_vector_mul(HiHi_eval, Hi_eval);
-  gsl_vector_memcpy(HiHiHi_eval, HiHi_eval);
+  gsl_vector_safe_memcpy(HiHiHi_eval, HiHi_eval);
   gsl_vector_mul(HiHiHi_eval, Hi_eval);
 
   gsl_vector_set_all(v_temp, 1.0);
@@ -929,13 +929,13 @@ void LogRL_dev12(double l, void *params, double *dev1, double *dev2) {
   *dev2 = 0.5 * trace_PKPK -
           0.5 * df * (2.0 * yPKPKPy * P_yy - yPKPy * yPKPy) / (P_yy * P_yy);
 
-  gsl_matrix_free(Pab);
-  gsl_matrix_free(PPab);
-  gsl_matrix_free(PPPab);
-  gsl_vector_free(Hi_eval);
-  gsl_vector_free(HiHi_eval);
-  gsl_vector_free(HiHiHi_eval);
-  gsl_vector_free(v_temp);
+  gsl_matrix_safe_free(Pab);  // FIXME
+  gsl_matrix_safe_free(PPab);
+  gsl_matrix_safe_free(PPPab);
+  gsl_vector_safe_free(Hi_eval);
+  gsl_vector_safe_free(HiHi_eval);
+  gsl_vector_safe_free(HiHiHi_eval);
+  gsl_vector_safe_free(v_temp);
 
   return;
 }
@@ -951,12 +951,12 @@ void LMM::CalcRLWald(const double &l, const FUNC_PARAM &params, double &beta,
   gsl_vector *Hi_eval = gsl_vector_safe_alloc(params.eval->size);
   gsl_vector *v_temp = gsl_vector_safe_alloc(params.eval->size);
 
-  gsl_vector_memcpy(v_temp, params.eval);
+  gsl_vector_safe_memcpy(v_temp, params.eval);
   gsl_vector_scale(v_temp, l);
   if (params.e_mode == 0) {
     gsl_vector_set_all(Hi_eval, 1.0);
   } else {
-    gsl_vector_memcpy(Hi_eval, v_temp);
+    gsl_vector_safe_memcpy(Hi_eval, v_temp);
   }
   gsl_vector_add_constant(v_temp, 1.0);
   gsl_vector_div(Hi_eval, v_temp);
@@ -976,9 +976,9 @@ void LMM::CalcRLWald(const double &l, const FUNC_PARAM &params, double &beta,
   se = sqrt(1.0 / (tau * P_xx));
   p_wald = gsl_cdf_fdist_Q((P_yy - Px_yy) * tau, 1.0, df);
 
-  gsl_matrix_free(Pab);
-  gsl_vector_free(Hi_eval);
-  gsl_vector_free(v_temp);
+  gsl_matrix_safe_free(Pab);
+  gsl_vector_safe_free(Hi_eval);
+  gsl_vector_safe_free(v_temp);
   return;
 }
 
@@ -993,12 +993,12 @@ void LMM::CalcRLScore(const double &l, const FUNC_PARAM &params, double &beta,
   gsl_vector *Hi_eval = gsl_vector_safe_alloc(params.eval->size);
   gsl_vector *v_temp = gsl_vector_safe_alloc(params.eval->size);
 
-  gsl_vector_memcpy(v_temp, params.eval);
+  gsl_vector_safe_memcpy(v_temp, params.eval);
   gsl_vector_scale(v_temp, l);
   if (params.e_mode == 0) {
     gsl_vector_set_all(Hi_eval, 1.0);
   } else {
-    gsl_vector_memcpy(Hi_eval, v_temp);
+    gsl_vector_safe_memcpy(Hi_eval, v_temp);
   }
   gsl_vector_add_constant(v_temp, 1.0);
   gsl_vector_div(Hi_eval, v_temp);
@@ -1020,9 +1020,9 @@ void LMM::CalcRLScore(const double &l, const FUNC_PARAM &params, double &beta,
   p_score =
       gsl_cdf_fdist_Q((double)ni_test * P_xy * P_xy / (P_yy * P_xx), 1.0, df);
 
-  gsl_matrix_free(Pab);
-  gsl_vector_free(Hi_eval);
-  gsl_vector_free(v_temp);
+  gsl_matrix_safe_free(Pab);
+  gsl_vector_safe_free(Hi_eval);
+  gsl_vector_safe_free(v_temp);
   return;
 }
 
@@ -1038,10 +1038,10 @@ void CalcUab(const gsl_matrix *UtW, const gsl_vector *Uty, gsl_matrix *Uab) {
     }
 
     if (a == n_cvt + 2) {
-      gsl_vector_memcpy(u_a, Uty);
+      gsl_vector_safe_memcpy(u_a, Uty);
     } else {
       gsl_vector_const_view UtW_col = gsl_matrix_const_column(UtW, a - 1);
-      gsl_vector_memcpy(u_a, &UtW_col.vector);
+      gsl_vector_safe_memcpy(u_a, &UtW_col.vector);
     }
 
     for (size_t b = a; b >= 1; --b) {
@@ -1053,17 +1053,17 @@ void CalcUab(const gsl_matrix *UtW, const gsl_vector *Uty, gsl_matrix *Uab) {
       gsl_vector_view Uab_col = gsl_matrix_column(Uab, index_ab);
 
       if (b == n_cvt + 2) {
-        gsl_vector_memcpy(&Uab_col.vector, Uty);
+        gsl_vector_safe_memcpy(&Uab_col.vector, Uty);
       } else {
         gsl_vector_const_view UtW_col = gsl_matrix_const_column(UtW, b - 1);
-        gsl_vector_memcpy(&Uab_col.vector, &UtW_col.vector);
+        gsl_vector_safe_memcpy(&Uab_col.vector, &UtW_col.vector);
       }
 
       gsl_vector_mul(&Uab_col.vector, u_a);
     }
   }
 
-  gsl_vector_free(u_a);
+  gsl_vector_safe_free(u_a);
   return;
 }
 
@@ -1077,12 +1077,12 @@ void CalcUab(const gsl_matrix *UtW, const gsl_vector *Uty,
     gsl_vector_view Uab_col = gsl_matrix_column(Uab, index_ab);
 
     if (b == n_cvt + 2) {
-      gsl_vector_memcpy(&Uab_col.vector, Uty);
+      gsl_vector_safe_memcpy(&Uab_col.vector, Uty);
     } else if (b == n_cvt + 1) {
-      gsl_vector_memcpy(&Uab_col.vector, Utx);
+      gsl_vector_safe_memcpy(&Uab_col.vector, Utx);
     } else {
       gsl_vector_const_view UtW_col = gsl_matrix_const_column(UtW, b - 1);
-      gsl_vector_memcpy(&Uab_col.vector, &UtW_col.vector);
+      gsl_vector_safe_memcpy(&Uab_col.vector, &UtW_col.vector);
     }
 
     gsl_vector_mul(&Uab_col.vector, Utx);
@@ -1105,10 +1105,10 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
     }
 
     if (a == n_cvt + 2) {
-      gsl_vector_memcpy(v_a, y);
+      gsl_vector_safe_memcpy(v_a, y);
     } else {
       gsl_vector_const_view W_col = gsl_matrix_const_column(W, a - 1);
-      gsl_vector_memcpy(v_a, &W_col.vector);
+      gsl_vector_safe_memcpy(v_a, &W_col.vector);
     }
 
     for (size_t b = a; b >= 1; --b) {
@@ -1119,10 +1119,10 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
       index_ab = GetabIndex(a, b, n_cvt);
 
       if (b == n_cvt + 2) {
-        gsl_vector_memcpy(v_b, y);
+        gsl_vector_safe_memcpy(v_b, y);
       } else {
         gsl_vector_const_view W_col = gsl_matrix_const_column(W, b - 1);
-        gsl_vector_memcpy(v_b, &W_col.vector);
+        gsl_vector_safe_memcpy(v_b, &W_col.vector);
       }
 
       gsl_blas_ddot(v_a, v_b, &d);
@@ -1130,8 +1130,8 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
     }
   }
 
-  gsl_vector_free(v_a);
-  gsl_vector_free(v_b);
+  gsl_vector_safe_free(v_a);
+  gsl_vector_safe_free(v_b);
   return;
 }
 
@@ -1147,19 +1147,19 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, const gsl_vector *x,
     index_ab = GetabIndex(n_cvt + 1, b, n_cvt);
 
     if (b == n_cvt + 2) {
-      gsl_vector_memcpy(v_b, y);
+      gsl_vector_safe_memcpy(v_b, y);
     } else if (b == n_cvt + 1) {
-      gsl_vector_memcpy(v_b, x);
+      gsl_vector_safe_memcpy(v_b, x);
     } else {
       gsl_vector_const_view W_col = gsl_matrix_const_column(W, b - 1);
-      gsl_vector_memcpy(v_b, &W_col.vector);
+      gsl_vector_safe_memcpy(v_b, &W_col.vector);
     }
 
     gsl_blas_ddot(x, v_b, &d);
     gsl_vector_set(ab, index_ab, d);
   }
 
-  gsl_vector_free(v_b);
+  gsl_vector_safe_free(v_b);
   return;
 }
 
@@ -1260,10 +1260,10 @@ void LMM::AnalyzeGene(const gsl_matrix *U, const gsl_vector *eval,
   }
   cout << endl;
 
-  gsl_vector_free(y);
-  gsl_vector_free(Uty);
-  gsl_matrix_free(Uab);
-  gsl_vector_free(ab);
+  gsl_vector_safe_free(y);
+  gsl_vector_safe_free(Uty);
+  gsl_matrix_safe_free(Uab);
+  gsl_vector_safe_free(ab);
 
   infile.close();
   infile.clear();
@@ -1324,7 +1324,7 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
     for (size_t i = 0; i < l; i++) {
       // for every batch...
       gsl_vector_view UtXlarge_col = gsl_matrix_column(UtXlarge, i);
-      gsl_vector_memcpy(Utx, &UtXlarge_col.vector);
+      gsl_vector_safe_memcpy(Utx, &UtXlarge_col.vector);
 
       CalcUab(UtW, Uty, Utx, Uab);
 
@@ -1425,7 +1425,7 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
 
     // copy genotype values for SNP into Xlarge cache
     gsl_vector_view Xlarge_col = gsl_matrix_column(Xlarge, c % msize);
-    gsl_vector_memcpy(&Xlarge_col.vector, x);
+    gsl_vector_safe_memcpy(&Xlarge_col.vector, x);
     c++; // count SNPs going in
 
     if (c % msize == 0)
@@ -1436,14 +1436,14 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
   // cout << "Counted SNPs " << c << " sumStat " << sumStat.size() << endl;
   cout << endl;
 
-  gsl_vector_free(x);
-  gsl_vector_free(x_miss);
-  gsl_vector_free(Utx);
-  gsl_matrix_free(Uab);
-  gsl_vector_free(ab);
+  gsl_vector_safe_free(x);
+  gsl_vector_safe_free(x_miss);
+  gsl_vector_safe_free(Utx);
+  gsl_matrix_safe_free(Uab);
+  gsl_vector_safe_free(ab);
 
-  gsl_matrix_free(Xlarge);
-  gsl_matrix_free(UtXlarge);
+  gsl_matrix_safe_free(Xlarge);
+  gsl_matrix_safe_free(UtXlarge);
 
 }
 
@@ -1598,10 +1598,10 @@ void MatrixCalcLR(const gsl_matrix *U, const gsl_matrix *UtX,
     pos_loglr.push_back(make_pair(i, log_lr));
   }
 
-  gsl_vector_free(w);
-  gsl_matrix_free(Utw);
-  gsl_matrix_free(Uab);
-  gsl_vector_free(ab);
+  gsl_vector_safe_free(w);
+  gsl_matrix_safe_free(Utw);
+  gsl_matrix_safe_free(Uab);
+  gsl_vector_safe_free(ab);
 
   return;
 }
@@ -1790,8 +1790,8 @@ void CalcLambda(const char func_name, const gsl_vector *eval,
 
   CalcLambda(func_name, param0, l_min, l_max, n_region, lambda, logl_H0);
 
-  gsl_matrix_free(Uab);
-  gsl_vector_free(ab);
+  gsl_matrix_safe_free(Uab);
+  gsl_vector_safe_free(ab);
 
   return;
 }
@@ -1816,8 +1816,8 @@ void CalcPve(const gsl_vector *eval, const gsl_matrix *UtW,
   pve = trace_G * lambda / (trace_G * lambda + 1.0);
   pve_se = trace_G / ((trace_G * lambda + 1.0) * (trace_G * lambda + 1.0)) * se;
 
-  gsl_matrix_free(Uab);
-  gsl_vector_free(ab);
+  gsl_matrix_safe_free(Uab);
+  gsl_vector_safe_free(ab);
   return;
 }
 
@@ -1843,14 +1843,14 @@ void CalcLmmVgVeBeta(const gsl_vector *eval, const gsl_matrix *UtW,
   gsl_matrix_set_zero(Uab);
   CalcUab(UtW, Uty, Uab);
 
-  gsl_vector_memcpy(v_temp, eval);
+  gsl_vector_safe_memcpy(v_temp, eval);
   gsl_vector_scale(v_temp, lambda);
   gsl_vector_set_all(Hi_eval, 1.0);
   gsl_vector_add_constant(v_temp, 1.0);
   gsl_vector_div(Hi_eval, v_temp);
 
   // Calculate beta.
-  gsl_matrix_memcpy(HiW, UtW);
+  gsl_matrix_safe_memcpy(HiW, UtW);
   for (size_t i = 0; i < UtW->size2; i++) {
     gsl_vector_view HiW_col = gsl_matrix_column(HiW, i);
     gsl_vector_mul(&HiW_col.vector, Hi_eval);
@@ -1881,15 +1881,15 @@ void CalcLmmVgVeBeta(const gsl_vector *eval, const gsl_matrix *UtW,
     gsl_vector_set(se_beta, i, sqrt(gsl_matrix_get(Vbeta, i, i)));
   }
 
-  gsl_matrix_free(Uab);
-  gsl_matrix_free(Pab);
-  gsl_vector_free(ab);
-  gsl_vector_free(Hi_eval);
-  gsl_vector_free(v_temp);
-  gsl_matrix_free(HiW);
-  gsl_matrix_free(WHiW);
-  gsl_vector_free(WHiy);
-  gsl_matrix_free(Vbeta);
+  gsl_matrix_safe_free(Uab);
+  gsl_matrix_safe_free(Pab);
+  gsl_vector_safe_free(ab);
+  gsl_vector_safe_free(Hi_eval);
+  gsl_vector_safe_free(v_temp);
+  gsl_matrix_safe_free(HiW);
+  gsl_matrix_safe_free(WHiW);
+  gsl_vector_safe_free(WHiy);
+  gsl_matrix_safe_free(Vbeta);
 
   gsl_permutation_free(pmt);
   return;
@@ -1929,7 +1929,7 @@ void LMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
   gsl_matrix *UtW_expand = gsl_matrix_safe_alloc(U->size1, UtW->size2 + 2);
   gsl_matrix_view UtW_expand_mat =
       gsl_matrix_submatrix(UtW_expand, 0, 0, U->size1, UtW->size2);
-  gsl_matrix_memcpy(&UtW_expand_mat.matrix, UtW);
+  gsl_matrix_safe_memcpy(&UtW_expand_mat.matrix, UtW);
   gsl_vector_view UtW_expand_env = gsl_matrix_column(UtW_expand, UtW->size2);
   gsl_blas_dgemv(CblasTrans, 1.0, U, env, 0.0, &UtW_expand_env.vector);
   gsl_vector_view UtW_expand_x = gsl_matrix_column(UtW_expand, UtW->size2 + 1);
@@ -2030,13 +2030,13 @@ void LMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
   }
   cout << endl;
 
-  gsl_vector_free(x);
-  gsl_vector_free(x_miss);
-  gsl_vector_free(Utx);
-  gsl_matrix_free(Uab);
-  gsl_vector_free(ab);
+  gsl_vector_safe_free(x);
+  gsl_vector_safe_free(x_miss);
+  gsl_vector_safe_free(Utx);
+  gsl_matrix_safe_free(Uab);
+  gsl_vector_safe_free(ab);
 
-  gsl_matrix_free(UtW_expand);
+  gsl_matrix_safe_free(UtW_expand);
 
   infile.close();
   infile.clear();
@@ -2078,7 +2078,7 @@ void LMM::AnalyzePlinkGXE(const gsl_matrix *U, const gsl_vector *eval,
   gsl_matrix *UtW_expand = gsl_matrix_safe_alloc(U->size1, UtW->size2 + 2);
   gsl_matrix_view UtW_expand_mat =
       gsl_matrix_submatrix(UtW_expand, 0, 0, U->size1, UtW->size2);
-  gsl_matrix_memcpy(&UtW_expand_mat.matrix, UtW);
+  gsl_matrix_safe_memcpy(&UtW_expand_mat.matrix, UtW);
   gsl_vector_view UtW_expand_env = gsl_matrix_column(UtW_expand, UtW->size2);
   gsl_blas_dgemv(CblasTrans, 1.0, U, env, 0.0, &UtW_expand_env.vector);
   gsl_vector_view UtW_expand_x = gsl_matrix_column(UtW_expand, UtW->size2 + 1);
@@ -2208,12 +2208,12 @@ void LMM::AnalyzePlinkGXE(const gsl_matrix *U, const gsl_vector *eval,
   }
   cout << endl;
 
-  gsl_vector_free(x);
-  gsl_vector_free(Utx);
-  gsl_matrix_free(Uab);
-  gsl_vector_free(ab);
+  gsl_vector_safe_free(x);
+  gsl_vector_safe_free(Utx);
+  gsl_matrix_safe_free(Uab);
+  gsl_vector_safe_free(ab);
 
-  gsl_matrix_free(UtW_expand);
+  gsl_matrix_safe_free(UtW_expand);
 
   infile.close();
   infile.clear();
diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index ba71b64..908a52d 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -64,6 +64,33 @@ bool has_nan(const vector<double> v) {
   return false;
 }
 
+bool has_nan(const gsl_vector *v) {
+  for (size_t i = 0; i < v->size; ++i)
+    if (gsl_isnan(gsl_vector_get(v,i))) return true;
+  return false;
+}
+bool has_inf(const gsl_vector *v) {
+  for (size_t i = 0; i < v->size; ++i) {
+    auto value = gsl_vector_get(v,i);
+    if (gsl_isinf(value) != 0) return true;
+  }
+  return false;
+}
+bool has_nan(const gsl_matrix *m) {
+  for (size_t i = 0; i < m->size1; ++i)
+    for (size_t j = 0; j < m->size2; ++j)
+      if (gsl_isnan(gsl_matrix_get(m,i,j))) return true;
+  return false;
+}
+bool has_inf(const gsl_matrix *m) {
+  for (size_t i = 0; i < m->size1; ++i)
+    for (size_t j = 0; j < m->size2; ++j) {
+      auto value = gsl_matrix_get(m,i,j);
+      if (gsl_isinf(value) != 0) return true;
+    }
+  return false;
+}
+
 // calculate variance of a vector
 double VectorVar(const gsl_vector *v) {
   double d, m = 0.0, m2 = 0.0;
@@ -96,8 +123,8 @@ void CenterMatrix(gsl_matrix *G) {
     }
   }
 
-  gsl_vector_free(w);
-  gsl_vector_free(Gw);
+  gsl_vector_safe_free(w);
+  gsl_vector_safe_free(Gw);
 
   return;
 }
@@ -120,7 +147,7 @@ void CenterMatrix(gsl_matrix *G, const gsl_vector *w) {
     }
   }
 
-  gsl_vector_free(Gw);
+  gsl_vector_safe_free(Gw);
 
   return;
 }
@@ -156,12 +183,12 @@ void CenterMatrix(gsl_matrix *G, const gsl_matrix *W) {
 
   gsl_matrix_add(G, Gtmp);
 
-  gsl_matrix_free(WtW);
-  gsl_matrix_free(WtWi);
-  gsl_matrix_free(WtWiWt);
-  gsl_matrix_free(GW);
-  gsl_matrix_free(WtGW);
-  gsl_matrix_free(Gtmp);
+  gsl_matrix_safe_free(WtW);
+  gsl_matrix_safe_free(WtWi);
+  gsl_matrix_safe_free(WtWiWt);
+  gsl_matrix_safe_free(GW);
+  gsl_matrix_safe_free(WtGW);
+  gsl_matrix_safe_free(Gtmp);
 
   return;
 }
@@ -228,7 +255,7 @@ bool isMatrixSymmetric(const gsl_matrix *G) {
 bool isMatrixPositiveDefinite(const gsl_matrix *G) {
   enforce(G->size1 == G->size2);
   auto G2 = gsl_matrix_safe_alloc(G->size1, G->size2);
-  enforce_gsl(gsl_matrix_memcpy(G2,G));
+  enforce_gsl(gsl_matrix_safe_memcpy(G2,G));
   auto handler = gsl_set_error_handler_off();
 #if GSL_MAJOR_VERSION >= 2 && GSL_MINOR_VERSION >= 3
   auto s = gsl_linalg_cholesky_decomp1(G2);
@@ -236,20 +263,20 @@ bool isMatrixPositiveDefinite(const gsl_matrix *G) {
   auto s = gsl_linalg_cholesky_decomp(G2);
 #endif
   gsl_set_error_handler(handler);
-  gsl_matrix_free(G2);
+  gsl_matrix_safe_free(G2);
   return (s == GSL_SUCCESS);
 }
 
 gsl_vector *getEigenValues(const gsl_matrix *G) {
   enforce(G->size1 == G->size2);
   auto G2 = gsl_matrix_safe_alloc(G->size1, G->size2);
-  enforce_gsl(gsl_matrix_memcpy(G2,G));
+  enforce_gsl(gsl_matrix_safe_memcpy(G2,G));
   auto eworkspace = gsl_eigen_symm_alloc(G->size1);
   enforce(eworkspace);
   gsl_vector *eigenvalues = gsl_vector_safe_alloc(G->size1);
   enforce_gsl(gsl_eigen_symm(G2, eigenvalues, eworkspace));
   gsl_eigen_symm_free(eworkspace);
-  gsl_matrix_free(G2);
+  gsl_matrix_safe_free(G2);
   return eigenvalues;
 }
 
@@ -359,9 +386,9 @@ void CenterVector(gsl_vector *y, const gsl_matrix *W) {
 
   gsl_blas_dgemv(CblasNoTrans, -1.0, W, WtWiWty, 1.0, y);
 
-  gsl_matrix_free(WtW);
-  gsl_vector_free(Wty);
-  gsl_vector_free(WtWiWty);
+  gsl_matrix_safe_free(WtW);
+  gsl_vector_safe_free(Wty);
+  gsl_vector_safe_free(WtWiWty);
 
   return;
 }
@@ -388,9 +415,9 @@ void StandardizeVector(gsl_vector *y) {
 // Calculate UtX.
 void CalcUtX(const gsl_matrix *U, gsl_matrix *UtX) {
   gsl_matrix *X = gsl_matrix_safe_alloc(UtX->size1, UtX->size2);
-  gsl_matrix_memcpy(X, UtX);
+  gsl_matrix_safe_memcpy(X, UtX);
   fast_dgemm("T", "N", 1.0, U, X, 0.0, UtX);
-  gsl_matrix_free(X);
+  gsl_matrix_safe_free(X);
 }
 
 void CalcUtX(const gsl_matrix *U, const gsl_matrix *X, gsl_matrix *UtX) {
@@ -407,7 +434,7 @@ void Kronecker(const gsl_matrix *K, const gsl_matrix *V, gsl_matrix *H) {
     for (size_t j = 0; j < K->size2; j++) {
       gsl_matrix_view H_sub = gsl_matrix_submatrix(
           H, i * V->size1, j * V->size2, V->size1, V->size2);
-      gsl_matrix_memcpy(&H_sub.matrix, V);
+      gsl_matrix_safe_memcpy(&H_sub.matrix, V);
       gsl_matrix_scale(&H_sub.matrix, gsl_matrix_get(K, i, j));
     }
   }
@@ -420,13 +447,13 @@ void KroneckerSym(const gsl_matrix *K, const gsl_matrix *V, gsl_matrix *H) {
     for (size_t j = i; j < K->size2; j++) {
       gsl_matrix_view H_sub = gsl_matrix_submatrix(
           H, i * V->size1, j * V->size2, V->size1, V->size2);
-      gsl_matrix_memcpy(&H_sub.matrix, V);
+      gsl_matrix_safe_memcpy(&H_sub.matrix, V);
       gsl_matrix_scale(&H_sub.matrix, gsl_matrix_get(K, i, j));
 
       if (i != j) {
         gsl_matrix_view H_sub_sym = gsl_matrix_submatrix(
             H, j * V->size1, i * V->size2, V->size1, V->size2);
-        gsl_matrix_memcpy(&H_sub_sym.matrix, &H_sub.matrix);
+        gsl_matrix_safe_memcpy(&H_sub_sym.matrix, &H_sub.matrix);
       }
     }
   }
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 1319a64..42e76e5 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -30,6 +30,10 @@ using namespace std;
 
 
 bool has_nan(const vector<double> v);
+bool has_nan(const gsl_vector *v);
+bool has_inf(const gsl_vector *v);
+bool has_nan(const gsl_matrix *m);
+bool has_inf(const gsl_matrix *m);
 
 double VectorVar(const gsl_vector *v);
 void CenterMatrix(gsl_matrix *G);