about summary refs log tree commit diff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile8
-rw-r--r--src/eigenlib.cpp107
-rw-r--r--src/eigenlib.h35
-rw-r--r--src/fastblas.cpp39
-rw-r--r--src/fastblas.h2
-rw-r--r--src/fastopenblas.h4
-rw-r--r--src/gemma.cpp22
-rw-r--r--src/gemma_io.cpp2
-rw-r--r--src/lmm.cpp4
-rw-r--r--src/param.cpp3
-rw-r--r--src/vc.cpp22
-rw-r--r--test/src/unittests-math.cpp8
12 files changed, 74 insertions, 182 deletions
diff --git a/Makefile b/Makefile
index 1376443..191285a 100644
--- a/Makefile
+++ b/Makefile
@@ -93,7 +93,7 @@ else
   ifdef GUIX
     # Effectively disable paths for GNU Guix
     OPENBLAS_INCLUDE_PATH = .
-    EIGEN_INCLUDE_PATH = $(GUIX)/include/eigen3
+    # EIGEN_INCLUDE_PATH = $(GUIX)/include/eigen3
     # RPATH = -Xlinker --rpath=$(GUIX)/lib
     ifdef FORCE_STATIC
       LIBS = -L$(GUIX)/lib
@@ -139,11 +139,11 @@ else
 endif
 
 ifneq ($(CXX), clang++)
-  # Clang does not like this switch
-  debug check fast-check: CPPFLAGS += -Og
+  # Clang does not like these switches
+  debug check fast-check: CPPFLAGS += -Og -Wfatal-errors
 endif
 
-debug check fast-check: CPPFLAGS += -g $(GCC_FLAGS) $(GSL_INCLUDE_PATH) -isystem$(EIGEN_INCLUDE_PATH) -Icontrib/catch-1.9.7 -Isrc
+debug check fast-check: CPPFLAGS += -g $(GCC_FLAGS) $(GSL_INCLUDE_PATH) -Icontrib/catch-1.9.7 -Isrc
 
 profile: CPPFLAGS += -pg
 
