From c1132606169875be6d07b54b30e8ae9446341bc2 Mon Sep 17 00:00:00 2001
From: Peter Carbonetto
Date: Sun, 4 Jun 2017 12:06:36 -0500
Subject: Removed FORCE_FLOAT from prdt.h/prdt.cpp.
---
src/eigenlib.cpp | 94 ++++++++++++++++++++++++++++----------------------------
1 file changed, 47 insertions(+), 47 deletions(-)
(limited to 'src/eigenlib.cpp')
diff --git a/src/eigenlib.cpp b/src/eigenlib.cpp
index 14ffbf1..7ad250f 100644
--- a/src/eigenlib.cpp
+++ b/src/eigenlib.cpp
@@ -1,6 +1,6 @@
/*
- Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011 Xiang Zhou
+ 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
@@ -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 .
+ along with this program. If not, see .
*/
#include
@@ -27,19 +27,23 @@
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, 0, OuterStride > A_mat(A->data, A->size1, A->size2, OuterStride(A->tda) );
- Map, 0, OuterStride > B_mat(B->data, B->size1, B->size2, OuterStride(B->tda) );
- Map, 0, OuterStride > C_mat(C->data, C->size1, C->size2, OuterStride(C->tda) );
+// 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, 0, OuterStride >
+ A_mat(A->data, A->size1, A->size2, OuterStride(A->tda) );
+ Map, 0, OuterStride >
+ B_mat(B->data, B->size1, B->size2, OuterStride(B->tda) );
+ Map, 0, OuterStride >
+ C_mat(C->data, C->size1, C->size2, OuterStride(C->tda) );
if (*TransA=='N' || *TransA=='n') {
if (*TransB=='N' || *TransB=='n') {
@@ -55,19 +59,18 @@ void eigenlib_dgemm (const char *TransA, const char *TransB, const double alpha,
}
}
- //gsl_matrix_view C_view = gsl_matrix_view_array (C_mat.data(), C->size1, C->size2);
- //gsl_matrix_memcpy (C, &C_view.matrix);
-
return;
}
-
-
-void eigenlib_dgemv (const char *TransA, const double alpha, const gsl_matrix *A, const gsl_vector *x, const double beta, gsl_vector *y)
-{
- Map, 0, OuterStride > A_mat(A->data, A->size1, A->size2, OuterStride(A->tda) );
- Map, 0, InnerStride > x_vec(x->data, x->size, InnerStride(x->stride) );
- Map, 0, InnerStride > y_vec(y->data, y->size, InnerStride(y->stride) );
+void eigenlib_dgemv (const char *TransA, const double alpha,
+ const gsl_matrix *A, const gsl_vector *x,
+ const double beta, gsl_vector *y) {
+ Map, 0, OuterStride >
+ A_mat(A->data, A->size1, A->size2, OuterStride(A->tda) );
+ Map, 0, InnerStride >
+ x_vec(x->data, x->size, InnerStride(x->stride) );
+ Map, 0, InnerStride >
+ y_vec(y->data, y->size, InnerStride(y->stride) );
if (*TransA=='N' || *TransA=='n') {
y_vec=alpha*A_mat*x_vec+beta*y_vec;
@@ -78,38 +81,35 @@ void eigenlib_dgemv (const char *TransA, const double alpha, const gsl_matrix *A
return;
}
-
-
-void eigenlib_invert(gsl_matrix *A)
-{
- Map > A_mat(A->data, A->size1, A->size2);
+void eigenlib_invert(gsl_matrix *A) {
+ Map >
+ A_mat(A->data, A->size1, A->size2);
A_mat=A_mat.inverse();
return;
}
-
-void eigenlib_dsyr (const double alpha, const gsl_vector *b, gsl_matrix *A)
-{
- Map > A_mat(A->data, A->size1, A->size2);
- Map, 0, OuterStride > b_vec(b->data, b->size, OuterStride(b->stride) );
-
+void eigenlib_dsyr (const double alpha, const gsl_vector *b, gsl_matrix *A) {
+ Map >
+ A_mat(A->data, A->size1, A->size2);
+ Map, 0, OuterStride >
+ b_vec(b->data, b->size, OuterStride(b->stride) );
A_mat=alpha*b_vec*b_vec.transpose()+A_mat;
-
return;
}
-
-void eigenlib_eigensymm (const gsl_matrix *G, gsl_matrix *U, gsl_vector *eval)
-{
- Map, 0, OuterStride > G_mat(G->data, G->size1, G->size2, OuterStride(G->tda) );
- Map, 0, OuterStride > U_mat(U->data, U->size1, U->size2, OuterStride(U->tda) );
- Map, 0, OuterStride > eval_vec(eval->data, eval->size, OuterStride(eval->stride) );
+void eigenlib_eigensymm (const gsl_matrix *G, gsl_matrix *U,
+ gsl_vector *eval) {
+ Map, 0, OuterStride >
+ G_mat(G->data, G->size1, G->size2, OuterStride(G->tda) );
+ Map, 0, OuterStride >
+ U_mat(U->data, U->size1, U->size2, OuterStride(U->tda) );
+ Map, 0, OuterStride >
+ eval_vec(eval->data, eval->size, OuterStride(eval->stride) );
SelfAdjointEigenSolver es(G_mat);
- if (es.info() != Success) abort();
-
+ if (es.info() != Success)
+ abort();
eval_vec=es.eigenvalues();
U_mat=es.eigenvectors();
-
return;
}
--
cgit v1.2.3