diff options
Diffstat (limited to 'src/lapack.cpp')
| -rw-r--r-- | src/lapack.cpp | 445 |
1 files changed, 264 insertions, 181 deletions
diff --git a/src/lapack.cpp b/src/lapack.cpp index d9e4149..18399a2 100644 --- a/src/lapack.cpp +++ b/src/lapack.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 <http://www.gnu.org/licenses/>. + along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include <iostream> @@ -25,82 +25,129 @@ using namespace std; -extern "C" void sgemm_(char *TRANSA, char *TRANSB, int *M, int *N, int *K, float *ALPHA, float *A, int *LDA, float *B, int *LDB, float *BETA, float *C, int *LDC); +extern "C" void sgemm_(char *TRANSA, char *TRANSB, int *M, int *N, int *K, + float *ALPHA, float *A, int *LDA, float *B, int *LDB, + float *BETA, float *C, int *LDC); extern "C" void spotrf_(char *UPLO, int *N, float *A, int *LDA, int *INFO); -extern "C" void spotrs_(char *UPLO, int *N, int *NRHS, float *A, int *LDA, float *B, int *LDB, int *INFO); -extern "C" void ssyev_(char* JOBZ, char* UPLO, int *N, float *A, int *LDA, float *W, float *WORK, int *LWORK, int *INFO); -extern "C" void ssyevr_(char* JOBZ, char *RANGE, char* UPLO, int *N, float *A, int *LDA, float *VL, float *VU, int *IL, int *IU, float *ABSTOL, int *M, float *W, float *Z, int *LDZ, int *ISUPPZ, float *WORK, int *LWORK, int *IWORK, int *LIWORK, int *INFO); +extern "C" void spotrs_(char *UPLO, int *N, int *NRHS, float *A, int *LDA, + float *B, int *LDB, int *INFO); +extern "C" void ssyev_(char* JOBZ, char* UPLO, int *N, float *A, int *LDA, + float *W, float *WORK, int *LWORK, int *INFO); +extern "C" void ssyevr_(char* JOBZ, char *RANGE, char* UPLO, int *N, + float *A, int *LDA, float *VL, float *VU, int *IL, + int *IU, float *ABSTOL, int *M, float *W, float *Z, + int *LDZ, int *ISUPPZ, float *WORK, int *LWORK, + int *IWORK, int *LIWORK, int *INFO); extern "C" double sdot_(int *N, float *DX, int *INCX, float *DY, int *INCY); -extern "C" void dgemm_(char *TRANSA, char *TRANSB, int *M, int *N, int *K, double *ALPHA, double *A, int *LDA, double *B, int *LDB, double *BETA, double *C, int *LDC); +extern "C" void dgemm_(char *TRANSA, char *TRANSB, int *M, int *N, int *K, + double *ALPHA, double *A, int *LDA, double *B, + int *LDB, double *BETA, double *C, int *LDC); extern "C" void dpotrf_(char *UPLO, int *N, double *A, int *LDA, int *INFO); -extern "C" void dpotrs_(char *UPLO, int *N, int *NRHS, double *A, int *LDA, double *B, int *LDB, int *INFO); -extern "C" void dsyev_(char* JOBZ, char* UPLO, int *N, double *A, int *LDA, double *W, double *WORK, int *LWORK, int *INFO); -extern "C" void dsyevr_(char* JOBZ, char *RANGE, char* UPLO, int *N, double *A, int *LDA, double *VL, double *VU, int *IL, int *IU, double *ABSTOL, int *M, double *W, double *Z, int *LDZ, int *ISUPPZ, double *WORK, int *LWORK, int *IWORK, int *LIWORK, int *INFO); +extern "C" void dpotrs_(char *UPLO, int *N, int *NRHS, double *A, int *LDA, + double *B, int *LDB, int *INFO); +extern "C" void dsyev_(char* JOBZ, char* UPLO, int *N, double *A, int *LDA, + double *W, double *WORK, int *LWORK, int *INFO); +extern "C" void dsyevr_(char* JOBZ, char *RANGE, char* UPLO, int *N, + double *A, int *LDA, double *VL, double *VU, + int *IL, int *IU, double *ABSTOL, int *M, + double *W, double *Z, int *LDZ, int *ISUPPZ, + double *WORK, int *LWORK, int *IWORK, + int *LIWORK, int *INFO); extern "C" double ddot_(int *N, double *DX, int *INCX, double *DY, int *INCY); - -//cholesky decomposition, A is distroyed -void lapack_float_cholesky_decomp (gsl_matrix_float *A) -{ +// Cholesky decomposition, A is destroyed. +void lapack_float_cholesky_decomp (gsl_matrix_float *A) { int N=A->size1, LDA=A->size1, INFO; char UPLO='L'; - if (N!=(int)A->size2) {cout<<"Matrix needs to be symmetric and same dimension in lapack_cholesky_decomp."<<endl; return;} + if (N!=(int)A->size2) { + cout << "Matrix needs to be symmetric and same dimension in" << + "lapack_cholesky_decomp." << endl; + return; + } spotrf_(&UPLO, &N, A->data, &LDA, &INFO); - if (INFO!=0) {cout<<"Cholesky decomposition unsuccessful in lapack_cholesky_decomp."<<endl; return;} + if (INFO!=0) { + cout << "Cholesky decomposition unsuccessful in" << + "lapack_cholesky_decomp." << endl; + return; + } return; } -//cholesky decomposition, A is distroyed -void lapack_cholesky_decomp (gsl_matrix *A) -{ +// Cholesky decomposition, A is destroyed. +void lapack_cholesky_decomp (gsl_matrix *A) { int N=A->size1, LDA=A->size1, INFO; char UPLO='L'; - if (N!=(int)A->size2) {cout<<"Matrix needs to be symmetric and same dimension in lapack_cholesky_decomp."<<endl; return;} + if (N!=(int)A->size2) { + cout << "Matrix needs to be symmetric and same dimension in" << + "lapack_cholesky_decomp." << endl; + return; + } dpotrf_(&UPLO, &N, A->data, &LDA, &INFO); - if (INFO!=0) {cout<<"Cholesky decomposition unsuccessful in lapack_cholesky_decomp."<<endl; return;} - + if (INFO!=0) { + cout << "Cholesky decomposition unsuccessful in" << + "lapack_cholesky_decomp."<<endl; + return; + } + return; } -//cholesky solve, A is decomposed, -void lapack_float_cholesky_solve (gsl_matrix_float *A, const gsl_vector_float *b, gsl_vector_float *x) -{ +// Cholesky solve, A is decomposed. +void lapack_float_cholesky_solve (gsl_matrix_float *A, + const gsl_vector_float *b, + gsl_vector_float *x) { int N=A->size1, NRHS=1, LDA=A->size1, LDB=b->size, INFO; char UPLO='L'; - if (N!=(int)A->size2 || N!=LDB) {cout<<"Matrix needs to be symmetric and same dimension in lapack_cholesky_solve."<<endl; return;} + if (N!=(int)A->size2 || N!=LDB) { + cout << "Matrix needs to be symmetric and same dimension in" << + "lapack_cholesky_solve." << endl; + return; + } gsl_vector_float_memcpy (x, b); spotrs_(&UPLO, &N, &NRHS, A->data, &LDA, x->data, &LDB, &INFO); - if (INFO!=0) {cout<<"Cholesky solve unsuccessful in lapack_cholesky_solve."<<endl; return;} + if (INFO!=0) { + cout << "Cholesky solve unsuccessful in lapack_cholesky_solve." << + endl; + return; + } return; } -//cholesky solve, A is decomposed, -void lapack_cholesky_solve (gsl_matrix *A, const gsl_vector *b, gsl_vector *x) -{ +// Cholesky solve, A is decomposed. +void lapack_cholesky_solve (gsl_matrix *A, const gsl_vector *b, + gsl_vector *x) { int N=A->size1, NRHS=1, LDA=A->size1, LDB=b->size, INFO; char UPLO='L'; - if (N!=(int)A->size2 || N!=LDB) {cout<<"Matrix needs to be symmetric and same dimension in lapack_cholesky_solve."<<endl; return;} + if (N!=(int)A->size2 || N!=LDB) { + cout << "Matrix needs to be symmetric and same dimension in" << + "lapack_cholesky_solve." << endl; + return; + } gsl_vector_memcpy (x, b); dpotrs_(&UPLO, &N, &NRHS, A->data, &LDA, x->data, &LDB, &INFO); - if (INFO!=0) {cout<<"Cholesky solve unsuccessful in lapack_cholesky_solve."<<endl; return;} + if (INFO!=0) { + cout << "Cholesky solve unsuccessful in lapack_cholesky_solve." << + endl; + return; + } return; } - -void lapack_sgemm (char *TransA, char *TransB, float alpha, const gsl_matrix_float *A, const gsl_matrix_float *B, float beta, gsl_matrix_float *C) -{ +void lapack_sgemm (char *TransA, char *TransB, float alpha, + const gsl_matrix_float *A, const gsl_matrix_float *B, + float beta, gsl_matrix_float *C) { int M, N, K1, K2, LDA=A->size1, LDB=B->size1, LDC=C->size2; if (*TransA=='N' || *TransA=='n') {M=A->size1; K1=A->size2;} @@ -111,8 +158,14 @@ void lapack_sgemm (char *TransA, char *TransB, float alpha, const gsl_matrix_flo else if (*TransB=='T' || *TransB=='t') {N=B->size1; K2=B->size2;} else {cout<<"need 'N' or 'T' in lapack_sgemm"<<endl; return;} - if (K1!=K2) {cout<<"A and B not compatible in lapack_sgemm"<<endl; return;} - if (C->size1!=(size_t)M || C->size2!=(size_t)N) {cout<<"C not compatible in lapack_sgemm"<<endl; return;} + if (K1!=K2) { + cout<<"A and B not compatible in lapack_sgemm"<<endl; + return; + } + if (C->size1!=(size_t)M || C->size2!=(size_t)N) { + cout<<"C not compatible in lapack_sgemm"<<endl; + return; + } gsl_matrix_float *A_t=gsl_matrix_float_alloc (A->size2, A->size1); gsl_matrix_float_transpose_memcpy (A_t, A); @@ -121,7 +174,8 @@ void lapack_sgemm (char *TransA, char *TransB, float alpha, const gsl_matrix_flo gsl_matrix_float *C_t=gsl_matrix_float_alloc (C->size2, C->size1); gsl_matrix_float_transpose_memcpy (C_t, C); - sgemm_(TransA, TransB, &M, &N, &K1, &alpha, A_t->data, &LDA, B_t->data, &LDB, &beta, C_t->data, &LDC); + sgemm_(TransA, TransB, &M, &N, &K1, &alpha, A_t->data, &LDA, + B_t->data, &LDB, &beta, C_t->data, &LDC); gsl_matrix_float_transpose_memcpy (C, C_t); gsl_matrix_float_free (A_t); @@ -132,8 +186,9 @@ void lapack_sgemm (char *TransA, char *TransB, float alpha, const gsl_matrix_flo -void lapack_dgemm (char *TransA, char *TransB, double alpha, const gsl_matrix *A, const gsl_matrix *B, double beta, gsl_matrix *C) -{ +void lapack_dgemm (char *TransA, char *TransB, double alpha, + const gsl_matrix *A, const gsl_matrix *B, + double beta, gsl_matrix *C) { int M, N, K1, K2, LDA=A->size1, LDB=B->size1, LDC=C->size2; if (*TransA=='N' || *TransA=='n') {M=A->size1; K1=A->size2;} @@ -144,8 +199,14 @@ void lapack_dgemm (char *TransA, char *TransB, double alpha, const gsl_matrix *A else if (*TransB=='T' || *TransB=='t') {N=B->size1; K2=B->size2;} else {cout<<"need 'N' or 'T' in lapack_dgemm"<<endl; return;} - if (K1!=K2) {cout<<"A and B not compatible in lapack_dgemm"<<endl; return;} - if (C->size1!=(size_t)M || C->size2!=(size_t)N) {cout<<"C not compatible in lapack_dgemm"<<endl; return;} + if (K1!=K2) { + cout << "A and B not compatible in lapack_dgemm"<<endl; + return; + } + if (C->size1!=(size_t)M || C->size2!=(size_t)N) { + cout<<"C not compatible in lapack_dgemm"<<endl; + return; + } gsl_matrix *A_t=gsl_matrix_alloc (A->size2, A->size1); gsl_matrix_transpose_memcpy (A_t, A); @@ -154,7 +215,8 @@ void lapack_dgemm (char *TransA, char *TransB, double alpha, const gsl_matrix *A gsl_matrix *C_t=gsl_matrix_alloc (C->size2, C->size1); gsl_matrix_transpose_memcpy (C_t, C); - dgemm_(TransA, TransB, &M, &N, &K1, &alpha, A_t->data, &LDA, B_t->data, &LDB, &beta, C_t->data, &LDC); + dgemm_(TransA, TransB, &M, &N, &K1, &alpha, A_t->data, &LDA, + B_t->data, &LDB, &beta, C_t->data, &LDC); gsl_matrix_transpose_memcpy (C, C_t); @@ -164,58 +226,79 @@ void lapack_dgemm (char *TransA, char *TransB, double alpha, const gsl_matrix *A return; } - - -//eigen value decomposition, matrix A is destroyed, float seems to have problem with large matrices (in mac) -void lapack_float_eigen_symmv (gsl_matrix_float *A, gsl_vector_float *eval, gsl_matrix_float *evec, const size_t flag_largematrix) -{ +// Eigen value decomposition, matrix A is destroyed, float seems to +// have problem with large matrices (in mac). +void lapack_float_eigen_symmv (gsl_matrix_float *A, gsl_vector_float *eval, + gsl_matrix_float *evec, + const size_t flag_largematrix) { if (flag_largematrix==1) { int N=A->size1, LDA=A->size1, INFO, LWORK=-1; char JOBZ='V', UPLO='L'; - if (N!=(int)A->size2 || N!=(int)eval->size) {cout<<"Matrix needs to be symmetric and same dimension in lapack_eigen_symmv."<<endl; return;} - - // float temp[1]; - // ssyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, temp, &LWORK, &INFO); - // if (INFO!=0) {cout<<"Work space estimate unsuccessful in lapack_eigen_symmv."<<endl; return;} - // LWORK=(int)temp[0]; + if (N!=(int)A->size2 || N!=(int)eval->size) { + cout << "Matrix needs to be symmetric and same" << + "dimension in lapack_eigen_symmv."<<endl; + return; + } LWORK=3*N; float *WORK=new float [LWORK]; - ssyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK, &LWORK, &INFO); - if (INFO!=0) {cout<<"Eigen decomposition unsuccessful in lapack_eigen_symmv."<<endl; return;} + ssyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK, + &LWORK, &INFO); + if (INFO!=0) { + cout << "Eigen decomposition unsuccessful in" << + "lapack_eigen_symmv."<<endl; + return; + } - gsl_matrix_float_view A_sub=gsl_matrix_float_submatrix(A, 0, 0, N, N); + gsl_matrix_float_view A_sub = + gsl_matrix_float_submatrix(A, 0, 0, N, N); gsl_matrix_float_memcpy (evec, &A_sub.matrix); gsl_matrix_float_transpose (evec); delete [] WORK; } else { - int N=A->size1, LDA=A->size1, LDZ=A->size1, INFO, LWORK=-1, LIWORK=-1; + int N=A->size1, LDA=A->size1, LDZ=A->size1, INFO, + LWORK=-1, LIWORK=-1; char JOBZ='V', UPLO='L', RANGE='A'; float ABSTOL=1.0E-7; - //VL, VU, IL, IU are not referenced; M equals N if RANGE='A' + // VL, VU, IL, IU are not referenced; M equals N if RANGE='A'. float VL=0.0, VU=0.0; int IL=0, IU=0, M; - if (N!=(int)A->size2 || N!=(int)eval->size) {cout<<"Matrix needs to be symmetric and same dimension in lapack_float_eigen_symmv."<<endl; return;} + if (N!=(int)A->size2 || N!=(int)eval->size) { + cout << "Matrix needs to be symmetric and same" << + "dimension in lapack_float_eigen_symmv." << endl; + return; + } int *ISUPPZ=new int [2*N]; float WORK_temp[1]; int IWORK_temp[1]; - ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK_temp, &LWORK, IWORK_temp, &LIWORK, &INFO); - if (INFO!=0) {cout<<"Work space estimate unsuccessful in lapack_float_eigen_symmv."<<endl; return;} + ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, + &VU, &IL, &IU, &ABSTOL, &M, eval->data, + evec->data, &LDZ, ISUPPZ, WORK_temp, &LWORK, + IWORK_temp, &LIWORK, &INFO); + if (INFO!=0) { + cout << "Work space estimate unsuccessful in" << + "lapack_float_eigen_symmv." << endl; + return; + } LWORK=(int)WORK_temp[0]; LIWORK=(int)IWORK_temp[0]; - //LWORK=26*N; - //LIWORK=10*N; float *WORK=new float [LWORK]; int *IWORK=new int [LIWORK]; - ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK, &LWORK, IWORK, &LIWORK, &INFO); - if (INFO!=0) {cout<<"Eigen decomposition unsuccessful in lapack_float_eigen_symmv."<<endl; return;} + ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, + &VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, + &LDZ, ISUPPZ, WORK, &LWORK, IWORK, &LIWORK, &INFO); + if (INFO!=0) { + cout << "Eigen decomposition unsuccessful in" << + "lapack_float_eigen_symmv." << endl; + return; + } gsl_matrix_float_transpose (evec); @@ -230,24 +313,28 @@ void lapack_float_eigen_symmv (gsl_matrix_float *A, gsl_vector_float *eval, gsl_ -//eigen value decomposition, matrix A is destroyed -void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, const size_t flag_largematrix) -{ +// Eigenvalue decomposition, matrix A is destroyed. +void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, + const size_t flag_largematrix) { if (flag_largematrix==1) { int N=A->size1, LDA=A->size1, INFO, LWORK=-1; char JOBZ='V', UPLO='L'; - if (N!=(int)A->size2 || N!=(int)eval->size) {cout<<"Matrix needs to be symmetric and same dimension in lapack_eigen_symmv."<<endl; return;} - - // double temp[1]; - // dsyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, temp, &LWORK, &INFO); - // if (INFO!=0) {cout<<"Work space estimate unsuccessful in lapack_eigen_symmv."<<endl; return;} - // LWORK=(int)temp[0]; + if (N!=(int)A->size2 || N!=(int)eval->size) { + cout << "Matrix needs to be symmetric and same" << + "dimension in lapack_eigen_symmv." << endl; + return; + } LWORK=3*N; double *WORK=new double [LWORK]; - dsyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK, &LWORK, &INFO); - if (INFO!=0) {cout<<"Eigen decomposition unsuccessful in lapack_eigen_symmv."<<endl; return;} + dsyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK, + &LWORK, &INFO); + if (INFO!=0) { + cout<<"Eigen decomposition unsuccessful in" << + "lapack_eigen_symmv." << endl; + return; + } gsl_matrix_view A_sub=gsl_matrix_submatrix(A, 0, 0, N, N); gsl_matrix_memcpy (evec, &A_sub.matrix); @@ -255,32 +342,48 @@ void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, cons delete [] WORK; } else { - int N=A->size1, LDA=A->size1, LDZ=A->size1, INFO, LWORK=-1, LIWORK=-1; + int N=A->size1, LDA=A->size1, LDZ=A->size1, INFO; + int LWORK=-1, LIWORK=-1; char JOBZ='V', UPLO='L', RANGE='A'; double ABSTOL=1.0E-7; - //VL, VU, IL, IU are not referenced; M equals N if RANGE='A' + // VL, VU, IL, IU are not referenced; M equals N if RANGE='A'. double VL=0.0, VU=0.0; int IL=0, IU=0, M; - if (N!=(int)A->size2 || N!=(int)eval->size) {cout<<"Matrix needs to be symmetric and same dimension in lapack_eigen_symmv."<<endl; return;} + if (N!=(int)A->size2 || N!=(int)eval->size) { + cout << "Matrix needs to be symmetric and same" << + "dimension in lapack_eigen_symmv." << endl; + return; + } int *ISUPPZ=new int [2*N]; double WORK_temp[1]; int IWORK_temp[1]; - dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK_temp, &LWORK, IWORK_temp, &LIWORK, &INFO); - if (INFO!=0) {cout<<"Work space estimate unsuccessful in lapack_eigen_symmv."<<endl; return;} + dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, + &IL, &IU, &ABSTOL, &M, eval->data, evec->data, + &LDZ, ISUPPZ, WORK_temp, &LWORK, IWORK_temp, + &LIWORK, &INFO); + if (INFO!=0) { + cout << "Work space estimate unsuccessful in" << + "lapack_eigen_symmv." << endl; + return; + } LWORK=(int)WORK_temp[0]; LIWORK=(int)IWORK_temp[0]; - //LWORK=26*N; - //LIWORK=10*N; double *WORK=new double [LWORK]; int *IWORK=new int [LIWORK]; - dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK, &LWORK, IWORK, &LIWORK, &INFO); - if (INFO!=0) {cout<<"Eigen decomposition unsuccessful in lapack_eigen_symmv."<<endl; return;} + dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, + &IL, &IU, &ABSTOL, &M, eval->data, evec->data, + &LDZ, ISUPPZ, WORK, &LWORK, IWORK, &LIWORK, &INFO); + if (INFO!=0) { + cout << "Eigen decomposition unsuccessful in" << + "lapack_eigen_symmv." << endl; + return; + } gsl_matrix_transpose (evec); @@ -292,9 +395,9 @@ void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, cons return; } -//DO NOT set eigen values to be positive -double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, const size_t flag_largematrix) -{ +// DO NOT set eigenvalues to be positive. +double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, + const size_t flag_largematrix) { #ifdef WITH_LAPACK lapack_eigen_symmv (G, eval, U, flag_largematrix); #else @@ -302,15 +405,8 @@ double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, const size_t gsl_eigen_symmv (G, eval, U, w); gsl_eigen_symmv_free (w); #endif - /* - for (size_t i=0; i<eval->size; ++i) { - if (gsl_vector_get (eval, i)<1e-10) { -// cout<<gsl_vector_get (eval, i)<<endl; - gsl_vector_set (eval, i, 0); - } - } - */ - //calculate track_G=mean(diag(G)) + + // Calculate track_G=mean(diag(G)). double d=0.0; for (size_t i=0; i<eval->size; ++i) { d+=gsl_vector_get(eval, i); @@ -321,55 +417,54 @@ double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, const size_t } -//DO NOT set eigen values to be positive -double EigenDecomp (gsl_matrix_float *G, gsl_matrix_float *U, gsl_vector_float *eval, const size_t flag_largematrix) -{ +// DO NOT set eigen values to be positive. +double EigenDecomp (gsl_matrix_float *G, gsl_matrix_float *U, + gsl_vector_float *eval, const size_t flag_largematrix) { #ifdef WITH_LAPACK lapack_float_eigen_symmv (G, eval, U, flag_largematrix); #else - //gsl doesn't provide float precision eigen decomposition; plus, float precision eigen decomposition in lapack may not work on OS 10.4 - //first change to double precision + + // GSL doesn't provide single precision eigen decomposition; + // plus, float precision eigen decomposition in LAPACK may not + // work on OS 10.4 first change to double precision. gsl_matrix *G_double=gsl_matrix_alloc (G->size1, G->size2); gsl_matrix *U_double=gsl_matrix_alloc (U->size1, U->size2); gsl_vector *eval_double=gsl_vector_alloc (eval->size); for (size_t i=0; i<G->size1; i++) { for (size_t j=0; j<G->size2; j++) { - gsl_matrix_set(G_double, i, j, gsl_matrix_float_get(G, i, j)); + gsl_matrix_set(G_double, i, j, + gsl_matrix_float_get(G, i, j)); } } gsl_eigen_symmv_workspace *w_space=gsl_eigen_symmv_alloc (G->size1); gsl_eigen_symmv (G_double, eval_double, U_double, w_space); gsl_eigen_symmv_free (w_space); - //change back to float precision + // Change back to float precision. for (size_t i=0; i<G->size1; i++) { for (size_t j=0; j<G->size2; j++) { - gsl_matrix_float_set(K, i, j, gsl_matrix_get(G_double, i, j)); + gsl_matrix_float_set(K, i, j, + gsl_matrix_get(G_double, i, j)); } } for (size_t i=0; i<U->size1; i++) { for (size_t j=0; j<U->size2; j++) { - gsl_matrix_float_set(U, i, j, gsl_matrix_get(U_double, i, j)); + gsl_matrix_float_set(U, i, j, + gsl_matrix_get(U_double, i, j)); } } for (size_t i=0; i<eval->size; i++) { gsl_vector_float_set(eval, i, gsl_vector_get(eval_double, i)); } - //delete double precision matrices + // Delete double-precision matrices. gsl_matrix_free (G_double); gsl_matrix_free (U_double); gsl_vector_free (eval_double); #endif - /* - for (size_t i=0; i<eval->size; ++i) { - if (gsl_vector_float_get (eval, i)<1e-10) { - gsl_vector_float_set (eval, i, 0); - } - } - */ - //calculate track_G=mean(diag(G)) - double d=0.0; + + // Calculate track_G=mean(diag(G)). + double d = 0.0; for (size_t i=0; i<eval->size; ++i) { d+=gsl_vector_float_get(eval, i); } @@ -379,8 +474,7 @@ double EigenDecomp (gsl_matrix_float *G, gsl_matrix_float *U, gsl_vector_float * } -double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty) -{ +double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty) { double logdet_O=0.0; #ifdef WITH_LAPACK @@ -394,7 +488,6 @@ double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty) int status = gsl_linalg_cholesky_decomp(Omega); if(status == GSL_EDOM) { cout << "## non-positive definite matrix" << endl; - // exit(0); } for (size_t i=0; i<Omega->size1; ++i) { @@ -404,16 +497,15 @@ double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty) gsl_vector_memcpy (OiXty, Xty); gsl_blas_dtrsv(CblasLower, CblasNoTrans, CblasNonUnit, Omega, OiXty); - gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, OiXty); - // gsl_linalg_cholesky_solve(XtX, Xty, iXty); + gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, OiXty); #endif return logdet_O; } -double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, gsl_vector_float *OiXty) -{ +double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, + gsl_vector_float *OiXty) { double logdet_O=0.0; #ifdef WITH_LAPACK @@ -436,7 +528,6 @@ double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, gsl_vector_ int status = gsl_linalg_cholesky_decomp(Omega_double); if(status == GSL_EDOM) { cout << "## non-positive definite matrix" << endl; - // exit(0); } for (size_t i=0; i<Omega->size1; ++i) { @@ -450,8 +541,7 @@ double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, gsl_vector_ gsl_vector_float_memcpy (OiXty, Xty); gsl_blas_strsv(CblasLower, CblasNoTrans, CblasNonUnit, Omega, OiXty); - gsl_blas_strsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, OiXty); - // gsl_linalg_cholesky_solve(XtX, Xty, iXty); + gsl_blas_strsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, OiXty); gsl_matrix_free (Omega_double); #endif @@ -460,129 +550,124 @@ double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, gsl_vector_ } -//LU decomposition -void LUDecomp (gsl_matrix *LU, gsl_permutation *p, int *signum) -{ +// LU decomposition. +void LUDecomp (gsl_matrix *LU, gsl_permutation *p, int *signum) { gsl_linalg_LU_decomp (LU, p, signum); return; } -void LUDecomp (gsl_matrix_float *LU, gsl_permutation *p, int *signum) -{ +void LUDecomp (gsl_matrix_float *LU, gsl_permutation *p, int *signum) { gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2); - //copy float matrix to double + // Copy float matrix to double. for (size_t i=0; i<LU->size1; i++) { for (size_t j=0; j<LU->size2; j++) { - gsl_matrix_set (LU_double, i, j, gsl_matrix_float_get(LU, i, j)); + gsl_matrix_set (LU_double, i, j, + gsl_matrix_float_get(LU, i, j)); } } - //LU decomposition + // LU decomposition. gsl_linalg_LU_decomp (LU_double, p, signum); - //copy float matrix to double + // Copy float matrix to double. for (size_t i=0; i<LU->size1; i++) { for (size_t j=0; j<LU->size2; j++) { - gsl_matrix_float_set (LU, i, j, gsl_matrix_get(LU_double, i, j)); + gsl_matrix_float_set (LU, i, j, + gsl_matrix_get(LU_double, i, j)); } } - //free matrix + // Free matrix. gsl_matrix_free (LU_double); return; } -//LU invert -void LUInvert (const gsl_matrix *LU, const gsl_permutation *p, gsl_matrix *inverse) -{ +// LU invert. +void LUInvert (const gsl_matrix *LU, const gsl_permutation *p, + gsl_matrix *inverse) { gsl_linalg_LU_invert (LU, p, inverse); return; } -void LUInvert (const gsl_matrix_float *LU, const gsl_permutation *p, gsl_matrix_float *inverse) -{ +void LUInvert (const gsl_matrix_float *LU, const gsl_permutation *p, + gsl_matrix_float *inverse) { gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2); - gsl_matrix *inverse_double=gsl_matrix_alloc (inverse->size1, inverse->size2); + gsl_matrix *inverse_double=gsl_matrix_alloc (inverse->size1, + inverse->size2); - //copy float matrix to double + // Copy float matrix to double. for (size_t i=0; i<LU->size1; i++) { for (size_t j=0; j<LU->size2; j++) { - gsl_matrix_set (LU_double, i, j, gsl_matrix_float_get(LU, i, j)); + gsl_matrix_set (LU_double, i, j, + gsl_matrix_float_get(LU, i, j)); } } - //LU decomposition + // LU decomposition. gsl_linalg_LU_invert (LU_double, p, inverse_double); - //copy float matrix to double + // Copy float matrix to double. for (size_t i=0; i<inverse->size1; i++) { for (size_t j=0; j<inverse->size2; j++) { - gsl_matrix_float_set (inverse, i, j, gsl_matrix_get(inverse_double, i, j)); + gsl_matrix_float_set (inverse, i, j, + gsl_matrix_get(inverse_double, + i, j)); } } - //free matrix + // Free matrix. gsl_matrix_free (LU_double); gsl_matrix_free (inverse_double); return; } -//LU lndet -double LULndet (gsl_matrix *LU) -{ +// LU lndet. +double LULndet (gsl_matrix *LU) { double d; d=gsl_linalg_LU_lndet (LU); return d; } -double LULndet (gsl_matrix_float *LU) -{ +double LULndet (gsl_matrix_float *LU) { gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2); double d; - //copy float matrix to double + // Copy float matrix to double. for (size_t i=0; i<LU->size1; i++) { for (size_t j=0; j<LU->size2; j++) { gsl_matrix_set (LU_double, i, j, gsl_matrix_float_get(LU, i, j)); } } - //LU decomposition + // LU decomposition. d=gsl_linalg_LU_lndet (LU_double); - //copy float matrix to double - /* - for (size_t i=0; i<LU->size1; i++) { - for (size_t j=0; j<LU->size2; j++) { - gsl_matrix_float_set (LU, i, j, gsl_matrix_get(LU_double, i, j)); - } - } - */ - //free matrix + // Free matrix gsl_matrix_free (LU_double); return d; } -//LU solve -void LUSolve (const gsl_matrix *LU, const gsl_permutation *p, const gsl_vector *b, gsl_vector *x) -{ +// LU solve. +void LUSolve (const gsl_matrix *LU, const gsl_permutation *p, + const gsl_vector *b, gsl_vector *x) { gsl_linalg_LU_solve (LU, p, b, x); return; } -void LUSolve (const gsl_matrix_float *LU, const gsl_permutation *p, const gsl_vector_float *b, gsl_vector_float *x) -{ +void LUSolve (const gsl_matrix_float *LU, const gsl_permutation *p, + const gsl_vector_float *b, gsl_vector_float *x) { gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2); gsl_vector *b_double=gsl_vector_alloc (b->size); gsl_vector *x_double=gsl_vector_alloc (x->size); - //copy float matrix to double + // Copy float matrix to double. for (size_t i=0; i<LU->size1; i++) { for (size_t j=0; j<LU->size2; j++) { - gsl_matrix_set (LU_double, i, j, gsl_matrix_float_get(LU, i, j)); + gsl_matrix_set (LU_double, i, j, + gsl_matrix_float_get(LU, i, j)); } } @@ -594,15 +679,15 @@ void LUSolve (const gsl_matrix_float *LU, const gsl_permutation *p, const gsl_ve gsl_vector_set (x_double, i, gsl_vector_float_get(x, i)); } - //LU decomposition + // LU decomposition. gsl_linalg_LU_solve (LU_double, p, b_double, x_double); - //copy float matrix to double + // Copy float matrix to double. for (size_t i=0; i<x->size; i++) { gsl_vector_float_set (x, i, gsl_vector_get(x_double, i)); } - //free matrix + // Free matrix. gsl_matrix_free (LU_double); gsl_vector_free (b_double); gsl_vector_free (x_double); @@ -610,8 +695,7 @@ void LUSolve (const gsl_matrix_float *LU, const gsl_permutation *p, const gsl_ve } -bool lapack_ddot(vector<double> &x, vector<double> &y, double &v) -{ +bool lapack_ddot(vector<double> &x, vector<double> &y, double &v) { bool flag=false; int incx=1; int incy=1; @@ -625,8 +709,7 @@ bool lapack_ddot(vector<double> &x, vector<double> &y, double &v) } -bool lapack_sdot(vector<float> &x, vector<float> &y, double &v) -{ +bool lapack_sdot(vector<float> &x, vector<float> &y, double &v) { bool flag=false; int incx=1; int incy=1; |