diff --git a/src/eigenlib.cpp b/src/eigenlib.cpp
deleted file mode 100644
index 470aa08..0000000
--- a/src/eigenlib.cpp
+++ /dev/null
@@ -1,107 +0,0 @@
-/*
-    Genome-wide Efficient Mixed Model Association (GEMMA)
-    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
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    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/>.
-*/
-
-#include "Eigen/Dense"
-// #include "gsl/gsl_linalg.h"
-#include "gsl/gsl_matrix.h"
-// #include "gsl/gsl_vector.h"
-#include <cmath>
-#include <iostream>
-#include <vector>
-// #include <cblas.h>
-
-using namespace std;
-using namespace Eigen;
-
-
-// On two different clusters, compare eigen vs lapack/gsl:
-//
-// dgemm, 5x or 0.5x faster or slower than lapack, 5x or 4x faster than gsl
-// dgemv, 20x or 4x faster than gsl,
-// eigen, 1x or 0.3x slower than lapack
-// invert, 20x or 10x faster than lapack
-//
-void eigenlib_dgemm(const char *TransA, const char *TransB, const double alpha,
-                    const gsl_matrix *A, const gsl_matrix *B, const double beta,
-                    gsl_matrix *C) {
-  Map<Matrix<double, Dynamic, Dynamic, RowMajor>, 0, OuterStride<Dynamic>>
-      A_mat(A->data, A->size1, A->size2, OuterStride<Dynamic>(A->tda));
-  Map<Matrix<double, Dynamic, Dynamic, RowMajor>, 0, OuterStride<Dynamic>>
-      B_mat(B->data, B->size1, B->size2, OuterStride<Dynamic>(B->tda));
-  Map<Matrix<double, Dynamic, Dynamic, RowMajor>, 0, OuterStride<Dynamic>>
-      C_mat(C->data, C->size1, C->size2, OuterStride<Dynamic>(C->tda));
-
-  if (*TransA == 'N' || *TransA == 'n') {
-    if (*TransB == 'N' || *TransB == 'n') {
-      C_mat = alpha * A_mat * B_mat + beta * C_mat;
-    } else {
-      C_mat = alpha * A_mat * B_mat.transpose() + beta * C_mat;
-    }
-  } else {
-    if (*TransB == 'N' || *TransB == 'n') {
-      C_mat = alpha * A_mat.transpose() * B_mat + beta * C_mat;
-    } else {
-      C_mat = alpha * A_mat.transpose() * B_mat.transpose() + beta * C_mat;
-    }
-  }
-}
-
-void eigenlib_dgemv(const char *TransA, const double alpha, const gsl_matrix *A,
-                    const gsl_vector *x, const double beta, gsl_vector *y) {
-  Map<Matrix<double, Dynamic, Dynamic, RowMajor>, 0, OuterStride<Dynamic>>
-      A_mat(A->data, A->size1, A->size2, OuterStride<Dynamic>(A->tda));
-  Map<Matrix<double, Dynamic, 1>, 0, InnerStride<Dynamic>> x_vec(
-      x->data, x->size, InnerStride<Dynamic>(x->stride));
-  Map<Matrix<double, Dynamic, 1>, 0, InnerStride<Dynamic>> y_vec(
-      y->data, y->size, InnerStride<Dynamic>(y->stride));
-
-  if (*TransA == 'N' || *TransA == 'n') {
-    y_vec = alpha * A_mat * x_vec + beta * y_vec;
-  } else {
-    y_vec = alpha * A_mat.transpose() * x_vec + beta * y_vec;
-  }
-}
-
-void eigenlib_invert(gsl_matrix *A) {
-  Map<Matrix<double, Dynamic, Dynamic, RowMajor>> A_mat(A->data, A->size1,
-                                                        A->size2);
-  A_mat = A_mat.inverse();
-}
-
-void eigenlib_dsyr(const double alpha, const gsl_vector *b, gsl_matrix *A) {
-  Map<Matrix<double, Dynamic, Dynamic, RowMajor>> A_mat(A->data, A->size1,
-                                                        A->size2);
-  Map<Matrix<double, Dynamic, 1>, 0, OuterStride<Dynamic>> b_vec(
-      b->data, b->size, OuterStride<Dynamic>(b->stride));
-  A_mat = alpha * b_vec * b_vec.transpose() + A_mat;
-}
-
-void eigenlib_eigensymm(const gsl_matrix *G, gsl_matrix *U, gsl_vector *eval) {
-  Map<Matrix<double, Dynamic, Dynamic, RowMajor>, 0, OuterStride<Dynamic>>
-      G_mat(G->data, G->size1, G->size2, OuterStride<Dynamic>(G->tda));
-  Map<Matrix<double, Dynamic, Dynamic, RowMajor>, 0, OuterStride<Dynamic>>
-      U_mat(U->data, U->size1, U->size2, OuterStride<Dynamic>(U->tda));
-  Map<Matrix<double, Dynamic, 1>, 0, OuterStride<Dynamic>> eval_vec(
-      eval->data, eval->size, OuterStride<Dynamic>(eval->stride));
-
-  SelfAdjointEigenSolver<MatrixXd> es(G_mat);
-  if (es.info() != Success)
-    abort();
-  eval_vec = es.eigenvalues();
-  U_mat = es.eigenvectors();
-}
diff --git a/src/eigenlib.h b/src/eigenlib.h
index 7fb69ad..e69de29 100644
--- a/src/eigenlib.h
+++ b/src/eigenlib.h
@@ -1,35 +0,0 @@
-/*
-    Genome-wide Efficient Mixed Model Association (GEMMA)
-    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
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    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/>.
-*/
-
-#ifndef __EIGENLIB_H__
-#define __EIGENLIB_H__
-
-// #include <vector>
-
-// using namespace std;
-
-void eigenlib_dgemm(const char *TransA, const char *TransB, const double alpha,
-                    const gsl_matrix *A, const gsl_matrix *B, const double beta,
-                    gsl_matrix *C);
-void eigenlib_dgemv(const char *TransA, const double alpha, const gsl_matrix *A,
-                    const gsl_vector *x, const double beta, gsl_vector *y);
-void eigenlib_invert(gsl_matrix *A);
-void eigenlib_dsyr(const double alpha, const gsl_vector *b, gsl_matrix *A);
-void eigenlib_eigensymm(const gsl_matrix *G, gsl_matrix *U, gsl_vector *eval);
-
-#endif
diff --git a/src/fastblas.cpp b/src/fastblas.cpp
index b3fcddb..e9e38d0 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>
@@ -235,8 +235,39 @@ 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>
+
+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);
+}
+
+  
diff --git a/src/fastblas.h b/src/fastblas.h
index 343a73a..cca3258 100644
--- a/src/fastblas.h
+++ b/src/fastblas.h
@@ -37,4 +37,6 @@ 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);
 
