Merge branch 'master' of github.com:genenetwork/GEMMA

1 files changed, 41 insertions, 4 deletions

diff --git a/src/fastblas.cpp b/src/fastblas.cpp
index b3fcddb..42c59ee 100644
--- a/src/fastblas.cpp
+++ b/src/fastblas.cpp
@@ -18,7 +18,7 @@
     along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include "gsl/gsl_matrix.h"
+// #include "gsl/gsl_matrix.h"
 #include <algorithm>    // std::min
 #include <cmath>
 #include <iomanip>
@@ -29,6 +29,7 @@
 #include "mathfunc.h"
 #include <string.h>
 #include "eigenlib.h"
+#include <cblas.h>
 
 const char *FastblasTrans = "T";
 const char *FastblasNoTrans = "N";
@@ -235,8 +236,44 @@ void fast_dgemm(const char *TransA, const char *TransB, const double alpha,
 void fast_eigen_dgemm(const char *TransA, const char *TransB, const double alpha,
                       const gsl_matrix *A, const gsl_matrix *B, const double beta,
                       gsl_matrix *C) {
-  if (is_legacy_mode())
-    eigenlib_dgemm(TransA,TransB,alpha,A,B,beta,C);
-  else
     fast_cblas_dgemm(TransA,TransB,alpha,A,B,beta,C);
 }
+
+/*
+ *  Inverse in place
+ */
+
+#include <gsl/gsl_permutation.h>
+// #include <gsl/gsl_linalg.h>
+
+extern "C" {
+  int gsl_linalg_LU_invert(const gsl_matrix * LU, const gsl_permutation * p, gsl_matrix * inverse);
+  int gsl_linalg_LU_decomp(gsl_matrix * A, gsl_permutation * p, int * signum);
+}
+
+void gsl_matrix_inv(gsl_matrix *m)  
+{  
+    size_t n=m->size1;  
+  
+    gsl_matrix *temp1=gsl_matrix_calloc(n,n);  
+    gsl_matrix_memcpy(temp1,m);  
+  
+    gsl_permutation *p=gsl_permutation_calloc(n);  
+    int sign=0;  
+    gsl_linalg_LU_decomp(temp1,p,&sign);  
+    gsl_matrix *inverse=gsl_matrix_calloc(n,n);  
+  
+    gsl_linalg_LU_invert(temp1,p,inverse);  
+    gsl_matrix_memcpy(m,inverse);  
+  
+    gsl_permutation_free(p);  
+    gsl_matrix_free(temp1);  
+    gsl_matrix_free(inverse);  
+  
+}
+
+void fast_inverse(gsl_matrix *m) {
+    gsl_matrix_inv(m);
+}
+
+