Refactored debug settings

Replaced eigenlib_dgemm with fast_dgemm - 10-30% speed gain for GEMMA

12 files changed, 128 insertions, 104 deletions

diff --git a/src/debug.cpp b/src/debug.cpp
index dacb89d..eb5c041 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -19,36 +19,46 @@
 #include "mathfunc.h"
 
 static bool debug_mode     = false;
-static bool debug_no_check = false;
-static bool debug_strict   = false;
+static bool debug_check    = false; // check data/algorithms
+static bool debug_strict   = false; // fail on error
+static bool debug_quiet    = false;
+static uint debug_issue    = 0;     // github issues
+static bool debug_legacy   = false; // legacy mode
 
 void debug_set_debug_mode(bool setting) { debug_mode = setting; }
-void debug_set_no_check_mode(bool setting) {debug_no_check = setting; }
+void debug_set_no_check_mode(bool setting) {debug_check = !setting; }
 void debug_set_strict_mode(bool setting) { debug_strict = setting; }
+void debug_set_quiet_mode(bool setting) { debug_quiet = setting; }
+void debug_set_issue(uint issue) { debug_issue = issue; }
+void debug_set_legacy_mode(bool setting) { debug_legacy = setting; }
 
 bool is_debug_mode() { return debug_mode; };
-bool is_no_check_mode() { return debug_no_check; };
+bool is_no_check_mode() { return !debug_check; };
+bool is_check_mode() { return debug_check; };
 bool is_strict_mode() { return debug_strict; };
+bool is_quiet_mode() { return debug_quiet; };
+bool is_issue(uint issue) { return issue == debug_issue; };
+bool is_legacy_mode() { return debug_legacy; };
 
 /*
   Helper function to make sure gsl allocations do their job because
   gsl_matrix_alloc does not initiatize values (behaviour that changed
   in GSL2) we introduced a 'strict mode' by initializing the buffer
-  with NaNs. This happens in STRICT mode without NO-CHECKS
-  (i.e. -strict option).
+  with NaNs. This happens when NO-CHECKS is not set
+  (i.e. -no-checks option).
 */
 gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols) {
   gsl_matrix *m = gsl_matrix_alloc(rows,cols);
   enforce_msg(m,"Not enough memory"); // just to be sure when there is no error handler set
-  if (debug_strict && !debug_no_check) {
+  if (is_check_mode()) {
     gsl_matrix_set_all(m, nan(""));
   }
   return m;
 }
 
 // Helper function called by macro validate_K(K, check)
-void do_validate_K(const gsl_matrix *K, bool do_check, bool strict, const char *__file, int __line) {
-  if (do_check) {
+void do_validate_K(const gsl_matrix *K, const char *__file, int __line) {
+  if (is_check_mode()) {
     // debug_msg("Validating K");
     auto eigenvalues = getEigenValues(K);
     const uint count_small = count_small_values(eigenvalues,EIGEN_MINVALUE);
@@ -61,13 +71,13 @@ void do_validate_K(const gsl_matrix *K, bool do_check, bool strict, const char *
     if (!isMatrixIllConditioned(eigenvalues))
       warning_at_msg(__file,__line,"K is ill conditioned!");
     if (!isMatrixSymmetric(K))
-      fail_at_msg(strict,__file,__line,"K is not symmetric!" );
+      fail_at_msg(is_strict_mode(),__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))
-      fail_at_msg(strict,__file,__line,"K is not positive definite!");
+      fail_at_msg(is_strict_mode(),__file,__line,"K is not positive definite!");
     gsl_vector_free(eigenvalues);
   }
 }
diff --git a/src/debug.h b/src/debug.h
index 0f7b38f..ee36887 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -13,18 +13,25 @@ void gemma_gsl_error_handler (const char * reason,
 void debug_set_debug_mode(bool setting);
 void debug_set_no_check_mode(bool setting);
 void debug_set_strict_mode(bool setting);
+void debug_set_quiet_mode(bool setting);
+void debug_set_issue(uint issue);
+void debug_set_legacy_mode(bool setting);
 
 bool is_debug_mode();
 bool is_no_check_mode();
+bool is_check_mode();
 bool is_strict_mode();
+bool is_quiet_mode();
+bool is_issue(uint issue);
+bool is_legacy_mode();
 
 gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols);
 
 // Validation routines
-void do_validate_K(const gsl_matrix *K, bool do_check, bool strict, const char *__file, int __line);
+void do_validate_K(const gsl_matrix *K, const char *__file, int __line);
 
 #define ROUND(f) round(f * 10000.)/10000
-#define validate_K(K,check,strict) do_validate_K(K,check,strict,__FILE__,__LINE__)
+#define validate_K(K) do_validate_K(K,__FILE__,__LINE__)
 
 #define warning_at_msg(__file,__line,msg) cerr << "**** WARNING: " << msg << " in " << __file << " at line " << __line << endl;
 
diff --git a/src/fastblas.cpp b/src/fastblas.cpp
index b25d964..0a4eba3 100644
--- a/src/fastblas.cpp
+++ b/src/fastblas.cpp
@@ -44,14 +44,16 @@ gsl_matrix *fast_copy(gsl_matrix *m, const double *mem) {
       }
     }
   } else { // faster goes by row
+    auto v = gsl_vector_calloc(cols);
+    enforce(v); // just to be sure
     for (auto r=0; r<rows; r++) {
-      auto v = gsl_vector_calloc(cols);
       assert(v->size == cols);
       assert(v->block->size == cols);
       assert(v->stride == 1);
       memcpy(v->block->data,&mem[r*cols],cols*sizeof(double));
       gsl_matrix_set_row(m,r,v);
     }
+    gsl_vector_free(v);
   }
   return m;
 }
