Removed WITH_LAPACK compile/build option.

9 files changed, 172 insertions, 564 deletions

diff --git a/src/bslmm.cpp b/src/bslmm.cpp
index 55a05ca..0bd4234 100644
--- a/src/bslmm.cpp
+++ b/src/bslmm.cpp
@@ -1,6 +1,6 @@
 /*
  Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011  Xiang Zhou
+ Copyright (C) 2011-2017 Xiang Zhou
  
  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
@@ -13,7 +13,7 @@
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
- along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
  */
 
 #include <iostream>
@@ -38,30 +38,15 @@
 #include "gsl/gsl_cdf.h"
 #include "gsl/gsl_roots.h"
 
-
-
-
 #include "lapack.h"
-
-#ifdef FORCE_FLOAT
-#include "param_float.h"
-#include "bslmm_float.h"
-#include "lmm_float.h"  //for class FUNC_PARAM and MatrixCalcLR
-#include "lm_float.h"
-#include "mathfunc_float.h"  //for function CenterVector
-#else
 #include "param.h"
 #include "bslmm.h"
 #include "lmm.h"
 #include "lm.h"
 #include "mathfunc.h"
-#endif
 
 using namespace std;
 
-
-
-
 void BSLMM::CopyFromParam (PARAM &cPar) 
 {
 	a_mode=cPar.a_mode;
@@ -330,11 +315,8 @@ double BSLMM::CalcPveLM (const gsl_matrix *UtXgamma, const gsl_vector *Uty, cons
 	gsl_matrix_set_identity (Omega);
 	gsl_matrix_scale (Omega, 1.0/sigma_a2); 
 
-#ifdef WITH_LAPACK
-	lapack_dgemm ((char *)"T", (char *)"N", 1.0, UtXgamma, UtXgamma, 1.0, Omega);
-#else
-	gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, UtXgamma, UtXgamma, 1.0, Omega);	
-#endif
+	lapack_dgemm ((char *)"T", (char *)"N", 1.0, UtXgamma, UtXgamma,
+		      1.0, Omega);
 	gsl_blas_dgemv (CblasTrans, 1.0, UtXgamma, Uty, 0.0, Xty);
 
 	CholeskySolve(Omega, Xty, OiXty);
@@ -517,11 +499,8 @@ double BSLMM::CalcPosterior (const gsl_matrix *UtXgamma, const gsl_vector *Uty,
 	gsl_matrix_set_identity (Omega);
 	
 	time_start=clock();
-#ifdef WITH_LAPACK
-	lapack_dgemm ((char *)"T", (char *)"N", sigma_a2, UtXgamma_eval, UtXgamma, 1.0, Omega);
-#else
-	gsl_blas_dgemm (CblasTrans, CblasNoTrans, sigma_a2, UtXgamma_eval, UtXgamma, 1.0, Omega);
-#endif	
+	lapack_dgemm ((char *)"T", (char *)"N", sigma_a2, UtXgamma_eval,
+		      UtXgamma, 1.0, Omega);
 	time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
 	
 	
@@ -1229,11 +1208,8 @@ void BSLMM::CalcXtX (const gsl_matrix *X, const gsl_vector *y, const size_t s_si
   gsl_matrix_view XtX_sub=gsl_matrix_submatrix(XtX, 0, 0, s_size, s_size);
   gsl_vector_view Xty_sub=gsl_vector_subvector(Xty, 0, s_size);
 
-#ifdef WITH_LAPACK
-  lapack_dgemm ((char *)"T", (char *)"N", 1.0, &X_sub.matrix, &X_sub.matrix, 0.0, &XtX_sub.matrix);
-#else
-  gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, &X_sub.matrix, &X_sub.matrix, 0.0, &XtX_sub.matrix);
-#endif
+  lapack_dgemm ((char *)"T", (char *)"N", 1.0, &X_sub.matrix,
+		&X_sub.matrix, 0.0, &XtX_sub.matrix);
   gsl_blas_dgemv(CblasTrans, 1.0, &X_sub.matrix, y, 0.0, &Xty_sub.vector);
 
   time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
@@ -1282,127 +1258,6 @@ void BSLMM::SetXgamma (const gsl_matrix *X, const gsl_matrix *X_old, const gsl_m
   gsl_vector_view Xtynew_sub=gsl_vector_subvector(Xty_new, 0, rank_new.size());
 
   //get X_new and calculate XtX_new
-  /*
-  if (rank_remove.size()==0 && rank_add.size()==0) {
-    gsl_matrix_memcpy(&Xnew_sub.matrix, &Xold_sub.matrix);
-    gsl_matrix_memcpy(&XtXnew_sub.matrix, &XtXold_sub.matrix);
-    gsl_vector_memcpy(&Xtynew_sub.vector, &Xtyold_sub.vector);
-  } else {
-    gsl_matrix *X_temp=gsl_matrix_alloc(X_old->size1, rank_old.size()-rank_remove.size() );
-    gsl_matrix *XtX_temp=gsl_matrix_alloc(X_temp->size2, X_temp->size2);
-    gsl_vector *Xty_temp=gsl_vector_alloc(X_temp->size2);
-    
-    if (rank_remove.size()==0) {
-      gsl_matrix_memcpy (X_temp, &Xold_sub.matrix);
-      gsl_matrix_memcpy (XtX_temp, &XtXold_sub.matrix);
-      gsl_vector_memcpy (Xty_temp, &Xtyold_sub.vector);
-    } else {
-      size_t i_temp=0, j_temp;
-      for (size_t i=0; i<rank_old.size(); i++) {
-	if (mapRank2in_remove.count(rank_old[i])!=0) {continue;}
-	gsl_vector_const_view Xold_col=gsl_matrix_const_column(X_old, i);	
-	gsl_vector_view Xtemp_col=gsl_matrix_column(X_temp, i_temp);
-	gsl_vector_memcpy (&Xtemp_col.vector, &Xold_col.vector);
-
-	d=gsl_vector_get (Xty_old, i);
-	gsl_vector_set (Xty_temp, i_temp, d);
-	
-	j_temp=i_temp;
-	for (size_t j=i; j<rank_old.size(); j++) {
-	  if (mapRank2in_remove.count(rank_old[j])!=0) {continue;}
-	  d=gsl_matrix_get (XtX_old, i, j);
-	  gsl_matrix_set (XtX_temp, i_temp, j_temp, d);
-	  if (i_temp!=j_temp) {gsl_matrix_set (XtX_temp, j_temp, i_temp, d);}
-	  j_temp++;
-	}
-	i_temp++;
-      }
-    }
-
-    if (rank_add.size()==0) {
-      gsl_matrix_memcpy (&Xnew_sub.matrix, X_temp);
-      gsl_matrix_memcpy (&XtXnew_sub.matrix, XtX_temp);
-      gsl_vector_memcpy (&Xtynew_sub.vector, Xty_temp);
-    } else {
-      gsl_matrix *X_add=gsl_matrix_alloc(X_old->size1, rank_add.size() );
-      gsl_matrix *XtX_aa=gsl_matrix_alloc(X_add->size2, X_add->size2);
-      gsl_matrix *XtX_at=gsl_matrix_alloc(X_add->size2, X_temp->size2);
-      gsl_vector *Xty_add=gsl_vector_alloc(X_add->size2);
-
-      //get X_add
-      SetXgamma (X_add, X, rank_add);
-
-      //get t(X_add)X_add and t(X_add)X_temp	
-      clock_t time_start=clock();
-      
-      //somehow the lapack_dgemm does not work here
-      //#ifdef WITH_LAPACK
-      //lapack_dgemm ((char *)"T", (char *)"N", 1.0, X_add, X_add, 0.0, XtX_aa);
-      //lapack_dgemm ((char *)"T", (char *)"N", 1.0, X_add, X_temp, 0.0, XtX_at);
-      
-      //#else
-      gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, X_add, X_add, 0.0, XtX_aa);
-      gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, X_add, X_temp, 0.0, XtX_at);
-      //#endif
-      gsl_blas_dgemv(CblasTrans, 1.0, X_add, y, 0.0, Xty_add);
-
-      time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
-
-      //save to X_new, XtX_new and Xty_new
-      size_t i_temp=0, j_temp, i_flag=0, j_flag=0;
-      for (size_t i=0; i<rank_new.size(); i++) {
-	if (mapRank2in_add.count(rank_new[i])!=0) {i_flag=1;} else {i_flag=0;}
-	gsl_vector_view Xnew_col=gsl_matrix_column(X_new, i); 
-	if (i_flag==1) {
-	  gsl_vector_view Xcopy_col=gsl_matrix_column(X_add, i-i_temp);
-	  gsl_vector_memcpy (&Xnew_col.vector, &Xcopy_col.vector);
-	} else {
-	  gsl_vector_view Xcopy_col=gsl_matrix_column(X_temp, i_temp);	  
-	  gsl_vector_memcpy (&Xnew_col.vector, &Xcopy_col.vector);
-	}	
-
-	if (i_flag==1) {
-          d=gsl_vector_get (Xty_add, i-i_temp);
-        } else {
-          d=gsl_vector_get (Xty_temp, i_temp);
-        }
-	gsl_vector_set (Xty_new, i, d);
-
-	j_temp=i_temp;
-	for (size_t j=i; j<rank_new.size(); j++) {
-          if (mapRank2in_add.count(rank_new[j])!=0) {j_flag=1;} else {j_flag=0;}
-
-	  if (i_flag==1 && j_flag==1) {
-            d=gsl_matrix_get(XtX_aa, i-i_temp, j-j_temp);
-	  } else if (i_flag==1) {
-	    d=gsl_matrix_get(XtX_at, i-i_temp, j_temp);
-	  } else if (j_flag==1) {
-	    d=gsl_matrix_get(XtX_at, j-j_temp, i_temp);
-	  } else {
-	    d=gsl_matrix_get(XtX_temp, i_temp, j_temp);
-	  }
-
-	  gsl_matrix_set (XtX_new, i, j, d);
-	  if (i!=j) {gsl_matrix_set (XtX_new, j, i, d);}
-
-	  if (j_flag==0) {j_temp++;}
-        }
-	if (i_flag==0) {i_temp++;}
-      }
-
-      gsl_matrix_free(X_add);
-      gsl_matrix_free(XtX_aa);
-      gsl_matrix_free(XtX_at);
-      gsl_vector_free(Xty_add);
-    }
-
-    gsl_matrix_free(X_temp);
-    gsl_matrix_free(XtX_temp);
-    gsl_vector_free(Xty_temp);
-  }
-  */
-
-
   if (rank_remove.size()==0 && rank_add.size()==0) {
     gsl_matrix_memcpy(&Xnew_sub.matrix, &Xold_sub.matrix);
     gsl_matrix_memcpy(&XtXnew_sub.matrix, &XtXold_sub.matrix);
@@ -1447,14 +1302,10 @@ void BSLMM::SetXgamma (const gsl_matrix *X, const gsl_matrix *X_old, const gsl_m
       clock_t time_start=clock();
       
       //somehow the lapack_dgemm does not work here
-      //#ifdef WITH_LAPACK
-      //lapack_dgemm ((char *)"T", (char *)"N", 1.0, X_add, X_add, 0.0, XtX_aa);
-      //lapack_dgemm ((char *)"T", (char *)"N", 1.0, X_add, X_temp, 0.0, XtX_at);
-      
-      //#else
-      gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, X_add, X_add, 0.0, XtX_aa);
-      gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, X_add, X_old, 0.0, XtX_ao);
-      //#endif
+      gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, X_add, X_add,
+		      0.0, XtX_aa);
+      gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, X_add, X_old,
+		      0.0, XtX_ao);
       gsl_blas_dgemv(CblasTrans, 1.0, X_add, y, 0.0, Xty_add);
 
       time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
diff --git a/src/bslmm.h b/src/bslmm.h
index 8b5edc7..863a22d 100644
--- a/src/bslmm.h
+++ b/src/bslmm.h
@@ -1,6 +1,6 @@
 /*
  Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011  Xiang Zhou
+ Copyright (C) 2011-2017 Xiang Zhou
  
  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
@@ -13,7 +13,7 @@
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
- along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
  */
 
 
@@ -25,24 +25,14 @@
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
 
-#ifdef FORCE_FLOAT
-#include "param_float.h"
-#else
 #include "param.h"
-#endif
-
 
 using namespace std;
 
-
-
-
-
-
 class BSLMM {
 
 public:	
-	// IO related parameters
+	// IO-related parameters.
 	int a_mode;	
 	size_t d_pace;
 	
diff --git a/src/bslmmdap.cpp b/src/bslmmdap.cpp
index f6e9e9c..c66189b 100644
--- a/src/bslmmdap.cpp
+++ b/src/bslmmdap.cpp
@@ -1,6 +1,6 @@
 /*
  Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011  Xiang Zhou
+ Copyright (C) 2011-2017 Xiang Zhou
 
  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
@@ -13,7 +13,7 @@
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
- along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
  */
 
 #include <iostream>
@@ -38,31 +38,17 @@
 #include "gsl/gsl_cdf.h"
 #include "gsl/gsl_roots.h"
 
-
-
 #include "logistic.h"
 #include "lapack.h"
 #include "io.h"
-
-#ifdef FORCE_FLOAT
-#include "param_float.h"
-#include "bslmmdap_float.h"
-#include "lmm_float.h"  //for class FUNC_PARAM and MatrixCalcLR
-#include "lm_float.h"
-#include "mathfunc_float.h"  //for function CenterVector
-#else
 #include "param.h"
 #include "bslmmdap.h"
 #include "lmm.h"
 #include "lm.h"
 #include "mathfunc.h"
-#endif
 
 using namespace std;
 
-
-
-
 void BSLMMDAP::CopyFromParam (PARAM &cPar)
 {
 	file_out=cPar.file_out;
@@ -506,11 +492,8 @@ double BSLMMDAP::CalcMarginal (const gsl_matrix *UtXgamma, const gsl_vector *Uty
 	gsl_matrix_set_identity (Omega);
 
 	time_start=clock();
-#ifdef WITH_LAPACK
-	lapack_dgemm ((char *)"T", (char *)"N", sigma_a2, UtXgamma_eval, UtXgamma, 1.0, Omega);
-#else
-	gsl_blas_dgemm (CblasTrans, CblasNoTrans, sigma_a2, UtXgamma_eval, UtXgamma, 1.0, Omega);
-#endif
+	lapack_dgemm ((char *)"T", (char *)"N", sigma_a2, UtXgamma_eval,
+		      UtXgamma, 1.0, Omega);
 	time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
 
 	//calculate beta_hat
@@ -819,198 +802,3 @@ void BSLMMDAP::DAP_EstimateHyper (const size_t kc, const size_t kd, const vector
 	gsl_vector_free(coef);
 	return;
 }
-
-/*
-//readin the estimated hyper-parameters and perform fine mapping for each region
-void BSLMM::DAP_FineMapping (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_matrix *A, const gsl_vector *Uty, const gsl_vector *K_eval, const gsl_vector *y, gsl_matrix *Hyper, gsl_vector *alpha, gsl_vector *pip) {
-	clock_t time_start;
-
-	//two priority sets: S_1 contains all candidate causal SNPs; S_2 contains the prioritized combintion of them
-	//two marginal probability sets: P_1 contains marginals for S_1; P_2 contains marginals for S_2;
-	set<size_t> S1set, S2set;
-	vector<size_t> S1vec;
-	vector<set<size_t> > S2vec;
-	vector<double> P1, P2;
-
-	//calculate P0 (null) and P1 (for every SNP)
-
-
-
-	//loop through the number of combinations
-	for (size_t s=0; s<p; s++) {
-	  //if (s==0), set up S_1: compute marginal of the null model, then compute P_1, then compute BF_1 and use them to select S_1; compute C_1
-
-
-
-	  //if (s==1), set up S_2: compute pair-wise P_2 and use them to select S_2; compute C_2
-
-	  //otherwise, match each combination of S_2 with each SNP from S_1, select into S_3; and replace S_2 with S_3; compute C_s
-
-
-	  //stop when the stopping critieria are reached (if S_2 is empty; if t; if kappa); add the residual component R
-
-	for (size_t t=0; t<total_step; ++t) {
-		if (t%d_pace==0 || t==total_step-1) {ProgressBar ("Running MCMC ", t, total_step-1, (double)n_accept/(double)(t*n_mh+1));}
-//		if (t>10) {break;}
-
-		if (a_mode==13) {
-			SampleZ (y, z_hat, z);
-			mean_z=CenterVector (z);
-
-			time_start=clock();
-			gsl_blas_dgemv (CblasTrans, 1.0, U, z, 0.0, Utz);
-			time_UtZ+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
-
-			//First proposal
-			if (cHyp_old.n_gamma==0 || cHyp_old.rho==0) {
-				logPost_old=CalcPosterior(Utz, K_eval, Utu_old, alpha_old, cHyp_old);
-				beta_old.clear();
-				for (size_t i=0; i<cHyp_old.n_gamma; ++i) {
-				  beta_old.push_back(0);
-				}
-			}
-			else {
-				gsl_matrix *UtXgamma=gsl_matrix_alloc (ni_test, cHyp_old.n_gamma);
-				gsl_vector *beta=gsl_vector_alloc (cHyp_old.n_gamma);
-				SetXgamma (UtXgamma, UtX, rank_old);
-				logPost_old=CalcPosterior(UtXgamma, Utz, K_eval, UtXb_old, Utu_old, alpha_old, beta, cHyp_old);
-
-				beta_old.clear();
-				for (size_t i=0; i<beta->size; ++i) {
-					beta_old.push_back(gsl_vector_get(beta, i));
-				}
-				gsl_matrix_free (UtXgamma);
-				gsl_vector_free (beta);
-			}
-		}
-
-
-	delete [] p_gamma;
-	beta_g.clear();
-
-	return;
-}
-
-*/
-
-
-
-
-
-
-/*
-//below fits MCMC for rho=1
-void BSLMM::CalcXtX (const gsl_matrix *X, const gsl_vector *y, const size_t s_size, gsl_matrix *XtX, gsl_vector *Xty)
-{
-  time_t time_start=clock();
-  gsl_matrix_const_view X_sub=gsl_matrix_const_submatrix(X, 0, 0, X->size1, s_size);
-  gsl_matrix_view XtX_sub=gsl_matrix_submatrix(XtX, 0, 0, s_size, s_size);
-  gsl_vector_view Xty_sub=gsl_vector_subvector(Xty, 0, s_size);
-
-#ifdef WITH_LAPACK
-  lapack_dgemm ((char *)"T", (char *)"N", 1.0, &X_sub.matrix, &X_sub.matrix, 0.0, &XtX_sub.matrix);
-#else
-  gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, &X_sub.matrix, &X_sub.matrix, 0.0, &XtX_sub.matrix);
-#endif
-  gsl_blas_dgemv(CblasTrans, 1.0, &X_sub.matrix, y, 0.0, &Xty_sub.vector);
-
-  time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
-
-  return;
-}
-
-
-
-double BSLMM::CalcPosterior (const double yty, class HYPBSLMM &cHyp)
-{
-	double logpost=0.0;
-
-	//for quantitative traits, calculate pve and pge
-	//pve and pge for case/control data are calculted in CalcCC_PVEnZ
-	if (a_mode==11) {
-		cHyp.pve=0.0;
-		cHyp.pge=1.0;
-	}
-
-	//calculate likelihood
-	if (a_mode==11) {logpost-=0.5*(double)ni_test*log(yty);}
-	else {logpost-=0.5*yty;}
-
-	logpost+=((double)cHyp.n_gamma-1.0)*cHyp.logp+((double)ns_test-(double)cHyp.n_gamma)*log(1-exp(cHyp.logp));
-
-	return logpost;
-}
-
-
-double BSLMM::CalcPosterior (const gsl_matrix *Xgamma, const gsl_matrix *XtX, const gsl_vector *Xty, const double yty, const size_t s_size, gsl_vector *Xb, gsl_vector *beta, class HYPBSLMM &cHyp)
-{
-	double sigma_a2=cHyp.h/( (1-cHyp.h)*exp(cHyp.logp)*(double)ns_test);
-	double logpost=0.0;
-	double d, P_yy=yty, logdet_O=0.0;
-
-	gsl_matrix_const_view Xgamma_sub=gsl_matrix_const_submatrix (Xgamma, 0, 0, Xgamma->size1, s_size);
-	gsl_matrix_const_view XtX_sub=gsl_matrix_const_submatrix (XtX, 0, 0, s_size, s_size);
-	gsl_vector_const_view Xty_sub=gsl_vector_const_subvector (Xty, 0, s_size);
-
-	gsl_matrix *Omega=gsl_matrix_alloc (s_size, s_size);
-	gsl_matrix *M_temp=gsl_matrix_alloc (s_size, s_size);
-	gsl_vector *beta_hat=gsl_vector_alloc (s_size);
-	gsl_vector *Xty_temp=gsl_vector_alloc (s_size);
-
-	gsl_vector_memcpy (Xty_temp, &Xty_sub.vector);
-
-	//calculate Omega
-	gsl_matrix_memcpy (Omega, &XtX_sub.matrix);
-	gsl_matrix_scale (Omega, sigma_a2);
-	gsl_matrix_set_identity (M_temp);
-	gsl_matrix_add (Omega, M_temp);
-
-	//calculate beta_hat
-	logdet_O=CholeskySolve(Omega, Xty_temp, beta_hat);
-	gsl_vector_scale (beta_hat, sigma_a2);
-
-	gsl_blas_ddot (Xty_temp, beta_hat, &d);
-	P_yy-=d;
-
-	//sample tau
-	double tau=1.0;
-	if (a_mode==11) {tau =gsl_ran_gamma (gsl_r, (double)ni_test/2.0,  2.0/P_yy); }
-
-	//sample beta
-	for (size_t i=0; i<s_size; i++)
-	{
-		d=gsl_ran_gaussian(gsl_r, 1);
-		gsl_vector_set(beta, i, d);
-	}
-	gsl_vector_view beta_sub=gsl_vector_subvector(beta, 0, s_size);
-	gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, &beta_sub.vector);
-
-	//it compuates inv(L^T(Omega)) %*% beta;
-	gsl_vector_scale(&beta_sub.vector, sqrt(sigma_a2/tau));
-	gsl_vector_add(&beta_sub.vector, beta_hat);
-	gsl_blas_dgemv (CblasNoTrans, 1.0, &Xgamma_sub.matrix, &beta_sub.vector, 0.0, Xb);
-
-	//for quantitative traits, calculate pve and pge
-	if (a_mode==11) {
-		gsl_blas_ddot (Xb, Xb, &d);
-		cHyp.pve=d/(double)ni_test;
-		cHyp.pve/=cHyp.pve+1.0/tau;
-		cHyp.pge=1.0;
-	}
-
-	logpost=-0.5*logdet_O;
-	if (a_mode==11) {logpost-=0.5*(double)ni_test*log(P_yy);}
-	else {logpost-=0.5*P_yy;}
-
-	logpost+=((double)cHyp.n_gamma-1.0)*cHyp.logp+((double)ns_test-(double)cHyp.n_gamma)*log(1.0-exp(cHyp.logp));
-
-	gsl_matrix_free (Omega);
-	gsl_matrix_free (M_temp);
-	gsl_vector_free (beta_hat);
-	gsl_vector_free (Xty_temp);
-
-	return logpost;
-}
-*/
-
-
diff --git a/src/bslmmdap.h b/src/bslmmdap.h
index ac78f97..a9e2d51 100644
--- a/src/bslmmdap.h
+++ b/src/bslmmdap.h
@@ -1,6 +1,6 @@
 /*
  Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011  Xiang Zhou
+ Copyright (C) 2011-2017 Xiang Zhou
 
  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
@@ -13,10 +13,9 @@
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
- along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
  */
 
-
 #ifndef __BSLMMDAP_H__
 #define __BSLMMDAP_H__
 
@@ -24,25 +23,14 @@
 #include <map>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
-
-#ifdef FORCE_FLOAT
-#include "param_float.h"
-#else
 #include "param.h"
-#endif
-
 
 using namespace std;
 
-
-
-
-
-
 class BSLMMDAP {
 
 public:
-	// IO related parameters
+	// IO-related parameters.
 	int a_mode;
 	size_t d_pace;
 
diff --git a/src/io.cpp b/src/io.cpp
index 1a4d466..e8f00ed 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -1,6 +1,6 @@
 /*
     Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011-2017 Xiang Zhou
+    Copyright (C) 2011-2017, Xiang Zhou
 
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -194,11 +194,13 @@ bool ReadFile_snps_header (const string &file_snps, set<string> &setSnps) {
 	return true;
 }
 
-//Read log file
-bool ReadFile_log (const string &file_log, double &pheno_mean)
-{
+// Read log file.
+bool ReadFile_log (const string &file_log, double &pheno_mean) {
 	ifstream infile (file_log.c_str(), ifstream::in);
-	if (!infile) {cout<<"error! fail to open log file: "<<file_log<<endl; return false;}
+	if (!infile) {
+	  cout << "error! fail to open log file: " << file_log << endl;
+	  return false;
+	}
 
 	string line;
 	char *ch_ptr;
@@ -229,15 +231,18 @@ bool ReadFile_log (const string &file_log, double &pheno_mean)
 	return true;
 }
 
-
-//Read bimbam annotation file
-bool ReadFile_anno (const string &file_anno, map<string, string> &mapRS2chr, map<string, long int> &mapRS2bp, map<string, double> &mapRS2cM)
-{
+// Read bimbam annotation file.
+bool ReadFile_anno (const string &file_anno, map<string, string> &mapRS2chr,
+		    map<string, long int> &mapRS2bp,
+		    map<string, double> &mapRS2cM) {
 	mapRS2chr.clear();
 	mapRS2bp.clear();
 
 	ifstream infile (file_anno.c_str(), ifstream::in);
-	if (!infile) {cout<<"error opening annotation file: "<<file_anno<<endl; return false;}
+	if (!infile) {
+	  cout << "error opening annotation file: " << file_anno << endl;
+	  return false;
+	}
 
 	string line;
 	char *ch_ptr;
@@ -251,11 +256,23 @@ bool ReadFile_anno (const string &file_anno, map<string, string> &mapRS2chr, map
 		ch_ptr=strtok ((char *)line.c_str(), " , \t");
 		rs=ch_ptr;
 		ch_ptr=strtok (NULL, " , \t");
-		if (strcmp(ch_ptr, "NA")==0) {b_pos=-9;} else {b_pos=atol(ch_ptr);}
+		if (strcmp(ch_ptr, "NA")==0) {
+		  b_pos=-9;
+		} else {
+		  b_pos=atol(ch_ptr);
+		}
 		ch_ptr=strtok (NULL, " , \t");
-		if (ch_ptr==NULL || strcmp(ch_ptr, "NA")==0) {chr="-9";} else {chr=ch_ptr;}
+		if (ch_ptr==NULL || strcmp(ch_ptr, "NA")==0) {
+		  chr="-9";
+		} else {
+		  chr=ch_ptr;
+		}
 		ch_ptr=strtok (NULL, " , \t");
-		if (ch_ptr==NULL || strcmp(ch_ptr, "NA")==0) {cM=-9;} else {cM=atof(ch_ptr);}
+		if (ch_ptr==NULL || strcmp(ch_ptr, "NA")==0) {
+		  cM=-9;
+		} else {
+		  cM=atof(ch_ptr);
+		}
 
 		mapRS2chr[rs]=chr;
 		mapRS2bp[rs]=b_pos;
@@ -268,15 +285,17 @@ bool ReadFile_anno (const string &file_anno, map<string, string> &mapRS2chr, map
 	return true;
 }
 
-//read one column of phenotype
-bool ReadFile_column (const string &file_pheno, vector<int> &indicator_idv, vector<double> &pheno, const int &p_column)
-{
+// Read 1 column of phenotype.
+bool ReadFile_column (const string &file_pheno, vector<int> &indicator_idv,
+		      vector<double> &pheno, const int &p_column) {
 	indicator_idv.clear();
 	pheno.clear();
 
 	igzstream infile (file_pheno.c_str(), igzstream::in);
-//	ifstream infile (file_pheno.c_str(), ifstream::in);
-	if (!infile) {cout<<"error! fail to open phenotype file: "<<file_pheno<<endl; return false;}
+	if (!infile) {
+	  cout << "error! fail to open phenotype file: " << file_pheno << endl;
+	  return false;
+	}
 
 	string line;
 	char *ch_ptr;
@@ -288,8 +307,15 @@ bool ReadFile_column (const string &file_pheno, vector<int> &indicator_idv, vect
 		for (int i=0; i<(p_column-1); ++i) {
 			ch_ptr=strtok (NULL, " , \t");
 		}
-		if (strcmp(ch_ptr, "NA")==0) {indicator_idv.push_back(0); pheno.push_back(-9);}		//pheno is different from pimass2
-		else {p=atof(ch_ptr); indicator_idv.push_back(1); pheno.push_back(p);}
+		if (strcmp(ch_ptr, "NA")==0) {
+		  indicator_idv.push_back(0);
+		  pheno.push_back(-9);
+		} // Pheno is different from pimass2.
+		else {
+		  p=atof(ch_ptr);
+		  indicator_idv.push_back(1);
+		  pheno.push_back(p);
+		}
 	}
 
 	infile.close();
@@ -298,17 +324,19 @@ bool ReadFile_column (const string &file_pheno, vector<int> &indicator_idv, vect
 	return true;
 }
 
-
-
-//Read bimbam phenotype file, p_column=1, 2 ...
-bool ReadFile_pheno (const string &file_pheno, vector<vector<int> > &indicator_pheno, vector<vector<double> > &pheno, const vector<size_t> &p_column)
-{
+// Read bimbam phenotype file, p_column=1, 2,...
+bool ReadFile_pheno (const string &file_pheno,
+		     vector<vector<int> > &indicator_pheno,
+		     vector<vector<double> > &pheno,
+		     const vector<size_t> &p_column) {
 	indicator_pheno.clear();
 	pheno.clear();
 
 	igzstream infile (file_pheno.c_str(), igzstream::in);
-//	ifstream infile (file_pheno.c_str(), ifstream::in);
-	if (!infile) {cout<<"error! fail to open phenotype file: "<<file_pheno<<endl; return false;}
+	if (!infile) {
+	  cout << "error! fail to open phenotype file: " << file_pheno << endl;
+	  return false;
+	}
 
 	string line;
 	char *ch_ptr;
@@ -333,8 +361,15 @@ bool ReadFile_pheno (const string &file_pheno, vector<vector<int> > &indicator_p
 		size_t i=0;
 		while (i<p_max ) {
 			if (mapP2c.count(i+1)!=0) {
-				if (strcmp(ch_ptr, "NA")==0) {ind_pheno_row[mapP2c[i+1]]=0; pheno_row[mapP2c[i+1]]=-9;}
-				else {p=atof(ch_ptr); ind_pheno_row[mapP2c[i+1]]=1; pheno_row[mapP2c[i+1]]=p;}
+				if (strcmp(ch_ptr, "NA")==0) {
+				  ind_pheno_row[mapP2c[i+1]]=0;
+				  pheno_row[mapP2c[i+1]]=-9;
+				}
+				else {
+				  p=atof(ch_ptr);
+				  ind_pheno_row[mapP2c[i+1]]=1;
+				  pheno_row[mapP2c[i+1]]=p;
+				}
 			}
 			i++;
 			ch_ptr=strtok (NULL, " , \t");
@@ -350,13 +385,15 @@ bool ReadFile_pheno (const string &file_pheno, vector<vector<int> > &indicator_p
 	return true;
 }
 
-
-bool ReadFile_cvt (const string &file_cvt, vector<int> &indicator_cvt, vector<vector<double> > &cvt, size_t &n_cvt)
-{
+bool ReadFile_cvt (const string &file_cvt, vector<int> &indicator_cvt,
+		   vector<vector<double> > &cvt, size_t &n_cvt) {
 	indicator_cvt.clear();
 
 	ifstream infile (file_cvt.c_str(), ifstream::in);
-	if (!infile) {cout<<"error! fail to open covariates file: "<<file_cvt<<endl; return false;}
+	if (!infile) {
+	  cout << "error! fail to open covariates file: " << file_cvt << endl;
+	  return false;
+	}
 
 	string line;
 	char *ch_ptr;
@@ -374,7 +411,11 @@ bool ReadFile_cvt (const string &file_cvt, vector<int> &indicator_cvt, vector<ve
 			v_d.push_back(d);
 			ch_ptr=strtok (NULL, " , \t");
 		}
-		if (flag_na==0) {indicator_cvt.push_back(1);} else {indicator_cvt.push_back(0);}
+		if (flag_na==0) {
+		  indicator_cvt.push_back(1);
+		} else {
+		  indicator_cvt.push_back(0);
+		}
 		cvt.push_back(v_d);
 	}
 
@@ -382,10 +423,16 @@ bool ReadFile_cvt (const string &file_cvt, vector<int> &indicator_cvt, vector<ve
 	else {
 		flag_na=0;
 		for (vector<int>::size_type i=0; i<indicator_cvt.size(); ++i) {
-			if (indicator_cvt[i]==0) {continue;}
+			if (indicator_cvt[i]==0) {
+			  continue;
+			}
 
 			if (flag_na==0) {flag_na=1; n_cvt=cvt[i].size();}
-			if (flag_na!=0 && n_cvt!=cvt[i].size()) {cout<<"error! number of covariates in row "<<i<<" do not match other rows."<<endl; return false;}
+			if (flag_na!=0 && n_cvt!=cvt[i].size()) {
+			  cout << "error! number of covariates in row " <<
+			    i << " do not match other rows." << endl;
+			  return false;
+			}
 		}
 	}
 
@@ -395,15 +442,15 @@ bool ReadFile_cvt (const string &file_cvt, vector<int> &indicator_cvt, vector<ve
 	return true;
 }
 
-
-
-//Read .bim file
-bool ReadFile_bim (const string &file_bim, vector<SNPINFO> &snpInfo)
-{
-  snpInfo.clear();
+// Read .bim file.
+bool ReadFile_bim (const string &file_bim, vector<SNPINFO> &snpInfo) {
+        snpInfo.clear();
 
 	ifstream infile (file_bim.c_str(), ifstream::in);
-	if (!infile) {cout<<"error opening .bim file: "<<file_bim<<endl; return false;}
+	if (!infile) {
+	  cout << "error opening .bim file: " << file_bim << endl;
+	  return false;
+	}
 
 	string line;
 	char *ch_ptr;
@@ -429,7 +476,8 @@ bool ReadFile_bim (const string &file_bim, vector<SNPINFO> &snpInfo)
 		ch_ptr=strtok (NULL, " \t");
 		major=ch_ptr;
 
-		SNPINFO sInfo={chr, rs, cM, b_pos, minor, major, 0, -9, -9, 0, 0, 0};
+		SNPINFO sInfo={chr, rs, cM, b_pos, minor, major,
+			       0, -9, -9, 0, 0, 0};
 		snpInfo.push_back(sInfo);
 	}
 
@@ -438,8 +486,7 @@ bool ReadFile_bim (const string &file_bim, vector<SNPINFO> &snpInfo)
 	return true;
 }
 
-
-//Read .fam file
+// Read .fam file.
 bool ReadFile_fam (const string &file_fam, vector<vector<int> > &indicator_pheno, vector<vector<double> > &pheno, map<string, int> &mapID2num, const vector<size_t> &p_column)
 {
 	indicator_pheno.clear();
diff --git a/src/io.h b/src/io.h
index 365127a..9d6f8cc 100644
--- a/src/io.h
+++ b/src/io.h
@@ -1,6 +1,6 @@
 /*
     Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011-2017 Xiang Zhou
+    Copyright (C) 2011-2017, Xiang Zhou
 
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -121,22 +121,60 @@ bool ReadFile_gene (const string &file_gene, vector<double> &vec_read,
 		    vector<SNPINFO> &snpInfo, size_t &ng_total);
 
 bool ReadHeader_io (const string &line, HEADER &header);
-bool ReadFile_cat (const string &file_cat, map<string, size_t> &mapRS2cat, size_t &n_vc);
-bool ReadFile_mcat (const string &file_mcat, map<string, size_t> &mapRS2cat, size_t &n_vc);
-
-bool ReadFile_catc (const string &file_cat, map<string, vector<double> > &mapRS2catc, size_t &n_cat);
-bool ReadFile_mcatc (const string &file_mcat, map<string, vector<double> > &mapRS2catc, size_t &n_cat);
-
-bool BimbamKin (const string &file_geno, const int display_pace, const vector<int> &indicator_idv, const vector<int> &indicator_snp, const map<string, double> &mapRS2weight, const map<string, size_t> &mapRS2cat, const vector<SNPINFO> &snpInfo, const gsl_matrix *W, gsl_matrix *matrix_kin, gsl_vector *vector_ns);
-bool PlinkKin (const string &file_bed, const int display_pace, const vector<int> &indicator_idv, const vector<int> &indicator_snp, const map<string, double> &mapRS2weight, const map<string, size_t> &mapRS2cat, const vector<SNPINFO> &snpInfo, const gsl_matrix *W, gsl_matrix *matrix_kin, gsl_vector *vector_ns);
-bool MFILEKin (const size_t mfile_mode, const string &file_mfile, const int display_pace, const vector<int> &indicator_idv, const vector<vector<int> > &mindicator_snp, const map<string, double> &mapRS2weight, const map<string, size_t> &mapRS2cat, const vector<vector<SNPINFO> > &msnpInfo, const gsl_matrix *W, gsl_matrix *matrix_kin, gsl_vector *vector_ns);
-
-bool ReadFile_wsnp (const string &file_wsnp, map<string, double> &mapRS2double);
-bool ReadFile_wsnp (const string &file_wcat, const size_t n_vc, map<string, vector<double> > &mapRS2vector);
-
-void ReadFile_beta (const string &file_beta, const map<string, size_t> &mapRS2cat, const map<string, double> &mapRS2wA, vector<size_t> &vec_cat, vector<size_t> &vec_ni, vector<double> &vec_weight, vector<double> &vec_z2, size_t &ni_total, size_t &ns_total, size_t &ns_test);
-void ReadFile_beta (const string &file_beta, const map<string, double> &mapRS2wA, map<string, string> &mapRS2A1, map<string, double> &mapRS2z);
-void Calcq (const size_t n_block, const vector<size_t> &vec_cat, const vector<size_t> &vec_ni, const vector<double> &vec_weight, const vector<double> &vec_z2, gsl_matrix *Vq, gsl_vector *q, gsl_vector *s);
+bool ReadFile_cat (const string &file_cat, map<string, size_t> &mapRS2cat,
+		   size_t &n_vc);
+bool ReadFile_mcat (const string &file_mcat, map<string, size_t> &mapRS2cat,
+		    size_t &n_vc);
+
+bool ReadFile_catc (const string &file_cat,
+		    map<string, vector<double> > &mapRS2catc,
+		    size_t &n_cat);
+bool ReadFile_mcatc (const string &file_mcat, map<string,
+		     vector<double> > &mapRS2catc, size_t &n_cat);
+
+bool BimbamKin (const string &file_geno, const int display_pace,
+		const vector<int> &indicator_idv,
+		const vector<int> &indicator_snp,
+		const map<string, double> &mapRS2weight,
+		const map<string, size_t> &mapRS2cat,
+		const vector<SNPINFO> &snpInfo, const gsl_matrix *W,
+		gsl_matrix *matrix_kin, gsl_vector *vector_ns);
+bool PlinkKin (const string &file_bed, const int display_pace,
+	       const vector<int> &indicator_idv,
+	       const vector<int> &indicator_snp,
+	       const map<string, double> &mapRS2weight,
+	       const map<string, size_t> &mapRS2cat,
+	       const vector<SNPINFO> &snpInfo,
+	       const gsl_matrix *W, gsl_matrix *matrix_kin,
+	       gsl_vector *vector_ns);
+bool MFILEKin (const size_t mfile_mode, const string &file_mfile,
+	       const int display_pace, const vector<int> &indicator_idv,
+	       const vector<vector<int> > &mindicator_snp,
+	       const map<string, double> &mapRS2weight,
+	       const map<string, size_t> &mapRS2cat,
+	       const vector<vector<SNPINFO> > &msnpInfo,
+	       const gsl_matrix *W, gsl_matrix *matrix_kin,
+	       gsl_vector *vector_ns);
+
+bool ReadFile_wsnp (const string &file_wsnp,
+		    map<string, double> &mapRS2double);
+bool ReadFile_wsnp (const string &file_wcat, const size_t n_vc,
+		    map<string, vector<double> > &mapRS2vector);
+
+void ReadFile_beta (const string &file_beta,
+		    const map<string, size_t> &mapRS2cat,
+		    const map<string, double> &mapRS2wA,
+		    vector<size_t> &vec_cat, vector<size_t> &vec_ni,
+		    vector<double> &vec_weight, vector<double> &vec_z2,
+		    size_t &ni_total, size_t &ns_total, size_t &ns_test);
+void ReadFile_beta (const string &file_beta,
+		    const map<string, double> &mapRS2wA,
+		    map<string, string> &mapRS2A1,
+		    map<string, double> &mapRS2z);
+void Calcq (const size_t n_block, const vector<size_t> &vec_cat,
+	    const vector<size_t> &vec_ni,
+	    const vector<double> &vec_weight, const vector<double> &vec_z2,
+	    gsl_matrix *Vq, gsl_vector *q, gsl_vector *s);
 
 void ReadFile_study (const string &file_study, gsl_matrix *Vq,
 		     gsl_vector *q_vec, gsl_vector *s_vec, size_t &ni);
diff --git a/src/lapack.cpp b/src/lapack.cpp
index 18399a2..2bbdf62 100644
--- a/src/lapack.cpp
+++ b/src/lapack.cpp
@@ -398,13 +398,7 @@ void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
 // DO NOT set eigenvalues to be positive.
 double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
 		    const size_t flag_largematrix) {
-#ifdef WITH_LAPACK
 	lapack_eigen_symmv (G, eval, U, flag_largematrix);
-#else
-	gsl_eigen_symmv_workspace *w=gsl_eigen_symmv_alloc (G->size1);
-	gsl_eigen_symmv (G, eval, U, w);
-	gsl_eigen_symmv_free (w);	
-#endif	
 
 	// Calculate track_G=mean(diag(G)).
 	double d=0.0;
@@ -420,48 +414,7 @@ double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
 // DO NOT set eigen values to be positive.
 double EigenDecomp (gsl_matrix_float *G, gsl_matrix_float *U,
 		    gsl_vector_float *eval, const size_t flag_largematrix) {
-#ifdef WITH_LAPACK
 	lapack_float_eigen_symmv (G, eval, U, flag_largematrix);
-#else
-	
-	// GSL doesn't provide single precision eigen decomposition;
-	// plus, float precision eigen decomposition in LAPACK may not
-	// work on OS 10.4 first change to double precision.
-	gsl_matrix *G_double=gsl_matrix_alloc (G->size1, G->size2);
-	gsl_matrix *U_double=gsl_matrix_alloc (U->size1, U->size2);
-	gsl_vector *eval_double=gsl_vector_alloc (eval->size);
-	for (size_t i=0; i<G->size1; i++) {
-		for (size_t j=0; j<G->size2; j++) {
-			gsl_matrix_set(G_double, i, j,
-				       gsl_matrix_float_get(G, i, j));
-		}
-	}	
-	gsl_eigen_symmv_workspace *w_space=gsl_eigen_symmv_alloc (G->size1);
-	gsl_eigen_symmv (G_double, eval_double, U_double, w_space);
-	gsl_eigen_symmv_free (w_space);	
-	
-	// Change back to float precision.
-	for (size_t i=0; i<G->size1; i++) {
-		for (size_t j=0; j<G->size2; j++) {
-			gsl_matrix_float_set(K, i, j,
-					     gsl_matrix_get(G_double, i, j));
-		}
-	}
-	for (size_t i=0; i<U->size1; i++) {
-		for (size_t j=0; j<U->size2; j++) {
-			gsl_matrix_float_set(U, i, j,
-					     gsl_matrix_get(U_double, i, j));
-		}
-	}
-	for (size_t i=0; i<eval->size; i++) {
-		gsl_vector_float_set(eval, i, gsl_vector_get(eval_double, i));
-	}	
-	
-	// Delete double-precision matrices.
-	gsl_matrix_free (G_double);
-	gsl_matrix_free (U_double);
-	gsl_vector_free (eval_double);
-#endif
 
 	// Calculate track_G=mean(diag(G)).
 	double d = 0.0;
@@ -477,28 +430,12 @@ double EigenDecomp (gsl_matrix_float *G, gsl_matrix_float *U,
 double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty) {
 	double logdet_O=0.0;
 	
-#ifdef WITH_LAPACK
 	lapack_cholesky_decomp(Omega);
 	for (size_t i=0; i<Omega->size1; ++i) {
 		logdet_O+=log(gsl_matrix_get (Omega, i, i));
 	}	
 	logdet_O*=2.0;	
 	lapack_cholesky_solve(Omega, Xty, OiXty);	
-#else	
-	int status = gsl_linalg_cholesky_decomp(Omega);
-	if(status == GSL_EDOM) {
-		cout << "## non-positive definite matrix" << endl; 
-	}
-	
-	for (size_t i=0; i<Omega->size1; ++i) {
-		logdet_O+=log(gsl_matrix_get (Omega, i, i));
-	}
-	logdet_O*=2.0;	
-	
-	gsl_vector_memcpy (OiXty, Xty);
-	gsl_blas_dtrsv(CblasLower, CblasNoTrans, CblasNonUnit, Omega, OiXty); 
-	gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, OiXty);
-#endif
 	
 	return logdet_O;
 }
@@ -508,43 +445,12 @@ double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty,
 		     gsl_vector_float *OiXty) {
 	double logdet_O=0.0;
 	
-#ifdef WITH_LAPACK
 	lapack_float_cholesky_decomp(Omega);
 	for (size_t i=0; i<Omega->size1; ++i) {
 		logdet_O+=log(gsl_matrix_float_get (Omega, i, i));
 	}	
 	logdet_O*=2.0;	
 	lapack_float_cholesky_solve(Omega, Xty, OiXty);	
-#else
-	gsl_matrix *Omega_double=gsl_matrix_alloc (Omega->size1, Omega->size2);
-	double d;
-	for (size_t i=0; i<Omega->size1; ++i) {
-		for (size_t j=0; j<Omega->size2; ++j) {
-			d=(double)gsl_matrix_float_get (Omega, i, j);
-			gsl_matrix_set (Omega_double, i, j, d);
-		}
-	}
-	
-	int status = gsl_linalg_cholesky_decomp(Omega_double);
-	if(status == GSL_EDOM) {
-		cout << "## non-positive definite matrix" << endl; 
-	}	
-	
-	for (size_t i=0; i<Omega->size1; ++i) {
-		for (size_t j=0; j<Omega->size2; ++j) {
-			d=gsl_matrix_get (Omega_double, i, j);
-			if (j==i) {logdet_O+=log(d);}
-			gsl_matrix_float_set (Omega, i, j, (float)d);
-		}
-	}
-	logdet_O*=2.0;	
-	
-	gsl_vector_float_memcpy (OiXty, Xty);
-	gsl_blas_strsv(CblasLower, CblasNoTrans, CblasNonUnit, Omega, OiXty); 
-	gsl_blas_strsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, OiXty);
-	
-	gsl_matrix_free (Omega_double);
-#endif
 	
 	return logdet_O;
 }	
diff --git a/src/param.cpp b/src/param.cpp
index 63bf349..0a1d721 100644
--- a/src/param.cpp
+++ b/src/param.cpp
@@ -1,6 +1,6 @@
 /*
     Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011-2017 Xiang Zhou
+    Copyright (C) 2011-2017, Xiang Zhou
 
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
diff --git a/src/param.h b/src/param.h
index 18d5e36..9707790 100644
--- a/src/param.h
+++ b/src/param.h
@@ -1,6 +1,6 @@
 /*
     Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011-2017 Xiang Zhou
+    Copyright (C) 2011-2017, Xiang Zhou
 
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by