Tests still pass with safe_alloc (which sets the buffers to NaNs on allocation)

8 files changed, 286 insertions, 269 deletions

diff --git a/src/debug.cpp b/src/debug.cpp
index eb5c041..4db3daf 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -19,7 +19,7 @@
 #include "mathfunc.h"
 
 static bool debug_mode     = false;
-static bool debug_check    = false; // check data/algorithms
+static bool debug_check    = true;  // check data/algorithms
 static bool debug_strict   = false; // fail on error
 static bool debug_quiet    = false;
 static uint debug_issue    = 0;     // github issues
@@ -56,6 +56,22 @@ gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols) {
   return m;
 }
 
+/*
+  Helper function to make sure gsl allocations do their job because
+  gsl_vector_alloc does not initiatize values (behaviour that changed
+  in GSL2) we introduced a 'strict mode' by initializing the buffer
+  with NaNs. This happens when NO-CHECKS is not set
+  (i.e. -no-checks option).
+*/
+gsl_vector *gsl_vector_safe_alloc(size_t n) {
+  gsl_vector *v = gsl_vector_alloc(n);
+  enforce_msg(v,"Not enough memory"); // just to be sure when there is no error handler set
+  if (is_check_mode()) {
+    gsl_vector_set_all(v, nan(""));
+  }
+  return v;
+}
+
 // Helper function called by macro validate_K(K, check)
 void do_validate_K(const gsl_matrix *K, const char *__file, int __line) {
   if (is_check_mode()) {
diff --git a/src/debug.h b/src/debug.h
index 7d2197d..c127630 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -26,6 +26,7 @@ bool is_issue(uint issue);
 bool is_legacy_mode();
 
 gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols);
+gsl_vector *gsl_vector_safe_alloc(size_t n);
 
 // Validation routines
 void do_validate_K(const gsl_matrix *K, const char *__file, int __line);
diff --git a/src/fastblas.cpp b/src/fastblas.cpp
index 35d19c1..de7e27b 100644
--- a/src/fastblas.cpp
+++ b/src/fastblas.cpp
@@ -197,7 +197,7 @@ void fast_dgemm(const char *TransA, const char *TransB, const double alpha,
   } else {
     fast_cblas_dgemm(TransA,TransB,alpha,A,B,beta,C);
 
-    #ifndef NDEBUG
+    #ifdef DISABLE
     if (is_check_mode()) {
       // ---- validate with original implementation
       gsl_matrix *C1 = gsl_matrix_alloc(C->size1,C->size2);
diff --git a/src/gemma.cpp b/src/gemma.cpp
index 75b1c33..f9a2fc9 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -1636,7 +1636,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
   if (cPar.a_mode == 41 || cPar.a_mode == 42) {
     gsl_vector *y_prdt;
 
-    y_prdt = gsl_vector_alloc(cPar.ni_total - cPar.ni_test);
+    y_prdt = gsl_vector_safe_alloc(cPar.ni_total - cPar.ni_test);
 
     // set to zero
     gsl_vector_set_zero(y_prdt);
@@ -1648,8 +1648,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
     if (!cPar.file_kin.empty() && !cPar.file_ebv.empty()) {
       cout << "Adding Breeding Values ... " << endl;
 
-      gsl_matrix *G = gsl_matrix_alloc(cPar.ni_total, cPar.ni_total);
-      gsl_vector *u_hat = gsl_vector_alloc(cPar.ni_test);
+      gsl_matrix *G = gsl_matrix_safe_alloc(cPar.ni_total, cPar.ni_total);
+      gsl_vector *u_hat = gsl_vector_safe_alloc(cPar.ni_test);
 
       // read kinship matrix and set u_hat
       vector<int> indicator_all;
@@ -1707,25 +1707,25 @@ void GEMMA::BatchRun(PARAM &cPar) {
   if (cPar.a_mode == 43) {
     // first, use individuals with full phenotypes to obtain estimates of Vg and
     // Ve
-    gsl_matrix *Y = gsl_matrix_alloc(cPar.ni_test, cPar.n_ph);
-    gsl_matrix *W = gsl_matrix_alloc(Y->size1, cPar.n_cvt);
-    gsl_matrix *G = gsl_matrix_alloc(Y->size1, Y->size1);
-    gsl_matrix *U = gsl_matrix_alloc(Y->size1, Y->size1);
-    gsl_matrix *UtW = gsl_matrix_alloc(Y->size1, W->size2);
-    gsl_matrix *UtY = gsl_matrix_alloc(Y->size1, Y->size2);
-    gsl_vector *eval = gsl_vector_alloc(Y->size1);
+    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);
+    gsl_matrix *G = gsl_matrix_safe_alloc(Y->size1, Y->size1);
+    gsl_matrix *U = gsl_matrix_safe_alloc(Y->size1, Y->size1);
+    gsl_matrix *UtW = gsl_matrix_safe_alloc(Y->size1, W->size2);
+    gsl_matrix *UtY = gsl_matrix_safe_alloc(Y->size1, Y->size2);
+    gsl_vector *eval = gsl_vector_safe_alloc(Y->size1);
 
-    gsl_matrix *Y_full = gsl_matrix_alloc(cPar.ni_cvt, cPar.n_ph);
-    gsl_matrix *W_full = gsl_matrix_alloc(Y_full->size1, cPar.n_cvt);
+    gsl_matrix *Y_full = gsl_matrix_safe_alloc(cPar.ni_cvt, cPar.n_ph);
+    gsl_matrix *W_full = gsl_matrix_safe_alloc(Y_full->size1, cPar.n_cvt);
 
     // set covariates matrix W and phenotype matrix Y
     // an intercept should be included in W,
     cPar.CopyCvtPhen(W, Y, 0);
     cPar.CopyCvtPhen(W_full, Y_full, 1);
 
-    gsl_matrix *Y_hat = gsl_matrix_alloc(Y_full->size1, cPar.n_ph);
-    gsl_matrix *G_full = gsl_matrix_alloc(Y_full->size1, Y_full->size1);
-    gsl_matrix *H_full = gsl_matrix_alloc(Y_full->size1 * Y_hat->size2,
+    gsl_matrix *Y_hat = gsl_matrix_safe_alloc(Y_full->size1, cPar.n_ph);
+    gsl_matrix *G_full = gsl_matrix_safe_alloc(Y_full->size1, Y_full->size1);
+    gsl_matrix *H_full = gsl_matrix_safe_alloc(Y_full->size1 * Y_hat->size2,
                                           Y_full->size1 * Y_hat->size2);
 
     // read relatedness matrix G, and matrix G_full
@@ -1761,8 +1761,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
     // calculate variance component and beta estimates
     // and then obtain predicted values
     if (cPar.n_ph == 1) {
-      gsl_vector *beta = gsl_vector_alloc(W->size2);
-      gsl_vector *se_beta = gsl_vector_alloc(W->size2);
+      gsl_vector *beta = gsl_vector_safe_alloc(W->size2);
+      gsl_vector *se_beta = gsl_vector_safe_alloc(W->size2);
 
       double lambda, logl, vg, ve;
       gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0);
@@ -1792,10 +1792,10 @@ void GEMMA::BatchRun(PARAM &cPar) {
       gsl_vector_free(beta);
       gsl_vector_free(se_beta);
     } else {
-      gsl_matrix *Vg = gsl_matrix_alloc(cPar.n_ph, cPar.n_ph);
-      gsl_matrix *Ve = gsl_matrix_alloc(cPar.n_ph, cPar.n_ph);
-      gsl_matrix *B = gsl_matrix_alloc(cPar.n_ph, W->size2);
-      gsl_matrix *se_B = gsl_matrix_alloc(cPar.n_ph, W->size2);
+      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);
+      gsl_matrix *B = gsl_matrix_safe_alloc(cPar.n_ph, W->size2);
+      gsl_matrix *se_B = gsl_matrix_safe_alloc(cPar.n_ph, W->size2);
 
       // obtain estimates
       CalcMvLmmVgVeBeta(eval, UtW, UtY, cPar.em_iter, cPar.nr_iter,
@@ -1873,7 +1873,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
   if (cPar.a_mode == 21 || cPar.a_mode == 22) {
     cout << "Calculating Relatedness Matrix ... " << endl;
 
-    gsl_matrix *G = gsl_matrix_alloc(cPar.ni_total, cPar.ni_total);
+    gsl_matrix *G = gsl_matrix_safe_alloc(cPar.ni_total, cPar.ni_total);
     enforce_msg(G, "allocate G"); // just to be sure
 
     time_start = clock();
@@ -1918,8 +1918,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
   if (cPar.a_mode == 25 || cPar.a_mode == 26) {
     cout << "Calculating the S Matrix ... " << endl;
 
-    gsl_matrix *S = gsl_matrix_alloc(cPar.n_vc * 2, cPar.n_vc);
-    gsl_vector *ns = gsl_vector_alloc(cPar.n_vc + 1);
+    gsl_matrix *S = gsl_matrix_safe_alloc(cPar.n_vc * 2, cPar.n_vc);
+    gsl_vector *ns = gsl_vector_safe_alloc(cPar.n_vc + 1);
     gsl_matrix_set_zero(S);
     gsl_vector_set_zero(ns);
 
@@ -1928,13 +1928,13 @@ void GEMMA::BatchRun(PARAM &cPar) {
         gsl_matrix_submatrix(S, cPar.n_vc, 0, cPar.n_vc, cPar.n_vc);
     gsl_vector_view ns_vec = gsl_vector_subvector(ns, 0, cPar.n_vc);
 
-    gsl_matrix *K = gsl_matrix_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test);
-    gsl_matrix *A = gsl_matrix_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test);
+    gsl_matrix *K = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test);
+    gsl_matrix *A = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test);
     gsl_matrix_set_zero(K);
     gsl_matrix_set_zero(A);
 
-    gsl_vector *y = gsl_vector_alloc(cPar.ni_test);
-    gsl_matrix *W = gsl_matrix_alloc(cPar.ni_test, cPar.n_cvt);
+    gsl_vector *y = gsl_vector_safe_alloc(cPar.ni_test);
+    gsl_matrix *W = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_cvt);
 
     cPar.CopyCvtPhen(W, y, 0);
 
@@ -1971,9 +1971,9 @@ void GEMMA::BatchRun(PARAM &cPar) {
   // Compute the q vector, that is used for variance component estimation using
   // summary statistics
   if (cPar.a_mode == 27 || cPar.a_mode == 28) {
-    gsl_matrix *Vq = gsl_matrix_alloc(cPar.n_vc, cPar.n_vc);
-    gsl_vector *q = gsl_vector_alloc(cPar.n_vc);
-    gsl_vector *s = gsl_vector_alloc(cPar.n_vc + 1);
+    gsl_matrix *Vq = gsl_matrix_safe_alloc(cPar.n_vc, cPar.n_vc);
+    gsl_vector *q = gsl_vector_safe_alloc(cPar.n_vc);
+    gsl_vector *s = gsl_vector_safe_alloc(cPar.n_vc + 1);
     gsl_vector_set_zero(q);
     gsl_vector_set_zero(s);
 
@@ -2029,8 +2029,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
   // LM.
   if (cPar.a_mode == 51 || cPar.a_mode == 52 || cPar.a_mode == 53 ||
       cPar.a_mode == 54) { // Fit LM
-    gsl_matrix *Y = gsl_matrix_alloc(cPar.ni_test, cPar.n_ph);
-    gsl_matrix *W = gsl_matrix_alloc(Y->size1, cPar.n_cvt);
+    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);
 
     // set covariates matrix W and phenotype matrix Y
     // an intercept should be included in W,
@@ -2082,16 +2082,16 @@ void GEMMA::BatchRun(PARAM &cPar) {
       cPar.UpdateSNP(mapRS2wK);
 
       // Setup matrices and vectors.
-      gsl_matrix *S = gsl_matrix_alloc(cPar.n_vc * 2, cPar.n_vc);
-      gsl_matrix *Vq = gsl_matrix_alloc(cPar.n_vc, cPar.n_vc);
-      gsl_vector *q = gsl_vector_alloc(cPar.n_vc);
-      gsl_vector *s = gsl_vector_alloc(cPar.n_vc + 1);
+      gsl_matrix *S = gsl_matrix_safe_alloc(cPar.n_vc * 2, cPar.n_vc);
+      gsl_matrix *Vq = gsl_matrix_safe_alloc(cPar.n_vc, cPar.n_vc);
+      gsl_vector *q = gsl_vector_safe_alloc(cPar.n_vc);
+      gsl_vector *s = gsl_vector_safe_alloc(cPar.n_vc + 1);
 
-      gsl_matrix *K = gsl_matrix_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test);
-      gsl_matrix *A = gsl_matrix_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test);
+      gsl_matrix *K = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test);
+      gsl_matrix *A = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test);
 
-      gsl_vector *y = gsl_vector_alloc(cPar.ni_test);
-      gsl_matrix *W = gsl_matrix_alloc(cPar.ni_test, cPar.n_cvt);
+      gsl_vector *y = gsl_vector_safe_alloc(cPar.ni_test);
+      gsl_matrix *W = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_cvt);
 
       gsl_matrix_set_zero(K);
       gsl_matrix_set_zero(A);
@@ -2218,16 +2218,16 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
       cPar.n_vc = cPar.n_vc - 1;
 
-      gsl_matrix *S = gsl_matrix_alloc(2 * cPar.n_vc, cPar.n_vc);
-      gsl_matrix *Vq = gsl_matrix_alloc(cPar.n_vc, cPar.n_vc);
-      // gsl_matrix *V=gsl_matrix_alloc (cPar.n_vc+1,
+      gsl_matrix *S = gsl_matrix_safe_alloc(2 * cPar.n_vc, cPar.n_vc);
+      gsl_matrix *Vq = gsl_matrix_safe_alloc(cPar.n_vc, cPar.n_vc);
+      // gsl_matrix *V=gsl_matrix_safe_alloc (cPar.n_vc+1,
       // (cPar.n_vc*(cPar.n_vc+1))/2*(cPar.n_vc+1) );
-      // gsl_matrix *Vslope=gsl_matrix_alloc (n_lines+1,
+      // gsl_matrix *Vslope=gsl_matrix_safe_alloc (n_lines+1,
       // (n_lines*(n_lines+1))/2*(n_lines+1) );
-      gsl_vector *q = gsl_vector_alloc(cPar.n_vc);
-      gsl_vector *s_study = gsl_vector_alloc(cPar.n_vc);
-      gsl_vector *s_ref = gsl_vector_alloc(cPar.n_vc);
-      gsl_vector *s = gsl_vector_alloc(cPar.n_vc + 1);
+      gsl_vector *q = gsl_vector_safe_alloc(cPar.n_vc);
+      gsl_vector *s_study = gsl_vector_safe_alloc(cPar.n_vc);
+      gsl_vector *s_ref = gsl_vector_safe_alloc(cPar.n_vc);
+      gsl_vector *s = gsl_vector_safe_alloc(cPar.n_vc + 1);
 
       gsl_matrix_set_zero(S);
       gsl_matrix_view S_mat =
@@ -2286,9 +2286,9 @@ void GEMMA::BatchRun(PARAM &cPar) {
       gsl_vector_free(s_ref);
       gsl_vector_free(s);
     } else {
-      gsl_matrix *Y = gsl_matrix_alloc(cPar.ni_test, cPar.n_ph);
-      gsl_matrix *W = gsl_matrix_alloc(Y->size1, cPar.n_cvt);
-      gsl_matrix *G = gsl_matrix_alloc(Y->size1, Y->size1 * cPar.n_vc);
+      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);
+      gsl_matrix *G = gsl_matrix_safe_alloc(Y->size1, Y->size1 * cPar.n_vc);
 
       // set covariates matrix W and phenotype matrix Y
       // an intercept should be included in W,
@@ -2365,9 +2365,9 @@ void GEMMA::BatchRun(PARAM &cPar) {
   // the genotypes
   if (cPar.a_mode == 66 || cPar.a_mode == 67) {
     // read reference file first
-    gsl_matrix *S = gsl_matrix_alloc(cPar.n_vc, cPar.n_vc);
-    gsl_matrix *Svar = gsl_matrix_alloc(cPar.n_vc, cPar.n_vc);
-    gsl_vector *s_ref = gsl_vector_alloc(cPar.n_vc);
+    gsl_matrix *S = gsl_matrix_safe_alloc(cPar.n_vc, cPar.n_vc);
+    gsl_matrix *Svar = gsl_matrix_safe_alloc(cPar.n_vc, cPar.n_vc);
+    gsl_vector *s_ref = gsl_vector_safe_alloc(cPar.n_vc);
 
     gsl_matrix_set_zero(S);
     gsl_matrix_set_zero(Svar);
@@ -2392,14 +2392,14 @@ void GEMMA::BatchRun(PARAM &cPar) {
     cPar.ObtainWeight(setSnps_beta, mapRS2wK);
 
     // set up matrices and vector
-    gsl_matrix *Xz = gsl_matrix_alloc(cPar.ni_test, cPar.n_vc);
-    gsl_matrix *XWz = gsl_matrix_alloc(cPar.ni_test, cPar.n_vc);
+    gsl_matrix *Xz = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_vc);
+    gsl_matrix *XWz = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_vc);
     gsl_matrix *XtXWz =
-        gsl_matrix_alloc(mapRS2wK.size(), cPar.n_vc * cPar.n_vc);
-    gsl_vector *w = gsl_vector_alloc(mapRS2wK.size());
-    gsl_vector *w1 = gsl_vector_alloc(mapRS2wK.size());
-    gsl_vector *z = gsl_vector_alloc(mapRS2wK.size());
-    gsl_vector *s_vec = gsl_vector_alloc(cPar.n_vc);
+        gsl_matrix_safe_alloc(mapRS2wK.size(), cPar.n_vc * cPar.n_vc);
+    gsl_vector *w = gsl_vector_safe_alloc(mapRS2wK.size());
+    gsl_vector *w1 = gsl_vector_safe_alloc(mapRS2wK.size());
+    gsl_vector *z = gsl_vector_safe_alloc(mapRS2wK.size());
+    gsl_vector *s_vec = gsl_vector_safe_alloc(cPar.n_vc);
 
     vector<size_t> vec_cat, vec_size;
     vector<double> vec_z;
@@ -2523,19 +2523,19 @@ void GEMMA::BatchRun(PARAM &cPar) {
   if (cPar.a_mode == 1 || cPar.a_mode == 2 || cPar.a_mode == 3 ||
       cPar.a_mode == 4 || cPar.a_mode == 5 ||
       cPar.a_mode == 31) { // Fit LMM or mvLMM or eigen
-    gsl_matrix *Y = gsl_matrix_alloc(cPar.ni_test, cPar.n_ph);
+    gsl_matrix *Y = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_ph);
     enforce_msg(Y, "allocate Y"); // just to be sure
-    gsl_matrix *W = gsl_matrix_alloc(Y->size1, cPar.n_cvt);
-    gsl_matrix *B = gsl_matrix_alloc(Y->size2, W->size2); // B is a d by c
+    gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt);
+    gsl_matrix *B = gsl_matrix_safe_alloc(Y->size2, W->size2); // B is a d by c
                                                           // matrix
-    gsl_matrix *se_B = gsl_matrix_alloc(Y->size2, W->size2);
-    gsl_matrix *G = gsl_matrix_alloc(Y->size1, Y->size1);
-    gsl_matrix *U = gsl_matrix_alloc(Y->size1, Y->size1);
+    gsl_matrix *se_B = gsl_matrix_safe_alloc(Y->size2, W->size2);
+    gsl_matrix *G = gsl_matrix_safe_alloc(Y->size1, Y->size1);
+    gsl_matrix *U = gsl_matrix_safe_alloc(Y->size1, Y->size1);
     gsl_matrix *UtW = gsl_matrix_calloc(Y->size1, W->size2);
     gsl_matrix *UtY = gsl_matrix_calloc(Y->size1, Y->size2);
     gsl_vector *eval = gsl_vector_calloc(Y->size1);
-    gsl_vector *env = gsl_vector_alloc(Y->size1);
-    gsl_vector *weight = gsl_vector_alloc(Y->size1);
+    gsl_vector *env = gsl_vector_safe_alloc(Y->size1);
+    gsl_vector *weight = gsl_vector_safe_alloc(Y->size1);
     assert_issue(is_issue(26), UtY->data[0] == 0.0);
 
     // set covariates matrix W and phenotype matrix Y
@@ -2667,26 +2667,22 @@ void GEMMA::BatchRun(PARAM &cPar) {
         CalcLambda('L', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max,
                    cPar.n_region, cPar.l_mle_null, cPar.logl_mle_H0);
         assert(!std::isnan(UtY->data[0]));
-        assert(!std::isnan(B->data[0]));
-        assert(!std::isnan(se_B->data[0]));
 
         CalcLmmVgVeBeta(eval, UtW, &UtY_col.vector, cPar.l_mle_null,
                         cPar.vg_mle_null, cPar.ve_mle_null, &beta.vector,
                         &se_beta.vector);
 
         assert(!std::isnan(UtY->data[0]));
-        assert(!std::isnan(B->data[0]));
-        assert(!std::isnan(se_B->data[0]));
 
         cPar.beta_mle_null.clear();
         cPar.se_beta_mle_null.clear();
+        assert(!std::isnan(B->data[0]));
+        assert(!std::isnan(se_B->data[0]));
         for (size_t i = 0; i < B->size2; i++) {
           cPar.beta_mle_null.push_back(gsl_matrix_get(B, 0, i));
           cPar.se_beta_mle_null.push_back(gsl_matrix_get(se_B, 0, i));
         }
         assert(!std::isnan(UtY->data[0]));
-        assert(!std::isnan(B->data[0]));
-        assert(!std::isnan(se_B->data[0]));
         assert(!std::isnan(cPar.beta_mle_null.front()));
         assert(!std::isnan(cPar.se_beta_mle_null.front()));
 
@@ -2698,6 +2694,9 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
         cPar.beta_remle_null.clear();
         cPar.se_beta_remle_null.clear();
+        assert(!std::isnan(B->data[0]));
+        assert(!std::isnan(se_B->data[0]));
+
         for (size_t i = 0; i < B->size2; i++) {
           cPar.beta_remle_null.push_back(gsl_matrix_get(B, 0, i));
           cPar.se_beta_remle_null.push_back(gsl_matrix_get(se_B, 0, i));
@@ -2709,11 +2708,11 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
         // calculate and output residuals
         if (cPar.a_mode == 5) {
-          gsl_vector *Utu_hat = gsl_vector_alloc(Y->size1);
-          gsl_vector *Ute_hat = gsl_vector_alloc(Y->size1);
-          gsl_vector *u_hat = gsl_vector_alloc(Y->size1);
-          gsl_vector *e_hat = gsl_vector_alloc(Y->size1);
-          gsl_vector *y_hat = gsl_vector_alloc(Y->size1);
+          gsl_vector *Utu_hat = gsl_vector_safe_alloc(Y->size1);
+          gsl_vector *Ute_hat = gsl_vector_safe_alloc(Y->size1);
+          gsl_vector *u_hat = gsl_vector_safe_alloc(Y->size1);
+          gsl_vector *e_hat = gsl_vector_safe_alloc(Y->size1);
+          gsl_vector *y_hat = gsl_vector_safe_alloc(Y->size1);
 
           // obtain Utu and Ute
           gsl_vector_memcpy(y_hat, &UtY_col.vector);
@@ -2816,10 +2815,10 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
   // BSLMM
   if (cPar.a_mode == 11 || cPar.a_mode == 12 || cPar.a_mode == 13) {
-    gsl_vector *y = gsl_vector_alloc(cPar.ni_test);
-    gsl_matrix *W = gsl_matrix_alloc(y->size, cPar.n_cvt);
-    gsl_matrix *G = gsl_matrix_alloc(y->size, y->size);
-    gsl_matrix *UtX = gsl_matrix_alloc(y->size, cPar.ns_test);
+    gsl_vector *y = gsl_vector_safe_alloc(cPar.ni_test);
+    gsl_matrix *W = gsl_matrix_safe_alloc(y->size, cPar.n_cvt);
+    gsl_matrix *G = gsl_matrix_safe_alloc(y->size, y->size);
+    gsl_matrix *UtX = gsl_matrix_safe_alloc(y->size, cPar.ns_test);
 
     // set covariates matrix W and phenotype vector y
     // an intercept should be included in W,
@@ -2842,10 +2841,10 @@ void GEMMA::BatchRun(PARAM &cPar) {
       cBslmm.CopyToParam(cPar);
       // else, if rho!=1
     } else {
-      gsl_matrix *U = gsl_matrix_alloc(y->size, y->size);
-      gsl_vector *eval = gsl_vector_alloc(y->size);
-      gsl_matrix *UtW = gsl_matrix_alloc(y->size, W->size2);
-      gsl_vector *Uty = gsl_vector_alloc(y->size);
+      gsl_matrix *U = gsl_matrix_safe_alloc(y->size, y->size);
+      gsl_vector *eval = gsl_vector_safe_alloc(y->size);
+      gsl_matrix *UtW = gsl_matrix_safe_alloc(y->size, W->size2);
+      gsl_vector *Uty = gsl_vector_safe_alloc(y->size);
 
       // read relatedness matrix G
       if (!(cPar.file_kin).empty()) {
@@ -2926,10 +2925,10 @@ void GEMMA::BatchRun(PARAM &cPar) {
   // BSLMM-DAP
   if (cPar.a_mode == 14 || cPar.a_mode == 15 || cPar.a_mode == 16) {
     if (cPar.a_mode == 14) {
-      gsl_vector *y = gsl_vector_alloc(cPar.ni_test);
-      gsl_matrix *W = gsl_matrix_alloc(y->size, cPar.n_cvt);
-      gsl_matrix *G = gsl_matrix_alloc(y->size, y->size);
-      gsl_matrix *UtX = gsl_matrix_alloc(y->size, cPar.ns_test);
+      gsl_vector *y = gsl_vector_safe_alloc(cPar.ni_test);
+      gsl_matrix *W = gsl_matrix_safe_alloc(y->size, cPar.n_cvt);
+      gsl_matrix *G = gsl_matrix_safe_alloc(y->size, y->size);
+      gsl_matrix *UtX = gsl_matrix_safe_alloc(y->size, cPar.ns_test);
 
       // set covariates matrix W and phenotype vector y
       // an intercept should be included in W,
@@ -2953,10 +2952,10 @@ void GEMMA::BatchRun(PARAM &cPar) {
         cBslmm.CopyToParam(cPar);
         // else, if rho!=1
       } else {
-        gsl_matrix *U = gsl_matrix_alloc(y->size, y->size);
-        gsl_vector *eval = gsl_vector_alloc(y->size);
-        gsl_matrix *UtW = gsl_matrix_alloc(y->size, W->size2);
-        gsl_vector *Uty = gsl_vector_alloc(y->size);
+        gsl_matrix *U = gsl_matrix_safe_alloc(y->size, y->size);
+        gsl_vector *eval = gsl_vector_safe_alloc(y->size);
+        gsl_matrix *UtW = gsl_matrix_safe_alloc(y->size, W->size2);
+        gsl_vector *Uty = gsl_vector_safe_alloc(y->size);
 
         // read relatedness matrix G
         if (!(cPar.file_kin).empty()) {
diff --git a/src/io.cpp b/src/io.cpp
index d1c4bc2..8abdeec 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -606,12 +606,12 @@ bool ReadFile_geno(const string &file_geno, const set<string> &setSnps,
     return false;
   }
 
-  gsl_vector *genotype = gsl_vector_alloc(W->size1);
-  gsl_vector *genotype_miss = gsl_vector_alloc(W->size1);
-  gsl_matrix *WtW = gsl_matrix_alloc(W->size2, W->size2);
-  gsl_matrix *WtWi = gsl_matrix_alloc(W->size2, W->size2);
-  gsl_vector *Wtx = gsl_vector_alloc(W->size2);
-  gsl_vector *WtWiWtx = gsl_vector_alloc(W->size2);
+  gsl_vector *genotype = gsl_vector_safe_alloc(W->size1);
+  gsl_vector *genotype_miss = gsl_vector_safe_alloc(W->size1);
+  gsl_matrix *WtW = gsl_matrix_safe_alloc(W->size2, W->size2);
+  gsl_matrix *WtWi = gsl_matrix_safe_alloc(W->size2, W->size2);
+  gsl_vector *Wtx = gsl_vector_safe_alloc(W->size2);
+  gsl_vector *WtWiWtx = gsl_vector_safe_alloc(W->size2);
   gsl_permutation *pmt = gsl_permutation_alloc(W->size2);
 
   gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW);
@@ -817,12 +817,12 @@ bool ReadFile_bed(const string &file_bed, const set<string> &setSnps,
     return false;
   }
 
-  gsl_vector *genotype = gsl_vector_alloc(W->size1);
-  gsl_vector *genotype_miss = gsl_vector_alloc(W->size1);
-  gsl_matrix *WtW = gsl_matrix_alloc(W->size2, W->size2);
-  gsl_matrix *WtWi = gsl_matrix_alloc(W->size2, W->size2);
-  gsl_vector *Wtx = gsl_vector_alloc(W->size2);
-  gsl_vector *WtWiWtx = gsl_vector_alloc(W->size2);
+  gsl_vector *genotype = gsl_vector_safe_alloc(W->size1);
+  gsl_vector *genotype_miss = gsl_vector_safe_alloc(W->size1);
+  gsl_matrix *WtW = gsl_matrix_safe_alloc(W->size2, W->size2);
+  gsl_matrix *WtWi = gsl_matrix_safe_alloc(W->size2, W->size2);
+  gsl_vector *Wtx = gsl_vector_safe_alloc(W->size2);
+  gsl_vector *WtWiWtx = gsl_vector_safe_alloc(W->size2);
   gsl_permutation *pmt = gsl_permutation_alloc(W->size2);
 
   gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW);
@@ -1366,12 +1366,12 @@ bool BimbamKin(const string file_geno, const set<string> ksnps,
   bool process_ksnps = ksnps.size();
 
   size_t ni_total = matrix_kin->size1;
-  gsl_vector *geno = gsl_vector_alloc(ni_total);
-  gsl_vector *geno_miss = gsl_vector_alloc(ni_total);
+  gsl_vector *geno = gsl_vector_safe_alloc(ni_total);
+  gsl_vector *geno_miss = gsl_vector_safe_alloc(ni_total);
 
   // Xlarge contains inds x markers
   const size_t msize = K_BATCH_SIZE;
-  gsl_matrix *Xlarge = gsl_matrix_alloc(ni_total, msize);
+  gsl_matrix *Xlarge = gsl_matrix_safe_alloc(ni_total, msize);
   enforce_msg(Xlarge, "allocate Xlarge");
 
   gsl_matrix_set_zero(Xlarge);
@@ -1506,14 +1506,14 @@ bool PlinkKin(const string &file_bed, vector<int> &indicator_snp,
   double d, geno_mean, geno_var;
 
   size_t ni_total = matrix_kin->size1;
-  gsl_vector *geno = gsl_vector_alloc(ni_total);
+  gsl_vector *geno = gsl_vector_safe_alloc(ni_total);
 
   size_t ns_test = 0;
   int n_bit;
 
   // Create a large matrix.
   const size_t msize = K_BATCH_SIZE;
-  gsl_matrix *Xlarge = gsl_matrix_alloc(ni_total, msize);
+  gsl_matrix *Xlarge = gsl_matrix_safe_alloc(ni_total, msize);
   gsl_matrix_set_zero(Xlarge);
 
   // Calculate n_bit and c, the number of bit for each snp.
@@ -1649,8 +1649,8 @@ bool ReadFile_geno(const string file_geno, vector<int> &indicator_idv,
     gsl_matrix_set_zero(K);
   }
 
-  gsl_vector *genotype = gsl_vector_alloc(UtX->size1);
-  gsl_vector *genotype_miss = gsl_vector_alloc(UtX->size1);
+  gsl_vector *genotype = gsl_vector_safe_alloc(UtX->size1);
+  gsl_vector *genotype_miss = gsl_vector_safe_alloc(UtX->size1);
   double geno, geno_mean;
   size_t n_miss;
 
@@ -1760,8 +1760,8 @@ bool ReadFile_geno(const string &file_geno, vector<int> &indicator_idv,
     gsl_matrix_set_zero(K);
   }
 
-  gsl_vector *genotype = gsl_vector_alloc(ni_test);
-  gsl_vector *genotype_miss = gsl_vector_alloc(ni_test);
+  gsl_vector *genotype = gsl_vector_safe_alloc(ni_test);
+  gsl_vector *genotype_miss = gsl_vector_safe_alloc(ni_test);
   double geno, geno_mean;
   size_t n_miss;
 
@@ -1879,7 +1879,7 @@ bool ReadFile_bed(const string &file_bed, vector<int> &indicator_idv,
     gsl_matrix_set_zero(K);
   }
 
-  gsl_vector *genotype = gsl_vector_alloc(UtX->size1);
+  gsl_vector *genotype = gsl_vector_safe_alloc(UtX->size1);
 
   double geno, geno_mean;
   size_t n_miss;
@@ -2015,7 +2015,7 @@ bool ReadFile_bed(const string &file_bed, vector<int> &indicator_idv,
     gsl_matrix_set_zero(K);
   }
 
-  gsl_vector *genotype = gsl_vector_alloc(ni_test);
+  gsl_vector *genotype = gsl_vector_safe_alloc(ni_test);
 
   double geno, geno_mean;
   size_t n_miss;
@@ -2658,13 +2658,13 @@ bool BimbamKinUncentered(const string &file_geno, const set<string> ksnps,
   double d, geno_mean, geno_var;
 
   size_t ni_test = matrix_kin->size1;
-  gsl_vector *geno = gsl_vector_alloc(ni_test);
-  gsl_vector *geno_miss = gsl_vector_alloc(ni_test);
+  gsl_vector *geno = gsl_vector_safe_alloc(ni_test);
+  gsl_vector *geno_miss = gsl_vector_safe_alloc(ni_test);
 
-  gsl_vector *Wtx = gsl_vector_alloc(W->size2);
-  gsl_matrix *WtW = gsl_matrix_alloc(W->size2, W->size2);
-  gsl_matrix *WtWi = gsl_matrix_alloc(W->size2, W->size2);
-  gsl_vector *WtWiWtx = gsl_vector_alloc(W->size2);
+  gsl_vector *Wtx = gsl_vector_safe_alloc(W->size2);
+  gsl_matrix *WtW = gsl_matrix_safe_alloc(W->size2, W->size2);
+  gsl_matrix *WtWi = gsl_matrix_safe_alloc(W->size2, W->size2);
+  gsl_vector *WtWiWtx = gsl_vector_safe_alloc(W->size2);
   gsl_permutation *pmt = gsl_permutation_alloc(W->size2);
 
   gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW);
@@ -2681,7 +2681,7 @@ bool BimbamKinUncentered(const string &file_geno, const set<string> ksnps,
 
   // Create a large matrix.
   const size_t msize = K_BATCH_SIZE;
-  gsl_matrix *Xlarge = gsl_matrix_alloc(ni_test, msize * n_vc);
+  gsl_matrix *Xlarge = gsl_matrix_safe_alloc(ni_test, msize * n_vc);
   gsl_matrix_set_zero(Xlarge);
 
   size_t ns_test = 0;
@@ -2850,12 +2850,12 @@ bool PlinkKin(const string &file_bed, const int display_pace,
 
   size_t ni_test = matrix_kin->size1;
   size_t ni_total = indicator_idv.size();
-  gsl_vector *geno = gsl_vector_alloc(ni_test);
+  gsl_vector *geno = gsl_vector_safe_alloc(ni_test);
 
-  gsl_vector *Wtx = gsl_vector_alloc(W->size2);
-  gsl_matrix *WtW = gsl_matrix_alloc(W->size2, W->size2);
-  gsl_matrix *WtWi = gsl_matrix_alloc(W->size2, W->size2);
-  gsl_vector *WtWiWtx = gsl_vector_alloc(W->size2);
+  gsl_vector *Wtx = gsl_vector_safe_alloc(W->size2);
+  gsl_matrix *WtW = gsl_matrix_safe_alloc(W->size2, W->size2);
+  gsl_matrix *WtWi = gsl_matrix_safe_alloc(W->size2, W->size2);
+  gsl_vector *WtWiWtx = gsl_vector_safe_alloc(W->size2);
   gsl_permutation *pmt = gsl_permutation_alloc(W->size2);
 
   gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW);
@@ -2875,7 +2875,7 @@ bool PlinkKin(const string &file_bed, const int display_pace,
 
   // Create a large matrix.
   const size_t msize = K_BATCH_SIZE;
-  gsl_matrix *Xlarge = gsl_matrix_alloc(ni_test, msize * n_vc);
+  gsl_matrix *Xlarge = gsl_matrix_safe_alloc(ni_test, msize * n_vc);
   gsl_matrix_set_zero(Xlarge);
 
   // Calculate n_bit and c, the number of bit for each SNP.
@@ -3074,8 +3074,8 @@ bool MFILEKin(const size_t mfile_mode, const string &file_mfile,
 
   string file_name;
 
-  gsl_matrix *kin_tmp = gsl_matrix_alloc(matrix_kin->size1, matrix_kin->size2);
-  gsl_vector *ns_tmp = gsl_vector_alloc(vector_ns->size);
+  gsl_matrix *kin_tmp = gsl_matrix_safe_alloc(matrix_kin->size1, matrix_kin->size2);
+  gsl_vector *ns_tmp = gsl_vector_safe_alloc(vector_ns->size);
 
   size_t l = 0;
   double d;
@@ -3843,7 +3843,7 @@ void ReadFile_study(const string &file_study, gsl_matrix *Vq_mat,
   string sfile = file_study + ".size.txt";
   string qfile = file_study + ".q.txt";
 
-  gsl_vector *s = gsl_vector_alloc(s_vec->size + 1);
+  gsl_vector *s = gsl_vector_safe_alloc(s_vec->size + 1);
 
   ReadFile_matrix(Vqfile, Vq_mat);
   ReadFile_vector(sfile, s);
@@ -3868,7 +3868,7 @@ void ReadFile_ref(const string &file_ref, gsl_matrix *S_mat,
   string sfile = file_ref + ".size.txt";
   string Sfile = file_ref + ".S.txt";
 
-  gsl_vector *s = gsl_vector_alloc(s_vec->size + 1);
+  gsl_vector *s = gsl_vector_safe_alloc(s_vec->size + 1);
 
   ReadFile_vector(sfile, s);
   ReadFile_matrix(Sfile, S_mat, Svar_mat);
@@ -3894,9 +3894,9 @@ void ReadFile_mstudy(const string &file_mstudy, gsl_matrix *Vq_mat,
   gsl_vector_set_zero(s_vec);
   ni = 0;
 
-  gsl_matrix *Vq_sub = gsl_matrix_alloc(Vq_mat->size1, Vq_mat->size2);
-  gsl_vector *q_sub = gsl_vector_alloc(q_vec->size);
-  gsl_vector *s = gsl_vector_alloc(s_vec->size + 1);
+  gsl_matrix *Vq_sub = gsl_matrix_safe_alloc(Vq_mat->size1, Vq_mat->size2);
+  gsl_vector *q_sub = gsl_vector_safe_alloc(q_vec->size);
+  gsl_vector *s = gsl_vector_safe_alloc(s_vec->size + 1);
 
   igzstream infile(file_mstudy.c_str(), igzstream::in);
   if (!infile) {
@@ -3985,9 +3985,9 @@ void ReadFile_mref(const string &file_mref, gsl_matrix *S_mat,
   gsl_vector_set_zero(s_vec);
   ni = 0;
 
-  gsl_matrix *S_sub = gsl_matrix_alloc(S_mat->size1, S_mat->size2);
-  gsl_matrix *Svar_sub = gsl_matrix_alloc(Svar_mat->size1, Svar_mat->size2);
-  gsl_vector *s = gsl_vector_alloc(s_vec->size + 1);
+  gsl_matrix *S_sub = gsl_matrix_safe_alloc(S_mat->size1, S_mat->size2);
+  gsl_matrix *Svar_sub = gsl_matrix_safe_alloc(Svar_mat->size1, Svar_mat->size2);
+  gsl_vector *s = gsl_vector_safe_alloc(s_vec->size + 1);
 
   igzstream infile(file_mref.c_str(), igzstream::in);
   if (!infile) {
diff --git a/src/lmm.cpp b/src/lmm.cpp
index 7b32330..e80cd76 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -363,9 +363,9 @@ double LogL_f(double l, void *params) {
   double f = 0.0, logdet_h = 0.0, d;
   size_t index_yy;
 
-  gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size);
+  gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  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_scale(v_temp, l);
@@ -413,11 +413,11 @@ double LogL_dev1(double l, void *params) {
   double dev1 = 0.0, trace_Hi = 0.0;
   size_t index_yy;
 
-  gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_matrix *PPab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *HiHi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size);
+  gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_matrix *PPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size);
+  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_scale(v_temp, l);
@@ -476,13 +476,13 @@ double LogL_dev2(double l, void *params) {
   double dev2 = 0.0, trace_Hi = 0.0, trace_HiHi = 0.0;
   size_t index_yy;
 
-  gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_matrix *PPab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_matrix *PPPab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *HiHi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *HiHiHi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size);
+  gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_matrix *PPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_matrix *PPPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size);
+  gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
+  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_scale(v_temp, l);
@@ -553,13 +553,13 @@ void LogL_dev12(double l, void *params, double *dev1, double *dev2) {
   double trace_Hi = 0.0, trace_HiHi = 0.0;
   size_t index_yy;
 
-  gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_matrix *PPab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_matrix *PPPab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *HiHi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *HiHiHi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size);
+  gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_matrix *PPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_matrix *PPPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size);
+  gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
+  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_scale(v_temp, l);
@@ -636,10 +636,10 @@ double LogRL_f(double l, void *params) {
   double f = 0.0, logdet_h = 0.0, logdet_hiw = 0.0, d;
   size_t index_ww;
 
-  gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_matrix *Iab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size);
+  gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_matrix *Iab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  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_scale(v_temp, l);
@@ -701,11 +701,11 @@ double LogRL_dev1(double l, void *params) {
   double dev1 = 0.0, trace_Hi = 0.0;
   size_t index_ww;
 
-  gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_matrix *PPab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *HiHi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size);
+  gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_matrix *PPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size);
+  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_scale(v_temp, l);
@@ -777,13 +777,13 @@ double LogRL_dev2(double l, void *params) {
   double dev2 = 0.0, trace_Hi = 0.0, trace_HiHi = 0.0;
   size_t index_ww;
 
-  gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_matrix *PPab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_matrix *PPPab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *HiHi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *HiHiHi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size);
+  gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_matrix *PPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_matrix *PPPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size);
+  gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
+  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_scale(v_temp, l);
@@ -867,13 +867,13 @@ void LogRL_dev12(double l, void *params, double *dev1, double *dev2) {
   double trace_Hi = 0.0, trace_HiHi = 0.0;
   size_t index_ww;
 
-  gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_matrix *PPab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_matrix *PPPab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *HiHi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *HiHiHi_eval = gsl_vector_alloc((p->eval)->size);
-  gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size);
+  gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_matrix *PPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_matrix *PPPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size);
+  gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
+  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_scale(v_temp, l);
@@ -947,9 +947,9 @@ void LMM::CalcRLWald(const double &l, const FUNC_PARAM &params, double &beta,
 
   int df = (int)ni_test - (int)n_cvt - 1;
 
-  gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_vector *Hi_eval = gsl_vector_alloc(params.eval->size);
-  gsl_vector *v_temp = gsl_vector_alloc(params.eval->size);
+  gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  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_scale(v_temp, l);
@@ -989,9 +989,9 @@ void LMM::CalcRLScore(const double &l, const FUNC_PARAM &params, double &beta,
 
   int df = (int)ni_test - (int)n_cvt - 1;
 
-  gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_vector *Hi_eval = gsl_vector_alloc(params.eval->size);
-  gsl_vector *v_temp = gsl_vector_alloc(params.eval->size);
+  gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  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_scale(v_temp, l);
@@ -1030,7 +1030,7 @@ void CalcUab(const gsl_matrix *UtW, const gsl_vector *Uty, gsl_matrix *Uab) {
   size_t index_ab;
   size_t n_cvt = UtW->size2;
 
-  gsl_vector *u_a = gsl_vector_alloc(Uty->size);
+  gsl_vector *u_a = gsl_vector_safe_alloc(Uty->size);
 
   for (size_t a = 1; a <= n_cvt + 2; ++a) {
     if (a == n_cvt + 1) {
@@ -1096,8 +1096,8 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
   size_t n_cvt = W->size2;
 
   double d;
-  gsl_vector *v_a = gsl_vector_alloc(y->size);
-  gsl_vector *v_b = gsl_vector_alloc(y->size);
+  gsl_vector *v_a = gsl_vector_safe_alloc(y->size);
+  gsl_vector *v_b = gsl_vector_safe_alloc(y->size);
 
   for (size_t a = 1; a <= n_cvt + 2; ++a) {
     if (a == n_cvt + 1) {
@@ -1141,7 +1141,7 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, const gsl_vector *x,
   size_t n_cvt = W->size2;
 
   double d;
-  gsl_vector *v_b = gsl_vector_alloc(y->size);
+  gsl_vector *v_b = gsl_vector_safe_alloc(y->size);
 
   for (size_t b = 1; b <= n_cvt + 2; ++b) {
     index_ab = GetabIndex(n_cvt + 1, b, n_cvt);
@@ -1188,10 +1188,10 @@ void LMM::AnalyzeGene(const gsl_matrix *U, const gsl_vector *eval,
   // Calculate basic quantities.
   size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2;
 
-  gsl_vector *y = gsl_vector_alloc(U->size1);
-  gsl_vector *Uty = gsl_vector_alloc(U->size2);
-  gsl_matrix *Uab = gsl_matrix_alloc(U->size2, n_index);
-  gsl_vector *ab = gsl_vector_alloc(n_index);
+  gsl_vector *y = gsl_vector_safe_alloc(U->size1);
+  gsl_vector *Uty = gsl_vector_safe_alloc(U->size2);
+  gsl_matrix *Uab = gsl_matrix_safe_alloc(U->size2, n_index);
+  gsl_vector *ab = gsl_vector_safe_alloc(n_index);
 
   // Header.
   getline(infile, line);
@@ -1289,17 +1289,17 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
 
   const size_t inds = U->size1;
   enforce(inds == ni_test);
-  gsl_vector *x = gsl_vector_alloc(inds); // #inds
-  gsl_vector *x_miss = gsl_vector_alloc(inds);
-  gsl_vector *Utx = gsl_vector_alloc(U->size2);
-  gsl_matrix *Uab = gsl_matrix_alloc(U->size2, n_index);
-  gsl_vector *ab = gsl_vector_alloc(n_index);
+  gsl_vector *x = gsl_vector_safe_alloc(inds); // #inds
+  gsl_vector *x_miss = gsl_vector_safe_alloc(inds);
+  gsl_vector *Utx = gsl_vector_safe_alloc(U->size2);
+  gsl_matrix *Uab = gsl_matrix_safe_alloc(U->size2, n_index);
+  gsl_vector *ab = gsl_vector_safe_alloc(n_index);
 
   // Create a large matrix with LMM_BATCH_SIZE columns for batched processing
   // const size_t msize=(process_gwasnps ? 1 : LMM_BATCH_SIZE);
   const size_t msize = LMM_BATCH_SIZE;
-  gsl_matrix *Xlarge = gsl_matrix_alloc(inds, msize);
-  gsl_matrix *UtXlarge = gsl_matrix_alloc(inds, msize);
+  gsl_matrix *Xlarge = gsl_matrix_safe_alloc(inds, msize);
+  gsl_matrix *UtXlarge = gsl_matrix_safe_alloc(inds, msize);
   enforce_msg(Xlarge && UtXlarge, "Xlarge memory check"); // just to be sure
   enforce(Xlarge->size1 == inds);
   gsl_matrix_set_zero(Xlarge);
@@ -1573,10 +1573,10 @@ void MatrixCalcLR(const gsl_matrix *U, const gsl_matrix *UtX,
                   vector<pair<size_t, double>> &pos_loglr) {
   double logl_H0, logl_H1, log_lr, lambda0, lambda1;
 
-  gsl_vector *w = gsl_vector_alloc(Uty->size);
-  gsl_matrix *Utw = gsl_matrix_alloc(Uty->size, 1);
-  gsl_matrix *Uab = gsl_matrix_alloc(Uty->size, 6);
-  gsl_vector *ab = gsl_vector_alloc(6);
+  gsl_vector *w = gsl_vector_safe_alloc(Uty->size);
+  gsl_matrix *Utw = gsl_matrix_safe_alloc(Uty->size, 1);
+  gsl_matrix *Uab = gsl_matrix_safe_alloc(Uty->size, 6);
+  gsl_vector *ab = gsl_vector_safe_alloc(6);
 
   gsl_vector_set_zero(ab);
   gsl_vector_set_all(w, 1.0);
@@ -1781,8 +1781,8 @@ void CalcLambda(const char func_name, const gsl_vector *eval,
   size_t n_cvt = UtW->size2, ni_test = UtW->size1;
   size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2;
 
-  gsl_matrix *Uab = gsl_matrix_alloc(ni_test, n_index);
-  gsl_vector *ab = gsl_vector_alloc(n_index);
+  gsl_matrix *Uab = gsl_matrix_safe_alloc(ni_test, n_index);
+  gsl_vector *ab = gsl_vector_safe_alloc(n_index);
 
   gsl_matrix_set_zero(Uab);
   CalcUab(UtW, Uty, Uab);
@@ -1804,8 +1804,8 @@ void CalcPve(const gsl_vector *eval, const gsl_matrix *UtW,
   size_t n_cvt = UtW->size2, ni_test = UtW->size1;
   size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2;
 
-  gsl_matrix *Uab = gsl_matrix_alloc(ni_test, n_index);
-  gsl_vector *ab = gsl_vector_alloc(n_index);
+  gsl_matrix *Uab = gsl_matrix_safe_alloc(ni_test, n_index);
+  gsl_vector *ab = gsl_vector_safe_alloc(n_index);
 
   gsl_matrix_set_zero(Uab);
   CalcUab(UtW, Uty, Uab);
@@ -1831,15 +1831,15 @@ void CalcLmmVgVeBeta(const gsl_vector *eval, const gsl_matrix *UtW,
   size_t n_cvt = UtW->size2, ni_test = UtW->size1;
   size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2;
 
-  gsl_matrix *Uab = gsl_matrix_alloc(ni_test, n_index);
-  gsl_vector *ab = gsl_vector_alloc(n_index);
-  gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index);
-  gsl_vector *Hi_eval = gsl_vector_alloc(eval->size);
-  gsl_vector *v_temp = gsl_vector_alloc(eval->size);
-  gsl_matrix *HiW = gsl_matrix_alloc(eval->size, UtW->size2);
-  gsl_matrix *WHiW = gsl_matrix_alloc(UtW->size2, UtW->size2);
-  gsl_vector *WHiy = gsl_vector_alloc(UtW->size2);
-  gsl_matrix *Vbeta = gsl_matrix_alloc(UtW->size2, UtW->size2);
+  gsl_matrix *Uab = gsl_matrix_safe_alloc(ni_test, n_index);
+  gsl_vector *ab = gsl_vector_safe_alloc(n_index);
+  gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
+  gsl_vector *Hi_eval = gsl_vector_safe_alloc(eval->size);
+  gsl_vector *v_temp = gsl_vector_safe_alloc(eval->size);
+  gsl_matrix *HiW = gsl_matrix_safe_alloc(eval->size, UtW->size2);
+  gsl_matrix *WHiW = gsl_matrix_safe_alloc(UtW->size2, UtW->size2);
+  gsl_vector *WHiy = gsl_vector_safe_alloc(UtW->size2);
+  gsl_matrix *Vbeta = gsl_matrix_safe_alloc(UtW->size2, UtW->size2);
 
   gsl_matrix_set_zero(Uab);
   CalcUab(UtW, Uty, Uab);
@@ -1921,13 +1921,13 @@ void LMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
   // Calculate basic quantities.
   size_t n_index = (n_cvt + 2 + 2 + 1) * (n_cvt + 2 + 2) / 2;
 
-  gsl_vector *x = gsl_vector_alloc(U->size1);
-  gsl_vector *x_miss = gsl_vector_alloc(U->size1);
-  gsl_vector *Utx = gsl_vector_alloc(U->size2);
-  gsl_matrix *Uab = gsl_matrix_alloc(U->size2, n_index);
-  gsl_vector *ab = gsl_vector_alloc(n_index);
+  gsl_vector *x = gsl_vector_safe_alloc(U->size1);
+  gsl_vector *x_miss = gsl_vector_safe_alloc(U->size1);
+  gsl_vector *Utx = gsl_vector_safe_alloc(U->size2);
+  gsl_matrix *Uab = gsl_matrix_safe_alloc(U->size2, n_index);
+  gsl_vector *ab = gsl_vector_safe_alloc(n_index);
 
-  gsl_matrix *UtW_expand = gsl_matrix_alloc(U->size1, UtW->size2 + 2);
+  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);
@@ -2071,12 +2071,12 @@ void LMM::AnalyzePlinkGXE(const gsl_matrix *U, const gsl_vector *eval,
   // Calculate basic quantities.
   size_t n_index = (n_cvt + 2 + 2 + 1) * (n_cvt + 2 + 2) / 2;
 
-  gsl_vector *x = gsl_vector_alloc(U->size1);
-  gsl_vector *Utx = gsl_vector_alloc(U->size2);
-  gsl_matrix *Uab = gsl_matrix_alloc(U->size2, n_index);
-  gsl_vector *ab = gsl_vector_alloc(n_index);
+  gsl_vector *x = gsl_vector_safe_alloc(U->size1);
+  gsl_vector *Utx = gsl_vector_safe_alloc(U->size2);
+  gsl_matrix *Uab = gsl_matrix_safe_alloc(U->size2, n_index);
+  gsl_vector *ab = gsl_vector_safe_alloc(n_index);
 
-  gsl_matrix *UtW_expand = gsl_matrix_alloc(U->size1, UtW->size2 + 2);
+  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);
diff --git a/src/logistic.cpp b/src/logistic.cpp
index 2dd0402..a936682 100644
--- a/src/logistic.cpp
+++ b/src/logistic.cpp
@@ -7,6 +7,7 @@
 #include <stdio.h>
 
 #include "logistic.h"
+#include "debug.h"
 
 // I need to bundle all the data that goes to the function to optimze
 // together.
@@ -135,7 +136,7 @@ void wgsl_mixed_optim_hessian(const gsl_vector *beta, void *params,
   int K = p->X->size2;
   int Kc = p->Xc->size2;
   int npar = beta->size;
-  gsl_vector *gn = gsl_vector_alloc(npar); // gn
+  gsl_vector *gn = gsl_vector_safe_alloc(npar); // gn
 
   // Intitialize Hessian out necessary ???
   gsl_matrix_set_zero(out);
@@ -226,11 +227,11 @@ int logistic_mixed_fit(gsl_vector *beta, gsl_matrix_int *X,
   // Initial fit.
   mLogLik = wgsl_mixed_optim_f(beta, &p);
 
-  gsl_matrix *myH = gsl_matrix_alloc(npar, npar); // Hessian matrix.
-  gsl_vector *stBeta = gsl_vector_alloc(npar);    // Direction to move.
+  gsl_matrix *myH = gsl_matrix_safe_alloc(npar, npar); // Hessian matrix.
+  gsl_vector *stBeta = gsl_vector_safe_alloc(npar);    // Direction to move.
 
-  gsl_vector *myG = gsl_vector_alloc(npar); // Gradient.
-  gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR.
+  gsl_vector *myG = gsl_vector_safe_alloc(npar); // Gradient.
+  gsl_vector *tau = gsl_vector_safe_alloc(npar); // tau for QR.
 
   for (iter = 0; iter < 100; iter++) {
     wgsl_mixed_optim_hessian(beta, &p, myH); // Calculate Hessian.
@@ -456,11 +457,11 @@ int logistic_cat_fit(gsl_vector *beta, gsl_matrix_int *X, gsl_vector_int *nlev,
   // Initial fit.
   mLogLik = wgsl_cat_optim_f(beta, &p);
 
-  gsl_matrix *myH = gsl_matrix_alloc(npar, npar); // Hessian matrix.
-  gsl_vector *stBeta = gsl_vector_alloc(npar);    // Direction to move.
+  gsl_matrix *myH = gsl_matrix_safe_alloc(npar, npar); // Hessian matrix.
+  gsl_vector *stBeta = gsl_vector_safe_alloc(npar);    // Direction to move.
 
-  gsl_vector *myG = gsl_vector_alloc(npar); // Gradient.
-  gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR.
+  gsl_vector *myG = gsl_vector_safe_alloc(npar); // Gradient.
+  gsl_vector *tau = gsl_vector_safe_alloc(npar); // tau for QR.
 
   for (iter = 0; iter < 100; iter++) {
     wgsl_cat_optim_hessian(beta, &p, myH); // Calculate Hessian.
@@ -596,7 +597,7 @@ void wgsl_cont_optim_hessian(const gsl_vector *beta, void *params,
   int n = p->y->size;
   int Kc = p->Xc->size2;
   int npar = beta->size;
-  gsl_vector *gn = gsl_vector_alloc(npar); // gn.
+  gsl_vector *gn = gsl_vector_safe_alloc(npar); // gn.
 
   // Intitialize Hessian out necessary ???
 
@@ -673,11 +674,11 @@ int logistic_cont_fit(gsl_vector *beta,
   // Initial fit.
   mLogLik = wgsl_cont_optim_f(beta, &p);
 
-  gsl_matrix *myH = gsl_matrix_alloc(npar, npar); // Hessian matrix.
-  gsl_vector *stBeta = gsl_vector_alloc(npar);    // Direction to move.
+  gsl_matrix *myH = gsl_matrix_safe_alloc(npar, npar); // Hessian matrix.
+  gsl_vector *stBeta = gsl_vector_safe_alloc(npar);    // Direction to move.
 
-  gsl_vector *myG = gsl_vector_alloc(npar); // Gradient.
-  gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR.
+  gsl_vector *myG = gsl_vector_safe_alloc(npar); // Gradient.
+  gsl_vector *tau = gsl_vector_safe_alloc(npar); // tau for QR.
 
   for (iter = 0; iter < 100; iter++) {
     wgsl_cont_optim_hessian(beta, &p, myH); // Calculate Hessian.
diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index f55f3c6..e7dff73 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -80,8 +80,8 @@ double VectorVar(const gsl_vector *v) {
 // Center the matrix G.
 void CenterMatrix(gsl_matrix *G) {
   double d;
-  gsl_vector *w = gsl_vector_alloc(G->size1);
-  gsl_vector *Gw = gsl_vector_alloc(G->size1);
+  gsl_vector *w = gsl_vector_safe_alloc(G->size1);
+  gsl_vector *Gw = gsl_vector_safe_alloc(G->size1);
   gsl_vector_set_all(w, 1.0);
 
   gsl_blas_dgemv(CblasNoTrans, 1.0, G, w, 0.0, Gw);
@@ -105,7 +105,7 @@ void CenterMatrix(gsl_matrix *G) {
 // Center the matrix G.
 void CenterMatrix(gsl_matrix *G, const gsl_vector *w) {
   double d, wtw;
-  gsl_vector *Gw = gsl_vector_alloc(G->size1);
+  gsl_vector *Gw = gsl_vector_safe_alloc(G->size1);
 
   gsl_blas_ddot(w, w, &wtw);
   gsl_blas_dgemv(CblasNoTrans, 1.0, G, w, 0.0, Gw);
@@ -127,12 +127,12 @@ void CenterMatrix(gsl_matrix *G, const gsl_vector *w) {
 
 // Center the matrix G.
 void CenterMatrix(gsl_matrix *G, const gsl_matrix *W) {
-  gsl_matrix *WtW = gsl_matrix_alloc(W->size2, W->size2);
-  gsl_matrix *WtWi = gsl_matrix_alloc(W->size2, W->size2);
-  gsl_matrix *WtWiWt = gsl_matrix_alloc(W->size2, G->size1);
-  gsl_matrix *GW = gsl_matrix_alloc(G->size1, W->size2);
-  gsl_matrix *WtGW = gsl_matrix_alloc(W->size2, W->size2);
-  gsl_matrix *Gtmp = gsl_matrix_alloc(G->size1, G->size1);
+  gsl_matrix *WtW = gsl_matrix_safe_alloc(W->size2, W->size2);
+  gsl_matrix *WtWi = gsl_matrix_safe_alloc(W->size2, W->size2);
+  gsl_matrix *WtWiWt = gsl_matrix_safe_alloc(W->size2, G->size1);
+  gsl_matrix *GW = gsl_matrix_safe_alloc(G->size1, W->size2);
+  gsl_matrix *WtGW = gsl_matrix_safe_alloc(W->size2, W->size2);
+  gsl_matrix *Gtmp = gsl_matrix_safe_alloc(G->size1, G->size1);
 
   gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW);
 
@@ -227,7 +227,7 @@ bool isMatrixSymmetric(const gsl_matrix *G) {
 
 bool isMatrixPositiveDefinite(const gsl_matrix *G) {
   enforce(G->size1 == G->size2);
-  auto G2 = gsl_matrix_alloc(G->size1, G->size2);
+  auto G2 = gsl_matrix_safe_alloc(G->size1, G->size2);
   enforce_gsl(gsl_matrix_memcpy(G2,G));
   auto handler = gsl_set_error_handler_off();
 #if GSL_MAJOR_VERSION >= 2 && GSL_MINOR_VERSION >= 3
@@ -242,11 +242,11 @@ bool isMatrixPositiveDefinite(const gsl_matrix *G) {
 
 gsl_vector *getEigenValues(const gsl_matrix *G) {
   enforce(G->size1 == G->size2);
-  auto G2 = gsl_matrix_alloc(G->size1, G->size2);
+  auto G2 = gsl_matrix_safe_alloc(G->size1, G->size2);
   enforce_gsl(gsl_matrix_memcpy(G2,G));
   auto eworkspace = gsl_eigen_symm_alloc(G->size1);
   enforce(eworkspace);
-  gsl_vector *eigenvalues = gsl_vector_alloc(G->size1);
+  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);
@@ -345,9 +345,9 @@ double CenterVector(gsl_vector *y) {
 
 // Center the vector y.
 void CenterVector(gsl_vector *y, const gsl_matrix *W) {
-  gsl_matrix *WtW = gsl_matrix_alloc(W->size2, W->size2);
-  gsl_vector *Wty = gsl_vector_alloc(W->size2);
-  gsl_vector *WtWiWty = gsl_vector_alloc(W->size2);
+  gsl_matrix *WtW = gsl_matrix_safe_alloc(W->size2, W->size2);
+  gsl_vector *Wty = gsl_vector_safe_alloc(W->size2);
+  gsl_vector *WtWiWty = gsl_vector_safe_alloc(W->size2);
 
   gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW);
   gsl_blas_dgemv(CblasTrans, 1.0, W, y, 0.0, Wty);
@@ -387,7 +387,7 @@ void StandardizeVector(gsl_vector *y) {
 
 // Calculate UtX.
 void CalcUtX(const gsl_matrix *U, gsl_matrix *UtX) {
-  gsl_matrix *X = gsl_matrix_alloc(UtX->size1, UtX->size2);
+  gsl_matrix *X = gsl_matrix_safe_alloc(UtX->size1, UtX->size2);
   gsl_matrix_memcpy(X, UtX);
   fast_dgemm("T", "N", 1.0, U, X, 0.0, UtX);
   gsl_matrix_free(X);