@@ -160,28 +162,25 @@ static void fast_cblas_dgemm(const char *TransA, const char *TransB, const doubl
   // C++ is row-major
   auto transA = (*TransA == 'N' || *TransA == 'n' ? CblasNoTrans : CblasTrans);
   auto transB = (*TransB == 'N' || *TransB == 'n' ? CblasNoTrans : CblasTrans);
-  // A(m x k) * B(k x n) = C(m x n))
-  auto rowsA = A->size1;
-  auto colsA = A->size2;
-  blasint M = A->size1;
-  blasint K = B->size1;
-  assert(K == colsA);
-  blasint N = B->size2;
-  // cout << M << "," << N "," << K << endl;
-  // Layout = CblasRowMajor: Trans: K , NoTrans M
-  blasint lda = (transA==CblasNoTrans ? K : M );
-  blasint ldb = (transB==CblasNoTrans ? N : K );
-  blasint ldc = N;
-
-  fast_cblas_dgemm(CblasRowMajor, transA, transB, M, N, K, alpha,
-              /* A */ A->data,
-              /* lda */ lda,
-              /* B */ B->data,
-              /* ldb */ ldb,
-              /* beta */ beta,
-              /* C */ C->data, ldc);
+  const size_t M = C->size1;
+  const size_t N = C->size2;
+  const size_t MA = (transA == CblasNoTrans) ? A->size1 : A->size2;
+  const size_t NA = (transA == CblasNoTrans) ? A->size2 : A->size1;
+  const size_t MB = (transB == CblasNoTrans) ? B->size1 : B->size2;
+  const size_t NB = (transB == CblasNoTrans) ? B->size2 : B->size1;
+
+  if (M == MA && N == NB && NA == MB) {  /* [MxN] = [MAxNA][MBxNB] */
+
+    cblas_dgemm (CblasRowMajor, transA, transB, M, N,NA,
+                 alpha, A->data, A->tda, B->data, B->tda, beta,
+                 C->data, C->tda);
+
+  } else {
+    throw invalid_argument("Range error in dgemm");
+  }
 }
 
+
 /*
    Use the fasted/supported way to call BLAS dgemm
 */
