aboutsummaryrefslogtreecommitdiff
path: root/src/vc.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/vc.cpp')
-rw-r--r--src/vc.cpp22
1 files changed, 11 insertions, 11 deletions
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;
}