Toying with lmdbi

1 files changed, 14 insertions, 7 deletions

diff --git a/src/gemma.cpp b/src/gemma.cpp
index fdb19fd..6b42763 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -1914,8 +1914,14 @@ void GEMMA::BatchRun(PARAM &cPar) {
   if (cPar.a_mode == M_KIN || cPar.a_mode == M_KIN2) {
     cout << "Calculating Relatedness Matrix ... " << endl;
 
-    gsl_matrix *G = gsl_matrix_safe_alloc(cPar.ni_total, cPar.ni_total);
-    enforce_msg(G, "allocate G"); // just to be sure
+    if (cPar.is_mdb) {
+      cPar.MdbCalcKin();
+      return;
+    }
+
+    enforce(cPar.ni_total > 0);
+    gsl_matrix *K = gsl_matrix_safe_alloc(cPar.ni_total, cPar.ni_total);
+    enforce_msg(K, "allocate G"); // just to be sure
 
     time_start = clock();
 
@@ -1923,7 +1929,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
       cout << "error! failed to prepare for calculating relatedness matrix. " << endl;
       return;
     }
-    cPar.CalcKin(G);
+    cPar.CalcKin(K);
 
     cPar.time_G = (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0);
     if (cPar.error == true) {
@@ -1932,15 +1938,15 @@ void GEMMA::BatchRun(PARAM &cPar) {
     }
 
     // Now we have the Kinship matrix test it
-    validate_K(G);
+    validate_K(K);
 
     if (cPar.a_mode == M_KIN) {
-      cPar.WriteMatrix(G, "cXX");
+      cPar.WriteMatrix(K, "cXX");
     } else { // M_KIN2
-      cPar.WriteMatrix(G, "sXX");
+      cPar.WriteMatrix(K, "sXX");
     }
 
-    gsl_matrix_safe_free(G);
+    gsl_matrix_safe_free(K);
   }
 
   // Compute the LDSC weights (not implemented yet)
@@ -2582,6 +2588,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     gsl_matrix *B = gsl_matrix_safe_alloc(Y->size2, W->size2);      // B is a d by c (effect size)
                                                           // matrix
     gsl_matrix *se_B = gsl_matrix_safe_alloc(Y->size2, W->size2);   // se
+    enforce(Y->size1 > 0);
     gsl_matrix *K = gsl_matrix_safe_alloc(Y->size1, Y->size1);      // Kinship aka as G/GRM in the equation ind x ind
     gsl_matrix *U = gsl_matrix_safe_alloc(Y->size1, Y->size1);
     gsl_matrix *UtW = gsl_matrix_calloc(Y->size1, W->size2);        // Transformed ind x covariates