@@ -189,14 +188,19 @@ static void fast_cblas_dgemm(const char *TransA, const char *TransB, const doubl
 void fast_dgemm(const char *TransA, const char *TransB, const double alpha,
                 const gsl_matrix *A, const gsl_matrix *B, const double beta,
                 gsl_matrix *C) {
-  fast_cblas_dgemm(TransA,TransB,alpha,A,B,beta,C);
-
-  #ifndef NDEBUG
-  if (is_strict_mode() && !is_no_check_mode()) {
-    // ---- validate with original implementation
-    gsl_matrix *C1 = gsl_matrix_alloc(C->size1,C->size2);
-    eigenlib_dgemm(TransA,TransB,alpha,A,B,beta,C1);
-    enforce_msg(gsl_matrix_equal(C,C1),"dgemm outcomes are not equal for fast & eigenlib");
+  if (is_legacy_mode()) {
+    eigenlib_dgemm(TransA,TransB,alpha,A,B,beta,C);
+  } else {
+    fast_cblas_dgemm(TransA,TransB,alpha,A,B,beta,C);
+
+    #ifndef NDEBUG
+    if (is_check_mode()) {
+      // ---- validate with original implementation
+      gsl_matrix *C1 = gsl_matrix_alloc(C->size1,C->size2);
+      eigenlib_dgemm(TransA,TransB,alpha,A,B,beta,C1);
+      enforce_msg(gsl_matrix_equal(C,C1),"dgemm outcomes are not equal for fast & eigenlib");
+      gsl_matrix_free(C1);
+    }
+    #endif
   }
-  #endif
 }
diff --git a/src/gemma.cpp b/src/gemma.cpp
index 95630c6..e2881a4 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -721,6 +721,7 @@ void GEMMA::PrintHelp(size_t option) {
     cout << " -debug                   debug output" << endl;
     cout << " -nind       [num]        read up to num individuals" << endl;
     cout << " -issue      [num]        enable tests relevant to issue tracker" << endl;
+    cout << " -legacy                  run gemma in legacy mode" << endl;
     cout << endl;
   }
 
@@ -766,7 +767,7 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) {
       str.assign(argv[i]);
       cPar.file_mbfile = str;
     } else if (strcmp(argv[i], "-silence") == 0) {
-      cPar.mode_silence = true;
+      debug_set_quiet_mode(true);
     } else if (strcmp(argv[i], "-g") == 0) {
       if (argv[i + 1] == NULL || argv[i + 1][0] == '-') {
         continue;
@@ -1367,8 +1368,9 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) {
       ++i;
       str.clear();
       str.assign(argv[i]);
-      cPar.issue = atoi(str.c_str()); // for testing purposes
-      enforce(cPar.issue > 0);
+      auto issue = atoi(str.c_str()); // for testing purposes
+      enforce(issue > 0);
+      debug_set_issue(issue);
     } else if (strcmp(argv[i], "-emp") == 0) {
       if (argv[i + 1] == NULL || argv[i + 1][0] == '-') {
         continue;
@@ -1588,14 +1590,16 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) {
       str.assign(argv[i]);
       cPar.window_ns = atoi(str.c_str());
     } else if (strcmp(argv[i], "-debug") == 0) {
-      cPar.mode_debug = true;
+      // cPar.mode_debug = true;
       debug_set_debug_mode(true);
     } else if (strcmp(argv[i], "-no-check") == 0) {
-      cPar.mode_check = false;
+      // cPar.mode_check = false;
       debug_set_no_check_mode(true);
     } else if (strcmp(argv[i], "-strict") == 0) {
-      cPar.mode_strict = true;
+      // cPar.mode_strict = true;
       debug_set_strict_mode(true);
+    } else if (strcmp(argv[i], "-legacy") == 0) {
+      debug_set_legacy_mode(true);
     } else {
       cout << "error! unrecognized option: " << argv[i] << endl;
       cPar.error = true;
@@ -1742,7 +1746,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     // center matrix G
     CenterMatrix(G);
     CenterMatrix(G_full);
-    validate_K(G,cPar.mode_check,cPar.mode_strict);
+    validate_K(G);
 
     // eigen-decomposition and calculate trace_G
     cout << "Start Eigen-Decomposition..." << endl;
@@ -1882,7 +1886,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     }
 
     // Now we have the Kinship matrix test it
-    validate_K(G,cPar.mode_check,cPar.mode_strict);
+    validate_K(G);
 
     if (cPar.a_mode == 21) {
       cPar.WriteMatrix(G, "cXX");
@@ -2323,7 +2327,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
         // center matrix G
         CenterMatrix(G);
-        validate_K(G,cPar.mode_check,cPar.mode_strict);
+        validate_K(G);
 
         (cPar.v_traceG).clear();
         double d = 0;
@@ -2532,7 +2536,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     gsl_vector *eval = gsl_vector_calloc(Y->size1);
     gsl_vector *env = gsl_vector_alloc(Y->size1);
     gsl_vector *weight = gsl_vector_alloc(Y->size1);
-    assert_issue(cPar.issue == 26, UtY->data[0] == 0.0);
+    assert_issue(is_issue(26), UtY->data[0] == 0.0);
 
     // set covariates matrix W and phenotype matrix Y
     // an intercept should be included in W,
@@ -2552,7 +2556,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
       // center matrix G
       CenterMatrix(G);
-      validate_K(G,cPar.mode_check,cPar.mode_strict);
+      validate_K(G);
 
       // is residual weights are provided, then
       if (!cPar.file_weight.empty()) {
@@ -2633,7 +2637,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
       CalcUtX(U, W, UtW);
       CalcUtX(U, Y, UtY);
 
-      assert_issue(cPar.issue == 26, ROUND(UtY->data[0]) == -16.6143);
+      assert_issue(is_issue(26), ROUND(UtY->data[0]) == -16.6143);
 
       LMM cLmm;
       cLmm.CopyFromParam(cPar);
@@ -2650,7 +2654,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
       // calculate UtW and Uty
       CalcUtX(U, W, UtW);
       CalcUtX(U, Y, UtY);
-      assert_issue(cPar.issue == 26, ROUND(UtY->data[0]) == -16.6143);
+      assert_issue(is_issue(26), ROUND(UtY->data[0]) == -16.6143);
 
       // calculate REMLE/MLE estimate and pve for univariate model
       if (cPar.n_ph == 1) { // one phenotype
@@ -2658,7 +2662,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
         gsl_vector_view se_beta = gsl_matrix_row(se_B, 0);
         gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0);
 
-        assert_issue(cPar.issue == 26, ROUND(UtY->data[0]) == -16.6143);
+        assert_issue(is_issue(26), ROUND(UtY->data[0]) == -16.6143);
 
         CalcLambda('L', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max,
                    cPar.n_region, cPar.l_mle_null, cPar.logl_mle_H0);
@@ -2857,7 +2861,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
         // center matrix G
         CenterMatrix(G);
-        validate_K(G,cPar.mode_check,cPar.mode_strict);
+        validate_K(G);
       } else {
         cPar.ReadGenotypes(UtX, G, true);
       }
@@ -2968,7 +2972,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
           // center matrix G
           CenterMatrix(G);
-          validate_K(G,cPar.mode_check,cPar.mode_strict);
+          validate_K(G);
 
         } else {
           cPar.ReadGenotypes(UtX, G, true);
diff --git a/src/io.cpp b/src/io.cpp
index 84367ed..cc0f064 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -41,6 +41,7 @@
 
 #include "debug.h"
 #include "eigenlib.h"
+#include "fastblas.h"
 #include "gzstream.h"
 #include "io.h"
 #include "lapack.h"
@@ -594,7 +595,7 @@ bool ReadFile_geno(const string &file_geno, const set<string> &setSnps,
                    const double &r2_level, map<string, string> &mapRS2chr,
                    map<string, long int> &mapRS2bp,
                    map<string, double> &mapRS2cM, vector<SNPINFO> &snpInfo,
-                   size_t &ns_test, bool debug) {
+                   size_t &ns_test) {
   debug_msg("entered");
   indicator_snp.clear();
   snpInfo.clear();
@@ -667,7 +668,7 @@ bool ReadFile_geno(const string &file_geno, const set<string> &setSnps,
     }
 
     if (mapRS2bp.count(rs) == 0) {
-      if (debug && count_warnings++ < 10) {
+      if (is_debug_mode() && count_warnings++ < 10) {
         std::string msg = "Can't figure out position for ";
         msg += rs;
         debug_msg(msg);
@@ -1454,12 +1455,12 @@ bool BimbamKin(const string file_geno, const set<string> ksnps,
 
     // compute kinship matrix and return in matrix_kin a SNP at a time
     if (ns_test % msize == 0) {
-      eigenlib_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
+      fast_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
       gsl_matrix_set_zero(Xlarge);
     }
   }
   if (ns_test % msize != 0) {
-    eigenlib_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
+    fast_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
   }
   cout << endl;
 
@@ -1600,13 +1601,13 @@ bool PlinkKin(const string &file_bed, vector<int> &indicator_snp,
     ns_test++;
 
     if (ns_test % msize == 0) {
-      eigenlib_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
+      fast_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
       gsl_matrix_set_zero(Xlarge);
     }
   }
 
   if (ns_test % msize != 0) {
-    eigenlib_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
+    fast_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
   }
 
   cout << endl;
@@ -1633,7 +1634,7 @@ bool PlinkKin(const string &file_bed, vector<int> &indicator_snp,
 // genotype and calculate K.
 bool ReadFile_geno(const string file_geno, vector<int> &indicator_idv,
                    vector<int> &indicator_snp, gsl_matrix *UtX, gsl_matrix *K,
-                   const bool calc_K, bool debug) {
+                   const bool calc_K) {
   debug_msg("entered");
   igzstream infile(file_geno.c_str(), igzstream::in);
   if (!infile) {
@@ -1738,7 +1739,7 @@ bool ReadFile_geno(const string &file_geno, vector<int> &indicator_idv,
                    vector<int> &indicator_snp,
                    vector<vector<unsigned char>> &Xt, gsl_matrix *K,
                    const bool calc_K, const size_t ni_test,
-                   const size_t ns_test, bool debug) {
+                   const size_t ns_test) {
   debug_msg("entered");
   igzstream infile(file_geno.c_str(), igzstream::in);
   if (!infile) {
@@ -2757,7 +2758,7 @@ bool BimbamKinUncentered(const string &file_geno, const set<string> ksnps,
         ns_vec[0]++;
 
         if (ns_vec[0] % msize == 0) {
-          eigenlib_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
+          fast_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
           gsl_matrix_set_zero(Xlarge);
         }
       } else if (mapRS2cat.count(rs) != 0) {
@@ -2774,7 +2775,7 @@ bool BimbamKinUncentered(const string &file_geno, const set<string> ksnps,
               gsl_matrix_submatrix(Xlarge, 0, msize * i_vc, ni_test, msize);
           gsl_matrix_view kin_sub = gsl_matrix_submatrix(
               matrix_kin, 0, ni_test * i_vc, ni_test, ni_test);
-          eigenlib_dgemm("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0,
+          fast_dgemm("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0,
                          &kin_sub.matrix);
 
           gsl_matrix_set_zero(&X_sub.matrix);
@@ -2790,7 +2791,7 @@ bool BimbamKinUncentered(const string &file_geno, const set<string> ksnps,
           gsl_matrix_submatrix(Xlarge, 0, msize * i_vc, ni_test, msize);
       gsl_matrix_view kin_sub =
           gsl_matrix_submatrix(matrix_kin, 0, ni_test * i_vc, ni_test, ni_test);
-      eigenlib_dgemm("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0,
+      fast_dgemm("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0,
                      &kin_sub.matrix);
     }
   }
@@ -2983,7 +2984,7 @@ bool PlinkKin(const string &file_bed, const int display_pace,
         ns_vec[0]++;
 
         if (ns_vec[0] % msize == 0) {
-          eigenlib_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
+          fast_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
           gsl_matrix_set_zero(Xlarge);
         }
       } else if (mapRS2cat.count(rs) != 0) {
@@ -3000,7 +3001,7 @@ bool PlinkKin(const string &file_bed, const int display_pace,
               gsl_matrix_submatrix(Xlarge, 0, msize * i_vc, ni_test, msize);
           gsl_matrix_view kin_sub = gsl_matrix_submatrix(
               matrix_kin, 0, ni_test * i_vc, ni_test, ni_test);
-          eigenlib_dgemm("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0,
+          fast_dgemm("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0,
                          &kin_sub.matrix);
 
           gsl_matrix_set_zero(&X_sub.matrix);
diff --git a/src/io.h b/src/io.h
index c65a2b0..215e8ba 100644
--- a/src/io.h
+++ b/src/io.h
@@ -64,7 +64,7 @@ bool ReadFile_geno(const string &file_geno, const set<string> &setSnps,
                    const double &r2_level, map<string, string> &mapRS2chr,
                    map<string, long int> &mapRS2bp,
                    map<string, double> &mapRS2cM, vector<SNPINFO> &snpInfo,
-                   size_t &ns_test, bool debug);
+                   size_t &ns_test);
 bool ReadFile_bed(const string &file_bed, const set<string> &setSnps,
                   const gsl_matrix *W, vector<int> &indicator_idv,
                   vector<int> &indicator_snp, vector<SNPINFO> &snpInfo,
@@ -94,7 +94,7 @@ bool PlinkKin(const string &file_bed, vector<int> &indicator_snp,
 
 bool ReadFile_geno(const string file_geno, vector<int> &indicator_idv,
                    vector<int> &indicator_snp, gsl_matrix *UtX, gsl_matrix *K,
-                   const bool calc_K, bool debug);
+                   const bool calc_K);
 bool ReadFile_bed(const string &file_bed, vector<int> &indicator_idv,
                   vector<int> &indicator_snp, gsl_matrix *UtX, gsl_matrix *K,
                   const bool calc_K);
@@ -102,7 +102,7 @@ bool ReadFile_geno(const string &file_geno, vector<int> &indicator_idv,
                    vector<int> &indicator_snp,
                    vector<vector<unsigned char>> &Xt, gsl_matrix *K,
                    const bool calc_K, const size_t ni_test,
-                   const size_t ns_test, bool debug);
+                   const size_t ns_test);
 bool ReadFile_bed(const string &file_bed, vector<int> &indicator_idv,
                   vector<int> &indicator_snp, vector<vector<unsigned char>> &Xt,
                   gsl_matrix *K, const bool calc_K, const size_t ni_test,
diff --git a/src/lmm.h b/src/lmm.h
index c02325d..9c46fae 100644
--- a/src/lmm.h
+++ b/src/lmm.h
@@ -28,7 +28,7 @@
 
 using namespace std;
 
-#define LMM_BATCH_SIZE 1000 // used for batch processing
+#define LMM_BATCH_SIZE 20000 // used for batch processing
 
 class FUNC_PARAM {
 
diff --git a/src/main.cpp b/src/main.cpp
index 92c4d90..d752a72 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -61,7 +61,7 @@ int main(int argc, char *argv[]) {
     return EXIT_FAILURE;
   }
 
-  if (cPar.mode_silence) {
+  if (is_quiet_mode()) {
     stringstream ss;
     cout.rdbuf(ss.rdbuf());
   }
diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index fbbc061..d0a9d65 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -49,6 +49,7 @@
 
 #include "debug.h"
 #include "eigenlib.h"
+#include "fastblas.h"
 #include "lapack.h"
 #include "mathfunc.h"
 
@@ -388,20 +389,16 @@ void StandardizeVector(gsl_vector *y) {
 void CalcUtX(const gsl_matrix *U, gsl_matrix *UtX) {
   gsl_matrix *X = gsl_matrix_alloc(UtX->size1, UtX->size2);
   gsl_matrix_memcpy(X, UtX);
-  eigenlib_dgemm("T", "N", 1.0, U, X, 0.0, UtX);
+  fast_dgemm("T", "N", 1.0, U, X, 0.0, UtX);
   gsl_matrix_free(X);
-
-  return;
 }
 
 void CalcUtX(const gsl_matrix *U, const gsl_matrix *X, gsl_matrix *UtX) {
-  eigenlib_dgemm("T", "N", 1.0, U, X, 0.0, UtX);
-  return;
+  fast_dgemm("T", "N", 1.0, U, X, 0.0, UtX);
 }
 
 void CalcUtX(const gsl_matrix *U, const gsl_vector *x, gsl_vector *Utx) {
   gsl_blas_dgemv(CblasTrans, 1.0, U, x, 0.0, Utx);
-  return;
 }
 
 // Kronecker product.
diff --git a/src/mvlmm.cpp b/src/mvlmm.cpp
index c0dc143..95ef14b 100644
--- a/src/mvlmm.cpp
+++ b/src/mvlmm.cpp
@@ -39,6 +39,7 @@
 #include "gsl/gsl_vector.h"
 
 #include "eigenlib.h"
+#include "fastblas.h"
 #include "gzstream.h"
 #include "io.h"
 #include "lapack.h"
@@ -3250,8 +3251,8 @@ void MVLMM::AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval,
           gsl_matrix_submatrix(UtXlarge, 0, 0, UtXlarge->size1, l);
 
       time_start = clock();
-      eigenlib_dgemm("T", "N", 1.0, U, &Xlarge_sub.matrix, 0.0,
-                     &UtXlarge_sub.matrix);
+      fast_dgemm("T", "N", 1.0, U, &Xlarge_sub.matrix, 0.0,
+                 &UtXlarge_sub.matrix);
       time_UtX += (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0);
 
       gsl_matrix_set_zero(Xlarge);
@@ -3717,7 +3718,7 @@ void MVLMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval,
           gsl_matrix_submatrix(UtXlarge, 0, 0, UtXlarge->size1, l);
 
       time_start = clock();
-      eigenlib_dgemm("T", "N", 1.0, U, &Xlarge_sub.matrix, 0.0,
+      fast_dgemm("T", "N", 1.0, U, &Xlarge_sub.matrix, 0.0,
                      &UtXlarge_sub.matrix);
       time_UtX += (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0);
 
diff --git a/src/param.cpp b/src/param.cpp
index d637481..8aa7a64 100644
--- a/src/param.cpp
+++ b/src/param.cpp
@@ -66,7 +66,7 @@ void LOCO_set_Snps(set<string> &ksnps, set<string> &gwasnps,
 // (indicator_idv[x] == 1). This should match indicator_cvt etc. If
 // this gives problems with certain sets we can simply trim to size.
 
-void trim_individuals(vector<int> &idvs, size_t ni_max, bool debug) {
+void trim_individuals(vector<int> &idvs, size_t ni_max) {
   if (ni_max) {
     size_t count = 0;
     for (auto ind = idvs.begin(); ind != idvs.end(); ++ind) {
@@ -76,7 +76,7 @@ void trim_individuals(vector<int> &idvs, size_t ni_max, bool debug) {
         break;
     }
     if (count != idvs.size()) {
-      if (debug)
+      if (is_debug_mode())
         cout << "**** TEST MODE: trim individuals from " << idvs.size()
              << " to " << count << endl;
       idvs.resize(count);
@@ -87,7 +87,7 @@ void trim_individuals(vector<int> &idvs, size_t ni_max, bool debug) {
 // ---- PARAM class implementation
 
 PARAM::PARAM(void)
-    : mode_silence(false), a_mode(0), k_mode(1), d_pace(DEFAULT_PACE),
+    : a_mode(0), k_mode(1), d_pace(DEFAULT_PACE),
       file_out("result"), path_out("./output/"), miss_level(0.05),
       maf_level(0.01), hwe_level(0), r2_level(0.9999), l_min(1e-5), l_max(1e5),
       n_region(10), p_nr(0.001), em_prec(0.0001), nr_prec(0.0001),
@@ -221,7 +221,7 @@ void PARAM::ReadFiles(void) {
   } else {
     n_cvt = 1;
   }
-  trim_individuals(indicator_cvt, ni_max, mode_debug);
+  trim_individuals(indicator_cvt, ni_max);
 
   if (!file_gxe.empty()) {
     if (ReadFile_column(file_gxe, indicator_gxe, gxe, 1) == false) {
@@ -234,7 +234,7 @@ void PARAM::ReadFiles(void) {
     }
   }
 
-  trim_individuals(indicator_idv, ni_max, mode_debug);
+  trim_individuals(indicator_idv, ni_max);
 
   // Read genotype and phenotype file for PLINK format.
   if (!file_bfile.empty()) {
@@ -302,11 +302,11 @@ void PARAM::ReadFiles(void) {
     gsl_matrix *W = gsl_matrix_alloc(ni_test, n_cvt);
     CopyCvt(W);
 
-    trim_individuals(indicator_idv, ni_max, mode_debug);
-    trim_individuals(indicator_cvt, ni_max, mode_debug);
+    trim_individuals(indicator_idv, ni_max);
+    trim_individuals(indicator_cvt, ni_max);
     if (ReadFile_geno(file_geno, setSnps, W, indicator_idv, indicator_snp,
                       maf_level, miss_level, hwe_level, r2_level, mapRS2chr,
-                      mapRS2bp, mapRS2cM, snpInfo, ns_test, mode_debug) == false) {
+                      mapRS2bp, mapRS2cM, snpInfo, ns_test) == false) {
       error = true;
     }
     gsl_matrix_free(W);
@@ -416,7 +416,7 @@ void PARAM::ReadFiles(void) {
     while (!safeGetline(infile, file_name).eof()) {
       if (ReadFile_geno(file_name, setSnps, W, indicator_idv, indicator_snp,
                         maf_level, miss_level, hwe_level, r2_level, mapRS2chr,
-                        mapRS2bp, mapRS2cM, snpInfo, ns_test_tmp, mode_debug) == false) {
+                        mapRS2bp, mapRS2cM, snpInfo, ns_test_tmp) == false) {
         error = true;
       }
 
@@ -1282,7 +1282,7 @@ void PARAM::ReadGenotypes(gsl_matrix *UtX, gsl_matrix *K, const bool calc_K) {
     }
   } else {
     if (ReadFile_geno(file_geno, indicator_idv, indicator_snp, UtX, K,
-                      calc_K, mode_debug) == false) {
+                      calc_K) == false) {
       error = true;
     }
   }
@@ -1302,7 +1302,7 @@ void PARAM::ReadGenotypes(vector<vector<unsigned char>> &Xt, gsl_matrix *K,
     }
   } else {
     if (ReadFile_geno(file_geno, indicator_idv, indicator_snp, Xt, K, calc_K,
-                      ni_test, ns_test, mode_debug) == false) {
+                      ni_test, ns_test) == false) {
       error = true;
     }
   }
diff --git a/src/param.h b/src/param.h
index 4b473c0..c4316bb 100644
--- a/src/param.h
+++ b/src/param.h
@@ -26,8 +26,8 @@
 #include <set>
 #include <vector>
 
-#define K_BATCH_SIZE 10000 // #snps used for batched K
-#define DEFAULT_PACE 1000
+#define K_BATCH_SIZE 20000 // #snps used for batched K
+#define DEFAULT_PACE 1000  // for display only
 
 using namespace std;
 
@@ -116,11 +116,11 @@ public:
 class PARAM {
 public:
   // IO-related parameters
-  bool mode_check = true;   // run data checks (slower)
-  bool mode_strict = false; // exit on some data checks
-  bool mode_silence;
-  bool mode_debug = false;
-  uint issue; // enable tests for issue on github tracker
+  // bool mode_check = true;   // run data checks (slower)
+  // bool mode_strict = false; // exit on some data checks
+  // bool mode_silence;
+  // bool mode_debug = false;
+  // uint issue; // enable tests for issue on github tracker
 
   uint a_mode; // Analysis mode, 1/2/3/4 for Frequentist tests
   int k_mode; // Kinship read mode: 1: n by n matrix, 2: id/id/k_value;