+void fast_inverse(gsl_matrix *m);
+
 #endif
diff --git a/src/fastopenblas.h b/src/fastopenblas.h
index e98e62f..2fc743e 100644
--- a/src/fastopenblas.h
+++ b/src/fastopenblas.h
@@ -25,9 +25,9 @@
 #include <iostream>
 extern "C"
 {
-   #include <cblas.h>   // For OpenBlas / Atlas
+  // #include <cblas.h>   // For OpenBlas / Atlas
 }
-#include "gsl/gsl_matrix.h"
+// #include "gsl/gsl_matrix.h"
 
 void fast_cblas_dgemm(const enum CBLAS_ORDER Order,
                       const enum CBLAS_TRANSPOSE TransA,
diff --git a/src/gemma.cpp b/src/gemma.cpp
index 0ced70c..6bc739b 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -2697,25 +2697,25 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
         CalcLambda('L', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max,
                    cPar.n_region, cPar.l_mle_null, cPar.logl_mle_H0);
-        assert(!std::isnan(UtY->data[0]));
+        assert(!isnan(UtY->data[0]));
 
         CalcLmmVgVeBeta(eval, UtW, &UtY_col.vector, cPar.l_mle_null,
                         cPar.vg_mle_null, cPar.ve_mle_null, &beta.vector,
                         &se_beta.vector);
 
-        assert(!std::isnan(UtY->data[0]));
+        assert(!isnan(UtY->data[0]));
 
         cPar.beta_mle_null.clear();
         cPar.se_beta_mle_null.clear();
-        assert(!std::isnan(B->data[0]));
-        assert(!std::isnan(se_B->data[0]));
+        assert(!isnan(B->data[0]));
+        assert(!isnan(se_B->data[0]));
         for (size_t i = 0; i < B->size2; i++) {
           cPar.beta_mle_null.push_back(gsl_matrix_get(B, 0, i));
           cPar.se_beta_mle_null.push_back(gsl_matrix_get(se_B, 0, i));
         }
-        assert(!std::isnan(UtY->data[0]));
-        assert(!std::isnan(cPar.beta_mle_null.front()));
-        assert(!std::isnan(cPar.se_beta_mle_null.front()));
+        assert(!isnan(UtY->data[0]));
+        assert(!isnan(cPar.beta_mle_null.front()));
+        assert(!isnan(cPar.se_beta_mle_null.front()));
 
         // the following functions do not modify eval
         CalcLambda('R', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max,
@@ -2726,8 +2726,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
         cPar.beta_remle_null.clear();
         cPar.se_beta_remle_null.clear();
-        assert(!std::isnan(B->data[0]));
-        assert(!std::isnan(se_B->data[0]));
+        assert(!isnan(B->data[0]));
+        assert(!isnan(se_B->data[0]));
 
         for (size_t i = 0; i < B->size2; i++) {
           cPar.beta_remle_null.push_back(gsl_matrix_get(B, 0, i));
@@ -3124,7 +3124,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
   return;
 }
 
-#include "Eigen/Dense"
+// #include "Eigen/Dense"
 
 void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) {
   string file_str;
@@ -3144,7 +3144,7 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) {
   outfile << "## GCC version      = " << __GNUC__ << "." << __GNUC_MINOR__ << "." << __GNUC_PATCHLEVEL__ << endl;
   #endif
   outfile << "## GSL Version      = " << GSL_VERSION << endl;
-  outfile << "## Eigen Version    = " << EIGEN_WORLD_VERSION << "." << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION << endl;
+  // outfile << "## Eigen Version    = " << EIGEN_WORLD_VERSION << "." << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION << endl;
 #ifdef OPENBLAS
 
   #ifndef OPENBLAS_LEGACY
diff --git a/src/gemma_io.cpp b/src/gemma_io.cpp
index 3c71ea6..2380d45 100644
--- a/src/gemma_io.cpp
+++ b/src/gemma_io.cpp
@@ -1470,7 +1470,7 @@ bool BimbamKin(const string file_geno, const set<string> ksnps,
         gsl_vector_set(geno_miss, i, 0);
         n_miss++;
       } else {
-        double d = stod(field);
+        double d = atof(field);
         if (is_strict_mode() && d == 0.0)
           enforce_is_float(std::string(field));  // rule out non NA and non-float fields
         gsl_vector_set(geno, i, d);
diff --git a/src/lmm.cpp b/src/lmm.cpp
index 16fc3c8..b811594 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -1582,7 +1582,7 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
         continue;
 
       double geno = gs[i];
-      if (std::isnan(geno)) {
+      if (isnan(geno)) {
         gsl_vector_set(x_miss, pos, 1.0);
         n_miss++;
       } else {
@@ -1829,7 +1829,7 @@ void LMM::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 de0c257..31b7382 100644
--- a/src/param.cpp
+++ b/src/param.cpp
@@ -39,6 +39,7 @@
 #include "gemma_io.h"
 #include "mathfunc.h"
 #include "param.h"
+#include "fastblas.h"
 
 using namespace std;
 
@@ -1412,7 +1413,7 @@ void compKtoV(const gsl_matrix *G, gsl_matrix *V) {
           gsl_matrix_const_submatrix(G, 0, j * ni_test, ni_test, ni_test);
       gsl_matrix_view KiKj_sub =
           gsl_matrix_submatrix(KiKj, 0, t * ni_test, ni_test, ni_test);
-      eigenlib_dgemm("N", "N", 1.0, &Ki.matrix, &Kj.matrix, 0.0,
+      fast_dgemm("N", "N", 1.0, &Ki.matrix, &Kj.matrix, 0.0,
                      &KiKj_sub.matrix);
       t++;
     }
diff --git a/src/vc.cpp b/src/vc.cpp
index 19590f7..22aaea9 100644
--- a/src/vc.cpp
+++ b/src/vc.cpp
@@ -43,7 +43,6 @@
 
 // #include "Eigen/Dense"
 
-#include "eigenlib.h"
 #include "gzstream.h"
 #include "gemma_io.h"
 #include "lapack.h"
@@ -51,6 +50,7 @@
 #include "mathfunc.h"
 #include "param.h"
 #include "vc.h"
+#include "fastblas.h"
 
 using namespace std;
 // using namespace Eigen;
@@ -198,16 +198,16 @@ void UpdateParam(const gsl_vector *log_sigma2, VC_PARAM *p) {
   }
 
   // Calculate H^{-1}.
-  eigenlib_invert(p->P);
+  fast_inverse(p->P);
 
-  eigenlib_dgemm("N", "N", 1.0, p->P, p->W, 0.0, HiW);
-  eigenlib_dgemm("T", "N", 1.0, p->W, HiW, 0.0, WtHiW);
+  fast_dgemm("N", "N", 1.0, p->P, p->W, 0.0, HiW);
+  fast_dgemm("T", "N", 1.0, p->W, HiW, 0.0, WtHiW);
 
-  eigenlib_invert(WtHiW);
+  fast_inverse(WtHiW);
   gsl_matrix_memcpy(WtHiWi, WtHiW);
 
-  eigenlib_dgemm("N", "T", 1.0, WtHiWi, HiW, 0.0, WtHiWiWtHi);
-  eigenlib_dgemm("N", "N", -1.0, HiW, WtHiWiWtHi, 1.0, p->P);
+  fast_dgemm("N", "T", 1.0, WtHiWi, HiW, 0.0, WtHiWiWtHi);
+  fast_dgemm("N", "N", -1.0, HiW, WtHiWiWtHi, 1.0, p->P);
 
   // Calculate Py, KPy, PKPy.
   gsl_blas_dgemv(CblasNoTrans, 1.0, p->P, p->y, 0.0, p->Py);
@@ -224,7 +224,7 @@ void UpdateParam(const gsl_vector *log_sigma2, VC_PARAM *p) {
           gsl_matrix_const_submatrix(p->K, 0, n1 * i, n1, n1);
 
       // Seems to be important to use gsl dgemv here instead of
-      // eigenlib_dgemv; otherwise.
+      // fast_dgemv; otherwise.
       gsl_blas_dgemv(CblasNoTrans, 1.0, &K_sub.matrix, p->Py, 0.0, &KPy.vector);
     }
 
@@ -232,15 +232,15 @@ void UpdateParam(const gsl_vector *log_sigma2, VC_PARAM *p) {
 
     // When phenotypes are not normalized well, then some values in
     // the following matrix maybe NaN; change that to 0; this seems to
-    // only happen when eigenlib_dgemv was used above.
+    // only happen when fast_dgemv was used above.
     for (size_t j = 0; j < p->KPy_mat->size1; j++) {
       d = gsl_matrix_get(p->KPy_mat, j, i);
-      if (std::isnan(d)) {
+      if (isnan(d)) {
         gsl_matrix_set(p->KPy_mat, j, i, 0);
         cout << "nan appears in " << i << " " << j << endl;
       }
       d = gsl_matrix_get(p->PKPy_mat, j, i);
-      if (std::isnan(d)) {
+      if (isnan(d)) {
         gsl_matrix_set(p->PKPy_mat, j, i, 0);
         cout << "nan appears in " << i << " " << j << endl;
       }
diff --git a/test/src/unittests-math.cpp b/test/src/unittests-math.cpp
index d44d8c4..d82f656 100644
--- a/test/src/unittests-math.cpp
+++ b/test/src/unittests-math.cpp
@@ -2,7 +2,7 @@
 #include <iostream>
 #include <fenv.h>
 #include "gsl/gsl_matrix.h"
-#include <cblas.h>
+// #include <cblas.h>
 
 #include <algorithm>
 #include <limits>
@@ -63,13 +63,13 @@ TEST_CASE( "Math functions", "[math]" ) {
 
   // ---- NaN checks
   vector<double> v = {1.0, 2.0};
-  REQUIRE (!std::isnan(std::accumulate(v.begin(), v.end(), 0)));
+  // REQUIRE (!isnan(std::accumulate(v.begin(), v.end(), 0)));
   vector<double> v2 = {1.0, 2.0, std::numeric_limits<double>::quiet_NaN()};
-  REQUIRE (std::isnan(v2[2]));
+  REQUIRE (isnan(v2[2]));
   REQUIRE(has_nan(v2));
   // test minus nan
   vector<double> v3 = {1.0, 2.0, -std::numeric_limits<double>::quiet_NaN()};
-  REQUIRE (std::isnan(v3[2]));
+  REQUIRE (isnan(v3[2]));
   REQUIRE(has_nan(v3));
 }
generated by cgit-pink 1.4.1 (git 2.36.1) at 2025-10-07 07:46:31 +0000