From 17deca2d54827a00df3ea4d98df700fc2b8ed777 Mon Sep 17 00:00:00 2001 From: xiangzhou Date: Sat, 20 Sep 2014 10:17:34 -0400 Subject: initial upload, version 0.95alpha --- lapack.cpp | 609 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 609 insertions(+) create mode 100644 lapack.cpp (limited to 'lapack.cpp') diff --git a/lapack.cpp b/lapack.cpp new file mode 100644 index 0000000..83d5290 --- /dev/null +++ b/lapack.cpp @@ -0,0 +1,609 @@ +/* + Genome-wide Efficient Mixed Model Association (GEMMA) + Copyright (C) 2011 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 . +*/ + +#include +#include +#include "gsl/gsl_vector.h" +#include "gsl/gsl_matrix.h" +#include "gsl/gsl_linalg.h" + +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 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 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); + + +//cholesky decomposition, A is distroyed +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."<data, &LDA, &INFO); + if (INFO!=0) {cout<<"Cholesky decomposition unsuccessful in lapack_cholesky_decomp."<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."<data, &LDA, &INFO); + if (INFO!=0) {cout<<"Cholesky decomposition unsuccessful in lapack_cholesky_decomp."<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."<data, &LDA, x->data, &LDB, &INFO); + if (INFO!=0) {cout<<"Cholesky solve unsuccessful in lapack_cholesky_solve."<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."<data, &LDA, x->data, &LDB, &INFO); + if (INFO!=0) {cout<<"Cholesky solve unsuccessful in lapack_cholesky_solve."<size1, LDB=B->size1, LDC=C->size2; + + if (*TransA=='N' || *TransA=='n') {M=A->size1; K1=A->size2;} + else if (*TransA=='T' || *TransA=='t') {M=A->size2; K1=A->size1;} + else {cout<<"need 'N' or 'T' in lapack_sgemm"<size2; K2=B->size1;} + else if (*TransB=='T' || *TransB=='t') {N=B->size1; K2=B->size2;} + else {cout<<"need 'N' or 'T' in lapack_sgemm"<size1!=(size_t)M || C->size2!=(size_t)N) {cout<<"C not compatible in lapack_sgemm"<size2, A->size1); + gsl_matrix_float_transpose_memcpy (A_t, A); + gsl_matrix_float *B_t=gsl_matrix_float_alloc (B->size2, B->size1); + gsl_matrix_float_transpose_memcpy (B_t, B); + 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); + gsl_matrix_float_transpose_memcpy (C, C_t); + + gsl_matrix_float_free (A_t); + gsl_matrix_float_free (B_t); + gsl_matrix_float_free (C_t); + return; +} + + + +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;} + else if (*TransA=='T' || *TransA=='t') {M=A->size2; K1=A->size1;} + else {cout<<"need 'N' or 'T' in lapack_dgemm"<size2; K2=B->size1;} + else if (*TransB=='T' || *TransB=='t') {N=B->size1; K2=B->size2;} + else {cout<<"need 'N' or 'T' in lapack_dgemm"<size1!=(size_t)M || C->size2!=(size_t)N) {cout<<"C not compatible in lapack_dgemm"<size2, A->size1); + gsl_matrix_transpose_memcpy (A_t, A); + gsl_matrix *B_t=gsl_matrix_alloc (B->size2, B->size1); + gsl_matrix_transpose_memcpy (B_t, B); + 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); + + gsl_matrix_transpose_memcpy (C, C_t); + + gsl_matrix_free (A_t); + gsl_matrix_free (B_t); + gsl_matrix_free (C_t); + 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) +{ + 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."<data, &LDA, eval->data, temp, &LWORK, &INFO); + // if (INFO!=0) {cout<<"Work space estimate unsuccessful in lapack_eigen_symmv."<data, &LDA, eval->data, WORK, &LWORK, &INFO); + if (INFO!=0) {cout<<"Eigen decomposition unsuccessful in lapack_eigen_symmv."<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' + 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."<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."<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."<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."<data, &LDA, eval->data, temp, &LWORK, &INFO); + // if (INFO!=0) {cout<<"Work space estimate unsuccessful in lapack_eigen_symmv."<data, &LDA, eval->data, WORK, &LWORK, &INFO); + if (INFO!=0) {cout<<"Eigen decomposition unsuccessful in lapack_eigen_symmv."<size1, LDA=A->size1, LDZ=A->size1, INFO, 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' + 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."<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."<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."<size1); + gsl_eigen_symmv (G, eval, U, w); + gsl_eigen_symmv_free (w); +#endif + /* + for (size_t i=0; isize; ++i) { + if (gsl_vector_get (eval, i)<1e-10) { +// cout<size; ++i) { + d+=gsl_vector_get(eval, i); + } + d/=(double)eval->size; + + return d; +} + + +//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_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; isize1; i++) { + for (size_t j=0; jsize2; 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 + for (size_t i=0; isize1; i++) { + for (size_t j=0; jsize2; j++) { + gsl_matrix_float_set(K, i, j, gsl_matrix_get(G_double, i, j)); + } + } + for (size_t i=0; isize1; i++) { + for (size_t j=0; jsize2; j++) { + gsl_matrix_float_set(U, i, j, gsl_matrix_get(U_double, i, j)); + } + } + for (size_t i=0; isize; i++) { + gsl_vector_float_set(eval, i, gsl_vector_get(eval_double, i)); + } + + //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; isize; ++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; + for (size_t i=0; isize; ++i) { + d+=gsl_vector_float_get(eval, i); + } + d/=(double)eval->size; + + return d; +} + + +double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty) +{ + double logdet_O=0.0; + +#ifdef WITH_LAPACK + lapack_cholesky_decomp(Omega); + for (size_t i=0; isize1; ++i) { + logdet_O+=log(gsl_matrix_get (Omega, i, i)); + } + logdet_O*=2.0; + lapack_cholesky_solve(Omega, Xty, OiXty); +#else + int status = gsl_linalg_cholesky_decomp(Omega); + if(status == GSL_EDOM) { + cout << "## non-positive definite matrix" << endl; + // exit(0); + } + + for (size_t i=0; isize1; ++i) { + logdet_O+=log(gsl_matrix_get (Omega, i, i)); + } + logdet_O*=2.0; + + 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); +#endif + + return logdet_O; +} + + +double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, gsl_vector_float *OiXty) +{ + double logdet_O=0.0; + +#ifdef WITH_LAPACK + lapack_float_cholesky_decomp(Omega); + for (size_t i=0; isize1; ++i) { + logdet_O+=log(gsl_matrix_float_get (Omega, i, i)); + } + logdet_O*=2.0; + lapack_float_cholesky_solve(Omega, Xty, OiXty); +#else + gsl_matrix *Omega_double=gsl_matrix_alloc (Omega->size1, Omega->size2); + double d; + for (size_t i=0; isize1; ++i) { + for (size_t j=0; jsize2; ++j) { + d=(double)gsl_matrix_float_get (Omega, i, j); + gsl_matrix_set (Omega_double, i, j, d); + } + } + + 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; isize1; ++i) { + for (size_t j=0; jsize2; ++j) { + d=gsl_matrix_get (Omega_double, i, j); + if (j==i) {logdet_O+=log(d);} + gsl_matrix_float_set (Omega, i, j, (float)d); + } + } + logdet_O*=2.0; + + 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_matrix_free (Omega_double); +#endif + + return logdet_O; +} + + +//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) +{ + gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2); + + //copy float matrix to double + for (size_t i=0; isize1; i++) { + for (size_t j=0; jsize2; j++) { + gsl_matrix_set (LU_double, i, j, gsl_matrix_float_get(LU, i, j)); + } + } + + //LU decomposition + gsl_linalg_LU_decomp (LU_double, p, signum); + + //copy float matrix to double + for (size_t i=0; isize1; i++) { + for (size_t j=0; jsize2; j++) { + gsl_matrix_float_set (LU, i, j, gsl_matrix_get(LU_double, i, j)); + } + } + + //free matrix + gsl_matrix_free (LU_double); + return; +} + + +//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) +{ + gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2); + gsl_matrix *inverse_double=gsl_matrix_alloc (inverse->size1, inverse->size2); + + //copy float matrix to double + for (size_t i=0; isize1; i++) { + for (size_t j=0; jsize2; j++) { + gsl_matrix_set (LU_double, i, j, gsl_matrix_float_get(LU, i, j)); + } + } + + //LU decomposition + gsl_linalg_LU_invert (LU_double, p, inverse_double); + + //copy float matrix to double + for (size_t i=0; isize1; i++) { + for (size_t j=0; jsize2; j++) { + gsl_matrix_float_set (inverse, i, j, gsl_matrix_get(inverse_double, i, j)); + } + } + + //free matrix + gsl_matrix_free (LU_double); + gsl_matrix_free (inverse_double); + return; +} + +//LU lndet +double LULndet (gsl_matrix *LU) +{ + double d; + d=gsl_linalg_LU_lndet (LU); + return d; +} + +double LULndet (gsl_matrix_float *LU) +{ + gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2); + double d; + + //copy float matrix to double + for (size_t i=0; isize1; i++) { + for (size_t j=0; jsize2; j++) { + gsl_matrix_set (LU_double, i, j, gsl_matrix_float_get(LU, i, j)); + } + } + + //LU decomposition + d=gsl_linalg_LU_lndet (LU_double); + + //copy float matrix to double + /* + for (size_t i=0; isize1; i++) { + for (size_t j=0; jsize2; j++) { + gsl_matrix_float_set (LU, i, j, gsl_matrix_get(LU_double, i, j)); + } + } + */ + //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) +{ + 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) +{ + 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 + for (size_t i=0; isize1; i++) { + for (size_t j=0; jsize2; j++) { + gsl_matrix_set (LU_double, i, j, gsl_matrix_float_get(LU, i, j)); + } + } + + for (size_t i=0; isize; i++) { + gsl_vector_set (b_double, i, gsl_vector_float_get(b, i)); + } + + for (size_t i=0; isize; i++) { + gsl_vector_set (x_double, i, gsl_vector_float_get(x, i)); + } + + //LU decomposition + gsl_linalg_LU_solve (LU_double, p, b_double, x_double); + + //copy float matrix to double + for (size_t i=0; isize; i++) { + gsl_vector_float_set (x, i, gsl_vector_get(x_double, i)); + } + + //free matrix + gsl_matrix_free (LU_double); + gsl_vector_free (b_double); + gsl_vector_free (x_double); + return; +} + + -- cgit v1.2.3