1 files changed, 31 insertions, 24 deletions

diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index 4203837..e7dff73 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -32,7 +32,7 @@
 #include <tuple>
 #include <vector>
 
-#include "Eigen/Dense"
+// #include "Eigen/Dense"
 
 #include "gsl/gsl_version.h"
 
@@ -49,11 +49,12 @@
 
 #include "debug.h"
 #include "eigenlib.h"
+#include "fastblas.h"
 #include "lapack.h"
 #include "mathfunc.h"
 
 using namespace std;
-using namespace Eigen;
+// using namespace Eigen;
 
 bool has_nan(const vector<double> v) {
   for (const auto& e: v) {
@@ -79,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);
@@ -104,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);
@@ -126,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);
 
@@ -226,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
@@ -241,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);
@@ -313,6 +314,13 @@ bool isMatrixIllConditioned(const gsl_vector *eigenvalues, double max_ratio) {
   return ret_valid;
 }
 
+double sum(const double *m, size_t rows, size_t cols) {
+  double sum = 0.0;
+  for (auto i = 0; i<rows*cols; i++)
+    sum += m[i];
+  return sum;
+}
+
 double SumVector(const gsl_vector *v) {
   double sum = 0;
   for (int i = 0; i < v->size; i++ ) {
@@ -337,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);
@@ -379,22 +387,18 @@ 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);
-  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.
@@ -520,6 +524,7 @@ unsigned char Double02ToUchar(const double dosage) {
   return (int)(dosage * 100);
 }
 
+/*
 void uchar_matrix_get_row(const vector<vector<unsigned char>> &X,
                           const size_t i_row, VectorXd &x_row) {
   if (i_row < X.size()) {
@@ -531,3 +536,5 @@ void uchar_matrix_get_row(const vector<vector<unsigned char>> &X,
     exit(1);
   }
 }
+
+*/